# How to create model with series of multinomials, each of which has observations on only a subset of components?

Hi!

I’m trying to infer the values of proportions `p`, where my observations are integer counts of subsets `q` of two or more of those proportions.

``````import arviz as az
import numpy as np
import pymc as pm

basic_model = pm.Model()

with basic_model:
conc = 1.0

p = pm.Dirichlet("p", a=3 * [conc])

q1 = pm.Dirichlet("q1", a=2 * [conc])
q2 = pm.Dirichlet("q2", a=2 * [conc])
q3 = pm.Dirichlet("q3", a=2 * [conc])

pm.Multinomial("m1", p=q1, n=15, observed=[5, 10])
pm.Multinomial("m2", p=q2, n=22, observed=[11, 11])
pm.Multinomial("m3", p=q3, n=40, observed=[32, 8])

eq = [
pm.math.eq(
q1,
p / (p + p),
),
pm.math.eq(
q1,
p / (p + p),
),
pm.math.eq(
q2,
p / (p + p),
),
pm.math.eq(
q2,
p / (p + p),
),
pm.math.eq(
q3,
(p + p) / (p + p + p),
),
pm.math.eq(
q3,
p / (p + p + p),
),
]

for i, e in enumerate(eq):
pm.Potential(
f"constraint{i}",
pm.math.switch(
e,
0,
-np.inf,
),
)
result = pm.sample(draws=1000, tune=1000, return_inferencedata=True)

``````

In this example, each of the subsets is modelled as a multinomial with its own Dirichlet prior. I try to provide the relationship between the overall `p` mixture and the `q` subsets using equality constraints; notably, the `q3` subset has its first component being a sum of two `p` components. This results in the following exception:

``````File "/data/jeff/mambaforge3/envs/happyjeff/lib/python3.8/site-packages/aesara/tensor/elemwise.py", line 621, in transform
if isinstance(r.type, (NullType, DisconnectedType)):
AttributeError: 'float' object has no attribute 'type'
``````

As I suspected, I was being dumb – the solution is simpler than I thought. Here’s working code that achieves what I wanted.

``````import arviz as az
import numpy as np
import pymc as pm

basic_model = pm.Model()

with basic_model:
conc = 1.0

p = pm.Dirichlet("p", a=3 * [conc])

q1 = pm.Dirichlet("q1", a=2 * [conc])
q2 = pm.Dirichlet("q2", a=2 * [conc])
q3 = pm.Dirichlet("q3", a=2 * [conc])

q1 = [
pm.Deterministic("q1_0", p / (p + p)),
pm.Deterministic("q1_1", p / (p + p)),
]
q2 = [
pm.Deterministic("q2_0", p / (p + p)),
pm.Deterministic("q2_1", p / (p + p)),
]

pm.Multinomial("m1", p=q1, n=150, observed=[50, 100])
pm.Multinomial("m2", p=q2, n=220, observed=[110, 110])
result = pm.sample(draws=1000, tune=1000, return_inferencedata=True)
print(az.summary(result)["mean"])
``````
1 Like