Reconciliation via Poisson Sampling

Hello, Ricardo.

>>> pm.__version__
# '5.10.2'

I changed as you proposed but it is still a up and running example.

>>> import pytensor.tensor as pt
>>> pt.matmul(S_2, l3)
# DropDims{axis=1}.0

And for the complete application

lvl3 = np.array([0,1,1,0,1,1,2,1])
lvl2 = np.array([1,1,2,3])
lvl1 = np.array([4,6])
lvl0 = np.array([10])

S_2 = np.array([[1,1,0,0,0,0,0,0],[0,0,1,1,0,0,0,0],[0,0,0,0,1,1,0,0],[0,0,0,0,0,0,1,1]])
S_1 = np.array([[1,1,0,0],[0,0,1,1]])
S_0 = np.array([1,1])

with pm.Model() as model:
    l3 = pm.Poisson('l3', mu=lvl3, shape=lvl3.shape)
    l2 = pm.Poisson('l2', mu=lvl2, shape=lvl2.shape, observed=pt.matmul(S_2, l3))
    l1 = pm.Poisson('l1', mu=lvl1, shape=lvl1.shape, observed=pt.matmul(S_1, l2))
    Y = pm.Poisson('Y', mu=lvl0[0], observed=pt.matmul(S_0, l1))
    trace = pm.sample(draws=5000, tune=3000)

trace.posterior.l3.mean(['chain', 'draw']).values
# array([0.     , 1.06295, 1.10495, 0.     , 1.0249 , 1.0223 , 2.1539 , 1.0237 ])