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

Code before:

    The test statistic is equal to the estimated proportion, which is simply
    ``3/15``:

    >>> result.statistic
    0.2

    We can use the `proportion_ci()` method of the result to compute the
    confidence interval of the estimate:

    >>> result.proportion_ci(confidence_level=0.95)
    ConfidenceInterval(low=0.05684686759024681, high=1.0)

    """
    k, n, p = np.broadcast_arrays(k, n, p)
    k_valid = (k >= 0) & (k <= n) & (k == np.floor(k))
    n_valid = (n >= 1) & (n == np.floor(n))
    p_valid = (p >= 0) & (p <= 1)
    valid = k_valid & n_valid & p_valid
    k = np.where(valid, k, np.nan)
    n = np.where(valid, n, np.nan)
    p = np.where(valid, p, np.nan)

    if alternative not in ('two-sided', 'less', 'greater'):
        raise ValueError(f"alternative ('{alternative}') not recognized; \n"
                         "must be 'two-sided', 'less' or 'greater'")

    B = _SimpleBinomial(n, p)
    if alternative == 'less':
        pval = B.cdf(k)
    elif alternative == 'greater':
        pval = B.sf(k - 1)
    else:
        # alternative is 'two-sided'
        d = B.pmf(k)
        rerr = 1 + 1e-7

        def k_lt_pn(d, k, p, n):
            B = _SimpleBinomial(n, p)
            ix = _binary_search_for_binom_tst(lambda x1: -B.pmf(x1),
                                              -d*rerr, np.ceil(p * n), n)
            # y is the number of terms between mode and n that are <= d*rerr.
            # ix gave us the first term where a(ix) <= d*rerr < a(ix-1)
            # if the first equality doesn't hold, y=n-ix. Otherwise, we
            # need to include ix as well as the equality holds. Note that
            # the equality will hold in very very rare situations due to rerr.
            y = n - ix + np.asarray(d*rerr == B.pmf(ix), dtype=ix.dtype)
            pval = B.cdf(k) + B.sf(n - y)
            return pval

        def k_gte_pn(d, k, p, n):

Review comment (suggestion):
Unfortunately, `k` and `n` are public attributes of the return object. This makes it unclear whether it's OK to NaN these out when there is, say, a problem with `p`. 

To keep it as simple as possible, I NaN-out all of `k`, `n`, `p` wherever there is a problem with one of them because the NaNs propagate naturally through the calculations. We can do whatever we want; it just might require more effort.

Code after:
    0.18406106910639114

    The null hypothesis cannot be rejected at the 5% level of significance
    because the returned p-value is greater than the critical value of 5%.

    The test statistic is equal to the estimated proportion, which is simply
    ``3/15``:

    >>> result.statistic
    0.2

    We can use the `proportion_ci()` method of the result to compute the
    confidence interval of the estimate:

    >>> result.proportion_ci(confidence_level=0.95)
    ConfidenceInterval(low=0.05684686759024681, high=1.0)

    """
    k, n, p = np.broadcast_arrays(k, n, p)
    k_valid = (k >= 0) & (k <= n) & (k == np.floor(k))
    n_valid = (n >= 1) & (n == np.floor(n))
    p_valid = (p >= 0) & (p <= 1)
    valid = k_valid & n_valid & p_valid
    k = np.where(valid, k, np.nan)
    n = np.where(valid, n, np.nan)
    p = np.where(valid, p, np.nan)

    if alternative not in ('two-sided', 'less', 'greater'):
        raise ValueError(f"alternative ('{alternative}') not recognized; \n"
                         "must be 'two-sided', 'less' or 'greater'")

    B = _SimpleBinomial(n, p)
    if alternative == 'less':
        pval = B.cdf(k)
    elif alternative == 'greater':
        pval = B.sf(k - 1)
    else:
        # alternative is 'two-sided'
        d = B.pmf(k)
        rerr = 1 + 1e-7

        def k_lt_pn(d, k, p, n):
            B = _SimpleBinomial(n, p)
            ix = _binary_search_for_binom_tst(lambda x1: -B.pmf(x1),
                                              -d*rerr, np.ceil(p * n), n)
            # y is the number of terms between mode and n that are <= d*rerr.
            # ix gave us the first term where a(ix) <= d*rerr < a(ix-1)
            # if the first equality doesn't hold, y=n-ix. Otherwise, we
            # need to include ix as well as the equality holds. Note that
            # the equality will hold in very very rare situations due to rerr.
            y = n - ix + np.asarray(d*rerr == B.pmf(ix), dtype=ix.dtype)