# 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");