# Using Multinomial with an underlying model for the probabilities

**URL:** <https://discourse.pymc.io/t/using-multinomial-with-an-underlying-model-for-the-probabilities/3986>\
**Category:** Questions\
**Created:** [October 18, 2019, 3:22pm UTC](https://discourse.pymc.io/t/using-multinomial-with-an-underlying-model-for-the-probabilities/3986 "2019-10-18T15:22:38Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Elijah\_Roberts](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/elijah_roberts/32/2260_2.png) [@Elijah\_Roberts](https://discourse.pymc.io/u/Elijah_Roberts)\
**Post date:** [October 18, 2019, 3:22pm UTC](https://discourse.pymc.io/t/using-multinomial-with-an-underlying-model-for-the-probabilities/3986/1 "2019-10-18T15:22:39Z")

</div>

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!

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

```

---

<div class="post-metadata">

**Author:** ![nze22](https://avatars.discourse-cdn.com/v4/letter/n/e9bcb4/32.png) [@nze22](https://discourse.pymc.io/u/nze22)\
**Post date:** [October 18, 2019, 5:18pm UTC](https://discourse.pymc.io/t/using-multinomial-with-an-underlying-model-for-the-probabilities/3986/2 "2019-10-18T17:18:42Z")

</div>

How about the following:

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

```

---

<div class="post-metadata">

**Author:** ![AlexAndorra](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/alexandorra/32/9142_2.png) [@AlexAndorra](https://discourse.pymc.io/u/AlexAndorra)\
**Post date:** [October 18, 2019, 6:46pm UTC](https://discourse.pymc.io/t/using-multinomial-with-an-underlying-model-for-the-probabilities/3986/3 "2019-10-18T18:46:56Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![clausherther](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/clausherther/32/414_2.png) [@clausherther](https://discourse.pymc.io/u/clausherther)\
**Post date:** [October 23, 2019, 1:00pm UTC](https://discourse.pymc.io/t/using-multinomial-with-an-underlying-model-for-the-probabilities/3986/4 "2019-10-23T13:00:18Z")

</div>

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:

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

```

?

---

<div class="post-metadata">

**Author:** ![Elijah\_Roberts](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/elijah_roberts/32/2260_2.png) [@Elijah\_Roberts](https://discourse.pymc.io/u/Elijah_Roberts)\
**Post date:** [October 23, 2019, 1:51pm UTC](https://discourse.pymc.io/t/using-multinomial-with-an-underlying-model-for-the-probabilities/3986/5 "2019-10-23T13:51:10Z")

</div>

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!

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

```

---

<div class="post-metadata">

**Author:** ![Elijah\_Roberts](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/elijah_roberts/32/2260_2.png) [@Elijah\_Roberts](https://discourse.pymc.io/u/Elijah_Roberts)\
**Post date:** [October 23, 2019, 1:52pm UTC](https://discourse.pymc.io/t/using-multinomial-with-an-underlying-model-for-the-probabilities/3986/6 "2019-10-23T13:52:56Z")

</div>

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