Probabilistic Programming

Across the earlier lessons in this domain, building a probabilistic model and computing with it were treated as two separate acts: specify a Bayesian network, hidden Markov model, or state-space model's structure, and then apply (or hand-derive) an inference algorithm suited to that specific structure - variable elimination for one kind of network, the forward algorithm for another, a bespoke Metropolis-Hastings sampler for a third. Probabilistic programming proposes to collapse that distinction. A probabilistic program is an ordinary-looking program extended with two additional primitives: a sample statement, which draws a random value from a specified distribution, and an observe (or condition) statement, which asserts that a particular random variable's value matches an observed data point. Together, these two primitives let a program describe an entire generative process - however it branches, loops, or recurses - implicitly defining a probability distribution over everything the program could have done, and a generic inference engine, rather than a hand-built algorithm specific to that one model, computes the posterior distribution over any of the program's random choices given the data asserted by its observe statements.

Goodman, Mansinghka, Roy, Bonawitz, and Tenenbaum's Church, an early and influential such language built on Lisp, made this idea concrete by treating ordinary program evaluation itself as sampling from the underlying generative process: evaluating a Church program once produces one sample from the process it describes, and a query operation generalizes plain evaluation to conditional sampling - drawing from the distribution over some expression's value given that some other expression (the evidence) is constrained to be true. Church also introduced stochastic memoization, letting the same random choice be reused consistently across multiple uses within a single execution, which matters for expressing nonparametric models where the number of underlying objects or parameters is itself unknown and can grow as the program runs, something that is awkward to express in fixed-dimensional graphical model notation, where the number of variables and their relationships must typically be fixed in advance before inference can even begin. Church's inference engine handled queries through generic Monte Carlo methods, including rejection sampling and Metropolis-Hastings, applied automatically to whatever program the user had written, without requiring a bespoke algorithm tailored to that program's particular structure.

The conceptual shift this represents is significant: a modeler writes down what the generative process is, in terms natural to the domain, rather than working out how to compute its posterior by hand for each new modeling choice. This genuinely repays the graphical-model investment from earlier in this domain, because most probabilistic programming inference engines, under the hood, still construct something resembling a graph or factor graph over the program's execution trace - the dependency structure a Bayesian network or factor graph would encode explicitly is implicit in which variables a program's control flow actually touches on a given run, and the inference engine has to reconstruct or approximate that structure automatically, on the fly, often differently across different executions of the same program when the program's control flow itself depends on random choices it makes along the way. A model expressed with loops or recursion can therefore define an effectively unbounded graphical model, something ordinary fixed-structure Bayesian networks were never designed to represent directly.

Stan, described by Carpenter and a large group of collaborators, takes a somewhat different approach from Church's sampling-based program semantics: a Stan program imperatively defines a log-probability (target) function over its parameters, conditioned on the data supplied to it, rather than describing a step-by-step generative sampling process directly. Stan performs full Bayesian inference over continuous parameters primarily via Hamiltonian Monte Carlo, and specifically its No-U-Turn Sampler (NUTS), an adaptive variant of HMC that automatically tunes the step size and trajectory length that plain HMC would otherwise require a user to hand-tune. Hamiltonian Monte Carlo's core advantage over a plain random-walk Metropolis sampler is that it uses the gradient of the log-target density to simulate a physically-inspired trajectory through parameter space, proposing distant moves with high acceptance probability rather than the small, undirected steps a random-walk proposal is forced to take - a difference that becomes decisive as the number of parameters grows, since random-walk proposals scale very poorly with dimension while gradient-guided proposals scale considerably better.

PyMC3, described by Salvatier, Wiecki, and Fonnesbeck, brings a similar capability into Python, again centered on NUTS as its default sampler, and again relying on automatic differentiation (originally implemented via the Theano library) to obtain the gradients that NUTS needs. This is the detail that makes the whole approach practical for arbitrary user-written models: rather than requiring a user to symbolically differentiate their particular log-posterior by hand - a task that becomes error-prone and tedious for any model of realistic complexity - automatic differentiation computes exact gradients mechanically by applying the chain rule through every operation the program actually performed while evaluating the log-probability, for essentially any model expressible in the language. Automatic differentiation is, in this sense, the general enabling technology underneath both Stan's and PyMC's gradient-based samplers, and its later adoption by deep-learning frameworks is part of what made combining probabilistic programming with large neural models feasible, as the next lesson takes up.

