Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Christopher Fonnesbeck - Probabilistic Programm...

Christopher Fonnesbeck - Probabilistic Programming with PyMC3

Bayesian statistics offers robust and flexible methods for data analysis that, because they are based on probability models, have the added benefit of being readily interpretable by non-statisticians. Until recently, however, the implementation of Bayesian models has been prohibitively complex for use by most analysts. But, the advent of probabilistic programming has served to abstract the complexity of Bayesian statistics, making such methods more broadly available. PyMC3 is a open-source Python module for probabilistic programming that implements several modern, computationally-intensive statistical algorithms for fitting Bayesian models, including Hamiltonian Monte Carlo (HMC) and variational inference. PyMC3’s intuitive syntax is helpful for new users, and the reliance on Theano for much of the computational work has allowed developers to keep the code base simple, making it easy to extend the software to meet analytic needs. PyMC3 itself extends Python's powerful "scientific stack" of development tools, which provide fast and efficient data structures, parallel processing, and interfaces for describing statistical models.

https://us.pycon.org/2017/schedule/presentation/473/

PyCon 2017

May 21, 2017
Tweet

More Decks by PyCon 2017

Other Decks in Programming

Transcript

  1. Stochastic language "primitives" Distribution over values: X ~ Normal(μ, σ)

    x = X.random(n=100) Distribution over functions: Y ~ GaussianProcess(mean_func(x), cov_func(x)) y = Y.predict(x2) Conditioning: p ~ Beta(1, 1) z ~ Bernoulli(p) # z|p
  2. Why Bayes? ❝The Bayesian approach is attractive because it is

    useful. Its usefulness derives in large measure from its simplicity. Its simplicity allows the investigation of far more complex models than can be handled by the tools in the classical toolbox.❞ —Link and Barker 2010
  3. model { for (j in 1:J){ y[j] ~ dnorm (theta[j],

    tau.y[j]) theta[j] ~ dnorm (mu.theta, tau.theta) tau.y[j] <- pow(sigma.y[j], -2) } mu.theta ~ dnorm (0.0, 1.0E-6) tau.theta <- pow(sigma.theta, -2) sigma.theta ~ dunif (0, 1000) }
  4. PyMC3 ☞ started in 2003 ☞ PP framework for fitting

    arbitrary probability models ☞ based on Theano ☞ implements "next generation" Bayesian inference methods ☞ NumFOCUS sponsored project !" github.com/pymc-devs/pymc3
  5. Calculating Gradients in Theano >>> from theano import function, tensor

    as tt >>> x = tt.dmatrix('x') >>> s = tt.sum(1 / (1 + tt.exp(-x))) >>> gs = tt.grad(s, x) >>> dlogistic = function([x], gs) >>> dlogistic([[3, -1],[0, 2]]) array([[ 0.04517666, 0.19661193], [ 0.25 , 0.10499359]])
  6. Calculating Gradients in Theano >>> from theano import function, tensor

    as tt >>> x = tt.dmatrix('x') >>> s = tt.sum(1 / (1 + tt.exp(-x))) >>> gs = tt.grad(s, x) >>> dlogistic = function([x], gs) >>> dlogistic([[3, -1],[0, 2]]) array([[ 0.04517666, 0.19661193], [ 0.25 , 0.10499359]])
  7. Priors with Model() as unpooled_model: α = Normal('α', 0, sd=1e5,

    shape=counties) β = Normal('β', 0, sd=1e5) σ = HalfCauchy('σ', 5)
  8. Bayesian approximation ☞ Maximum a posteriori (MAP) estimate ☞ Laplace

    (normal) approximation ☞ Rejection sampling ☞ Importance sampling ☞ Sampling importance resampling (SIR) ☞ Approximate Bayesian Computing (ABC)
  9. MCMC Markov chain Monte Carlo simulates a Markov chain for

    which some function of interest is the unique, invariant, limiting distribution.
  10. MCMC Markov chain Monte Carlo simulates a Markov chain for

    which some function of interest is the unique, invariant, limiting distribution. This is guaranteed when the Markov chain is constructed that satisfies the detailed balance equation:
  11. Hamiltonian Monte Carlo Uses a physical analogy of a frictionless

    particle moving on a hyper-surface Requires an auxiliary variable to be specified ☞ position (unknown variable value) ☞ momentum (auxiliary)
  12. Hamiltonian MC ➀ Sample a new velocity from univariate Gaussian

    ➁ Perform n leapfrog steps to obtain new state ➂ Perform accept/reject move of
  13. Variational Inference Variational inference minimizes the Kullback-Leibler divergence from approximate

    distributions, but we can't calculate the true posterior distribution.
  14. ADVI* * Kucukelbir, A., Tran, D., Ranganath, R., Gelman, A.,

    & Blei, D. M. (2016, March 2). Automatic Differentiation Variational Inference. arXiv.org.
  15. with partial_pooling: approx = fit(n=100000) Average Loss = 1,115.5: 100%|██████████|

    100000/100000 [00:13<00:00, 7690.51it/s] Finished [100%]: Average Loss = 1,115.5
  16. with partial_pooling: approx = fit(n=100000) Average Loss = 1,115.5: 100%|██████████|

    100000/100000 [00:13<00:00, 7690.51it/s] Finished [100%]: Average Loss = 1,115.5 >>> approx <pymc3.variational.approximations.MeanField at 0x119aa7c18>
  17. New PyMC3 Features ☞ Gaussian processes ☞ Elliptical slice sampling

    ☞ Sequential Monte Carlo methods ☞ Time series models
  18. Convolutional variational autoencoder♌ import keras class Decoder: ... def decode(self,

    zs): keras.backend.theano_backend._LEARNING_PHASE.set_value(np.uint8(0)) return self._get_dec_func()(zs) ♌ Taku Yoshioka (c) 2016
  19. with pm.Model() as model: # Hidden variables zs = pm.Normal('zs',

    mu=0, sd=1, shape=(minibatch_size, dim_hidden), dtype='float32') # Decoder and its parameters dec = Decoder(zs, net=cnn_dec) # Observation model xs_ = pm.Normal('xs_', mu=dec.out.ravel(), sd=0.1, observed=xs_t.ravel(), dtype='float32')
  20. The Future ☞ Operator Variational Inference ☞ Normalizing Flows ☞

    Riemannian Manifold HMC ☞ Stein variational gradient descent ☞ ODE solvers
  21. The PyMC3 Team ☞ Colin Carroll ☞ Peadar Coyle ☞

    Bill Engels ☞ Maxim Kochurov ☞ Junpeng Lao ☞ Osvaldo Martin ☞ Kyle Meyer ☞ Austin Rochford ☞ John Salvatier ☞ Adrian Seyboldt ☞ Thomas Wiecki ☞ Taku Yoshioka