Approximate Bayesian Computation (ABC) for the Ising/Potts model

ABC_intro2.R

set.seed(13)
y <- rnorm(n=5, mean=1, sd=2/3)
n <- length(y)
s_sq <- sum((y-1)^2)/n
# weakly informative priors:
nu0 <- 1
s0 <- 1
post_nu <- nu0 + n
post_ssd <- (nu0 * s0^2 + n*s_sq)/2
post_tau <- rgamma(10000, post_nu/2, post_ssd)
hist(post_tau, main="", xlab=expression(paste(pi,"(",tau, "|y)")), col="lightgreen", freq=F,xlim=c(0,7))
lines(density(post_tau), col=4, lwd=3, lty=2)
abline(v=(2/3)^{-2}, lty=3, col=2, lwd=3)
dev.copy2eps(file="postTau_exact.eps")

# now with ABC:
prop_tau <- rgamma(10000, nu0/2, 0.5*nu0*s0^2)
pseudoMx <- matrix(nrow=10000,ncol=n)
for (i in 1:n) {
  pseudoMx[,i] <- rnorm(10000, 1, 1/sqrt(prop_tau))
}
pseudoSSD <- rowSums((pseudoMx - 1)^2)/n
ps_norm <- abs(pseudoSSD - s_sq)
epsilon <- sort(ps_norm)[2000]
prop_keep <- prop_tau[ps_norm <= epsilon]

hist(prop_tau,breaks=50,freq=FALSE,col="lightgreen",main="", xlab=expression(pi(tau)),ylab="",xlim=c(0,7))
curve(dgamma(x, nu0/2, 0.5*nu0*s0^2), col="darkblue", lwd=3, lty=2, add=TRUE)
abline(v=(2/3)^{-2}, lty=3, col=2, lwd=3)
dev.copy2eps(file="priorTau.eps")

hist(prop_keep,breaks=50,freq=FALSE,col="lightgreen",main="", xlab=expression(paste(pi[epsilon],"(",tau, " | ",delta(s(w),s(y)) < epsilon,")")),ylab="",xlim=c(0,7))
curve(dgamma(x, post_nu/2, post_ssd), col="darkblue", lwd=2, add=TRUE)
abline(v=(2/3)^{-2}, lty=3, col=2, lwd=3)
dev.copy2eps(file="abc_postTau.eps")

# effect of epsilon
epsilon <- sort(ps_norm)[7500]
prop_keep <- prop_tau[ps_norm <= epsilon]
hist(prop_keep,breaks=20,freq=FALSE,col="lightgreen",main="", xlab="",ylab="",xlim=c(0,7))
curve(dgamma(x, post_nu/2, post_ssd), col="darkblue", lwd=2, add=TRUE)
abline(v=(2/3)^{-2}, lty=3, col=2, lwd=3)
dev.copy2eps(file="abcTau_post5000.eps")

epsilon <- sort(ps_norm)[5000]
prop_keep <- prop_tau[ps_norm <= epsilon]
hist(prop_keep,breaks=20,freq=FALSE,col="lightgreen",main="", xlab="",ylab="",xlim=c(0,7))
curve(dgamma(x, post_nu/2, post_ssd), col="darkblue", lwd=2, add=TRUE)
abline(v=(2/3)^{-2}, lty=3, col=2, lwd=3)
dev.copy2eps(file="abcTau_post3000.eps")

epsilon <- sort(ps_norm)[100]
prop_keep <- prop_tau[ps_norm <= epsilon]
hist(prop_keep,breaks=20,freq=FALSE,col="lightgreen",main="", xlab="",ylab="",xlim=c(0,7))
curve(dgamma(x, post_nu/2, post_ssd), col="darkblue", lwd=2, add=TRUE)
abline(v=(2/3)^{-2}, lty=3, col=2, lwd=3)
dev.copy2eps(file="abcTau_post20.eps")

epsilon <- sort(ps_norm)[10]
prop_keep <- prop_tau[ps_norm <= epsilon]
hist(prop_keep,breaks=20,freq=FALSE,col="lightgreen",main="", xlab="",ylab="",xlim=c(0,7))
curve(dgamma(x, post_nu/2, post_ssd), col="darkblue", lwd=2, add=TRUE)
abline(v=(2/3)^{-2}, lty=3, col=2, lwd=3)
dev.copy2eps(file="abcTau_post5.eps")

