Mitzi Morris sends along the paper trail (I love archaic idioms) for her data spelunking in advance of her all-day introductory Stan tutorial next week at StanCon (we’ll see you in Uppsala). She was trying to reproduce the results from my decade-old case study on hierarchical modeling when she ran into a couple discrepancies.

The case study was designed to provide a Bayesian replication of the results in Efron and Morris’s five decade old paper on hierarchical modeling (aka Stein’s estimator, aka population regularization, aka “empirical” Bayes), which is still under a paywall courtesy of our “friends” at the American Statistical Association. If your organization isn’t paying ASA for access to a paper that an academic donated for free 50 years ago, I’ll leave you to find your own pirated copy in good conscience, or you can follow the link and let Google hoist the Jolly Roger for you (now that’s an even more obscure and archaic reference).

• Bradley Efron and Carl Morris. 1975. Data Analysis Using Stein’s Estimator and Its Generalizations. JASA 70(350).

My case study’s been out for ten years. Mitzi found a problem when matching the data provided in the R package pscl against that in Efron and Morris’s paper. In particular, the data for the player named “Williams” was wrong. Efron and Morris manually transcribed the data from a newspaper with the goal of finding a bunch of players with the exact same number of at bats on a given day (45, it turns out). They did so accurately.

Sadly, I imported an error in the 18-item data set from the R package. Here’s an edited report from Claude after it tracked down what went wrong.

1. The 1975 paper is correct, although it only gives the surname “Williams.”

2. The R package pscl expanded that name incorrectly to Billy Williams, while retaining Walt’s team and remaining-season data. Published pscl data.

3. The Bean Machine tutorial () contains the exact mixed row found in your CSV:

     Billy,Williams,...,591,0.3299,636,205,0.251
   

4. That row combines Billy Williams’s at-bats and hits with Walt Williams’s .251 average. Thus, ID 15 in your data should be corrected throughout to:

     Walt,Williams,45,10,0.222,270,0.2556,315,79,0.251
   

But wait, that’s not all. There was another typo in the Bean Machine tutorial. They mistranscibed Ron Swoboda‘s surname as “Swaboda.” I should have caught that as I was a huge baseball fan in the early 1970s, am an excellent speller, and grew up around a lot of Poles and Czechs in Detroit and its western suburbs. Mitzi then had the following conversation with Claude after she noticed the typo. This may be the first time I’ve seen Claude crack a spontaneous joke.

Mitzi: did that [the typo] come also from the R file?

