one more thing for the Rstan 3 specification

[Andrew Gelman]

In my specs, I don’t think I had much on R-hat, n_eff, and these related post-processing things. Along with these, we’ll also want waic and cross-validation (which will only operate if certain quantities are defined as generated quantities in the model, as discussed in the Waic-Stan paper: http://www.stat.columbia.edu/~gelman/research/unpublished/waic_stan.pdf

Anyway, I just wanted to mention this so we don’t forget it when it becomes time to implement all this.
A

[Bob Carpenter]

Nothing's going to get implemented without a design.
So that's the first step and the question to you is how you
want these interfaces to look from a user's perspective.

- Bob

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

[Andrew Gelman]

OK, here’s a shot. It’s two functions; see below. I also welcome the suggestions of Ben, Jiqiang, and others.
A

P.S. I saw the wiki page and I would’ve put my suggestion there but I couldn’t figure out how to edit it.

foo.summary <- stan_summary(foo.fit.multiple): returns the matrix of mean, se_mean, sd, quantiles, n_eff, R_hat (what’s currently returned in the stan print statement)

foo.waic <- stan_waic_xval(foo.fit.multiple): returns a list of waic, elpd_waic, p_waic, elpd_loo, p_loo, pointwise, total, and se. This function only works if one of the returned parameters in foo.fit.multiple is something called log_lik. If there is no log_lik, it returns NULL. We could start with the function on page 14-15 of this paper: http://www.stat.columbia.edu/~gelman/research/unpublished/waic_stan.pdf

[Bob Carpenter]

On Jun 27, 2014, at 6:03 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

> OK, here's a shot.

That's a sketch of a spec, but it's a step in the right direction
for making things precise.

> It's two functions; see below. I also welcome the suggestions of Ben, Jiqiang, and others.
> A
>
> P.S. I saw the wiki page and I would've put my suggestion there but I couldn't figure out how to edit it.
>

Click on the "edit" button. You probably need to
log in to GitHub.

> foo.summary <- stan_summary(foo.fit.multiple): returns the matrix of mean, se_mean, sd, quantiles, n_eff, R_hat

What do you mean by "matrix of"? Are there row and column labels?
Don't you like structured objects here? So if "Sigma" is an M x N
parameter, how do I access the summary stats for Sigma[m,n]?

> (what's currently returned in the stan print statement)

The way it works now is that stan() returns a fit object, and
that defines the behavior of the generic function print(), which
you won't find direct doc as a function for in RStan because it's
a method on the fit object.

> foo.waic <- stan_waic_xval(foo.fit.multiple): returns a list of waic, elpd_waic, p_waic, elpd_loo, p_loo, pointwise, total, and se. This function only works if one of the returned parameters in foo.fit.multiple is something called log_lik. If there is no log_lik, it returns NULL. We could start with the function on page 14-15 of this paper: http://www.stat.columbia.edu/~gelman/research/unpublished/waic_stan.pdf
>

So you want a list rather than an object so you're OK with users
getting p_loo using the syntax foo.waic$p_loo or foo.waic[["p_loo"]]?

Unless p_loo corresponds to a leave-one-out cross-validation, this
seems like a confusing name.

We'll need doc for what each of these elements is. What I really want
in all of our RStan and PyStan doc is:

argument types: what kind of R object do you pass in?

argument content: what is the argument?

argument constraints: what constraints on the arguments passed in
are there, either singly or in conjunction.

return type: what type of object does it return
what is the content of that object field by field

Do you want to give them standalone names? Names like "total" seem
a bit general.

And speaking of names, looks like you're torn between using something
like "_summarize()" that's a verb and "_summary()" which is the current
function name in R. I'm OK with whatever we decide on naming conventions,
but we should probably stick to it. It'd be nice if we could go with
"summarize" and "optimize" or "summary" and "optimimum" --- either describe
what it's doing or what it returns.

- Bob

[Andrew Gelman]

Hi, see responses below.

On Jun 27, 2014, at 8:58 PM, Bob Carpenter <ca...@alias-i.com> wrote:

> Click on the "edit" button. You probably need to
> log in to GitHub.

I did log into Github but I could not find an edit button.

>> foo.summary <- stan_summary(foo.fit.multiple): returns the matrix of mean, se_mean, sd, quantiles, n_eff, R_hat
>
> What do you mean by "matrix of"? Are there row and column labels?

It’s the matrix that is output here, for example:

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu 7.9 0.2 4.9 -2.1 4.5 7.9 11.0 17.8 422 1
tau 6.3 0.3 5.0 0.2 2.5 5.2 8.9 18.7 214 1
eta[1] 0.4 0.0 0.9 -1.5 -0.2 0.4 1.0 2.1 928 1
eta[2] 0.0 0.0 0.9 -1.8 -0.6 0.0 0.5 1.8 1640 1
eta[3] -0.2 0.0 1.0 -2.1 -0.8 -0.2 0.4 1.8 1243 1
eta[4] 0.0 0.0 0.9 -1.7 -0.6 0.0 0.6 1.7 1421 1
eta[5] -0.3 0.0 0.9 -2.0 -1.0 -0.4 0.3 1.5 883 1
eta[6] -0.2 0.0 0.9 -2.0 -0.8 -0.2 0.4 1.6 926 1
eta[7] 0.4 0.0 0.9 -1.4 -0.2 0.4 0.9 2.1 969 1
eta[8] 0.1 0.0 1.0 -1.8 -0.6 0.1 0.7 2.0 1365 1
theta[1] 11.4 0.3 8.1 -1.4 5.9 10.3 15.2 30.6 574 1
theta[2] 7.7 0.2 6.1 -3.7 3.9 7.8 11.4 19.5 762 1
theta[3] 5.8 0.3 7.9 -12.1 1.8 6.5 10.5 19.9 715 1
theta[4] 8.0 0.2 6.5 -5.4 3.9 8.1 12.3 20.2 977 1
theta[5] 5.0 0.3 6.7 -10.3 1.3 5.7 9.5 16.5 667 1
theta[6] 6.0 0.2 6.6 -8.4 2.0 6.2 10.2 18.6 976 1
theta[7] 10.8 0.3 6.8 -1.1 6.2 10.2 14.9 26.0 596 1
theta[8] 8.6 0.3 7.9 -6.1 4.0 8.1 12.6 27.7 629 1
lp__ -5.0 0.1 2.6 -10.7 -6.6 -4.8 -3.1 -0.5 367 1

Row names are “mean”, “se_mean”, etc.; column names are scalar parameters.

> Don't you like structured objects here? So if "Sigma" is an M x N
> parameter, how do I access the summary stats for Sigma[m,n]?

Sure, we could do that, but that seems like a lot of work to me and it’s not clear that there are many cases where people will want to access the summaries in that way. That said, perhaps it would be logical to have a primitive function that computes mean, se_mean, sd, quantiles, n_eff, and Rhat for a single scalar parameter. Also I’m not actually quite sure how we are computing n_eff, Rhat, etc. for constrained parameters such as cov matrices.

>> foo.waic <- stan_waic_xval(foo.fit.multiple): returns a list of waic, elpd_waic, p_waic, elpd_loo, p_loo, pointwise, total, and se. This function only works if one of the returned parameters in foo.fit.multiple is something called log_lik. If there is no log_lik, it returns NULL. We could start with the function on page 14-15 of this paper: http://www.stat.columbia.edu/~gelman/research/unpublished/waic_stan.pdf
>>
>
> So you want a list rather than an object so you're OK with users
> getting p_loo using the syntax foo.waic$p_loo or foo.waic[[“p_loo"]]?

I don’t have a strong feeling here. It somehow seems like overkill to have extractor functions for all of these, but maybe that would be simpler, to just have functions waic(), elpd_waic(), p_waic(), elpd_loo(), p_loo(), etc. This might be easiest for the user. But if we do it that way, we will have to relabel “pointwise”, “total”, and “se”. Also we have to decide if we want the “stan_” prefix for all these functions.

> Unless p_loo corresponds to a leave-one-out cross-validation, this
> seems like a confusing name.

Yes, loo is leave-one-out cross-validation, as described in the paper.

> We'll need doc for what each of these elements is. What I really want
> in all of our RStan and PyStan doc is:
>
> argument types: what kind of R object do you pass in?
>
> argument content: what is the argument?
>
> argument constraints: what constraints on the arguments passed in
> are there, either singly or in conjunction.
>
> return type: what type of object does it return
> what is the content of that object field by field
>

> And speaking of names, looks like you're torn between using something
> like "_summarize()" that's a verb and "_summary()" which is the current
> function name in R.

I don’t recall using “summarize”. I was trying to use “summary” to be roughly consistent with R.

[Bob Carpenter]

The edit button is to the right of the page title on the top
next to the green button that says "new page".

- Bob

On Jun 27, 2014, at 9:54 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

...

> I did log into Github but I could not find an edit button.

...

[Bob Carpenter]

I think the confusion is that there's not a matrix per-se that's output here.
It's the behavior of the print() function on the fit object. The fit object
separates the samples by chain and structures them by parameter, I believe.
But I'm not familiar with the internals.

We really want to nail down the details rather than have them trickled
out as you realize there are things you don't like. Unless you literally
mean exactly this output is OK for you.

Should the quantiles be progammable?

Would you rather defaults other than 2.5%, 25%, 50%, etc? If so, should
we change CmdStan's print to match (it's now just doing 90% intervals and
50% quantile, I think).

Is the camel-case "Rhat" OK along with the C++-like "n_eff"?

>> Don't you like structured objects here? So if "Sigma" is an M x N
>> parameter, how do I access the summary stats for Sigma[m,n]?
>
> Sure, we could do that, but that seems like a lot of work to me and it's not clear that there are many cases where people will want to access the summaries in that way. That said, perhaps it would be logical to have a primitive function that computes mean, se_mean, sd, quantiles, n_eff, and Rhat for a single scalar parameter. Also I'm not actually quite sure how we are computing n_eff, Rhat, etc. for constrained parameters such as cov matrices.

I'm not sure which is easier given the existing structure.
But it is a question of what you want. With only a flat
structure, it's a pain to pull out a summary for just an array.

Should we be able to create sub-summary objects? That is,
take a summary, give it a sequence of parameter names and
get a new summary back? I don't know if that's how print works
or whether print gets more arguments for the things to print.

The point is that all of this needs to be specified and documented.
You can look at the existing doc now, which is pretty thorough, to decide
if things behave the way you want them to.

I'm pretty sure we're computing n_eff and Rhat with the constrained
values directly. If that's wrong, we need a deeper rethink of what's
going on. They are all (inverse) transforms of the unconstrained
parameters.

>>> foo.waic <- stan_waic_xval(foo.fit.multiple): returns a list of waic, elpd_waic, p_waic, elpd_loo, p_loo, pointwise, total, and se. This function only works if one of the returned parameters in foo.fit.multiple is something called log_lik. If there is no log_lik, it returns NULL. We could start with the function on page 14-15 of this paper: http://www.stat.columbia.edu/~gelman/research/unpublished/waic_stan.pdf
>>>
>>
>> So you want a list rather than an object so you're OK with users
>> getting p_loo using the syntax foo.waic$p_loo or foo.waic[["p_loo"]]?
>
> I don't have a strong feeling here. It somehow seems like overkill to have extractor functions for all of these, but maybe that would be simpler, to just have functions waic(), elpd_waic(), p_waic(), elpd_loo(), p_loo(), etc. This might be easiest for the user. But if we do it that way, we will have to relabel "pointwise", "total", and "se". Also we have to decide if we want the "stan_" prefix for all these functions.

OK --- I was just reacting to you saying you didn't want to see
the $ signs.

I think the point is that they'd be methods on an object, so you'd
doc them along with the WAIC return object, not as standalone functions.

For example, you don't see "print" documented directly when you look at the
RStan index --- you see print.stanfit. So it'd be something like
total.waic in the index if we added a total() method to the WAIC object.
Assuming you want the return type of doing WAIC to be an object.

I have no idea what people usually do in R, what users are comfortable
with, and what's the recommended practice.

>
>> Unless p_loo corresponds to a leave-one-out cross-validation, this
>> seems like a confusing name.
>
> Yes, loo is leave-one-out cross-validation, as described in the paper.

I see --- it's "truncated importance-sampling LOO". That's a bit of
a mouthful in a variable name.

>> We'll need doc for what each of these elements is. What I really want
>> in all of our RStan and PyStan doc is:
>>
>> argument types: what kind of R object do you pass in?
>>
>> argument content: what is the argument?
>>
>> argument constraints: what constraints on the arguments passed in
>> are there, either singly or in conjunction.
>>
>> return type: what type of object does it return
>> what is the content of that object field by field
>>
>> And speaking of names, looks like you're torn between using something
>> like "_summarize()" that's a verb and "_summary()" which is the current
>> function name in R.
>
> I don't recall using "summarize". I was trying to use "summary" to be roughly consistent with R.

The point is that R's not consistent. It uses "optimize" and "summary", which
are a verb and a noun. The optimize function creates an list object and
the summary function returns a value that depends on what it's summarizing.

This is the same issue as logLik, which is what R calls it and what we'd need
to call it as a method in order for R's general logLik() function to work on
our fit objects the way they do on a glm() output. We can either go with what R uses for
compatibility with R or we can try to make Stan internally consistent.

It's not a huge deal --- I'm resigned to the fact that naming can never
be made consistent because of inconsistent existing standards (at least
that's the case in Java, C++, R, and everything else I've ever used).

- Bob

[Andrew Gelman]

I see no green button. Perhaps I’m looking at the wrong thing? I’m here:
https://github.com/stan-dev/rstan/wiki/Interfaces-3.0-Spec

[Marco Inacio]

An owner (Daniel or Bob, I think) need to add you as a member in GitHub.

[Bob Carpenter]

I added Andrew.

Andrew and Marco --- you should both have write access
everywhere on stan-dev on GitHub.

Let me know if you still can't edit the wiki page.

- Bob

[Daniel Lee]

Thanks, Bob.

I thought everyone was on a team that had write permission to wikis. If that's not true, I'll need to revisit permissions. Please let me know if you can't edit the wiki pages.