# Histogram of the MLE of the probability in binomial distribution and the plot of beta distributions

I have data with columns "y" and "n", which for this example can be "y" count of heads out of "n" coin flips. There are "i" rows ie "i" different coins: $$P(y_{i} \mid p_{i}, n_{i}) = \binom{n_{i}}{y_{i}}\: p_{i}^{y_{i}} \: (1-p_{i})^{n_{i}-y_{i}}$$ When I calculate the MLE of $$p_{i}$$, this comes out equal to $$\frac{y_{i}}{n_{i}}$$. Calculating the MLE of $$p_{i}$$ for each row of data, and plotting a histogram I get a "U" shape.

Separately, I believe binomial distributions of the kind of situation I have described can be modelled as beta distributions of choosing values of $$p_{i}$$. The form of the beta distribution is: $${p_{i}^{a - 1}(1 - p_{i})^{b-1} \over B(a, b)}$$ With $$a-1 = \sum_{i} y$$ and $$b-1 = \sum_{i} (n_{i} - y_{i})$$, the shape of this beta distribution does not match the "U" shaped histogram of the MLE of $$p_{i}$$. Why is this? If the beta distribution is the probability density of $$p$$ given a, the sum of all positive cases and b, the sum of all negative cases, then should it not have the same "U" shape as the histogram of the MLE of $$p$$?

EDIT: more info on the data, this can be implied from the histogram of $$\hat p_{i}$$ however the data seems to form two clusters, one with a high frequency of low $$\hat p_{i}$$ (so low probability of getting heads with $$y_{i}$$ being low compared to $$n_{i}$$ for that particular row) and a high frequency of high $$\hat p_{i}$$ (so high probability of getting heads with $$y_{i}$$ being high compared to $$n_{i}$$ for that particular row). And then a low frequency of $$\hat p_{i}$$ between the two extremes. This corresponds to the coins being either heavily biased towards heads or heavily biased towards tails with few coins in between.

• You are using the wrong Beta distribution (because you are applying a formula to solve a different problem). To model a U shape correctly, both $a$ and $b$ must be less than $1.$ You can estimate them with a full MLE.
– whuber
Commented Feb 7, 2022 at 18:21
• @whuber I thought there was only 1 beta distribution - this one: ${p_{i}^{a - 1}(1 - p_{i})^{b-1} \over B(a, b)}$? Also, how can $a$ and $b$ both be less than 1, if $a$ is the number of positive occurances and $b$ is $n-a$? Commented Feb 7, 2022 at 18:23
• That's the misapplication to which I was referring. You describe a situation in which the underlying values of $p$ are not fixed: they appear to vary randomly according to some distribution, and you have elected to model that distribution with one of the infinitely many Betas parameterized by some unknown $(a,b).$ This is very different than updating an estimate of a fixed value of $p.$
– whuber
Commented Feb 7, 2022 at 18:26
• @whuber Well $p_{i}$ is bounded between 0 and 1 as it is the probability of coin "i" landing on heads. In this situation it is essentially a measure of how biased the coin is. Also from plotting $\hat p_{i}$ I can see the distribution of p, so I guess my question then is how do I choose a beta distribution to model this? I know a and b should be < 1, however I thought a and b were $a=\sum_{i} y_{i}$ and $b=n-a$ respectively. What other values could they take? Commented Feb 7, 2022 at 18:33
• They will take whatever values maximize the likelihood of your observations. Because you describe a U-shaped histogram, that strongly indicates those values will both be less than $1.$ (Unless, perhaps, you mean an upside-down U?) If this isn't perfectly clear, then please include a small but characteristic subset of your data in your question to make sure everyone understands exactly what you are working with and that people who answer can use as an example.
– whuber
Commented Feb 7, 2022 at 18:36

Notice that $$\hat p_i = y_i/n_i$$ is the MLE of $$p_i$$ given the single observation $$y_i$$. It doesn't make sense to compare it to the posterior probability of $$p$$ given all the data.

The posterior probability of $$p_i$$ given that single observation is the Beta distribution you wrote with $$a = \alpha + y_i - 1$$ , $$b = \beta + n_i - y_i - 1$$ , where $$\alpha,\beta$$ are the prior parameters (which you can take to be 1 if you want a uniform distribution).

Or conversely you could compare the posterior distribution you wrote to the MLE of $$p$$ given all the data, $$\hat p = \sum y_i / \sum n_i$$.

It's also worth noting that in general the posterior probablilty and the distribution of the MLE are not the same thing, although they are related in some cases.

• Would comparing the posterior distribution to the MLE of p given all the data not result in a $\hat p$ value of around 0.5 seen as the data is clustered towards p=0.1 and p=0.9? Commented Feb 7, 2022 at 19:18
• Do you mean that the data come from a mixture of two populations with different binomial parameters ? Commented Feb 7, 2022 at 19:25

One approach is to "integrate out" the random term ($$p$$) and then maximize the resulting likelihood function. The joint density of $$Y_i$$ and $$p_i$$ (assuming you've decided upon a beta distribution for the random values of $$p$$) is

$$\frac{\binom{n_i}{y_i} p_i^{a+y_i-1} (1-p_i)^{b+n_i-y_i-1}}{B(a,b)}$$

Integrating over $$p_i$$ we have

$$\int_0^1 \frac{\binom{n_i}{y_i} p_i^{a+y_i-1} (1-p_i)^{b+n_i-y_i-1}}{B(a,b)} \, dp_i=\frac{\binom{n_i}{y_i} \Gamma (a+y_i) \Gamma (b+n_i-y_i)}{B(a,b) \Gamma (a+b+n_i)}$$

and the log of the likelihood becomes

$$\log(L)=\sum _{i=1}^n \left(\log \binom{n_i}{y_i} + \log\Gamma (a+y_i)+\log\Gamma(b+n_i-y_i) -\log (B(a,b)) - \log\Gamma(a+b+n_i)\right)$$

We can drop the $$\log \binom{n_i}{y_i}$$ term as it is a constant for this purpose. Using R one can perform the above steps:

# Generate some data
a <- 0.23
b <- 0.40
nsim <- 1000
set.seed(12345)
n <- rpois(nsim, 20)
p <- rbeta(nsim, a, b)
y <- rep(0, nsim)
for (i in 1:nsim) {y[i] <- rbinom(1, n[i], p[i])}

# Log of likelihood
logL <- function(par, n, y) {
a <- par[1]
b <- par[2]
sum(lgamma(a + y) + lgamma(b + n - y) - lbeta(a, b) - lgamma(a + b + n))
}

# Obtain maximum likelihood estimates
mle <- optim(c(0.25,0.25), logL, n=n, y=y, control=list(fnscale=-1), hessian=TRUE)

# Maximum likelihood estimates
mle$par # [1] 0.2443108 0.3808030 # Estimated standard errors cov <- solve(-mle$hessian)
(se.a <- cov[1,1]^0.5) # 0.01307452
(se.b <- cov[2,2]^0.5) # 0.02134923

# Plot of histogram and estimated density
hist(y/n, freq=FALSE)
lines(c(0:100)/100, dbeta(c(0:100)/100, mle$$par[1], mle$$par[2]))


• +1 -- but please fix the discrepancies between how your code computes $\log L$ (which is correct) and your formula (which is incorrect).
– whuber
Commented Feb 8, 2022 at 15:28
• @whuber Thanks for catching that error. I believe it's now fixed.
– JimB
Commented Feb 8, 2022 at 16:24