R help - Oct 2012 - Coda, HPDinterval and multiple chains

Dear all,

I'm not 100% sure if this question is best directed at the r-list, or a
mailing list concerned with Bayesian analysis, so please accept my apologies if
another audience may be more appropriate.

I have been using the rjags package to run Jags models with multiple chains and
store the results in a Coda based mcmc list.  For instance, having created a
jags model and done initial adapting and updating, I run the following command:

coda_odp_gini_only <- coda.samples(odp_gini_only,
                                  
variable.names=c("beta","sigma2.u2", "deviance"),
                                   n.iter=itercount, thin=thincount)

This create an object with 4 separate chains as requested in my initial call to
Jags which created the object "odp_gini_only".  I then use the  Coda
package to look at the results stored in "coda_odp_gini_only".  Prior
to running any analysis with Coda, I used the command,

coda.options(combine.plots=TRUE, combine.stats=TRUE)

to ask for results that combine the four separate chains.  Sure enough, if I
enter "summary(coda_odp_gini_only)", I am given a single set of output
combining the four chains.  However, if I enter
"HPDinterval(coda_odp_gini_only)" I receive 4 sets of HPD figures, one
for each chain.  Is it possible to combine the four chains together to receive a
single set of HPD estimates?

In a similar vein, is it possible to use the Coda object to estimate the
proportion of a given parameters distribution which is above (or below) a given
value, for instance, the proportion of the distribution of beta[1] greater than
zero?  Again, in doing so, is it possible to combine the results of the four
chains into a single estimate?

Kind regards, and many thanks in advance for any advice anyone can offer me.

Paul
Dr Paul Norris
Lecturer in Social Policy
University of Edinburgh
-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.

How about

require(coda)
data(line)
str(line)
cline <- as.mcmc(do.call(rbind,line))
str(cline)

# Thus
HPDinterval(cline) # (or any FUN.mcmc)
sum(cline[,'alpha'] >= median(cline[,'alpha']))/nrow(cline)

Cheers

On Thu, Oct 4, 2012 at 12:22 PM, NORRIS Paul <P.Norris@ed.ac.uk> wrote:
> Dear all,
>
> I'm not 100% sure if this question is best directed at the r-list, or a
> mailing list concerned with Bayesian analysis, so please accept my
> apologies if another audience may be more appropriate.
>
> I have been using the rjags package to run Jags models with multiple
> chains and store the results in a Coda based mcmc list.  For instance,
> having created a jags model and done initial adapting and updating, I run
> the following command:
>
> coda_odp_gini_only <- coda.samples(odp_gini_only,
>                                   
variable.names=c("beta","sigma2.u2",
> "deviance"),
>                                    n.iter=itercount, thin=thincount)
>
> This create an object with 4 separate chains as requested in my initial
> call to Jags which created the object "odp_gini_only".  I then
use the
>  Coda package to look at the results stored in
"coda_odp_gini_only".  Prior
> to running any analysis with Coda, I used the command,
>
> coda.options(combine.plots=TRUE, combine.stats=TRUE)
>
> to ask for results that combine the four separate chains.  Sure enough, if
> I enter "summary(coda_odp_gini_only)", I am given a single set of
output
> combining the four chains.  However, if I enter
> "HPDinterval(coda_odp_gini_only)" I receive 4 sets of HPD
figures, one for
> each chain.  Is it possible to combine the four chains together to receive
> a single set of HPD estimates?
>
> In a similar vein, is it possible to use the Coda object to estimate the
> proportion of a given parameters distribution which is above (or below) a
> given value, for instance, the proportion of the distribution of beta[1]
> greater than zero?  Again, in doing so, is it possible to combine the
> results of the four chains into a single estimate?
>
> Kind regards, and many thanks in advance for any advice anyone can offer
> me.
>
> Paul
> Dr Paul Norris
> Lecturer in Social Policy
> University of Edinburgh
> --
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
	[[alternative HTML version deleted]]