Empty stanfit Object For Computing MCMC Summaries

Interfaces RStan
1

I’m hoping to use stan_ess to compute the effective sample size of a Markov chain created externally to RStan. What’s the best practice for creating an empty stanfit object and then populating the samples with an external array in order to call the official Stan effective sample size estimator? Ideally this wouldn’t involve having to compile any intermediate Stan programs…

1 Like
2

In the R-hat paper I used rstan::monitor (which just takes a 3 dimensional array with dims specified in the help) and extracted the appropriate thing, so I imagine stan_ess works the same.

Edit:
File is here: https://github.com/avehtari/rhat_ess/blob/master/code/simple_example.R

Relevant bit of code:

for(chain in 1:nchains) {
    for (par in 1:2) {
      sims[,chain,par] = arima.sim(n=niter, list(ar=ar_param),sd=sqrt(1-ar_param^2))*sd[chain,par]
    }
    sims[,chain,3] = (arima.sim(n=niter, list(ar=ar_param),sd=sqrt(1-ar_param^2)))/arima.sim(n=niter, list(ar=ar_param),sd=sqrt(1-ar_param^2)) + means[chain]
    sims[,chain,4] = arima.sim(n=niter, list(ar=ar_param),sd=sqrt(1-ar_param^2))/arima.sim(n=niter, list(ar=ar_param),sd=sqrt(1-ar_param^2))
  }
  index_start =  2*npars*(rep-1) + 1;
  index_end = 2*npars*rep;
  nothing <- capture.output(monitor_old <- rstan::monitor(sims,warmup=0))
  nothing <- capture.output(monitor_new <- monitor(sims,warmup=0))
  par_names <- c("Finite mean, different variance", "Finite mean, same variance", "Infinite mean, different location","Infinite mean, same location")
  
  Rhat$rhat[index_start:index_end] = c(monitor_old[,10],monitor_new$Rhat)
  Rhat$version[index_start:index_end] = c(rep("old",npars),rep("new",npars))
  Rhat$par[index_start:index_end] = c(par_names,par_names)
  Rhat$neff[index_start:index_end] = c(monitor_old[,9],monitor_new$Bulk_ESS)
  Rhat$neff_bulk[index_start:index_end] = c(rep(NA,npars),monitor_new$Bulk_ESS)
  Rhat$neff_tail[index_start:index_end] = c(rep(NA,npars), monitor_new$Tail_ESS)
3

Save to csv and populate header to looklike cmdstan output? Then read with rstan (or what tools there are on R side)

4

This works well, thanks.

5

Is there any functionality like monitor in PyStan?

6

PyStan functions only work against fit.

ArviZ works with any mcmc result. (e.g az.from_dict is quite flexible). Rank based methods PR is not yet merged.

7

Unfortunately arviZ is much to bloated for my use case here.