Category Archives: Stan

A slew of improvements to NUTS

Posted on by
6

This post is from Bob.

Hold onto your hats, because 2026 promises to bring a whole slew of improvements to MCMC for continuously differentiable densities. The really awesome part is that all of these improvements are orthogonal, so they stack. I’m going to list the ones we know work first, followed by a couple in which we have high hopes.

I’m currently working on implementing all this in a C++ reference implementation, which I plan to roll out with an interface like that of Nutpie. You can follow here:

You’ll find the work in progress on branches. We are, of course, always happy to get feedback on the code, and I’m happy to get feedback on these or other ideas for sampling.

Fisher divergence for mass matrix adaptation

Adrian Seyboldt developed much faster and more robust and better targeted adaptation using Fisher divergence for his sampler Nutpie. I’m helping Adrian and Eliot Carlsen finish up a paper on this, which we hope to release in a week or so. In addition to better diagonal and dense estimators, it contains a really nice low-rank plus diagonal preconditioner that seems very effective (the risk is getting stuck too much in the subspace defined by the low rank structure). It also contains all the mathematical proofs, thanks to Adrian. The basic idea is that Fisher divergence adaptation targets getting the preconditioned target as close to a standard normal in terms of gradients as possible (not as close as possible in terms of density—that’s KL divergence). With some surprising mathematical magic, you can solve the exact optimization of Fisher divergence by taking the geometric mean in the affine-invariant manifold of positive-definite matrices of two quantities: the inverse covariance of the draws and the covariance of the scores (where the score is the gradient of the log density). In a multivariate normal, the covariance of the scores is the inverse covariance of the target density, has zero expectation, and thus acts like a control variate. Who knew? (That was a rhetorical question. Ben Goodrich, of course, knew—he knows everything.)

Adaptive step size on the fly

This was developed by Nawaf Bou-Rabee, Tore Kleppe, Sifan Liu, and me over the course of a few papers expanding on our Gibbs self tuning (GIST) ideas, culminating in our WALNUTS paper (on arXiv). The basic idea here is at each leapfrog step to try to take a step, and if the Hamiltonian diverges past a tolerance, try a smaller step size (half the original). There’s a bit of adjustment to do for reversibility (very much like in the MALT sampler of Lionel Riou-Durand et al. [the paper has an all star cast]), but the expectation is that the adjustment will be close to unity and that seems to hold in practice. There’s some subtlety here in that the WALNUTS sampler we’ve produced can be slower than NUTS (on a per gradient basis) for densities that NUTS can fit well, but if it’s well tuned in terms of minimal numbers of micro steps per macro step in the WALNUTS algorithm, it will be more efficient on a per-iteration basis. This shows we really are fixing the integrator, even if that doesn’t always give us a better sampler on a per-gradient basis. This same kind of issue comes up with higher-order leapfrog algorithms—they’re more precise, but often not worth the effort. On the plus side, when you have highly varying scales, like in the funnel density, WALNUTS can sample it effectively where NUTS will just fail silently reporting overly optimistic R-hat and ESS values.

Isokinetic sampler

As the name implies, the idea of isokinetic sampler is that you keep the kinetic energy constant. This is carried out by designing a kinetic energy function that’s different than the usual quadratic model derived from Newtonian physics (decomposed a la Hamilton, of course). This was developed decades ago by Mark Tuckerman at NYU’s Courant Institute for molecular dynamics, written about in Leimkuhler and Matthews Molecular Dynamics textbook, and reinvented recently by Jakob Robink and Uroš Seljak under the name “microcanonical sampling.” Isokinetic sampling has the remarkable property that you can fix a single energy level set and the isokinetic sampler can be restricted to that, yet remain ergodic for the distribution. Although they’re singing the praises of unadjusted methods, Reuben Cohn-Gordon joined Jakob and Uroš and they wrote a joint paper on Metropolis-adjusted methods (on arXiv). Of course, you could extend this to something like jittered HMC, multinomial HMC, or NUTS (we discuss these alternatives for standard HMC in the various GIST papers in detail). Recently my colleagues Tore Kleppe and Nawaf Bou-Rabee (my partners in crime on GIST) have verified the results presented by Robnik et al. and it seems to give a pretty clean factor of 1.5 to 3 over the standard NUTS kinetic energy model.

Unadjusted sampling, at least for warmup

In statistics, we’re used to thinking of unbiased MCMC with the correct ergodicity properties—you run them longer and you get a better answer. But in estimation in stats, we’re often OK with introducing a bit of bias if it gives us a large enough reduction in variance that we get better answers. If you view sampling the same way, then removing the Metropolis step from HMC gives you a biased sampler that can be a lot faster at moving toward the typical set where we want to sample than an adjusted sampler. And while Jakob, Uroš and Reuben are mainly working on unadjusted samplers and trying to prove error bounds for them under pretty strong assumptions, they have also suggested the reasonable tactic of using unadjusted sampling during warmup then switching to adjusted later so that everything works as expected for longer runs.

Smoother adaptation of Nutpie

Nutpie, like NUTS before it, works in blocks. It evaluates a given number of MCMC iterations, then updates its mass matrix estimate. I’ve developed a smooth alternative that simply exponentially discounts the past at a decreasing rate to mimic the block structure of Stan’s warmup. This seems to work very well and doesn’t suffer the problem of Nutpie of never converging due to a finite cap on block size (Adrian’s working on changing the version in Nutpie in some way to get around the initial design of 100-length blocks; Stan works in a sequence of blocks that doubles in size).

Adam replacement for dual averaging

Stan used the dual averaging stochastic gradient descent algorithm to set step size. It’s not spelled out in the short NUTS paper that dual averaging is an SGD algorithm because the NUTS paper never gives you the objective and its gradient. If you work it out backwards from the dual averaging algorithm, you can deduce that Hoffman and Gelman used a normal target (i.e., squared error), the gradient of which is the negative difference between the observation and the target. This gives you simple stochastic gradients you can plug into any SGD optimizer. I’ve just done this this week, but I’ve already found that Adam is much faster to converge, less highly variable after convergence, and more monotonic getting to convergence in the short tests I’ve done so far. I had to add an update discounting factor to Adam or you get the terrible oscillations rather than convergence for which both dual averaging and Adam are known.

Concurrent adaptation and convergence monitoring

This is a big one that Andrew’s been asking for for years. I’ve finished the online R-hat monitor (it’s in a branch of that name in the Walnuts repo on flatironinstitute). It uses the C++11 threads library to build an asynchronous, non-blocking monitor. The chains roll along as usual, each within its own thread, but they accumulate their within chain means and variances via Welford’s algorithm. They publish these in a per-chain Atomic using relaxed memory guarantees (hence the lack of blocking). Then there is also a monitor thread running that continually reads the per-chain means and variances and computes R-hat (the original one, not the split or ranked versions that are currently used in posterior (R) and ArviZ (Python)). I am able to run at least 100 R-hat checks per second without slowing down the chains noticeably, so it can detect convergence within a few iterations of when it happens even for very fast models. This all works blazingly fast on my new Mac Studio Server running 16 chains in parallel(*). Next up, I have to monitor adaptation. For that, there’s not a pre-built solution, but I’ll be using something like R-hat to monitor whether the mass matrix and step sizes have converged (on a log scale, which makes the monitoring respect the positive-definite manifold structure and operate scale free). And I’ll have to use a double-buffered store rather than an atomic because the vectors required for monitoring convergence of the mass matrix are D-dimensional, not default copyable.

(*) Mac ARM hardware

Not an algorithmic change, but … if you haven’t gotten the memo yet, the new Mac ARM chips are very well suited to parallel sampling like we do in Stan. Their memory architecture is much more tightly integrated into the CPU and more multiplexed than Intel architecture. This is great for sampling or other tasks where data and parameters are being hit asynchronously in memory in parallel.