Claude: If you mean the Swaboda misspelling: no. The pscl R dataset correctly uses Ron Swoboda. The misspelling appears in the Bean Machine tutorial
(https://beanmachine.org/docs/overview/tutorials/Hierarchical_modeling/HierarchicalModeling/), which contains Ron,Swaboda.

Mitzi: are you sure? I’ve heard it both ways.

Claude: Yes—despite what Shawn Spencer might say, it is definitively Ron Swoboda, with an o.

Mitzi: who’s Shawn Spencer? in 2017, Bob Carpenter published this dataset in a Stan case study.

Claude: Shawn Spencer is the fictional detective in the TV comedy Psych, played by James Roday Rodriguez. “I’ve heard it both ways” is his recurring response when someone corrects him—hence my reference.

Mitzi and I use “I’ve heard it both ways” all the time, having watched umpteen seasons of Psych. I would recommend if you want a lighthearted detective story that doesn’t take itself too seriously with lovable characters and you can put up with the ubiquitous mid-00s background sexism and racism.

We’re very excited about this book. It’s the result of several years of effort. You can order from the publisher or from Amazon.

Here’s the book’s webpage, which includes the data and code for the book’s examples and case studies, of which there are many.

Here’s the table of contents:

Part 1: From Bayesian inference to Bayesian workflow
1. Bayesian theory and Bayesian practice
2. Statistical modeling and workflow
3. Computational tools
4. Introduction to workflow: Modeling performance on a multiple choice exam

Part 2: Statistical workflow
5. Building statistical models
6. Using simulations to capture uncertainty
7. Prediction, generalization, and causal inference
8. Visualizing and checking fitted models
9. Comparing and improving models
10. Statistical inference and scientific inference

Part 3: Computational workflow
11. Fitting statistical models
12. Diagnosing and fixing problems with fitting
13. Approximate algorithms and approximate models
14. Simulation-based calibration checking
15. Statistical modeling as software development

Part 4. Case studies
16. Coding a series of models: Simulated data of movie ratings
17. Prior specification for regression models: Reanalysis of a sleep study
18. Predictive model checking and comparison: Clinical trial
19. Building up to a hierarchical model: Coronavirus testing
20. Using a fitted model for decision analysis: Classification competition
21. Posterior predictive checking: Stochastic learning in dogs
22. Incremental development and testing: Black cat adoptions
23. Debugging a model: World Cup football
24. Leave-one-out cross validation model checking and comparison: Roaches
25. Model building and expansion: Golf putting
26. Model building with latent variables: Markov models for animal movement
27. Model building: Time-series decomposition for birthdays
28. Models for regression coefficients and variable selection: Student grades
29. Sampling problems with latent variables: No vehicles in the park
30. Challenge of multimodality: Differential equation for planetary motion
31. Simulation-based calibration checking in model development workflow

Appendices
A. Statistical and computational workflow for Bayesians and non-Bayesians
B. How to get the most out of Bayesian Data Analysis

One way to think of the book is that it’s all the things missing from BDA, like how to set up an informative prior, what to do when your computations aren’t converging, how to work through a series of models fit to the same data, how to design and perform simulated-data experiments . . . and all sorts of other things too.

The core of the book–parts 1 through 3–clock in under 200 pages, and then we have another 300 pages full of case studies demonstrating different aspects of Bayesian statistical and computational workflow. The appendices should be useful to you too, first because the workflow ideas in this book apply to non-Bayesian inference too, and second because BDA still has lots of valuable material in it, so it’s good to know where to look.

This new Bayesian Workflow book could change your life (we hope), and I thank my coauthors, Aki Vehtari and Richard McElreath, with Daniel Simpson, Charles C. Margossian, Yuling Yao, Lauren Kennedy, Jonah Gabry, Paul-Christian Bürkner, Martin Modrák, Vianey Leos Barajas, for all their care and effort. We thank our employers and various funding agencies for giving us the resources to be able to write this book as a side project along with all our daily responsibilities. And we thank many people for their input on earlier versions of the book, along with the Stan developers making so much of this work possible and the Stan community of users for supplying a continuing series of challenges that have motivated many of the ideas and methods discussed in the book.

I posted this already on the blog and you can see answers to some questions in the comments there. I’m posting it again here because, hey, we don’t come out with a new book every day!

I hope you find the book readable, interesting, and useful.

Lu Zhang writes:

As a member of the BayesComp 2027 conference committee, I would like to share the announcement of the call for invited session proposals for the upcoming BayesComp conference, which will be held in College Station, Texas, on May 18–20, 2027.

The scientific committee is currently soliciting proposals for invited sessions. Each invited session will consist of three speakers, and proposals should focus on timely, important, and broadly engaging topics in Bayesian computation and related areas.

The submission deadline (as of now) for invited session proposals is August 15, 2026.

Proposal form: https://forms.gle/wpYvkkjKGZ5vHqhF6

Additional details are available in the official announcement:

The LOC for BayesComp 2027 is pleased to announce that the next edition of BayesComp will take place in College Station, TX during May 18–20, 2027. The scientific committee is now opening calls for invited sessions. Each invited session will consist of 3 speakers. Proposals should highlight timely, important, and broadly engaging topics in Bayesian computation and related areas. Each speaker may be listed as a speaker in only one invited or contributed session proposal.

Lu is the first author on the Pathfinder paper and continues to do interesting work on Bayesian statistics and computing. Based on what I’ve heard about past BayesComps, the conference should be really interesting.

This post is from Bob

Nutpie laps NUTS

The short story is that the Nutpie sampler is usually twice as fast as Stan’s sampler for Stan models and twice as fast as PyMC’s samplers for PyMC models. The speed improvement is due entirely to better mass-matrix adaptation. Everything else is the same—dual averaging for step size adaptation and the biased-progressive, multinomial no-U-turn sampler for Hamiltonian Monte Carlo. Note that the speed is not coming from Python or Rust; we have a C++ implementation in the works we plan to interface with Python and then R.

The backstory

Adrian Seyboldt, a member of PyMC Labs, designed the sampler and released the open-source implementation linked above. He built it so he could try out Rust! Adrian’s also behind Stan’s game-changing new sum-to-zero parameterization along with Sean Pinkney. He’s a physicist by training and has amazing geometric intuition.

I tried Nutpie and it didn’t disappoint. So I wanted to understand the algorithm. I couldn’t understand the Rust code and Adrian hadn’t written any pseudocode or algorithm description. Eliot Carlson, with a recently minted MA in stats from Columbia, showed up at my office on the recommendation of Philip Greengard to look for a computational stats project. I asked him what he thought about Rust and like any great grad student, he wasn’t daunted by the fact that he knew nothing about Rust and was a statistician by training. With serious digging and some help from Adrian, Eliot extracted pseudocode descriptions and we were able to formulate the algorithm in what I hope is an understandable way. Eliot now works full time at PyMC Labs.

The paper

The arXiv paper just came out yesterday.

• Adrian Seyboldt, Eliot L. Carlson, Bob Carpenter. 2026. Preconditioning Hamiltonian Monte Carlo by minimizing Fisher Divergence. arXiv 2603.18845.

We’d love to hear any comments, corrections, or suggestions people might have.

My role in this was just helping to formulate the whole thing as a stats paper in a language that computational statisticians could understand. Adrian and Eliot worked out the proofs, from which I learned a lot.

Why is Nutpie better?

There are two basic reasons why Nutpie’s diagonal mass-matrix adaptation is better than Stan’s:

1. Better initialization: Use the gradients at the initial point to estimate mass matrix and make the whole sampler scale invariant.

2. Two sources of information: Use the position and the gradient of the log density at the position as two sources of information with which to estimate a mass matrix

The motivation for (1) is that the outer product of the scores is a one-sample Monte Carlo estimate of the expected negative Hessian of the target log density. The one-sample estimate is regularized by taking a geometric average with an identity matrix—using a few more draws here may be more stable.

Will wonders never cease, but it turns out to be relatively simple to take the geometric mean as the midpoint of a geodesic in the affine-invariant Riemannian manifold (AIRM) of positive-definite matrices. It works out to regular geometric mean in the diagonal case, making it simple to implement.

The motivation for (2) is that in a multivariate normal distribution with covariance Sigma, the covariance of the draws is Sigma and the covariance of the scores (gradients of log density of the draws) is Sigma inverse. These are then geometrically averaged to produce an estimate. This can converge very fast—many fewer iterations than the dimensionality of the problem.

As with Stan, the default in Nutpie is to use a diagonal mass matrix. Unlike Stan, Nutpie also lets you use a low-rank plus diagonal approximation (as in L-BFGS optimization), which can go a long way to decorrelating hierarchical models. Using a rank K approximation, what would otherwise be N^3 operations like solves or matrix products become K^2 * N, which is manageable for K of order 10 or less.

Mathematical motivation

The geometric average of the covariance of the draws and the inverse covariance of the scores minimizes the Fisher divergence from a standard normal distribution to the preconditioned target density. Fisher divergence is like Kullback-Leibler divergence, but instead measures the expected distance between gradients rather than the difference between log densities. This turns out in practice to lead to a better diagonal preconditioner. The paper provides proofs and an appendix with geometric motivations.

What we do not prove and what I do not believe anyone understands is how to select an optimal static preconditioner when the target log density has varying curvature, in any form: diagonal, low rank plus diagonal, or dense. I think what we really need is dynamic adaptation that isn’t cubic like Riemannian and doesn’t require implicit integrators.

Background on Stan Playground is here.

And now Brian Ward reports:

There is now a new tab in the data generation section which lets you add files like csvs that can be accessed by the data scripts

image1719×773 47.4 KB

 

Generating Stan data and code with brms

You can now use BRMS to generate the Stan code and data for use in Stan Playground

 

image1228×377 26.5 KB

 

Quick share option

The export window now has an option to create a sharable link where all the data for the project is encoded in the URL. This can lead to some very long URLs, but it’s a very quick way to share your models without needing to authenticate for Github Gists.

JOSS Paper

Our paper on Stan Playground was published in JOSS this morning! Journal of Open Source Software: stan-playground: Run Stan models directly in your browser

How cool is that??

This post is from Bob.

Cole Monnahan (who came up with the idea and did the heavy lifting), two of his colleagues from the Template Model Builder (TMB) project, Kasper Kristensen and James T. Thorson, and I just put a paper up on arXiv.

• Monnahan et al. 2026. Leveraging sparsity to improve no-U-turn sampling efficiency for hierarchical Bayesian models. arXiv:2603.02437
   
• R package on GitHub (MIT licensed): SparseNUTS

If you have any feedback on either, we’d love to hear it.

The method

It’s quite simple if you know how the components work. The tricky part is getting them all to play together nicely, efficiently, and robustly.

1. Take a max marginal mode of a hierarchical model, of the kind produced by the linear mixed effects package (lme4) or Template Model Builder (TMB). Both packages are implemented in C++ under the hood and distributed with an R interface.

2. Center a second-order Taylor series approximation at the marginal mode using a precision matrix (inverse covariance) rather than covariance matrix. This is like a Laplace approximation, but it’s not centered at a global mode.

3. Take the resulting sparse precision matrix and use it to precondition the target density. This is equivalent to using the precision as a mass matrix in HMC, as Radford Neal showed in his introduction in the MCMC Handbook. This approach is necessitated by Stan’s lack of sparse mass matrix support.

4. Use TMB’s R interface and Andrew Johnson’s StanEstimators package to make the model available to Stan’s samplers and other tooling.

There’s nothing Stan-specific about this technique. It could be rolled into PyMC and NumPyro with a lot of work if JAX’s experimental sparse library is up to the task of duplicating what TMB is doing.

Empirical evaluation

The paper is evaluated on 15+ realistic models, some of which can be scaled. The results show that with this preconditioning, Stan can scale to 10K+ parameter hierarchical models when there is sparseness and high correlations.

The paper demonstrate how much better SparseNUTS is than trying to use Stan’s built-in diagonal or dense mass matrix estimators. Obviously, the dense mass matrix won’t scale to 10K parameters as the mass matrix will have 100M entries and that becomes nearly impossible to estimate or manage computationally in Stan. Stan is very slow to estimate mass matrices in these cases (WALNUTS should be better, which we hope to roll out soon), but the real problem is that diagonal preconditioning is insufficient for hierarchical models.

Should you use it?

Yes! The form of hierarchical models addressed by SparseNUTS are widely used and Stan by itself simply cannot fit them. SparseNUTS also provides much more flexibility than INLA in writing custom likelihoods and priors.

Why not code this method in Stan?

We do not have the sparse Hessian tooling within Stan’s automatic differentiation library to implement this directly. Nor are we likely to roll out something for specific model classes, as we try to stick to black-box techniques in the core of Stan. This isn’t because we don’t like model-class specific methods. We are happy to implement them in special packages like brms. It’s just because one has to limit one’s scope somehow to keep a software project manageable with limited developer hours.

Some alternatives on the horizon

We are evaluating how far we’ll be able to scale Charles Margossian’s black-box approach to marginalization, which Steve Bronder and Charles have gotten into shape so it will be in the next Stan release. As INLA demonstrates, marginalization can be an efficient alternative to sampling the whole model if you can get the marginalization right.

We’re also evaluating how far we can get with Adrian Seyboldt’s Nutpie sampler’s approach to estimating low-rank plus diagonal mass matrices based on Fisher divergence, which could also potentially solve this problem as most of the structure in hierarchical models is low rank. Adrian and Eliot Carlsen (both of PyMC Labs) and I are writing that up now and we will release the paper soon. So far, Nutpie’s low-rank plus diagonal adaptation works incredibly well for some problems, but fails spectacularly on others for reasons we do not yet understand.

David Zoltowski, Skyler Wu, Xavier Gonzalez, Leo Kozachkov, and Scott Linderman write:

We propose algorithms to evaluate MCMC samplers in parallel across the chain length. To do this, we build on recent methods for parallel evaluation of nonlinear recursions that formulate the state sequence as a solution to a fixed-point problem and solve for the fixed-point using a parallel form of Newton’s method. We show how this approach can be used to parallelize Gibbs, Metropolis-adjusted Langevin, and Hamiltonian Monte Carlo sampling across the sequence length.

Bob Carpenter presented their paper at the computation seminar last week. Bob was enthusiastic about this work, and his enthusiasm was contagious.

The above graph shows what’s happening.

“Sequential Sampling” is when they run HMC on this example. The gray line shows the trajectory or sequence of simulation draws. (Here they just do one chain to demonstrate the idea, but later in the paper they talk about the benefits of simulating multiple chains.) The trajectory is a random vector, but it can be thought of as deterministic, conditional on the sequence of random numbers that is used to do the jumping and make the accept/reject decisions.

Here’s the idea: Given the input sequence of random numbers, the goal is to compute that trajectory. But rather than doing it the usual way, by running HMC, they do it by starting with an initial guess–here, the flat or constant trajectory shown as the horizontal line of the top graph above in the “Parallel Sampling” column–and then moving the points around until they reach the correct trajectory. This can be done because the correct trajectory is Markovian (in this case, first-order Markov), and the value at each time point from t=1 through t=1000 is a known deterministic function of the value at time t-1 and the set of input random numbers corresponding to that iteration. They write this as s_t = f_t(s_{t-1}), where f_t is that deterministic function, and the random numbers for time t are implicitly included in f_t.

Each iteration of the parallel algorithm is a step in which the values of the trajectory at each point are jiggled up and down in an attempt to make each s_t as close as possible to f_t(s_{t-1}), for each t. The beauty of this algorithm is that these updates can be done locally, so you can do them in parallel. It turns out that this can be done in O(log(T)) ticks of wall time. The update requires O(T) computations but with enough processors the parallelism saves lots of time.

When does this make sense? I know next to nothing about parallel computing, but apparently it can be done on GPUs, and there are now programming languages that can distribute the computing steps.

There are lots of details that make this work, including using autodiff on the functions f_t.

Discarding bad approximations in the middle of the chain?

Here’s something I noticed about their example.

Look at Figure B above: it shows that after 147 iterations of the parallel updating, the computation has reached the correct values of the chain for this set of input random numbers. But what about after 10 iterations? It’s pretty damn close already! There seem to just be 4 or 5 points that are in the wrong place.

So here’s my question. Could there be a way, if you want to stop early, to identify those bad points and throw them out?

Usually when we’re doing HMC, you have to keep the whole chain or throw it all out, or maybe you keep running it and you throw away the initial part that hasn’t completely warmed up. But in this new algorithm, the entire chain is a static target, so it seems that it might be possible to just throw away a few bad iterations in the middle.

At first this idea might seem hopeless: if you don’t know the final answer, how would you know which points to throw out? But we can compare each s_t with f_t(s_{t-1}), so maybe just take the points where the discrepancy is large and toss them?

As a start, I’d be interested in how things work with the example shown above. Do the bad points in parallel iteration 10 correspond to large values of |s_t – f_t(s_{t-1})|? If so, maybe there’s something we can do here! Nothing stopping us from continuing to run all the way to parallel iteration 147 or whatever, but in the meantime maybe we can do better along the way.

It goes like this. We start with a simple linear regression. In the Stan code we first simulate data given the true parameter values and then express the model:

data {
  int N;
  real a_true, b_true, sigma_true;
}
transformed data {
  vector[N] x;
  vector[N] y;
   for (n in 1:N) {
    x[n] = uniform_rng(0, 50);
    y[n] = normal_rng(a_true + b_true * x[n], sigma_true);
  }
}
parameters {
  real a, b
  real<lower=0> sigma;
}
model {
  for (n in 1:N) {
    y[n] ~ normal(a + b * x[n], sigma);
  }
}

And here are the input data:

{
  "N": 100,
  "a_true": 20,
  "b_true": 0.8,
  "sigma_true": 15
}

We compile and run it in Stan Playground, setting the seed to 123, and here are the results:

This all seems reasonable. The true values of a, b, and sigma are 20, 0.8, and 15, and the estimates are 19.0 +/- 2.9, 0.85 +/- 0.10, and 15.5 +/- 1.2.

Now let’s try throwing in a prior. I’m gonna use my default prior for parameters on unit scale:

  a ~ normal(0, 1);
  b ~ normal(0, 1);
  sigma ~ exponential(1);

This goes into the model block of the above Stan program.

Before going on, let me point out that the parameters of this problem are not on unit scale. By “unit scale,” I mean, “being drawn from a distribution with scale that is of order 1,” and a_true = 20 here. Not unit scale at all! Also, sigma_true = 15, but this will be less of an issue because the exponential prior is much weaker than the normal in the tail.

So, yeah, this is a terrible prior. But, y’know, that happens. Not on purpose (usually), but we’re all busy and sometimes we apply models that don’t make sense or don’t fit the data, or both.

So let’s see what happens with this strong prior.

Actually, we kind of know what to expect. The likelihood is centered at (a, b, c) = (19.0, 0.85, 15.5), and the prior mode is at (0, 0, 0), so we expect the parameter estimates all to be pulled toward zero. Especially the intercept, a, since the prior for a is much stronger than the likelihood. And it’s a linear normal model so the posterior should be pulled smoothly toward the prior, not the sort of funny stuff that would happen with a Cauchy prior, say.

But let’s run it and see what happens. Here goes:

Let’s look at each parameter:
– The estimate for “a” has gone from 19.0 to 1.7: So, yeah, the prior was strong, it pulled the estimate almost all the way to the prior mode of 0.
– The estimate for “b” has gone from 0.85 to 1.4: Wha….?? It got pulled away from the prior!
– The estimate for “sigma” has gone from 15.5 to 16.5: Again, it got pulled in the wrong direction!

What’s going on here?

The key to understanding this problem is to realize that, although the three parameters a, b, and sigma are independent in the prior distribution here, they’re strongly dependent in the likelihood.

So let’s talk through what’s happening.

First, the prior for “a” is very strong, so the intercept is pulled almost all the way to 0.

Second, look at the data. The predictor x ranges from 0 to 50, and the true regression line is y = 20 + 0.8x, so if you grab the line at x=0 and pull it down from y=20 to y=0, while still trying to go through the data, you’ll increase the slope. Just draw the scatterplot and you’ll see what I mean.

So the strong prior on “a” causes the estimated slope to go up. Also there’s a prior on “b” that pulls the slope toward 0, but it’s a weak prior relative to the likelihood (recall that the likelihood-based estimate of b had a standard error of 0.1), so it doesn’t have much effect. The net effect of the joint prior is to increase the slope of the fitted line, thus pull the estimate of b away from 0 in this case.

Finally, the new line doesn’t fit the data so well–that’s what happens when you move away from the maximum likelihood estimate, the fit gets worse!–and so the estimate of sigma goes up. Again, there is a prior on sigma that pulls its estimate toward zero, but it’s not as strong as the signal from the data telling us that the fit is getting worse.

We can see this even more clearly by removing the priors on b and sigma, so that all we have is the prior on the intercept, a:

data {
  int N;
  real a_true, b_true, sigma_true;
}
transformed data {
  vector[N] x;
  vector[N] y;
   for (n in 1:N) {
    x[n] = uniform_rng(0, 50);
    y[n] = normal_rng(a_true + b_true * x[n], sigma_true);
  }
}
parameters {
  real a, b;
  real<lower=0> sigma;
}
model {
  a ~ normal(0, 1);
 // b ~ normal(0, 1);
 // sigma ~ exponential(1);
  for (n in 1:N) {
    y[n] ~ normal(a + b * x[n], sigma);
  }
}

And here’s what happens:

So, yeah, the inference for “b” is about the same, “sigma” gets pulled up even more, and “a” gets pulled even more toward zero. That last bit makes sense too: in the model with informative priors on all three parameters, the priors on b and sigma have the effect of somewhat decreasing the amount that the intercept can be pulled down to zero.

Lessons from this example

First, put your parameters on unit scale so that you won’t accidentally use priors that are so contrary to the data. At this point you might say that you don’t know anything about your parameters so you can’t put them on unit scale. But I disagree! In every application I’ve seen, we have some sense of the scale of the effects and adjustments in the model. If you don’t have this information–if you’re flying blind, as happens when you’re developing default software–then there should be a way to scale the parameters based on the data. That’s how we do things in rstanarm, as discussed in Section 9.5 of Regression and Other Stories.

Second, plot the data and fitted model. Check the fit of the model to the data! You’re under no obligation to stick with a model that clearly doesn’t fit.

Third, experiment! As shown above, you can learn a lot by perturbing the model and re-fitting. No need to just stare at your code and your computer output to try to figure out what is happening. You can and should use experimentation as a regular tool.

Nikolas Siccha writes:

I got inspired by Bob’s latest blog post to bring your attention to this code generation tool. It’s currently supremely buggy and only really usable by myself – but it may still be interesting for you to see, as it IMO solves many of Stan’s “problems.”

“Supremely buggy and only really usable by myself” . . . that describes lots of the software that I write!

Follow the link and here’s what you see:

Brings Julia syntax to Stan models by implementing a (limited) Julia to Stan transpilation with many caveats. See test/slic.jl for implementations of a few simple posteriordb models and see src/slic_stan/builtin.jlfor a list of built-in functions and examples of user defined functions.

Current features include

• activity analysis (automatically determines what is data, transformed_data, parameters, transformed parameters, model, or generated quantities),
• automatically inferred types, shapes and constraints – including for user defined functions (including the function arguments, function body, and function return type),
• automatic posterior pointwise likelihood and predictive generation,
• (variadic) user defined functions,
• higher order (user defined) functions (such as map, broadcasted, sum and more),
• sub models,
• post-hoc model adjustment,
• named tuples,
• (approximate) automatic code formatting à la Blue,
• and more.

Upcoming features include, in order of priority and estimated arrival,

• easy runtime assertions (like Julia’s @assert) – and support for other macros,
• model docstrings,
• custom types (for method dispatch – this would help with more “Julia-style” broadcasting, e.g. via Ref),
• closures via Julia’s Do-Block Syntax (to make within chain parallelization via reduce_sum less painful),
• lower transpilation runtimes (currently, transpilation can sometimes take longer than compilation – there is currently at least one algorithmic inefficiency on top of the systemic implementation inefficiency),
• a much better user experience,
• more and better tests,
• keyword arguments,
• default arguments,
• inlining (to reduce potential runtime overhead),
• easier custom parameter transformations (going from sampler parametrization to user parametrization – aka as constraining),
• array comprehensions,
• a more complete (and more correct) coverage of built-in Stan functions,
• better name resolution (currently user defined functions or sub models have to be defined in Main),
• functions that mutate their arguments (solved via inlining),
• and more.

Almost anything that’s possible in Julia should be possible to be transpiled to Stan. Of course, unless Stan is much faster than Julia (+Mooncake or Enzyme) for the model in question, just sticking to Julia comes with many advantages.

Features which I am on the fence about, but currently not planning to implement:

• a Julia backend,
• target += statements,
• top level control flow,
• top level mutability,
• getting rid of superfluous parentheses.

Features which are NOT planned:

• (automatically) transpiling Julia functions which have not been defined via @deffun.

The earn_height.stan model below becomes

I don’t quite get that code at the end, but then again I don’t speak Julia.

One reason I think this is cool is that Julia and Stan started at the same time! Both projects were originally supported by a Department of Energy grant that funded Alan Edelman and me to hire postdocs and programmers. Alan’s group built Julia; we built Stan. Thanks, Bob, Matt, Daniel, and everybody else who did this!

Stepping back, it’s good to be able to go between languages and not be stuck in just one framework. We wrote Stan for two reasons:

1. To be able to easily and flexibly express Bayesian models and fit them to data.

2. To be able to perform statistical workflow, which includes fitting multiple models, postprocessing inference, plotting data and fitted models, etc.

For step 2, it was absolutely necessary for Stan to link smoothly to general-purpose computing environments such as R and Python, which is why we built rstan and pystan (which have since mostly been replaced by the lighter interfaces cmdstanr and cmdstanpy; thanks, Jonah, Mitzi, and everybody else who did this!). You can also run Stan from Julia. What Niko is doing here is going beyond that in allowing some sort of improved version of Stan that can be written in Julia, which seems like it could be really useful.

Beyond all this, as a scholar I’ve always tried to avoid any territorial attitude. I’ve just seen this too often, with researchers who have some method or another and are out there dissing the competition. OK, don’t get me wrong here, I’m happy to diss, when it’s appropriate, but not just cos it’s competition. Science isn’t a competition; we’re all in this together.

So, yeah, Stan is great. I love Stan, and I use it. I eat my own dogfood. But the purpose of Stan is Bayesian workflow. And if Stan can be helpful to Bayesian workflow in other ways, that’s great.

When building models, it’s helpful to check our understanding by simulating fake data and seeing if the fitted model can recover the true underlying parameters. And it’s so easy to do in Stan!

See above–I just did a simple example in Stan Playground. (Thanks, Brian!) The code’s kind of ugly because I specified the true parameter values in the transformed data block Stan code rather than calling them as data; on the upside the whole thing is then in one program.

This is not the same as full simulation-based calibration checking (SBC). In SBC you draw the true parameter values from the prior, repeat the entire process in parallel many times, and then check the average coverage of posterior inferences with respect to the true parameter values, averaging over the prior. Here I’m just running one simulation and setting the true parameter values just once. It’s a kind of quick-and-dirty SBC which can still be useful in catching problems such as nonidentified models or poor mixing.

And it’s so easy to do! Whenever you fit a model, you should be checking it on fake data.

I was doing the above example because I wanted to quickly check for one of the exercises I’m writing for the forthcoming Bayesian Workflow book.

P.S. Here’s the code for the above example if you want to try it yourself:
Continue reading →

This post is from Bob.

Drift-diffusion models

Whew. The cdf function for the seven-parameter drift-diffusion model was just merged. The pdf was merged a few months ago. This is a big deal. These pdfs and cdfs are used for in decision-time models in cognitive psychology. There’s a really nice illustration through NLM on nih.gov. The basic idea is that you have a binary task like deciding if an image is red or blue. The data being recorded is time to decision and the decision being made. The underlying generative model is a continuous Wiener diffusion process that has a lag time to get started before drifting with some bias toward opposing decision boundaries. The decision is determined by when it crosses a boundary and which one it crosses (see the illustration). The cdf is important when the task ends before a decision is made, giving you censored observations, which require cdfs or truncated pdfs to implement.

The first time I saw this model being applied was by Bruno Nicenboim and Shravan Vasishth (psycholinguists at Potsdam at the time, though Bruno has since moved to Tilburg) about six or eight years ago. At that point, it took Stan a month or so to fit the model (yes, that’s a month, not a typo)—you may know them as two of the three authors of the really wonderful book, Introduction to Bayesian Data Analysis for Cognitive Science (2025, CRC), which, in its final chapter, covers accumulator models of which the drift-diffusion model is one form. Now these models are very fast in Stan with the new built-in functions.

The pull requests and engineering challenge

Hats off to Franziska Henrich, a cognitive psychologist and Stan developer at the University of Freiburg, aka GitHub user Franzi2114, for writing the code and bearing with Steve Bronder’s hundreds of comments and fixes and my final round of a hundred or so change requests. You can see all the gory details in the discussions around the pull requests and in the code itself.

• GitHub: pdf pull request
• GitHub: cdf pull request

These pair of functions were perhaps the two hardest functions to get into Stan for a myriad of reasons. The most challenging obstacle beyond the inherent complication of the functions themselves is that our testing framework for densities can’t handle seven-parameter densities. So all the tests had to be projected into subsets of parameters (and seven choose four or five or whatever it was led to a lot of tests). A further difficulty is that to make the arithmetic stable, the code branches all over the place (see the Hartmann and Klaeur article linked below), which also complicates testing.

Some academic background

In addition to Vasishth et al.’s book chapter, there is a vast literature on drift-diffusion models in cognitive psychology and elsewhere. Most relevantly, Franziska wrote an open-access article about the model and the Stan implementation.

• Henrich, F., Hartmann, R., Pratz, V., Voss, A., & Klauer, K. C. 2024. The Seven-parameter Diffusion Model: an Implementation in Stan for Bayesian Analyses. Behavior Research Methods.

Luckily, Hartmann and Klauer provided the derivatives in a previous (closed-access) article.

• Raphael Hartmann and Karl Christoph Klauer. 2021. Partial derivatives for the first-passage time distribution in Wiener diffusion models. Journal of Mathematical Psychology.

As is often the case, you can find a pdf through Google Scholar. It’s a daunting pile of mathematics that puts the “M” in “mathematical psychology.” Luckily for us, the authors published an R implementation in package WienR on CRAN (the name is because it’s the Wiener diffusion model underlying the process), which Franziska could use for testing.

Coming to the Stan language next release

We just put out a new Stan release, so we have plenty of time to get the language wrappers around the math library functions before the next release of Stan. Ideally, we’ll also have a User’s Guide chapter with examples of how to use them. We’re always open to new *User’s Guide* chapters about models or methodologies in wide use, and as you can see from this example, we take pull requests, which go down much more easily with the User’s Guide.