"""
Tests and Confidence Intervals for Binomial Proportions

Created on Fri Mar 01 00:23:07 2013

Author: Josef Perktold
License: BSD-3
"""

from statsmodels.compat.python import lzip

from typing import Callable

import numpy as np
import pandas as pd
from scipy import optimize, stats

from statsmodels.stats.base import AllPairsResults, HolderTuple
from statsmodels.stats.weightstats import _zstat_generic2
from statsmodels.tools.sm_exceptions import HypothesisTestWarning
from statsmodels.tools.testing import Holder
from statsmodels.tools.validation import array_like

FLOAT_INFO = np.finfo(float)


def _bound_proportion_confint(
    func: Callable[[float], float], qi: float, lower: bool = True
) -> float:
    """
    Try hard to find a bound different from eps/1 - eps in proportion_confint

    Parameters
    ----------
    func : callable
        Callable function to use as the objective of the search
    qi : float
        The empirical success rate
    lower : bool
        Whether to fund a lower bound for the left side of the CI

    Returns
    -------
    float
        The coarse bound
    """
    default = FLOAT_INFO.eps if lower else 1.0 - FLOAT_INFO.eps

    def step(v):
        return v / 8 if lower else v + (1.0 - v) / 8

    x = step(qi)
    w = func(x)
    cnt = 1
    while w > 0 and cnt < 10:
        x = step(x)
        w = func(x)
        cnt += 1
    return x if cnt < 10 else default


def _bisection_search_conservative(
    func: Callable[[float], float], lb: float, ub: float, steps: int = 27
) -> tuple[float, float]:
    """
    Private function used as a fallback by proportion_confint

    Used when brentq returns a non-conservative bound for the CI

    Parameters
    ----------
    func : callable
        Callable function to use as the objective of the search
    lb : float
        Lower bound
    ub : float
        Upper bound
    steps : int
        Number of steps to use in the bisection

    Returns
    -------
    est : float
        The estimated value.  Will always produce a negative value of func
    func_val : float
        The value of the function at the estimate
    """
    upper = func(ub)
    lower = func(lb)
    best = upper if upper < 0 else lower
    best_pt = ub if upper < 0 else lb
    if np.sign(lower) == np.sign(upper):
        raise ValueError("problem with signs")
    mp = (ub + lb) / 2
    mid = func(mp)
    if (mid < 0) and (mid > best):
        best = mid
        best_pt = mp
    for _ in range(steps):
        if np.sign(mid) == np.sign(upper):
            ub = mp
            upper = mid
        else:
            lb = mp
        mp = (ub + lb) / 2
        mid = func(mp)
        if (mid < 0) and (mid > best):
            best = mid
            best_pt = mp
    return best_pt, best


