Analysis and visualization of samples from stan

[Douglas Bates]

What tools are popular for the analysis and visualization of samples created by stan?  In particular, I would be interested in tools that use both lp__ and the parameter values.

The manual mentions calculation of ESS, the effective sample size.  Are there convenient ways of assessing that, say in R?

[Ben Goodrich]

Hi Doug,

On Monday, October 8, 2012 11:32:24 AM UTC-4, Douglas Bates wrote:

What tools are popular for the analysis and visualization of samples created by stan?  In particular, I would be interested in tools that use both lp__ and the parameter values.

The manual mentions calculation of ESS, the effective sample size.  Are there convenient ways of assessing that, say in R?

 
If you use the rstan package, there are print and plot methods for the resulting stanfit object that do such things. As in

stopifnot(require(rstan))

schools_dat <- list(J = 8,
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))

eight_schools <- stan(file = "http://wiki.stan.googlecode.com/git/rstangettingstarted/8schools.stan",
                      data = schools_dat)
print(eight_schools)
plot(eight_schools)

If you use command-line Stan, it just writes a .csv file, but at the moment, I think there is no convenience function to construct a stanfit object from the .csv files.

You can load library(mcmc), make a mcmc object from the .csv file, and and analyze it with their functions with

stopifnot(require(mcmc))
posterior <- as.mcmc(read.csv("samples.csv", comment.char = "#"))

Possibly also the library(boa).

Ben
 

[Bob Carpenter]

Oops -- forgot to send this until I saw Ben's response.

RStan derives from Coda for analysis and R2WinBUGS and R2jags for
graphics.

The main difference is that Coda uses a slightly different
ESS estimator and a slightly different version of the Gelman
and Rubin R-hat convergence statistic. Another minor diff
is that RStan shows the warmup samples as well as the
actual samples in the traceplots.

Another alternative that also provides Coda-like analysis
and R2WinBUGS-like graphing is mcmcplots by Curtis and Goldin.

- Bob

[Andrew Gelman]

In addition to the default summaries and graphics from print() and plot(), rstan also has extractor functions to pull out the simulation draws directly.  So you can extract any parameters of interest and also lp__

A

--
 
 

[Xavier Fernández i Marín]

Douglas Bates vas escriure el dia dl, 08 oct 2012:

> What tools are popular for the analysis and visualization of samples
> created by stan? In particular, I would be interested in tools that use
> both lp__ and the parameter values.

Hi all,

Although this was posted several days ago, I would like to add that there
is a new package in R called "ggmcmc" that is intended to do precisely
that: visual inspection of convergence. Now it only shows the usual plots
(densities, traceplots, running means, comparisons of full chain and last
part, autocorrelations, crosscorrelations, histograms and even
caterpillarplots), but I am planning to extend it to show visually the
formal diagnostics.

To use it in stan:
install.packages("ggmcmc", dependencies=TRUE)

Assume that 'm' is the result of a call to stan(), with the default 4
chains, then:
-----8<---------------
library(ggmcmc)
s <- extract(m, inc_warmup=FALSE)
s <- mcmc.list(mcmc(s[,1,-4]), mcmc(s[,2,-4]), mcmc(s[,3,-4]), mcmc(s[,4,-4]))
S <- ggs(s)
-----8<---------------

And now 'S' is an object ready to be plotted using ggmcmc
(ggs_traceplot(S), ggs_running(S), ...)

Or simply the wrapper function "ggmcmc(S)".


If stan is used without 'rstan', and the samples have been stored in
s_1.csv and s_2.csv:
-----8<---------------
library(coda)
s1 <- read.table("s_1.csv", sep=",", header=TRUE)
s2 <- read.table("s_2.csv", sep=",", header=TRUE)
s.mcmc <- mcmc.list(mcmc(s1[,-c(1:3)]), mcmc(s2[,-c(1:3)]))
S <- ggs(s.mcmc)
-----8<---------------


I have written a detailed example on how to process Stan samples in ggmcmc:
http://blog.xavier-fim.net/2012/10/processing-stan-samples-with-ggmcmc

And there is a general introduction to ggmcmc:
http://xavier-fim.net/packages/ggmcmc/

Hope this helps,

--
- Xavier -

[Sebastian Weber]

Just found a neat way to get posterior samples from stan to coda:

posterior <- extract(se$stanFit, inc_warmup=FALSE)
posterior_coda <- do.call(mcmc.list, alply(posterior, 2, mcmc) )

S <- ggs(posterior_coda)


… but this version does not exclude element 4 in the third dimension which I did not understand. Of course, you need plyr loaded.

> --
>
>

--
Sebastian Weber
Group of Cell Communication and Control
Freiburg Institute for Advanced Studies - FRIAS
School of Life Sciences - LIFENET
Albert-Ludwigs-Universität Freiburg
Albertstr. 19
79104 Freiburg
T.: +49-761-203-97237
Fax:+49-761-203-97334

[Andrew Gelman]

Xavier:
In this case, the following extract from BDA3 might be helpful in setting up your formal convergence diagnostics. We have some new twists that go beyond the earlier R-hat stuff.
Andrew

[Xavier Fernández i Marín]

Sebastian Weber vas escriure el dia dj, 18 oct 2012:

> Just found a neat way to get posterior samples from stan to coda:
>
> posterior <- extract(se$stanFit, inc_warmup=FALSE) posterior_coda <-
> do.call(mcmc.list, alply(posterior, 2, mcmc) )
>
> S <- ggs(posterior_coda)
>
>
> … but this version does not exclude element 4 in the third dimension
> which I did not understand. Of course, you need plyr loaded.

Thank you Sebastian,

Excluding element 4 is trivial:
posterior_coda <- do.call(mcmc.list, alply(posterior[,,-4], 2, mcmc) )

But of course your solution is much cleaner.

I have updated the post:
http://blog.xavier-fim.net/2012/10/processing-stan-samples-with-ggmcmc-2/

Best,

--
- Xavier -

[Xavier Fernández i Marín]

Andrew Gelman vas escriure el dia dj, 18 oct 2012:

> Xavier: In this case, the following extract from BDA3 might be helpful in
> setting up your formal convergence diagnostics. We have some new twists
> that go beyond the earlier R-hat stuff.

Thank you very much Andrew. I have the Geweke test plot almost ready, and
will continue working with the Rhat.

I am planning to do the following four: Geweke, Gelman & Rubin,
Heidelberger-Welch and Raftery-Lewis (see on github
https://github.com/xfim/ggmcmc/issues/1). Am I missing anything.

Will keep you posted when it is ready, so that you may want to announce in
your blog.

[Stéphane Laurent]

Hello Xavier,

 Thanks for this great post on your blog. 

However the output of the extract() function is a list for me and then your code doesn't work, I don't understand why you get an array.

[Stéphane Laurent]

Ah ok I have find how to get the array : one has to use the option permuted=FALSE in extract().