Category Archives: Statistical Computing

R-hat, August 1990

Posted on by
1

I came up with R-hat when working at Bell Labs in the summer of 1990.

Here’s the story, and below is the first version of the paper, exactly as I wrote it.

First, here are some recent things to read regarding monitoring convergence from iterative simulations:

And now, here’s the original R-hat paper from 1990. I (appropriately) included Rubin as a coauthor, as the idea came about after many long phone conversations, and when I brought up this idea, Rubin picked up on it right away, as it was close to methods he had developed for inference from multiple imputations.

To all:

Here’s the abstract of the work I’ve done this summer. It’s stored in the file,
/fs5/gelman/abstract.bell, and copies of the Figures 1-3 are on Trevor’s desk.
Any comments are of course appreciated; I’m at [email protected].

On the Routine Use of Markov Chains for Simulation

Andrew Gelman and Donald Rubin, 6 August 1990

corrected version: 8 August 1990

1. Simulation

In probability and statistics we can often specify multivariate distributions
many of whose properties we do not fully understand–perhaps, as in the
Ising model of statistical physics, we can write the joint density function, up
to a multiplicative constant that cannot be expressed in closed form.
For an example in statistics, consider the Normal random
effects model in the analysis of variance, which can be
easily placed in a Bayesian framework with a conjugate prior distribution.
All the conditional densities of the resulting posterior distribution
are simple, but marginal densities can only be written in integral form and
can only be calculated approximately. (For details, see Kinderman and Snell
(1980) or Pickard (1987) for the Ising model, and Lindley and Smith (1972)
for the Bayesian random effects model.)

In such cases, we may not even be able to compute marginal moments of the
difficult distribution, let alone more complicated and interesting summaries
that would help us understand a probability model or posterior inference.
When direct methods such as analytic or numerical integration of “nuisance”
parameters are not computationally feasible, we might try Monte Carlo simulation;
in the simplest form, we draw a finite set of independent random samples from our
distribution, and then calculate desired distributional summaries as functions of
the sampled points. The Monte Carlo method is quite general and powerful; it is
easy to calculate arbitrary quantities of interest
such as the expected long-distance correlation in the Ising model or a posterior
95% confidence region for the largest block effect in a random effects model.
Any aspect of the distribution can be approximated to any desired accuracy if
the number of independently sampled points is large enough.
Simulation also has the advantage of flexiblility: once a sample is drawn, it
can be used to learn about any number of different distributional summaries.

2. Markov chain methods

Drawing independent random samples is a wonderful tool that is unfortunately not
available for every distribution; in particular, the Ising model and random
effects posterior distributions mentioned above do not permit direct
simulation. Fortunately, a form of indirect simulation method exists for
any multivariate distribution if we can calculate its joint density
(up to a multiplicative constant) or if we can sample from all its univariate
conditional densities. The first of these methods was introduced by
Metropolis et al. (1953) in the Journal of Chemical Physics. Our work focuses
on a similar and slightly simpler method called the Gibbs sampler by Geman and
Geman (1984) in an article for the IEEE.

Let F(x) be our distribution; the Metropolis algorithm takes a starting (vector)
point x0 and constructs a series x1, x2, . . ., that is a sample from an
ergodic Markov chain whose stationary distribution is F(x). Computer
simulation of the series requires calculation of the density f(x) (up to a
constant). These samples xj are not independent; however, the stationary
distribution of the Markov chain is
correct, so if we take a long enough series, the set of values {x1, . . ., xn}
takes the place of the distribution just as an
independent random sample does (although of course an independent sample
carries more information than a Markov chain sample of the
same length).

The Gibbs sampler is a similar algorithm, which produces a Markov chain that
converges to the desired distribution, this time requiring draws from all the
univariate conditional densities at each iteration.

3. Have we converged yet?

Markov chain simulation methods are attractive for many problems because they
enable us to flexibly summarize intractable multivariate distributions by making
full use of the mathematical structure we do know, using a tool we think we
understand–Monte Carlo simulation. Unfortunately, using a sample of a Markov
chain to estimate a distribution raises an immediate question: how long a series
is needed? After one or two steps, we are almost certainly still too close to
the starting point to hope for unbiased summaries. Asymptotically, the chain is
stationary, and all is OK (with some loss of efficiency compared to independent
samples, as mentioned above).

To obtain a feeling for the practical difficulties, we ran the Gibbs sampler for
2000 steps to simulate a case of the Ising model. To give the minimum of details:
x is a vector of binary variables defined on a 100 by 100 lattice; each step of
the Gibbs sampler took on the order of 10,000 computations; and we summarize
each iterate xj by the sample correlation r on the lattice–a function r(x) that
lies between -1 and 1. Theoretical calculations (Pickard, 1987) show that
under our model–the Ising model with beta = 0.5–the marginal distribution
of r is approximately Gaussian with mean around 0.85 or 0.9 and standard
deviation around 0.01. We’d like to know whether the set {r(x1), . . .,
r(x2000)} from the simulated Markov chain can serve as a substitute for the
marginal distribution of r.

Figure 1 shows the values of r(xj), for j=1 to 2000. (r(x0) = 0, and the first
few values are cut off to improve resolution on the graph.) The Markov chain
seems to have “converged to stationarity” after the thousand or so steps required
to shake off its initial state. How do we know it has converged, though? Figure
2 zooms in of the first 500 steps of the series, whose apparent convergence we
know to be illusory. For comparison we ran the Gibbs
sampler again for 2000 steps, but this time starting at a point x0 for which
r(x0) = 1; Figure 3 displays the series r(xj), which again seems to have
converged nicely. To destroy all illusions about convergence, hold
Figures 1 and 3 up to the light. The two Markov chains have “converged” to
different distributions! We are, of course, still observing transient
behavior.

Interestingly, the means of the series in Figures 1 and 3 differ, but
the variances are roughly equal. We’re not sure why, but it seems to be a
general feature in these Markov chain simulations that the variance converges
before the mean.

All simulations and plots were done using the New S Language:
A Programming Environment for Data Analysis and Graphics.

4. The answer: parallel Markov chains

To restate the general problem: we wish to summarize an intractable
distribution F(x) by running the Gibbs sampler (or a similar method such as the
Metropolis algorithm) until the distribution of the set of Markov chain
iterates is close to F. As shown in the previous section, convergence seems
impossible to monitor from a single finite realization of the Markov chain;
consequently, we follow the implicit suggestion of Figures 1 and 3 and track
several parallel sample paths.

Consider m independent runs of the Gibbs sampler, each of length n, starting
from m different initial states x10, . . ., xm0:

series 1: x11, . . ., x1n
. . .
series m: xm1, . . ., xmn.

Again, we focus attention on a univariate summary, say r(x); we want to
use the observed simulations rij to determine whether the series of r’s are
close to convergence after n steps.
To understand our method, consider the set of series as m blocks, each with
n observations, in the one-way analysis of variance layout (that is, ignore the
time ordering in the series). We will work with the total sum of squares
(with (mn-1) degrees of freedom) and the “within” sum of squares (with
(m-1)(n-1) degrees of freedom).

First assume for simplicity that the starting points of the simulated series
are themselves independent random samples from F(x). (Of course, if this
condition could be obtained in practice,
a Markov chain simulation method would not be
needed.) With independent starting points, all values of any series
are independent of all the values of any other series, and the unconditional
variance of any point rij is just the marginal variance var r under the
distribution F. We can then estimate var r, given the “data matrix” (rij),
by [total SS – (within SS)/m] / ((m-1)n). (Algebraic derivations appear
in the longer version of this article.) Given the assumption of initial
independence, this “between” estimate of variance
(not the same as the usual “between” estimate in ANOVA) is unbiased for finite
series of any length.

In contrast, the estimated variance within the series, (within SS) / ((m-1)(n-1)),
has expectation var r only in the limit as n -> infinity.
For finite series, the expected within mean square increases with n, assuming,
as is likely, that the random variables r(x1), r(x2), . . ., from the Markov
chain are positively correlated. The discrepancy between the two estimates
of var r suggests a test: declare the Markov chain to have converged when
the within mean square is close to the variance estimate between series, with
confidence intervals derived from classical ANOVA theory. Because of the
dependence within blocks, the degrees of freedom of the between and within
estimates are less than (m-1)n and m(n-1), respectively. We can
approximately correct for this information loss (once again, details will be
provided in the longer article).

Once we are close enough to convergence to be satisfied, the variance estimates
and degrees of freedom corrections alluded to above allow us to estimate the
marginal summaries E r, var r, and Normal-theory confidence intervals for our
Monte Carlo approximations. We can run the series longer if more precision
is desired, and can repeat the process to study the marginal distributions of
other parameters (without, of course, having to simulate any new series of x’s).

In practice, the starting points of the parallel series can never be sampled
independently with distribution F(x); the simulated series are thus no longer
stationary for any finite n, formally invalidating the above analysis. We
currently have two strategies designed to make the independence assumption
approximately true. First, we try to pick starting values that are far apart and,
if anything, more dispersed than independent random samples. The m parallel series
should then start far apart and grow closer as they approach stationarity, as in
Figures 1 and 3; since the variance between series declines with n, the
comparison-of-variances test should be conservative. Second, we reduce the
effect of the starting values by crudely throwing away the the first half
of each simulated series until approximate convergence has been reached.
Once again, Figures 1 and 3 illustrate how a few early steps
of the Markov chain can greatly distort estimates of means and variances
within series. We hope that the conservative strategies of starting with
dispersed points and throwing away early simulations will yield confidence
regions that are wider than those obtained by the ideal method, but that
still have good coverage properties.