def proportion_confint(
    count, nobs, alpha: float = 0.05, method="normal", alternative: str = "two-sided"
):
    """
    Confidence interval for a binomial proportion

    Parameters
    ----------
    count : {int or float, array_like}
        number of successes, can be pandas Series or DataFrame. Arrays
        must contain integer values if method is "binom_test".
    nobs : {int or float, array_like}
        total number of trials.  Arrays must contain integer values if method
        is "binom_test".
    alpha : float
        Significance level, default 0.05. Must be in (0, 1)
    method : {"normal", "agresti_coull", "beta", "wilson", "binom_test"}
        default: "normal"
        method to use for confidence interval. Supported methods:

         - `normal` : asymptotic normal approximation
         - `agresti_coull` : Agresti-Coull interval
         - `beta` : Clopper-Pearson interval based on Beta distribution
         - `wilson` : Wilson Score interval
         - `jeffreys` : Jeffreys Bayesian Interval
         - `binom_test` : Numerical inversion of binom_test

    alternative : {"two-sided", "larger", "smaller"}
        default: "two-sided"
        specifies whether to calculate a two-sided or one-sided confidence interval.

    Returns
    -------
    ci_low, ci_upp : {float, ndarray, Series DataFrame}
        larger and smaller confidence level with coverage (approximately) 1-alpha.
        When a pandas object is returned, then the index is taken from `count`.
        When side is not "two-sided", lower or upper bound is set to 0 or 1 respectively.

    Notes
    -----
    Beta, the Clopper-Pearson exact interval has coverage at least 1-alpha,
    but is in general conservative. Most of the other methods have average
    coverage equal to 1-alpha, but will have smaller coverage in some cases.

    The "beta" and "jeffreys" interval are central, they use alpha/2 in each
    tail, and alpha is not adjusted at the boundaries. In the extreme case
    when `count` is zero or equal to `nobs`, then the coverage will be only
    1 - alpha/2 in the case of "beta".

    The confidence intervals are clipped to be in the [0, 1] interval in the
    case of "normal" and "agresti_coull".

    Method "binom_test" directly inverts the binomial test in scipy.stats.
    which has discrete steps.

    TODO: binom_test intervals raise an exception in small samples if one
       interval bound is close to zero or one.

    References
    ----------
    .. [*] https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval

    .. [*] Brown, Lawrence D.; Cai, T. Tony; DasGupta, Anirban (2001).
       "Interval Estimation for a Binomial Proportion", Statistical
       Science 16 (2): 101-133. doi:10.1214/ss/1009213286.
    """
    is_scalar = np.isscalar(count) and np.isscalar(nobs)
    is_pandas = isinstance(count, (pd.Series, pd.DataFrame))
    count_a = array_like(count, "count", optional=False, ndim=None)
    nobs_a = array_like(nobs, "nobs", optional=False, ndim=None)

    def _check(x: np.ndarray, name: str) -> np.ndarray:
        if np.issubdtype(x.dtype, np.integer):
            return x
        y = x.astype(np.int64, casting="unsafe")
        if np.any(y != x):
            raise ValueError(
                f"{name} must have an integral dtype. Found data with "
                f"dtype {x.dtype}"
            )
        return y

    if method == "binom_test":
        count_a = _check(np.asarray(count_a), "count")
        nobs_a = _check(np.asarray(nobs_a), "count")

    q_ = count_a / nobs_a

    if alternative == "two-sided":
        if method != "binom_test":
            alpha = alpha / 2.0
    elif alternative not in ["larger", "smaller"]:
        raise NotImplementedError(f"alternative {alternative} is not available")

    if method == "normal":
        std_ = np.sqrt(q_ * (1 - q_) / nobs_a)
        dist = stats.norm.isf(alpha) * std_
        ci_low = q_ - dist
        ci_upp = q_ + dist
    elif method == "binom_test" and alternative == "two-sided":

        def func_factory(count: int, nobs: int) -> Callable[[float], float]:
            if hasattr(stats, "binomtest"):

                def func(qi):
                    return stats.binomtest(count, nobs, p=qi).pvalue - alpha

            else:
                # Remove after min SciPy >= 1.7
                def func(qi):
                    return stats.binom_test(count, nobs, p=qi) - alpha

            return func

        bcast = np.broadcast(count_a, nobs_a)
        ci_low = np.zeros(bcast.shape)
        ci_upp = np.zeros(bcast.shape)
        index = bcast.index
        for c, n in bcast:
            # Enforce symmetry
            reverse = False
            _q = q_.flat[index]
            if c > n // 2:
                c = n - c
                reverse = True
                _q = 1 - _q
            func = func_factory(c, n)
            if c == 0:
                ci_low.flat[index] = 0.0
            else:
                lower_bnd = _bound_proportion_confint(func, _q, lower=True)
                val, _z = optimize.brentq(func, lower_bnd, _q, full_output=True)
                if func(val) > 0:
                    power = 10
                    new_lb = val - (val - lower_bnd) / 2**power
                    while func(new_lb) > 0 and power >= 0:
                        power -= 1
                        new_lb = val - (val - lower_bnd) / 2**power
                    val, _ = _bisection_search_conservative(func, new_lb, _q)
                ci_low.flat[index] = val
            if c == n:
                ci_upp.flat[index] = 1.0
            else:
                upper_bnd = _bound_proportion_confint(func, _q, lower=False)
                val, _z = optimize.brentq(func, _q, upper_bnd, full_output=True)
                if func(val) > 0:
                    power = 10
                    new_ub = val + (upper_bnd - val) / 2**power
                    while func(new_ub) > 0 and power >= 0:
                        power -= 1
                        new_ub = val - (upper_bnd - val) / 2**power
                    val, _ = _bisection_search_conservative(func, _q, new_ub)
                ci_upp.flat[index] = val
            if reverse:
                temp = ci_upp.flat[index]
                ci_upp.flat[index] = 1 - ci_low.flat[index]
                ci_low.flat[index] = 1 - temp
            index = bcast.index
    elif method == "beta" or (method == "binom_test" and alternative != "two-sided"):
        ci_low = stats.beta.ppf(alpha, count_a, nobs_a - count_a + 1)
        ci_upp = stats.beta.isf(alpha, count_a + 1, nobs_a - count_a)

        if np.ndim(ci_low) > 0:
            ci_low.flat[q_.flat == 0] = 0
            ci_upp.flat[q_.flat == 1] = 1
        else:
            ci_low = 0 if q_ == 0 else ci_low
            ci_upp = 1 if q_ == 1 else ci_upp
    elif method == "agresti_coull":
        crit = stats.norm.isf(alpha)
        nobs_c = nobs_a + crit**2
        q_c = (count_a + crit**2 / 2.0) / nobs_c
        std_c = np.sqrt(q_c * (1.0 - q_c) / nobs_c)
        dist = crit * std_c
        ci_low = q_c - dist
        ci_upp = q_c + dist
    elif method == "wilson":
        crit = stats.norm.isf(alpha)
        crit2 = crit**2
        denom = 1 + crit2 / nobs_a
        center = (q_ + crit2 / (2 * nobs_a)) / denom
        dist = crit * np.sqrt(q_ * (1.0 - q_) / nobs_a + crit2 / (4.0 * nobs_a**2))
        dist /= denom
        ci_low = center - dist
        ci_upp = center + dist
    # method adjusted to be more forgiving of misspellings or incorrect option name
    elif method[:4] == "jeff":
        ci_low = stats.beta.ppf(alpha, count_a + 0.5, nobs_a - count_a + 0.5)
        ci_upp = stats.beta.isf(alpha, count_a + 0.5, nobs_a - count_a + 0.5)
    else:
        raise NotImplementedError(f"method {method} is not available")
    if method in ["normal", "agresti_coull"]:
        ci_low = np.clip(ci_low, 0, 1)
        ci_upp = np.clip(ci_upp, 0, 1)
    if is_pandas:
        container = pd.Series if isinstance(count, pd.Series) else pd.DataFrame
        ci_low = container(ci_low, index=count.index)
        ci_upp = container(ci_upp, index=count.index)
    if alternative == "larger":
        ci_low = 0
    elif alternative == "smaller":
        ci_upp = 1
    if is_scalar:
        return float(ci_low), float(ci_upp)
    return ci_low, ci_upp


