I’m new to Pymc3 and it seemed like a great tool for a problem I am encountering, but I am having trouble with the last bit of code.
Briefly, I have two multinomial distributions representing disease models (each disease has the same three objects present in imaging data, but in different frequencies). I would like to model an instance where both diseases are present in imaging data as mixture. This would be straight-forward if each observation (patch of the image) came from from disease A or B, but it is more realistic that each observation will be a mixture itself (some convolution of disease A and B).
Currently what I have is below:
Setting up our prior and underlying parameters:
#Number of image patches N=100 #Number of objects in each image patch objPerPatch = 10 #Initial guess for the fraction of the two diseases present beta = np.array([0.5, 0.5]) #True frequency of our diseases in image patches beta_true = np.array([0.7, 0.3]) #Frequency for three objects in each disease alphas_true12 = np.array([[0, 0.2, 0.8],[0.8, 0.2, 0]])
Simulating the data:
observed_data12 = np.zeros([N,3]) #Each image patch for simplicity is the same combination of the two multinomials # and will always the same object count for m in range(N): observed_data1 = np.random.multinomial(objPerPatch*beta_true, alphas_true12, size=1) observed_data2 = np.random.multinomial(objPerPatch*beta_true, alphas_true12, size=1) observed_data12[m] = observed_data1+observed_data2
Creating the model:
with pm.Model() as simple_comp_model: #Underlying probability we want to infer p1 = pm.Uniform('p', 0, 1) p2 = 1 - p1 p = T.stack([p1, p2]) #Get distribution of our multinomials obj_dist = [pm.Multinomial.dist(p=alphas_true12[i], n=objPerPatch, shape=3) for i in range(2)] #*This would be the approach if each image patch was assigned to one disease #disease_mix = pm.Mixture('disease_mix', w = p, comp_dists = obj_dist, observed=observed_data12) #Now get the distribution of this mixture? Unclear disease_mix = pm.Mixture.dist( w = p, comp_dists = obj_dist) #Treat problem as mixture of mixture (no bueno) x = pm.Mixture('x',w = [1,0], comp_dists = disease_mix,observed=observed_data12)
When I attempt this, I get TypeError: ‘Mixture’ object is not iterable
with simple_comp_model: trace = pm.sample(draws=3000,chains=1, tune=2000,discard_tuned_samples=True)
I’ve looked at pm.Potential and pm.DensityDist rather than pm.Mixture, but to no avail. I feel that I’m missing one last line to wrap it up, but I can’t figure it out.
Apologies if this problem has already been answered. I looked through other questions but none seemed to be quite the same.