PTsampler(ltarget,dim.theta,mcmc=NULL,support=NULL,pts.options=NULL,
status=TRUE,state=NULL)nburn giving the number of burn-in
scans, nsave giving
the total numbernlevel (an integer giving the number
of levels of the finite Polya tree approximation; default=5),
TRUE) or the
continuation of a previous analysis (FALSE). In the latter case
the current value of the parameters must be specifietheta (a vector of dimension dim.theta of parameters),
u (a Polya tree decomposition matrix), uinv (a matrix giving
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.
PTdensity###############################
# 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)Run the code above in your browser using DataLab