def multinomial_proportions_confint(counts, alpha=0.05, method="goodman"):
    """
    Confidence intervals for multinomial proportions.

    Parameters
    ----------
    counts : array_like of int, 1-D
        Number of observations in each category.
    alpha : float in (0, 1), optional
        Significance level, defaults to 0.05.
    method : {'goodman', 'sison-glaz'}, optional
        Method to use to compute the confidence intervals; available methods
        are:

         - `goodman`: based on a chi-squared approximation, valid if all
           values in `counts` are greater or equal to 5 [2]_
         - `sison-glaz`: less conservative than `goodman`, but only valid if
           `counts` has 7 or more categories (``len(counts) >= 7``) [3]_

    Returns
    -------
    confint : ndarray, 2-D
        Array of [lower, upper] confidence levels for each category, such that
        overall coverage is (approximately) `1-alpha`.

    Raises
    ------
    ValueError
        If `alpha` is not in `(0, 1)` (bounds excluded), or if the values in
        `counts` are not all positive or null.
    NotImplementedError
        If `method` is not kown.
    Exception
        When ``method == 'sison-glaz'``, if for some reason `c` cannot be
        computed; this signals a bug and should be reported.

    Notes
    -----
    The `goodman` method [2]_ is based on approximating a statistic based on
    the multinomial as a chi-squared random variable. The usual recommendation
    is that this is valid if all the values in `counts` are greater than or
    equal to 5. There is no condition on the number of categories for this
    method.

    The `sison-glaz` method [3]_ approximates the multinomial probabilities,
    and evaluates that with a maximum-likelihood estimator. The first
    approximation is an Edgeworth expansion that converges when the number of
    categories goes to infinity, and the maximum-likelihood estimator converges
    when the number of observations (``sum(counts)``) goes to infinity. In
    their paper, Sison & Glaz demo their method with at least 7 categories, so
    ``len(counts) >= 7`` with all values in `counts` at or above 5 can be used
    as a rule of thumb for the validity of this method. This method is less
    conservative than the `goodman` method (i.e. it will yield confidence
    intervals closer to the desired significance level), but produces
    confidence intervals of uniform width over all categories (except when the
    intervals reach 0 or 1, in which case they are truncated), which makes it
    most useful when proportions are of similar magnitude.

    Aside from the original sources ([1]_, [2]_, and [3]_), the implementation
    uses the formulas (though not the code) presented in [4]_ and [5]_.

    References
    ----------
    .. [1] Levin, Bruce, "A representation for multinomial cumulative
           distribution functions," The Annals of Statistics, Vol. 9, No. 5,
           1981, pp. 1123-1126.

    .. [2] Goodman, L.A., "On simultaneous confidence intervals for multinomial
           proportions," Technometrics, Vol. 7, No. 2, 1965, pp. 247-254.

    .. [3] Sison, Cristina P., and Joseph Glaz, "Simultaneous Confidence
           Intervals and Sample Size Determination for Multinomial
           Proportions," Journal of the American Statistical Association,
           Vol. 90, No. 429, 1995, pp. 366-369.

    .. [4] May, Warren L., and William D. Johnson, "A SAS® macro for
           constructing simultaneous confidence intervals  for multinomial
           proportions," Computer methods and programs in Biomedicine, Vol. 53,
           No. 3, 1997, pp. 153-162.

    .. [5] May, Warren L., and William D. Johnson, "Constructing two-sided
           simultaneous confidence intervals for multinomial proportions for
           small counts in a large number of cells," Journal of Statistical
           Software, Vol. 5, No. 6, 2000, pp. 1-24.
    """
    if alpha <= 0 or alpha >= 1:
        raise ValueError("alpha must be in (0, 1), bounds excluded")
    counts = np.array(counts, dtype=float)
    if (counts < 0).any():
        raise ValueError("counts must be >= 0")

    n = counts.sum()
    k = len(counts)
    proportions = counts / n
    if method == "goodman":
        chi2 = stats.chi2.ppf(1 - alpha / k, 1)
        delta = chi2**2 + (4 * n * proportions * chi2 * (1 - proportions))
        region = (
            (2 * n * proportions + chi2 + np.array([-np.sqrt(delta), np.sqrt(delta)]))
            / (2 * (chi2 + n))
        ).T
    elif method[:5] == "sison":  # We accept any name starting with 'sison'
        # Define a few functions we'll use a lot.
        def poisson_interval(interval, p):
            """
            Compute P(b <= Z <= a) where Z ~ Poisson(p) and
            `interval = (b, a)`.
            """
            b, a = interval
            prob = stats.poisson.cdf(a, p) - stats.poisson.cdf(b - 1, p)
            return prob

        def truncated_poisson_factorial_moment(interval, r, p):
            """
            Compute mu_r, the r-th factorial moment of a poisson random
            variable of parameter `p` truncated to `interval = (b, a)`.
            """
            b, a = interval
            return p**r * (
                1
                - (
                    (
                        poisson_interval((a - r + 1, a), p)
                        - poisson_interval((b - r, b - 1), p)
                    )
                    / poisson_interval((b, a), p)
                )
            )

        def edgeworth(intervals):
            """
            Compute the Edgeworth expansion term of Sison & Glaz's formula
            (1) (approximated probability for multinomial proportions in a
            given box).
            """
            # Compute means and central moments of the truncated poisson
            # variables.
            mu_r1, mu_r2, mu_r3, mu_r4 = (
                np.array(
                    [
                        truncated_poisson_factorial_moment(interval, r, p)
                        for (interval, p) in zip(intervals, counts)
                    ]
                )
                for r in range(1, 5)
            )
            mu = mu_r1
            mu2 = mu_r2 + mu - mu**2
            mu3 = mu_r3 + mu_r2 * (3 - 3 * mu) + mu - 3 * mu**2 + 2 * mu**3
            mu4 = (
                mu_r4
                + mu_r3 * (6 - 4 * mu)
                + mu_r2 * (7 - 12 * mu + 6 * mu**2)
                + mu
                - 4 * mu**2
                + 6 * mu**3
                - 3 * mu**4
            )

            # Compute expansion factors, gamma_1 and gamma_2.
            g1 = mu3.sum() / mu2.sum() ** 1.5
            g2 = (mu4.sum() - 3 * (mu2**2).sum()) / mu2.sum() ** 2

            # Compute the expansion itself.
            x = (n - mu.sum()) / np.sqrt(mu2.sum())
            phi = np.exp(-(x**2) / 2) / np.sqrt(2 * np.pi)
            H3 = x**3 - 3 * x
            H4 = x**4 - 6 * x**2 + 3
            H6 = x**6 - 15 * x**4 + 45 * x**2 - 15
            f = phi * (1 + g1 * H3 / 6 + g2 * H4 / 24 + g1**2 * H6 / 72)
            return f / np.sqrt(mu2.sum())

        def approximated_multinomial_interval(intervals):
            """
            Compute approximated probability for Multinomial(n, proportions)
            to be in `intervals` (Sison & Glaz's formula (1)).
            """
            return np.exp(
                np.sum(
                    np.log(
                        [
                            poisson_interval(interval, p)
                            for (interval, p) in zip(intervals, counts)
                        ]
                    )
                )
                + np.log(edgeworth(intervals))
                - np.log(stats.poisson._pmf(n, n))
            )

        def nu(c):
            """
            Compute interval coverage for a given `c` (Sison & Glaz's
            formula (7)).
            """
            return approximated_multinomial_interval(
                [
                    (np.maximum(count - c, 0), np.minimum(count + c, n))
                    for count in counts
                ]
            )

        # Find the value of `c` that will give us the confidence intervals
        # (solving nu(c) <= 1 - alpha < nu(c + 1).
        c = 1.0
        nuc = nu(c)
        nucp1 = nu(c + 1)
        while not (nuc <= (1 - alpha) < nucp1):
            if c > n:
                raise Exception(
                    "Couldn't find a value for `c` that "
                    "solves nu(c) <= 1 - alpha < nu(c + 1)"
                )
            c += 1
            nuc = nucp1
            nucp1 = nu(c + 1)

        # Compute gamma and the corresponding confidence intervals.
        g = (1 - alpha - nuc) / (nucp1 - nuc)
        ci_lower = np.maximum(proportions - c / n, 0)
        ci_upper = np.minimum(proportions + (c + 2 * g) / n, 1)
        region = np.array([ci_lower, ci_upper]).T
    else:
        raise NotImplementedError('method "%s" is not available' % method)
    return region


