Bayesian modeling with PyMC3 and ArviZ






Córdoba 2018

Bayesian Modeling *


* AKA Bayesian statistics, Probabilistic modeling, probabilistic programming, probabilistic machine learning...

A Gaussian Model


Why do we need software for Bayesian modeling?


Bayesian modeling is conceptually simple :-)

  • (Data, Likelihood, Priors) ---Bayes' theorem--- > Posterior Distribution

Bayesian modeling is mathematically difficult :-(

  • Only a small subset of interesting models can be solved analytically --> we must use numerical methods
  • We need tools to define, solve, interpret and critize models

Enter Probabilistic Programming


  • It is all about abstracting the inference procedure and providing tools for transparent, easy and iterative model building
  • Models are written in code and then compiled to get the posterior



The promise: clear separation of modeling and inference. Practitioners should focus on modeling, not computational details

PyMC3: Probabilistic programming with Python


  • Model building
    • A large collections of probability distributions
    • A clear and powerful syntax
    • Integration with the PyData-stack
  • Inference
    • Markov Chain Monte Carlo
    • Variational Inference
  • Computational backend:
    • Theano --> Speed, automatic differenciation, mathematical optimizations, GPU Support
    • PyMC4 will probably use Tensorflow
In [5]:
with pm.Model() as model_g:
    μ = pm.Normal('μ', 0, 10)  # Prior
    σ = pm.HalfNormal('σ', 25)  # Prior
    y = pm.Normal('y', μ, σ, observed=data)  # Likelihood
    trace_g = pm.sample(1000)
az.plot_joint(trace_g, kind='kde', fill_last=False, textsize=16);

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [σ, μ]
Sampling 2 chains: 100%|██████████| 3000/3000 [00:01<00:00, 2684.83draws/s]

Coal Mining Disasters

with pm.Model() as model_am:
    # Priors
    cp = pm.DiscreteUniform('cp', lower=years.min(), upper=years.max(), testval=1900)

    rate_0 = pm.Exponential('rate_0', 1)
    rate_1 = pm.Exponential('rate_1', 1)

    rates = pm.math.switch(pc >= years, rate_0, rate_1)
    # Likelihood
    dissasters = pm.Poisson('dissasters', rates, observed=data)
    # Inference
    trace_am = pm.sample(1000)

Linear regression




with pm.Model() as linear_regression:
    # Priors
    α = pm.Normal('α', mu=0, sd=10)
    β = pm.Normal('β', mu=0, sd=1, shape=K)
    ϵ = pm.HalfNormal('ϵ', 10)
    μ = α + pm.math.dot(X, β) # linear model
    #Likelihood
    y_pred = pm.Normal('y_pred', mu=μ, sd=ϵ, observed=y_obs)
    # Inference
    trace = pm.Sample(1000)

Hierarchical models

with pm.Model() as model_gh:
    # hyper-priors
    μ_μ = pm.Normal('μ_μ', 0, 10)
    σ_μ = pm.HalfNormal('σ_μ', 25)
    # Priors
    μ = pm.Normal('μ', μ_μ, σ_μ, shape=K) 
    σ = pm.HalfNormal('σ', 25)
    # Likelihood
    y = pm.Normal('y', μ, σ, observed=data)
    # Inference
    trace_gh = pm.sample(1000)

and many more...



PyMC3 documentation

Hierarchical Bayesian neural networks





https://twiecki.github.io/blog/2018/08/13/hierarchical_bayesian_neural_network/

Physical and statistical models



ArviZ: Exploratory analysis of Bayesian models



  • Includes numerical and visual functions for posterior analysis, sample diagnostics, model checking, and comparison.
  • The goal is to provide backend-agnostic tools for exploratory analysis of Bayesian models in Python, by first converting inference data into xarray objects.
  • It is currently in heavy development (0.1.0 version released a week ago).

Xarray



  • N-dimensional labeled arrays and datasets in Python

  • xarray.Dataset is an in-memory representation of a netCDF file

In [6]:
inference_data = az.load_arviz_data('centered_eight')
inference_data
Out[6]:
Inference data with groups:
	> posterior
	> sample_stats
	> posterior_predictive
	> prior
	> observed_data
In [7]:
az.plot_forest(inference_data, var_names=['mu', 'theta'], combined=True, figsize=(10, 8));

Model Comparison



  • WAIC: A way to approximate the generalization error. Measures model-fit and penalizes model complexity
  • LOO: A way to approximate leave-one-out cross validation without re-fitting the model



  • Styles

    az.style.use('arviz-darkgrid')

  • Support for PyMC3 and PyStan objects

  • Python >= 3.4
  • A unified datastructure for probabilistic programming?
  • An opportunity to think deeper about Exploratory analysis of Bayesian models?

PyMC and Arviz future




  • Integrate Approximate Bayesian Computation (part of a GSoC)
  • Discontinuous Hamiltonian Monte Carlo
  • ODE solvers
  • BARTS (Bayesian additive regression trees)
  • Integrate ArviZ into PyMC3
  • More than one plotting-backend for ArviZ (altair is a good candidate)
  • Reach out to other libraries that may be interested to use ArviZ









Image Markov Chain Monte Carlo https://github.com/ColCarroll/imcmc