Thanks. I do not fully understand mixture models or how to construct marginalized models. Is there a tutorial that I should follow?
By the way, I end up changing Bernoulli
skills = pm.Bernoulli(‘skills’, p=0.5*np.ones((22, 7)), shape=(22,7))
To Beta distribution
skills = pm.Beta(‘skills’, alpha=2, beta=2, shape=(22,7))
And model is providing better estimates. I also felt perhaps I should have used beta in the first place as probability of skill could be between 0 or 1, rather than being binary.