As I say, it’s not clear to me what the problem was. I fixed it with manually broadcasting two numpy arrays which were 2nd and 3rd arguments to tt.where. Everything worked correctly before this fix when I set parameters like jm.B1 = 1 and call jm.image().eval(), but somehow failed when using this with pymc. After this fix I managed to sample this model and get correct results.