def samplesize_confint_proportion(proportion, half_length, alpha=0.05, method="normal"):
    """
    Find sample size to get desired confidence interval length

    Parameters
    ----------
    proportion : float in (0, 1)
        proportion or quantile
    half_length : float in (0, 1)
        desired half length of the confidence interval
    alpha : float in (0, 1)
        significance level, default 0.05,
        coverage of the two-sided interval is (approximately) ``1 - alpha``
    method : str in ['normal']
        method to use for confidence interval,
        currently only normal approximation

    Returns
    -------
    n : float
        sample size to get the desired half length of the confidence interval

    Notes
    -----
    this is mainly to store the formula.
    possible application: number of replications in bootstrap samples

    """
    q_ = proportion
    if method == "normal":
        n = q_ * (1 - q_) / (half_length / stats.norm.isf(alpha / 2.0)) ** 2
    else:
        raise NotImplementedError('only "normal" is available')

    return n


def proportion_effectsize(prop1, prop2, method="normal"):
    """
    Effect size for a test comparing two proportions

    for use in power function

    Parameters
    ----------
    prop1, prop2 : float or array_like
        The proportion value(s).

    Returns
    -------
    es : float or ndarray
        effect size for (transformed) prop1 - prop2

    Notes
    -----
    only method='normal' is implemented to match pwr.p2.test
    see http://www.statmethods.net/stats/power.html

    Effect size for `normal` is defined as ::

        2 * (arcsin(sqrt(prop1)) - arcsin(sqrt(prop2)))

    I think other conversions to normality can be used, but I need to check.

    Examples
    --------
    >>> import statsmodels.api as sm
    >>> sm.stats.proportion_effectsize(0.5, 0.4)
    0.20135792079033088
    >>> sm.stats.proportion_effectsize([0.3, 0.4, 0.5], 0.4)
    array([-0.21015893,  0.        ,  0.20135792])

    """
    if method != "normal":
        raise ValueError('only "normal" is implemented')

    es = 2 * (np.arcsin(np.sqrt(prop1)) - np.arcsin(np.sqrt(prop2)))
    return es


