trans.dens: Functions for transdimensional MCMC

In HI: Simulation from Distributions Supported by Nested Hyperplanes

Functions for transdimensional MCMC

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

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 [e-mail], 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

#### ==> 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 April 30, 2022, 1:06 a.m.