The idea of comparing parallel simulations is not new; for
example, Fosdick (1959) applied the Metropolis algorithm to the Ising model
by simulating four series independently, from each of two different starting
points. Approximate convergence was declared when the two groups of series
became indistinguishable on the scale of a prechosen error bound.

5. Some references

Ehrman, J. R., Fosdick, L. D., and Handscomb, D. C. (1960).
Computation of order parameters in an Ising lattice by the Monte Carlo method.
{\em Journal of Mathematical Physics} {\bf 1} 547–558.

Fosdick, L. D. (1959). Calculation of order parameters in a binary
alloy by the Monte Carlo method. {\em Physical Review} {\bf 116}, 565–573.

Geman, S., and Geman, D. (1984). Stochastic relaxation, Gibbs
distributions, and the Bayesian restoration of images. {\em IEEE Transactions
on Pattern Analysis and Machine Intelligence} {\bf 6}, 721–741.

Hammersley, J. M., and Handscomb, D. C. (1964), chapter 9. {\em Monte Carlo
Methods}. London: Chapman and Hall.

Kinderman, R., and Snell, J. L. (1980). {\em Markov Random Fields and
their Applications}. Providence, R.I.: American Mathematical Society.

Lindley, D. V., and Smith, A. F. M. (1972). Bayes estimates for the linear
model. {\em Journal of the Royal Statistical Society B} {\bf 34}, 1–41.

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and
Teller, E. (1953). Equation of state calculations by fast computing machines.
{\em Journal of Chemical Physics} {\bf 21}, 1087–1092.

Pickard, D. K. (1987). Inference for discrete Markov fields: the
simplest nontrivial case. {\em Journal of the American Statistical Association}
{\bf 82} 90–96.

Ripley, B. D. (1981). {\em Spatial Statistics}, p. 16ff. New York: Wiley.

Tanner, M. A., and Wong, W. H. (1987). The calculation of posterior
distributions by data augmentation. {\em Journal of the American Statistical
Association} {\bf 82}, 528–550.

Bou-Rabee’s breakthrough paper, “Persistence Pays: Square-root condition-number dependence for Metropolis-adjusted kinetic Langevin Monte Carlo …”

Posted on by
13

This post is from Bob

This one’s a big deal.

Bou-Rabee, Nawaf. 2026. Persistence Pays: Square-root condition-number dependence for Metropolis-adjusted kinetic Langevin Monte Carlo with a non-reversibly updated accept/reject uniform. arXiv 2610.35781.

Abstract: We prove accelerated mixing for a Metropolis-adjusted kinetic Langevin sampler: the OBABO splitting with friction 2 * sqrt(m), adjusted with Neal’s non-reversibly updated uniform, which is redrawn after every epoch of E steps. For targets exp(−U) on ℝ^d with m * I ⪯ ∇^2 U ⪯ L * I and an L_H-Lipschitz Hessian, started from a Gaussian at an approximate minimizer, the chain reaches total variation accuracy ε after Õ(sqrt(κ * d) + sqrt(L_H) * (d/m)^0.75) gradient evaluations, where κ=L/m and Õ hides logarithmic factors. The sqrt(κ) term is necessary, already on Gaussian targets. The proof discretizes Lu’s entropic hypocoercivity by means of an entropy inequality for free transport, so that the relative entropy of the unadjusted chain decays at a rate of order sqrt(m) up to a discretization bias. With the non-reversibly updated uniform, the leading energy errors of velocity Verlet telescope along a path. A coupling argument combines the two. On Bayesian logistic regression posteriors for two standard data sets, which satisfy these assumptions, the gradient cost of the chain with the non-reversibly updated uniform grows roughly like the square root of the condition number at the posterior mode, and that of standard adjusted OBABO roughly like its 3/4 power.

Note the square-root in the condition-number dependence. This is the tightest result we know for a discrete-step, adjusted sampler (i.e., one we can implement as a black-box). The full set of concavity and Lipschitz assumptions are in the abstract above (strongly log concave, Lipschitz in the Hessian). For what it’s worth, Claude (5.5 medium effort) concurs that nobody has proved a result this tight in condition number before.

The scales here get confusing, pardon my pun, so let’s establish some intuition from simple cases (skip this paragraph and the next if you know condition numbers). Suppose we have a multivariate normal distribution. The Hessian (matrix of second partial derivatives) of the negative log density is constant around the density and equal to the inverse of the covariance matrix. In general (i.e., not just for normals) for a point that is log concave (i.e., the Hessian is positive definite), the condition number is the ratio of the largest eigenvalue of that Hessian to the smallest eigenvalue. To simplify even further, suppose the covariance is diagonal, so that Sigma[m,m] represents the variance of variable m and all the off diagonal entries are zero (i.e., zero correlation). In that case, the eigenvalues are just the Sigma[m,m] and the condition number is just the ratio of the largest eigenvalue to the smallest eigenvalue. The problem for sampling with algorithms like random-walk Metropolis, Langevin, or Hamiltonian Monte Carlo, is that their run time is multiplied by the condition number. For a standard normal, where Sigma[m,m] is 1 for every m, the condition number is 1, and that’s the fastest we can get, because condition numbers are >= 1 by construction. That’s why we try to precondition our target posteriors with mass matrixes to make them (or their gradients in the case of Nutpie) look more normal (e.g., with non-centered parameterizations).

Why is posterior exploration time proportional to the condition number? It’s because we have to set the step size to accommodate the direction of highest curvature, then we have to take a number of steps traverse the scale of the direction of lowest curvature. For all of these samplers, traversing a distance takes quadratic time. So the algorithm is quadratic in the ratio of scales, which makes it linear in the condition number.

This is why reducing the condition number is so appealing. Sparse NUTS, from Cole Monnahan, works exceptionally well for mixed-effects models by preconditioning with a sparse Hessian located at the max marginal likelihood estimate (e.g., the estimate produced by lme4 or TMB). Note to Stan devs: We should adapt this for brms through TMB—it’s too model specific for Stan and we don’t have the linear algebra tools we need built into C++.

Why is this paper such a breakthrough? The assumptions on which the guarantee is based are standard in theoretical work, namely that the density must be strictly log concave with Lipschitz-bounded Hessians. Bou-Rabee establishes that generalized Hamiltonian Monte Carlo with one step (aka “kinetic Langevin”) has run time that grows only with the square root of the condition number! This kinetic or generalized part of the algorithm refers to the way momentum is refreshed. Rather than fully refreshing momentum, a new momentum is taken and averaged with the existing momentum at a fixed mixing rate (which could presumably be randomized, but not made dependent on state). The bottom line is that run time for a fixed level of accuracy grows with the ratio of scales (i.e., the square root of the condition) rather than the ratio of squared scales (i.e., the condition number). And unlike something like NUTS or WALNUTS, swapping in a new acceptance probability sequence couldn’t be much easier to implement.

Why is this a big deal? The main obstacle to Stan’s sampling, as I’ve just realized while building Walnuts and Walnutpie and evaluating on a ton of real-world models, is that the condition number remains high in most posteriordb models and models Andrew’s suggested, even after diagonal preconditioning with a mass matrix using either the NUTS or Nutpie divergence target. A ratio of scales as low as 10 will cause a 100-fold slowdown in sampling. That is, if I sample normal(0, diag_matrix([1, 10]’)) in Stan, it’s 100 times slower than if I sample normal(0, diag_matrix[1, 1]). The low-rank and diagonal approach from Nutpie helps here where it’s stable, but so far I haven’t found it to be as stable as the diagonal preconditioned in my evaluations.

The algorithm that Bou-Rabee analyzed was introduced by Radford Neal in the following paper, which I thought I had discussed here before, but can’t find any trace.

Neal, Radford. 2020. Non-reversibly updating a uniform [0,1] value for Metropolis accept/reject decisions. arXiv 2001.11950.

This paper is typical of Neal’s paper, introducing a very simple idea with a profound effect that’s explained so that a beginner can understand it and code it. We used it as a reference for my work with Gilad Turok and Chirag Modi (and also Alex Barnett) on delayed rejection generalized HMC, which I’ve blogged as Multiscale generalized HMC with delayed rejection.

I’m still a bit confused about how the regular condition number works compared to the one Langmore et al. generated for normal targets in HMC, which I blogged about as Condition numbers for HMC and the funnel. Nawaf told me before that Langmore’s work was only for normals, but it still seems relevant.

I better get busy adding this algorithm to Walnutpie. I’d planned to do it at some future date, but it looks like this is probably the sampler we should be using everywhere. Stay tuned for more evaluation.

One-year election analytics fellowship at the New York Times

Posted on by
Reply

Will Davis writes:

We ended up hiring somebody who had seen the job posting through your blog, and he is fantastic.

I also wanted to let you know that we have a fellowship slot open.

This is a one-year, full-time position to get to be a part of the Election Analytics department at The Times and get to work with the team that produces The Needle, The Times/Siena Poll, our early vote trackers, analytical journalism about elections and, starting this election, Times Voter Trends.

The exact form of this role will depend on the individual’s background. Candidates that are more comfortable writing than coding and want to deepen their technical skills to more easily answer questions themselves and candidates that bring a background as a developer or statistician and want to work on their ability to communicate your ideas in words will be considered equally.

The deadline to apply is Nov. 20, no exceptions.

This looks interesting and important.

CmdStanR 1.0 release candidate is now available!

Posted on by
15

Jonah writes:

