PTsampler: Polya Tree sampler function

In DPpackage: Bayesian Nonparametric Modeling in R

Description

This function allows a user to generate a sample from a user-defined unormalized continuos distribution using the Polya tree sampler algorithm.

Usage

PTsampler(ltarget,dim.theta,mcmc=NULL,support=NULL,pts.options=NULL,
	  status=TRUE,state=NULL)

Arguments

ltarget	a function giving the log of the target density.
dim.theta	an integer indicating the dimension of the target density.
mcmc	an optional list giving the MCMC parameters. The list must include the following integers: nburn giving the number of burn-in scans, nsave giving the total number of scans to be saved, and ndisplay giving the number of saved scans to be displayed on screen (the function reports on the screen when every ndisplay iterations have been carried out). Default values are 1000, 1000, and 100 for nburn, nsave, and ndisplay, respectively.
support	an optional matrix, of dimension dim.theta * npoints, giving the initial support points. By default the function generates 400 support points from a dim.theta normal distribution with mean 0 and diagonal covariance matrix with 1000 in the diagonal.
pts.options	an optional list of giving the parameters needed for the PTsampler algorithm. The list must include: nlevel (an integer giving the number of levels of the finite Polya tree approximation; default=5), tune1 (a double precision variable representing the standard deviation of the log-normal candidate distribution for the precision parameter of the Polya tree; default=1), delta (a double precision number indicating the maximum distance between the target and the approximation; default=0.2), max.warmup (an integer giving the maximum number of steps allowed for the warm-up phase; default=50000), minc (a double precision variable giving the minimum value allowed for the precision parameter of the Polya tree approximation; default=1), cpar0 (a double precision variable giving the initial value for the precision parameter of the Polya tree approximation; default=1000), and nadd (an integer variable giving the number of warm-up steps after convergence; default=1000).
status	a logical variable indicating whether this run is new (TRUE) or the continuation of a previous analysis (FALSE). In the latter case the current value of the parameters must be specified in the object state.
state	a list giving the starting points for the MCMC algorithm. The list must include: theta (a vector of dimension dim.theta of parameters), u (a Polya tree decomposition matrix), uinv (a matrix giving the inverse of the decompositon matrix), cpar (giving the value of the Polya tree precision parameter), support (a matrix giving the final support points), dim.theta (an integer giving the dimension of the problem), and L1 (a double precision number giving the final convergence criterion value).

Details

PTsampler produces a sample from a user-defined multivariate distribution using the Polya tree sampler algorithm. The algorithm constructs an independent proposal based on an approximation of the target density. The approximation is built from a set of support points and the predictive density of a finite multivariate Polya tree. In an initial warm-up phase, the support points are iteratively relocated to regions of higher support under the target distribution to minimize the distance between the target distribution and the Polya tree predictive distribution. In the sampling phase, samples from the final approximating mixture of finite Polya trees are used as candidates which are accepted with a standard Metropolis-Hastings acceptance probability. We refer to Hanson, Monteiro, and Jara (2011) for more details on the Polya tree sampler.

Value

An object of class PTsampler representing the MCMC sampler. Generic functions such as print, plot, and summary have methods to show the results of the fit.

The list state in the output object contains the current value of the parameters necessary to restart the analysis. If you want to specify different starting values to run multiple chains set status=TRUE and create the list state based on this starting values.

The object thetsave in the output list save.state contains the samples from the target density.

Author(s)

Alejandro Jara <[e-mail]>

Tim Hanson <[e-mail]>

References

Hanson, T., Monteiro, J.V.D, and Jara, A. (2011) The Polya Tree Sampler: Toward Efficient and Automatic Independent Metropolis-Hastings Proposals. Journal of Computational and Graphical Statistics, 20: 41-62.

See Also

PTdensity

Examples

## Not run: 

###############################
# EXAMPLE 1 (Dog Bowl)
###############################

# Target density

  target <- function(x,y)
  {
     out <- (-3/2)*log(2*pi)-0.5*(sqrt(x^2+y^2)-10)^2-
            0.5*log(x^2+y^2)
     exp(out)
  }	

  ltarget <- function(x)
  {
     out <- -0.5*((sqrt(x[1]^2+x[2]^2)-10)^2)-
             0.5*log(x[1]^2+x[2]^2)
     out
  }	

# MCMC

  mcmc <- list(nburn=5000,
               nsave=10000,
               ndisplay=500)

# Initial support points (optional)

  support <- cbind(rnorm(300,15,1),rnorm(300,15,1))

