In that case, it is probably better to model it with a Binomial distribution, maybe try this:
for i in range(4):
for j in range(4):
pm.Binomial('obs_{}_{}'.format(i,j), n=Ntrial[i,j], p=p[i,j], observed=data_raw[i,j])
In that case, it is probably better to model it with a Binomial distribution, maybe try this:
for i in range(4):
for j in range(4):
pm.Binomial('obs_{}_{}'.format(i,j), n=Ntrial[i,j], p=p[i,j], observed=data_raw[i,j])