I’m trying to understand how to fit a joint distribution in pymc3.

Here’s a toy model to illustrate my question. Suppose I have data y[2,N] drawn from a mixture model, with some mixing proportion, where the mixture component y_i = k determines two things:

- The mean of the Normal distribution that y[0,i] is sampled from (either theta1 or theta2).
- The rate of the Poisson distribution that y[1,i] is sampled from (either rate1 or rate2).

How do I implement such a model in pymc3?

My attempt is below (marginalizing out cluster assignments). I may be going down the wrong track since I don’t understand how to dealing with data where each observation has elements drawn from multiple distributions. The current code gives the error “too many indices for array” at line “normal_data = y[0,i]”.

def likelihood(theta, rate): def logp(y): normal_data = y[0,:] poisson_data = y[1,:] return (pm.Poisson.dist(rate).logp(poisson_data) + pm.Normal.dist(theta, 1).logp(normal_data)) return logp basic_model = pm.Model() with basic_model: theta1 = pm.Uniform('theta1', 0, 1) theta2 = pm.Uniform('theta2', 0, 1) rate1 = pm.Exponential('rate1', 1) rate2 = pm.Exponential('rate2', 1) mix_proportion = pm.Uniform('lambda', 0, 1) custom = pm.DensityDist('custom', likelihood(theta1, theta2)) y_obs = pm.Mixture('y_obs', mix_proportion, [custom(theta1, rate1), custom(theta2, rate2)], observed=data)

Thanks!