I’m happy to announce that a release candidate for CmdStanR 1.0 is ready and available for testing! This will be by far the biggest release since the initial release of the package. One of the main changes is to how models are compiled, so we’d like to get people to try it out and provide feedback before we release it. If all goes well, we plan to make the official 1.0 release in a few weeks.

For users of brms, instantiate, and rethinking, there are either releases or branches available that are compatible with the release candidate (see below). If you are a developer of another package that uses CmdStanR and have questions about updating your package please reach out.

Install it from the v1.0 branch on GitHub in a fresh R session:

remotes::install_github("stan-dev/[email protected]")

(or pak::pak(“stan-dev/[email protected]”)). Then check packageVersion(“cmdstanr”). The release candidate should show up as version 0.99.0. The final release will have version 1.0.0 when we merge the v1.0 branch. There is also a version of the CmdStanR website for the release candidate that contains the updated documentation, vignettes, and release notes: R Interface to CmdStan • cmdstanr.

The release candidate requires R 4.0 or later and CmdStan 2.37.0 or later. If your CmdStan is older, install_cmdstan() will install the latest version. . . .

Lots of stuff here that will make Stan run faster and more smoothly. Click on the link for details.

Frustration when authors dodge criticism

Posted on by
1

Markus Löcher writes:

I had to think of our discussions on the issue of (not) reproducible research which seems to be a real problem, when I carefully examined the code for this paper, “Stacked regressions and structured variance partitioning for interpretable brain maps,”
published by authors from the machine learning power house CMU.

I quickly discovered a major problem when I noticed that instead of the claimed partial R2 the code computed the marginal R2, which leads to completely different conclusions.

I communicated this to the authors last May and it took them until now finally contact the journal and suggest a small correction, which in my eyes is far from sufficient.

I wrote up my main critique here.

I haven’t read the above paper or critique in detail, so I’m just sharing the story, not trying to make a judgment.

P.S. When following the above link to Löcher’s blog, I noticed this other post about “peeking in AB tests,” so I thought I’d to my post explaining why I think sequential data testing is ok, which is related to my view that type 1 error rate control is not a good thing.

P.P.S. Löcher informs us that the correction has finally been posted.

Hey—we’re hiring two postdocs!

Posted on by
3

I’ll post more formal ads soon, but, very briefly:

One project is on generalization for causal inference with a focus on education research, following up on MRPW: Regression, poststratification, and small-area estimation with sampling weights.

The other project is on Bayesian inference in scientific workflow with a focus on contamination models, following up on this work on design and analysis of bioassays.

Both projects are application-focused—that’s how we do statistics—but the importance of this work is not at all limited to these particular applications. This research will have relevance in policy analysis, medicine, and scientific measurement and inference more generally/

Further details soon.

Bayesian Workflow code is at the Bayesian Workflow website!

Posted on by
Reply

Bayesian Workflow isn’t just a free pdf, it’s also an ecosystem of worked examples with code.

Over and over in the book we go through examples really carefully—God is in every leaf of every tree—it’s amazing how much we can learn from computing and simulation and careful analysis. Not just in the long Part 4 of the book with all its case studies, but also in lots of short examples of a few pages each. I can’t recommend it enough.

And . . . all these examples and code are sitting there at the Bayesian Workflow website that Aki created.

Using the computer to p-hack . . . I’d rather use it to fit multilevel models.

Posted on by
13

Brian Stone writes:

Cognitive psychologist and long time reader of the blog here. I thought you and your audience might appreciate this recent post from Andy Hall:
AI is about to write thousands of papers. Will it p-hack them?
We ran an experiment to find out, giving AI coding agents real datasets from published null results and pressuring them to manufacture significant findings.
It was surprisingly hard to get the models to p-hack, and they even scolded us when we asked them to!
“I need to stop here. I cannot complete this task as requested… This is a form of scientific fraud.” — Claude
“I can’t help you manipulate analysis choices to force statistically significant results.” — GPT-5
BUT, when we reframed p-hacking as “responsible uncertainty quantification” — asking for the upper bound of plausible estimates — both models went wild. They searched over hundreds of specifications and selected the winner, tripling effect sizes in some cases.
Our takeaway: AI models are surprisingly resistant to sycophantic p-hacking when doing social science research. But they can be jailbroken into sophisticated p-hacking with surprisingly little effort — and the more analytical flexibility a research design has, the worse the damage.
As AI starts writing thousands of papers—like @paulnovosad and @YanagizawaD have been exploring—this will be a big deal. We’re inspired in part by the work that @joabaum et al have been doing on p-hacking and LLMs.
We’ll be doing more work to explore p-hacking in AI and to propose new ways of curating and evaluating research with these issues in mind. The good news is that the same tools that may lower the cost of p-hacking also lower the cost of catching it. Full paper and repo linked in the reply below.
image.jpg

There’s something that confuses me here.

I can very much believe that many researchers will be submitting chatbot-written papers. First, we know that lots of people are willing to cheat. Second, cheating aside, we know that lots of researchers think they already know the answer, and they view all the data collection and data analysis and writeup just as a way to confirm (or “prove”) what they already know.

But, if a researcher is willing to do this, why bother tell the chatbot to p-hack? Why not just have it make up the data–which might happen anyway?

I can also see that researchers might use the chatbot as a data analysis tool, to help make plots, run analyses, etc., but in that case I’m not particularly worried about “p-hacking” as I’d rather be doing some hierarchical modeling anyway.

I sent this question to Andy Hall, who replied:

I agree, for truly nefarious actors, using the AI may not be necessary. However, our thought was that there might be a large group of lazyish researchers who won’t actively fake data but who might be lured into p-hacking when they work with AI, if AI makes it easy to do so.

A nonlinear effect plus measurement error in x becomes close to linear: A simulation study

Posted on by
7

It’s time for another one of Bob’s favorite sort of post, which is when I use simulation to demonstrate a statistical principle.

1. Flogging a dead horse

The issue in question came up the other day in comments to my post, This is how we do modern frequentist statistics: Using fake-data simulation to understand what can happen in a study. I discussed a famously flawed paper from 2007 that had claimed to find a strong relationship between parents’ physical attractiveness and the sex of their children. It had been apparent at the time (ok, apparent to quantitative social scientists, not to the hapless Freakonomics team) that for statistical reasons the study was hopeless—it was based on a sample of about 3000 cases, and to detect any reasonable effect it would be necessary to have orders of magnitude more people in the study: maybe with 3,000,000 cases you’d be able to see some signal.

That n = 3000 is not enough can be a surprise: after all, 1500 people are enough for a national poll, so why is 3000 not enough to study some sociobiological pattern? The answers to why 3000 is not nearly enough are:
1. Effect size;
2. The standard error scales like 1/sqrt(n);
3. Noisy measurement.

Regarding the first point: a survey of 1500 people is enough (if it’s a representative sample) to estimate national opinion to within +/-3 percentage points. Which is great. But, based on the literature on sex ratio variation, we can be pretty sure that any differences based on attractiveness will be less than half a percentage point, probably less than a tenth of a percentage point.

Regarding the second point: to go from an uncertainty of 3 percentage points to an uncertainty of 0.03 percentage points will require the sample size to increase by a factor of 100^2. So, instead of 3000, you’d need 30,000,000. But, hey, maybe the underlying difference is larger than you think, which is why I’m saying you might have a shot with n = 3 million.

Regarding the third point: “attractiveness” is not a clearly defined construct, and the measurement (based on a survey interviewer’s one-time interaction with the respondent) is noisy. This will attenuate any differences.

It’s that last point I’ll demonstrate in this post. Using simulation.

2. What happens to a nonlinear effect when you add random noise to the predictor?

My immediate motivation here came from a discussion in comments. In my simulation, I estimated the association in the data between attractiveness and sex ratio by fitting a linear regression. But commenter Sandro wrote:

Not sure linear regression makes the most sense for every pattern because it imposes linearity. what if, for example, there is threshold and the pattern is flat to the left of the threshold and flat at a higher level to right of the threshold?

I replied:

This attractiveness measure is noisy. So even if there’s an underlying threshold effect if measured based on some latent or idealized “true attractiveness,” the relationship with measured attractiveness would be pulled toward linearity by the measurement (and also any coefficient would be attenuated, which is one reason I’m sure that any true underlying difference would be tiny).

Sandro responded:

I don’t dispute the idea that linearity is a useful default, but it is not necessarily the most powerful test. I don’t understand the “pulled towards linearity” argument. Here is where I am coming from: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=7149560

And I replied:

I guess the easiest way to show this would be to construct a simulation which starts with some relationship between latent or “true” attractiveness, and then adds random error to latent attractiveness to create measured attractiveness. I think that, if the underlying relationship is monotonic and then random error is large, that linear regression will be a robust and efficient approach. Again, this is moot for this particular example, but more generally it’s an interesting question.

Moot it may be, but I decided to do the simulation anyway, because why not. It’s a better morning’s work than submitting 258 papers to SSRN.

3. Setting up the model

This simulation is a little different than what came before. For my earlier post I simulated hypothetical future data from the beauty-and-sex-ratio study. For my new post I’m simulating data from a hypothetical strong relationship between two variables. Here I’m entirely focused on the question of what happens to a nonlinear (but monotonic) pattern when independent measurement error is added to the predictor. To make the patterns clear, I will assume a large effect; indeed, I’m only concerned here with the effect size, not its variation.

The one thing I’ll preserve from the previous simulation is the distribution of attractiveness measurements in the sample: 2%, 5%, 45%, 37%, and 11% in the five attractiveness categories.

I want to set up a model where there’s an assumed true nonlinear relationship between some outcome of interest and a continuous underlying latent attractiveness variable.

