Probabilistic Programming Languages

Building a custom Bayesian model the traditional way means writing down a prior, a likelihood, deriving the posterior's exact or approximate form by hand, and then hand-coding a sampler (Gibbs, Metropolis-Hastings, or otherwise) specific to that exact model, real work, redone from scratch for every new model. A probabilistic programming language (PPL), Stan, PyMC, NumPyro, removes almost all of that: you write the generative story as a short program (this parameter is drawn from a Normal, that observation is drawn from a Poisson whose rate depends on the parameter), and the language's built-in engine automatically produces posterior samples, no custom derivation or sampler-coding required.

The analogy: describing a recipe versus writing the kitchen's operating manual

Building a custom sampler by hand is like being asked, every time you want a new dish, to also write out the complete operating manual for a fully automated kitchen capable of cooking exactly that dish, an enormous amount of engineering for one specific outcome. A probabilistic programming language is like handing a general-purpose chef (the inference engine) a recipe written in plain steps, "start with 2 cups flour drawn from this distribution, then combine with...", and trusting the chef to figure out how to execute it correctly, no matter what specific recipe you hand over next. You describe what the model is; the engine figures out how to fit it.

What the code actually specifies

A PPL program specifies a generative model: a sequence of random variable declarations, each with a distribution, often depending on variables declared earlier. Schematically, in Stan-like pseudocode:

mu ~ Normal(0, 10)          # prior on the mean
sigma ~ HalfNormal(5)       # prior on the standard deviation
y ~ Normal(mu, sigma)       # likelihood: observed data y

The engine reads this declaration, automatically constructs the joint density $p(\mu, \sigma, y) = p(\mu)\,p(\sigma)\,p(y\mid \mu,\sigma)$ implied by the code, exactly the Bayesian network factorization idea, just written as executable code instead of drawn as a graph, and then runs a general-purpose sampler (typically Hamiltonian Monte Carlo) against that joint density to produce samples from the posterior $p(\mu,\sigma\mid y)$, with no model-specific math derived by the user at all.

Worked example 1: fitting a simple mean-estimation model

Suppose you observe 5 data points, $y = [5.1, 4.8, 5.3, 5.0, 4.9]$, and want the posterior for the true mean $\mu$, assuming known standard deviation $\sigma=0.3$ and prior $\mu \sim \text{Normal}(0, 10)$. The full derivation by hand: with a conjugate Normal-Normal setup, the exact posterior is $\text{Normal}(\mu_{\text{post}}, \sigma_{\text{post}}^2)$ where $\sigma_{\text{post}}^2 = \left(\frac{1}{10^2} + \frac{5}{0.3^2}\right)^{-1} \approx \left(0.01 + 55.6\right)^{-1} \approx 0.018$, and $\mu_{\text{post}} = \sigma_{\text{post}}^2\left(\frac{0}{10^2} + \frac{5\times 5.02}{0.3^2}\right) \approx 0.018 \times 278.9 \approx 5.02$ (using $\bar y = 5.02$). A PPL user never derives this: they write the three-line model above, call the sampler, and get back thousands of posterior samples for $\mu$ whose mean and spread numerically match $5.02$ and $\sqrt{0.018}\approx 0.135$ almost exactly, validating that the automatic engine reproduces what the hand-derivation would have given, without the user ever doing the algebra.

Worked example 2: the payoff shows up on a harder model

Now add a twist: $\sigma$ is also unknown, with its own prior, and the model has a hierarchical structure, say, 5 different desks each with their own mean $\mu_k$, all drawn from a shared "meta-mean." No clean, closed-form posterior exists for this model at all; deriving a custom Gibbs sampler by hand would require working out several conditional distributions analytically, a multi-page derivation prone to errors. In a PPL, the change from the simple model is just a few more lines declaring the hierarchical structure, the sampler code itself does not change at all, because it was never written for a specific model in the first place; it works against the joint density mechanically, regardless of how many layers or how much structure the model has.

The prior-to-posterior update shown above is exactly the computation a PPL automates at scale, instead of one Beta-Binomial pair worked by hand, a general engine runs the same logic for models with dozens of interacting parameters.

What this means in practice

PPLs are the practical tool of choice for real-world Bayesian modeling in finance: hierarchical models pooling information across many similar but distinct groups (desks, sectors, strategies), custom risk models with non-standard likelihoods that have no closed-form conjugate solution, and any model complex enough that hand-deriving a sampler would be impractical. The tradeoff is computational cost and diagnostic burden: general-purpose samplers like Hamiltonian Monte Carlo can still be slow on very large datasets or badly-specified models, and every fit still needs convergence diagnostics to confirm the sampler actually explored the posterior properly, the engine automates the mechanics, not the judgment.

Practice

1. Redo the mean-estimation posterior calculation in worked example 1 with $\sigma = 1.0$ instead of $0.3$, and note how much wider the resulting posterior becomes.
2. Explain in one sentence why a PPL's sampler code does not need to change when the model becomes hierarchical.
3. Name one situation where deriving a model by hand (rather than using a PPL) would still be the better choice.