def std_prop(prop, nobs):
    """
    Standard error for the estimate of a proportion

    This is just ``np.sqrt(p * (1. - p) / nobs)``

    Parameters
    ----------
    prop : array_like
        proportion
    nobs : int, array_like
        number of observations

    Returns
    -------
    std : array_like
        standard error for a proportion of nobs independent observations
    """
    return np.sqrt(prop * (1.0 - prop) / nobs)


def _std_diff_prop(p1, p2, ratio=1):
    return np.sqrt(p1 * (1 - p1) + p2 * (1 - p2) / ratio)


def _power_ztost(
    mean_low,
    var_low,
    mean_upp,
    var_upp,
    mean_alt,
    var_alt,
    alpha=0.05,
    discrete=True,
    dist="norm",
    nobs=None,
    continuity=0,
    critval_continuity=0,
):
    """
    Generic statistical power function for normal based equivalence test

    This includes options to adjust the normal approximation and can use
    the binomial to evaluate the probability of the rejection region

    see power_ztost_prob for a description of the options
    """
    # TODO: refactor structure, separate norm and binom better
    if not isinstance(continuity, tuple):
        continuity = (continuity, continuity)
    crit = stats.norm.isf(alpha)
    k_low = mean_low + np.sqrt(var_low) * crit
    k_upp = mean_upp - np.sqrt(var_upp) * crit
    if discrete or dist == "binom":
        k_low = np.ceil(k_low * nobs + 0.5 * critval_continuity)
        k_upp = np.trunc(k_upp * nobs - 0.5 * critval_continuity)
        if dist == "norm":
            # need proportion
            k_low = (k_low) * 1.0 / nobs  # -1 to match PASS
            k_upp = k_upp * 1.0 / nobs
    #    else:
    #        if dist == 'binom':
    #            # need counts
    #            k_low *= nobs
    #            k_upp *= nobs
    # print mean_low, np.sqrt(var_low), crit, var_low
    # print mean_upp, np.sqrt(var_upp), crit, var_upp
    if np.any(k_low > k_upp):  # vectorize
        import warnings

        warnings.warn("no overlap, power is zero", HypothesisTestWarning, stacklevel=2)
    std_alt = np.sqrt(var_alt)
    z_low = (k_low - mean_alt - continuity[0] * 0.5 / nobs) / std_alt
    z_upp = (k_upp - mean_alt + continuity[1] * 0.5 / nobs) / std_alt
    if dist == "norm":
        power = stats.norm.cdf(z_upp) - stats.norm.cdf(z_low)
    elif dist == "binom":
        power = stats.binom.cdf(k_upp, nobs, mean_alt) - stats.binom.cdf(
            k_low - 1, nobs, mean_alt
        )
    return power, (k_low, k_upp, z_low, z_upp)


