MLE vs Binomial - Beta conjugate pairs

Question

I'm trying to understand how to use the conjugate relationship between the Binomial and Beta distributions to update the parameters of the Beta distribution. Specifically, I am imagining a series of binomial trials with a prior $\text{Beta}(2,2)$ on the probability of success. As I understand it, I add the number of successes to the first parameter of the Beta distribution ($\alpha$), and the number of failures to the second ($\beta$).

However, this seems to give an answer that is far too precise, and also wrong! Maximum likelihood fit to the proportions gives a much more sensible answer. What am I missing?

Simple R code to reproduce

set.seed(1234)
# simulate a variable number of sample sizes and successes
ns <- rpois(45, 200)
xs <- rbinom(45, ns, 0.1)

# fit a beta distribution with maximum likelihood
beta_fit <- fitdistrplus::fitdist(xs/ns, "beta", start = list(shape1 = 2, shape2 = 2))
beta_fit
#> Fitting of the distribution ' beta ' by maximum likelihood 
#> Parameters:
#>         estimate Std. Error
#> shape1  21.36964   4.470611
#> shape2 176.69459  37.348557

hist(xs/ns, freq = FALSE, ylim = c(0,20), xlim = c(0, .3), col = "lightgrey")

# curve from maximum likelihood
curve(dbeta(x, beta_fit$estimate[[1]], beta_fit$estimate[[2]]), add = TRUE, col = "orange", lwd = 3)
# curve from conjugate priors
curve(dbeta(x, sum(xs) + 2, sum(ns - xs) + 2), add = TRUE, col = "darkgreen", lwd = 3)

legend(.2, 16, legend=c("MLE", "conjugate"),
       col=c("orange", "darkgreen"), 
       lty=1, cex=1, box.lty=0, lwd = 3)

^{Created on 2019-11-11 by the reprex package (v0.3.0)}

Answer 1

Beta is a conjugate prior for binomial, not the other way around.

In beta-binomial model, you assume beta prior for the probability of success $p$. The distribution of $p$ a posteriori follows beta distribution. The posterior predictive distribution of the data follows beta-binomial distribution.

If you wanted to compare the Bayesian estimate to the frequentist estimate, you should rather compare it to the confidence interval around the parameter (so, something roughly $\sqrt{N}$ times smaller).

Moreover, in your example you are using varying binomial sample sizes, while in binomial distribution this is a fixed parameter, so beta-binomial model in the simplest form is not appropriate in here.

Answer 2

There is a saying in computer science that the hardest computer language to learn is your second one. I believe that is true in statistics too. The most challenging probability method to learn is your second.

Let’s split these problems apart into a null hypothesis solution and a Bayesian method and, instead, we will try and mentally work through what the two types of tools are trying to do.

Let’s begin with the method of maximum likelihood and with the minimum variance unbiased estimator. For the binomial, the MLE is also the MVUE. I bring this up because it will create a potential issue with how we think about the Bayesian problem.

If you are randomizing the sample size, then neither method you are using is appropriate because the choice of sample size is also a random process with a parameter that needs to be estimated. You are adding a gigantic problem into what you are doing called “optional stopping.” You do not want to do that. If you add randomness anywhere into a process that would not have that randomness, then you need to account for that in some way.

For both cases, we are going to fix the sample size at 10 each. We are also going to imagine four parallel universes out of the infinite number possible in the parameter space of $\Theta=[0,1]$. While the integrations will consider all possible worlds, it is impossible to visualize an infinite number of options.

In universe one, $\theta=.5$ and successes equal 5; in universe two, $\theta=.41$ and successes equal 5; in universe three $\theta=.5$ and successes equal four; while in universe four, $\theta=.41$ and success equal four.

In all universes, you believe that the coin is roughly fair, and on the Bayesian calculations, you assign a prior distribution of $\Pr(\theta)=B(2,2)$. What you are looking at is not a graph of heads or tails, but of possible parameter locations. This is not about successes or failures, it is about parameter location. What you are doing is underweighting any world that is distant from .5 if you are really in one of those worlds.

In all worlds, this graph is your belief system before seeing any data.

In universe one, the true long-run mass function for ten draws is as follows. Five successes were observed. $\hat{\theta}_{MLE/MVUE}=.5$ Your 95% confidence interval is $[.19,.81]$. If you would “plug-in” the parameter estimate as if true, then your mass would match the true long-run distribution below.

There are two logical Bayesian point estimators based on your use of the MLE or the MVUE. The MLE minimizes the all-or-nothing loss function, which would map to the Bayesian posterior mode if the loss function were preserved. The MVUE minimizes quadratic loss. The Bayesian mapping of that loss function is the posterior mean. The split in the Bayesian result will require you to choose a loss function if you need a point estimator of $\theta$.

The Bayesian posterior mean is .5, and the posterior mode is .5. While the Frequentist confidence interval depends on the type of point estimator that you choose, the Bayesian one does not. The Bayesian interval is $[.251,.749]$. The Bayesian density function for your beliefs is graphed below.

In world two, the true long-run mass function is below. Nonetheless, because the sample is the same, the parameter estimates remain the same as do the intervals and the graphs.

This is a good place to look at the Bayesian predictive distribution, which is what you were trying for above. It is the predicted number of successes based on the data and prior. It is below.

The prediction isn’t considering a point estimate as valid. The prediction is the weighted average binomial distribution over the set of all possible parameter values weighted by the posterior. Note the density values on the left. The Bayesian prediction is far more agnostic about the future than is the one created from plugging in the MLE/MVUE.

The Bayesian estimator is still giving weight to unlikely parameter values. It is possible that the actual value of the parameter is .6 with a 5:5 success to failure ratio. The Bayesian estimate converges to the asymptotic distribution as the sample size goes to infinity.

Now let us consider the universes where we see four successes. The Frequentist point estimator is $\hat{\theta}_{MLE/MVU}=.4$ with a 95% confidence interval based on the normal approximation of $[.0964,.7036]$. The Bayesian posterior mean is .4286, while the posterior mode is .4167. The Bayesian 95% highest density region credible interval is $[.192,.684]$. The graph of your new beliefs is below.

Note that the mode is not co-located with the MLE/MVUE. You can see that in the graph of posterior beliefs.

The Bayesian predictive distribution is below.

Now let’s summarize what is above. The Frequentist binomial distribution is valid only when $\frac{k}{n}=\theta$ although it will tend to be close. Assuming that the estimator is also the real value creates steep predictions. What if the actual value had been $\theta=.5$, but there had been 9 observations. It would be impossible for the Frequentist plot to match reality.

The Bayesian strength is the Bayesian weakness. If your beliefs about reality are actually stated in the prior and you have reasonable beliefs, then you will get gains in estimation because it improves precision. If your beliefs are nonsense, then your posterior will be, at least in part, based on nonsense.

What if $\theta=.9$, then it would take quite a bit of time to learn the true value of the parameter.

If you saw 9 successes, your posterior beliefs would be in the graph below. The graph falls below 0 because of the smoothing method I used in geom_smooth. Sorry, it is a bit of an effort to fix it so I left it.

The Bayesian prediction on the ten next values is below.

The weakness of a Frequentist method is if there are zero or ten successes. For zero successes, the estimator is zero. The interval is the point containing zero. The prediction would be that every single value will be failures. It gives no weight at all to any other possible parameter value.