I just got a Mac Studio with 20 performance cores on the M3 Ultra chip and 256GB of memory. I would highly recommend this specific machine if you have US$6K to spend (your cost may vary—electronics seem to be way more expensive outside of the U.S.). They also have a US$2K starter model, which should also be great compared to just about anything else. Even the MacBooks are good—the inexpensive Air I got a few years ago crushed my mega-expensive (albeit 8 year old) iMac running Intel Xeon chips. People have told me Stan’s broken on Windows after comparing Windows and Mac, but really, it’s just the memory architecture. I didn’t do any controlled tests, but I swear that tests that were taking 3s now take 1-2s after upgrading to the Tahoe OS (bunch of small UI changes that make pretty much no difference to my experience); Apple says they’re continuing to integrate more ARM goodness into their OS, so maybe that was it? The release notes don’t mention anything about heavy thread performance.

I don’t know how long this is going to be useful. I plan to develop a lot of these algorithms in JAX to run on GPU, which is going to be way faster than anything you can do on the desktop.

TENTATIVE IMPROVEMENTS

These haven’t been well tested by us yet.

Generalized HMC

NUTS is not good for GPUs. Matt Hoffman et al. have a new paper out that’s going into the next edition of the MCMC handbook about MCMC on modern hardware that explains why (on arXiv for now). Hugh Dance et al. just wrote another paper showing how to code something like NUTS on GPU, but it’s still expensive. Something like generalized HMC promises to give much of the advantage of NUTS implicitly without the elaborate recursive doubling structure. This is the kind of sampling to which Matt Hoffman and his former colleague Pavel Sountsov developed (Matt left Google and doesn’t work on sampling any more as far as I know). This is also what Gilad Turok, Chirag Modi, and I found with delayed rejection. It also sidesteps the problem of having to chain all the adjustments for local step-size adaptativity which can add up in a bad way for longer chains. It does add some additional tuning, which is how much to partially refresh momentum, which is a kind of proxy for path length for U-turns. It also has the advantage of being implicitly Rao-Blackwellized compared to HMC (you can average all the steps on an HMC path when computing expectations, but it’s just usually not worth the storage—this may turn out to be the case for G-HMC too).

Local mass matrix adaptation

I think of this as the final frontier. If we could locally condition, we wouldn’t need variable step sizes. We don’t have anything that works for this in complex cases. In Stan, you get a choice of a diagonal or dense preconditioned, and in Nutpie you also get low rank plus diagonal, but they’re all global preconditions, not local.

In very simple cases, we can use GIST to generate a mass matrix by taking an inverse Wishart sample around the negative Hessian. If you set the degrees of freedom correctly to get low variance, this pretty much perfectly preconditions a multivariate normal and should be a sound strategy in any log concave density. It’s just expensive, even with everything Cholesky factored. But that’s not enough for a system like Stan. The problem in going to more general models is that the Hessian’s no longer guaranteed to be positive definite (as it would be in a log concave density). This means we have to use something like Michael Betancourt’s softabs technique to condition the Hessian back to positive definite (e.g., eigendecompose, move negative eigenvalues up to positive, put back together, which is prohibitive because of the cubic cost). As an alternative, Nawaf and Tore have been looking at some explicit Riemannian integrators that only require a few Hessian-vector products, which are cheap with autodiff (linear in dimension rather than quadratic). This avoids the problem with the implicit integrator in Riemannian HMC, which is an additional obstacle beyond the need for positive-definite metrics. We might go back to implicit midpoint, as Arya Pourzanjani explored in his thesis, because we can use GIST to set a step size where we can guarantee stability; even so, the implicit nature of the algorithm is a real challenge for reversibility.

StanCon 2026 registration and abstract submission are now open

Posted on by
2

Registration for StanCon 2026, 17-21 August, in Uppsala, Sweden, is now open! You can already register and submit abstracts for contributions. Our first keynote speaker will be announced soon.

New this StanCon: In addition to tutorials, we are introducing workshops. These are structured similarly to workshops at other conferences. We will provide three-hour rooms on the Friday of StanCon 2026 for contributor-led mini-conferences to empower the community to shape the meeting by diving deep into emerging areas and shared interests.

Submit your abstract at https://www.stancon2026.org/abstracts/

Register for StanCon 2026 at Registration – Stan Conference 2026

Deadlines

  • 28 January: Workshop and tutorial proposal deadlines; early poster acceptance deadline (for participants needing early confirmation)
  • 25 February: Oral presentation abstracts deadline
  • 27 May: Final poster submission deadline
  • 10 June: Early-bird registration deadline

Please share this announcement with anyone who might be interested!

— The Local Organising Committee
Måns Magnusson, Sara Hamis, Aki Vehtari

From Bayesian inference to LLMs (Steve Bronder’s 2025 CppCon talk)

Posted on by
1

This post is from Bob.

Steve is now a C++ celebrity!

The production values for these are top notch. They’re also highly Google ranked, so you can’t miss them when you search for advanced C++ topics.

Steve’s been writing C++ code for Stan since before the pandemic. He now works with me at Flatiron Institute as a software engineer. He joined the Stan team after taking a GPU class and asking us if we had any use for GPUs. We said, “Yes, please.”

I say Steve’s a celebrity because of the prestige of CppCon—it’s where the language developers, top educators, and infrastructure tool builders hang out. For example, Steve had lunch with Matt Goldbolt. I use Goldbolt’s app all the time and didn’t realize it was a person’s name until Steve mentioned the lunch. To get some sense of CppCon, you can see Godbolt himself going full fanboy over a CppCon selfie with Laurie Kirk.

Effective sample size depends on the quantity

Posted on by
6

This post is by Aki.

I recently blogged about the term effective sample size and also commented “A further important point is that the effective sample size depends also on which expectation is estimated.” In the comments Visruth asked “Why does E(X^2) net a different ESS than E(X)?” and Kyurae Kim provided correct short answer. I decided to provide additional illustration

The No-U-Turn-Sampling (NUTS) variant of Hamiltonian Monte Carlo (HMC) aims to maximize the expected jump distance. The jump distance is not exactly maximized as there is randomness to keep the Markov chain reversible (which is a useful property to prove that the stationary distribution is the desired target distribution) and due to some algorithm efficiency choices. In some cases, this maximization of jump distance can lead to odd lag negative autocorrelations and higher effective sample size than the total number of draws. Let’s consider a case with theta being normally distributed. When the Markov chain is in a tail, a large jump distance will tend to take the Markov chain on the opposite side of the distribution, and if the next jump has also a large jump distance we go again to the opposite side and close to the first point. Then lag 1 and other odd lag autocorrelations are negative and lag 2 and other even lag correlations are positive. This will make effective sample size for E[theta] to be larger than the number of draws. Next it’s easiest to consider a normal distribution with mean 0 (vector of 0’s for multivariate normal) but this generalizes to non-zero means. If we now consider the absolute value of abs(theta) when we jump from tail to tail, the jump distance is likely to small and both odd and even lag autocorrelations are positive. For abs(theta) to maximize the jump distance, it would be better to jump between the tail and near mode, but NUTS is not designed to do that. theta^2 behaves as abs(theta).

We can test this with Stan and sampling from a multivariate normal. I sample from 16-dimensional unit normal, as I know that with 16 dimensions the algorithm details on how the Hamiltonian simulation is extended and U-turn decided, happen to be such that we get strong negative autocorrelations for theta.


data {
  int<lower=0> D; // number of dimensions
}
parameters {
  vector[D] theta;
}
model {
  theta ~ normal(0, 1);
}

