[R] Computing ergodic mean with CODA

[Ben Bolker]

Raquel Rangel de Meireles Guimarães <raquelrguimadem <at> gmail.com> writes:

> 
> Hi all,
> 
> I would like to compute ergodic mean using MCMC output from WinBUGS. I 
> tried using CODA package, but it seems that it is not implemented yet.
> 
> Could anyone help me to compute this? Attached to this email are my 
> output and index files.
> 
> Kind regards,
> 
> Raquel
> 
> --

[the data were not in the most convenient form; 
one of the preferences is for a *small* reproducible
example ...  I cut & pasted them into a file to
recreate the original mcmc object].

  Based on a quick look at [Ntzoufras, Ioannis. 2009. Bayesian modeling using
WinBUGS. John Wiley and Sons.], I think cumsum(x)/(1:length(x)) is the
ergodic mean (essentially a running cumulative mean).


x <- read.table("~/R/misc/ergod.dat")
library(coda)
v <- as.mcmc(matrix(x[,2],nrow=5000))
em <- sweep(apply(v,2,cumsum),1,(1:nrow(v)),"/")
library(reshape)
m <- melt(em)
xyplot(value~X1|X2,type="l",data=m,
       scales=list(y=list(relation="free")))