Using Multinomial with an underlying model for the probabilities

Questions
1

Hi all,
I’m trying to formulate a simple example in PyMC3 of using a Multinomial distribution to infer the probabilities of a die throw, except where the side probabilities are linked by a specific model with some lesser number of parameters. In the example below it is just a linearly decreasing probability with some characteristic slope m, which is the parameter to infer. I’ve got it working using the code below by setting up a Uniform random variable for the prior and then making a Deterministic variable for each side probability. But I don’t really like having all the terms written out individually, especially if one would extend the number of possible states it would get very ugly. I was thinking of using a custom DensityDist to specify the probabilities, but I can’t quite figure out if that is feasible or not. Or possibly I could write a custom distribution, but also it is not readily apparent how to do that. Does anyone have any insight on how to improve the code below? Thanks!

import numpy as np
import pymc3 as mc

def randomDataLinearDiceModel(p, n):
    
    # Calculate the cumulative probabilities.
    c = np.cumsum(p)

    # Generate a random sequence of dice rolls.
    r = np.random.rand(n)
    x = np.zeros((n,), dtype=np.uint8)
    for i in range(n):
        x[i] = np.argmax(r[i]<c)

    # Get the counts for each die side.
    (counts,edges) = np.histogram(x,np.arange(0-0.5,len(p)+0.5,1.0))
    centers = (edges[:-1]+edges[1:])/2
    
    return centers,counts

def bayesianInferenceLinearDiceModel(p, n, counts):
    
    with mc.Model() as model:

        # Prior distributions.
        m = mc.Uniform('m', lower=0, upper=0.05)
        p0 = mc.Deterministic('p0', (1/6) + 3*m)
        p1 = mc.Deterministic('p1', (1/6) + 2*m)
        p2 = mc.Deterministic('p2', (1/6) + 1*m)
        p3 = mc.Deterministic('p3', (1/6) - 1*m)
        p4 = mc.Deterministic('p4', (1/6) - 2*m)
        p5 = mc.Deterministic('p5', (1/6) - 3*m)

        # Likelihood.
        counts_likelihood = mc.Multinomial('counts', n=n, p=[p0,p1,p2,p3,p4,p5], observed=counts)

        # Perform the sampling.
        trace = mc.sample(2000, tune=1000)
        
        return trace


m=0.045
p = np.array([1/6+3*m,1/6+2*m,1/6+1*m,1/6-1*m,1/6-2*m,1/6-3*m])
n=10000
sides,rolls = randomDataLinearDiceModel(p, n)
trace = bayesianInferenceLinearDiceModel(p, n, rolls)
print("Posterior distribution summary:")
print(mc.summary(trace))
print("Finished.")
2

How about the following:

def bayesianInferenceLinearDiceModel(p, n, counts):

    lst = np.concatenate((np.arange(3,0,-1), np.arange(-1,-4,-1)))
    with mc.Model() as model:

        # Prior distributions.
        m = mc.Uniform('m', lower=0, upper=0.05)
        p = mc.Deterministic('p', 1/6 + lst*m)

        # Likelihood.
        counts_likelihood = mc.Multinomial('counts', n=n, p=p, observed=counts)

        # Perform the sampling.
        trace = mc.sample(2000, tune=1000)
        
        return trace


m=0.045
p = np.array([1/6+3*m,1/6+2*m,1/6+1*m,1/6-1*m,1/6-2*m,1/6-3*m])
n=10000
sides,rolls = randomDataLinearDiceModel(p, n)
trace = bayesianInferenceLinearDiceModel(p, n, rolls)
print("Posterior distribution summary:")
print(mc.summary(trace))
print("Finished.")
1 Like
3

@nze22 solution looks good to me!
You might also wanna look into softmax regression if you need something more generalizable - e.g. you have to estimate the intercepts also, or the intercept and/or slopes vary by category, etc.

4

Do we need to pass in the previously calculated values for m and p to the model, since they’re RVs in the model?
So, for clarity maybe:

def bayesianInferenceLinearDiceModel(n, counts):

    lst = np.concatenate((np.arange(3,0,-1), np.arange(-1,-4,-1)))
    with mc.Model() as model:

        # Prior distributions.
        m = mc.Uniform('m', lower=0, upper=0.05)
        p = mc.Deterministic('p', 1/6 + lst*m)

        # Likelihood.
        counts_likelihood = mc.Multinomial('counts', n=n, p=p, observed=counts)

        # Perform the sampling.
        trace = mc.sample(2000, tune=1000)
        
        return trace


m=0.045
p = np.array([1/6+3*m,1/6+2*m,1/6+1*m,1/6-1*m,1/6-2*m,1/6-3*m])
n=10000
_, counts = randomDataLinearDiceModel(p, n)
trace = bayesianInferenceLinearDiceModel(n, counts)
print("Posterior distribution summary:")
print(mc.summary(trace))
print("Finished.")

?

5

First, thanks for the reply. I didn’t realize that the deterministic RV would automatically be expanded to the correct dimensions. That said, I now have the two following versions, one of which uses Deterministic variables and one which just does Theano math directly, which both work correctly. But I don’t have a good feel for if one approach is better for performance in PyMC3. As in if I were using a 1000 side dice would making 1000 deterministic RVs be bad compared to manipulating the probabilities using Theano. And in that case would it be better to make a custom distribution class, similar to how the Dirichlet can be used to specify the p values for Multinomial, or does it not matter? Thanks for any insight anyone might have!

def linearDiceProbabilities_1(m):
    return 1/6 + np.array([3,2,1,-1,-2,-3])*m

def linearDiceProbabilities_2(m):
    return mc.Deterministic('p', 1/6 + np.array([3,2,1,-1,-2,-3])*m)

def bayesianInferenceLinearDiceModel(p, n, counts):

    with mc.Model() as model:

        # Prior distributions.
        m = mc.Uniform('m', lower=0, upper=0.05)
        p = linearDiceProbabilities_1(m)
        #p = linearDiceProbabilities_2(m)

        # Likelihood.
        counts_likelihood = mc.Multinomial('counts', n=n, p=p, observed=counts)

        # Perform the sampling.
        trace = mc.sample(2000, tune=1000)
        
        return trace
6

@clausherther Yes, that’s a good point. That was just leftover code.

1 Like