The solution to 1m21 seems wrong to me. The solution code lists an array of all ones, which would indicate k trials that each resulted in one success. But that is not what the question tells us: we observe one trial with each of Y total successes. The tricky part is that we have to guarantee that the quantity we are attempting to model, n (i.e. the number of trials), must be at least the number of observed successes. Here is my solution for this problem:
obs = [0, 5, 10]
thetas = [0.2, 0.5]
fig, axes = plt.subplots(
nrows=len(obs), ncols=len(thetas), figsize=(10, 8), sharex=False
)
for i, y in enumerate(obs):
for j, theta in enumerate(thetas):
with pm.Model() as model:
# Shifted Poisson formulation: n = n_failures + y
n_failures = pm.Poisson("n_failures", mu=4.5)
n = pm.Deterministic("n", n_failures + y)
# Likelihood
Y = pm.Binomial("Y", n=n, p=theta, observed=y)
# Sampling
step = pm.Metropolis(vars=[n_failures])
trace = pm.sample(
draws=3000,
tune=1000,
step=step,
random_seed=42,
progressbar=False,
return_inferencedata=True,
)
# Plot posterior for 'n' (kind="hist" renders clean discrete bars)
ax = axes[i, j]
az.plot_posterior(trace, var_names=["n"], ax=ax, kind="hist")
ax.set_title(f"Observed Y = {y}, θ = {theta}")
plt.style.use('seaborn-v0_8-whitegrid')
plt.tight_layout()
plt.show()
The solution to 1m21 seems wrong to me. The solution code lists an array of all ones, which would indicate k trials that each resulted in one success. But that is not what the question tells us: we observe one trial with each of Y total successes. The tricky part is that we have to guarantee that the quantity we are attempting to model, n (i.e. the number of trials), must be at least the number of observed successes. Here is my solution for this problem: