Implementation in R of algorithms from chapter 3 of Perfect Sampling (M. Huber, 2016, pp. 43-52): coupling from the past; doubling CFTP; and monotonic CFTP for the Ising model.

Huber2016_ch3_part1.R

set.seed(1234)

# Exercise 3.1: simple symmetric random walk on {0,1,2}
upd <- function(x, U) {
  x + ifelse(x < 2 && U > 0.5, 1, 0) - ifelse(x > 0 && U <= 0.5, 1, 0)
}

x <- 1
upd(x,0.25)
upd(x,0.75)
upd(upd(x,0.75),0.25)

# complete coupling: simulating *forwards*
x = 0
y = 1
z = 2
t = 0
while (x != y || y != z) {
  t <- t+1
  u <- runif(1)
  print(paste(t,x,y,z,u))
  x <- upd(x,u)
  y <- upd(y,u)
  z <- upd(z,u)
}
identical(x,y,z)
# we actually need the time-reversal of this chain...

# Coupling_from_the_past (pg. 46)
cftp <- function(t, x0, A, upd) {
  U <- runif(t)
  x <- x0
  print(U)
  if (!all(U > A[1,], U < A[2,])) {
    x <- cftp(t, x0, A, upd)
  }
  print(x)
  for (i in 1:t) {
    x <- upd(x,U[i])
  }
  return(x)
}

# Example 3.2 (pg. 47)
cftp(2, sample(0:2,1), matrix(c(0,0.5,0,0.5),nrow=2), upd)

# Doubling_Coupling_from_the_past (pg. 49)
# NOTE: since t is now variable-length, need to search for a subsequence
# that satisfies A. Therefore, A is limited to an open interval
# on [0,1]^t where the lower bound (and upper bound) is identical for
# every dimension. This limit is imposed for computational convenience,
# rather than any mathematical necessity.
double_cftp <- function(t, x0, A, upd) {
  U <- runif(t)
  x <- x0
  print(t)
  idx <- which(U > A[1] & U < A[2]) # search for subsequence
  if (sum(diff(idx) == 1) == 0) {
    x <- double_cftp(2*t, x0, A, upd)
  }
  for (i in 1:t) {
    x <- upd(x,U[i])
  }
  return(x)
}

double_cftp(1, sample(0:2,1), c(0,0.5), upd)

library(bayesImageS)
Xmax <- matrix(1,nrow=2,ncol=2)
Xmin <- matrix(-1,nrow=2,ncol=2)
Xneigh <- getNeighbors(Xmax, neiStruc=c(2,2,0,0))

# single site update using random scan (pg. 18)
Ising_Gibbs <- function(x, u, idx, beta, n) {
  neigh <- n[idx,]
  n1 <- sum(x[neigh] == 1)
  n2 <- sum(x[neigh] == -1)
  if (log(u) < beta*n1 - log(exp(beta*n1) + exp(beta*n2))) {
    x[idx] <- 1
  } else {
    x[idx] <- -1
  }
  return(x);
}

Ising_Gibbs(c(Xmax,0), 0.9, 2, 0.5, Xneigh)
Ising_Gibbs(c(Xmin,0), 0.9, 2, 0.5, Xneigh)

Ising_Gibbs(c(Xmax,0), 0.02, 2, 0.5, Xneigh)
Ising_Gibbs(c(Xmin,0), 0.02, 2, 0.5, Xneigh)

# Monotonic_coupling_from_the_past (pg. 52)
mono_cftp_Ising <- function(t, x0, y0, beta, neigh) {
  print(t)
  U <- runif(t)
  I <- sample(1:nrow(neigh), t, replace = TRUE)
  for (i in 1:t) {
    x0 <- Ising_Gibbs(x0, U[i], I[i], beta, neigh)
    y0 <- Ising_Gibbs(y0, U[i], I[i], beta, neigh)
  }
  if (identical(x0, y0)) {
    return(x0)
  } else {
    x0 <- mono_cftp_Ising(2*t, x0, y0, beta, neigh)
    # read twice
    for (i in 1:t) {
      x0 <- Ising_Gibbs(x0, U[i], I[i], beta, neigh)
    }
    return(x0)
  }
}

system.time(X <- mono_cftp_Ising(length(Xmax), c(Xmin,0), c(Xmax,0), 0.5, Xneigh))
matrix(X[1:length(Xmin)],nrow=nrow(Xmin)) # drop the dummy state

Huber2016_ch3_part2.R

