R/GelmanRubin.R

In Scopia/StatComp18053: generating and monitoring convergence of Markov chain and permutation test for equal distribution

#' @title function for monitoring convergence of Markov chain
#' @description This function use the Gelman-Rubin method to monitor convergence of Markov chains.
#' @usage GelmanRubin(psi)
#' @param psi the statistic for chain of X
#' @details user should provide the statistic psi(X)
#' @references Statistic computing with R. Maria L. Rizzo
#' @examples
#' \dontrun{
#' lden=dcauchy
#' sigma=3
#' M4=cumsum((MC4=Metropolis(dist=lden,N=n,sigma=sigma,print_acc=F)))/1:length(MC4)
#' M3=cumsum((MC3=Metropolis(dist=lden,N=n,sigma=sigma,print_acc=F)))/1:length(MC3)
#' M2=cumsum((MC2=Metropolis(dist=lden,N=n,sigma=sigma,print_acc=F)))/1:length(MC2)
#' M1=cumsum((MC1=Metropolis(dist=lden,N=n,sigma=sigma,print_acc=F)))/1:length(MC1)
#' psi=rbind(M1,M2,M3,M4)
#' plot((R=sapply(1:2000,function(i)GelmanRubin(psi[,1:i]))),
#'       main="R value of Gelman-Rubin method",ylim=c(1,2))
#' }
#' @export
GelmanRubin <- function(psi) {
  psi <- as.matrix(psi)
  n <- ncol(psi)
  k <- nrow(psi)
  psi.means <- rowMeans(psi)
  B <- n * var(psi.means)
  psi.w <- apply(psi, 1, "var")
  W <- mean(psi.w)
  v.hat <- W*(n-1)/n + (B/n)
  r.hat <- v.hat / W
  return(r.hat)
}