Tag Archives: pymc3

MCMC in Python: Power-posterior calculation with PyMC3

I’ve got a fun little project that has let me check in on the PyMC project after a long time away. Long-time readers of Healthy Algorithms might remember my obsession with PyMC2 from my DisMod days nearly ten years ago, but for those of you joining us more recently… there is a great way to build Bayesian statistical models with Python, and it is the PyMC package. It was completely rewritten for version 3, which is now available, and using it looks like this:

import numpy as np
import matplotlib.pyplot as plt

# Initialize random number generator

# True parameter values
alpha, sigma = 1, 1
beta = [1, 2.5]

# Size of dataset
size = 100

# Predictor variable
X1 = np.random.randn(size)
X2 = np.random.randn(size) * 0.2

# Simulate outcome variable
Y = alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(size)*sigma

from pymc3 import Model, Normal, HalfNormal
basic_model = Model()

with basic_model:

    # Priors for unknown model parameters
    alpha = Normal('alpha', mu=0, sd=10)
    beta = Normal('beta', mu=0, sd=10, shape=2)
    sigma = HalfNormal('sigma', sd=1)

    # Expected value of outcome
    mu = alpha + beta[0]*X1 + beta[1]*X2

    # Likelihood (sampling distribution) of observations
    Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)


What I need is a way to compute the “power-posterior” of a model, which is to say p(\theta\mid y,t) = \frac{p(y\mid \theta)^t p(\theta)}{\mathcal{Z}_t(y)}  .

I’ve figured out a pretty cool way to compute the prior and the likelihood of the model:

from theano import theano, tensor as tt

def logpriort(self):
    """Theano scalar of log-prior of the model"""
    factors = [var.logpt for var in self.free_RVs]
    return tt.add(*map(tt.sum, factors))
def logprior(self):
    """Compiled log probability density function"""
    return self.model.fn(logpriort(self))
def logliket(self):
    """Theano scalar of log-prior of the model"""
    factors = [var.logpt for var in self.observed_RVs] + self.potentials
    return tt.add(*map(tt.sum, factors))
def loglike(self):
    """Compiled log probability density function"""
    return self.model.fn(logliket(self))

Now I just need to put them together, to get the log of the power-posterior. There must be an easy way to do it, but this is not it:

def logpowerpt(self, t):
    return t*loglike(self) + logprior(self)
def logpowerp(self, t):
    return self.model.fn(logpowerpt(self, t))

vals = {'alpha':1,
logpowerp(basic_model, 0)(vals)

Any suggestions?

Comments Off on MCMC in Python: Power-posterior calculation with PyMC3

Filed under MCMC

Laplace approximation in PyMC3, revisited

I found an even better example of the value of Laplace approximation, and its just a small tweak to the example I did a few weeks ago: http://nbviewer.ipython.org/gist/aflaxman/6d0a9ff2441348f3a130

Comments Off on Laplace approximation in PyMC3, revisited

Filed under statistics

Porting PyMC2 models to PyMC3

It is time for me to start doing this. Here is a StackOverflow question and answer about it, if it is time for you, too: http://stackoverflow.com/questions/30798447/porting-pymc2-code-to-pymc3-hierarchical-model-for-sports-analytics/30853077#30853077

Comments Off on Porting PyMC2 models to PyMC3

Filed under software engineering

MCMC in Python: Gaussian mixture model in PyMC3

PyMC3 is really coming along. I tried it out on a Gaussian mixture model that was the subject of some discussion on GitHub: https://github.com/pymc-devs/pymc3/issues/443#issuecomment-109813012 http://nbviewer.ipython.org/gist/aflaxman/64f22d07256f67396d3a

1 Comment

Filed under MCMC, software engineering, statistics

Laplace approximation in Python: another cool trick with PyMC3

I admit that I’ve been skeptical of the complete rewrite of PyMC that underlies version 3. It seemed to me motivated by an interest in using unproven new step methods that require knowing the derivative of the posterior distribution. But, it is really coming together, and regardless of whether or not the Hamiltonian Monte Carlo stuff pays off, there are some cool tricks you can do when you can get derivatives without a hassle.

Exhibit 1: A Laplace approximation approach to fitting mixed effect models (as described in http://www.seanet.com/~bradbell/tmb.htm)


Comments Off on Laplace approximation in Python: another cool trick with PyMC3

Filed under software engineering, statistics