Stable Bernoulli and Binomial Log-Likelihoods

Why Stable Likelihood Evaluation Matters

For Bernoulli and binomial data, the likelihood is simple to write down, but a direct computer evaluation can still fail. The usual expression contains products such as \(0\log(0)\), which are undefined as ordinary arithmetic, and terms such as \(\log(1-p)\), which can lose accuracy when \(p\) is close to one. These are not merely cosmetic issues: a likelihood calculation that returns a NaN or a substantially inaccurate value can disrupt an optimizer's comparisons between candidate parameter values.

This tutorial starts with probabilities \(p\) already computed. Its goal is to evaluate the Bernoulli and grouped-binomial log-likelihood contributions safely, including when a count is zero or a probability is exactly zero or one. It does not change the statistical model; it changes how its log-likelihood is represented and evaluated.

Bernoulli and Grouped-Binomial Contributions

Let \(Y\in\{0,1\}\) be a Bernoulli observation with success probability \(p\). Its probability mass is \(p^Y(1-p)^{1-Y}\), so its log-likelihood contribution is \(\operatorname{xlogy}(Y,p)+\operatorname{xlogy}(1-Y,1-p)\). For a grouped observation, let \(m\) be the number of trials and \(y\) the number of successes, where \(m\) is a nonnegative integer and \(y\in\{0,\ldots,m\}\). The binomial probability mass is \(\binom{m}{y}p^y(1-p)^{m-y}\), giving the log-likelihood contribution \[ \log\binom{m}{y}+\operatorname{xlogy}(y,p)+\operatorname{xlogy}(m-y,1-p). \]

The binomial coefficient is part of the full probability mass. It does not depend on \(p\), so it can be dropped when optimizing \(p\) for fixed \(m\) and \(y\). Keep it when the goal is the full log probability or a likelihood comparison for which that constant matters. For ordinary count sizes, the coefficient's logarithm can be computed using log-gamma functions; for very large counts, a specialized log-binomial-coefficient routine may be preferable.

Here \(\log1p(z)\) denotes the function that evaluates \(\log(1+z)\) accurately when \(z\) is near zero. Thus \(\log1p(-p)\) is useful when \(p\) is close to zero as well as when \(p\) is close to one; in particular, it avoids first rounding \(1-p\) and then taking its logarithm. For the success term, use \(\log(p)\) when \(p\) is not close to one. When \(p\) is close to one, \(\log1p(p-1)\) is an alternative that can preserve accuracy in the small value of \(\log(p)\). The expression \(p-1\) is exact in ordinary floating-point arithmetic for representable \(p\in[1/2,1]\), by Sterbenz's lemma. A robust implementation must also respect the domain \(p\in[0,1]\); a value outside that interval is an invalid probability, not a value to repair silently.

Worked Evaluations

Why the Stable Formula Is the Correct Log-Likelihood

Proof. If \(0<p<1\), both category probabilities are positive. Applying the logarithm to the binomial mass and using the logarithm product rule gives \[ \log\left\{\binom{m}{y}p^y(1-p)^{m-y}\right\} =\log\binom{m}{y}+y\log(p)+(m-y)\log(1-p). \] The xlogy definition gives exactly these two count-weighted terms, including a zero term when either count is zero.

At \(p=0\), if \(y=0\), the binomial mass is \(\binom{m}{0}1^m=1\), whose log is zero; the formula gives \(\log\binom{m}{0}=0\), a zero success term, and \(m\log(1)=0\). If \(y>0\), the mass is zero because it includes \(0^y\); the formula assigns \(-\infty\) to the positive success count paired with probability zero. At \(p=1\), if \(y=m\), the mass is one and the formula gives zero. If \(y<m\), the mass is zero and the positive failure count paired with probability zero gives \(-\infty\). These cases include \(m=0,y=0\), for which the mass and log-likelihood are both one and zero, respectively. Thus the expression agrees with the log mass throughout the stated domain. \(\square\)

Likelihood Maximization and Boundary Counts

The same representation makes the boundary cases explicit. The next result describes the maximum over \(p\) for a single grouped-binomial observation. It is useful as a basic check on computed likelihoods, not as a substitute for evaluating them stably.

Proof. The term \(\log\binom{m}{y}\) is constant in \(p\). First suppose \(0<y<m\). For \(0<p<1\), differentiating the remaining terms gives \[ \ell'(p)=\frac{y}{p}-\frac{m-y}{1-p} =\frac{y-mp}{p(1-p)}. \] The denominator is positive. Therefore \(\ell'(p)>0\) when \(p<y/m\), equals zero at \(p=y/m\), and is negative when \(p>y/m\). The log-likelihood increases up to \(y/m\) and decreases after it. The endpoint log-likelihoods are \(-\infty\), since each endpoint assigns zero probability to an outcome with both successes and failures. Hence \(y/m\) is the unique maximizer.

If \(y=0\), for \(0\leq p<1\) the likelihood-dependent part is \(m\log(1-p)\), which is maximized at \(p=0\), where it equals zero; at \(p=1\) it is \(-\infty\). If \(y=m\), the corresponding part is \(m\log(p)\), maximized at \(p=1\); at \(p=0\) it is \(-\infty\). This proves all cases. \(\square\)

Common Mistake: Computing Every Term Literally

A tempting implementation is to calculate \(y\log(p)+(m-y)\log(1-p)\) unconditionally. At \(y=0,p=0\), for example, the first product is \(0\cdot(-\infty)\). Real-number algebra says the count contributes zero, but floating-point arithmetic commonly returns NaN for this product. A different mistake is to compute \(1-p\), then take its logarithm, and assume the result is reliable for every probability near one.

The correction is to branch on the counts before computing their log terms: skip a term whose count is zero; for a positive failure count use \(\log1p(-p)\); and for a positive success count use a stable evaluation of \(\log(p)\). At exact boundaries, a positive count paired with probability zero must return \(-\infty\). If the probability has already rounded to zero or one, no later log-likelihood routine can recover information lost when it was formed; the stable evaluation can only work with the supplied value.