Beta-Binomial Likelihood for Grouped Success Counts

From Grouped Counts to a Likelihood

A grouped binomial observation records a number of successes and the number of trials that produced it. In the beta-binomial model, each group has its own success probability, drawn from a beta distribution; conditional on that probability, its count is binomial. Integrating out the group-specific probability gives a beta-binomial count. The likelihood therefore uses one beta-binomial probability for each observed group, not one binomial probability for the total number of successes pooled across groups.

Suppose there are $m$ independent groups. For group $i$, the fixed trial total is $n_i$, and the observed success count is $y_i$, where $n_i$ is a nonnegative integer and $y_i\in\{0,\ldots,n_i\}$. Assume a common pair of shape parameters $\alpha>0$ and $\beta>0$, with group-specific probabilities $P_i$ independently distributed as $\operatorname{Beta}(\alpha,\beta)$. Conditional on $P_i$, assume $Y_i\mid P_i\sim\operatorname{Binomial}(n_i,P_i)$.

Definition: The beta-binomial probability mass function for $y_i$ successes in $n_i$ trials is $$ P(Y_i=y_i)=\binom{n_i}{y_i}\frac{B(y_i+\alpha,n_i-y_i+\beta)}{B(\alpha,\beta)},\qquad y_i=0,\ldots,n_i, $$ where $B(a,b)=\int_0^1 p^{a-1}(1-p)^{b-1}\,dp$ for $a,b>0$. The grouped likelihood is the product of these probabilities over the independent groups.

The beta-binomial formula follows directly from the mixture. The binomial probability conditional on $P_i=p$ is $\binom{n_i}{y_i}p^{y_i}(1-p)^{n_i-y_i}$, and the beta density is $p^{\alpha-1}(1-p)^{\beta-1}/B(\alpha,\beta)$. Multiplying and integrating over $0<p<1$ gives

\begin{aligned} P(Y_i=y_i) &=\int_0^1 \binom{n_i}{y_i}p^{y_i}(1-p)^{n_i-y_i} \frac{p^{\alpha-1}(1-p)^{\beta-1}}{B(\alpha,\beta)}\,dp\\ &=\binom{n_i}{y_i} \frac{\int_0^1p^{y_i+\alpha-1}(1-p)^{n_i-y_i+\beta-1}\,dp} {B(\alpha,\beta)}\\ &=\binom{n_i}{y_i}\frac{B(y_i+\alpha,n_i-y_i+\beta)}{B(\alpha,\beta)}. \end{aligned}

Independence across groups now gives the likelihood. The trial totals are treated as fixed, so they enter the combinatorial factors and beta-function arguments but are not themselves modeled as random.

Likelihood: For observed counts $(y_1,\ldots,y_m)$ and fixed totals $(n_1,\ldots,n_m)$, $$ L(\alpha,\beta)=\prod_{i=1}^{m}\binom{n_i}{y_i} \frac{B(y_i+\alpha,n_i-y_i+\beta)}{B(\alpha,\beta)},\qquad \alpha,\beta>0. $$ The maximum-likelihood estimates maximize this function over the positive shape-parameter space.

Worked Example: Writing and Evaluating a Grouped Likelihood

Two independent groups each have four trials. The observed success counts are $y_1=0$ and $y_2=3$. Write their likelihood and evaluate it at $\alpha=\beta=1$.

The two group probabilities multiply, so $$ L(\alpha,\beta)= \left[\binom40\frac{B(\alpha,\beta+4)}{B(\alpha,\beta)}\right] \left[\binom43\frac{B(\alpha+3,\beta+1)}{B(\alpha,\beta)}\right]. $$

At $\alpha=\beta=1$, the beta density is uniform on $(0,1)$. For any $y\in\{0,\ldots,n\}$, the resulting beta-binomial probability is $$ \binom ny B(y+1,n-y+1) =\binom ny\frac{y!(n-y)!}{(n+1)!} =\frac{1}{n+1}. $$ Thus each of these two counts has probability $1/5$, and the joint likelihood is $L(1,1)=(1/5)(1/5)=1/25=0.04$. This is a likelihood value for the observed data at the specified parameter pair, not a probability distribution over parameter values.

Log-Likelihood and Shape-Parameter Scores