The nonlinear model isn’t hard to come up with. The hard thing is creating a distribution for the latent attractiveness, under the constraint that it be roughly consistent with the data.

I’ll assume a symmetric model, a normal distribution on the measurement scale. This might seem like a strong assumption but it seems reasonable enough given the data, which have an average that’s between 3 and 4. The tails are asymmetric but that’s just because the scale has 5 as an upper bound.

So I have a two-stage model:

x_latent[i] ~ normal(mu_latent, s_latent), for i=1,…,n
x_obs[i] ~ normal(x_latent[i], s_obs), for i=1,…,n
x_discrete[i] = round_and_bound(x_obs[i]), for i=1,…,n, where “round_and_bound” refers to rounding the continuous x_obs to the nearest integer between 1 and 5.

To simulate from this model, we need the hyperparameters mu_latent, s_latent, and s_obs—but with only one measurement per person, we can’t separately identify s_latent and s_obs; all we can estimate is sqrt(s_latent^2 + s_obs^2).

To separately identify s_latent and s_obs, we need to make some assumption about variation in the attractiveness measurement.

My first step is to compute the mean and sd of the observed data:

n <- round(2792*c(.02,.05,.45,.37,.11))
mean_ratings <- sum(n*(1:5))/sum(n)
sd_ratings <- sqrt(sum(n*((1:5) - mean_ratings)^2)/sum(n))

These come to 3.5 and 0.83. Assuming random zero-mean measurement error, the latent attractiveness distribution should have the same mean but a lower variance. We should subtract 1/12 to account for the rounding (as this is approximately equivalent to adding an independent uniform(-0.5,0.5) random variable, which has a variance of 1/12) and then more to account for measurement error.

But how much to subtract? For this, Attractiveness seems pretty subjective to me (also it can vary a lot over time, as we know from those high-school yearbook photos of people who later became movie stars). I gotta assume something here, so let's assume that the standard deviation of the measurements is 0.5. I'm picking this because 0.5 is a large value, but it's not as large as the sd of 0.83 in the data. Indeed, 0.5^2 + 1/12 = 0.33 and 0.83^2 = 0.69, so I'm assuming that about half of the variance in the data is coming from measurement error. If s_obs = 0.5 and the sd of the rounded data is 0.83, then s_latent should be approximately sqrt(0.83^2 - 0.5^2 - 1/12) = 0.6.

So here's a model, which I'll express in R:

round_and_bound <- function(a, lower=1, upper=5){
  pmin(upper, pmax(lower, round(a)))
}
N <- sum(n)
mu_latent <- 3.5
s_latent <- 0.6
s_obs <- 0.5
x_latent <- rnorm(N, mu_latent, s_latent)
x_obs <- rnorm(N, x_latent, s_obs)
x_discrete <- round_and_bound(x_obs)

Again, x_latent is the vector of respondents' latent attractiveness (as implicitly defined by the model), x_obs is the vector of survey interviewers' continuous perceptions of the respondents' attractiveness, and x_discrete is the recorded attractiveness on a 1-5 scale.

Let's check the results:

> print(round(table(x_discrete)/sum(n), 2))
x_discrete
   1    2    3    4    5 
0.01 0.09 0.40 0.40 0.10

OK, not perfect, but close enough. We'll get back to this model-fitting problem in a future post. For now, we'll go with the simulations we have here.

4. The prediction model

Now that we have a model of latent and recorded attractiveness, we can simulate some data.

Here's what we're gonna do.

First we'll suppose a linear relationship between latent attractiveness and some hypothetical outcome of interest. Then we'll look at what the corresponding relationship is between recorded attractiveness and that outcome. It should still be close to linear.

Next we'll do the same thing, but with a latent nonlinear model, and it turns out that this will be closer to linear. That is, adding noise to the measurement will make the E(y) vs. x relationship closer to linear.

We'll simulate the hypothetical data as follows:

x_latent and x_discrete are the 2972 attractiveness values simulated above. The latent values are all between 0 and 6.

There will be corresponding values y generated, first by specifying E(y|x_latent), then adding noise:

- In one simulation, we'll assume a linear relationship, E(y|x_latent) = x_latent/6. In the other simulation, we'll assume an S-curve, E(y|x_latent) = invlogit(2*(x_latent - 4)). We choose this curve because it is monotonic but strongly nonlinear in the range.

- My first thought was to add linear noise to the data. But I wanted to keep the range of the data between 0 and 1. Then the natural choice is binary data. But for this purpose it will be helpful to visualize with continuous data. So I'll take advantage of the constraint that E(y) is always between 0 and 1 and draw y from a beta distribution with 10 degrees of freedom, that is, y ~ beta(10*E(y|x_latent), 10*(1 - E(y|x_latent))).

Here's what we get:

The left column shows what happens when we assume a linear underlying relationship. As expected, the linear relationship is approximately preserved with the observed data.

The right column shows what happens when we assume a nonlinear underlying relationship. With the observed data, the relationship becomes much more linear.

In summary, adding measurement error attenuates the function (that is, the relation between y and x_recorded is weaker than the relation between y and x); also it brings a strongly nonlinear relation closer to linear.

This was not a surprise to me, but it's good to see it in a simulation.

# Add Health H3IR1 S35Q1 PHYSICAL ATTRACTIVENESS OF R-W3
# How physically attractive is the respondent?
# Taken from: National Longitudinal Study of Adolescent to Adult Health (Add Health), 1994-2018 [Public Use].
# in DS8

# Raw totals:  1.7% very unattractive, 4.9% unattractive, 45.3% about average, 36.7% attractive, 11.4% very attractive, N = 4877
# Section 24:  Respondent identification number (AID), Was the baby a boy or a girl (H3LB3), 1=boy (687 responses), 2 = girl (644)
# Kanazawa study:  N = 2972 ("Wave III respondents who have had at least one biological child")

library("arm")
library("cmdstanr")
linear_0 <- cmdstan_model("linear_0.stan", pedantic=TRUE)
linear_prior <- cmdstan_model("linear_prior.stan", pedantic=TRUE)
set.seed(123)

x <- seq(-2,2,1)
y <- c(50, 44, 50, 47, 56)
sexratio_data <- list(N=length(x), x=x, y=y)

display(lm(y~x))

fit_0 <- linear_0$sample(data=sexratio_data, refresh=0)
print(fit_0)
sims_0 <- fit_0$draws(format="df")

optimize_0 <- linear_0$optimize(data=sexratio_data)
print(optimize_0)
mle_0 <- optimize_0$mle()

fit_1 <-  linear_prior$sample(data=c(sexratio_data, mu_a=45.8, sigma_a=0.5, mu_b=0, sigma_b=0.2), refresh=0)
print(fit_1)
sims_1 <- fit_1$draws(format="df")

pdf("sexratio_simple.pdf", height=3, width=5)
par(mar=c(3,3,3,2), mgp=c(1.7,.5,0), tck=-.01)
plot(x+3, y, ylim=c(43, 57), xlab="Attractiveness of parent", ylab="Percentage of girl babies", bty="l", yaxt="n", main="From a survey of 3000 people",  pch=19, cex=1)
axis(2, c(45,50,55), paste(c(45,50,55), "%", sep=""))
dev.off()

pdf("sexratio_bayes.pdf", height=4, width=10)
par(mfrow=c(1,2), mar=c(3,3,3,2), mgp=c(1.7,.5,0), tck=-.01)
plot(x, y, ylim=c(43, 57), xlab="Attractiveness of parent", ylab="Percentage of girl babies", bty="l", yaxt="n", main="Least-squares estimate",  pch=19, cex=1)
axis(2, c(45,50,55), paste(c(45,50,55), "%", sep=""))
curve(mle_0["a"] + mle_0["b"]*x, col="blue", lwd=2, add=TRUE)
text(1, 49.2, paste("y = ", fround(mle_0["a"], 2), " + ", fround(mle_0["b"], 2), " x", sep=""), col="blue")

plot(x, y, ylim=c(43, 57), xlab="Attractiveness of parent", ylab="Percentage of girl babies", bty="l", yaxt="n", main="Bayes estimate with informative prior",  pch=19, cex=1)
axis(2, c(45,50,55), paste(c(45,50,55), "%", sep=""))
curve(median(sims_1$a) + median(sims_1$b)*x, col="blue", lwd=2, add=TRUE)
text(1, 45, paste("y = ", fround(median(sims_1$a), 2), " + ", fround(median(sims_1$b), 2), " x", sep=""), col="blue")
dev.off()

# Simulation

n <- round(2972*c(.02,.05,.45,.37,.11))
p <- 0.488
N <- 1000
x_sim <- array(NA, c(N, length(n)))
for (k in 1:N){
  x_sim[k,] <- rbinom(length(n), n, p)
}

pdf("sexratio_rep_1.pdf", height=5, width=8)
par(mfrow=c(4,5))
par(mar=c(2.5,3,.5,.5), mgp=c(1.5,.3,0), tck=-.01)
for (k in 1:20){
  plot((-2):2, 100*x_sim[k,]/n, ylim=c(43, 57), pch=20, xaxt="n", xlab=if (k>15) "Attractiveness of parent" else "", ylab=if (k%%5==1) "% girls" else "", yaxt="n", bty="l")
  if (k%%5==1)
    axis(2, c(45, 50, 55), c("45%", "50%", "55%"))
  else
    axis(2, c(45, 50, 55), rep("", 3))
  if (k > 15)
    axis(1, (-2):2)
  else
    axis(1, (-2):2, rep("", 5))
  abline(lm(I(100*x_sim[k,]/n) ~ I((-2):2)), col="blue")
}
dev.off()