# concentration of measure
n <- 25
y <- rnorm(n=n, mean=1, sd=2/3)
s_sq <- sum((y-1)^2)/n
post_nu <- nu0 + n
post_ssd <- (nu0 * s0^2 + n*s_sq)/2
prop_tau <- rgamma(10000, nu0/2, 0.5*nu0*s0^2)
pseudoMx <- matrix(nrow=10000,ncol=n)
for (i in 1:n) {
  pseudoMx[,i] <- rnorm(10000, 1, 1/sqrt(prop_tau))
}
pseudoSSD <- rowSums((pseudoMx - 1)^2)/n
ps_norm <- abs(pseudoSSD - s_sq)
epsilon <- sort(ps_norm)[1000]
prop_keep <- prop_tau[ps_norm <= epsilon]
hist(prop_keep,breaks=20,freq=FALSE,col="lightgreen",main="", xlab="",ylab="",xlim=c(0,4),ylim=c(0,1))
curve(dgamma(x, post_nu/2, post_ssd), col="darkblue", lwd=2, add=TRUE)
abline(v=(2/3)^{-2}, lty=3, col=2, lwd=3)
dev.copy2eps(file="abcTau_n25.eps")

n <- 500
y <- rnorm(n=n, mean=1, sd=2/3)
s_sq <- sum((y-1)^2)/n
post_nu <- nu0 + n
post_ssd <- (nu0 * s0^2 + n*s_sq)/2
prop_tau <- rgamma(10000, nu0/2, 0.5*nu0*s0^2)
pseudoMx <- matrix(nrow=10000,ncol=n)
for (i in 1:n) {
  pseudoMx[,i] <- rnorm(10000, 1, 1/sqrt(prop_tau))
}
pseudoSSD <- rowSums((pseudoMx - 1)^2)/n
ps_norm <- abs(pseudoSSD - s_sq)
epsilon <- sort(ps_norm)[200]
prop_keep <- prop_tau[ps_norm <= epsilon]
hist(prop_keep,breaks=20,freq=FALSE,col="lightgreen",main="", xlab="",ylab="",xlim=c(0,4),ylim=c(0,3))
curve(dgamma(x, post_nu/2, post_ssd), col="darkblue", lwd=2, add=TRUE)
abline(v=(2/3)^{-2}, lty=3, col=2, lwd=3)
dev.copy2eps(file="abcTau_n500.eps")

n <- 10000
y <- rnorm(n=n, mean=1, sd=2/3)
s_sq <- sum((y-1)^2)/n
post_nu <- nu0 + n
post_ssd <- (nu0 * s0^2 + n*s_sq)/2
prop_tau <- rgamma(10000, nu0/2, 0.5*nu0*s0^2)
pseudoMx <- matrix(nrow=10000,ncol=n)
for (i in 1:n) {
  pseudoMx[,i] <- rnorm(10000, 1, 1/sqrt(prop_tau))
}
pseudoSSD <- rowSums((pseudoMx - 1)^2)/n
ps_norm <- abs(pseudoSSD - s_sq)
epsilon <- sort(ps_norm)[100]
prop_keep <- prop_tau[ps_norm <= epsilon]
hist(prop_keep,breaks=5,freq=FALSE,col="lightgreen",main="", xlab="",ylab="",xlim=c(0,4),ylim=c(0,9))
curve(dgamma(x, post_nu/2, post_ssd), col="darkblue", lwd=2, add=TRUE)
abline(v=(2/3)^{-2}, lty=3, col=2, lwd=3)
dev.copy2eps(file="abcTau_n10k.eps")

# for 1 million observations: multiple threads!
n <- 1e6
y <- rnorm(n=n, mean=1, sd=2/3)
s_sq <- sum((y-1)^2)/n
post_nu <- nu0 + n
post_ssd <- (nu0 * s0^2 + n*s_sq)/2
prop_tau <- rgamma(10000, nu0/2, 0.5*nu0*s0^2)

