[R] initial points for arms in package HI

[plummer@iarc.fr]

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. ...{{dropped}}