Model not reproducing observed data?

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])