Quadratic Approximation

Modern Bayesian data analysis is based on efficient Markov Chain Monte Carlo (MCMC) techniques.

Learning Bayesian statistic can be hard for a Frequentist! The Statistical Rethinking book by Richard McElreath tries to avoid cognitive overload for its readers by not requireing them to learn Bayesian statistics and MCMC at the same time. The posterior distribution is in most interesting cases analytically intractable and hence MCMC is used to numerically determine it.

As a simpler alternative to MCMC one can you Quadratic Approximation¹. A lot of people has ported the R code examples from Statistical Rethinking to Python using frameworks like PYMC, Pyro, NumPyro, and TensorFlow Probability. Apparently there is no Quadratic Approximation solution available for PYMC². Numpyro has AutoLaplaceApproximation³ but it is not clear to me that this is the same as Quadratic Approximation.

quad5 leverages PYMC and works by adding a custom step for the sample method on PYMC models. By doing so it benefits from standard PYMC functionality that automatically adds constant, deterministic, and observed data to the inferencedata output from the sample method. This allows for easy usage of the Arviz visualization library for posterior distributions and more general the ecosystem around PYMC.

Word of warning

This package is not production-grade code. It is primarily meant for educational purposes. Secondary for the author to learn more about the internals of the PYMC library.

I am quite sure there will be valid PYMC models where this package not is able to produce a quadratic approximation for the posterior distribution. You are more than welcome to submit either a PR or create a issue for such cases.

Example

        import arviz as az
        import numpy as np
        import pymc as pm
        import quad5 as quad5

        y = np.array([2642, 3503, 4358], dtype=np.float64)

        with pm.Model() as m:
            logsigma = pm.Uniform("logsigma", 1, 100)
            mu = pm.Uniform("mu", 0, 10000)
            _ = pm.Normal("y", mu=mu, sigma=pm.math.exp(logsigma), observed=y)
            custom_step = quad5.QuadraticApproximation([mu, logsigma], m)
            trace = pm.sample(draws=1000, chains=4, tune=100, step=custom_step)

        az.plot_posterior(trace)    

See more examples in this notebook with examples from chapter 4 in Statistical Rethinking.