Extract Sampler Parameters in RStan

[Michael Betancourt]

I’ve said it before and I’ll say it again, not having sampler parameters
(accept stat, n divergent, etc) included in extract by default is dangerous.
The user should be forced to look at these Clockwork Orange-style.

But I also can’t figure out how to access them from a stanfit object at
all (there’s no official doc on how to do this and the archived users
list is no help). Seriously, why is this so hard?

[Bob Carpenter]

The magic key you're looking for is:

get_adaptation_info

It IS in the doc.

> On Mar 2, 2015, at 11:58 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
> I’ve said it before and I’ll say it again, not having sampler parameters
> (accept stat, n divergent, etc) included in extract by default is dangerous.
> The user should be forced to look at these Clockwork Orange-style.

I agree that these should be included in extract, just like lp__.

> But I also can’t figure out how to access them from a stanfit object at
> all (there’s no official doc on how to do this and the archived users
> list is no help). Seriously, why is this so hard?

In the doc for stanfit-class, you see:

get_sampler_params
signature(object = "stanfit"): obtain the parameters used for
the sampler such as stepsizeand treedepth. The results are returned
as a list, each element of which is an array for a chain. The array
has number of columns corresponding to the number of parameters used
in the sampler and its column names provide the parameter names. Optional
parameter inc_warmup indicates whether to include the warmup period.

Here's an example:

> fit <- stan(model_code = "parameters { real y; } model { y ~ normal(0,1); }")

> samp_params <- get_sampler_params(fit);

At this point, samp_params[[1]] is the first chain, samp_params[[2]] is
the second chain and so on. Each chain has everyting in order.

> head(sps[[1]])
accept_stat__ stepsize__ treedepth__ n_leapfrog__ n_divergent__
[1,] 4.781267e-01 2.0000000 1 1 0
[2,] 1.860020e-218 5.5698093 1 1 0
[3,] 9.682160e-01 0.7105644 3 5 0
[4,] 8.290335e-01 0.7877146 2 3 0
[5,] 9.444573e-01 0.7123420 2 3 0
[6,] 9.966285e-01 0.9768503 1 1 0

If you want to get rid of warmup, use

get_sampler_params(fit, inc_warmup=FALSE)

but you probably don't for these since you've said we should be
paying attention to n_divergent__ in warmup.

It's in the R doc, but not under what you may be looking for,