def binom_tost(count, nobs, low, upp):
    """
    Exact TOST test for one proportion using binomial distribution

    Parameters
    ----------
    count : {int, array_like}
        the number of successes in nobs trials.
    nobs : int
        the number of trials or observations.
    low, upp : floats
        lower and upper limit of equivalence region

    Returns
    -------
    pvalue : float
        p-value of equivalence test
    pval_low, pval_upp : floats
        p-values of lower and upper one-sided tests

    """
    # binom_test_stat only returns pval
    tt1 = binom_test(count, nobs, alternative="larger", prop=low)
    tt2 = binom_test(count, nobs, alternative="smaller", prop=upp)
    return (
        np.maximum(tt1, tt2),
        tt1,
        tt2,
    )


def binom_tost_reject_interval(low, upp, nobs, alpha=0.05):
    """
    Rejection region for binomial TOST

    The interval includes the end points,
    `reject` if and only if `r_low <= x <= r_upp`.

    The interval might be empty with `r_upp < r_low`.

    Parameters
    ----------
    low, upp : floats
        lower and upper limit of equivalence region
    nobs : int
        the number of trials or observations.

    Returns
    -------
    x_low, x_upp : float
        lower and upper bound of rejection region

    """
    x_low = stats.binom.isf(alpha, nobs, low) + 1
    x_upp = stats.binom.ppf(alpha, nobs, upp) - 1
    return x_low, x_upp


