Last data update: 2014.03.03

R: Functions for transdimensional MCMC
 trans.dens R Documentation

## 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 GPetris@uark.edu, 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)

```

### Results

```
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(HI)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/HI/trans.dens.Rd_%03d_medium.png", width=480, height=480)
> ### Name: trans.dens
> ### Title: Functions for transdimensional MCMC
> ### Aliases: trans.dens trans.up trans2 transUp2 transBack2
> ### Keywords: distribution
>
> ### ** 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

0     1     2
0.330 0.348 0.322
> 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

FALSE  TRUE
0.496 0.504
> 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)
>
>
>
>
>
>
> dev.off()
null device
1
>

```