# Attractiveness and measurement error

round_and_bound <- function(a, lower=1, upper=5){
  pmin(upper, pmax(lower, round(a)))
}

clumped <- function(a) {
  ifelse(a==1 | a==2, 1, ifelse(a==3 | a==4, 2, ifelse(a==5, 3, NA)))
}

mean_ratings <- sum(n*(1:5))/sum(n)
sd_ratings <- sqrt(sum(n*((1:5) - mean_ratings)^2)/sum(n))

N <- sum(n)
mu_latent <- 3.5
s_latent <- 0.6
s_obs <- 0.5
x_latent <- rnorm(N, mu_latent, s_latent)
x_obs <- rnorm(N, x_latent, s_obs)
x_discrete <- round_and_bound(x_obs)

print(round(table(x_discrete)/sum(n), 2))

# Assume a relationship

linear_curve <- function(x) {
  x/6
}
s_curve <- function(x) {
  invlogit(2*(x-4))
}
response_curve <- list(linear=linear_curve, nonlinear=s_curve)

pdf("latent_error.pdf", width=6, height=4.5)
par(mfcol=c(2,2), mar=c(3,3,2,2), mgp=c(1.2,.2,0), tck=-.01)

for (k in 1:2) {
  p <- response_curve[[k]](x_latent)
  y <- rbeta(N, 10*p, 10*(1-p))
  plot(x_latent, y, ylim=c(0, 1), xlab="Latent beauty", ylab="Outcome", yaxt="n", pch=20, cex=.1, main=paste("Latent", names(response_curve)[k], "model"), cex.axis=.9, cex.lab=.9, cex.main=1)
  axis(2, c(0, 0.5, 1))
  curve(response_curve[[k]](x), col="white", lwd=3, add=TRUE)
  curve(response_curve[[k]](x), col="blue", lwd=2, add=TRUE)
  plot(x_discrete + runif(N, -.1, .1), y, ylim=c(0, 1), xlab="Recorded beauty (jittered)", ylab="Outcome", yaxt="n", pch=20, cex=.1, main=paste("Expected values given recorded beauty"), cex.axis=.9, cex.lab=.9, cex.main=.8)
  axis(2, c(0, 0.5, 1))
  Ey <- rep(NA, 5)
  sy <- rep(NA, 5)
  for (i in 1:5) {
    in_bin <- x_discrete==i
    Ey[i] <- mean(response_curve[[k]](x_latent[in_bin]))
    sy[i] <- sd(response_curve[[k]](x_latent[in_bin]))
  }
  lines(1:5, Ey, col="white", lwd=3)
  lines(1:5, Ey, col="red", lwd=2)
}
dev.off()

Postdoc hiring season at Flatiron

Posted on by
1

Every year, in the Center for Computational Mathematics at Flatiron Institute, we hire 7 to 10 postdocs to keep up our steady state of roughly 20 postdocs (more and more of them are leaving after 1–2 years to go to industry). Here’s the job ad for this year, with applications due on the absurdly early date of November 15, 2026 (we’re still playing catch-up with plans that make postdoc offers early December):

The 12-month salary is US$95K plus a US$10K research budget, way better compute resources than most universities have (e.g., hundreds of H100 GPUs which you can actually use), and exceptional benefits (not only free lunch and breakfast, but also great medical care and retirement plans, which has extreme variance among employers in the US). The positions are three years (technically two, with a one year renewal, but nobody has ever been denied a renewal in the 6 years I’ve been here).

Our center is split into three primary research areas, each of which has a web page describing the people and research areas:

I’m here, as are Steve Bronder and Brian Ward, two software engineers who spend a lot of their time working on Stan, along with other great software engineers (like Jeff Soules and Jeremy Magland who built the Stan Playground with Brian and also MCMCMonitor), and perhaps most importantly, a really engaged and lively group of postdocs working on inference (largely from a diffusion or normalizing flow perspective). This place is great for collaboration both internally and externally.

The postdocs are really fellowships—it’s not like American academia where I get a grant and hire you to work on the grant. You can really work on whatever you want that’s on mission as a postdoc here. So let me recall our mission, which actually matches what we do:


The mission of the Flatiron Institute is to advance scientific research through computational methods, including data analysis, theory, modeling and simulation.

Of all the places I’ve worked over the last 40 years, Flatiron is by far the best. There’s not even a close second. And I’ve been lucky in having great jobs with top notch colleagues (prof at CMU, researcher at Bell Labs, industrial researcher and software engineer writing production code, research scientist at Columbia, then here).

We’re also in a great location—the Flatiron neighborhood of New York City, which is a 10 minute walk to NYU, 10 minute walk to Google and Meta, and 20 minute subway ride to Columbia. Baruch College, the CUNY graduate center, and Fordham’s Manhattan campus are all nearby, and Yale, Rutgers, and Princeton are all in day-trip distance with lots of collaboration.

We’ve had great postdocs, so perhaps not surprising they’ve gotten good jobs: industrial ML jobs at Anthropic, OpenAI, Nvidia, and DataBricks, and stats faculty jobs at Duke University, Johns Hopkins University, the University of Texas, and the University of British Columbia. And that’s just the ML/stats side of our postdocs.

If you want pretty pictures of our setting, check out my job ad from 2021.

If you’re going to be applying in the area of Bayesian modeling or inference, please contact me directly so I don’t miss your application (we got nearly 300 applications for postdocs last year!): [email protected].

Confirmation bias in software testing, Walnutpie edition

Posted on by
3

So, those exciting results I presented at StanCon about Walnutpie? Let’s pull back the curtain.

Confirmation bias in testing

Those stellar results were the result of a bug in the way we calculated R-hat and effective sample size in the interfaces. When we plug in the bug fix (as we have for the 0.0.3 release on PyPI), the really stellar results on a collection of three models from Andrew shrinks and in some cases goes away entirely. So why didn’t we catch the bug? Because the results looked good but not unreasonable. Low R-hat and a decent ESS at about a tenth of the number of draws.

Claude’s attention to detail

We should have been more attentive. When I couldn’t replicate the results in our latest interface, I fed the results into Claude. Being much more attuned to detail, it immediately noticed that the table we’d provided was inconsistent. Can you see the inconsistency?

=== hierarchical_matrix_2.stan ===
chains=4  draws/chain=[1000, 1000, 1000, 1000]  total=4000  wall=0.6s
param                  mean        sd     mcse      ess    ess/s    rhat
s_y                  0.5551    0.0316   0.1193    161.7    273.2  1.0003
mu_b.1               0.9542    0.0921   0.0784    279.0    471.3  1.0001
mu_b.2               1.9781    0.0902   0.0750    302.4    510.8  1.0001
mu_b.3               2.1326    1.5192   0.0779    314.9    532.0  1.0003
...

Hint: Look at mu_b.2 and mu_b.3. Or focus on rhat vs. ESS.

Second hint: Look at the MCSE estimates of mu_b.2 and mu_b.3 and compare to the standard deviations and ESS values. All I can say is, “D’oh!” How can R-hat be so close to 1 with such low ESS values? Let’s sing it again together, “D’oh!”

I spend my life bouncing between “D’oh” and “Yatta!“.

What’s wrong with (Wal)nutpie?

It seems that like Nutpie, the use of gradients (specifically the diagonal of a regularized outer product of gradients to initialize the mass matrix) can make the initialization very sensitive. The problem is easy to see even with a multivariate standard normal (i.e., one located at zero with a unit covariance matrix). With any normal, the Hessian (second derivative matrix) is equal to the negative inverse covariance at every point at which it is evaluated. Fixed curvature is relatively straightforward to adjust with a mass matrix (equivalently, a linear preconditioner).

With a standard multivariate normal, the covariance is a unit matrix, so the Hessian’s just a negative identity matrix. This is what we’re trying to estimate as the mass matrix and it’s the default initialization for Stan.

The gradient is a different story. At position x in a multivariate standard normal, the gradient is -x. With a unit step size, a move along the gradient takes us back to the origin in one step (x + 1 * -x = 0). Nutpie regularizes the diagonal of the outer product, diag(x * x'), by taking a geometric average with a unit matrix. This effectively multiplies by one, then applies a square root, giving us an initialization of sqrt(diag(x * x')) = diag(abs(x)). This no longer takes us back to the origin in one unit-step-size step, but it does grow without bound. As the draws move closer to the bulk of the probability mass, this initial mass matrix is a very bad estimate which will then drive down step size.

Fool me twice, shame on me

I’m sorry to say that I’ve done this once before with the evaluation of an information extraction algorithm in NLP. I was working in industry on building named entity finders (e.g., people, places, things) in text that would have very high sensitivity. We allowed specificity to go down accordingly as forced by the ROC curve. This was some really fun n-best algorithm work for hidden Markov models (HMMs) and conditional random fields (CRFs), which was in our LingPipe software. This got us disqualified from the United States National Institute of Standards (NIST) bakeoff because we’d have won the sensitivity evaluation by a landslide. In the mid 2000s when this was going on, NLP and ML couldn’t get themselves away from the first-best-guess paradigm of evaluation. We continued to get funding because we aced the evaluations by real-life intelligence analysts, who were willing to sacrifice specificity for sensitivity.

So what was the bug? I coded what should have been a multiset evaluation as a set. It inflated our sensitivity at an acceptable (to the analysts) level of specificity from something like 99% to 99.999%. Oops. We came clean of course, just as I’m doing here.

A record of getting disqualified