And we run NUTS with CmdStanR

library(cmdstanr)
mod <- cmdstan_model("normal.stan")
fit <- mod$sample(data=list(D=16), refresh=0)

In the following, I examine just the first dimension of theta[1]

The autocorrelations for theta and theta^2 clearly show what I described above.

Negative autocorrelations lead to super-efficiency and the estimated ESS for theta is 10457, which is 2.6 times bigger than the total number of draws (4000). The estimated ESS for theta^2 is 1534, which is just 38% of the total number of draws! Even if we are not interested in theta^2 directly, we may be interested in sd of the posterior of theta, and computation of sd uses theta^2.

The following plot shows the ESS for the quantiles 0.05,…,0.95 (see Rank-normalization, folding, and localization: An improved Rhat for assessing convergence of MCMC for the details)

We see that if we want to report e.g. 90% posterior interval for theta, the accuracy of the interval end-points are based on about ESS of 2900. Thus, if we want to report posterior sd or posterior interval, there is not much benefit from having very high ESS for theta. In fact, as maximizing the jump distance in theta space is decreasing ESS for theta^2 and quantiles, we are wasting some computation time. However, multivariate isotropic normal is not that common posterior for most interesting models and data, and most of the time we don’t see negative autocorrelations and superefficient ESS, and thus there is no big need to change the NUTS algorithm related to this.

ESS is used to estimate Monte Carlo standard error and for any quantity of interest when MCSE is computed using, e.g., posterior package, the appropriate ESS is used. For convenience, posterior package report also Bulk-ESS and Tail-ESS which provide two summary values on the sampling efficiency in general (see the Rhat paper for the details). These are useful as quick summaries as they are scale free unlike MCSE which needs to be interpreted in the context of the scale of the quantity of interest (see more in the Digits case study).

EDIT: added 0 break for y axis

Predictive Modelling for Football Analytics is available!

Posted on by
5

This post is by Leo.
After a long and exciting journey, the book I co-authored with Dimitris Karlis, and Ioannis Ntzoufras Predictive Modelling for Football Analytics edited by CRC Press is available!

The book discusses the most well-known classical and Bayesian models, along with the main computational tools used in the football analytics domain. It also introduces the footBayes R package (built on Stan and CmdStan), which accompanies the reader through all the examples proposed in the book. It aims to be both a practical guide and a theoretical foundation for students, data scientists, sports analysts, and football professionals who wish to understand and apply predictive modelling in a football context.

This text is primarily for senior undergraduates, graduate students, and academic researchers in mathematics, statistics, and computer science who are interested in learning about football analytics. For sure, we really enjoyed writing this book.

Here’s the table of contents:

Chapter 1 – A short introduction to football analytics

Chapter 2 – Methods, algorithms and computational tools

Chapter 3 – Tournament and game prediction via simulation

Chapter 4 – Implementation of basic models in R via footBayes

Chapter 5 – Additional statistical models for the scores

Chapter 6 – Modelling international matches: the Euro and World cup experience

Chapter 7 – Compare statistical models’ performance with the bookmakers

You can order the book here.

And you can download all the data and reproducible R code of the book here.

P.S. I still remember fitting my first Stan model on football data when I was a visiting scholar at Columbia, in the Department of Statistics, back in 2016. I was sitting in Andrew’s office, hoping for a decent fit, and discussing with Jonah how to improve it. Almost ten years later, we now have some proof that we can always improve our models;)

P.S.2 As the great coach José Mourinho once said, “He who knows only about football, knows nothing about football”. For me, football is simply a tool to make statistics more accessible and engaging especially for those outside the field.

StanCon 2026 in Uppsala, Sweden

Posted on by
5

StanCon 2026 will take place in Uppsala, Sweden, from August 17th to August 21st, 2026.

The conference brings together researchers and practitioners passionate about Bayesian inference and probabilistic programming in one of Sweden’s most historic and vibrant university cities. Attendees will enjoy a week of conference talks, workshops, and tutorials spanning both foundational methods and real-world applications.

Stay tuned for updates! More information about registration, abstract submission, the program, and deadlines will be posted at https://www.stancon2026.org/.

The conference is sponsored by Uppsala University, The Swedish Excellence Centre for Computational Social Science at Linköping University, eSSENCE, and Beijerstiftelsen.

With kind regards from the organising committee
Måns Magnusson (Uppsala University)
Sara Hamis (Uppsala University)
Aki Vehtari (Aalto University)

Aki looking for a doctoral student to develop Bayesian workflow

Posted on by
1

I (Aki) am looking for a doctoral student with Bayesian background to work on Bayesian workflow and cross-validation (see my publication list for my recent work) at Aalto University, Finland (the world’s happiest country). You would also collaborate with Andrew and Stan and ArviZ developers. You can apply through the ELLIS PhD program (dl October 31)

The Dodgers are hiring

Posted on by
Reply

Brendan Cooley writes:

I’m an analyst with the Los Angeles Dodgers. We are looking for budding Bayesian statisticians and deep learners that are interested in joining us for an internship next summer. We have a soft spot for jax and numpyro but Stan and PyMC folks are obviously always of interest. Job description is here.

Students can get in touch at [email protected] with any questions.

It looks like the Dodgers might be facing the Mets in the first round of the playoffs this year. So, if you do take this job, please do me a favor and hold off for a couple weeks on giving the Dodgers any useful ideas!

The Miami Marlins are hiring

Posted on by
Reply

This post is from Daniel.

Bryant Davis at the Miami Marlins forwarded me a job posting. They’re looking for someone Bayesian to join them!

The Miami Marlins Baseball Research and Baseball Solutions departments are seeking entry- and senior-level Data Scientists and Baseball Analysts. In particular, they are seeking individuals who’ve worked in probabilistic programming languages – anyone with experience in Stan or PyMC who has a passion for working in sports is encouraged to apply! For more details, check out the job posting here.

Show, don’t tell: ChatGPT 5 marginalizing Gelman’s measurement error model in Stan

Posted on by
48

This post is from Bob

Even though the post is from me, the story is partly about Andrew, but mostly about how LLMs are getting better at math. Rather than telling everyone about LLMs and having a theoretical discussion, I thought it’d be useful to show you another example of the kind of interaction with LLMs that I find super useful.

The backstory

Andrew frequently mentions that he’s impatient with Stan’s default sampling times and would like something to give him a rough answer faster. So I and several colleagues have been working on finding faster adaptation (for now, I would recommend Adrian Seyboldt’s fast adapting Nutpie sampler, which can already be run with models coded in Stan or PyMC). That means evaluation. So I finally asked Andrew for an example.

Gelman’s measurement error model

Here’s the measurement error model that Andrew gave me:

data {
  int<lower=0> N;
  vector[N] y;
  vector[N] x_star;
  real<lower=0> sigma_x_star;
}
parameters {
  real a, b, mu_x;
  real<lower=0> sigma, sigma_x;
  vector[N] x;
}
model {
  x ~ normal(mu_x, sigma_x);
  y ~ normal(a + b*x, sigma);
  x_star ~ normal(x, sigma_x_star);
}

As an aside, I love that there’s an appropriate lower bound on N.

This is a classic measurement error model—if you’re not familiar, there’s a short introduction to exactly this flavor of model in the Stan User’s Guide introduction to measurement error models.

Bob’s revised model

This looks like a BDA type model—there are improper flat priors. And the latent true values of the covariates x have a centered parameterization. Both of these can cause problem fitting, so I swithced to a non-centered parameterization using offset and multiplers, then added some weakly informative priors for all the parameters.