For differentiation, write $B(a,b)=\Gamma(a)\Gamma(b)/\Gamma(a+b)$, where $\Gamma$ is the gamma function, and define the digamma function $\psi(x)=d\log\Gamma(x)/dx$ for $x>0$. Taking logarithms turns the product over groups into a sum. The binomial coefficients do not depend on either shape parameter, but they are retained in the log-likelihood because they are part of the likelihood.

\begin{aligned} \ell(\alpha,\beta) =\sum_{i=1}^{m}\Big[& \log\binom{n_i}{y_i} +\log\Gamma(y_i+\alpha) +\log\Gamma(n_i-y_i+\beta)\\ &-\log\Gamma(n_i+\alpha+\beta) +\log\Gamma(\alpha+\beta) -\log\Gamma(\alpha)-\log\Gamma(\beta)\Big]. \end{aligned}

Differentiating each term gives the two score components. The shared $\log\Gamma(\alpha+\beta)$ term contributes to both derivatives, as does the term involving $\log\Gamma(n_i+\alpha+\beta)$.

Score equations: With $\psi$ the digamma function, the likelihood scores are $$ \begin{aligned} U_\alpha(\alpha,\beta) &=\frac{\partial\ell}{\partial\alpha} =\sum_{i=1}^{m} \left[\psi(y_i+\alpha)-\psi(n_i+\alpha+\beta) +\psi(\alpha+\beta)-\psi(\alpha)\right],\\ U_\beta(\alpha,\beta) &=\frac{\partial\ell}{\partial\beta} =\sum_{i=1}^{m} \left[\psi(n_i-y_i+\beta)-\psi(n_i+\alpha+\beta) +\psi(\alpha+\beta)-\psi(\beta)\right]. \end{aligned} $$ An interior maximum must satisfy $U_\alpha=0$ and $U_\beta=0$.

These equations are obtained by differentiating a finite sum of smooth functions. For every $\alpha,\beta>0$, all gamma-function arguments in the log-likelihood are positive, so these derivatives exist. This establishes the score expressions without requiring an interchange of differentiation and integration: the mixture integral has already been evaluated as a beta function.

Worked Example: Checking the Scores at a Specific Parameter Pair

Consider two independent groups, each with $n_i=2$, with observed counts $y_1=0$ and $y_2=2$. Evaluate both score components at $\alpha=\beta=1$.

For positive integers, $\psi(k+1)-\psi(k)=1/k$. In either group, $\psi(n_i+\alpha+\beta)=\psi(4)$ and $\psi(\alpha+\beta)=\psi(2)$, so $\psi(4)-\psi(2)=1/2+1/3=5/6$. For the first group, the $\alpha$-contribution is $\psi(1)-\psi(1)-5/6=-5/6$, and the $\beta$-contribution is $\psi(3)-\psi(1)-5/6=3/2-5/6=2/3$.

For the second group, the $\alpha$-contribution is $\psi(3)-\psi(1)-5/6=2/3$, and the $\beta$-contribution is $\psi(1)-\psi(1)-5/6=-5/6$. Adding group contributions gives $$ U_\alpha(1,1)=-\frac56+\frac23=-\frac16,\qquad U_\beta(1,1)=\frac23-\frac56=-\frac16. $$ The score is not zero, so $(1,1)$ is not an interior stationary point for these observations.

Why Maximum Likelihood Usually Requires Numerical Optimization

The score equations are coupled sums of digamma-function terms. In general they do not simplify to algebraic equations with a closed-form solution for $\alpha$ and $\beta$. An interior maximum-likelihood estimate is therefore usually found numerically by maximizing the log-likelihood or solving the two score equations. Solving the score equations alone is not enough: a stationary point must be checked as a maximum, and the parameter space also has boundaries approached as either shape parameter tends to zero or infinity.

A convenient way to enforce positivity is to optimize in unconstrained coordinates $\eta=\log\alpha$ and $\xi=\log\beta$, then set $\alpha=e^\eta$ and $\beta=e^\xi$. The transformed log-likelihood is $\ell(e^\eta,e^\xi)$, and its gradient follows from the chain rule:

