# Reproducing a didactic example of Lindley (1993)

Lindley (1993) discusses the following mixed discrete and continuous prior for the tea tasting lady experiment, where $$\pi$$ is probability of a correct classification: $$p(\pi=0.5) = 0.8$$ (discrete part) and $$p(\pi | \pi > 0.5) = 1.6 (1-\pi)$$ (continuous part). The prior reflects a skeptic stance towards Muriel Bristol's tea tasting ability, but not entirely discards the possibility pf $$\pi > 0.5$$. Lindley states the posterior probability for $$P(\pi = 0.5)$$ for several combinations of correct and wrong classifications, e.g., 5 right, one wrong: $$P(\pi = 0.5 | 5R, 1W) = 0.59$$.

The math for this uncommon mixed prior is straightforward: By Bayes' theorem $$P(\pi = 0.5 | 5R, 1W) = \frac{P(5R,1W|\pi=0.5)p(\pi = 0.5)}{P(5R,1W)}$$ and the denominator leads to an integral of the continuous part of the posterior, which I evaluate numerically. However, I failed to reproduce the numbers, coming up with $$P(\pi = 0.5 | 5R, 1W) = 0.605$$. The discrepancy for $$P(\pi = 0.5 | 6R)$$ is even higher. Also for the other example in Lindley (1993) [wine tasting] I find the maximum-a-posteriori estimate closer to 0.8 than to 0.75 (what Lindley states).

prior_tea_discrete <- function(p){(p==0.5) * 0.8}
prior_tea_cont <- function(p){1.6*(1-p)}
y1 <- c(1,1,1,1,1,0)
lik1 <- function(p){choose(length(y1),sum(y1)) * p^sum(y1) * (1-p)^{length(y1)-sum(y1)}}
lik1 <- Vectorize(lik1)
y2 <- c(1,1,1,1,1,1)
lik2 <- function(p){choose(length(y2),sum(y2)) * p^sum(y2) * (1-p)^{length(y2)-sum(y2)}}
lik2 <- Vectorize(lik2)

post_tea_cont_1 <- function(p)lik1(p)*prior_tea_cont(p)
v1 <- integrate(post_tea_cont_1, 0.5, 1)$value post_tea_discrete_1 <- lik1(0.5)* prior_tea_discrete(0.5) post_prob_discr1 <- post_tea_discrete_1/(post_tea_discrete_1 + v1) # 0.605, slightly higher than Lindley's (1993) 0.59 post_tea_cont_2 <- function(p)lik2(p)*prior_tea_cont(p) v2 <- integrate(post_tea_cont_2, 0.5, 1)$value
post_tea_discrete_2 <- lik2(0.5)* prior_tea_discrete(0.5)
post_prob_discr2 <- post_tea_discrete_2/(post_tea_discrete_2 + v2)
# 0.312, substantially higher than Lindley's 0.23


## Reference

Lindley, D.V. (1993). The analysis of experimental data: The appreciation of tea and wine. Teaching Statistics, 15(1), 22-25.

Yes, and like you, for five right and one wrong, I get (using Lindley's $$P$$ rather than your $$\pi$$) $$\mathbb P(P=\tfrac12 \mid RRRRRW) = \frac{ 0.8 P^5(1-P)\Big|_{P=1/2}}{0.8 P^5(1-P)\Big|_{P=1/2} + \int\limits_{1/2}^1 1.6(1-P)P^5(1-P)\, dP} \\= \frac{112}{185} \approx 0.6054$$
and for six right and none wrong, I get $$\mathbb P(P=\tfrac12 \mid RRRRRR) = \frac{ 0.8 P^6\Big|_{P=1/2}}{ 0.8 P^6\Big|_{P=1/2} + \int\limits_{1/2}^1 1.6(1-P)P^6\, dP} \\= \frac{112}{359} \approx 0.3120.$$
For the wine tasting, Lindley used a prior of $$48(1-P)(P-\frac12)$$ on $$[\frac12,1]$$. By assigning a zero prior probability to $$P<\frac12$$ he seems to have forgotten his separate advocacy of "Cromwell's rule".
In the 5R,1W case, I get a posterior density of $$\frac{258048}{275}(1-P)^2(P-\frac12)P^5$$ on $$[\frac12,1]$$ with a maximum, i.e. a mode or MAP, at $$P=\frac{\sqrt{41}+19}{32} \approx 0.794$$, presumably the same as your calculation, and a mean of $$\frac{1058}{1375} \approx 0.769$$.
In the 6R,0W case, I get a posterior density of $$\frac{258048}{1291}(1-P)(P-\frac12)P^6$$ on $$[\frac12,1]$$ with a maximum at $$P=\frac{\sqrt{57}+21}{32} \approx 0.892$$, higher than the roughly $$0.87$$ suggested in Lindley's Figure 1, and a mean of $$\frac{5397}{6455} \approx 0.836$$.