Polya Tree sampler function

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

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 <atjara@uc.cl>

Tim Hanson <hansont@stat.sc.edu>

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)