# How do I compute the estimated values of x for a beta-binomial distribution?

I understand how to set up a binomial probability distribution. I'm trying to extend my understanding to the beta-binomial. On Wikipedia, there is a beta-binomial example given at https://en.wikipedia.org/wiki/Beta-binomial_distribution. I can follow the computations to get m1 and m2 and from there (6.23 and 42.31 respectively, n=12), hat-alpha and hat-beta (using the method of moments -- 34.1350 and 31.6085 respectively). I can reproduce in the example table the values for the row "Predicted (Binomial p = 0.519215)", but I can't figure out how to reproduce the values for the row "Predicted (Beta-Binomial)". How were those values computed?

Here's the data from the table (the number of male children among the first 12 children of family size 13 in 6115 families taken from hospital records in 19th century Saxony, so x ranges from 0-12):

Number of males:

0 1 2 3 4 5 6 7 8 9 10 11 12


Observed:

3 24 104 286 670 1033 1343 1112 829 478 181 45 7


Beta:

2.3 22.6 104.8 310.9 655.7 1036.2 1257.9 1182.1 853.6 461.9 177.9 43.8 5.2


Your example data is contained in the R package vcd, so I will use it from there. First, the betabinomial distribution has two parameters which we can estimate by maximum likelihood. Let us do that in R:

data(Saxony, package="vcd")

library(bbmle)
# negative loglikelihood function:

bb_nloglik0 <- function(x) {
function(alpha, beta) {
by=1), length(x)-1, alpha, beta, log=TRUE)*x)
}
}

bb_nloglik <- bb_nloglik0(Saxony)


Then finding the maximum likelihood estimates (by minimizing the negative log likelihood defined above)

mod <- bbmle::mle2(bb_nloglik, list(alpha=50, beta=50))
mod.prof <-bbmle::profile(mod)
coef(mod)
alpha     beta
34.10599 31.58112


For good measure we can also find confidence intervals by likelihood profiling, even if you didn't ask:

confint(mod.prof)
2.5 %   97.5 %
alpha 27.38399 44.53472
beta  25.35719 41.23744


And, finally, we can calculate the expected number of boys under the estimated beta-binomial model, by calculating the point probabilities and multiplying by N:

    round(rbind(Obs=Saxony,  Exp=extraDistr::dbbinom(0:12, 12,
coef(mod)[1], coef(mod)[2])*N), 1)
0    1     2     3     4      5    6      7     8     9    10   11  12
Obs 3.0 24.0 104.0 286.0 670.0 1033.0 1343 1112.0 829.0 478.0 181.0 45.0 7.0
Exp 2.3 22.6 104.8 310.9 655.7 1036.2 1258 1182.2 853.6 461.9 177.9 43.8 5.2