arms: Function to perform Adaptive Rejection Metropolis Sampling In HI: Simulation from distributions supported by nested hyperplanes

Description

Generates a sequence of random variables using ARMS. For multivariate densities, ARMS is used along randomly selected straight lines through the current point.

Usage

 `1` ```arms(y.start, myldens, indFunc, n.sample, ...) ```

Arguments

 `y.start` Initial point `myldens` Univariate or multivariate log target density `indFunc` Indicator function of the convex support of the target density `n.sample` Desired sample size `...` Parameters passed to `myldens` and `indFunc`

Details

Strictly speaking, the support of the target density must be a bounded convex set. When this is not the case, the following tricks usually work. If the support is not bounded, restrict it to a bounded set having probability practically one. A workaround, if the support is not convex, is to consider the convex set generated by the support and define `myldens` to return `log(.Machine\$double.xmin)` outside the true support (see the last example.)

The next point is generated along a randomly selected line through the current point using arms.

Make sure the value returned by `myldens` is never smaller than `log(.Machine\$double.xmin)`, to avoid divisions by zero.

Value

An `n.sample` by `length(y.start)` matrix, whose rows are the sampled points.

Note

The function is based on original C code by Wally Gilks for the univariate case, http://www.mrc-bsu.cam.ac.uk/pub/methodology/adaptive_rejection/.

Author(s)

Giovanni Petris [email protected]

References

Gilks, W.R., Best, N.G. and Tan, K.K.C. (1995) Adaptive rejection Metropolis sampling within Gibbs sampling (Corr: 97V46 p541-542 with Neal, R.M.), Applied Statistics 44:455–472.

Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94``` ```#### ==> Warning: running the examples may take a few minutes! <== #### ### Univariate densities ## Unif(-r,r) y <- arms(runif(1,-1,1), function(x,r) 1, function(x,r) (x>-r)*(x-7)*((x-mean)<7), 5000, mean=10) summary(y); hist(y,prob=TRUE,main="Gaussian(m,1); m=10") curve(dnorm(x,mean=10),3,17,add=TRUE) ## Exponential(1) y <- arms(5, function(x) -x, function(x) (x>0)*(x<70), 5000) summary(y); hist(y,prob=TRUE,main="Exponential(1)") curve(exp(-x),0,8,add=TRUE) ## Gamma(4.5,1) y <- arms(runif(1,1e-4,20), function(x) 3.5*log(x)-x, function(x) (x>1e-4)*(x<20), 5000) summary(y); hist(y,prob=TRUE,main="Gamma(4.5,1)") curve(dgamma(x,shape=4.5,scale=1),1e-4,20,add=TRUE) ## Gamma(0.5,1) (this one is not log-concave) y <- arms(runif(1,1e-8,10), function(x) -0.5*log(x)-x, function(x) (x>1e-8)*(x<10), 5000) summary(y); hist(y,prob=TRUE,main="Gamma(0.5,1)") curve(dgamma(x,shape=0.5,scale=1),1e-8,10,add=TRUE) ## Beta(.2,.2) (this one neither) y <- arms(runif(1), function(x) (0.2-1)*log(x)+(0.2-1)*log(1-x), function(x) (x>1e-5)*(x<1-1e-5), 5000) summary(y); hist(y,prob=TRUE,main="Beta(0.2,0.2)") curve(dbeta(x,0.2,0.2),1e-5,1-1e-5,add=TRUE) ## Triangular y <- arms(runif(1,-1,1), function(x) log(1-abs(x)), function(x) abs(x)<1, 5000) summary(y); hist(y,prob=TRUE,ylim=c(0,1),main="Triangular") curve(1-abs(x),-1,1,add=TRUE) ## Multimodal examples (Mixture of normals) lmixnorm <- function(x,weights,means,sds) { log(crossprod(weights, exp(-0.5*((x-means)/sds)^2 - log(sds)))) } y <- arms(0, lmixnorm, function(x,...) (x>(-100))*(x<100), 5000, weights=c(1,3,2), means=c(-10,0,10), sds=c(1.5,3,1.5)) summary(y); hist(y,prob=TRUE,main="Mixture of Normals") curve(colSums(c(1,3,2)/6*dnorm(matrix(x,3,length(x),byrow=TRUE),c(-10,0,10),c(1.5,3,1.5))), par("usr")[1], par("usr")[2], add=TRUE) ### Bivariate densities ## Bivariate standard normal y <- arms(c(0,2), function(x) -crossprod(x)/2, function(x) (min(x)>-5)*(max(x)<5), 500) plot(y, main="Bivariate standard normal", asp=1) ## Uniform in the unit square y <- arms(c(0.2,.6), function(x) 1, function(x) (min(x)>0)*(max(x)<1), 500) plot(y, main="Uniform in the unit square", asp=1) polygon(c(0,1,1,0),c(0,0,1,1)) ## Uniform in the circle of radius r y <- arms(c(0.2,0), function(x,...) 1, function(x,r2) sum(x^2)1) ) 0 else 1 y <- arms(c(0.2,0.2), function(x) 1, simp, 500) plot(y, xlim=c(0,1), ylim=c(0,1), main="Uniform in the simplex", asp=1) polygon(c(0,1,0), c(0,0,1)) ## A bimodal distribution (mixture of normals) bimodal <- function(x) { log(prod(dnorm(x,mean=3))+prod(dnorm(x,mean=-3))) } y <- arms(c(-2,2), bimodal, function(x) all(x>(-10))*all(x<(10)), 500) plot(y, main="Mixture of bivariate Normals", asp=1) ## A bivariate distribution with non-convex support support <- function(x) { return(as.numeric( -1 < x[2] && x[2] < 1 && -2 < x[1] && ( x[1] < 1 || crossprod(x-c(1,0)) < 1 ) ) ) } Min.log <- log(.Machine\$double.xmin) + 10 logf <- function(x) { if ( x[1] < 0 ) return(log(1/4)) else if (crossprod(x-c(1,0)) < 1 ) return(log(1/pi)) return(Min.log) } x <- as.matrix(expand.grid(seq(-2.2,2.2,length=40),seq(-1.1,1.1,length=40))) y <- sapply(1:nrow(x), function(i) support(x[i,])) plot(x,type='n',asp=1) points(x[y==1,],pch=1,cex=1,col='green') z <- arms(c(0,0), logf, support, 1000) points(z,pch=20,cex=0.5,col='blue') polygon(c(-2,0,0,-2),c(-1,-1,1,1)) curve(-sqrt(1-(x-1)^2),0,2,add=TRUE) curve(sqrt(1-(x-1)^2),0,2,add=TRUE) sum( z[,1] < 0 ) # sampled points in the square sum( apply(t(z)-c(1,0),2,crossprod) < 1 ) # sampled points in the circle ```

HI documentation built on May 29, 2017, 5:38 p.m.