# file "Rdiag" -- very simple convergence diagnostic

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

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

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

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

for (therun in 1:numruns) {

	X = runif(1,-10,20)  # initial Markov chain value ~ Uniform[-10,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
	}

	# update the convergence diagnostic values
	singlemean[therun] = mean(hlist)

}

# do the convergence diagnostic analysis
ESTSD = sd(singlemean)
print(singlemean)
cat("estimator standard error (ESTSD) =", ESTSD, "\n")