Bayes' theorem applied to a binomial variable

Question

I just got an unexpected result, and wanted to check if I made any mistakes in my calculations, or if you can point me to posts/literature on the subject I describe below.

Suppose you sample $n$ random individuals from a large population (so you can assume replacement although you don't actually replace), and test them for a given feature that can be present or absent.
You find that the feature is present in $s$ individuals out of the $n$ tested.

Normally I would conclude that the expected fraction ($p_0$) of individuals with the feature in the overall population is $\frac s n$.

Concerning the distribution of $p_0$, I used the binomial distribution for the probability density:

$P(s \ out \ of \ n | p_0 = p) = \binom n s p^s (1-p)^{n-s}$

and Bayes' theorem with a uniform prior:

$P(p_0 = p | s \ out \ of \ n) = \frac {P(s \ out \ of \ n | p_0 = p)}{P(s \ out \ of \ n)} = \frac {P(s \ out \ of \ n | p_0 = p)}{\int_0^1 P(s \ out \ of \ n | p_0 = p) \cdot dp} $

Given that:

$\int_0^1 P(s \ out \ of \ n | p_0 = p) dp = \int_0^1 \binom n s p^s (1-p)^{n-s} \cdot dp = \frac 1 {n+1}$

it follows that:

$P(p_0 = p | s \ out \ of \ n) = (n+1) \binom n s p^s (1-p)^{n-s} $

which is the distribution I was looking for.

I tried plotting it for $p \in [0,1]$ with some numerical values of $n, s$, and the curve seemed OK. It was narrower when $n$ was larger, and it looked more symmetrical when $\frac s n$ was far from the extremes $0,1$.
Differentiating it w.r.t. $p$ and solving for $p$ gave $p = \frac s n$, which I took as an indication that the 'most frequent' value of $p_0$ was indeed the assumed one.

However, when I applied the 'expected value' method:

$E[p_0] = \int_0^1 P(p_0 = p | s \ out \ of \ n) \cdot p \cdot dp = \int_0^1 (n+1) \binom n s p^s (1-p)^{n-s} \cdot p \cdot dp = \frac {s+1}{n+2}$

This is the same as $\frac s n$ when $s = \frac n 2$, whereas in all other cases it gives a 'skewed' result.
E.g. for $s=2, n=6$, it gives $\frac 3 8$ instead of the 'assumed' $\frac 2 6 = \frac 1 3$.

Do you think this makes sense, or is it instead a consequence of my (wrong?) calculations or Bayesian prior?

---

EDIT: prior resulting in an expected value equal to the observed % ($\frac s n$)

Starting from the concept that a uniform Bayesian prior might not be appropriate, as suggested by user Jarle Tufto in the comments, I examined a different prior.

In particular, suppose that percentages (here the fraction $p_0$) are not uniformly distributed, whereas their logit ($u$) is:

$u=ln(\frac p {1-p})$

This implies that:

$P(p_0=p) = \frac 1 {p (1-p)}$

(the way I obtained this is quite nebulous even to me, and I will need to post a new question for that).

Reapplying Bayes' theorem with this new prior gives:

$P(p_0 = p | s \ out \ of \ n) = \frac {P(s \ out \ of \ n | p_0 = p) P(p_0=p)}{\int_0^1 P(s \ out \ of \ n | p_0 = p) P(p_0=p) \cdot dp} $

The denominator integrates to:

$\int_0^1 P(s \ out \ of \ n | p_0 = p) P(p_0=p) \cdot dp = \int_0^1 \binom n s p^s (1-p)^{n-s} \frac 1 {p (1-p)} \cdot dp = \frac n {s (n-s)} $

Thus:

$P(p_0 = p | s \ out \ of \ n) = \frac {\binom n s p^s (1-p)^{n-s} \frac 1 {p (1-p)}}{\frac n {s (n-s)}} = \binom {n-1} {s-1} (n-s) p^{s-1} (1-p)^{n-s-1}$

And the expected value is:

$E[p_0] = \int_0^1 P(p_0 = p | s \ out \ of \ n) \cdot p \cdot dp = \int_0^1 \binom {n-1} {s-1} (n-s) p^{s-1} (1-p)^{n-s-1} \cdot p \cdot dp = \frac s n$

However now it's the maximum value that does not seem to make much sense: by differentiating and solving for $p$, I get:

$p = \frac {s-1}{n-2}$

:/

Answer 1

You start using an Uniform prior, equivalent to a $Be(1,1)$. This two 1's can be understood as pseudo observations of "previous" experiments, let's say, having obtained 1 Yes, 1 No previously.

Now let's make an experiment where we obtain 2 Yes and 4 No, as in your example. The posterior is now a $Be(1+2,1+4)=Be(3,5)$.

The mode of the posterior is equivalent to the MLE in the frequentist method, this is, 2/6. Because you used an Uniform prior.

E.g. for $s=2, n=6$, it gives $\frac 3 8$ instead of the 'assumed' $\frac 2 6 = \frac 1 3$.

The expected value of the posterior is the MLE if we are taking into account the pseudo observations, this is, (2+1)/(6+1+1)=2/8.

In the second case I would say you are using the Haldane prior, $Be(0,0)$, where there are no "previous" observations. So, the expected value of the posterior is now the frequentist MLE.