[R] initial points for arms in package HI

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


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

## 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



Christoph Buser <buser@stat.math.ethz.ch> Seminar fuer Statistik, LEO C13
ETH (Federal Inst. Technology)	8092 Zurich	 SWITZERLAND
phone: x-41-44-632-4673		fax: 632-1228

http://stat.ethz.ch/~buser/

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 03:07:38 2005

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