None of this automation removes the modeler's responsibility for judging whether the resulting posterior is actually trustworthy. Posterior predictive checks (simulating data from the fitted model and comparing it to the real observed data), sensitivity analysis to the choice of prior, and convergence diagnostics for the MCMC chain itself - effective sample size and measures like the potential scale reduction factor computed across multiple independent chains started from different initial values - remain necessary practice regardless of how automated the sampling mechanics have become. A model can be syntactically valid, run without errors, and still produce a posterior that is a poor fit to the data or a sampler that never actually converged to that posterior, and no amount of automatic differentiation or adaptive step-size tuning changes that. Probabilistic programming automates the computational mechanics of inference; it does not automate away the modeler's need to check that a model's assumptions are reasonable or that its sampler actually converged before trusting the resulting conclusions.

The combination of probabilistic programming with deep learning motivated a further generation of languages designed for models with many more parameters and neural-network components embedded directly inside the generative process, rather than treating deep learning and probabilistic modeling as two separate toolkits that a practitioner has to stitch together by hand. Pyro, described by Bingham and a large group of collaborators as "deep universal probabilistic programming," is built on PyTorch, giving it access to GPU-accelerated computation and, importantly, to the same automatic differentiation infrastructure that PyTorch uses for training ordinary neural networks - letting a single Pyro model mix conventional neural network layers with explicit sample and observe statements in the same program, so that a neural network can, for instance, define the parameters of a distribution that a later part of the program samples from, or a probabilistic program's structured output can feed into a downstream neural component. This tight integration is the point: a purely neural model has no natural way to represent calibrated uncertainty over its own predictions, while a purely classical probabilistic model has no natural way to absorb the representational power of learned features from raw, high-dimensional data such as images or text - combining the two in one language lets a single model draw on both strengths at once.

At the scale these models reach - potentially millions of parameters, when a probabilistic program's generative process includes a full neural network - plain MCMC, including gradient-based methods like NUTS, generally becomes impractical: even a single HMC trajectory requires computing gradients through the entire model repeatedly, and a Markov chain that must be run for a large number of iterations to adequately explore a high-dimensional space simply does not scale to that regime in reasonable time. This is where the variational inference lesson from earlier in this domain becomes directly practical rather than just theoretically appealing: rather than sampling, Pyro's default approach is stochastic variational inference, which posits a tractable family of approximate posteriors (as Jordan, Ghahramani, Jaakkola, and Saul's synthesis described for graphical models generally) and optimizes that family's parameters via stochastic gradient ascent on a lower bound on the log-evidence, using mini-batches of data the same way ordinary neural network training does. This trades sampling's asymptotic exactness for the same kind of scalable, gradient-based optimization that made deep learning itself practical at scale.

It is worth being honest about where this leaves the field. Exact inference in simple, small-to-moderate graphical models, and basic sampling methods like Metropolis-Hastings and Gibbs sampling, are by now settled, textbook technology with well-understood guarantees. Scaling probabilistic programming to very large, very expressive models - deep generative models, models mixing discrete and continuous latent structure, models that need calibrated, trustworthy uncertainty estimates in a production system rather than just a point prediction - remains active research and engineering territory. Questions like how to reliably diagnose convergence for stochastic variational inference (which lacks MCMC's asymptotic exactness guarantee), how to combine discrete latent variables with gradient-based inference cleanly, and how to make automated inference scale without silently degrading the quality of the resulting uncertainty estimates, do not yet have fully settled, universally reliable answers. Students entering this area should therefore treat the mechanics covered in the first two lessons of this module - the sample/observe abstraction, and gradient-based samplers backed by automatic differentiation - as durable foundations to build on, while treating claims about any particular large-scale system's speed, scalability, or uncertainty calibration as findings to verify against current evidence rather than settled fact, since the state of the art in this corner of the field continues to move.