Gaurav Sood writes:

I was thinking about Bob’s post, and I think we can do something with Wasserstein Barycenter Priors

I have the Gaussian thing going here.

But the right approach is the barycenter priors ….

I haven’t been following the details here, but I think this discussion is about a problem that I’ve seen many times in applications, which is that you want to update your posterior as new data come in, and it’s expensive to re-fit the model with all your data each time a new little batch of data come in. We’ve solved this problem in different ways at different times: sometimes we approximate the posterior distribution by a normal or mixture of normals; other times we approximate the posterior for the hyperparameters and use it as a prior for the model we’re fitting to the new data; sometimes the problem can be simplified by factorization in some way or another. There’s also particle filtering, which is a particular class of iterative simulation algorithm using multiple chains.

Naming the problem

As with any problem that comes up often, there are many solutions. I’m excited that Bob recently implemented a version in Stan. Bob called it “chaining Bayesian inference” . . . I’m not thrilled with that term. “Sequential Bayesian updating” sounds better to me, so that’s what I used in the title of this post. I’m open to other suggestions.

Gaurav’s approach

Here’s what Gaurav writes in his new post:

The Problem

Bayesian inference becomes challenging when you need to chain analyses across sequential datasets but lack conjugate priors with closed-form posteriors. While textbook solutions work beautifully for simple cases like binomial-beta models, real-world scenarios often involve complex posteriors that resist analytical treatment.

Chenyang Zhong has an elegant solution using kernel density estimation to approximate posteriors as priors for subsequent analyses. Zhong approximates the posterior from the first dataset using a mixture of normals: p(θ | y₁, x₁) ≈ (1/M) Σᵢ Normal(θ | θᵢ⁽¹⁾, h·I)

where θᵢ⁽¹⁾ are posterior draws from the first analysis, and h controls the kernel bandwidth. This creates an empirical prior for analyzing the second dataset.

The approach is clever—it preserves the shape of the posterior while allowing efficient MCMC sampling. Zhong even develops a sophisticated nearest-neighbor method for fast Metropolis sampling.

A Wasserstein Alternative

The Gaussian distribution that minimizes Wasserstein-2 distance to any empirical distribution is simply the one matching the first two moments. This means we can replace Zhong’s kernel density approximation with:

fit_wasserstein_prior(posterior_draws):
    """Wasserstein-optimal Gaussian approximation"""
    mu = np.mean(posterior_draws, axis=0)
    Sigma = np.cov(posterior_draws.T)
    return mu, Sigma

That’s it. No optimization, no bandwidth selection, no nearest-neighbor graphs.

Implementation Comparison

Zhong’s approach requires:

Choosing bandwidth h Implementing mixture model in Stan Complex log-sum-exp calculations

Wasserstein approach:

Automatic parameter selection Simple multivariate normal prior Standard Stan implementation Comparable or faster performance

Here’s the key Stan difference:

stan// Zhong's mixture approach
model {
  y ~ bernoulli_logit(x * beta);
  vector[B] lp;
  for (b in 1:B) {
    lp[b] = normal_lpdf(beta | beta0[b], h);
  }
  target += log_sum_exp(lp) - log(B);
}

// Wasserstein approach  
model {
  y ~ bernoulli_logit(x * beta);
  beta ~ multi_normal(mu_prior, Sigma_prior);
}

When Does This Work?

The Wasserstein approximation is optimal when:

The true posterior is approximately Gaussian You want to preserve the covariance structure Computational simplicity matters

It may underperform when posteriors have:

Strong multimodality Heavy tails Complex dependence structures

Gaurav includes a bunch more Stan and Python code in his post.