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

Code before:

    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)
        return mode[()]

    def _moment_raw_formula(self, order, *, n, p, **kwargs):

Review comment (suggestion):
```suggestion
            lambda *args: np.log1p(-scu._binom_sf(*args))
```

Code after:

    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)
        return mode[()]

    def _moment_raw_formula(self, order, *, n, p, **kwargs):