#' @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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.