data {
  int<lower=0> N;
  vector[N] y;
  vector[N] x_star;
  real sigma_x_star;
}
parameters {
  real a, b, mu_x;
  real<lower=0> sigma, sigma_x;
  vector[N] x;
}
model {
  a ~ normal(0, 1);
  b ~ normal(0, 1);
  mu_x ~ normal(0, 1);
  sigma ~ lognormal(0, 0.5);
  sigma_x ~ lognormal(0, 0.5);
  
  x ~ normal(mu_x, sigma_x);
  y ~ normal(a + b * x, sigma);
  x_star ~ normal(x, sigma_x_star);
}

It fits better than Andrew’s but still not well. With the default 1000 warmup and 1000 sampling iterations on four chains, there were still R-hat values above 1.2. If we crank that up to an order of magnitude more iterations, everything looks OK, but it still has an integrated autocorrelation time of around 100 (meaning ESS is 1/100 times the sample size).

GPT 5’s revised model

Looking at the model, I thought like it looked like a prime candidate for marginalizing out the x parameter—it has a normal prior and it’s used in a normal likelihood. But I’m terrible at this kind of math and there weren’t any math stats folks around like Philip Greengard, so I turned to GPT 5, which had already done some math for me I couldn’t do myself (calculating the second and fourth moments of Neal’s funnel). Here’s my interaction with GPT 5.0 set to “auto”. GPT 3.5 could barely write a syntactically correct Stan program for linear regression, so let me just say it’s come a long way since then.

[Bob] I have the following Stan program for measurement error in x and I’m wondering if there is some way to marginalize out some of the variables to make this sample more efficiently. I figured since everything is normal there would be some chance of this. In particular, I would like to marginalize the x out if possible.

```stan
[code from above]
```

I would be OK slightly modifying the priors on the sigma to be inverse gammas on sigma^2 if that would help grease the wheels of the marginalization.

And lo and behold, the genie returned a correctly marginalized Stan program:

data {
  int<lower=0> N;
  vector[N] y;
  vector[N] x_star;
  real<lower=0> sigma_x_star;
}
parameters {
  real a, b, mu_x;
  real<lower=0> sigma, sigma_x;
}
transformed parameters {
  real inv_var_x  = inv_square(sigma_x);
  real inv_var_xs = inv_square(sigma_x_star);
  real tilde_v    = 1.0 / (inv_var_x + inv_var_xs);              // Var(x | x*)
  real<lower=0> sd_xstar = sqrt(square(sigma_x) + square(sigma_x_star));
  real<lower=0> sd_y_cond = sqrt(square(sigma) + square(b) * tilde_v);
  vector[N] tilde_mu = tilde_v * ( inv_var_x * rep_vector(mu_x, N)
                                 + inv_var_xs * x_star );
}
model {
  // hyperpriors (use whatever you prefer; half-normal/lognormal work fine)
  a ~ normal(0, 1);
  b ~ normal(0, 1);
  mu_x ~ normal(0, 1);
  sigma ~ lognormal(0, 0.5);
  sigma_x ~ lognormal(0, 0.5);

  // marginalized likelihood
  x_star ~ normal(mu_x, sd_xstar);                 // p(x*)
  y ~ normal(a + b * tilde_mu, sd_y_cond);         // p(y | x*)
}

Followed by this advice:

  • [GPT5] Integrating out μx as well is possible with a normal prior, but it induces dense cross-observation covariance (rank-1 update) and is typically worse for HMC/NUTS. Keeping μx as a parameter is usually better.

Did it work?

Yes! Not only did it work, it worked on the first try. And the results matched what I got from a much longer run of the full model, but with a great ESS.

GPT, can you explain to our users at home how you did that?

I was thinking this kind of operation would be useful to include in the efficiency chapter of the Stan User’s Guide. So I asked the genie to explain itself.

