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

Questions v5
1

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[0],
            p[0] / (p[0] + p[1]),
        ),
        pm.math.eq(
            q1[1],
            p[1] / (p[0] + p[1]),
        ),
        pm.math.eq(
            q2[0],
            p[1] / (p[1] + p[2]),
        ),
        pm.math.eq(
            q2[1],
            p[2] / (p[1] + p[2]),
        ),
        pm.math.eq(
            q3[0],
            (p[0] + p[1]) / (p[0] + p[1] + p[2]),
        ),
        pm.math.eq(
            q3[1],
            p[2] / (p[0] + p[1] + p[2]),
        ),
    ]

    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'

I suspect I’m going about this in completely the wrong way. What is a better approach?

Thanks in advance!

2

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[0] / (p[0] + p[1])),
        pm.Deterministic("q1_1", p[1] / (p[0] + p[1])),
    ]
    q2 = [
        pm.Deterministic("q2_0", p[1] / (p[1] + p[2])),
        pm.Deterministic("q2_1", p[2] / (p[1] + p[2])),
    ]

    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