PTsampler {DPpackage} | R Documentation |
This function allows a user to generate a sample from a user-defined unormalized continuos distribution using the Polya tree sampler algorithm.
PTsampler(ltarget,dim.theta,mcmc=NULL,support=NULL,pts.options=NULL, status=TRUE,state=NULL)
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). |
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 (2009) for more details on the Polya tree sampler.
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.
Alejandro Jara <ajarav@udec.cl>
Tim Hanson <hanson@biostat.umn.edu>
Hanson, T., Monteiro, J.V.D, and Jara, A. (2009) The Polya Tree Sampler: Towards Efficient and Automatic Independent Metropolis Proposals. Technical Report.
## 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)