Solved – Using R2WinBUGS, how to extract information from each chain

bugsmarkov-chain-montecarlor

I am new to Gibbs Sampling, and have been using WinBUGS, but I find that it is not well-suited towards storing/presenting results, so I have been calling it from R using the R2WinBUGS package. The data is apparently stored as a "bug" class.

I converted it to coda to run diagnostics, and it displays each of the chains, but I am confused as to the $ extension for each individual chain. I cannot find any good documentation for the "coda" class (the cran instructions are not helpful).

My code is below:

> bugs.sim <- bugs(data, inits, parameter, "gl.bug", n.chains = 5,
                   codaPkg = TRUE, DIC = FALSE, n.iter = 5000)
> codaobject <- read.bugs(bugs.sim)

As you can see, I have 5 chains, and I would like to take the mean and standard deviation of each. How do I go about doing this? I can use the codaobject to take the Geweke diagnostic of each chain, it displays each chain as [[i]] ($i=1,\dots,5$).

Thanks in advance. And any references to a detailed documentation for R2Winbugs, would also be greatly appreciated.

Best Answer

The object returned by read.bugs is an object of S3 class mcmc.list. You can use the double brackets [[ to access the separate chains, i.e. the different mcmc-objects that make up the larger mcmc.list object, which really is simply a list of mcmc-objects that inherits some information about thinning and chain length from its components.

More to the point, s.th. like lapply(codaobject, function(x){ colMeans(x) }) should return the posterior means for each parameter in each chain and lapply(codaobject, function(x){ apply(x, 2, sd) }) should give chain- and parameter-specific posterior sd's, since each chain is essentially just a numeric matrix with rows corresponding to the (saved) iterations and columns corresponding to the different params.

EDIT: I think Gelman/Hill's "Bayesian Data Analysis" contains some worked examples using R2WinBUGS.

Related Solutions

Solved – How to use draws from a Markov chain (Monte Carlo) for (bayesian) inference

Bayesian Statistics does not depend on Markov chains (well in theory), Markov chain Monte Carlo is a method for making the computations in Bayesian Statistics easier (doable). Generally we want to generate data from the posterior distribution which we can easily compute parts of, but not always all of (the normalizing coefficient is usually the hard part), but McMC methods will generate data from the posterior without needing to know some of the parts that are harder to find, that is the benifit of McMC with Bayesian stats (though there are some problems where the posterior is known exactly and McMC is not needed at all). It is very rare that the posterior is uniform, so your example in not realistic, but can still help with understanding.

The short answer to length of burn-in and how long to run the chain is "We Don't Know". The little longer (and more useful) answer is to run the chain for a while, then look at it to see if it looks good (yes that is subjective, but experiance helps, and if we knew the exact answer we would know enough not to need McMC), then often run longer to be sure. I have seen cases where the chain was run for quite a few iterations, it looked good, but the researchers decided to run it for 4 times longer to be sure, and near the end sudenly it jumped to a different part of the distribution that had not been covered before.

More dependency does mean that the chain needs to be run for longer (both burn-in and total runs), but as long as the chain has run long enough the dependency no longer matters. Think of a simple random sample with all values independent, now sort the values into order statistics, they are no longer independent, but the mean and sd are the same and estimate the population just as well. The fact that McMC iterations are dependent does not matter so long as there are enough iterations to fully represent the distribution. Your uniform example shows the risks from running for too short of a time.

Solved – ARIMA: extract date/time information from ARIMA model

Did you see the help file for the ts function? At the bottom there are examples for quarterly and monthly time series. Perhaps you could directly extend this to second-wise series? You only need to specify the start date, not the end date.

For example, instead of using start = c(1954, 7), frequency = 12 for monthly data starting in July 1954 you could use start = c(d, s), frequency = 60*24 for second-wise data starting on day d, second s. The time attribute of your data would be given in days.

set.seed(1)
x=rnorm(10^4) # generate random data
d=1
s=1
x=ts(x,start=c(d,s),frequency=60*24)
plot(x)
print(time(x)) # extract the 'time' attribute of the 'ts' object

There might be a better solution, but perhaps this one is sufficient?

Best Answer

Related Solutions

Solved – How to use draws from a Markov chain (Monte Carlo) for (bayesian) inference

Solved – ARIMA: extract date/time information from ARIMA model

Related Question