# simple slice sampler (eg. 8.1 in Robert & Casella, 2003, pg. 323)
f_x <- function(x) {0.5 * exp(- sqrt(x))}
inv_fx <- function(y) {log(2*y)^2}

sliceSamp <- function(niter, x0, f_x, inv_fx) {
  samp <- samp_y <- samp_u <- numeric(niter)
  samp[1] <- x0
  samp_y[1] <- 0
  samp_u[1] <- 0
  for (it in 2:niter) {
    samp_u[it] <- runif(1)
    samp_y[it] <- samp_u[it] * f_x(samp[it-1])
    u2 <- runif(1)
    samp[it] <- u2 * inv_fx(samp_y[it])
  }
  list(x=samp, y=samp_y, u=samp_u)
}

# see Fig. 8.1 on pg. 324
x_0 <- rexp(1)
s_samp <- sliceSamp(50000, x_0, f_x, inv_fx)
hist(s_samp$x,col="lightgreen",breaks=100,freq=FALSE,main="",xlab="x")

# eg. 8.2: truncated Gaussian
f_x <- function(x) {(x >= 0) * exp(-0.5 * (x+3)^2)  * (x <= 1)}
inv_fx <- function(y) {
  x <- sqrt(-2 * log(y)) - 3
  ifelse(x >= 0 && x <= 1, x, 1)
}

# see Fig. 8.2 on pg. 325
x_0 <- runif(1)
nit <- 11
s_samp <- sliceSamp(nit, x_0, f_x, inv_fx)
plot(s_samp$x, s_samp$y, ylim=c(0,f_x(0)), xlim=c(0,1), xlab="x", ylab="y")
curve(exp(-0.5 * (x+3)^2), add=TRUE, col=4, lwd=2, lty=2)
lines(rep(s_samp$x,each=2), c(0, rep(s_samp$y[2:nit], each=2), s_samp$y[nit]))
abline(h=0,lty=2)

# perfect slice sampler (Huber 2016, pg. 57)
x_1 <- 0 # unique maximum value of f_x
x_0 <- 1 # unique minimum

coupledSliceSamp <- function(x_0, x_1, u, u2, f_x, inv_fx) {
  y_0 <- u * f_x(x_0)
  x_0 <- u2 * inv_fx(y_0)

  y_1 <- u * f_x(x_1)
  if (f_x(x_0) >= y_1) {
    x_1 <- x_0
  } else {
    x_1 <- u2 * inv_fx(y_1)
  }
  list(x_0=x_0, y_0=y_0, x_1=x_1, y_1=y_1)
}

coupledSliceSamp(x_0, x_1, runif(1), runif(1), f_x, inv_fx)

perfectSliceSamp <- function(x_0, x_1, t, f_x, inv_fx) {
  u <- runif(t)
  u2 <- runif(t)
  for (i in 1:t) {
    res <- coupledSliceSamp(x_0, x_1, u[i], u2[i], f_x, inv_fx)
    x_0 <- res$x_0
    x_1 <- res$x_1
  }
  if (x_0 == x_1) {
    return(x_1)
  } else {
    print(paste("T = ", log(t)/log(2) + 1, "; t =", 2*t))
    x_1 <- perfectSliceSamp(x_0, x_1, 2*t, f_x, inv_fx)
    # read twice
    for (i in t:1) {
      res <- coupledSliceSamp(x_1, x_1, u[i], u2[i], f_x, inv_fx)
      x_0 <- res$x_0
      x_1 <- res$x_1
    }
    return(x_1)
  }
}

x_1 <- 0 # unique maximum value of f_x
x_0 <- 1 # unique minimum
perfectSliceSamp(x_0, x_1, 1, f_x, inv_fx)

library(truncnorm)
nit <- 1000
sam <- numeric(nit)
for (it in 1:nit) {
  sam[it] <- perfectSliceSamp(1, 0, 1, f_x, inv_fx)
}
hist(sam, freq=FALSE, main="Truncated Gaussian", xlab="x", ylim=c(0,dtruncnorm(0,0,1,-3)),breaks=20)
curve(dtruncnorm(x,a=0,b=1,mean=-3,sd=1), add=TRUE, col=4, lwd=2, lty=2)

sam2 <- sliceSamp(1000, runif(1), f_x, inv_fx)
hist(sam2$x, freq=FALSE, main="Truncated Gaussian", xlab="x", ylim=c(0,dtruncnorm(0,0,1,-3)),breaks=30)
curve(dtruncnorm(x,a=0,b=1,mean=-3,sd=1), add=TRUE, col=4, lwd=2, lty=2)
range(sam2$x)
which(sam2$x > 1)