[Bob] This is great and it worked exactly as advertised and had an effective sample size for estimating parameter means of about half the total number of sampling iterations and agreed with a long run of the original model. Now I would like to write this up as an example for our users in the efficiency chapter of the [*Stan User’s Guide*](https://mc-stan.org/docs/stan-users-guide/efficiency-tuning.html). So I want an explanation pitched for a similar audience. Given that we’re going with the factored conditional form, could you go into a bit more detail about how you derived the forms of the four variables you have, the constants $\tilde{v}$ and $\mu_i$, and the distributions for $\x_i^\star$ and $y_i \mid x_i^\star$? You can break it down into known properties of normals, though please remind me of what these are because I’m not very good with classical math stats manipulations of exponential families, nor will many of our users be.

And, of course, with an LLM, Bob’s your uncle (from the French idiom et voilà):

Of course, I will have to double check all the math until I understand it thoroughly, then convert to the Stan User’s Guide style.

They’re looking for businesses that want to use their Bayesian inference software, I think?

Posted on by
1

Kirill Parinov writes:

I recently co-founded Boston Bayes, a Dutch venture studio, together with the founders who created RxInfer — a powerful open-source probabilistic programming framework enabling organizations to make data-driven decisions under uncertainty. https://lazydynamics.com/

Boston Bayes is now seeking independent founders (for CEO positions) to launch new ventures powered by this technology — startups that turn uncertainty into a competitive edge through real-time probabilistic inference.

Boston Bayes will provide the tools, deep tech, and ongoing support needed to build the future of probabilistic AI, while the founders maintain independence and ownership of their ventures. You can read more here: https://bostonbayes.com

Do you happen to know anyone in your network who might be interested in founding a spin-off with Boston Bayes’ support?

All of this is baffling to me, starting with that I don’t know what a “venture studio” is, also it’s named Boston but it’s in the Netherlands, also I don’t understand what it means to be seeking “independent founders for CEO positions,” also I’m not quite sure what is the real-time probabilistic inference they are talking about. Also I don’t get what’s up with RxInfer, but Bayesian inference is cool, and anything we put in Stan and our workflow book and our research articles is open-source, so anyone is free to use these ideas in whatever computer program they’re writing.

In any case, maybe someone reading this post is doing Bayesian workflow in the business world and would find this to be useful?

Uber could use your statistical analysis.

Posted on by
Reply

Johannes Hallermeier writes:

After CU, I went on to work at Uber, as an applied scientist on the policy research team.
I’m reaching out because Uber is hiring for a number of interesting applied science & engineering roles across policy, marketing and marketplace.
—————————————-
Applied science / engineering roles at Uber
– two roles on policy research (#1, #3)
– two roles on marketing applied science (#2, #4, plus a third role coming online soon)
– several roles on marketplace (pricing, matching, incentives, etc., see rest of list below)
– before applying, please ping [email protected], including your CV and preferred role (for questions and/or referrals)

Applied Scientist II, Policy & Consumer Research

Applied Scientist II, Brand Marketing

Applied Scientist II, Earnings Policy

Applied Scientist II, Marketing Applied Scientist II

Machine Learning Engineer

Optimization Engineer (Operations Research)

Backend Engineer, Rider Pricing & Incentives

Senior Optimization Engineer

Senior Software Engineer, Rider Pricing and Incentives

Senior Software Engineer, Dynamic Pricing

Staff Machine Learning Engineer, Causal Inference

Staff Software Engineer

Staff Software Engineer, Rider Pricing Platform

Staff Machine Learning Engineer, Pricing and Incentives

Staff Machine Learning Engineer, Dynamic Pricing

Staff Optimization Engineer, Dynamic Pricing

Senior Staff Machine Learning Engineer, Marketplace Pricing & Incentives

Science Manager – Dynamic Pricing

Manager, Science

Senior Staff Engineer – Marketplace Competitive Intelligence

Scientist II – Competitive Intelligence

Sr. Scientist – Competitive Intelligence

Sr. Scientist – Competitive Intelligence

Scientist II, Pricing and Incentives

Senior Scientist, Pricing and Incentives

I guess maybe if you interview for a job there, they’ll pick you up directly from the airport?

In any case, I can only assume that expertise in Bayes and Stan will come in handy.

(1) Fitting hierarchical models in genetics, (2) A Stan model that runs faster with 400,000 latent parameters, (3) Super-scalable penalized maximum likelihood inference for biome problems, (4) “In the end, I basically gave up working on biology because of the politics.”

Posted on by
4

Our recent post elicited this comment from Bob that was so detailed I thought it deserved its own post.

This one comment from Bob has enough material for 4 standalone posts:

1. Fitting hierarchical models in genetics: why the full Bayesian approach gives the right answer, whereas an existing shortcut method using approximations and likelihood ratio tests does not.

2. Stan’s HMC scales well in high dimension, but not so well in bad geometry, so for this problem it was better to use a Poisson/gamma model with 400,000 latent parameters than a much lower dimensional but computationally awkward negative binomial model. (Recall the dictum that mixture models all have computational problems and that all computational problems are essentially mixture models.)

3. “Super-scalable penalized maximum likelihood inference for biome problems with some collaborators in optimization. With 150K biomes averaging 3 megabases each and arbitrary amounts of data, it fits in about 20m on a couple decent servers (doesn’t need GPU).”

4. The struggles of interdisciplinary, or interlab, collaboration within certain subfields of biology.

Bob’s summary of that last point: “It’s trivial to rewrite a statistically sound, fully Bayesian version of DESeq2, but biologists are too conservative to use new tools, so I haven’t even tried to fix this. In the end, I basically gave up working on biology because of the politics.”

And here’s Bob’s comment in full:

There are already hierarchical models in wide use in genomics. For example, rMATS is a system for inferring splice variants. The statistical computation in that package is a mess, which led them to the ridiculous conclusion that fewer variants from a population with deeper sequencing is better for estimating the population values than more variants with less deep sequencing (basic stats will show you the latter is better if your goal is estimating the population). I rewrote it all as a Bayesian model. It gets the right answer statistically, namely that it’s better to have more samples at lower sequence depth, unlike rMATS. The problem with rMATS is their Laplace approximation for max marginal likelihood and then their likelihood ratio tests, which together give them the wrong results statistically, even for their model.

Negative binomials induce really bad geometry, as does just about anything that gives you multiple ways to explain the same data (i.e., either higher mean or higher dispersion can explain a large value). So I reparameterized in two ways. First, I coded the differential expression with a gamma-Poisson rather than negative binomial, which with 20 sequencing runs over 20K splice variants led to 400K latent parameters. Second, I use an overall mean plus difference parameterization on Gelman’s advice. Stan’s HMC scales well in high dimension, but not so well in bad geometry, so fitting is OK.

For political reasons, I haven’t been able to publish. As an example, my original collaborator’s Ph.D. supervisor literally won’t let her work on anything other than his project and the person whose lab she rotated through is only interested if it produces “breakthrough biology” we can publish in Nature, which is not something I can deliver. I haven’t been able to find anyone who can help me evaluate this and write it up in a way that the biologists would be OK with. I’d be happy to share the code. I even hired an intern specifically to work on this, but the intern showed up and refused to work on the problem!

More recently, I worked on super-scalable penalized maximum likelihood inference for biome problems with some collaborators in optimization. With 150K biomes averaging 3 megabases each and arbitrary amounts of data, it fits in about 20m on a couple decent servers (doesn’t need GPU). It’s like Rob Patro et al.’s salmon, but I figured out how to do all the bias adjustments for things like hexamer and positional bias (just push the uncertainty through Bayes style—not hard), which is why Rob told me that he abandoned the approach. I’ve had the same problem getting a biologist to help evaluate and write up. Even my collaborators on the autism project where I have a ton of gut biome data from Simons Foundation’s Autism Research Initiative (SFARI) and we did publish a paper, but they say they have zero time to do this. Everyone says they know someone, but when I contact them, they’re too busy. If I could come up with half a year funding for a postdoc, they’d loan me one (they really do assign their postdocs work—it’s not at all like CS and stats). It’s amazing to me how constrained biologists are in working with other people. The code we produced was really great, because Robert Gower (an optimization specialist) helped add a mirror descent optimization algorithm, which is about an order of magnitude faster than the L-BFGS default in Stan. And Robert Blackwell (it’s the 3 Roberts project), a high-performace compute specialist, added some high performance sparse matrix computation in C++ and figured out how to scale it out on a (CPU) cluster. What’s cool is that time doesn’t depend on the data size other than the really fast conversion to k-mer counts.

The autism gut biome paper for which I helped with the Stan modeling is Multi-level analysis of the gut–brain axis shows autism spectrum disorder-associated molecular and microbial profiles. This project was led by Jamie Morton, who was an amazing postdoc at Flatiron. Jamie wrote a really nice overview of how to think about compositional analyses like differential expression in Establishing microbial composition measurement standards with reference frames. But he has zero time to help build a state-of-the-art differential expression tool.

I would suggest looking at DESeq2 to see what a mess current differential expression inference is. It’s pretty easy to see why none of their p-values are well calibrated. It’s trivial to rewrite a statistically sound, fully Bayesian version of DESeq2, but biologists are too conservative to use new tools, so I haven’t even tried to fix this.

In the end, I basically gave up working on biology because of the politics.

Using hierarchical modeling to get more stable rankings of gene expression

Posted on by
33

Will Macnair writes:

I am a computational biologist (with a maths/stats background) working with genomic data. One task that comes up often is estimating changes in expression level between two conditions, for many thousands of genes, then ranking the genes in order of some measure of “confidence.”

My question is: In a Bayesian regression setting, what is a good way to rank the genes? If I calculate Bayes factors for “full model” vs “covariates-only model” and use these to rank the genes, is this ok? Are there better options?

To give a bit more detail, this task is typically referred to as “differential expression”, and we usually do this for all genes that pass some filter on very low / very boring expression. In the frequentist setting, the ranking is done with an adjusted p-value, and there are various popular methods for doing this (e.g. edgeR and DESeq2). For some recent analysis I have been using rstanarm, and I’ve been wondering what would be the best way to achieve a confidence-type ranking in this setting. The reason we do ranking is that we might do follow-up experiments for some genes, however the experiments are a lot of work. So we would like to prioritise genes where we have high confidence that something interesting is going on.

Some other thoughts:

– I searched previous blog entries and found this one. There is some nice discussion in the comments, in particular the thread started by Bob Carpenter where he talks about framing it as ranking. If I’ve missed something more recent, please point me to that!

– There are also some comments from Daniel Lakeland on how there has been a lot of work on the FDR approach, despite it being a bad framing of the problem (i.e. “effect yes/no” rather than “ranking under cost constraints”). The blog was posted in 2016, and the approach taken in the field has not really changed much since then, so your input would be valuable.

– In the blog post, and also in BDA3, you say that you don’t really like Bayes factors, as having a lump of probability at 0 doesn’t really make any sense. I agree with that, however in this case I don’t really want to make statements about the value of the parameter, just how confident we are that something is going on. Could Bayes factors here be ok?

– I wonder if an alternative approach would be to decide on a minimum interesting effect size, then ask which genes have the highest posterior probability of an effect size larger than this.

– I think elsewhere you have suggested addressing multiple testing-type problems by putting everything into one large hierarchical model. I have over 10k genes, and a few hundred samples, so this doesn’t seem practical (at least not for me!).

– In my analysis, I have calculated Bayes factors (for “full model” vs “only covariates”), but also tried ranking by mean(case_vs_control) / sd(case_vs_control). The Bayes factors seem to have some advantages, e.g. they are especially low for genes that are known to have strong sex effects, and sex is one of the covariates.

My reply: I know very little about genomics–ok, actually I know nothing at all about genomics, although I’ve had the occasional conversation of this sort where it seems that I’ve been able to offer helpful advice. So, rather than attempt any sort of solution, I will just share some scattered thoughts:

1. Why do you want to rank the genes? What are you going to do with the ranking?

My general thought is that, even if you’re doing ranking, there are lots of ways to do this, and it’s hard for me to think of a realistic example in which tail-area probabilities are the right way to do the ranking–see discussion here–except in some very simple symmetric problems where all analyses lead to the same result.

2. My problem with the Bayes factor is not so much with the lump of probability at zero, but rather with the dependence on the prior distribution for the unconstrained parameters in the model. Take K well-identified parameters on unit scale and change their weak prior from independent normal(0,100) to independent normal(0,1000), and you’ve decreased the model’s Bayes factor by a factor of 10^K without changing the inferences conditional on the model. We discuss this sort of thing further in Chapter 7 of BDA3 (chapter 6 of the earlier editions). When comparing models, you can establish conventions with weak priors so as to obtain stable Bayes factors, but that’s what these are–conventions–and there’s reason to expect the resulting Bayes factors to make much sense (that’s the point that I made in my 1995 article with Rubin on the topic), so you’ll just have to evaluate these as one more possible decision rule with no clearly-defined theoretical basis.

3. You’re right that I have suggested addressing multiple testing-type problems by putting everything into one large hierarchical model! Here’s my paper with Jennifer and Masanao explaining this, published in the Journal of Research on Educational Effectiveness in 2012.

4. You write that you over 10k genes, and a few hundred samples, so this doesn’t seem practical. But it is practical! We routinely fit multilevel regressions in Stan with hundreds of thousands of data points and thousands of predictors. Try it and see how it goes! If it’s slow, then my recommendation is to set the group-level variance parameters to reasonable values, then it’s just a simple regression model with a prior, and you can fit it fast using an optimizer, just a simple mode-finder or you can use ADVI or Pathfinder if you’d like.

5. I can see that ranking by mean(case_vs_control) / sd(case_vs_control) might not work so well because the sd can be noisy. You can probably do better by fitting a hierarchical model to the sd’s.

6. Again speaking generally, you should be able to use simulated-data experimentation to get a sense of how any of these methods is working. Simulated-data experimentation has two advantages here: first, by construction you know the true values of all the parameters so you can directly compare the performance of different methods; second, the very act of setting up the simulation forces you to figure out exactly what you are trying to learn, which brings us back to my item #1 above.

Stan for multimodal mixtures—from exponential CPS to linear DP

Posted on by
13

This post is from Bob

I’ve been thinking about evaluation recently because I’ve been working with colleagues on new samplers, which means evaluating how well they work (more on that soon). This in turn means coming up with target densities on which to evaluate them.

A combinatorial multimodal test case

I wanted something clearly multimodal and hence not log concave. I remember somebody’s paper (help with citation?) used a mixture of four two-dimensional isotropic normals, separated enough to make transition possible, but still a bit difficult. Not to give the game away, but here’s a posterior plot of a sample drawn from Stan—the imbalance in component weights is intentional, as I’ll describe below.



I’m still working on posteriordb with the Stan gang (see the authors of the linked paper) and Inference Gym with Reuben Cohn-Gordon (another linguist by training and programming language geek turned to MCMC), and thought it’d be nice to have something a little more general than just the 2D example. So I got out my notebook, and realized the generalization to D dimensions involves 2^D mixture components that are normal with unit covariance located at the points in {-r, r}^D. 

p(y | r) = SUM_{mu in {-r, r}^D} 1/2^D normal(y | mu, I).

I then generalized to allow setting the probability that Y[d] > 0 to be p in (0, 1) to get a non-uniform mixture. This leads to a slightly more complex density because of the non-uniformity.

p(y | r) = SUM_{mu in {-r, r}^D} binomial(sum(mu == r) | D, p) * normal(y | mu, I).

Coding in Stan with continuation-passing style

So how do we code this in Stan? Obviously it needs to be recursive or at least iterative to deal with the D being unknown at compilation time. Whenever I see recursion, I immediately think of continuation passing style (CPS). So I came up with this Stan program to code a generalization in D dimensions.

functions {
  real mm(vector y, real r, real p, int d, real lp) {
    if (d == 0) {
      return lp;
    }
    real lp1 = mm(y, r, p, d - 1, lp + normal_lpdf(y[d] | r, 1));
    real lp2 = mm(y, r, p, d - 1, lp + normal_lpdf(y[d] | -r, 1));
    return log_mix(p, lp1, lp2);
  }

  real mm_lpdf(vector y, real r, real p, int D) {
    return mm(y, r, p, D, 0);
  }
}
data {
  int D;   // number of dimensions
  real r;  // modes in {-r, r}^D
  real p;  // p = Pr[Y[d] > 0]
}
parameters {
  vector[D] y;
}
model {
  y ~ mm(r, p, D);
}

The log_mix function is defined as follows, but implemented in a more stable way.

log_mix(p, lp1, lp2)
    = log_sum_exp(log(p) + lp1, log(1 - p) + lp2)
    = log(exp(log(p) + lp1) + exp(log(1 - p) + lp2))
    = log(p * exp(lp1) + (1 - p) * exp(lp2)).

If you unfold the recursion manually, the leaves wind up being the log densities and the weights wind up percolating as described in the definition. If you’re having trouble seeing this, manually expanding the D = 1 and then D = 2 cases will help. It’s compact, but it’s still exponential in cost to evaluate a log density and gradient (i.e., O(2^D)). Although it’s slow in higher dimensions, it works.

Python scripts

The plot above is from the following Python code that sets

D = 2, r = 2.5, and p = 2.0/3.0.

For those of you considering a move to Python, having a clone of data frames (pandas) and ggplot2 (plotnine) is a godsend. And yes, of course the LLMs know how to code pandas and plotnine.

import cmdstanpy as csp
import pandas as pd
import plotnine as pn

model = csp.CmdStanModel(stan_file='mm.stan')
D = 2
r = 2.5
p = 2.0 / 3.0
data = {'D': D, 'r': r, 'p': p}
fit = model.sample(data = data, iter_sampling=5_000)
print(fit.summary(sig_figs=2))

y = fit.stan_variable('y')
df = pd.DataFrame({'y1': y[:, 0], 'y2': y[:, 1]})
plot = (
    pn.ggplot(df, pn.aes(x='y1', y='y2'))
    + pn.geom_vline(xintercept=[-r, r], color='red', linetype='dashed')
    + pn.geom_hline(yintercept=[-r, r], color='red', linetype='dashed')
    + pn.geom_point(alpha=0.1)
    + pn.scale_x_continuous(breaks=[-r, 0, r])
    + pn.scale_y_continuous(breaks=[-r, 0, r])
    + pn.coord_fixed()
    + pn.theme_minimal()
)
plot.save('mm.jpg', dpi=300)

The knockoff of data frames in pandas and ggplot2 in plotnine are a godsend if you’re transitioning to Python from R (which I would highly recommend).

Dynamic programming to the rescue

Because it involved CPS, I mailed it off to Brian Ward around midnight last night. I’m a decent programmer, but Brian’s next level. By the time I arrived today at 10 am, he had rewritten the target density as follows.

  real mm_lpdf(vector y, real r, real p, int d) {
    if (d == 0) {
      return 0;
    }
    real lower_mixture = mm_lpdf(y | r, p, d - 1);
    real lp1 = lower_mixture + normal_lpdf(y[d] | r, 1);
    real lp2 = lower_mixture + normal_lpdf(y[d] | -r, 1);
    return log_mix(p, lp1, lp2);
  }

[Edit: Switched everything to lpdf from a mix of lpdf and lupdf.]

He saw that the recursions were doing the same thing in each branch and could be shared. Because there’s only one recursive call, Brian’s code is linear (i.e., O(D)). It achieves this speedup using dynamic programming (DP). DP calculates partial solutions that can be combined into larger solutions rather than recomputing them. DP’s the technique that you need to solve the harder L33T-code quizzes you’ll get during technical interviews these days. Other examples where DP can be helpful for statistical models include the fast Fourier transform (FFT), the forward algorithm for hidden Markov models (HMMs), and the Poisson-binomial distribution. The first two are coded efficiently in Stan and the latter I showed how to code in a Stan forum post on Poisson-binomial.

Try it yourself in the Stan Playground

If you want to play with this yourself, Brian built a version using the Stan Playground that you can run in the browser.

Here’s what it looks like after setting D = 3, running sampling, and then viewing a histogram with all three dimensions selected.



It’s a live demo, so you can edit the data to set r, D, and p. And it’s really fast due to the DP. Just like in ShinyStan and especially like its generic in-the-browser version MCMCMonitor (from many of the same developers as Stan Playground), you can view 3D projections of the higher-dimensional draws and rotate them to see it making 8 balls in 3D, 7 of which are visible in the screen grab. Or you can go to higher dimensions and view projections down to two or three dimensions. You might want to increase the number of draws per chain to get cleaner delineation of the posterior densities in the visualizations.

The Dodgers are hiring

Posted on by
Reply

Richard Anderson writes:

I manage the data science team with the Los Angeles Dodgers and have a job post that may be of interest to students or readers of your blog. If you know anyone who may be interested, they are welcome to reach out to me directly.

From the job description:

We are especially interested in candidates with either (1) demonstrated strength in deep learning with applications to spatiotemporal data or (2) demonstrated strength in Bayesian hierarchical modeling & probabilistic forecasting.

And, in the list of areas of preferred expertise:

You have written probabilistic models in NumPyro, PyMC, or Stan with custom likelihoods and priors.

I can only assume that Leo Durocher would be spinning in his grave. “Quantitative Analysis,” indeed!

StanBio

Posted on by
Reply

This is Eric.

On 30 May, we held the first virtual StanBio Connect conference. One large video recording is now available on the conference website and the Stan YouTube channel under the Live tab. We have individual videos, but we have not processed them yet. I was impressed with the depth and the variety of the talks, and I hope we will be able to do it again next year.

In the meantime, I am considering turning stanbio.org into a destination for all things Stan in biomedical research. This would involve restructuring the website (it’s built with Quarto, so that’s not too difficult) and hosting both the original and linked content. If you think it’s a good idea and want to volunteer and contribute, please send an email to [email protected] with the subject line ‘StanBio’.

loo R package 10 years!

Posted on by
5

This post is by Aki.

The loo R package has its 10 year anniversary today! Jonah Gabry made the first loo package release (v0.1.0) 10 years ago on June 26.

  • loo has been downloaded more than 4 million times from the RStudio CRAN mirror (there are more than 80 mirrors, but the RStudio mirror is likely to be one of the most popular ones)
  • R-universe counts 304 other R packages using loo
  • based on R-universe scores, loo is in top 100 among 26,819 packages

Here’s a blog post about the background, advances during the years, and a bit about the future.

Stan

When I (Aki) got involved with Stan project, there was a need for a model selection criterion that would be fast, robust, and easy to compute. I had used cross-validation (CV) a lot, but it did require some expertise to know which computational approach to use in which case, and how to diagnose the reliability. Brute-force leave-one-out cross-validation with repeated model inference is slow. Gelfand et al. (1992) had proposed importance sampling leave-one-out (LOO) CV, but 1) that estimate may have infinite variance (e.g. Peruggia, 1997), 2) due to skewed distribution the estimate would be over-optimistic with high probability, and 3) there was no good diagnostic for reliability.

