# Using Multinomial with an underlying model for the probabilities

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.")
``````

``````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

@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.

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.")
``````

?

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
``````

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

1 Like