# file "Rdiag2" -- simple convergence diagnostic for fast-mixing chain

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

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

M = 500  # run length
B = 50  # amount of burn-in
sigma = 10  # proposal scaling (for fast mixing)

# variables for convergence diagnostic
numruns = 50
singlevar = rep(0,numruns)
B1vals = NULL

for (therun in 1:numruns) {

	X = runif(1,0,20)  # initial Markov chain value, chosen Uniform[0,20]

	hlist = NULL  # for keeping track of h function values

	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 values, after burn-in
	    if (i == B+1)
		B1vals = c(B1vals, h(X))
	}

	# update the convergence diagnostic values
	singlevar[therun] = var(hlist)

}

# do the convergence diagnostic analysis
WITH = mean(singlevar)
INTER = var(B1vals)
cat("WITH =", WITH, "... INTER =", INTER, "... ratio =", INTER/WITH, "\n");