DIC

Andrew had been using DIC (Spiegelhalter et al., 2002), which was simple and fast, but 1) it assumes the predictions are made using posterior mean of the parameters (and not by integrating over the posterior), 2) it’s not invariant to parameter transformations, and 3) it was known to fail for multimodal posteriors and flexible models.

WAIC

The Widely-applicable information criterion (WAIC; Watanabe, 2010) seemed promising at first, 1) being simple and fast, 2) assuming predictions using posterior predictive distribution, 3) being invariant to parameter transformations, and 4) it works with multimodal and singular posteriors. Alas, more testing of WAIC revealed it also fails with more flexible models without a warning.

Diagnostics

It was now clear that there was a pressing need for a diagnostic for the model comparison criterion computation.

I started to investigate a diagnostic for WAIC. WAIC can be presented as a truncated Taylor series approximation. In difficult cases, the higher order functional cumulant terms are not small. Can we diagnose whether MCMC estimates of higher order terms have finite variance?

Koopman et al. (2009) had proposed a diagnostic for importance sampling: 1) fit generalized Pareto distribution (GPD) to the tail of the importance ratio distribution, 2) use shape parameter k to diagnose whether the variance is finite (k <= 1/2), and 3) reject the estimate if the variance is not finite. The problem was that there was no advice what to do if the estimate is rejected and this seemed to happen a lot with WAIC and importance sampling LOO. Furthermore, finite but high variance is also problematic.

