Is there function to convert from Stan to Coda?

[Suyin Chang]

Hi all, I looked through the past messages in the group looking for good ways to convert from Stan to Coda in R, and did found some one suggesting:

#s <- extract(model2.stanfit, permuted=FALSE,inc_warmup=FALSE)
#s <- do.call(mcmc.list, alply(s[,,], 2, mcmc)) #not excluding l_ph
#s <- do.call(mcmc.list, alply(s[,,-(length(s[1,1,]))], 2, mcmc)) #excluding l_ph

With a few changes it 's been working for me, but I'd appreaciate to know how other people do this and there is a finnished package-function to do such convertion?

Thanks for the help

Chang

[Ben Goodrich]

If you don't need the warmup, you can do

library(coda)
s <- mcmc.list(lapply(1:ncol(fit), function(x) mcmc(as.array(fit)[,x,])))

But, we're always interested to know what people are using coda for that isn't already available in rstan?

Ben

[Suyin Chang]

Hi, thank you Ben, this is a really efficient solution.

About functions in coda, I am interested in doing geweke test,calculating effective sample sizes and ploting running means. Is it possible to do such tests with rstan directly? I didn't find any mention to them through the the rstan functions listed in R help of rstan. Specially geweke test (which I do not like, I think it to be not accurate, but I was asked to use anyway).

[Suyin Chang]

Oh, I just wrote it wrong: where I wrote "effective sample sizes" (which of course stan gives directly), I meant cross-corrlation between variables and its plot.
But Geweke and running means are more important to me anyway.

[Bob Carpenter]

I think it might be nice to have an option to visualize the running
average of parameter samples over the traceplot. It's on the same scale.
I'd restrict it to after warmup, because then the value represents the
estimated mean if sampling had stopped at that iteration.

I don't have any opinion about the Geweke tests.
On the one hand, it would be nice to have Stan satisfy all of
our user's requirements --- if we don't have these things, people will
still need to use Coda, etc., if they need them. On the other
hand, we want to encourage what we think are good practices.

- Bob

[Stéphane Laurent]

Hello,

What is fit in your code ? For me the output of the extract() function is a list, hence this code doesn't work.

[Bob Carpenter]

I created an RStan feature request:

https://github.com/stan-dev/rstan/issues/52

We would love to get some more help on implementing these
things in RStan. If you have a GitHub account, you should
be able to comment on the feature request. If you don't, you
can comment here.

What I want is a way to map a Coda object (or whatever BUGS
and JAGS and others output) to a (reduced) Stan fit object (or to
a new object that the Stan fit object holds onto).

And then all the features that are in Coda. This should stop
all these messages in the future and also help us run the
evaluations vs. BUGS and JAGS that Andrew's been calling for.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.

[Stéphane Laurent]

I have found how to do: https://groups.google.com/forum/#!topic/stan-users/5BZwhLM5OfE (use permuted=FALSE)

[Jeffrey Arnold]

If you mean the R objects that coda uses (as opposed to the text file output), use as.mcmc.list

On Sun, Feb 23, 2014 at 11:55 AM, Stéphane Laurent <lauren...@yahoo.fr> wrote:

I have found how to do: https://groups.google.com/forum/#!topic/stan-users/5BZwhLM5OfE (use permuted=FALSE)

[Bob Carpenter]

That's good news for people who want to use Coda.

I'm more interested in going the other way around and
letting us analyze BUGS and JAGS output using RStan's posterior
analysis tools so we can do the comparisons Andrew wants to do
using the new split R-hat and n_eff calcs. Is that also easy?

- Bob

[Ben Goodrich]

On Sunday, February 23, 2014 12:53:29 PM UTC-5, Bob Carpenter wrote:

That's good news for people who want to use Coda.

I'm more interested in going the other way around and
letting us analyze BUGS and JAGS output using RStan's posterior
analysis tools so we can do the comparisons Andrew wants to do
using the new split R-hat and n_eff calcs.  Is that also easy?

 
The monitor() function already exists in RStan. You just need an array whose dimensions are

simulations x chains x parameters

but coda (and the rest of the world) is inclined to think in terms of

simulations x parameters x chains

so you need

mcmc.list2array <-
  function(x) {
    stopifnot(class(x)[1] != "mcmc.list")
    arr <- array(NA_real_,
                 dim = c(nrow(x[[1]]), length(x), ncol(x[[1]])))
    dimnames(arr)[[3]] <- colnames(x[[1]])
    for(i in seq_along(x)) arr[,i,] <- as.matrix(x[[i]])
    return(arr)
  }

Ben