Interesting. You’re right but I’m not quite sure how to get what you want. For anyone else as confused as me, you can reproduce the desired result with the following code to see that the resulting distributions are different.
import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt
prior_draws = np.random.randn(100)
traces = []
z_observation = 2.0
for mu in prior_draws:
with pm.Model() as individual_model:
y = pm.Normal('y', mu=mu)
z = pm.Normal('z',mu=y,observed=z_observation)
traces.append(pm.sample(progressbar=False))
individual_trace = np.concatenate([trace['y'].ravel() for trace in traces])
with pm.Model() as combined_model:
x = pm.Normal('x')
y = pm.Normal('y', mu=x)
z = pm.Normal('z',mu=y,observed=z_observation)
combined_trace = pm.sample(draws=20000)
plt.hist(individual_trace,bins=100, density=True,alpha=0.5)
plt.hist(combined_trace['y'],bins=100, density=True,alpha=0.5);