def binom_test_reject_interval(value, nobs, alpha=0.05, alternative="two-sided"):
    """
    Rejection region for binomial test for one sample proportion

    The interval includes the end points of the rejection region.

    Parameters
    ----------
    value : float
        proportion under the Null hypothesis
    nobs : int
        the number of trials or observations.

    Returns
    -------
    x_low, x_upp : int
        lower and upper bound of rejection region
    """
    if alternative in ["2s", "two-sided"]:
        alternative = "2s"  # normalize alternative name
        alpha = alpha / 2

    if alternative in ["2s", "smaller"]:
        x_low = stats.binom.ppf(alpha, nobs, value) - 1
    else:
        x_low = 0
    if alternative in ["2s", "larger"]:
        x_upp = stats.binom.isf(alpha, nobs, value) + 1
    else:
        x_upp = nobs

    return int(x_low), int(x_upp)


def binom_test(count, nobs, prop=0.5, alternative="two-sided"):
    """
    Perform a test that the probability of success is p.

    This is an exact, two-sided test of the null hypothesis
    that the probability of success in a Bernoulli experiment
    is `p`.

    Parameters
    ----------
    count : {int, array_like}
        the number of successes in nobs trials.
    nobs : int
        the number of trials or observations.
    prop : float, optional
        The probability of success under the null hypothesis,
        `0 <= prop <= 1`. The default value is `prop = 0.5`
    alternative : str in ['two-sided', 'smaller', 'larger']
        alternative hypothesis, which can be two-sided or either one of the
        one-sided tests.

    Returns
    -------
    p-value : float
        The p-value of the hypothesis test

    Notes
    -----
    This uses scipy.stats.binom_test for the two-sided alternative.
    """

    if np.any(prop > 1.0) or np.any(prop < 0.0):
        raise ValueError("p must be in range [0,1]")
    if alternative in ["2s", "two-sided"]:
        try:
            pval = stats.binomtest(count, n=nobs, p=prop).pvalue
        except AttributeError:
            # Remove after min SciPy >= 1.7
            pval = stats.binom_test(count, n=nobs, p=prop)
    elif alternative in ["l", "larger"]:
        pval = stats.binom.sf(count - 1, nobs, prop)
    elif alternative in ["s", "smaller"]:
        pval = stats.binom.cdf(count, nobs, prop)
    else:
        raise ValueError(
            "alternative not recognized\nshould be two-sided, larger or smaller"
        )
    return pval