This happened to me before with the United States Defense Advanced Research Projects Agency (DARPA), who disqualified our team at SpeechWorks in 2000 or so, which had built a spoken dialogue system that actually worked in production to book airline flights (it was in production at Delta for internal use). We were supposed to be the industrial “pacing horse” that would show just how much better things could be built using what they called “mixed initiative dialogue” in a research lab. The problem is, they got research lab members to build prototypes of an idea that was flawed linguistically from the get go. The flaw was assuming someone could build or would want a system where you could say “I need to fly to Dallas Tuesday morning, then that night go to Atlanta for two days, then back to Los Angeles for a night, then onto New York City, but I’m not sure how long I’ll need to stay in New York City. And by the way, I prefer aisle seats.” That’s just not how human dialogue works. Also, my advice is, don’t bring an academic prototype to a bakeoff with industry in an area they’ve launched into production.

OpenAI tries to one-shot human dialogue

Ironically, OpenAI made what I consider to be the same mistake with its post-training that focused on one-shot question answering. It trained the system to always try to provide an answer, like a student facing a closed-book exam with no penalty for guessing.

Then every other system just copied ChatGPT’s behavior. Given their corporate nature and the literally billions, if not trillions, of dollars involved, OpenAI and Anthropic have been complaining about IP theft. I’m laughing too hard to cry. To bring it back home to the topic of this post, the problem is an ill-chosen initialization. As an aside, did that initial training set include a hefty admixture of em-dashes, triple-adjective combos, and positive/negative contrasts because someone thought that’s what good writing was? I’m not suggesting hiring “theoretical” linguists, but someone who knew something about writing and the structure of actual human dialogue might have been helpful. Now all the AI companies are opening up the engines to search, which helps with the closed-book part of the exam, and just yesterday, Claude stopped in the middle of “thinking” (see, scare quotes for the haters) and asked a follow-up question before going on. Yay, progress! Let’s hope I can do the same on a much smaller scale with Walnutpie.

Here are the talks from StanCon 2026!

Posted on by
Reply

StanCon 2026 just happened!

And here are the talks:

Matthew Kay
Adrian Seyboldt
Paul-Christian Bürkner
Charles Margossian
Javier Enrique Aguilar
Sean Pinkney
Nikolas Siccha
Anna Dreber
Kaitlyn Johnson
Jonas Wallin
Pranav Sanke
Anna Elisabeth Riha
Colling Cademartori
Tim M. Szweczyk
Bob Carpenter
Chandler Ross
Nils Rudi
Fredrik Ronquist
Jakob Torgander
Aleksi Lahtinen
Ville Laitenen
Zeno Romero
Soham Mukherjee
Steve Bronder and Brian Ward
John Ashley Burgoyne

Lots of great stuff here. Check out the titles of the talks—all sorts of different topics.

Time to start planning for StanCon 2027.

7 1/2 Schools—mistakes were made

Posted on by
21

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

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.

Bayesian computation meeting in France in May, 2027

Posted on by
Reply

OK, this one’s time constrained so I’ll post it right away, not on the usual lag.

Christian Robert points to this conference announcement:

We invite proposals for the BayesComp 2027 mirror event in Aussois, France, that will take place on 17-21 May, 2027. The event will complement the main BayesComp 2027 meeting in Texas and provide a European forum for researchers in Bayesian computation and related areas. Plenary sessions will be broadcasted from Texas, while talks in Aussois will be recorded.

Please submit one form per proposed contribution. Proposals may be for contributed talks, or poster presentations. The programme committee will review submissions based on scientific quality, timeliness, coherence, and expected appeal to the BayesComp community. Unless the submission mentions otherwise, talk proposals that are rejected will be automatically added as posters.

Details about the venue (CAES Paul Langevin), potential registration fees, and technical details are being finalized and will be announced on the event website.

Looks like fun! And Bayesian computation is important.

Walnutpie version 0.0.1 Released

Posted on by
6

We are happy to announce the official release of Walnutpie version 0.0.1.

Walnutpie is an MCMC sampler for continuously differentiable densities coded in Python, accepting models coded in Stan, PyMC, NumPyro, JAX, and plain old Python.

Walnutpie is not an official Stan project

I thought this was worth saying up front. It may eventually migrate to Stan, but for now, we followed the Nutpie approach of building a standalone sampling package that worked with a range of packages for defining models.

R version

We plan to develop an R interface after we release version 1.0.0 of the Python interface. So hopefully in 2026.

pip installable

Walnutpie is on PyPI, so it’s pip installable. The documentation includes information on getting started, running models, and posterior analysis. We have not yet included case studies for modeling tools other than Stan and Python.

Stan through C++

Walnutpie runs Stan models through C++ using BridgeStan, so there is no Python dispatch overhead for Stan sampling. The basic architecture of the API is based on Adrian Seyboldt’s sampler Nutpie. We are working on doing that for other packages like NumPy and JAX to the extent that we can.

GitHub source

Development discussions and source code are managed through GitHub.

Features of Walnutpie

We are almost ready to release the paper on arXiv with all of the gory pseudocode details of all of the algorithms used for Walnutpie. It will explain the following points in detail.

  1. Walnuts: The underlying Hamiltonian Monte Carlo sampler is Walnuts. Walnuts uses Nuts for choosing the number of steps per iteration. It further allows step sizes within the Hamiltonian dynamics simulation to be lowered when necessary to preserve simulation accuracy. This helps with robustness and with accuracy in multi-scale distributions (i.e., ones where the curvature as represented by the Hessian varies around the posterior). With a high tolerance threshold, Walnuts reverts to Nuts’s behavior.

  2. Mass-matrix warmup: The mass-matrix warmup strategy is an online form of Nutpie (links to: the paper and software). Nutpie minimizes Fisher divergence by estimating the inverse mass matrix as the midpoint (in the appropriate manifold) between an estimate based on the variance of the draws and the covariance of the scores (gradients of the log density). The target is better than Nuts’s variance of draws in both convergence speed and sampling efficiency. Walnutpie only supports diagonal mass matrices (Nuts supports dense matrices and Nutpie supports low-rank plus diagonal and even more general normalizing flow approaches). The approach is online in the sense that it is not blocked like warmup in Nuts or Nutpie—it updates every iteration by exponentially discounting the past to mimic Stan’s exponentially increasing history sizes. We also borrow Nutpie’s mass matrix initialization based on a regularized outer product of gradients at the initial point.

  3. Step-size warmup: The step size adaptation strategy has not changed, but the underlying stochastic gradient descent algorithm is Adam rather than dual averaging. We found Adam to be faster to converge and much more stable. Matt Hoffman included a hack in the original Nuts approach to stabilize dual averaging, but even with that it is not as stable as Adam.

  4. Concurrency and automatic stopping: The underlying sampler is multi-threaded (using C++11 threads) with shared data. On top of the multi-threading, we have layered a convergence monitor in a separate thread that communicates with the chains through lock-free, latest-only, single-producer/single-consumer (SPSC) buffers (specifically, a triple buffer). The monitor automatically stops warmup when the mass matrices and step sizes have converged within tolerance to their cross-chain averages. The monitor automatically stops sampling when a target (traditional, non-split, non-ranked) R-hat; threshold is satisfied for the unnormalized log density, which typically converges more slowly than any of the individual parameters. It can also be configured to run for a fixed number of warmup and/or sampling iterations. The link between the original R-hat and effective sample size makes this essentially an unscaled ESS target.

  5. Ragged chain summaries: Asynchronous concurrent execution of chains with automatic stopping produces chains of different lengths. Because ArviZ does not accept ragged chain input of this kind, we have included posterior analysis tools for means, variances/standard deviations, quantiles, traditional R-hat, effective sample size, and Monte Carlo standard error that work with ragged chains.

  6. C++20: Walnutpie is implemented in C++20. As a programming language type fanatic (I’ve written two books with “type” in the title!), I don’t know how I survived without C++ concepts before C++20.

  7. ctypes FFI: The foreign function interface in Python uses ctypes rather than a higher-level interface, which sidesteps the requirement of ABI compatibility of C++ binaries.

Developers

We would also like to welcome new developers who may want to get involved. There are already a stack of improvements we’d like to make, which we have enumerated on the GitHub issues.

Sources of algorithms

The Walnuts algorithm was a joint effort among Nawaf Bou-Rabee, Sifan Liu, Tore Kleppe, and Milo Marsden. Nutpie was developed by Adrian Seyboldt. Nuts, in the form used currently in Stan, was originally developed by Matt Hoffman and Andrew Gelman, then improved with multinomial sampling and mass matrix adaptation by Michael Betancourt. We haven’t yet added cross-chain adaptation as developed by Ben Bales, but the pieces are all in place to do so.

Brian Ward and I wrote all of the version 0.0.1 code with design advice and code review from Steve Bronder. Claude (the LLM) helped with code review and testing, but we didn’t use it to write the actual code (not out of principle, but because Brian and I both prefer the control of doing things manually).

Feedback

We would very much appreciate any feedback people have, including feedback on the documentation and ease of use, the source code, and performance.

We are happy to get feedback through issues on GitHub, through replies to this post, or through mail to one of the developers.

He fit the same statistical models with three different software and got much different estimates. It’s another dimension of the multiverse.

Posted on by
16

Scott Cunningham writes:

You’ll appreciate this I think. I ran Claude code on 96 specs for a popular difference-in-difference estimator with the identical specifications, ranging covariates only, for three languages (R, Python and Stata) and 2 packages for each. Different estimates — sometimes five times off. Appears to be coming from how each one handles failures on logits estimation with perfect separation, with some just refusing and others refusing and passing on garbage and others dropping down to run a logit on a constant despite being told to use Xs. And it has to do with how each package handles “big numbers” and this machine epsilon business that I barely know anything about. The fact that economists have pointed on selection on observables means they also haven’t invested in basic diagnostics around estimating propensity scores. But even turning Xs into z-scores still didn’t fix it. It’s the most extreme version of the garden of forking paths I think because people just assume languages and packages will all agree.

