Dirichlet vs Binomial in pymc3

Question

I am having trouble sampling from a Dirichlet/Multinomial distribution with pymc3. I tried to create a simple test-case to recreate a Beta/Binomial using Dirichlet/Multinomial with n=2, but I can't get it to work.

Below I have some code that works for Binomial but fails for Multinomial. One of the obvious differences is that the Multinomial model is more constrained: i.e. to start, rating is set to 10 in the Binomial model, and [10,10] in the Multinomial. The pymc3 Dirichlet code does say "Only the first k-1 elements of x are expected" but only arrays of shape 2 seem to work in my code.

The output shows that num_friends and rating are being sampled in the Binomial case, but not in the Multinomial case. friends_ratings is being sampled in both. Thanks!

Oh, also Dirichlet('d', np.array([1,1])) crashes with "Floating point error 8". It only appears to fail when two integers of value 1 are passed in. np.array([1.,1.]) works.

import pymc as pm
import numpy as np

print "TEST BINOMIAL"
with pm.Model() as model:
    friends_ratings = pm.Beta('friends_ratings', alpha=1, beta=2)
    num_friends = pm.DiscreteUniform('num_friends', lower=0, upper=100)
    rating = pm.Binomial('rating', n=num_friends, p=friends_ratings)

    step = pm.Metropolis([num_friends, friends_ratings, rating])
    start = {"friends_ratings":.5, "num_friends":20, 'rating':10}

    tr = pm.sample(5, step, start=start, progressbar=False)    
    print "friends", [tr[i]['num_friends'] for i in range(len(tr))]
    print "friends_ratings", [tr[i]['friends_ratings'] for i in range(len(tr))]
    print "rating", [tr[i]['rating'] for i in range(len(tr))]

print "TEST DIRICHLET"
with pm.Model() as model:
    friends_ratings = pm.Dirichlet('friends_ratings', np.array([1.,1.]), shape=2)
    num_friends = pm.DiscreteUniform('num_friends', lower=0, upper=100)    
    rating = pm.Multinomial('rating', n=num_friends, p=friends_ratings, shape=2)

    step = pm.Metropolis([num_friends, friends_ratings, rating])
    start = {'friends_ratings': np.array([0.5,0.5]), 'num_friends': 20, 'rating': [10,10]}

    tr = pm.sample(5, step, start=start, progressbar=False)    
    print "friends", [tr[i]['num_friends'] for i in range(len(tr))]
    print "friends_ratings", [tr[i]['friends_ratings'] for i in range(len(tr))]
    print "rating", [tr[i]['rating'] for i in range(len(tr))]

Output:

TEST BINOMIAL
friends [22.0, 24.0, 24.0, 23.0, 23.0]
friends_ratings [0.5, 0.5, 0.41, 0.41, 0.41]
ratingf [10.0, 11.0, 11.0, 11.0, 11.0]
TEST DIRICHLET
friends [20.0, 20.0, 20.0, 20.0, 20.0]
friends_ratings [array([ 0.51369621,  1.490608  ]), ... ]
rating [array([ 10.,  10.]), array([ 10.,  10.]), ... ]

Answer 1

PyMC3 does not automatically normalize the Dirichlet. So far you have to do this explicitly using simplextransform. See here for an example.

There is an issue of making this transform automatic though: https://github.com/pymc-devs/pymc3/issues/315

EDIT (9/14/2015): PyMC3 now automatically transforms the dirichlet distribution (as any other distribution). So you don't need to specify that manually anymore.