Hi,
Since eval statements are working fine for the likelihood (which means model would have a valid structure at that what ever random point eval choose) but fails to evaluate at the initial point (as suggested by the error stack involving point = model.initial_point()), something in your model does not make sense in the initial point. Indeed if you look at
pll[0].eval( {'delta': np.array([[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])})
it gives a value of [0]. Later on you use its first element p2 as
U[i, :p2-1]
which is a zero dimensional array which likely is the root of your problems (since the error also says something about not being able to broadcast zero dimensional array with sth else).
You can for instance try setting something like
initvals={"delta":np.array([[0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1,
0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1,
0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0]])}
and later
delta = pm.Bernoulli("delta",p = p, dims = ("sample", "specie"),
initval=initvals["delta"])
Now the model starts running though very slowly on my end. Don’t know why maybe bad initial point. But I also see that you have a peculiar structure here where
for i, p2 in enumerate(pll):
tem = pm.math.sum(pm.math.exp(U[i, :p2-1]))
concatenated_values = pm.math.concatenate((pm.math.exp(U[i, :p2-1]) / (1 + tem), [1 / (1 + tem)]))
if pm.draw(p2) == 0:
expend_con = pt.set_subtensor(expend_con[i,:], concatenated_values)
else:
expend_con = pt.set_subtensor(expend_con[i, ind_tem[i]], concatenated_values)#concatenated_values[0]
you update some of your values based on pm.draw(p). I wonder if there is a better way to achieve this (why not just p2==0 for instance? Also probably use switch statement to make it vectorized). Or maybe it is fine I don’t know but it is (in my limited experience) first time seeing a model structure with pm.draw. I will let more experienced folks comment on that.
ps: I have also noticed you are using slice sampling for everything but if you let pymc choose with trace = pm.sample(), it actually runs much faster.