Announcing gsim 0.1

[Lionel Henry]

Dear Stanimals,

gsim is a new R package to post-process posterior simulations. Its goal is to make it easy to create numerical and graphical summaries involving posterior parameters. In particular, it should  be helpful to perform posterior predictive checking.

In a nutshell, gsim abstracts away all loops over simulations so that the code reflects the math and not technicalities like array instantiation, subsetting, etc. It also provides tools to easily plot fitted models while displaying posterior uncertainty.

The project page is on github: http://github.com/lionelgit/gsim

You will find a tutorial here: https://dl.dropboxusercontent.com/u/239966/gsim.html

I've put in place unit tests to make sure that gsim outputs are reliable, but for the moment you should of course keep an eye on the results and watch out for bugs. I am interested in your comments on the design, potential features, and any issues.

Lionel

[Bob Carpenter]

This reminds me very much of Jouni Kerman's rv (random variable)
package for R:

http://cran.r-project.org/web/packages/rv/index.html

The idea is that random vars are represented as simulations and can
be combined by various operations that respects the random variableness,
such as (alpha + beta) doing the right thing (adding simulations
iteration by iteration).

I'm not enough of an R expert to comment on the rest. For what it's
worth, I don't find the loops that onerous. A

Within Stan, we're always torn between

1. doing the simulations within the generated quantities block, and