def power_binom_tost(low, upp, nobs, p_alt=None, alpha=0.05):
    if p_alt is None:
        p_alt = 0.5 * (low + upp)
    x_low, x_upp = binom_tost_reject_interval(low, upp, nobs, alpha=alpha)
    power = stats.binom.cdf(x_upp, nobs, p_alt) - stats.binom.cdf(
        x_low - 1, nobs, p_alt
    )
    return power


def power_ztost_prop(
    low,
    upp,
    nobs,
    p_alt,
    alpha=0.05,
    dist="norm",
    variance_prop=None,
    discrete=True,
    continuity=0,
    critval_continuity=0,
):
    """
    Power of proportions equivalence test based on normal distribution

    Parameters
    ----------
    low, upp : floats
        lower and upper limit of equivalence region
    nobs : int
        number of observations
    p_alt : float in (0,1)
        proportion under the alternative
    alpha : float in (0,1)
        significance level of the test
    dist : str in ['norm', 'binom']
        This defines the distribution to evaluate the power of the test. The
        critical values of the TOST test are always based on the normal
        approximation, but the distribution for the power can be either the
        normal (default) or the binomial (exact) distribution.
    variance_prop : None or float in (0,1)
        If this is None, then the variances for the two one sided tests are
        based on the proportions equal to the equivalence limits.
        If variance_prop is given, then it is used to calculate the variance
        for the TOST statistics. If this is based on an sample, then the
        estimated proportion can be used.
    discrete : bool
        If true, then the critical values of the rejection region are converted
        to integers. If dist is "binom", this is automatically assumed.
        If discrete is false, then the TOST critical values are used as
        floating point numbers, and the power is calculated based on the
        rejection region that is not discretized.
    continuity : bool or float
        adjust the rejection region for the normal power probability. This has
        and effect only if ``dist='norm'``
    critval_continuity : bool or float
        If this is non-zero, then the critical values of the tost rejection
        region are adjusted before converting to integers. This affects both
        distributions, ``dist='norm'`` and ``dist='binom'``.

    Returns
    -------
    power : float
        statistical power of the equivalence test.
    (k_low, k_upp, z_low, z_upp) : tuple of floats
        critical limits in intermediate steps
        temporary return, will be changed

    Notes
    -----
    In small samples the power for the ``discrete`` version, has a sawtooth
    pattern as a function of the number of observations. As a consequence,
    small changes in the number of observations or in the normal approximation
    can have a large effect on the power.

    ``continuity`` and ``critval_continuity`` are added to match some results
    of PASS, and are mainly to investigate the sensitivity of the ztost power
    to small changes in the rejection region. From my interpretation of the
    equations in the SAS manual, both are zero in SAS.

    works vectorized

    **verification:**

    The ``dist='binom'`` results match PASS,
    The ``dist='norm'`` results look reasonable, but no benchmark is available.

    References
    ----------
    SAS Manual: Chapter 68: The Power Procedure, Computational Resources
    PASS Chapter 110: Equivalence Tests for One Proportion.

    """
    mean_low = low
    var_low = std_prop(low, nobs) ** 2
    mean_upp = upp
    var_upp = std_prop(upp, nobs) ** 2
    mean_alt = p_alt
    var_alt = std_prop(p_alt, nobs) ** 2
    if variance_prop is not None:
        var_low = var_upp = std_prop(variance_prop, nobs) ** 2
    power = _power_ztost(
        mean_low,
        var_low,
        mean_upp,
        var_upp,
        mean_alt,
        var_alt,
        alpha=alpha,
