Mixture modeling example in pymc

finite_mixture.py

import numpy as np
from pymc import Model, Gamma, Normal, Dirichlet
from pymc import Categorical
from pymc import sample, Metropolis

k = 3
ndata = 500

v = np.random.randint(0, k, ndata)
data = ((v == 0)*(50 + np.random.randn(ndata))
        + (v == 1)*(-50 + np.random.randn(ndata))
        + (v == 2)*np.random.randn(ndata))

with Model() as model:
    dd = Dirichlet('dd', a=np.array([1., 1., 1.]), shape=k)
    precs = Gamma('precs', alpha=0.1, beta=0.1, shape=k)
    means = Normal('means', 0, 0.001, shape=k)
    category = Categorical('category',
                           p=dd)

    points = Normal('obs',
                     means[category],
                     precs[category],
                     shape=ndata,
                     observed=data)
    tr = sample(3000, step=Metropolis(model.vars))