How to marginalize latent discrete variables in CDMs (cognitive diagnosis models)?

Questions v3
6

I would suggest you try and write the correct model for a single participant and a single question. That should be easy to check if you are getting the right thing. Then gradually expand to multiple questions for a single participant, and finally to the whole dataset.

It will also be easier for folks here to provide feedback on the intermediate versions.

An advantage of the smaller scenario is that you can include the discrete variables before any marginalization. You should be able to sample such a model.

If you write them as separate variables (instead of as a vector): You can even use MarginalModel — pymc_experimental 0.0.14 documentation. This allows you to check if the logp you are getting with your manual marginalization matches with the one from the MarginalModel.

Here is a simplified scenario notebook: Google Colab

import pymc as pm
import pymc_experimental as pmx
import numpy as np

q_skills = np.array([
    [1, 0, 1],
    [1, 1, 0],
    [0, 0, 1],
    [1, 1, 1],
])  # Skills needed to answer q1-q4 correctly

with pmx.MarginalModel(coords={"skills": range(3), "q": range(4)}) as m:
    # We need to create separate alpha's so that we can use `sum` and still marginalize the variables
    alphas = [pm.Bernoulli(f"alpha_{i}", p=0.3) for i in range(3)]
    s = pm.Beta("s", 1, 5)
    g = pm.Beta("g", 2, 1)

    # shapes: sum((q, skills) * (skills,), axis=-1) -> (q,)
    all_skills = pm.math.sum(pm.math.stack([alphas]) * q_skills, axis=-1)
    prob_correct = pm.math.where(all_skills, s, g)
    obs = pm.Bernoulli("obs", p=prob_correct, observed=[1, 1, 0, 1])

# Unmarginalized probability
print(m.point_logps())
# {'alpha_0': -0.36, 'alpha_1': -0.36, 'alpha_2': -0.36, 's': -1.09, 'g': -1.22, 'obs': -2.32}

m.marginalize(alphas)

# Marginalized probability: alphas are gone
print(m.point_logps()) 
# {'s': -1.09, 'g': -1.22, 'obs': -3.27}

MarginalModel will not be able to cope with the full model here, because of the interaction between the distinct Bernoulli variables (that’s why we had to create them separately). However it can be useful to check you are defining the right model

1 Like