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