\frac{\partial\ell(e^\eta,e^\xi)}{\partial\eta} =\alpha U_\alpha(\alpha,\beta),\qquad \frac{\partial\ell(e^\eta,e^\xi)}{\partial\xi} =\beta U_\beta(\alpha,\beta).

A numerical routine can use these gradients, with suitable step control, to search for a maximum. In practice, starting from more than one positive parameter pair and comparing attained log-likelihoods helps detect poor local searches. The fitted point should be checked for a small score, a locally negative-definite Hessian when an interior maximum is claimed, and an adequate comparison with relevant boundary limits.

Worked Example: A Symmetric Likelihood and a Constrained Search

Suppose three independent groups each have two trials, with observed counts $0$, $1$, and $2$. To illustrate a one-dimensional constrained search, impose the restriction $\alpha=\beta=t$, where $t>0$. For $n=2$, the beta-binomial probabilities under this restriction are $$ P(Y=0)=P(Y=2)=\frac{t+1}{2(2t+1)},\qquad P(Y=1)=\frac{t}{2t+1}. $$

Therefore the likelihood along the restriction is $$ L(t)=\left(\frac{t+1}{2(2t+1)}\right)^2\left(\frac{t}{2t+1}\right), $$ and its log-derivative is $$ \frac{d}{dt}\log L(t) =2\left(\frac{1}{t+1}-\frac{2}{2t+1}\right) +\frac{1}{t}-\frac{2}{2t+1} =\frac{1-t}{t(t+1)(2t+1)}. $$ The denominator is positive for $t>0$. The derivative is positive when $0<t<1$, zero at $t=1$, and negative when $t>1$, so the constrained likelihood has its unique maximum at $t=1$. Its value is $L(1)=(1/3)^3=1/27$. This simple symmetric example reduces to a scalar equation; with unrestricted $\alpha$ and $\beta$, the two score equations generally need numerical solution. The constrained result does not by itself establish the unrestricted maximum.

What the Likelihood Can and Cannot Identify

The group-level structure is essential: the beta-binomial likelihood represents variation in the latent success probability between groups. Replacing the counts by a single pooled total would discard how successes are distributed across groups, which is information about the mixture. Conversely, if every group has only one trial, the data cannot identify both shape parameters. In that case, $$ P(Y_i=1)=\frac{B(1+\alpha,\beta)}{B(\alpha,\beta)} =\frac{\alpha}{\alpha+\beta},\qquad P(Y_i=0)=\frac{\beta}{\alpha+\beta}. $$ The likelihood then depends only on the ratio $\alpha/(\alpha+\beta)$, not on the two parameters separately. This is an identifiability issue, not a numerical failure.

Another important qualification is that the score equations characterize interior stationary points only. If the likelihood is maximized only as a parameter approaches the edge of its allowed range, an interior score solution need not exist. Numerical fitting should therefore report whether an interior estimate was found and should not treat failure of a solver as proof that the data contain no information.

Common Mistakes

Pooling counts before forming the likelihood. The groups are independent beta-binomial observations, so their probabilities multiply. A pooled total generally does not preserve the likelihood.
Using a binomial likelihood after mixing. A binomial probability conditions on a fixed success probability. Here that probability is random across groups and has been integrated out.
Dropping the shared beta-function denominator. The term $B(\alpha,\beta)$ appears once for every group, so its log contributes to both score sums for every $i$.
Assuming a zero score guarantees the MLE. It identifies an interior stationary point, which must still be checked against other stationary points and boundary behavior.
Expecting both shapes to be estimated from one-trial groups. With $n_i=1$ for all groups, the likelihood contains information only about $\alpha/(\alpha+\beta)$.

Key takeaway: For independent grouped success counts, construct one beta-binomial probability per group, sum their log-probabilities, and differentiate with respect to both positive shape parameters. The resulting digamma score equations generally call for constrained numerical optimization and careful checks of identifiability and boundaries.

Check Your Understanding

Use the model assumptions and likelihood expressions above to answer these questions.

Why does the grouped likelihood multiply one beta-binomial probability for each group?
Which terms in the log-likelihood contribute to the score with respect to $\alpha$?
How can a log transformation enforce positivity when optimizing over $\alpha$ and $\beta$?
Why can groups with one trial each fail to identify both shape parameters?
What additional checks are needed before calling an interior score solution a maximum-likelihood estimate?