# file "Rgel2" -- Gelman-Rubin convergence diagnostic for stationary chain

g = function(x) { (1/2) * dnorm(x,0) + (1/2) * dnorm(x,10) }

h = function(x) { return( x ) }

M = 5000  # run length
B = 50  # amount of burn-in
sigma = 10  # proposal scaling

# variables for convergence diagnostic
numruns = 10
singlemean = rep(0,numruns)
singlevar = rep(0,numruns)

varfact <- function(xxx) { 2 * sum(acf(xxx, plot=FALSE)$acf) - 1 }

for (therun in 1:numruns) {

	# choose initial X from stationary:
	X = rnorm(1)
	if (runif(1) < 0.5)
	    X = X + 10

	hlist = NULL  # for keeping track of h function values
	varfactlist = NULL

	for (i in 1:M) {
	    Y = X + sigma * rnorm(1)  # proposal value (dim=1)
	    U = runif(1)  # for accept/reject
	    alpha = g(Y) / g(X)  # for accept/reject
	    if (U < alpha)
		X = Y  # accept proposal
	    if (i>B)
		hlist = c(hlist, h(X))  # keep h value if burn-in over
	}

	# update the convergence diagnostic values
	singlemean[therun] = mean(hlist)
	singlevar[therun] = var(hlist)
	varfactlist = c(varfactlist, varfact(hlist))

}

# do the convergence diagnostic analysis
withinvar = mean(singlevar)
betweenvar = (M-B) * var(singlemean)
cat("withinvar =", withinvar, "    betweenvar =", betweenvar, "\n");
cat("ratio =", betweenvar/withinvar, "  varfact =", mean(varfactlist),
  "  bound =", 1 + 0.44*(M-B), "\n");