[R] initial points for arms in package HI

[Christoph Buser]

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 at iarc.fr writes:
 > Quoting Christoph Buser <buser at 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
 > -----------------------------------------------------------------------