trans.dens: Functions for transdimensional MCMC In HI: Simulation from distributions supported by nested hyperplanes

Description

Computes the value of the 'g' density at a given point and, optionally, returns the backtransformed point and the model to which the point belongs.

Usage

 ```1 2 3 4 5``` ```trans.dens(y, ldens.list, which.models, ..., back.transform=F) trans.up(x, ldens.list, which.models, ...) trans2(y, ldens.list, k, ...) transUp2(y, ldens.list, k, ...) transBack2(y, ldens.list, k, ...) ```

Arguments

 `y` Vector or matrix of points (by row) at which the density of the absolutely continuous auxiliary distribution has to be evaluated `x` Vector or matrix of points corresponding to `y` (see examples) `ldens.list` List of densities (of submodels) `which.models` List of integers, in increasing order, giving the number of components to be dropped when evaluating the density in `which.models` in the corresponding position. A first element equal 0 (full model) is added if not already present `back.transform` Logical that determines the output `k` Difference between the dimension of the larger model and the dimension of the smaller model `...` Other arguments passed to the functions in `ldens.list`

Details

See the reference for details. The functions with the 2 in the name operate on pairs of models only.

Value

If `back.transform=F`, `trans.dens` returns the density of the absolutely continuous auxiliary distribution evaluated at the point(s) `y`.\ If `back.transform=T`, `trans.dens` returns in addition the point `x` corresponding to `y` in the original space and the index of the subspace to which `x` belongs.\ `trans.up` is a (stochastic) right inverse of the correspondence between `y` and `x`

Author(s)

Giovanni Petris [email protected], Luca Tardella

References

Petris \& Tardella, A geometric approach to transdimensional Markov chain Monte Carlo. The Canadian Journal of Statistics, vol.31, n.4, (2003).

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``` ```#### ==> Warning: running the examples may take a few minutes! <== #### ### Generate a sample from a mixture of 0,1,2-dim standard normals ldens.list <- list(f0 = function(x) sum(dnorm(x,log=TRUE)), f1 = function(x) dnorm(x,log=TRUE), f2 = function() 0) trans.mix <- function(y) { trans.dens(y, ldens.list=ldens.list, which.models=0:2) } trans.rmix <- arms(c(0,0), trans.mix, function(x) crossprod(x)<1e4, 500) rmix <- trans.dens(y=trans.rmix, ldens.list=ldens.list, which.models=0:2, back.transform = TRUE) table(rmix[,2])/nrow(rmix) # should be about equally distributed plot(trans.rmix,col=rmix[,2]+3,asp=1, xlab="y.1", ylab="y.2", main="A sample from the auxiliary continuous distribution") x <- rmix[,-(1:2)] plot(x, col=rmix[,2]+3, asp=1, main="The sample transformed back to the original space") ### trans.up as a right inverse of trans.dens set.seed(6324) y <- trans.up(x, ldens.list, 0:2) stopifnot(all.equal(x, trans.dens(y, ldens.list, 0:2, back.transform=TRUE)[,-(1:2)])) ### More trans.up z <- trans.up(matrix(0,1000,2), ldens.list, 0:2) plot(z,asp=1,col=5) # should look uniform in a circle corresponding to model 2 z <- trans.up(cbind(runif(1000,-3,3),0), ldens.list, 0:2) plot(z,asp=1,col=4) # should look uniform in a region corresponding to model 1 ### trans2, transBack2 ldens.list <- list(f0 = function(x) sum(dnorm(x,log=TRUE)), f1 = function(x) dnorm(x,log=TRUE)) trans.mix <- function(y) { trans2(y, ldens.list=ldens.list, k=1)[-2] } trans.rmix <- arms(c(0,0), trans.mix, function(x) crossprod(x)<1e2, 1000) rmix <- transBack2(y=trans.rmix, ldens.list=ldens.list, k=1) table(rmix[,2]==0)/nrow(rmix) # should be about equally distributed plot(trans.rmix,col=(rmix[,2]==0)+3,asp=1, xlab="y.1", ylab="y.2", main="A sample from the auxiliary continuous distribution") plot(rmix, col=(rmix[,2]==0)+3, asp=1, main="The sample transformed back to the original space") ### trunsUp2 z <- t(sapply(1:1000, function(i) transUp2(c(-2+0.004*i,0), ldens.list, 1))) plot(z,asp=1,col=2) ```

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