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.