Re: [R] initial points for arms in package HI

From: Christoph Buser <buser_at_stat.math.ethz.ch>
Date: Wed 20 Jul 2005 - 21:23:53 EST

Dear Martyn

Thank you for your fast and helpful answer. It showed me to my wrong thinking.

I confused the previous value of the Markov Chain with the initial values to construct the envelope.

In the original C code by W. Gilks there are the arguments "xinit" and "xprev". The first is used to construct the envelope, the second is the previous value of the Markov chain. There are comments from the author how to set the initial values "xinit". It is important to choose these values independent from the current parameter, being updated (in the case of not log-concave functions).

In the R code and the C code behind the function arms the argument "xinit" is missing, since it is set in the function itself and the user can not change it.

Since I confused the two arguments I didn't realize that "y.start" is the value "xprev", the state of the Markov chain and set this value 1 or runif(...). But for an implementation of a Gibbs Sampler this makes no sense.

My example was a little bit misleading, since I simulated 1000 points from the same distribution. In reality my distribution is changing at each step and I always need only 1 new value in my Gibbs Sampler.
But by implementing it correctly with the previous value of the chain it should work in my context.

Best regards,

Christoph

plummer@iarc.fr writes:
> Quoting Christoph Buser <buser@stat.math.ethz.ch>:
>
> > Dear R-users
> >
> > I have a problem choosing initial points for the function arms()
> > in the package HI
> > I intend to implement a Gibbs sampler and one of my conditional
> > distributions is nonstandard and not logconcave.
> > Therefore I'd like to use arms.
> >
> > But there seem to be a strong influence of the initial point
> > y.start. To show the effect I constructed a demonstration
> > example. It is reproducible without further information.
> >
> > Please note that my target density is not logconcave.
> >
> > Thanks for all comments or ideas.
> >
> > Christoph Buser
>
> Dear Christoph,
>
> There is a Metropolis step at each iteration of the ARMS sampler, in
> which it may choose to reject the proposed move to a new point and stick
> at the current point (This is what the "M" in "ARMS" stands for) If you
> do repeated calls to arms with the same starting point, then the
> iterations where the Metropolis step rejects a move will create a spike
> in the sample density at your initial value. If you use a uniform random
> starting point, then your sample density will be a mixture of the
> target distribution (Metropolis accepts move) and a uniform distribution
> (Metropolis rejects move).
>
> You should be doing something like this:
>
> res1 <- arms(runif(1,0,100), logDichteGam, function(x) (x>0)&(x<100), 1)
> for(i in 2:1000)
> res1[i] <- arms(res1[i-1], logDichteGam, function(x) (x>0)&(x<100), 1)
>
> i.e., using each sampled point as the starting value for the next
> iteration. The sequence of values in res1 will then be a correlated
> sample from the given distribution:
>
> acf(res1)
>
> The bottom line is that you can't use ARMS to draw a single sample
> from a non-log-concave density.
>
> If you are still worried about using ARMS, you can verify your results
> using the random walk Metropolis sampler (MCMCmetrop1R) in the package
> MCMCpack.
>
> Martyn
>
> > ## R Code:
> >
> > library(HI)
> > ## parameter for the distribution
> > para <- 0.1
> >
> > ## logdensity
> > logDichteGam <- function(x, u = para, v = para) {
> > -(u*x + v*1/x) - log(x)
> > }
> > ## density except for the constant
> > propDichteGam <- function(x, u = para, v = para) {
> > exp(-(u*x + v*1/x) - log(x))
> > }
> > ## calculating the constant
> > (c <- integrate(propDichteGam, 0, 1000, rel.tol = 10^(-12))$value)
> > ## density
> > DichteGam <- function(x, u = para, v = para) {
> > exp(-(u*x + v*1/x) - log(x))/c
> > }
> >
> > ## calculating 1000 values by repeating a call of arms (this would
> > ## be the situation in an Gibbs Sample. Of course in a Gibbs sampler
> > ## the distribution would change. This is only for demonstration
> > res1 <- NULL
> > for(i in 1:1000)
> > res1[i] <- arms(runif(1,0,100), logDichteGam, function(x) (x>0)&(x<100), 1)
> >
> > ## Generating a sample of thousand observations with 1 call of arms
> > res2 <- arms(runif(1,0,100), logDichteGam, function(x) (x>0)&(x<100), 1000)
> >
> > ## Plot of the samples
> > mult.fig(4)
> > plot(res1, log = "y")
> > plot(res2, log = "y")
> > hist(res1, freq = FALSE, xlim = c(0,4), breaks = seq(0,100,by = 0.1),
> > ylim = c(0,1))
> > curve(DichteGam, 0,4, add = TRUE, col = 2)
> > hist(res2, freq = FALSE, xlim = c(0,4), breaks = seq(0,100,by = 0.1),
> > ylim = c(0,1))
> > curve(DichteGam, 0,4, add = TRUE, col = 2)
> >
> >
> > ## If we repeat the procedure, using the fix intial value 1,
> > ## the situation is even worse
> > res3 <- NULL
> > for(i in 1:1000)
> > res3[i] <- arms(1, logDichteGam, function(x) (x>0)&(x<100), 1)
> >
> > ## Generating a sample of thousand observations with 1 call of arms
> > res4 <- arms(1, logDichteGam, function(x) (x>0)&(x<100), 1000)
> >
> > ## Plot of the samples
> > par(mfrow = c(2,2))
> > plot(res3, log = "y")
> > plot(res4, log = "y")
> > hist(res3, freq = FALSE, xlim = c(0,4), breaks = seq(0,100,by = 0.1),
> > ylim = c(0,1))
> > curve(DichteGam, 0,4, add = TRUE, col = 2)
> > hist(res4, freq = FALSE, xlim = c(0,4), breaks = seq(0,100,by = 0.1),
> > ylim = c(0,1))
> > curve(DichteGam, 0,4, add = TRUE, col = 2)
> >
> >
> > ## If I generate the sample in a for-loop (one by one) I do not
> > ## get the correct density. But this is exactly the situation in
> > ## my Gibbs Sampler. Therfore I am concerned about the correct
> > ## application of arms
> >
>
>
>
> -----------------------------------------------------------------------
> This message and its attachments are strictly confidential. If you are
> not the intended recipient of this message, please immediately notify
> the sender and delete it. Since its integrity cannot be guaranteed,
> its content cannot involve the sender's responsibility. Any misuse,
> any disclosure or publication of its content, either whole or partial,
> is prohibited, exception made of formally approved use
> -----------------------------------------------------------------------



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Wed Jul 20 21:28:33 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:33:51 EST