# file "Rgel4" -- convergence diagnostics for slow-mixing 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 = 1 # proposal scaling (for slow mixing) # variables for convergence diagnostic numruns = 10 singlemean = rep(0,numruns) singlevar = rep(0,numruns) 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 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");