HI (version 0.5)

trans.dens: 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

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

Details

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

References

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

Examples

# NOT RUN {
#### ==> 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)

# }