I then realized that if we assume the tail of the ratio distribution is close to a generalized Pareto distribution, and we fit the generalized Pareto distribution, we can use that as a model for the tail. To make the computation practical, this model is used to replace the raw ratios with modelled (smoothed) ratios, which are then used in further computations. In theory, modeling should reduce variability (and if the model is good, the bias can be negligible) and this was observed in practice; smoothed importance sampling LOO performed better than WAIC (or had similar performance for simple boring models).

Jonah implemented the method in loo package while we were making the experiments and writing the paper, and loo v0.1.0 was released on Github on June 26 before the papers were public.

Main papers

As theory and experiments did take a lot of pages, we decided to split the paper into “Very Good Importance Sampling” (arXived July 9, the name inspired by WAIC) and “Efficient implementation of leave-one-out cross-validation and WAIC for evaluating fitted Bayesian models” (arXived July 16). Eventually the paper names were changed to “Pareto smoothed importance sampling” (PSIS) and “Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC”. The latter did get published in Statistics and Computing in less than 14 months, while the PSIS paper (final version 58 pages) took more than 9 years to get published eventually in JMLR.

While the PSIS algorithm stayed practically the same the whole time, the theoretical justifications and diagnostics did improve over the years. Dan Simpson and Yuling Yao provided help with the theory and joined as co-authors.

The PSIS estimator always has finite variance with a cost of some bias. If Pareto-k<0.7, the bias and variance are small (see details in the PSIS paper; Vehtari et al., 2024). Although 9 years felt too long time and we had a nasty case of reviewer 2, I’m happy about the improved theoretical understanding of the pre-asymptotic behaviour.

Outliers and influential observations

The loo package did gain additional practical advice on how examining the number of parameters, the effective number of parameters (p_loo), and the number of observations can provide information on whether the high Pareto-k values are likely due to a) badly misspecified model with outliers or b) well specified but flexible model.

As we get Pareto-k diagnostic for each LOO-fold, that can be used as identifying influential or problematic observations, but also to focus the additional computation only for the specific LOO-folds.

Additional computation for high Pareto-k cases

The simplest approach is to just re-run MCMC for the folds with high Pareto-k. rstanarm and brms packages know enough about the data and model, so that they can provide automated approach for this.

To speed-up computation, we developed moment matching approach (Paananen et al., 2021) that can adjust the posterior draws faster than what re-running MCMC would take, to better match the proposal and target. The needed functionality and a vignette was added to loo. brms makes it easy to use moment matching LOO with one option.