# Scanning the posterior

  fit <- PTsampler(ltarget,dim.theta=2,mcmc=mcmc,support=support)

  fit
  summary(fit)
  plot(fit,ask=FALSE)	

# Samples saved in 
# fit$save.state$thetasave
# Here is an example of how to use them
	
  par(mfrow=c(1,2))
  plot(acf(fit$save.state$thetasave[,1],lag=100))
  plot(acf(fit$save.state$thetasave[,1],lag=100))
	  
# Plotting resulting support points

  x1 <- seq(-15,15,0.2)	
  x2 <- seq(-15,15,0.2)	
  z <- outer(x1,x2,FUN="target")
  par(mfrow=c(1,1))
  image(x1,x2,z,xlab=expression(theta[1]),ylab=expression(theta[2]))
  points(fit$state$support,pch=19,cex=0.25)

# Plotting the samples from the target density

  par(mfrow=c(1,1))
  image(x1,x2,z,xlab=expression(theta[1]),ylab=expression(theta[2]))
  points(fit$save.state$thetasave,pch=19,cex=0.25)

# Re-starting the chain from the last sample

  state <- fit$state
  fit <- PTsampler(ltarget,dim.theta=2,mcmc=mcmc,
                   state=state,status=FALSE)


###############################
# EXAMPLE 2 (Ping Pong Paddle)
###############################

  bivnorm1 <- function(x1,x2)
  {
       eval <- (x1)^2+(x2)^2 
       logDET <-  0
       logPDF <- -(2*log(2*pi)+logDET+eval)/2
       out <- exp(logPDF)
       out
  }

  bivnorm2 <- function(x1,x2)
  {
       mu <- c(-3,-3)
       sigmaInv <- matrix(c(5.263158,-4.736842,
                           -4.736842,5.263158),
                            nrow=2,ncol=2)
       eval <- (x1-mu[1])^2*sigmaInv[1,1]+
               2*(x1-mu[1])*(x2-mu[2])*sigmaInv[1,2]+
               (x2-mu[2])^2*sigmaInv[2,2] 
       logDET <-  -1.660731
       logPDF <- -(2*log(2*pi)+logDET+eval)/2
       out <- exp(logPDF)
       out
  }

  bivnorm3 <- function(x1,x2)
  {
       mu <- c(2,2)
       sigmaInv <- matrix(c(5.263158,4.736842,
                            4.736842,5.263158),
                            nrow=2,ncol=2)
       eval <- (x1-mu[1])^2*sigmaInv[1,1]+
               2*(x1-mu[1])*(x2-mu[2])*sigmaInv[1,2]+
               (x2-mu[2])^2*sigmaInv[2,2] 
       logDET <-  -1.660731
       logPDF <- -(2*log(2*pi)+logDET+eval)/2
       out <- exp(logPDF)
       out
  }

  target <- function(x,y)
  {
     out <- 0.34*bivnorm1(x,y)+
	    0.33*bivnorm2(x,y)+
	    0.33*bivnorm3(x,y)
     out
  }	

  ltarget <- function(theta)
  {
     out <- 0.34*bivnorm1(x1=theta[1],x2=theta[2])+
	    0.33*bivnorm2(x1=theta[1],x2=theta[2])+
	    0.33*bivnorm3(x1=theta[1],x2=theta[2])
     log(out)
  }	


# MCMC

  mcmc <- list(nburn=5000,
               nsave=10000,
               ndisplay=500)

# Initial support points (optional)

  support <- cbind(rnorm(300,6,1),rnorm(300,6,1))

# Scanning the posterior

  fit <- PTsampler(ltarget,dim.theta=2,mcmc=mcmc,support=support)

  fit
  summary(fit)
  plot(fit,ask=FALSE)	

# Samples saved in 
# fit$save.state$thetasave
# Here is an example of how to use them
	
  par(mfrow=c(1,2))
  plot(acf(fit$save.state$thetasave[,1],lag=100))
  plot(acf(fit$save.state$thetasave[,1],lag=100))
	  
# Plotting resulting support points

  x1 <- seq(-6,6,0.05)	
  x2 <- seq(-6,6,0.05)	
  z <- outer(x1,x2,FUN="target")
  par(mfrow=c(1,1))
  image(x1,x2,z,xlab=expression(theta[1]),ylab=expression(theta[2]))
  points(fit$state$support,pch=19,cex=0.25)

# Plotting the samples from the target density

  par(mfrow=c(1,1))
  image(x1,x2,z,xlab=expression(theta[1]),ylab=expression(theta[2]))
  points(fit$save.state$thetasave,pch=19,cex=0.25)



## End(Not run)

DPpackage documentation built on May 1, 2019, 10:23 p.m.