Manuscript Error: Conceptual Error with Problem 1M21 ?

[ikarosilva]

Description

Problem 1M21 asks us to estimate the number of customer "n" visiting a store, given the number of purchases made at the store, "Y", and the purchasing rate theta.

Both my solution and the solution provide the book (https://bayesiancomputationbook.com/solutions/chp_01.html#m21), give the following graph:

https://bayesiancomputationbook.com/_images/chp_01_71_1.png

However, if you look careful (for any case with Y>0), the HDI lower bound should be Y. for example at case Y=5, theta=0.2 the HDI lower bound is lower than Y. Unless I am missing something, this does not make sense. Because for Y=5 for example, we know that at least 5 customers made purchases and thus visited the store, so that n>=5.

Maybe the problem is with the choice of like-likelihood (instead of Binomial should be something like discrete Pareto, or constrained it so that it is 0 for values<Y)?

Does this makes sense, or am I missing something obvious?
Thanks.

Reference

Author Section

Do not close until

• Added to Errata
• Fixed in Open Access
• Fixed in latex source

[ikarosilva]

Update:
I think the solution below makes more intuitive sense.

#Solution
ys = [[0], [5], [10]] #Note: Incorrect in previous solution because we were only passing ones
θs = [0.2, 0.5]
traces = []
params = []
for θ in θs:
for y in ys:
with pm.Model() as model:
n = pm.Poisson("n", 10) #Note: 4.5 will give sampling error, makes sense if the purchases are higher than trials(?)
#n = pm.DiscreteUniform("n",lower=0,upper=20)
Y = pm.Binomial("Y", n, θ, observed=y)
trace = pm.sample()
traces.append(trace)
params.append((y, θ))

_, axes = plt.subplots(2, 3, figsize=(12, 8))
axes = axes.ravel()

for trace, ax, param in zip(traces, axes, params):
az.plot_posterior(trace, ax=ax)
ax.set_title(f"Y={np.sum(param[0])}, θ={param[1]}")

for trace, ax, param in zip(traces, axes, params):
az.plot_posterior(trace, ax=ax)
ax.set_title(f"Y={np.sum(param[0])}, θ={param[1]}")

The two key changes are (commented above):

1. We require the prior to be pm.Poisson(10)
   #Avoids "SamplingError: Initial evaluation of model at starting point failed!" because our data in some cases our successes are much higher than the mean number of trials (?).

2. Changed "ys" from:
   ys = [0, np.ones(5), np.ones(10)] #I think this a mistake in the original solution?
   to:
   ys = [0, [5], [10]]