This kind of doc issue (not knowing how the stanfit object was constructed
so not knowing what's in it by design) is one reason I don't like having single
classes or functions overloaded for different types of return types.
There's no way to provide clean doc for the get, because the columns
will vary based on the type of sampler (e.g., NUTS has treedepth__ but
vanilla HMC doesn't). Similarly, stan() is used for the call for both vanilla
HMC and NUTS. Just like in the command line.

- Bob

[Bob Carpenter]

Sorry --- those first two lines are wrong --- it's get_sampler_params,
not get_adaptation_info. I have a hard time remembering which is which,
since adaptation info is related to the samplers.

- 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/d/optout.
>
>

[Michael Betancourt]

How did you get the doc for the stanfit class? I tried
all of the method()-like functions I could google but
the best I could get was

"> methods(class='stanfit')
[1] as.array.stanfit* as.data.frame.stanfit* as.matrix.stanfit* dim.stanfit*
[5] dimnames.stanfit* is.array.stanfit* pairs.stanfit* print.stanfit*

Non-visible functions are asterisked"

And then was does get_sampler_params return? It’s
a list of what, exactly? typeof() claims that it’s a double,
but it’s single-indexed and it’s not at all clear how to
extract a single sampler parameter.

[Bob Carpenter]

I opened the R GUI. Then

> help(rstan)

and that pops up a page the bottom of which has an
index. Click on that, and you get and you get an
index, one element of which is stanfit-class. Click on
that and you get all the doc.

get_sampler_params returns a list the same size as
the number of chains. I find R's types totally messed
up --- lists are like Python dictionaries, so they can
be heterogeneous, but in this case, each element is
itself a matrix (or array or data frame or something --- I
always find it hard to tell since they all have labeled columns).
The column labels have the names. I'm not even sure
how to pull elements out.

I find R very confusing. But a little trial and error and
you get this:

> sps[[1]][1,5]
n_divergent__
0

So that's the first iteration's n_divergent__. No idea
why it has a label. If you do

> sps[[1]][1:1000,5]

you get all the n_divergent__
for chain 1, but for reasons I don't even begin to understand,
the label disappears.

> sps[[1]][1:1000,5]
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[71] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[141] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[211] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[281] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[351] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[421] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[491] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[561] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[631] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[701] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[771] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[841] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[911] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[981] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Good luck.

- Bob

[Ben Goodrich]

On Monday, March 2, 2015 at 9:39:09 AM UTC-5, Michael Betancourt wrote:

How did you get  the doc for the stanfit class?  I tried
all of the method()-like functions I could google but
the best I could get was

"> methods(class='stanfit')
[1] as.array.stanfit*      as.data.frame.stanfit* as.matrix.stanfit*     dim.stanfit*          
[5] dimnames.stanfit*      is.array.stanfit*      pairs.stanfit*         print.stanfit*        

   Non-visible functions are asterisked"

Those aren't the methods you are looking for. More specifically, those are the S3 methods. You want the S4 methods, but you can't get them by doing

showMethods("stanfit-class")

because S4 methods are defined for generic functions with S4 classes as signatures instead of being part of the classes definition (like every other object-oriented language and the S4.5 classes in R, a.k.a. reference classes, which are not currently used by RStan but I am in favor of converting to). The shortest path to what you want probably would have been to look at all the functions exported by the rstan package and then guess what they do.

> ls("package:rstan")
 [1] "constrain_pars"      "extract"             "get_adaptation_info"
 [4] "get_cppcode"         "get_cppo_mode"       "get_cxxflags"
 [7] "get_inits"           "get_logposterior"    "get_num_upars"
[10] "get_posterior_mean"  "get_sampler_params"  "get_seed"
[13] "get_seeds"           "get_stancode"        "get_stanmodel"
[16] "grad_log_prob"       "log_prob"            "lookup"
[19] "makeconf_path"       "monitor"             "optimizing"
[22] "plot"                "read_rdump"          "read_stan_csv"
[25] "rstan_options"       "sampling"            "sflist2stanfit"
[28] "show"                "stan"                "stanc"
[31] "stan_model"          "stan_rdump"          "stan_version"
[34] "summary"             "testify"             "traceplot"
[37] "unconstrain_pars"

And then was does get_sampler_params return?  It’s
a list of what, exactly?  typeof() claims that it’s a double,
but it’s single-indexed and it’s not at all clear how to
extract a single sampler parameter.

It is a list whose length is equal to the number of chains. Each list element is an iterations x 4 matrix. It works almost like a std::vector except that you have to index it with double brackets. So,

my_list <- get_sampler_params(posterior)
my_list[[1]][2,3]

extracts the first list element and then the element in the second row and third column of that matrix. Typically, you would want to apply the same function to each chain, in which case it would be something like

sapply(my_list, FUN = colMeans)

Ben

[Krzysztof Sakrejda]

On Monday, March 2, 2015 at 10:17:35 AM UTC-5, Ben Goodrich wrote:

On Monday, March 2, 2015 at 9:39:09 AM UTC-5, Michael Betancourt wrote:

Those aren't the methods you are looking for. More specifically, those are the S3 methods. You want the S4 methods, but you can't get them by doing

[snip]

With condolences to Michael, this is by far the most entertaining response I've seen on the whole list, now we just need to record Obi Wan saying it.  It's like the
"ya can't get there from here" joke for computational statistics, and it's funny for the same reason. 

Krzysztof
 

[Krzysztof Sakrejda]

btw....

On Monday, March 2, 2015 at 9:51:07 AM UTC-5, Bob Carpenter wrote:

the number of chains.  I find R's types totally messed
up --- lists are like Python dictionaries, so they can
be heterogeneous, but in this case, each element is
itself a matrix (or array or data frame or something --- I
always find it hard to tell since they all have labeled columns).  
The column labels have the names.  I'm not even sure
how to pull elements out.

 
Unless you know what a function will do already, the only way to find out are the 'str' and 'class' functions.

I find R very confusing.  But a little trial and error and
you get this:

> sps[[1]][1,5]
n_divergent__
            0

So that's the first iteration's n_divergent__.  No idea
why it has a label.  If you do
 
> sps[[1]][1:1000,5]

you get all the n_divergent__
for chain 1, but for reasons I don't even begin to understand,
the label disappears.

The "reason" is that R helpfully does helpful conversions for you (like dropping the 'dimnames' attribute), and the bracket operator only behaves like something reasonable if you add "drop=FALSE" as an argument:

> str(sps[[1]][1:10,5] )
 num [1:10] 0 1 0 0 0 0 0 0 0 0
> str(sps[[1]][1:10,5,drop=FALSE] )
 num [1:10, 1] 0 1 0 0 0 0 0 0 0 0
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr "n_divergent__"

 

[Jonah]

A summary table of the min, max, and avg values of all the sampler parameters is also included automatically in shinyStan under the "Diagnose" tab. And it's probably much quicker to just  do launch_shinystan(stanfit) than try to figure out how to use get_sampler_params and then deal with the list object it gives you. In fact the entire motivation for making shinyStan was so I would never have to type 'get_sampler_params' ever again. But now I've gone and done it twice in one message. 

[Jonah]

Here's a function that will return the max treedepth for each chain. If you change "treedepth__" to any of the other sampler params and change "max" whatever function you want then you can pretty much get any sampler param info you want using this template. 

max_treedepth <- function(stanfit, inc_warmup = TRUE) {

  max_td <- sapply(get_sampler_params(stanfit, inc_warmup = inc_warmup),

                   function(x) max(x[,"treedepth__"]))

  names(max_td) <- paste0("chain", 1:length(max_td))

  max_td

}

If you don't feel like writing each of those functions, I've already done most of them for you and whole bunch of others. There's a handy little convenience function in the shinyStan package called stan_get that you can call from R.  See help("stan_get", "shinyStan"). 

[Michael Betancourt]

You really want to look at the entire distribution of each diagnostic
variable aggregated across the chains, so in the end I had to do

sampler_params <- get_sampler_params(fit, inc_warmup=FALSE)
sapply(sampler_params, function(x) c(x[,'treedepth__']))[,1]

In PyStan it’s a similar problem but reduce makes it slightly cleaner,

sampler_params = fit.get_sampler_params(inc_warmup=False)
reduce(lambda x, y: x + y['n_divergent__'].tolist(), sampler_params, [])

[Jonah]

On Wednesday, March 4, 2015 at 3:48:07 AM UTC-5, Michael Betancourt wrote:

You really want to look at the entire distribution of each diagnostic
variable aggregated across the chains, so in the end I had to do

sampler_params <- get_sampler_params(fit, inc_warmup=FALSE)
sapply(sampler_params, function(x) c(x[,'treedepth__']))[,1]

Is there a reason for the [,1] index? That should just limit the return value to the treedepths for the first chain. If you want to aggregate then 

c(sapply(sampler_params, function(x) (x[,'treedepth__']))) 

should do the trick, unless I misunderstood what you meant. 

[Michael Betancourt]

sapply just returned one of those weird R arrays with zero

length in dimension and the closing [,1] cleaned it up.  I’m

sure concatenation would also be fine.

[Steve Bellan]

another option:

sampler_params <- get_sampler_params(fit)
rbindlist(lapply(sampler_params, function(x) data.table(x)))

is a clean way to aggregate across chains