Overflow during the estimation of empirical accuracy

[fc-jian]

Issue:
At line 268 of alignparse/consensus.py, the binomial coefficients are precomputed for subsequent likelihood calculations. However, for large values of n and k, directly computing them using scipy.special.binom results in inf values, which can cause errors during likelihood optimization.

Proposed solution:
To address this, we can estimate the logarithm of the binomial coefficients using Stirling’s approximation instead of computing them directly. The subsequent likelihood calculation can then be divided into two cases:

• When n = u: the binomial coefficients remain small, so direct computation can be retained.
• When n ≠ u: the binomial coefficients can become extremely large, but since the additional p ** n term is absent, the computation can be performed entirely in log space.

Example implementation:

def log_comb(n, k):
    """
    Robustly compute log[C(n, k)].
    For small/moderate n, use scipy.special.binom.
    For large n, use Stirling approximation.
    """
    try:
        # Try direct computation
        val = scipy.special.binom(n, k)
        # if finite and not huge
        if numpy.isfinite(val) and val > 0 and val < numpy.inf:
            return numpy.log(val)
        else:
            raise OverflowError
    except Exception:
        # Stirling approximation for large n,k
        # Use: log C(n,k) = n*log(n) - k*log(k) - (n-k)*log(n-k)
        #       + 0.5*log(n / (2π k (n-k)))
        log2pi = numpy.log(2 * numpy.pi)
        term = (
            n * numpy.log(n)
            - k * numpy.log(k)
            - (n - k) * numpy.log(n - k)
            + 0.5 * (numpy.log(n) - numpy.log(k) - numpy.log(n - k) - log2pi)
        )
        return term

class _LnL_error_rate:
    def __init__(self, df, *, n_col, u_col, count_col):
        """See main class docstring."""
        _df = df.assign(
            n=lambda x: x[n_col],
            u=lambda x: x[u_col],
            count=lambda x: x[count_col],
        )
        self._df_un = _df[['n', 'u', 'count']].query("n == u").assign(
            binom=lambda x: scipy.special.binom(x["n"], x["u"] - 1)
        )
        self._df_other = _df[['n', 'u', 'count']].query("n != u").assign(
            log_binom=lambda x: [log_comb(n, k) for n, k in zip(x["n"], x["u"] - 1)]
        )

    def lnlik(self, eps):
        """Log likelihood for error rate `eps`."""
        return sum(
            self._df_other["count"]
            * (
                self._df_other["log_binom"]
                + (self._df_other["n"] - self._df_other["u"] + 1) * numpy.log(1 - eps)
                + (self._df_other["u"] - 1) * numpy.log(eps)
            )
        ) + sum(
            self._df_un["count"]
            * numpy.log(
                self._df_un["binom"]
                * (1 - eps) ** (self._df_un["n"] - self._df_un["u"] + 1)
                * eps ** (self._df_un["u"] - 1)
                + eps ** self._df_un["n"]
            )
        )

    def neg_lnlik(self, eps):
        """Negative log likelihood for error rate `epsilon`."""
        return -self.lnlik(eps)

    def maxlik_eps(self):
        """Maximum likelihood value of error rate `epsilon`."""
        res = scipy.optimize.minimize_scalar(
            self.neg_lnlik, bounds=(1e-8, 1 - 1e-8), method="bounded"
        )
        if not res.success:
            raise RuntimeError(f"optimization failed:\n{res}")
        return res.x

Would you be open to this approach? If it seems reasonable, I’d be happy to open a PR to address this issue.

[jbloom]

This looks correct to me, and it would be great if you are interested in making a pull request @fc-jian. I can review more carefully after that, but the only suggestion I have right now is that in log_comb it is generally not a good idea to have an except with Exception as that will catch any error, and we want it to be more specific (https://jerrynsh.com/stop-using-exceptions-like-this-in-python/).

But I'm not clear why you need a try / except at all, could you not just move the part in the except block to the else block and get rid of the try / except altogether? Or does val = scipy.special.binom(n, k) sometimes raise an error itself?

[jbloom]

This issue is being closed by #109.

The key point is that the computational of the empirical accuracy (consensus.empirical_accuracy) ran into numerical issues if the numbers were large as it involved computing very large numbers and then taking their logs. Above in this issue, @fc-jian originally proposed using Stirling's approximation, and made a draft pull request #108 to fix that.

However, in looking more I discovered that the built-in python gammaln function is even a better way to do this. I also updated the docs to describe the new math being done, which is also pasted below:

This pull request therefore solves #106 and is in lieu of #108, as I think it is a better solution.

@fc-jian, thanks so much for noting and flagging all of this!