if (file.exists("pseudoSSD.rda")) {
  load("pseudoSSD.rda")
  load("data_n1mil.rda")
  s_sq <- sum((y-1)^2)/n
  post_nu <- nu0 + n
  post_ssd <- (nu0 * s0^2 + n*s_sq)/2
} else {
  library(doParallel)
  registerDoParallel(cores=min(4, detectCores()))

  pseudoSSD <- foreach (i=1:length(prop_tau), .combine=c, .multicombine=TRUE) %dopar% {
    pseudo <- rnorm(n, 1, 1/sqrt(prop_tau[i]))
    sum((pseudo - 1)^2)/n
  }
  save(pseudoSSD, file="pseudoSSD.rda")
}
ps_norm <- abs(pseudoSSD - s_sq)
epsilon <- sort(ps_norm)[15]
prop_keep <- prop_tau[ps_norm <= epsilon]
hist(prop_keep,breaks=5,freq=FALSE,col="lightgreen",main="", xlab="",ylab="",xlim=c(0,4))
curve(dgamma(x, post_nu/2, post_ssd), col="darkblue", lwd=2, add=TRUE)
abline(v=(2/3)^{-2}, lty=3, col=2, lwd=3)
dev.copy2eps(file="abcTau_n1mil.eps")

library(bayess)
data(Menteith)
image(as.matrix(Menteith), col=gray(1:256/256), useRaster=T)
#dev.copy2eps(file="Menteith.eps")

library(bayesImageS)
y <- as.matrix(Menteith)
mask <- matrix(1, nrow=nrow(y), ncol=ncol(y))
neigh <- getNeighbors(mask, c(2,2,0,0))
block <- getBlocks(mask, 2)
# weakly-informative priors for the mixture components
priors <- list()
priors$k <- 6
priors$mu <- rep(256/2, priors$k)
priors$mu.sd <- rep(256/6,priors$k)
priors$sigma <- rep(256/6,priors$k)
priors$sigma.nu <- rep(6, priors$k)
priors$beta <- c(0,2)
# ABC-MCMC
#system.time(mcmc_res <- mcmcPotts(as.vector(y), neigh, block, c(0,2), niter=2000,nburn=1000,
 #                                 priors, mh = list(algorithm="abc", auxiliary=200, aux_swap=TRUE)))

# ABC-SMC
if (file.exists("resSMC-ABC.rda")) {
  load("resSMC-ABC.rda")
} else {
  system.time(smc_res <- smcPotts(as.vector(y), neigh, block, param=list(npart=2000,nstat=5),priors=priors))
  smc_res$priors <- priors
  save(smc_res, file="resSMC-ABC.rda")  
}  

plot.ts(smc_res$epsilon, xlab="SMC iteration", ylab=expression(epsilon[t]))
abline(h=0,lty=2)
dev.copy2eps(file="ABC-SMCepsilon.eps")

smc_res$wt[is.na(smc_res$wt)] <- 0
hist(smc_res$beta, freq=F, col="lightgreen", main="", xlab=expression(beta))
lines(density(smc_res$beta, weights=smc_res$wt), col=4, lty=2, lwd=3)
dev.copy2eps(file="ABC-SMCposterior.eps")

plot.ts(smc_res$ess, xlab="SMC iteration", ylab="ESS")
abline(h=1000, col=4, lty=2)
dev.copy2eps(file="ABC-SMC_ESS.eps")

plot.ts(smc_res$variance, xlab="SMC iteration", ylab=expression(sigma[MH]))
abline(h=0,lty=2)

plot.ts(smc_res$accept/length(smc_res$wt), xlab="SMC iteration", ylab="M-H acceptance",ylim=c(0,1))

colMeans(smc_res$mu)
pred <- smc_res$alloc/rowSums(smc_res$alloc)
predMx <- array(dim=c(100,100,3))
predMx[,,3] <- pred[,4]
predMx[,,1] <- pred[,6]
predMx[,,2] <- (pred[,1] + pred[,2] + pred[,3] + pred[,5])/4
predMx[,,1] <- t(predMx[,,1])
predMx[,,2] <- t(predMx[,,2])
predMx[,,3] <- t(predMx[,,3])
plot(c(0,1),c(0,1),type='n',xaxt='n',yaxt='n',xlab="",ylab="",asp=1)
rasterImage(as.raster(predMx[100:1,,]), 0, 0, 1, 1, interpolate = FALSE)
dev.copy2eps(file="ABC-SMCsegmentation.eps")

plot(c(0,1),c(0,1),type='n',xaxt='n',yaxt='n',xlab="",ylab="",asp=1)
rasterImage(as.raster(t(y[,100:1])/max(y)), 0, 0, 1, 1, interpolate = FALSE)
dev.copy2eps(file="ABC-SMC_Menteith.eps")

2/diff(range(smc_res$beta)) * 2000