mattjj/pybasicbayes

This library provides objects that model probability distributions and the related operations that are common in generative Bayesian modeling and Bayesian inference, including Gibbs sampling and variational mean field algorithms. The file abstractions.py describes the queries a distribution must support to be used in each algorithm, as well as an API for models, which compose the distribution objects.

Example

The file models.py shows how to construct mixture models building on the distribution objects in this library. For example, to generate data from a Gaussian mixture model, we might set some hyperparameters, construct a Mixture object, and then ask it to randomly generate some data from the prior:

import numpy as np
from pybasicbayes import models, distributions

# hyperparameters
alpha_0=5.0
obs_hypparams = dict(mu_0=np.zeros(2),sigma_0=np.eye(2),kappa_0=0.05,nu_0=5)

# create the model
priormodel = models.Mixture(alpha_0=alpha_0,
        components=[distributions.Gaussian(**obs_hypparams) for itr in range(30)])

# generate some data
data = priormodel.rvs(400)

# delete the model
del priormodel

If we throw away the prior model at the end, we're left just with the data, which look like this:

from matplotlib import pyplot as plt
plt.figure()
plt.plot(data[:,0],data[:,1],'kx')
plt.title('data')

Imagine we loaded these data from some measurements file and we wanted to fit a mixture model to it. We can create a new Mixture and run inference to get a representation of the posterior distribution over mixture models conditioned on observing these data:

posteriormodel = models.Mixture(alpha_0=alpha_0,
        components=[distributions.Gaussian(**obs_hypparams) for itr in range(30)])

posteriormodel.add_data(data)

Since pybasicbayes implements both Gibbs sampling and variational mean field inference algorithms, we can use both together in a hybrid algorithm.

import copy
from pybasicbayes.util.text import progprint_xrange

allscores = [] # variational lower bounds on the marginal data log likelihood
allmodels = []
for superitr in range(5):
    # Gibbs sampling to wander around the posterior
    print 'Gibbs Sampling'
    for itr in progprint_xrange(100):
        posteriormodel.resample_model()

    # mean field to lock onto a mode
    print 'Mean Field'
    scores = [posteriormodel.meanfield_coordinate_descent_step()
                for itr in progprint_xrange(100)]

    allscores.append(scores)
    allmodels.append(copy.deepcopy(posteriormodel))

import operator
models_and_scores = sorted([(m,s[-1]) for m,s
    in zip(allmodels,allscores)],key=operator.itemgetter(1),reverse=True)

Now we can plot the score trajectories:

plt.figure()
for scores in allscores:
    plt.plot(scores)
plt.title('model vlb scores vs iteration')

And show the point estimate of the best model by calling the convenient Mixture.plot():

models_and_scores[0][0].plot()
plt.title('best model')

Since these are Bayesian methods, we have much more than just a point estimate for plotting: we have fit entire distributions, so we can query any confidence or marginal that we need.

See the file demo.py for the code for this demo.

Authors

Matt Johnson, Alex Wiltschko, Yarden Katz, Nick Foti, and Scott Linderman.