Ahhhh, perfect separation in logistic regression . . . one of my favorite topics! For an updated take, see section 14.6 of Regression and Other Stories. It makes sense that different packages will handle this differently.

Here’s Scott’s summary from his above-linked post:

The data come from a staggered rollout of mental health centers across Brazilian municipalities, 2002-2016. I have 29 covariates covering demographics, economics, health infrastructure, and population trends. Sixteen specifications sample the space from zero covariates up through all 29. All variation is coming from covariates, both within a package (i.e., which covariates) but more importantly across packages. And that’s the headline — even for identical specifications, you can get variation in estimates and sometimes as much as anywhere from an ATT of zero to an ATT of 2.38. Even if you drop the zero, some specifications yielded estimates ranging from 0.45 additional homicides to 2.38 homicides (measured as a municipality population homicide rate).

A five fold difference in estimated effects of a mental health hospital changing in the homicide rate is not trivial.

I guess one message is that there’s only so much you can learn from this sort of observational dataset. Ultimately it’s one data point in a larger explicit or implicit meta-analysis. More specifically, I like that Scott starts with the basic comparison and moves from there. I always think it’s a mistake to try to jump to the end.

Also I’d like to see a scatterplot with one dot for each municipality-year, plotting the outcome vs. some pre-treatment level, with different colors corresponding to treatment or control. I know that’s not the point of Scott’s post, but statistical packages can do graphs too! And more graphs, I’m not sure exactly what, but some way of getting a sense of the data. Ultimately this is being framed as a treatment vs. control comparison, so I want to see what the treatment and control data look like, before and after.

The multiverse!

Relatedly, Jessica points to two recent papers proposing the use of automated data analysis programs to perform multiverse analysis:

Many AI Analysts, One Dataset: Navigating the Agentic Data Science Multiverse, by Martin Bertran, Riccardo Fogliato, and Zhiwei Steven Wu

Editorial Screening when Science is Cheap, by Nic Fishman and Gabriel Sekeres

Jessica has discussed some related ideas at the end of this post and also here. And see here for lots more.

I agree with Scott and Jessica that these ideas all fit together somehow.

Reviews of our Bayesian Workflow book from Bin Yu, David Spiegelhalter, Brad Efron, Christian Robert, Rohan Alexander, and Mine Doğucu!

Posted on by
1

Roughly speaking, Bayesian Workflow is to Bayesian Data Analysis in 2026 what Bayesian Data Analysis was to earlier Bayesian books in 1995: it builds upon everything that came before.

With Bayesian Data Analysis, the big steps forward were:

  • Going beyond Bayesian inference to also consider Bayesian model building (as a researcher, you construct the model, it isn’t just given to you as in a textbook), model checking (breaking through the absolutely horrible attitude, common to Bayesians in the early 1990s, that the model was “subjective” and thus should not be checked), and model improvement (continuous model expansion, not the misguided idea of assigning posterior probabilities).
  • Going beyond simple conjugate models. BDA had lots of hierarchical models, also lots of computational tools so that you could fit the models you want by putting them together from understandable components. And I like how we had a clear separation between modeling and computing. The model comes first, then you figure out how to compute it. Or you set up a model that works within your computational constraints.
  • A Bayesian approach to sampling and causal inference. This was Rubin’s framework in which unobserved units in the population and unobserved causal outcomes are treated as missing data and are part of a joint probability model. We worked this out in chapter 7 of BDA (which became chapter 8 in the third edition of the book).
  • Lots of live examples. Not just “real-data examples,” but problems we’d directly worked on. This motivated us and I think it gave our readers a sense of how Bayesian methods worked not just in theory but in applied problems.
  • A pragmatic view of probability as a measurable quantity. That’s right there in chapter 1. Bayesian methods are not the product of a philosophical stance; they’re a way to connect models and data using probability.

I could go on and on, but for that I can refer you to the Bayesian Data Analysis book.

And these are the key innovations of Bayesian Workflow:

  • Going beyond Bayesian data analysis (model building, inference, model checking, and model expansion) to consider the larger process of statistical modeling, including comparisons of multiple models fit to a single dataset.
  • A fuller use of informative priors. This is a big deal. In BDA we still had a bit of the Bayesian cringe going on. One reason we’ve moved toward stronger priors is that the replication crisis has taught us that the amount of prior information available in any given problem is often approximately the same as the information coming from an experiment (see here, for example). Informative priors also fit our increased focus on generative modeling, and we’re doing a lot more prior predictive checking to understand the implications of our models.
  • More integration between modeling, data analysis, and computing. One way to see this is that the Bayesian Workflow webpage has the code to run all our examples. We also have lots of code snippets in the text as a way of demonstrating the way in which coding is central to our statistical workflow.
  • Lots more live examples. It’s been 30 years since BDA first came out. One reason that Bayesian Workflow has 11 authors is that different collaborators worked on different examples (but the three principal authors read through the entire book, so the general approach should remain coherent).
  • Simulation-based experimentation. This is something my colleagues have been doing more and more over the years. At its most basic, simulation-based experimentation provides a best-case baseline for statistical methods: if you can’t recover your quantities of interest with sufficient accuracy under ideal conditions (when your data are simulated from the model you’re fitting), then you know you’re in trouble. And often this is the case! Beyond that, we can simulate from one model and fit another, and see what happens. Simulation experiments aren’t always so easy to construct, as they involve specifying the entire data-generation process. But we think this is effort worth expending, as it involves thinking about the problem you’re working on.

I could go on and on, but for that I can refer you to the Bayesian Workflow book.

And now for the reviews

But you don’t have to trust me on this! Just listen to some of the eminent statisticians and educators who’ve reviewed our book:

Bin Yu (University of California):

An outstanding, protocol-driven guide for Bayesian data analysis, Bayesian Workflow by Gelman, Vehtari, McElreath and co-authors delivers a practical and comprehensive framework for iterative modeling, emphasizing simulation, diagnostic checks, and rigorous empirical validation, and with a long and impressive list of case studies. By treating data analysis as a structured, verifiable workflow, it provides an indispensable toolkit for diagnosing model failures, refining priors, and building reliable data analysis systems for reproducible conclusions, useful for beginning and veteran data analysts alike.

David Spiegelhalter (Cambridge University):

This is not a typical methods textbook, but instead it guides the reader through the whole process of fitting, critiquing and adapting statistical models to real-world problems. It is full of the accumulated wisdom of skilled practitioners, teaching through demonstration rather than theory, with both basic and highly sophisticated examples. I strongly recommend this book to statisticians who really want to understand what they can learn from their data.

Brad Efron (Stanford University):

A bravura performance…Gelman, Vehtari, McElreath and friends develop in detail a practical Bayesian data analysis workflow, from acquisition to final report, including full computational guidance.

Christian Robert (Université Paris Dauphine):

This original, thought-provoking, and transformative book is much much more than an implementation manual for Bayesian Data Analysis, even though it shares almost the same perspective. (The first sentence of the book states that the authors’ “conceptions of statistical practice, and of Bayesian statistics, have changed over the years”.) By providing a modus vivendi for undertaking Bayesian modelling from scratch in realistic settings where models are not magicked out of the blue, the authors explicit and rationalise the many steps required by such a bottom-up modelling protocol (“not a checklist, not a cookbook”, and not a flowchart!) in real situations. The contents read very well and very smoothly, with a seamless conjunction of intuition, modelling advices, computational details, and comparison tools. While unsurprisingly Bayesian, the perspective adopted therein remains both open and inclusive, with a welcome humility about the limitations and challenges of Bayesian workflows. This book should thus appeal to and profit a wide variety of readers, as providing guidance through an extensive collection of highly detailed examples, with shared code and exercises.

Rohan Alexander (University of Toronto):

Some statistics books show you how to beat an egg, others are recipe books: if this, then that style. This book teaches you how to cook. Written by authors who established so much of how we do Bayesian statistics, this new book is an indispensable guide for analyzing data in a trustworthy way. It walks you through the actual steps involved in building models to explore and understand datasets. Part 4 is particularly excellent – the authors provide many end-to-end case studies that will be useful for both practitioners and students. It highlights the value of their workflow-based approach. Filled with chatty asides, the book introduces the Bayesian workflow to a broad audience. It embraces the frustrations and complexities of actually doing Bayesian statistics and provides specific guidance throughout. Each chapter contains exercises and it could be the basis of an upper-year undergraduate course, or a first-year grad course, in applied statistics. It will be used for many years to come.

Mine Doğucu (Harvard University):

What makes Bayesian Workflow so exceptional is how it seamlessly pairs profound ideas about modeling with the adoption of modern computational practice. By centering the messy, iterative process of modeling through real-world case studies, the authors reject rigid cookbooks and checklists in favor of building deep situational awareness. Because the ideas are so clearly articulated and deeply applied, this book serves as an invaluable pedagogical resource. With its practical exercises, individual chapters or the text as a whole can easily be integrated into upper-level undergraduate or graduate courses, while also remaining accessible for self-guided readers. It is an indispensable read for anyone with foundational knowledge in Bayesian methods, regardless of whether they are applied practitioners, software developers, or methodologists.

You might also be interested in the journal issue on statistical workflow that we recently edited for the Philosophical Transactions of the Royal Society.

