R/MurraySchemePotts.R

In ick003/SpTMixture: Mixture models with spatio-temporal dependencies

#' Simulate the parameters of a Potts model using the Murray Scheme.
#'
#' This function allows you to simulate the parameters of a Potts model using the Murray Scheme.
#' @param beta.init Starting parameter values for the parameters of the Potts model.
#' @param range.beta Range of values for the parameters of the Potts model.
#' @param n.colors number of different "spins".
#' @param Nb.list List of neighbours.
#' @param col.obs Binary matrix of colors.
#' @param N.run Number of MCMC runs to perform.
#' @keywords MCMC Murray Scheme
#' @references Murray, I. and Ghahramani, Z. and MacKay, D. (2006) MCMC for doubly-intractable distributions. \emph{UAI Proceedings} \bold{22} 359--366.
#' Everitt, R. (2012) Bayesian Parameter Estimation for Latent Markov Random Fields and Social Networks. \emph{Journal of Computational and Graphical Statistics} \bold{21} 940--960.
#' @export
#' @examples
#' beta = c(1.2, 1.2,1.2,1.2)
#' Obs.loc = expand.grid(1:5,1:5)
#' id.pos = sample(1:4, nrow(Obs.loc), replace = TRUE)
#' Res = SwendsenWangAlg(beta, Coord = as.matrix(Obs.loc), col.obs = id.pos, Nrun = 20)
#' par(mfrow = c(1,1))
#' col = grey.colors(4)
#' image(matrix(Res$col.new, ncol =5), col = col)
#' beta.init = runif(1,0,1)
#' set.seed(1)
#' beta.mcmc = MurraySchemePotts(beta.init = beta.init,Coords = Obs.loc ,n.colors = 4, col.obs = Res$col.new, N.run = 5000, range.beta = c(0,3))
#' Nbs = getBonds(as.matrix(Obs.loc), NN = 4, th = 2)
#' beta.mcmc.nb = MurraySchemePotts(beta.init = beta.init,Nb = Nbs,n.colors = 4, col.obs = Res$col.new, N.run = 5000, range.beta = c(0,3))
#'
#' print(beta.mcmc$rate)
#' par(mfrow = c(2,2))
#' hist(beta.mcmc$beta.sample, breaks = 50, xlim = c(0,1.75))
#' plot(beta.mcmc$beta.sample, type = "l", ylim = c(0,1.75))
#' hist(beta.mcmc.nb$beta.sample, breaks = 50, xlim = c(0,1.75))
#' plot(beta.mcmc.nb$beta.sample, type = "l", ylim = c(0,1.75))
MurraySchemePotts = function(beta.init, range.beta = c(0, 1), n.colors=2, Coords = NULL, Nb = NULL, col.obs, N.run=100, print = FALSE){

  beta.hist = matrix(0, nrow = N.run, ncol = length(beta.init))
  beta.hist[1,] = beta.init

    zz = as.matrix(col.obs)

    n.obs = length(zz)

  rate = 0

  beta = as.matrix(beta.init)

  if(!is.null(Coords)){Coords = as.matrix(Coords)}

  if(is.null(Nb)){bds = getBonds(Coords, NN = 4, th = 2)
  }else{bds = Nb}

  col.sum = getCanStat(bds, zz, n.colors)

  for(t in 2:N.run){

    if(print == TRUE){if(t/(N.run/10) == round(t/(N.run/10))){print(t)}}

    # Simulate beta.p ~ Q(.|beta)

    if(length(beta)>1){
      beta.p = runif(length(beta), sapply(beta, function(x) max(9*x/10,range.beta[1])),
                   sapply(beta, function(x) min(11*x/10,range.beta[2])))
    } else {
      beta.p = runif(1, max(9*beta/10,range.beta[1]),min(11*beta/10,range.beta[2]))
    }

    # Simulating z.p ~ f(.|beta.p) using Swendsen-Wang Alg. Need to reach stationary distribution

    beta.temp = rep(beta.p, n.colors)

    res = SwendsenWangAlg(beta.temp, Neigh_Bonds = bds, col.obs = zz, Nrun = 20, sampling = "uniform")

    z.p = res$col.new


    # Calculate probability of acceptance


    col.sum.p = getCanStat(bds, z.p, n.colors)

    p.num = sapply(11/10*beta.p, function(x) min(x,range.beta[2])) - 9 * beta.p/10

    p.denom = sapply(11/10*beta, function(x) min(x,range.beta[2])) - 9 * beta/10

    p.move =  exp(sum((beta.p - beta) * (col.sum - col.sum.p))) * prod(p.num / p.denom)

    p.move = min(p.move, 1, na.rm=T)

    # Exchange with probability p.move

    if(runif(1) < p.move){
      beta = beta.p
      rate = rate + 1
    }

    beta.hist[t,] = beta

  }

  RET = list(beta.old = beta.init, beta.new = beta, beta.sample = beta.hist,
             rate = rate)
  return(RET)
}