2. doing them on the outside (doing the loop version of what gsim's doing).

A drawback to (1) is that you have to figure out what you need ahead of time.
But that comes with a benefit that you get all the convergence monitoring
automatically, and things like R-hat and n_eff are sensitive to non-linear transforms.

(2) is complementary to (1) in that you can figure out what you need later,
but then you have to work to get R-hat and n_eff and what not.

The major drawback to (2) [from a Stan perspective] is that it's done in an
outside language (R, Matlab, Python, etc.), and then you're doing many of the same
calculations as in the original Stan model, but with different syntaxes and perhaps
different parameterizations.

Of course, the benefit of doing something like (2) in an external R package
like gsim is that it'll work for other forms of output than Stan's.

We're having the same issue with how to organize posterior visualizations and
analysis (similar to what Coda does). We have lots of R developers who want
to write gsim-like packages external to Stan, the most recent example being the
Shiny visualizations discussed on stan-users a week or two ago.

- 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.

[Lionel Henry]

Thanks for the feedback Bob.

This reminds me very much of Jouni Kerman's rv (random variable) package for R

Yes, gsim was first implemented using a similar object oriented approach. However, this was not flexible and rather slow for complex operations on many simulations. It is not flexible because a generic method is needed for every function and operator to be used with posterior objects. The reactive functions in particular were cumbersome to implement and really slow (they had to precompute a grid of possible values over some range defined ahead of time).

So after reading Wickham's book on R programming I reimplemented everything using a 'computing on the language' approach. This makes gsim fully flexible for any array operations (with some current limitations such as subsetting, but that should work in the future), and much faster. I also may rewrite it in Rcpp in the long term, which would make it faster than using R loops.

For what it's worth, I don't find the loops that onerous.

Project files quickly become a mess though. I think it's much better when each line of code is a direct translation of a mathematical idea. Also, an important part of gsim is to save snippets of code to be reused with different models. So this package aims to make project files neater and easier to read, as well as more organized.

Within Stan, we're always torn between 
1.  doing the simulations within the generated quantities block, and
2.  doing them on the outside (doing the loop version of what gsim's doing).

I think it's not so much a tension than an emulation. Model exploration and model checking (which is also exploratory in its nature) are a critical part of bayesian analysis. I don't think that Stan can (nor it should) offer a framework to perform these exploratory tasks.

The major drawback to (2)  [from a Stan perspective] is that it's done in an outside language (R, Matlab, Python, etc.), and then you're doing many of the same calculations as in the original Stan model, but with different syntaxes and perhaps different parameterizations.

One great thing with Stan is its compatibility with many platforms in which post-processing of simulations and exploration of models take place. gsim is merely an addition to the R platform so that working with simulations is seamless. It does not compete with Stan's generated quantities. Once the exploration is done, the model can be refitted in Stan to obtain convergence diagnostics on critical quantities.

We're having the same issue

I don't think this is an 'issue' though. Why strive for Stan's compatibility with many frameworks and then be worried that these frameworks develop tools to work with Stan's outputs? As I see it, ShinyStan is only a way of displaying information provided by Stan.

Lionel

[Tamas Papp]

Lionel, thanks for writing this! I have been looking for approaches to
posterior predictive checks [1] for a while, since I found rv cumbersome
in practice.

[1] https://groups.google.com/d/topic/stan-users/eOf6oOm9xIY/discussion

I looked at gsim briefly and I think your approach is preferable to
overloading the standard operators (like rv does). However, I have to
say that I find the whole non-standard evaluation concept very brittle
in R -- I know that Wickham uses it extensively, but it becomes tricky
when you are trying to do something nontrivial (reuse code, etc).

So for the moment I will continue to use loops -- I don't find them
cumbersome per se, but I think that more abstraction is less
error-prone, but I will follows gsim with interest.

IMO posterior predictive checks can be best conceptualized as mapping
tuples to tuples (lists in R) with a function. A library would help in
extracting the tuples from the posterior for each draw, and then
organizing the results in containers. In R, for efficiency, posteriors
are lists of arrays instead of vectors of lists, so the key is going
back and forth between these two representations, and making sure
variable names are mapped correctly, which is where computing on the
language comes in. I am wondering if the latter could be more
transparent and composable with a bit of extra work on the part of the
user: I imagine something like

map_posterior(function, posterior, inargs, outargs)

where inargs and outargs would be alists (list with names elements in R)
specifying variable names and the corresponding columns, and the
function would do the following:

1. extract a draw from the posterior,
2. arrange the arguments and do_call the function,
3. for the first call, guess dimensions from the list that was returned
and create the containers, for subsequent calls, just check consistency,
4. save the results for that draw,
5. return them at the end.

The benefit of this is that there is no non-standard evaluation, just a
function application. I wrote this in R but it turns out to be
hopelessly slow (in my naive implementation, I am not an R wizard). If
you have a similar core function in gsim, or something that can be
abstracted like this, maybe you could consider exposing it so people
could contribute an alternative interface with standard evaluation.

Best,

Tamas

[Lionel Henry]

Hello Tamas,

I am not sure I undertand well your idea of mapping tuples to tuples. I also tried to think about a way of abstracting ppc in functions, but I had trouble finding a solution that would be both flexible, and easier to read than simply spelling out the checking.

At one point I thought of something like this:

  bernoulli_check(y = switched, p = fitted, stat = sd(y - p))

Well I just went ahead and implemented it :-) See in https://github.com/lionelgit/gsim/blob/master/R/summaries.r.

It only required a handful of tweaks to the expression parser. Besides `y`, `p` and `stat`, it takes a `p_value` argument to indicate whether it should return n_sims p values (the default) or a list containing y and y_rep (returning a list is not 100% clean at the moment because the code that handles listed outputs has not yet been adapted).

This function uses both SE and NSE. The SE part computes replicates. The NSE part takes the 'stat' argument to evaluate it on both y and y_rep, for each simulation.

Here it is in action:

  sims <- gsim(wells_fit, wells_data, n_sims = 200)

  p <- sims({

    X <- cbind(intercept(), c_dist100, c_arsenic, c_dist100 * c_arsenic)

    fitted <- inv_logit(X %*% col_vector(beta))

    bernoulli_check(y = switched, p = fitted, stat = sd(y - p))

  })

It could easily be adapted for other distributions. However I'm not sure how this approach would play out with more complicated checks for hierarchical models etc.

> I find the whole non-standard evaluation concept very brittle in R

At first I also thought that computing on the language to wrap everything into loops before execution would be brittle. However, when I tested a rough implementation, it worked surprisingly well. One thing that helps is that every input is an array or a numeric vector, and that most function used inside gsim operate only on this type of inputs.

> it becomes tricky when you are trying to do something nontrivial (reuse code, etc).

I'm not sure what you mean by something nontrivial and by reusing code. Would quoted snippets fit the use cases you have in mind?

Lionel

[Bob Carpenter]

Thanks for these detailed discussions. It's very helpful.
I hope Andrew reads this because he's thought a lot about these
issues.

I've never used rv, but we also overload all the functions/operators in
Stan to support autodiff, and it's definitely a lot of work to do right.
I had no idea you could do that kind of thing in R for anything other than
some of the OO methods. Can you overload addition for an R object type?

The portability issue for Stan is whether things will work across
interfaces in R, Python, MATLAB, Julia, etc. (as well as across C++
compilers and operating system and hardware configs). The main advantages
to working directly in our generated quantities block is that:

* it's portable across Stan interfaces, and

* you get vectorization across samples for free.

The latter is what makes Stan variables (other than data) feel like random
variables and thus makes Stan feel like a probabilistic programming language.
This is the same feel I get from rv and gsim (though I haven't used either),
as well as from BUGS and JAGS, the obvious influences on our language.

The downside for generated quantities is that

* it only works for Stan, and

* it requires the posterior analysis to be done along with the fit.

The latter's mainly an issue for flexibility, though it has the advantage
of making things like convergence monitoring integrated.

When we add all of the MLE functions lm(), glm(), MML functions lmer(), glmer(),
and Bayesian MML functions blmer(), and bglmer() directly to Stan (with
wrappers for the lm/lmer expression language in Stan's R package), there'll be
less motivation for anyone to use any other package :-) We already do MLE
and I'm guessing MML will land in the next six or nine months now that we
have both Rob and Marcus working on it.

What Stan users have been asking us for is something that makes it possible
to use one model for both estimation and simulation. Unlike
Stan isn't a pure directed graphical modeling language, and it compiles
code dedicated for its inferential task. This makes the code a lot faster,
but it has the enormous drawback of requiring users to write
two models where BUGS/JAGS would only require one. There is the bright
side that each model is easier to read for intent, but that's a separate
matter.

- Bob

[Andrew Gelman]

Just to add to what Bob said: One thing we’ve been talking about is R or Python wrapper for Stan that does graphical models. The idea is that we’d create some sort of graphical modeling language and then the wrapper would take that model, along with a choice of which random variables are being specified as data and which are unknown, and then it would write a Stan program, putting the appropriate code into the model and generated quantities blocks. This would facilitate fake-data simulation and posterior predictive checking, which are two things we want to be doing all the time.

A

P.S. For an example of what is and is not a graphical model, consider the following 2 lines of Stan code:

theta ~ normal(0,1);
theta ~ normal(0,1);

[Or, how I’d prefer it written:
normal(theta|0,1);
normal(theta|0,1);]

This code, although it looks like two specifications of a prior for theta, is actually in Stan just two contributions to the log-posterior, and it is mathematically equivalent to:

theta ~ normal(0,1/sqrt(2));
[or, as I’d prefer, normal(theta|0,1/sqrt(2));]

[Bob Carpenter]

I didn't want to subvert the stan-users thread, so I started
a new one.

It'd be trivial to add

foo(y | theta1, ..., thetaK)

as an alternative syntax for

y ~ foo(theta1, ..., thetaK);

And we can deprecate the latter, though it'll be a bit of
work changing everything, since the syntax is involved in almost
every single model we've written or a user's written.

I think of the foo(y|theta) notation as an expression (specifically
a function application) and thus evaluating to a value without
side effects, so it may be confusing just using it as a bare statement
that implicitly increments the total log density.

- Bob

> On Nov 26, 2014, at 9:09 PM, Andrew Gelman <Gel...@stat.columbia.edu> wrote:
>
> P.S. For an example of what is and is not a graphical model, consider the following 2 lines of Stan code:
>
> theta ~ normal(0,1);
> theta ~ normal(0,1);
>
> [Or, how I’d prefer it written:
> normal(theta|0,1);
> normal(theta|0,1);]

...

[Tamas Papp]

I find the ~ syntax much better. My eyes parse `a | b` as a single
argument to a function (probably the influence of R's formula syntax).

Sure, I could get used to | in Stan if everyone prefers that, but I
don't think ~ is confusing: I think it is reasonable to expect that a
Stan user has some minimal level of understanding about Bayesian
inference and how Stan works, eg from reading at least the very
beginning of the manual --- it is pretty much hopeless to use Stan
without that level of understanding.

Best,

Tamas

PS.: didn't cross-post to stan-dev as I don't have posting privileges
there.

[Joshua Pritikin]

On Thursday, November 27, 2014 2:21:31 AM UTC-5, Bob Carpenter wrote:

It'd be trivial to add

  foo(y | theta1, ..., thetaK)

as an alternative syntax for

  y ~ foo(theta1, ..., thetaK);

I think of the foo(y|theta) notation as an expression (specifically
a function application) and thus evaluating to a value without
side effects, so it may be confusing just using it as a bare statement
that implicitly increments the total log density.

I don't follow how foo(y|theta) is an expression and y ~ foo(theta) is

not an expression. Neither of them mention the total log density.

To me, the distinction is that the tilde implies simulating draws

from foo(theta). The way to generate a tilde in latex to use $\sim$

suggesting simulation. And in stan models, we are not doing

simulation, but we are really saying that the model should account

for the density foo(y|theta).

I'm not sure what to do about the implicit accumulation of these

densities. But that's a separate question. I think the question is

what can go inside of increment_log_prob(...).

increment_log_prob(y ~ normal(mu,sigma));

increment_log_prob(normal_log(y,mu,sigma));

increment_log_prob(normal_log(y | mu,sigma));

[Tamas Papp]

On Wed, Nov 26 2014, Lionel Henry <lione...@gmail.com> wrote:

> I am not sure I undertand well your idea of mapping tuples to tuples. I
> also tried to think about a way of abstracting ppc in functions, but I had
> trouble finding a solution that would be both flexible, and easier to read
> than simply spelling out the checking.

Hi Lionel,

Sorry for not being clear. Here is stylized code for what I have in
mind, using the speed of light example from BDA (2nd ed, Chapter 6):

--8<---------------cut here---------------start------------->8---
y <- is_the_data

my_statistic <- function(y, mu) {
y_sorted <- sort(y)
abs(y_sorted[61]-mu)-abs(y_sorted[6]-mu)
}

my_rep <- function(mu, sigma) {
rnorm(length(y), mu, sigma)
}

## with a loop
M <- length(posterior$mu)
T_data <- numeric(M)
T_sim <- numeric(M)
for (i in 1:M) {
T_data[i] <- my_statistic(y, posterior$mu[i])
T_sim[i] <- my_statistic(my_rep(posterior$mu[i], posterior$sigma[i]),
posterior$mu[i])
}

## with a hypothetical function map_posterior
map_posterior(function(mu, sigma) {
list(T_data = my_statistic(y, mu),
T_sim = my_statistic(my_rep(mu, sigma), mu))
},
posterior,
list(mu = "mu", sigma = "sigma")) # input mapping, may not be necessary
## result is a data frame with two columns.
--8<---------------cut here---------------end--------------->8---

map_posterior would take care of calling the function with each draw in
the posterior, and collecting the results etc. I find functional
approaches very clean.

>> I find the whole non-standard evaluation concept very brittle in R
>
> At first I also thought that computing on the language to wrap everything
> into loops before execution would be brittle. However, when I tested a
> rough implementation, it worked surprisingly well. One thing that helps is
> that every input is an array or a numeric vector, and that most function
> used inside gsim operate only on this type of inputs.
>
>
>> it becomes tricky when you are trying to do something nontrivial (reuse
> code, etc).
>
> I'm not sure what you mean by something nontrivial and by reusing code.
> Would quoted snippets fit the use cases you have in mind?

As you mention in the documentation, control flow does not work for
gsim, but I don't know if this is an inherent limitation.

I usually run into problems where I have functions that call functions
which have nonstandard evaluation, since it is hard to parametrize the
latter (Wickham has, of course, packages that work around this, like
lazyeval). But this is not super-relevant for posterior predictive
checks I guess.

Best,

Tamas

[Lionel Henry]

Hi Bob,

> I had no idea you could do that kind of thing in R for anything

> other than some of the OO methods. Can you overload addition for an

> R object type?

Yes, most operators are generic functions so you can specify methods

for them. Actually, the main OO systems in R (S3 and S4) are based on

overloading generics. The methods do not belong to classes but to

generic functions. Dispatching then takes place based on the class of

the first argument (S3), or based on the classes of the arguments of

the method (S4).

You can also entirely overwrite addition. In R everything is a

function, even operators, so if your new addition function has scoping

precedence, it will prevail.

>  We already do MLE and I'm guessing MML will land in the next six

> or nine months now that we have both Rob and Marcus working on it.

Awesome news!

[Lionel Henry]

Hello Andrew,

> Just to add to what Bob said: One thing we’ve been talking about is

> R or Python wrapper for Stan that does graphical models.

Having such a tool to generate model files would be great and would

indeed make life easier to get replicates and test stats from Stan.

However, I still think that there is a need for explorative tools that

Stan cannot fulfill on its own. It's disruptive to go back and refit a

model in Stan whenever you think of a new aspect of the model to

check, especially if the sampling is slow.

There are also other explorative tasks that gsim tries to make

easier. In particular I spent a lot of effort to implement the

reactive outputs (they may need better names), to help creating

functions to display fitted lines with posterior uncertainty.

Lionel

[Lionel Henry]

Hi Tamas,

Thanks for the example. My function bernoulli_check is not very

different from your idea. The interface uses NSE for the statistic

function, but you can still define and use your own external statistic

functions. So maybe that meets your idea of a 'functional' solution?

By the way, I uploaded normal_check on Github. I also tweaked both

checking functions to make them independent of gsim. This was easy

because gsim uses the same data structure than Stan's extracted

parameters. So now you can directly use them with simulations

extracted from a stanfit.

I also changed the output of these functions. They now return a list

with `y_stat` and `y_rep_stat` with the p-value recorded as an

attribute. I think it's better to have the list returned by default

because checking the test graphically requires the statistics. And the

p-value is still easy to access because I provide a getter function

called `p`: p_value <- p(my_check).

So, for the speed of light example, you would use it as follows:

sims <- extract(speed_stanfit)

symmetry_check <- normal_check(

  y = newcomb, mu = sims$theta, sigma = sims$sigma,

  stat = {

    sorted_y <- sort(y)

    abs(sorted_y[61] - mu) - abs(sorted_y[6] - mu)

  })

p(symmetry_check)

Of course you can just use stat = symmetry_stat(y, mu) if you have

already defined one. And to make everything a bit nicer, you can use

the %$% operator from the package magrittr to expose the data and the

parameters to normal_check (equivalent to use the `with`

function). Then it would be:

symmetry_check <- c(sims, list(y = newcomb)) %$%

  normal_check(y, mu = theta, sigma, stat = symmetry_stat(y, mu))

The problem with the current implementation of these checking

functions is that you can only operate on the parameters of the

data distribution. So, for normal_check, the stat can only depend on

mu and sigma. Does someone has an example of a checking stat that

depends on other parameters than those of the data distribution?

I thought of this issue while I was thinking of a more generic

model_check() that would take y and y_rep as inputs. I thought this

checking function would be useful for the common case where the

replicates are generated with Stan. But a checking function that

cannot operate on any parameters is actually not so practical :-).

However I could make it work inside gsim. Having:

sims(

  model_check(y, y_rep, my_stat(y, theta))

)

is a bit nicer than:

sims({

  y_stat <- my_stat(y, theta)

  y_rep_stat <- my_stat(y_rep, theta)

  list(y_stat, y_rep_stat)

})

Lionel

[Bob Carpenter]

Syntactically, in a programming language, expressions denote
values and statements denote actions.

So

sin(x + 3)

is an expression, whereas

y = sin(x + 3);

is a statement.

Languages like C like to overload things, so that the above statement
may be used as an expression with value sin(x + 3). Similarly, the expression
may be used as a no-op statement, though no reason to do that if there are no
side effects.

You can see how the distinction's made in Stan in the reference
section of the manual, which has separate chapters on statements
and expressions, as well as the BNF syntax for the language which
has the precise syntactic definition.

All I was trying to say is that foo(y|theta) "feels" more like
an expression to me.

- Bob

[Michael Betancourt]

I agree with this feeling. We’re in an awkward spot here because of
Andrew’s example where you can have multiple “sampling” statements,

y ~ p(theta)
y ~ p(theta)

in which case “y is distributed as p(theta)” doesn’t make as much sense.
Personally I still think

y ~= p(theta)
y ~= p(theta)

gets the point across more, but at some point we may be too burdened
by the BUGS tradition.

[Lionel Henry]

> The problem with the current implementation of these checking

> functions is that you can only operate on the parameters of the data

> distribution.

Ok, this was actually easy enough to make the checking functions able

to operate on any parameters. I had to add a `params` argument. The

tricky part internally (and also the speed bottleneck) was to find a

way of subsetting arrays of varying dimensions. With this I could

implement `check_model()`. Note that I renamed all checking functions

from *_check to check_* to make them verbs.

So now the global idea behind the checking functions is that you apply

a statistic on both y and y_rep. If you already have y_rep, you call

check_model. If you need to sample replicates, you call the

specialised functions like check_normal or check_bernoulli. Actually,

the distribution-specific functions just sample replicates and then

call check_model themselves.

These functions work outside gsim on extracted stanfit

simulations. (Note that if you need to get simulations in this form,

you can use `as.mclist()` on any object supported in gsim.) They don't

work inside gsim yet, but when they do the only difference will be

that they won't have a params argument.

Example with check_model:

  newcomb_sims <- extract(stanfit)

  out <- check_model(y = radon$y, y_rep, params = newcomb_sims,

    stat = {

      y_sorted <- sort(y)

      abs(y_sorted[61] - y_hat) - abs(y_sorted[6] - y_hat)

    })

Example with check_bernoulli:

  wells_sims <- extract(stanfit)

  out <- check_bernoulli(y = wells$switched, p = y_hat,

    params = wells_sims, stat = sd(y - p)

  )

check_normal is similar to check_bernoulli, only with mu and sigma

arguments.

When params is supplied, NSE is used for the random variables (y_rep,

mu, sigma, p, etc). So if you specify mu = y_hat, they will look

inside params for y_hat. On the other hand, when params is not

supplied, the random variables are obtained from the calling

environment. In both cases, the statistic operates on 'official' param

names (mu, sigma, p), not on supplied names.

I may change that interface if it is confusing. I would love some

feedback from experimented ppc users. One remaining issue is speed,

although it's fast enough for practical purposes. Computing a bernoulli

check on 2000 simulations takes 1.6 seconds on my computer.

Rewriting check_model in c++ would solve this.

Anyway, thanks Tamas. I don't know if these functions would be a good

replacement of map_posterior for you but our discussion helped me see

things more clearly!

Lionel

[Lionel Henry]

I agree with Bob that foo(x | theta) feels like an expression that evaluates to something. Namely, I think of normal(x | mu, sigma) as a density function.

But what if the `model` block gets renamed to something suggestive like `posterior contributions` or `probability factors`? Then it would feel like the objects in this block are compounded in the posterior.

factors {

  normal(theta | 0, 1)

  normal(theta | 0, 1)

}

Lionel