Stan

Stan is a DSL for implementing Bayesian models. The model is composed in Stan, and then compiled and called from an interface in another language—in R, the rstan package. The Stan code compiles to C++, and rstan makes parallelization across chains simple.

It has its own shirts.

greta

greta is an R package by Nick Golding that implements MCMC sampling in TensorFlow, and can thus run on GPUs and scale well, provided the appropriate hardware and drivers. It is built on top of RStudio’s tensorflow package.

It does not have shirts, as far as I am aware, but does have a nice logo that would make a handsome hex sticker.

Stan

The Stan model has to be composed in, well, Stan. This can be in a separate file or a string, but I find the former usually too disconnected from the rest of the code and the latter a magnet for typos. The better solution in my attempts thus far is to implement the model in an RMarkdown code chunk with an output.var parameter set. When the chunk is run, the model is compiled and assigned to an R object with the supplied name. Code highlighting, completion, and so on work as usual.

Stan is declarative imperative¹, with chunks to define input data, parameters, the model, and more as necessary. Keeping it very simple, we will supply

• n, a number of observations,
• x, our single feature (a vector of length n), and
• y, our output (another vector of length n).

The model will calculate

• beta0, the y intercept,
• beta1, the slope, and
• sigma, the standard deviation of y.

While we could supply informative priors by adding lines like beta0 = normal(30, 2);, we will instead leave it blank, so Stan will use a uniform prior, determining a plausible range by running lots of warm-up iterations.

Our model is just a line. We assume our outputs are normally distributed, with a mean defined by the line, and a standard deviation parameterized by sigma.

All together, and set to output to an R variable called model_stan:

data {
    int<lower=0> n;
    vector[n] x;
    vector[n] y;
}

parameters {
    real beta0;
    real beta1;
    real<lower=0> sigma;
}

model {
    y ~ normal(beta0 + beta1*x, sigma);
}

Compiling takes about 40 seconds on my computer. It also produces 15 warnings that

warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas] #pragma clang diagnostic pop

but apparently they’re harmless, so I set message=FALSE.

To run the model, we just load rstan and call sampling on the compiled model, supplying the data. I set the algorithm to Hamiltonian Monte Carlo and the number of chains to 1 so as to compare more equally with greta.

library(rstan)

draws_stan <- sampling(model_stan, 
                       data = list(n = nrow(mtcars),
                                   x = mtcars$wt, 
                                   y = mtcars$mpg),
                       seed = 47,
                       algorithm = "HMC",
                       chains = 1)
#> 
#> SAMPLING FOR MODEL 'mtcars' NOW (CHAIN 1).
#> 
#> Gradient evaluation took 2.3e-05 seconds
#> 1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds.
#> Adjust your expectations accordingly!
#> 
#> 
#> Iteration:    1 / 2000 [  0%]  (Warmup)
#> Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Iteration: 2000 / 2000 [100%]  (Sampling)
#> 
#>  Elapsed Time: 0.23787 seconds (Warm-up)
#>                0.034739 seconds (Sampling)
#>                0.272609 seconds (Total)

draws_stan
#> Inference for Stan model: mtcars.
#> 1 chains, each with iter=2000; warmup=1000; thin=1; 
#> post-warmup draws per chain=1000, total post-warmup draws=1000.
#> 
#>         mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
#> beta0  37.02    0.08 1.53  34.08  36.02  36.96  37.99  40.10   402 1.00
#> beta1  -5.27    0.02 0.46  -6.15  -5.57  -5.25  -4.96  -4.43   389 1.00
#> sigma   3.10    0.01 0.36   2.46   2.86   3.08   3.31   3.94  1000 1.00
#> lp__  -51.08    0.10 1.28 -54.51 -51.61 -50.73 -50.17 -49.65   179 1.01
#> 
#> Samples were drawn using HMC(diag_e) at Fri May  4 19:52:44 2018.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at 
#> convergence, Rhat=1).

Since this is a very small model, it runs nearly instantaneously (0.26 seconds). The print methods of the call and resulting model are nicely informative. It looks good, but let’s check the diagnostic plots:

plot(draws_stan)
#> ci_level: 0.8 (80% intervals)
#> outer_level: 0.95 (95% intervals)

traceplot(draws_stan)

The intercept is moving more than the slope or standard deviation, which is reasonable considering the dimensions. The data is not totally linear (it could be better fit with a quadratic function), but the traceplots are fuzzy caterpillars, so it seems all is well.

Let’s extract the data and plot:

draws_stan_df <- as.data.frame(draws_stan)

ggplot(draws_stan_df, aes(beta0, beta1, color = sigma)) + 
    geom_point(alpha = 0.3) + 
    geom_density2d(color = "gray30") + 
    scale_color_viridis_c(option = "plasma") + 
    labs(title = "Stan parameter space")

ggplot(mtcars, aes(wt, mpg)) + 
    geom_point() + 
    geom_abline(aes(intercept = beta0, slope = beta1), draws_stan_df, 
                alpha = 0.1, color = "gray50") + 
    geom_abline(slope = mean(draws_stan_df$beta1), 
                intercept = mean(draws_stan_df$beta0), 
                color = "blue", size = 1) + 
    labs(title = "Bayesian regression on mtcars via Stan")

Looks quite plausible for the data at hand.

Stan

Stan is actually another language, although it is not particularly complicated or hard to learn. Since it is declarative imperative, it is easy to see what goes where in the model, but not so easy to figure out where bugs are when compilation fails. Error messages (aside from the annoying warnings) are pretty good, though.

Written in an RMarkdown chunk, Stan is really easy to integrate into a larger workflow without breaking stride. There is also the rstanarm package, which assembles up typical models into R functions. While this seems like it should be simpler, a quick attempt showed it to actually be somewhat more finicky, as the resulting functions have a lot of parameters and are fairly complicated compared to the smaller, simpler individual pieces of stan itself.

The biggest problem I ran into is that since Stan is another language, the documentation is not available with ?. In fact, it is a 637-page PDF. The PDF is in fact very well-written and organized, but still, it’s a PDF, which is annoying as it does not integrate into the web well. (I can’t put a link to the section on linear regression here.) Hopefully the documentation will eventually be better integrated into the otherwise-great Stan website. In an ideal world, it could be integrated into R’s help, too, though that would require some work.

greta

greta avoids a lot of the issues of Stan. Documentation is available as usual. Intermediary products are easily viewable, and an individual problematic call will fail, making debugging simpler. It is not quite as easy to look at the code and quickly see the model structure, but the models plot method is a pretty solution to this problem.

In terms of organization, greta almost unavoidably leaves a lot of objects in the global environment with no inherent organization. The objects will not fit neatly in a data frame, so there is no natural organizational approach. Everything could (and perhaps should) be kept in a list or environment. A more robust approach would be to make a class with slots for the parts of the model and a suitable print method. While more work to build and apply, such an approach would allow quicker manipulations of multiple models without mess or confusion.

The biggest drawback in my short experiment is speed. I ran simulations on slightly larger (not “big”, in any sense) data, and it had no effect on run time, which seems to be purely a function of the number of iterations. I am not sure how it scales, but given that it was explicitly built for such a purpose, hopefully what I experienced is just overhead.

Overall, greta is a bit less mature than Stan, which has a bigger, funded team working on it. Given that greta leverages other projects that have such an infrastructure, it may not require such dedicated resources, but it does have room to grow yet. Hopefully it will!

Bayesian simulation

For a relative beginner, both interfaces proved surprisingly simple and easily adaptable. I—as ever—have much yet to learn, but these frameworks will doubtlessly prove to be great assets with which to further explore the potential of Bayesian statistics.