# Multinomial probabilities don't add to 1

**URL:** <https://discourse.pymc.io/t/multinomial-probabilities-dont-add-to-1/13609>\
**Category:** v5\
**Tags:** modeling\
**Created:** [January 10, 2024, 4:35am UTC](https://discourse.pymc.io/t/multinomial-probabilities-dont-add-to-1/13609 "2024-01-10T04:35:16Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![yahrMason](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/yahrmason/32/5826_2.png) [@yahrMason](https://discourse.pymc.io/u/yahrMason)\
**Post date:** [January 10, 2024, 4:35am UTC](https://discourse.pymc.io/t/multinomial-probabilities-dont-add-to-1/13609/1 "2024-01-10T04:35:16Z")

</div>

Hello,  
I am working on converting [this notebook](https://github.com/pymc-devs/pymc-resources/blob/main/Rethinking_2/Chp_11.ipynb) to PyMC v5 (from v4). I am running into an issue with code block 11.59. This is the code:

```auto
N = 500

# simulate family incomes for each individual
family_income = np.random.rand(N)

# assign a unique coefficient for each type of event
b = np.array([-2.0, 0.0, 2.0])

p = softmax(np.array([0.5, 1.0, 1.5])[:, None] + np.outer(b, family_income), axis=0).T

career = np.asarray([np.random.multinomial(1, pp) for pp in p])
career = np.where(career == 1)[1]

with pm.Model() as m11_14:
    a = pm.Normal("a", 0.0, 1.5, shape=2) # intercepts
    b = pm.Normal("b", 0.0, 1.0, shape=2) # coefficients on family income

    s0 = a[0] + b[0] * family_income
    s1 = a[1] + b[1] * family_income
    s2 = np.zeros(N) # pivot
    s = pm.math.stack([s0, s1, s2]).T

    p_ = pt.special.softmax(s)
    career_obs = pm.Categorical("career", p=p_, observed=career)

    idata_11_14 = pm.sample(1000, tune=2000, target_accept=0.9, random_seed=RANDOM_SEED)

```

However, the model fails to sample.

```auto
Logp initial evaluation results:
{'a': -2.79, 'b': -1.9, 'career': -inf}
You can call `model.debug()` for more details.

```

When I call `m11_14.debug()` I can see that `p_` doesn’t add up to 1

```auto
point={'a': array([0., 0.]), 'b': array([0., 0.])}

The variable career has the following parameters:
0: Softmax{axis=None} [id A] <Matrix(float64, shape=(500, 3))>
 └─ Transpose{axes=[1, 0]} [id B] <Matrix(float64, shape=(500, 3))>
    └─ Join [id C] <Matrix(float64, shape=(3, 500))>
       ├─ 0 [id D] <Scalar(int8, shape=())>
       ├─ Add [id E] <Matrix(float64, shape=(1, 500))>
       │ ├─ ExpandDims{axes=[0, 1]} [id F] <Matrix(float64, shape=(1, 1))>
       │ │ └─ Subtensor{i} [id G] <Scalar(float64, shape=())>
       │ │ ├─ a [id H] <Vector(float64, shape=(2,))>
       │ │ └─ 0 [id I] <uint8>
       │ └─ Mul [id J] <Matrix(float64, shape=(1, 500))>
       │ ├─ ExpandDims{axes=[0, 1]} [id K] <Matrix(float64, shape=(1, 1))>
       │ │ └─ Subtensor{i} [id L] <Scalar(float64, shape=())>
       │ │ ├─ b [id M] <Vector(float64, shape=(2,))>
       │ │ └─ 0 [id I] <uint8>
       │ └─ [[0.145023 ... 61251361]] [id N] <Matrix(float64, shape=(1, 500))>
       ├─ Add [id O] <Matrix(float64, shape=(1, 500))>
       │ ├─ ExpandDims{axes=[0, 1]} [id P] <Matrix(float64, shape=(1, 1))>
       │ │ └─ Subtensor{i} [id Q] <Scalar(float64, shape=())>
       │ │ ├─ a [id H] <Vector(float64, shape=(2,))>
       │ │ └─ 1 [id R] <uint8>
       │ └─ Mul [id S] <Matrix(float64, shape=(1, 500))>
       │ ├─ ExpandDims{axes=[0, 1]} [id T] <Matrix(float64, shape=(1, 1))>
       │ │ └─ Subtensor{i} [id U] <Scalar(float64, shape=())>
       │ │ ├─ b [id M] <Vector(float64, shape=(2,))>
       │ │ └─ 1 [id R] <uint8>
       │ └─ [[0.145023 ... 61251361]] [id N] <Matrix(float64, shape=(1, 500))>
       └─ [[0. 0. 0. ... 0. 0. 0.]] [id V] <Matrix(float64, shape=(1, 500))>
The parameters evaluate to:
0: [[0.00066667 0.00066667 0.00066667]
 [0.00066667 0.00066667 0.00066667]
 [0.00066667 0.00066667 0.00066667]
 ...
 [0.00066667 0.00066667 0.00066667]
 [0.00066667 0.00066667 0.00066667]
 [0.00066667 0.00066667 0.00066667]]
This does not respect one of the following constraints: 0 <= p <=1, sum(p) = 1

```

```auto
Last updated: Tue Jan 09 2024

Python implementation: CPython
Python version : 3.11.7
IPython version : 8.19.0

numpy : 1.26.3
matplotlib: 3.8.2
scipy : 1.11.4
pandas : 2.1.4
pytensor : 2.18.4
arviz : 0.17.0
pymc : 5.10.3

Watermark: 2.4.3

```

I suspect that I am treating the Pytensor variables incorrectly. Any guidance on this would be greatly appreciated!

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 10, 2024, 7:14am UTC](https://discourse.pymc.io/t/multinomial-probabilities-dont-add-to-1/13609/2 "2024-01-10T07:14:15Z")

</div>

By default Softmax works across all axis. You probably need `axes=-1`

---

<div class="post-metadata">

**Author:** ![yahrMason](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/yahrmason/32/5826_2.png) [@yahrMason](https://discourse.pymc.io/u/yahrMason)\
**Post date:** [January 10, 2024, 10:07pm UTC](https://discourse.pymc.io/t/multinomial-probabilities-dont-add-to-1/13609/3 "2024-01-10T22:07:05Z")

</div>

> [@ricardoV94](#):
>
> `axes=-1`

Thanks, Ricardo. This is exactly what I needed. Silly oversight on my part.