Repository: scipy/scipy
File: scipy/stats/_new_distributions.py (Python)

Code before:

    _parameterizations = [_Parameterization(_n_param, _p_param)]
    _variable = _x_param

    def __init__(self, *, n, p, **kwargs):
        super().__init__(n=n, p=p, **kwargs)

    def _pmf_formula(self, x, *, n, p, **kwargs):
        return scu._binom_pmf(x, n, p)

    def _logpmf_formula(self, x, *, n, p, **kwargs):
        # This implementation is from the ``scipy.stats.binom`` and could be improved
        # by using a more numerically sound implementation of the absolute value of
        # the binomial coefficient.
        combiln = (
            special.gammaln(n+1) - (special.gammaln(x+1) + special.gammaln(n-x+1))
        )
        return combiln + special.xlogy(x, p) + special.xlog1py(n-x, -p)

    def _cdf_formula(self, x, *, n, p, **kwargs):
        return scu._binom_cdf(x, n, p)

    def _logcdf_formula(self, x, *, n, p, **kwargs):
        # todo: add this strategy to infrastructure more generally while avoiding
        # expensive calculations to figure out which tail we're in
        mode = self._mode_formula(n=n, p=p)
        return xpx.apply_where(x < mode, (x, n, p),
            lambda *args: np.log(scu._binom_cdf(*args)),
            lambda *args: scu._log1mexp(np.log(scu._binom_sf(*args)))
        )

    def _ccdf_formula(self, x, *, n, p, **kwargs):
        return scu._binom_sf(x, n, p)

    def _logccdf_formula(self, x, *, n, p, **kwargs):
        mode = self._mode_formula(n=n, p=p)
        return xpx.apply_where(x > mode, (x, n, p),
            lambda *args: np.log(scu._binom_sf(*args)),
            lambda *args: scu._log1mexp(np.log(scu._binom_cdf(*args)))
        )

    def _icdf_formula(self, x, *, n, p, **kwargs):
        return scu._binom_ppf(x, n, p)

    def _iccdf_formula(self, x, *, n, p, **kwargs):
        return scu._binom_isf(x, n, p)

    def _mode_formula(self, *, n, p, **kwargs):
        # https://en.wikipedia.org/wiki/Binomial_distribution#Mode
        mode = np.floor((n+1)*p)
        mode = np.where(p == 1, mode - 1, mode)

Review comment (security):
Here, I've used the mode to determine which tail to use ~~because it is faster than other methods~~(not true, see https://github.com/scipy/scipy/pull/23716#discussion_r2395250383 [^1]. 

Ideally, we'd have this sort of logic in the infrastructure itself rather than in individual distributions. We could *always* make the decision based on the median, but sometimes the median is very expensive to calculate. So I think I'll allow distribution authors to specify a method that identifies the tail in which the CDF/CCDF should be evaluated. Something similar could be useful to determine where the probability mass is concentrated. That will help the integrator out when there are large, finite integration intervals where the integrand is essentially zero except for within a small sub-interval.

[^1]: I should think a bit about whether this can be off by one compared to the optimal choice for extreme parameters, but this is still an improvement over status quo.

Code after:

    _parameterizations = [_Parameterization(_n_param, _p_param)]
    _variable = _x_param

    def __init__(self, *, n, p, **kwargs):
        super().__init__(n=n, p=p, **kwargs)

    def _pmf_formula(self, x, *, n, p, **kwargs):
        return scu._binom_pmf(x, n, p)

    def _logpmf_formula(self, x, *, n, p, **kwargs):
        # This implementation is from the ``scipy.stats.binom`` and could be improved
        # by using a more numerically sound implementation of the absolute value of
        # the binomial coefficient.
        combiln = (
            special.gammaln(n+1) - (special.gammaln(x+1) + special.gammaln(n-x+1))
        )
        return combiln + special.xlogy(x, p) + special.xlog1py(n-x, -p)

    def _cdf_formula(self, x, *, n, p, **kwargs):
        return scu._binom_cdf(x, n, p)

    def _logcdf_formula(self, x, *, n, p, **kwargs):
        # todo: add this strategy to infrastructure more generally, but allow dist
        #   author to specify threshold other than median in case median is expensive
        median = self._icdf_formula(0.5, n=n, p=p)
        return xpx.apply_where(x < median, (x, n, p),
            lambda *args: np.log(scu._binom_cdf(*args)),
            lambda *args: np.log1p(-scu._binom_sf(*args))
        )

    def _ccdf_formula(self, x, *, n, p, **kwargs):
        return scu._binom_sf(x, n, p)

    def _logccdf_formula(self, x, *, n, p, **kwargs):
        median = self._icdf_formula(0.5, n=n, p=p)
        return xpx.apply_where(x < median, (x, n, p),
            lambda *args: np.log1p(-scu._binom_cdf(*args)),
            lambda *args: np.log(scu._binom_sf(*args))
        )

    def _icdf_formula(self, x, *, n, p, **kwargs):
        return scu._binom_ppf(x, n, p)

    def _iccdf_formula(self, x, *, n, p, **kwargs):
        return scu._binom_isf(x, n, p)

    def _mode_formula(self, *, n, p, **kwargs):
        # https://en.wikipedia.org/wiki/Binomial_distribution#Mode
        mode = np.floor((n+1)*p)
        mode = np.where(p == 1, mode - 1, mode)