The loo package also has a vignette and support for K-fold-CV which is robust and relatively fast if, e.g. K=10.

Large data

Even though PSIS-LOO is generally fast, with big enough data it can be slow. We developed a sub-sampling LOO approach (Magnusson et al., 2019, 2022) with a vignette in the loo package. In the sub-sampling approach we use a faster but biased estimate for all LOO-folds and a slower but (almost) unbiased estimate for a subset of LOO-folds. The biggest benefits of this approach can be seen in the projpred package (github version), where PSIS-LOO given the full data search path is fast, but doing the search for each LOO-fold is slow. Using subsampling and the difference-estimator we can get more than a 10-fold speedup as demonstrated in one of the case studies.

Predictive checking

Cross-validation can be used to improve predictive checking. Posterior predictive checking can fail with flexible models as the same data are used for fitting and checking. We added useful functions to loo package to support LOO predictive checks and LOO probability integral transformation (LOO-PIT) calibration checks to bayesplot (Gabry et al., 2019; Säilynoja et al., 2022, 2025).

Beyond LOO-CV

The package is named loo as it started as an implementation of the PSIS-LOO algorithm (and we had only US and Finnish people thinking about the name). But it was natural to extend it beyond LOO.

Leave-one-group-out (LOGO) is useful if we want to predict for new groups. LOGO is challenging for importance sampling as the posterior for group specific parameters changes a lot if we leave out all the group specific observations. We can use K-fold-CV (loo vignette) or integrate out the group specific parameters (Roaches case study).

While LOO is valid for analysing the observation model in time series models, we may sometimes prefer leave-future-out cross-validation (LFO-CV), as it has smaller bias if the prediction task is in the future. The loo package has a vignette demonstrating how PSIS and occasional re-fits can be used for fast LFO-CV (Bürkner, Gabry and Vehtari, 2020).

The downside of LFO-CV is that it uses only a small part of the data for fitting the model and making predictions for the future, and thus has high variance. If the focus is model comparison, it is better to use K-fold-CV or hv-block-CV with joint log score as these have smaller variance, which leads to higher model selection efficiency for time series models (Cooper et al., 2025a).

Leaving out more than one observation and using joint log score improves the model selection efficiency also in the case of spatial models (Cooper et al., 2025b).

Often temporal and spatial models are presented as non-factorized normal (or t) models, in which case we need to use properties of multivariate normal (or t) to compute LOO (Bürkner, Gabry and Vehtari, 2020). The loo package has a vignette and brms uses this approach for non-factorized models.

Comparing models

Instead of just using point estimates of the predictive performance, we can also quantify the related uncertainty, which is especially useful when doing model comparison. From the beginning, the loo package was reporting the log score (elpd) difference and the related standard error based on the recommendation by Vehtari and Lampinen (2002). We (Sivula et al., 2025) have investigated in more detail the conditions when the standard error and related normal approximation are accurate. The normal approximation can be used to estimate probability that model B has better predictive performance than model A.

Model averaging

The LOO computed log score (elpd_loo) has a connection to information criteria (see, e.g. Vehtari and Ojanen, 2012). Inspired by information criteria based model weights and stacking, we developed Bayesian stacking (Yao et al., 2018), which is also implemented in the loo package. Furthermore, we extended Bayesian stacking to Bayesian hierarchical stacking (Yao et al., 2022a) and stacking for non-mixing computations (Yao et al., 2022b). There is a loo package vignette for stacking.

Other scores and metrics

The loo package supports also (S)CRPS, MAE, RMSE, MSE, ACC, and BACC, although not as nicely as log score (see below for future plans).

ArviZ

loo is an R package, but PyStan and PyMC users needed fast cross-validation, too. The Python and Julia ArviZ libraries were subsequently developed to include most of the methods that are in the loo R package,

Future

Because PSIS is useful in many other cases beyond just LOO, PSIS functionality has now been implemented in the posterior package, so that packages that would like to use just the Pareto smoothing and diagnostics do not need to depend on the loo package. We’re also in the progress of changing loo to use more of these functions from posterior.

We’re also refactoring loo to improve modularity and usability (as part of Google Summer of Code) focusing on:

  • easier use of different scores and metrics (e.g. RMSE, R^2, CRPS)
  • easier use of different cross-validation variants
  • easier use of joint log score
  • more diagnostic information provided to the user, e.g. for the uncertainty normal approximation

The ArviZ team is also working to update ArviZ.

CV-FAQ

The loo package users have been asking many questions, and eventually I wrote CV-FAQ answering the most frequently asked questions (and setting straight some common misconceptions).

Thanks

Very big thanks to Jonah Gabry for writing the loo package (see also other contributors) and to the ArviZ team for implementing the methods in Python and Julia! All the methods developed in the papers would not be widely used without these packages (The Practical LOO-CV paper has been cited more than 5900 times). While the methods in the papers made certain things possible, it is the software that makes using the methods easy!

Stanford Human Trafficking Data Lab Hiring for a Full-time Postdoctoral Scholar

Posted on by
1

Ben Seiler sends along this opportunity:

The Stanford Human Trafficking Data Lab is accepting applications for a postdoctoral fellowship position to join a project investigating trafficking risks in charcoal supply chains in Brazil. The position is open to recent graduates of PhD programs in statistics, economics, computer science, operations research, or related data science fields. The position provides opportunities to participate in rigorous, quantitative research on human trafficking, including supply chain network analysis and geospatial modeling. The successful candidate will have strong data science skills, including experience working with large, complex data from varied sources, and machine learning methodologies. The underlying data are complex and will require sophisticated data management and integration skills. A candidate should have proficiency with GIS software and Python, strong written and interpersonal communication skills, and a demonstrated interest in addressing social justice issues through data-driven research. The postdoc will work in partnership with PI Grant Miller (Stanford University School of Medicine) and other research team members, and will contribute to study design, participate in field research, conduct data analysis, and disseminate findings through academic publications and presentations. The postdoctoral fellow will be expected to focus mainly on this project, but may spend up to 20% of their time on independent research. For more on the Lab’s other ongoing projects see https://htdatalab.stanford.edu/projects/. The postdoc will be based at Stanford University in the Department of Health Policy, located near to the Departments of Economics and Political Science and the Graduate School of Business.

TO APPLY: Please email a single PDF document named
“Lastname_Firstname_HTDL_Postdoc_2025” containing the following materials to Lydia Aletrais at [email protected]:

1. A cover letter describing your interest in this position, your relevant training and experience, and our earliest and preferred start dates
2. A current CV
3. A transcript (unofficial is fine)
4. Names, e-mail addresses, and phone numbers of 2-3 references.

Applications will be considered on a rolling basis. Short-listed applicants will be asked to complete a technical exercise and may be called for an interview.

Given that they’re doing network analysis and geospatial modeling, I expect that knowledge of Stan would be useful too.

Election analytics positions available at the New York Times

Posted on by
Reply

Will Davis writes:

I oversee the Election Analytics department (aka The Needle and The Times/Siena Poll) at The Times.

We’re lucky enough to have two great roles open on the team right now for people excited to make a career in the field of election analytics.

* Election analyst: This is a mostly technical role, but it comes with the opportunity to write. This person would be responsible for some most essential elements of our election polling and modeling — modeling unit-level turnout and vote share, creating baselines for The Needle and helping us continue to innovate the design of The Times/Siena Poll.

* Election researcher: This is a great job for somebody with less of a technical skillset who’s eager to be around a team they can learn a ton from. It’s primarily focused on manual research and data entry.

These sound like excellent an opportunity to combine statistical modeling, computation, and graphics.