inst/doc/mcmcPottsNoData.R

In bayesImageS: Bayesian Methods for Image Segmentation using a Potts Model

## ----setup, include=FALSE-----------------------------------------------------
library(bayesImageS)
set.seed(12)

## ----fig.cap="Chequerboard Gibbs sampling for a regular 8x8 lattice with k=3 labels."----
k <- 3
n <- 8*8
mask <- matrix(1,nrow=sqrt(n),ncol=sqrt(n))
neigh <- getNeighbors(mask, c(2,2,0,0))
block <- getBlocks(mask, 2)
edges <- getEdges(mask, c(2,2,0,0))
print(paste(sum(mask),"pixels"))
print(paste("maximum sufficient statistic S(z) =",nrow(edges)))
res.Gibbs <- mcmcPottsNoData(beta=5, k=3, neigh, block, niter=100)
ts.plot(res.Gibbs$sum, ylim=c(nrow(edges)/3, nrow(edges)))
abline(h=nrow(edges), col=2, lty=3)
summary(res.Gibbs$sum[51:100])
var(res.Gibbs$sum[51:100])

## ----fig.cap="Chequerboard Gibbs sampling for a regular 125x125 lattice with k=3 labels."----
n <- 125*125
mask <- matrix(1,nrow=sqrt(n),ncol=sqrt(n))
neigh <- getNeighbors(mask, c(2,2,0,0))
block <- getBlocks(mask, 2)
edges <- getEdges(mask, c(2,2,0,0))
print(paste(sum(mask),"pixels"))
print(paste("maximum sufficient statistic S(z) =",nrow(edges)))
res.Gibbs <- mcmcPottsNoData(beta=5, k=3, neigh, block, niter=2000)
ts.plot(res.Gibbs$sum, ylim=c(nrow(edges)/3, nrow(edges)))
abline(h=nrow(edges), col=2, lty=3)
summary(res.Gibbs$sum[1001:2000])
var(res.Gibbs$sum[1001:2000])

## -----------------------------------------------------------------------------
res.SW <- swNoData(beta=5, k=3, neigh, block, niter=50)
ts.plot(res.SW$sum, ylim=c(nrow(edges)/3, nrow(edges)))
abline(h=nrow(edges), col=2, lty=3)
summary(res.SW$sum[26:50])
var(res.SW$sum[26:50])

## -----------------------------------------------------------------------------
n <- 1000*1000
mask <- matrix(1,nrow=sqrt(n),ncol=sqrt(n))
neigh <- getNeighbors(mask, c(2,2,0,0))
block <- getBlocks(mask, 2)
edges <- getEdges(mask, c(2,2,0,0))
print(paste(sum(mask),"pixels"))
print(paste("maximum sufficient statistic S(z) =",nrow(edges)))
res.SW <- swNoData(beta=5, k=3, neigh, block, niter=100)
ts.plot(res.SW$sum, ylim=c(nrow(edges)/3, nrow(edges)))
abline(h=nrow(edges), col=2, lty=3)
summary(res.SW$sum[81:100])
var(res.SW$sum[81:100])

## -----------------------------------------------------------------------------
n <- 125*125
mask <- matrix(1,nrow=sqrt(n),ncol=sqrt(n))
neigh <- getNeighbors(mask, c(2,2,0,0))
block <- getBlocks(mask, 2)
edges <- getEdges(mask, c(2,2,0,0))
system.time(res2.Gibbs <- mcmcPottsNoData(beta=0, k=3, neigh, block, niter=100, random=FALSE))
ts.plot(res2.Gibbs$sum, ylim=range(c(res2.Gibbs$sum, nrow(edges))))
abline(h=nrow(edges), col=2, lty=3)
abline(h=nrow(edges)/3, col=4, lty=3)
summary(res2.Gibbs$sum[51:100])
var(res2.Gibbs$sum[51:100])

## -----------------------------------------------------------------------------
muSz <- nrow(edges)/3
sdSz <- sqrt(nrow(edges)*(1/3)*(2/3))
hist(res2.Gibbs$sum[51:100], freq=FALSE, breaks=20, col=3,
     xlim=range(c(res2.Gibbs$sum,muSz - 3*sdSz, muSz + 3*sdSz)))
abline(v=nrow(edges)/3, col=4, lty=3, lwd=3)
curve(dnorm(x, mean=nrow(edges)/3, sd=sqrt(nrow(edges)*(1/3)*(2/3))), 
          col="darkblue", lwd=2, lty=2, add=TRUE, yaxt="n")

## -----------------------------------------------------------------------------
system.time(res2.SW <- swNoData(beta=0, k=3, neigh, block, niter=100, random=FALSE))
ts.plot(res2.SW$sum, ylim=range(c(res2.SW$sum, nrow(edges))))
abline(h=nrow(edges), col=2, lty=3)
abline(h=nrow(edges)/3, col=4, lty=3)
summary(res2.SW$sum[51:100])
var(res2.SW$sum[51:100])

## -----------------------------------------------------------------------------
hist(res2.SW$sum[51:100], freq=FALSE, breaks=20, col=3,
     xlim=range(c(res2.Gibbs$sum,muSz - 3*sdSz, muSz + 3*sdSz)))
abline(v=nrow(edges)/3, col=4, lty=3, lwd=3)
curve(dnorm(x, mean=muSz, sd=sdSz), 
          col="darkblue", lwd=2, lty=2, add=TRUE, yaxt="n")