Problem with pm.Categorical

I’m getting strange behavior with a very simple model with a categorical likelihood. I expect that the issue is my lack of experience with this particular kind of model.

Here a minimal example that illustrates the problem:

import numpy as np
import pymc3 as pm
import theano.tensor as tt
import matplotlib.pyplot as plt


with pm.Model():

    y = np.random.choice(4, 1000, p=[0.1, 0.2, 0.6, 0.1])
    print(y)

    probs = []

    for k in range(4):

        _p = pm.Beta(name='p%i' % k, alpha=1, beta=1)
        probs.append(_p)

    p = tt.stack(probs)
    pm.Categorical(name='y', p=p, observed=y)

    trace = pm.sample(draws=10000, tune=2000)
    pm.traceplot(trace)
    plt.savefig('traceplot.png')

I think this model should converge on the correct values of p, but it doesn’t:

Any ideas what’s going wrong here?

pm.Categorical normalized the input vector p so that the p.sum(-1) == 1

This is more transparent if you do:

preal = [0.1, 0.2, 0.6, 0.1]
y = np.random.choice(4, 1000, p=preal)
print(y)
with pm.Model():
    probs = []
    for k in range(4):
        _p = pm.Beta(name='p%i' % k, alpha=1, beta=1)
        probs.append(_p)

    p = tt.stack(probs)
    p1 = pm.Deterministic('p', p/p.sum())
    pm.Categorical(name='y', p=p1, observed=y)

    trace = pm.sample(draws=10000, tune=2000)

pm.traceplot(trace, varnames=['p'], lines=dict(p=preal))

Thanks! Looks like I was being completely dumb!

lol. BTW you might want to use a Dirichlet prior so it generates unit vector p

For future reference, here is the codification of junpenglao Dirichlet prior.

preal = [0.1, 0.2, 0.6, 0.1]
y = np.random.choice(4, 1000, p=preal)

with pm.Model():
    
    p = pm.Dirichlet('p', a=np.ones(4))

    pm.Categorical(name='y', p=p, observed=y)

    trace = pm.sample(draws=1000, tune=200)

pm.traceplot(trace, varnames=['p'], lines=dict(p=preal))
1 Like