From: <plummer_at_iarc.fr>

Date: Wed 20 Jul 2005 - 07:18:27 EST

This message and its attachments are strictly confidential. ...{{dropped}}

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 07:22:31 2005

Date: Wed 20 Jul 2005 - 07:18:27 EST

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

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 07:22:31 2005

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