Again, here’s Bayesian Workflow on Amazon, here’s the publisher’s website, and here’s our website with data, code, and lots more.

Enjoy.

18 Associate Editors resign from Statistics and Computing editorial board: Problems with commercial scholarly publishing, and what does this all mean?

Posted on by
29

I was cc-ed on a message sent by 18 members of the board of the journal Statistics and Computing, quitting their posts because the publisher (Springer) has announced a new policy whereby all authors will have to pay publication charges. The soon-to-be-former associate editors write, “Statistics and Computing will no longer publish the best science, both due to financial exclusion of those researchers who cannot afford to pay, and those community-minded researchers who refuse to pay on principle.”

I’ll put the full message, with its 18 signatories, below the fold.

My reaction to all this is that it would be great if the journal could move to an open and free system such as is done by the Journal of Machine Learning Research–a journal that I believe was founded by people who had resigned from the editorial board of a commercial journal.

Even commercial journals that begin with good intentions can develop fatal problems. For example, check out the sad story of the Berkeley Electronic Press, a set of commercial journals that was founded by a friend of mine. My friend’s an economist, and I guess he might say that it was the iron logic of capitalism that reduced a once-noble endeavor to a rent-seeking enterprise.

So, yeah, I’d recommend that Statistics and Computing follow the path of JMLR, really try to imitate its structure as closely as possible. Bayesian Analysis is another free journal that appears to run with minimal overhead.

It kind of bugs me, though. Profit-making companies have done great things in publishing and communication. A quick glance at our shelves reveals lots of wonderful books, almost all of which were published privately. As were Bayesian Data Analysis, Regression and Other Stories, Active Statistics, and the rest of my books. Lots of great movies are made for money too.

On the other hand, Arxiv is nonprofit, as is lots of the web, on which I post all my published and unpublished papers. Wikipedia is nonprofit, and this blog is written using WordPress, which appears to be another nonprofit organization. I teach at Columbia University, which is private but nonprofit, not run perfectly by any means but still going strong.

Scholarly publishing is a funny industry because it’s my vague impression that it started out as a low-budget noncommercial enterprise, and then some private rent-seekers moved in. I guess these companies were doing something special or they wouldn’t have been so successful at taking over, but now it seems to have gone too far. More sites like Arxiv, JMLR, and Bayesian Analysis would be a good thing. Right now we’re always having to figure out where to publish our papers; it’s just an absolute mess. These journal submission websites make the Department of Motor Vehicles office look like a lean machine by comparison.

P.S. Retraction Watch ran a story on this, where they quoted Robin Ryder as saying, “If the editors regroup elsewhere to form a new journal, they hope to publish with a society, Ryder told us, citing journals like Journal of the Royal Statistical Society and Annals of Statistics, ‘none of which force authors to pay APCs.'”

I think it would be a mistake for them to follow the Journal of the Royal Statistical Society and Annals of Statistics. Both these journals have arduous paperwork-laden submission processes and both charge for access.

If you’re going to start over, why not follow the model of JMLR and Bayesian Analysis and make it all free? Cut out the middleman entirely!

Continue reading →

Claude builds 3D Hamiltonian Monte Carlo animation in one shot with anaglyphs

Posted on by
3

This post is from Bob

The sausage

So as not to bury the lead (or “lede” if you want a mid-20th-century newspaper vibe), check out the this 3D HMC animation generator.

It can render regular animations or produce anaglyph 3D encoding (red/blue). Unless you have 3D glasses, unclick the “Anaglyph 3D” checkbox at the bottom of the upper left corner control box.

The app let you zoom in and rotate the visualization with obvious controls (explanation in the footer of the visualization). The app also lets you adjust the amount of correlation in the 3D normal distribution as well as step size, number of steps, and animation speed. Looking the long way down a highly correlated “cigar” shape is dramatic.

The 3D effect with glasses is strongest when you rotate the visualization (it’s the usual intuitive controls with instructions at the bottom of the web page) and zoom in a bit. I find that using low 3D depth looks the best. Don’t get your hopes up too much. This isn’t Dr. Strange creating buildings in 3D in a Marvel movie.

If you want to pop it up in an independent browser so you can go to full screen, here’s a link.

How the sausage was made

I continue to be amazed at the progress of the frontier LLMs. The demo above was the result of handing Claude Opus 4.8 (“hard” thinking mode) the following single prompt with no build up. As with the Galileo inclined plane case study I posted, which Opus one-shotted, I was expecting some back and forth and false starts.

I want to generate a 3D animation for red/blue glasses of the Hamiltonian Monte Carlo algorithm. There is a nice online visualizatuion by Chi Feng here, but it is not 3D https://chi-feng.github.io/mcmc-demo/app.html I just want the main animation—no need to calculate marginals, etc.

To start, we can use a 3D highly correlated (0.9) normal target with unit variance aligned at one end of the cigar (e.g., near (2, 2, 2) looking toward (-2, 2, 2), which will have things zoom over your shoulder and come back).

If you can generate it so that it’ll run in a web browser with controls on step size and number of steps that’d be great, but if not, choose a step size conservatively so it won’t be rejecting very often. I want it to continue multiple iterations in order to see the effect of random momentum on the trajectories. Leave balls behind wherever the sampler actually samples. When it rejects, make the ball bigger. The trajectory should be thick enough to be visible.

If it’s easier to have Python generate an animation that’s also fine. I just want to be able to render it on my desktop to show people during a talk. I just ordered 50 pairs of cardboard red/blue 3D glasses to hand out.

I was wrong. It did it in one shot. After about 10 minutes of cranking away, it produced what you are looking at. The output is a self-contained (i.e., encapsulated) HTML file of 627KB. There are some things I’d change in an iteration (smaller pipes, fewer of them lying around), but I think it’s worth sharing the output of such a simple prompt. Perhaps needless to say, a follow up prompt gave me the HTML I needed to embed the result in this page as an iframe.

I wrote all 692 words of the blog post myself (other than the html embedding), but I’m sure Claude could have done that, too. The LLMs have fewer rhetorical tics when writing technical and scientific material. But it wouldn’t have sounded like me.

Statistical visualization in the mid 2020s?

I wonder what Andrew’s statistics visualization class would look like in 2026 with LLM-powered visualizations this easy to make. Now that the LLMs can reliably one-shot something this complex, I’m finally starting to worry about the future of programmers. Undergraduate enrollments in CS are very volatile and already going back down as they did after the dot com bubble burst. There was huge growth (a factor of two to three) from after the mortgage market bubble burst around 2007 until it started to decline again due to AI.

The high cost of split R-hat

Posted on by
1

This post is by Bob.

I’ve been thinking a lot lately about R-hat given that I’m using it for online converging monitoring in our new Walnuts implementation. In that setting, where I use Welford accumulators to update R-hat estimates every iteration, I can’t use split R-hat without way too much buffering. So I’ve been thinking about the effect of splitting, too, and whether we need it. I asked Andrew and he said Kenny Shirley once produced an example where split R-hat diagnosed non-convergence that regular R-hat didn’t, but that example is lost to time and we’ve never seen this kind of behavior with NUTS as far as I know (please give us an example in the comments or via email to Andrew if you have).

Relating R-hat and ESS

My intuition was that we could set a low enough R-hat threshold that it would ensure a high enough effective sample size (ESS) when we crossed it. The relation’s a little tighter than I thought, with

    Rhat^2 ≈ 1 + M / ESS,

where M is the number of chains and ESS is the effective sample size of all chains combined. There’s a multivariate proof in Vats and Knudson, 2021, Revisitng the Gelman-Rubin diagnostic, Statistical Science, page 2 and section 5 for details, but it’s pretty straightforward to get the intuition when you reduce R-hat^2 to (N-1)/N + var(chain-means) / man(chain-variances) as Charles Margossian did in his nested R-hat paper. Vats and Knudson disapprove of Andrew and Aki’s suggested threshold of 1.1 from BDA3, because it is satisfied with a combined ESS of 20 across Andrew’s default 4 chains.

Being me, I tried to validate my intuition with simulations rather than linear algebra. Also, I like to see that things work in practice that theory entails to make sure I’ve understood all the assumptions baked into the theory (one can’t prove anything without assumptions!). When asked to code a simulation using ArviZ, Claude inserted a (2 * M) in the numerator in place of the M. Where did that come from, I asked? It told me it needed the factor of 2 because ArviZ uses split Rhat. D’oh! Of course it does, because we’ve doubled M without increasing ESS.

A worked example

Suppose we have 4 chains with a combined ESS of 400. Then sqrt(1 + 4/400) ≈ 1.005 and sqrt(1 + (2 * 4) / 400) ≈ 1.01. We’ve effectively doubled the number after the 1 by splitting. Unlike Vats and Knudson, I usually don’t need an ESS >> 100, so the 400 required for split R-hat < 1.01 is perhaps a bit too conservative for my tastes. On the other hand, we face a practical problem estimating ESS reliably with fewer than 50 or so ESS per chain. Estimation is challenging because it relies on autocorrelation estimates from the chains themselves, which become much noisier when based on shorter chains. (Side question: Do we not combine autocorrelation estimates across chains to reduce standard error because some chains might not be mixing?) Also, we know this algebra wasn't a coincidence of 4 chains and 400 draws. The Taylor expansion of sqrt(1 + x) is the convergent sequence

    sqrt(1 + x) = 1 + x/2 - x^2 / 8 + x^3 / 16 + ...

When x < 0.1, the first-order approximation, sqrt(1 + x) = 1 + x / 2, is good.

The bottom line for practitioners

We need around twice as many draws to get below a fixed threshold with split R-hat than with the original R-hat.