# 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)
# }
Run the code above in your browser using DataLab