B
    dT                 @   s\  d Z ddlZddlZddlmZ ddlmZ dddd	d
dddddddddddgZd	gZ	dgdgdZ
deded  Zde Zdd Zedddd?d!dZedddd@d#d	Zd$d% Zed&dddAd(d
ZdBd*dZdCd+dZdDd,dZdEd.dZd/d0 Zd1d2 Zd3d4 Zd5d Zd6d7 d8fd9dZd:d Zd;d Zd<d Zd=d Z d>d Z!dS )Fa2  
This module contains simple statistical algorithms that are
straightforwardly implemented as a single python function (or family of
functions).

This module should generally not be used directly.  Everything in
`__all__` is imported into `astropy.stats`, and hence that package
should be used for access.
    N)deprecated_renamed_argument   )_statsgaussian_fwhm_to_sigmagaussian_sigma_to_fwhmbinom_conf_intervalbinned_binom_proportionpoisson_conf_intervalmedian_absolute_deviationmad_stdsignal_to_noise_oir_ccd	bootstrapkuiper
kuiper_two!kuiper_false_positive_probabilitycdf_from_intervalsinterval_overlap_lengthhistogram_intervalsfold_intervalsscipy)r   r	   g       @      ?c                s   t | tjrt| } n
t| } t  ttfs6 f t | j }tj	j
 | t| j fddt|D }| |S )aa  
    Expand the shape of an array.

    Insert a new axis that will appear at the `axis` position in the
    expanded array shape.

    This function allows for tuple axis arguments.
    ``numpy.expand_dims`` currently does not allow that, but it will in
    numpy v1.18 (https://github.com/numpy/numpy/pull/14051).
    ``_expand_dims`` can be replaced with ``numpy.expand_dims`` when the
    minimum support numpy version is v1.18.

    Parameters
    ----------
    data : array-like
        Input array.
    axis : int or tuple of int
        Position in the expanded axes where the new axis (or axes) is
        placed.  A tuple of axes is now supported.  Out of range axes as
        described above are now forbidden and raise an `AxisError`.

    Returns
    -------
    result : ndarray
        View of ``data`` with the number of dimensions increased.
    c                s    g | ]}| krd nt qS )r   )next).0Zax)axisshape_it `/work/yifan.wang/ringdown/master-ringdown-env/lib/python3.7/site-packages/astropy/stats/funcs.py
<listcomp>U   s    z _expand_dims.<locals>.<listcomp>)
isinstancenpmatrixasarray
asanyarraytuplelistlenndimcorenumericZnormalize_axis_tupleitershaperangeZreshape)datar   Zout_ndimr*   r   )r   r   r   _expand_dims-   s    

r-   confconfidence_levelz4.0鲘?wilsonc             C   s  |dk s|dkrt dd| }t| t} t|t}|dk rTt d| dk  sl| |k rtt d|dks|dkrdd	lm} td
t	||d }| t
} |t
}| | }|dkrX| |d d
  ||d   }|t| ||d   t|d|  |d d|    }	t||	 ||	 g}
d|
|
dk < d|
|
dk< n4|}|t|d|  |  }	t||	 ||	 g}
n|dks|dkrddlm} |dkr|| d ||  d d| }|| d ||  d dd|  }n<|| d ||  d d| }|| d ||  d dd|  }|jdkr`| dkrPd}n| |krxd}nd|| dk< d|| |k< t||g}
nt d|d|
S )a!  Binomial proportion confidence interval given k successes,
    n trials.

    Parameters
    ----------
    k : int or numpy.ndarray
        Number of successes (0 <= ``k`` <= ``n``).
    n : int or numpy.ndarray
        Number of trials (``n`` > 0).  If both ``k`` and ``n`` are arrays,
        they must have the same shape.
    confidence_level : float, optional
        Desired probability content of interval. Default is 0.68269,
        corresponding to 1 sigma in a 1-dimensional Gaussian distribution.
        Confidence level must be in range [0, 1].
    interval : {'wilson', 'jeffreys', 'flat', 'wald'}, optional
        Formula used for confidence interval. See notes for details.  The
        ``'wilson'`` and ``'jeffreys'`` intervals generally give similar
        results, while 'flat' is somewhat different, especially for small
        values of ``n``.  ``'wilson'`` should be somewhat faster than
        ``'flat'`` or ``'jeffreys'``.  The 'wald' interval is generally not
        recommended.  It is provided for comparison purposes.  Default is
        ``'wilson'``.

    Returns
    -------
    conf_interval : ndarray
        ``conf_interval[0]`` and ``conf_interval[1]`` correspond to the lower
        and upper limits, respectively, for each element in ``k``, ``n``.

    Notes
    -----
    In situations where a probability of success is not known, it can
    be estimated from a number of trials (n) and number of
    observed successes (k). For example, this is done in Monte
    Carlo experiments designed to estimate a detection efficiency. It
    is simple to take the sample proportion of successes (k/n)
    as a reasonable best estimate of the true probability
    :math:`\epsilon`. However, deriving an accurate confidence
    interval on :math:`\epsilon` is non-trivial. There are several
    formulas for this interval (see [1]_). Four intervals are implemented
    here:

    **1. The Wilson Interval.** This interval, attributed to Wilson [2]_,
    is given by

    .. math::

        CI_{\rm Wilson} = \frac{k + \kappa^2/2}{n + \kappa^2}
        \pm \frac{\kappa n^{1/2}}{n + \kappa^2}
        ((\hat{\epsilon}(1 - \hat{\epsilon}) + \kappa^2/(4n))^{1/2}

    where :math:`\hat{\epsilon} = k / n` and :math:`\kappa` is the
    number of standard deviations corresponding to the desired
    confidence interval for a *normal* distribution (for example,
    1.0 for a confidence interval of 68.269%). For a
    confidence interval of 100(1 - :math:`\alpha`)%,

    .. math::

        \kappa = \Phi^{-1}(1-\alpha/2) = \sqrt{2}{\rm erf}^{-1}(1-\alpha).

    **2. The Jeffreys Interval.** This interval is derived by applying
    Bayes' theorem to the binomial distribution with the
    noninformative Jeffreys prior [3]_, [4]_. The noninformative Jeffreys
    prior is the Beta distribution, Beta(1/2, 1/2), which has the density
    function

    .. math::

        f(\epsilon) = \pi^{-1} \epsilon^{-1/2}(1-\epsilon)^{-1/2}.

    The justification for this prior is that it is invariant under
    reparameterizations of the binomial proportion.
    The posterior density function is also a Beta distribution: Beta(k
    + 1/2, n - k + 1/2). The interval is then chosen so that it is
    *equal-tailed*: Each tail (outside the interval) contains
    :math:`\alpha`/2 of the posterior probability, and the interval
    itself contains 1 - :math:`\alpha`. This interval must be
    calculated numerically. Additionally, when k = 0 the lower limit
    is set to 0 and when k = n the upper limit is set to 1, so that in
    these cases, there is only one tail containing :math:`\alpha`/2
    and the interval itself contains 1 - :math:`\alpha`/2 rather than
    the nominal 1 - :math:`\alpha`.

    **3. A Flat prior.** This is similar to the Jeffreys interval,
    but uses a flat (uniform) prior on the binomial proportion
    over the range 0 to 1 rather than the reparametrization-invariant
    Jeffreys prior.  The posterior density function is a Beta distribution:
    Beta(k + 1, n - k + 1).  The same comments about the nature of the
    interval (equal-tailed, etc.) also apply to this option.

    **4. The Wald Interval.** This interval is given by

    .. math::

       CI_{\rm Wald} = \hat{\epsilon} \pm
       \kappa \sqrt{\frac{\hat{\epsilon}(1-\hat{\epsilon})}{n}}

    The Wald interval gives acceptable results in some limiting
    cases. Particularly, when n is very large, and the true proportion
    :math:`\epsilon` is not "too close" to 0 or 1. However, as the
    later is not verifiable when trying to estimate :math:`\epsilon`,
    this is not very helpful. Its use is not recommended, but it is
    provided here for comparison purposes due to its prevalence in
    everyday practical statistics.

    This function requires ``scipy`` for all interval types.

    References
    ----------
    .. [1] 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

    .. [2] Wilson, E. B. (1927). "Probable inference, the law of
       succession, and statistical inference". Journal of the American
       Statistical Association 22: 209-212.

    .. [3] Jeffreys, Harold (1946). "An Invariant Form for the Prior
       Probability in Estimation Problems". Proc. R. Soc. Lond.. A 24 186
       (1007): 453-461. doi:10.1098/rspa.1946.0056

    .. [4] Jeffreys, Harold (1998). Theory of Probability. Oxford
       University Press, 3rd edition. ISBN 978-0198503682

    Examples
    --------
    Integer inputs return an array with shape (2,):

    >>> binom_conf_interval(4, 5, interval='wilson')  # doctest: +FLOAT_CMP
    array([0.57921724, 0.92078259])

    Arrays of arbitrary dimension are supported. The Wilson and Jeffreys
    intervals give similar results, even for small k, n:

    >>> binom_conf_interval([1, 2], 5, interval='wilson')  # doctest: +FLOAT_CMP
    array([[0.07921741, 0.21597328],
           [0.42078276, 0.61736012]])

    >>> binom_conf_interval([1, 2,], 5, interval='jeffreys')  # doctest: +FLOAT_CMP
    array([[0.0842525 , 0.21789949],
           [0.42218001, 0.61753691]])

    >>> binom_conf_interval([1, 2], 5, interval='flat')  # doctest: +FLOAT_CMP
    array([[0.12139799, 0.24309021],
           [0.45401727, 0.61535699]])

    In contrast, the Wald interval gives poor results for small k, n.
    For k = 0 or k = n, the interval always has zero length.

    >>> binom_conf_interval([1, 2], 5, interval='wald')  # doctest: +FLOAT_CMP
    array([[0.02111437, 0.18091075],
           [0.37888563, 0.61908925]])

    For confidence intervals approaching 1, the Wald interval for
    0 < k < n can give intervals that extend outside [0, 1]:

    >>> binom_conf_interval([1, 2], 5, interval='wald', confidence_level=0.99)  # doctest: +FLOAT_CMP
    array([[-0.26077835, -0.16433593],
           [ 0.66077835,  0.96433593]])

    g        g      ?z*confidence_level must be between 0. and 1.r   zn must be positivezk must be in {0, 1, .., n}r1   Zwald)erfinvg       @g    _B   r      ZjeffreysZflat)
betaincinvg      ?zUnrecognized interval: s)
ValueErrorr   r!   astypeintanyscipy.specialr2   sqrtminfloatarrayr5   r&   )knr/   intervalalphar2   kappapZmidpointZ
halflengthconf_intervalr5   Z
lowerboundZ
upperboundr   r   r   r   Z   sZ     %


$


" 


   c             C   s   t | } t |t}| j|jkr.tdt j| ||d\}}t j| | |d\}}|dd |dd  d }	|	|dd  }
|dk}|	| }	|
| }
|| }|| }|| }t||||d	}t || }|	|
||fS )
a   Binomial proportion and confidence interval in bins of a continuous
    variable ``x``.

    Given a set of datapoint pairs where the ``x`` values are
    continuously distributed and the ``success`` values are binomial
    ("success / failure" or "true / false"), place the pairs into
    bins according to ``x`` value and calculate the binomial proportion
    (fraction of successes) and confidence interval in each bin.

    Parameters
    ----------
    x : sequence
        Values.
    success : sequence of bool
        Success (`True`) or failure (`False`) corresponding to each value
        in ``x``.  Must be same length as ``x``.
    bins : int or sequence of scalar, optional
        If bins is an int, it defines the number of equal-width bins
        in the given range (10, by default). If bins is a sequence, it
        defines the bin edges, including the rightmost edge, allowing
        for non-uniform bin widths (in this case, 'range' is ignored).
    range : (float, float), optional
        The lower and upper range of the bins. If `None` (default),
        the range is set to ``(x.min(), x.max())``. Values outside the
        range are ignored.
    confidence_level : float, optional
        Must be in range [0, 1].
        Desired probability content in the confidence
        interval ``(p - perr[0], p + perr[1])`` in each bin. Default is
        0.68269.
    interval : {'wilson', 'jeffreys', 'flat', 'wald'}, optional
        Formula used to calculate confidence interval on the
        binomial proportion in each bin. See `binom_conf_interval` for
        definition of the intervals.  The 'wilson', 'jeffreys',
        and 'flat' intervals generally give similar results.  'wilson'
        should be somewhat faster, while 'jeffreys' and 'flat' are
        marginally superior, but differ in the assumed prior.
        The 'wald' interval is generally not recommended.
        It is provided for comparison purposes. Default is 'wilson'.

    Returns
    -------
    bin_ctr : ndarray
        Central value of bins. Bins without any entries are not returned.
    bin_halfwidth : ndarray
        Half-width of each bin such that ``bin_ctr - bin_halfwidth`` and
        ``bin_ctr + bins_halfwidth`` give the left and right side of each bin,
        respectively.
    p : ndarray
        Efficiency in each bin.
    perr : ndarray
        2-d array of shape (2, len(p)) representing the upper and lower
        uncertainty on p in each bin.

    Notes
    -----
    This function requires ``scipy`` for all interval types.

    See Also
    --------
    binom_conf_interval : Function used to estimate confidence interval in
                          each bin.

    Examples
    --------
    Suppose we wish to estimate the efficiency of a survey in
    detecting astronomical sources as a function of magnitude (i.e.,
    the probability of detecting a source given its magnitude). In a
    realistic case, we might prepare a large number of sources with
    randomly selected magnitudes, inject them into simulated images,
    and then record which were detected at the end of the reduction
    pipeline. As a toy example, we generate 100 data points with
    randomly selected magnitudes between 20 and 30 and "observe" them
    with a known detection function (here, the error function, with
    50% detection probability at magnitude 25):

    >>> from scipy.special import erf
    >>> from scipy.stats.distributions import binom
    >>> def true_efficiency(x):
    ...     return 0.5 - 0.5 * erf((x - 25.) / 2.)
    >>> mag = 20. + 10. * np.random.rand(100)
    >>> detected = binom.rvs(1, true_efficiency(mag))
    >>> bins, binshw, p, perr = binned_binom_proportion(mag, detected, bins=20)
    >>> plt.errorbar(bins, p, xerr=binshw, yerr=perr, ls='none', marker='o',
    ...              label='estimate')

    .. plot::

       import numpy as np
       from scipy.special import erf
       from scipy.stats.distributions import binom
       import matplotlib.pyplot as plt
       from astropy.stats import binned_binom_proportion
       def true_efficiency(x):
           return 0.5 - 0.5 * erf((x - 25.) / 2.)
       np.random.seed(400)
       mag = 20. + 10. * np.random.rand(100)
       np.random.seed(600)
       detected = binom.rvs(1, true_efficiency(mag))
       bins, binshw, p, perr = binned_binom_proportion(mag, detected, bins=20)
       plt.errorbar(bins, p, xerr=binshw, yerr=perr, ls='none', marker='o',
                    label='estimate')
       X = np.linspace(20., 30., 1000)
       plt.plot(X, true_efficiency(X), label='true efficiency')
       plt.ylim(0., 1.)
       plt.title('Detection efficiency vs magnitude')
       plt.xlabel('Magnitude')
       plt.ylabel('Detection efficiency')
       plt.legend()
       plt.show()

    The above example uses the Wilson confidence interval to calculate
    the uncertainty ``perr`` in each bin (see the definition of various
    confidence intervals in `binom_conf_interval`). A commonly used
    alternative is the Wald interval. However, the Wald interval can
    give nonsensical uncertainties when the efficiency is near 0 or 1,
    and is therefore **not** recommended. As an illustration, the
    following example shows the same data as above but uses the Wald
    interval rather than the Wilson interval to calculate ``perr``:

    >>> bins, binshw, p, perr = binned_binom_proportion(mag, detected, bins=20,
    ...                                                 interval='wald')
    >>> plt.errorbar(bins, p, xerr=binshw, yerr=perr, ls='none', marker='o',
    ...              label='estimate')

    .. plot::

       import numpy as np
       from scipy.special import erf
       from scipy.stats.distributions import binom
       import matplotlib.pyplot as plt
       from astropy.stats import binned_binom_proportion
       def true_efficiency(x):
           return 0.5 - 0.5 * erf((x - 25.) / 2.)
       np.random.seed(400)
       mag = 20. + 10. * np.random.rand(100)
       np.random.seed(600)
       detected = binom.rvs(1, true_efficiency(mag))
       bins, binshw, p, perr = binned_binom_proportion(mag, detected, bins=20,
                                                       interval='wald')
       plt.errorbar(bins, p, xerr=binshw, yerr=perr, ls='none', marker='o',
                    label='estimate')
       X = np.linspace(20., 30., 1000)
       plt.plot(X, true_efficiency(X), label='true efficiency')
       plt.ylim(0., 1.)
       plt.title('The Wald interval can give nonsensical uncertainties')
       plt.xlabel('Magnitude')
       plt.ylabel('Detection efficiency')
       plt.legend()
       plt.show()

    z!sizes of x and success must match)binsr+   )rH   Nr   g       @r   )r/   rB   )	r   Zravelr8   boolr*   r7   Z	histogramr   abs)xsuccessrH   r+   r/   rB   rA   Z	bin_edgesr@   Zbin_ctrZbin_halfwidthZvalidrE   ZboundsZperrr   r   r   r   >  s$     
c             C   sF   | dkrt d| |dkr,t d| |d k	rBt d| d S )Nr   z$Only sigma=1 supported for interval r   z&background not supported for interval z,confidence_level not supported for interval )r7   )sigma
backgroundr/   namer   r   r   _check_poisson_conf_inputs  s    rQ   Z	conflevelroot-nc             C   s  t | st | } |dkrPt|||| t | t |  | t |  g}n|dkrt|||| t | t |  | t |  g}t | r| dkrd|d< nd|d| dkf< n|dkrt|||| t | d t | d  | d t | d  g}n|dkrTt|||| t | d t | d	  | d t | d	  g}n|d
krtd||| ddl}|jj	|}t d|j
d|  | d|j
d|  d | g}t | r| dkrd|d< nd|d| dkf< n|dkrt | rt| ts8tdnt| jjt js8td|dkrPtd|t |}t |dkszt |dkrtdt |}t |dk rtdt jtdd| ||}t |}ntd| |S )a  Poisson parameter confidence interval given observed counts

    Parameters
    ----------
    n : int or numpy.ndarray
        Number of counts (0 <= ``n``).
    interval : {'root-n','root-n-0','pearson','sherpagehrels','frequentist-confidence', 'kraft-burrows-nousek'}, optional
        Formula used for confidence interval. See notes for details.
        Default is ``'root-n'``.
    sigma : float, optional
        Number of sigma for confidence interval; only supported for
        the 'frequentist-confidence' mode.
    background : float, optional
        Number of counts expected from the background; only supported for
        the 'kraft-burrows-nousek' mode. This number is assumed to be determined
        from a large region so that the uncertainty on its value is negligible.
    confidence_level : float, optional
        Confidence level between 0 and 1; only supported for the
        'kraft-burrows-nousek' mode.

    Returns
    -------
    conf_interval : ndarray
        ``conf_interval[0]`` and ``conf_interval[1]`` correspond to the lower
        and upper limits, respectively, for each element in ``n``.

    Notes
    -----

    The "right" confidence interval to use for Poisson data is a
    matter of debate. The CDF working group `recommends
    <https://www-cdf.fnal.gov/physics/statistics/notes/pois_eb.txt>`_
    using root-n throughout, largely in the interest of
    comprehensibility, but discusses other possibilities. The ATLAS
    group also `discusses
    <http://www.pp.rhul.ac.uk/~cowan/atlas/ErrorBars.pdf>`_  several
    possibilities but concludes that no single representation is
    suitable for all cases.  The suggestion has also been `floated
    <https://ui.adsabs.harvard.edu/abs/2012EPJP..127...24A>`_ that error
    bars should be attached to theoretical predictions instead of
    observed data, which this function will not help with (but it's
    easy; then you really should use the square root of the theoretical
    prediction).

    The intervals implemented here are:

    **1. 'root-n'** This is a very widely used standard rule derived
    from the maximum-likelihood estimator for the mean of the Poisson
    process. While it produces questionable results for small n and
    outright wrong results for n=0, it is standard enough that people are
    (supposedly) used to interpreting these wonky values. The interval is

    .. math::

        CI = (n-\sqrt{n}, n+\sqrt{n})

    **2. 'root-n-0'** This is identical to the above except that where
    n is zero the interval returned is (0,1).

    **3. 'pearson'** This is an only-slightly-more-complicated rule
    based on Pearson's chi-squared rule (as `explained
    <https://www-cdf.fnal.gov/physics/statistics/notes/pois_eb.txt>`_ by
    the CDF working group). It also has the nice feature that if your
    theory curve touches an endpoint of the interval, then your data
    point is indeed one sigma away. The interval is

    .. math::

        CI = (n+0.5-\sqrt{n+0.25}, n+0.5+\sqrt{n+0.25})

    **4. 'sherpagehrels'** This rule is used by default in the fitting
    package 'sherpa'. The `documentation
    <https://cxc.harvard.edu/sherpa4.4/statistics/#chigehrels>`_ claims
    it is based on a numerical approximation published in `Gehrels
    (1986) <https://ui.adsabs.harvard.edu/abs/1986ApJ...303..336G>`_ but it
    does not actually appear there.  It is symmetrical, and while the
    upper limits are within about 1% of those given by
    'frequentist-confidence', the lower limits can be badly wrong. The
    interval is

    .. math::

        CI = (n-1-\sqrt{n+0.75}, n+1+\sqrt{n+0.75})

    **5. 'frequentist-confidence'** These are frequentist central
    confidence intervals:

    .. math::

        CI = (0.5 F_{\chi^2}^{-1}(\alpha;2n),
              0.5 F_{\chi^2}^{-1}(1-\alpha;2(n+1)))

    where :math:`F_{\chi^2}^{-1}` is the quantile of the chi-square
    distribution with the indicated number of degrees of freedom and
    :math:`\alpha` is the one-tailed probability of the normal
    distribution (at the point given by the parameter 'sigma'). See
    `Maxwell (2011)
    <https://ui.adsabs.harvard.edu/abs/2011arXiv1102.0822M>`_ for further
    details.

    **6. 'kraft-burrows-nousek'** This is a Bayesian approach which allows
    for the presence of a known background :math:`B` in the source signal
    :math:`N`.
    For a given confidence level :math:`CL` the confidence interval
    :math:`[S_\mathrm{min}, S_\mathrm{max}]` is given by:

    .. math::

       CL = \int^{S_\mathrm{max}}_{S_\mathrm{min}} f_{N,B}(S)dS

    where the function :math:`f_{N,B}` is:

    .. math::

       f_{N,B}(S) = C \frac{e^{-(S+B)}(S+B)^N}{N!}

    and the normalization constant :math:`C`:

    .. math::

       C = \left[ \int_0^\infty \frac{e^{-(S+B)}(S+B)^N}{N!} dS \right] ^{-1}
       = \left( \sum^N_{n=0} \frac{e^{-B}B^n}{n!}  \right)^{-1}

    See `Kraft, Burrows, and Nousek (1991)
    <https://ui.adsabs.harvard.edu/abs/1991ApJ...374..344K>`_ for further
    details.

    These formulas implement a positive, uniform prior.
    `Kraft, Burrows, and Nousek (1991)
    <https://ui.adsabs.harvard.edu/abs/1991ApJ...374..344K>`_ discuss this
    choice in more detail and show that the problem is relatively
    insensitive to the choice of prior.

    This function has an optional dependency: Either `Scipy
    <https://www.scipy.org/>`_ or `mpmath <http://mpmath.org/>`_  need
    to be available (Scipy works only for N < 100).
    This code is very intense numerically, which makes it much slower than
    the other methods, in particular for large count numbers (above 1000
    even with ``mpmath``). Fortunately, some of the other methods or a
    Gaussian approximation usually work well in this regime.

    Examples
    --------
    >>> poisson_conf_interval(np.arange(10), interval='root-n').T
    array([[  0.        ,   0.        ],
           [  0.        ,   2.        ],
           [  0.58578644,   3.41421356],
           [  1.26794919,   4.73205081],
           [  2.        ,   6.        ],
           [  2.76393202,   7.23606798],
           [  3.55051026,   8.44948974],
           [  4.35424869,   9.64575131],
           [  5.17157288,  10.82842712],
           [  6.        ,  12.        ]])

    >>> poisson_conf_interval(np.arange(10), interval='root-n-0').T
    array([[  0.        ,   1.        ],
           [  0.        ,   2.        ],
           [  0.58578644,   3.41421356],
           [  1.26794919,   4.73205081],
           [  2.        ,   6.        ],
           [  2.76393202,   7.23606798],
           [  3.55051026,   8.44948974],
           [  4.35424869,   9.64575131],
           [  5.17157288,  10.82842712],
           [  6.        ,  12.        ]])

    >>> poisson_conf_interval(np.arange(10), interval='pearson').T
    array([[  0.        ,   1.        ],
           [  0.38196601,   2.61803399],
           [  1.        ,   4.        ],
           [  1.69722436,   5.30277564],
           [  2.43844719,   6.56155281],
           [  3.20871215,   7.79128785],
           [  4.        ,   9.        ],
           [  4.8074176 ,  10.1925824 ],
           [  5.62771868,  11.37228132],
           [  6.45861873,  12.54138127]])

    >>> poisson_conf_interval(
    ...     np.arange(10), interval='frequentist-confidence').T
    array([[  0.        ,   1.84102165],
           [  0.17275378,   3.29952656],
           [  0.70818544,   4.63785962],
           [  1.36729531,   5.91818583],
           [  2.08566081,   7.16275317],
           [  2.84030886,   8.38247265],
           [  3.62006862,   9.58364155],
           [  4.41852954,  10.77028072],
           [  5.23161394,  11.94514152],
           [  6.05653896,  13.11020414]])

    >>> poisson_conf_interval(
    ...     7, interval='frequentist-confidence').T
    array([  4.41852954,  10.77028072])

    >>> poisson_conf_interval(
    ...     10, background=1.5, confidence_level=0.95,
    ...     interval='kraft-burrows-nousek').T  # doctest: +FLOAT_CMP
    array([[ 3.47894005, 16.113329533]])

    zroot-nzroot-n-0r   r   Zpearsong      ?g      ?Zsherpagehrelsg      ?zfrequentist-confidenceg      ?Nr3   zkraft-burrows-nousekz!Number of counts must be integer.z7Set confidence_level for method {}. (sigma is ignored.)z2confidence_level must be a number between 0 and 1.zBackground must be >= 0.T)cachez1Invalid method for Poisson confidence intervals: )r   Zisscalarr"   rQ   r?   r<   Zscipy.statsstatsZnormZsfZchi2ZppfZisfr   r9   	TypeError
issubclassdtypetypeintegerr7   formatr:   Z	vectorize_kraft_burrows_nousekZvstack)rA   rB   rN   rO   r/   rF   r   rC   r   r   r   r	     sj     O






"





 
Fc             C   s   |dkr\t | tjjr@d}tjj}|rZtjjt| | dd} q`|rPd}tj}q`d}tj}nd}t| } || |d}|dk	rt	||d}|t
| | |dd}|dkrtj|r| }ntj|r|s|jtjd}|S )a  
    Calculate the median absolute deviation (MAD).

    The MAD is defined as ``median(abs(a - median(a)))``.

    Parameters
    ----------
    data : array-like
        Input array or object that can be converted to an array.
    axis : None, int, or tuple of int, optional
        The axis or axes along which the MADs are computed.  The default
        (`None`) is to compute the MAD of the flattened array.
    func : callable, optional
        The function used to compute the median. Defaults to `numpy.ma.median`
        for masked arrays, otherwise to `numpy.median`.
    ignore_nan : bool
        Ignore NaN values (treat them as if they are not in the array) when
        computing the median.  This will use `numpy.ma.median` if ``axis`` is
        specified, or `numpy.nanmedian` if ``axis==None`` and numpy's version
        is >1.10 because nanmedian is slightly faster in this case.

    Returns
    -------
    mad : float or `~numpy.ndarray`
        The median absolute deviation of the input array.  If ``axis``
        is `None` then a scalar will be returned, otherwise a
        `~numpy.ndarray` will be returned.

    Examples
    --------
    Generate random variates from a Gaussian distribution and return the
    median absolute deviation for that distribution::

        >>> import numpy as np
        >>> from astropy.stats import median_absolute_deviation
        >>> rand = np.random.default_rng(12345)
        >>> from numpy.random import randn
        >>> mad = median_absolute_deviation(rand.standard_normal(1000))
        >>> print(mad)    # doctest: +FLOAT_CMP
        0.6829504282771885

    See Also
    --------
    mad_std
    NT)copyF)r   )r   Zoverwrite_input)Z
fill_value)r   r   maZMaskedArrayZmedianZmasked_whereisnanZ	nanmedianr"   r-   rK   ZisMaskedArrayitemZfillednan)r,   r   func
ignore_nanZ	is_maskedZdata_medianresultr   r   r   r
     s,    /

c             C   s   t | |||d}|d S )as  
    Calculate a robust standard deviation using the `median absolute
    deviation (MAD)
    <https://en.wikipedia.org/wiki/Median_absolute_deviation>`_.

    The standard deviation estimator is given by:

    .. math::

        \sigma \approx \frac{\textrm{MAD}}{\Phi^{-1}(3/4)}
            \approx 1.4826 \ \textrm{MAD}

    where :math:`\Phi^{-1}(P)` is the normal inverse cumulative
    distribution function evaluated at probability :math:`P = 3/4`.

    Parameters
    ----------
    data : array-like
        Data array or object that can be converted to an array.
    axis : None, int, or tuple of int, optional
        The axis or axes along which the robust standard deviations are
        computed.  The default (`None`) is to compute the robust
        standard deviation of the flattened array.
    func : callable, optional
        The function used to compute the median. Defaults to `numpy.ma.median`
        for masked arrays, otherwise to `numpy.median`.
    ignore_nan : bool
        Ignore NaN values (treat them as if they are not in the array) when
        computing the median.  This will use `numpy.ma.median` if ``axis`` is
        specified, or `numpy.nanmedian` if ``axis=None`` and numpy's version is
        >1.10 because nanmedian is slightly faster in this case.

    Returns
    -------
    mad_std : float or `~numpy.ndarray`
        The robust standard deviation of the input data.  If ``axis`` is
        `None` then a scalar will be returned, otherwise a
        `~numpy.ndarray` will be returned.

    Examples
    --------
    >>> import numpy as np
    >>> from astropy.stats import mad_std
    >>> rand = np.random.default_rng(12345)
    >>> madstd = mad_std(rand.normal(5, 2, (100, 100)))
    >>> print(madstd)    # doctest: +FLOAT_CMP
    1.984147963351707

    See Also
    --------
    biweight_midvariance, biweight_midcovariance, median_absolute_deviation
    )r   ra   rb   gtV?)r
   )r,   r   ra   rb   ZMADr   r   r   r   _  s    7c       	      C   sB   | | | }t | || ||| |    ||d   }|| S )a  Computes the signal to noise ratio for source being observed in the
    optical/IR using a CCD.

    Parameters
    ----------
    t : float or numpy.ndarray
        CCD integration time in seconds
    source_eps : float
        Number of electrons (photons) or DN per second in the aperture from the
        source. Note that this should already have been scaled by the filter
        transmission and the quantum efficiency of the CCD. If the input is in
        DN, then be sure to set the gain to the proper value for the CCD.
        If the input is in electrons per second, then keep the gain as its
        default of 1.0.
    sky_eps : float
        Number of electrons (photons) or DN per second per pixel from the sky
        background. Should already be scaled by filter transmission and QE.
        This must be in the same units as source_eps for the calculation to
        make sense.
    dark_eps : float
        Number of thermal electrons per second per pixel. If this is given in
        DN or ADU, then multiply by the gain to get the value in electrons.
    rd : float
        Read noise of the CCD in electrons. If this is given in
        DN or ADU, then multiply by the gain to get the value in electrons.
    npix : float
        Size of the aperture in pixels
    gain : float, optional
        Gain of the CCD. In units of electrons per DN.

    Returns
    ----------
    SNR : float or numpy.ndarray
        Signal to noise ratio calculated from the inputs
    r3   )r   r<   )	tZ
source_epsZsky_epsZdark_epsrdZnpixZgainsignalnoiser   r   r   r     s    % d   c             C   s   |dkr| j d }|dk s"|dk r*td|dkrN|f|f | j dd  }n0y|t|| f}W n tk
r|   |f}Y nX t|}xNt|D ]B}tjjd| j d |d}|dkr| | ||< q|| | ||< qW |S )a
  Performs bootstrap resampling on numpy arrays.

    Bootstrap resampling is used to understand confidence intervals of sample
    estimates. This function returns versions of the dataset resampled with
    replacement ("case bootstrapping"). These can all be run through a function
    or statistic to produce a distribution of values which can then be used to
    find the confidence intervals.

    Parameters
    ----------
    data : ndarray
        N-D array. The bootstrap resampling will be performed on the first
        index, so the first index should access the relevant information
        to be bootstrapped.
    bootnum : int, optional
        Number of bootstrap resamples
    samples : int, optional
        Number of samples in each resample. The default `None` sets samples to
        the number of datapoints
    bootfunc : function, optional
        Function to reduce the resampled data. Each bootstrap resample will
        be put through this function and the results returned. If `None`, the
        bootstrapped data will be returned

    Returns
    -------
    boot : ndarray

        If bootfunc is None, then each row is a bootstrap resample of the data.
        If bootfunc is specified, then the columns will correspond to the
        outputs of bootfunc.

    Examples
    --------
    Obtain a twice resampled array:

    >>> from astropy.stats import bootstrap
    >>> import numpy as np
    >>> from astropy.utils import NumpyRNGContext
    >>> bootarr = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 0])
    >>> with NumpyRNGContext(1):
    ...     bootresult = bootstrap(bootarr, 2)
    ...
    >>> bootresult  # doctest: +FLOAT_CMP
    array([[6., 9., 0., 6., 1., 1., 2., 8., 7., 0.],
           [3., 5., 6., 3., 5., 3., 5., 8., 8., 0.]])
    >>> bootresult.shape
    (2, 10)

    Obtain a statistic on the array

    >>> with NumpyRNGContext(1):
    ...     bootresult = bootstrap(bootarr, 2, bootfunc=np.mean)
    ...
    >>> bootresult  # doctest: +FLOAT_CMP
    array([4. , 4.6])

    Obtain a statistic with two outputs on the array

    >>> test_statistic = lambda x: (np.sum(x), np.mean(x))
    >>> with NumpyRNGContext(1):
    ...     bootresult = bootstrap(bootarr, 3, bootfunc=test_statistic)
    >>> bootresult  # doctest: +FLOAT_CMP
    array([[40. ,  4. ],
           [46. ,  4.6],
           [35. ,  3.5]])
    >>> bootresult.shape
    (3, 2)

    Obtain a statistic with two outputs on the array, keeping only the first
    output

    >>> bootfunc = lambda x:test_statistic(x)[0]
    >>> with NumpyRNGContext(1):
    ...     bootresult = bootstrap(bootarr, 3, bootfunc=bootfunc)
    ...
    >>> bootresult  # doctest: +FLOAT_CMP
    array([40., 46., 35.])
    >>> bootresult.shape
    (3,)

    Nr   r   z3neither 'samples' nor 'bootnum' can be less than 1.)lowhighsize)	r*   r7   r%   rU   r   emptyr+   randomrandint)r,   ZbootnumZsamplesZbootfuncZ
resultdimsbootiZbootarrr   r   r   r     s"    S

c                s   ddl m ddlm ddlm ddlm fdd}| tt		fdd	fd
dfdd
 
fdd}|  d}
| }||fS )a  Upper limit on a poisson count rate

    The implementation is based on Kraft, Burrows and Nousek
    `ApJ 374, 344 (1991) <https://ui.adsabs.harvard.edu/abs/1991ApJ...374..344K>`_.
    The XMM-Newton upper limit server uses the same formalism.

    Parameters
    ----------
    N : int or np.int32/np.int64
        Total observed count number
    B : float or np.float32/np.float64
        Background count rate (assumed to be known with negligible error
        from a large background area).
    CL : float or np.float32/np.float64
       Confidence level (number between 0 and 1)

    Returns
    -------
    S : source count limit

    Notes
    -----
    Requires :mod:`~scipy`. This implementation will cause Overflow Errors for
    about N > 100 (the exact limit depends on details of how scipy was
    compiled). See `~astropy.stats.mpmath_poisson_upper_limit` for an
    implementation that is slower, but can deal with arbitrarily high numbers
    since it is based on the `mpmath <http://mpmath.org/>`_ library.
    r   )brentq)quad)	factorial)expc                s<   t j| d t jd}d | t t |||   S )Nr   )rW   g      ?)r   arangeZfloat64sumpower)NBrA   )rt   rs   r   r   eqn8[  s    z)_scipy_kraft_burrows_nousek.<locals>.eqn8c                s"   | | } | ||    S )Nr   )Srx   ry   SpB)eqn8_resrt   factorial_Nr   r   eqn7f  s    z)_scipy_kraft_burrows_nousek.<locals>.eqn7c                s    | |||fddS )Ni  )argslimitr   )S_minS_maxrx   ry   )r   rr   r   r   	eqn9_leftj  s    z._scipy_kraft_burrows_nousek.<locals>.eqn9_leftc                sB   |  d kr dS  fddd  S dS )a  
        Kraft, Burrows and Nousek suggest to integrate from N-B in both
        directions at once, so that S_min and S_max move similarly (see
        the article for details). Here, this is implemented differently:
        Treat S_max as the optimization parameters in func and then
        calculate the matching s_min that has has eqn7(S_max) =
        eqn7(S_min) here.
        r   g        c                s   |   S )Nr   )rL   )ry   rx   r   y_S_maxr   r   <lambda>z      zA_scipy_kraft_burrows_nousek.<locals>.find_s_min.<locals>.<lambda>Nr   )r   rx   ry   )rq   r   )ry   rx   r   r   
find_s_minm  s    	z/_scipy_kraft_burrows_nousek.<locals>.find_s_minc                s&   |  }||  }|d  S )Nr   r   )r6   s_minout)ry   CLrx   r   r   r   r   ra   |  s    z)_scipy_kraft_burrows_nousek.<locals>.funcrh   )	Zscipy.optimizerq   Zscipy.integraterr   r;   rs   mathrt   r>   )rx   ry   r   rz   ra   r   r   r   )ry   r   rx   rq   r   r}   r   rt   rs   r~   r   rr   r   _scipy_kraft_burrows_nousek7  s    
r   c       	         s  ddl m}mm
mmmm |t|t  |tdfdd}| fddfdd	
fd
d	 	fdd}t	  d}x||dk r|d7 }qW 
||d |gdd}	| }t|t|fS )aq  Upper limit on a poisson count rate

    The implementation is based on Kraft, Burrows and Nousek in
    `ApJ 374, 344 (1991) <https://ui.adsabs.harvard.edu/abs/1991ApJ...374..344K>`_.
    The XMM-Newton upper limit server used the same formalism.

    Parameters
    ----------
    N : int or np.int32/np.int64
        Total observed count number
    B : float or np.float32/np.float64
        Background count rate (assumed to be known with negligible error
        from a large background area).
    CL : float or np.float32/np.float64
       Confidence level (number between 0 and 1)

    Returns
    -------
    S : source count limit

    Notes
    -----
    Requires the `mpmath <http://mpmath.org/>`_ library.  See
    `~astropy.stats.scipy_poisson_upper_limit` for an implementation
    that is based on scipy and evaluates faster, but runs only to about
    N = 100.
    r   )mpfrs   findrootfsumrw   rt   rr   g-C6?c                s8    fddt t| d D }d  |  S )Nc                s   g | ]} || qS r   r   )r   rA   )ry   rs   rw   r   r   r     s    z>_mpmath_kraft_burrows_nousek.<locals>.eqn8.<locals>.<listcomp>r   g      ?)r+   r9   )rx   ry   Zsumterms)rt   rs   r   rw   )ry   r   rz     s    "z*_mpmath_kraft_burrows_nousek.<locals>.eqn8c                s"   | | } | ||    S )Nr   )r{   rx   ry   r|   )r}   rt   r~   r   r   r     s    z*_mpmath_kraft_burrows_nousek.<locals>.eqn7c                s    fdd}|| |gS )Nc                s   |  S )Nr   )r{   )ry   rx   r   r   r   eqn7NB  s    z?_mpmath_kraft_burrows_nousek.<locals>.eqn9_left.<locals>.eqn7NBr   )r   r   rx   ry   r   )r   rr   )ry   rx   r   r     s    z/_mpmath_kraft_burrows_nousek.<locals>.eqn9_leftc                sV   |   ks$d kr(dS  fdd}|d  gddS dS )a  
        Kraft, Burrows and Nousek suggest to integrate from N-B in both
        directions at once, so that S_min and S_max move similarly (see
        the article for details). Here, this is implemented differently:
        Treat S_max as the optimization parameters in func and then
        calculate the matching s_min that has has eqn7(S_max) =
        eqn7(S_min) here.
        r   g        c                s   |   S )Nr   )rL   )ry   rx   r   r   r   r   	eqn7ysmax  s    zC_mpmath_kraft_burrows_nousek.<locals>.find_s_min.<locals>.eqn7ysmaxridder)solvertolNr   )r   rx   ry   r   )r   r   r   )ry   rx   r   r   r     s    	z0_mpmath_kraft_burrows_nousek.<locals>.find_s_minc                s"   |  }||  }| S )Nr   )r6   r   r   )ry   r   rx   r   r   r   r   ra     s    z*_mpmath_kraft_burrows_nousek.<locals>.funcg      ?r   r   )r   r   )
Zmpmathr   rs   r   r   rw   rt   rr   r>   max)	rx   ry   r   r   rz   ra   Zs_max_guessr   r   r   )ry   r   rx   r   r}   r   rt   rs   r~   r   r   r   rw   rr   r   r   _mpmath_kraft_burrows_nousek  s&    $
	r   c             C   sf   ddl m}m} |rJ| dkrJyt| ||S  tk
rH   |sDtdY nX |rZt| ||S tddS )a=  Upper limit on a poisson count rate

    The implementation is based on Kraft, Burrows and Nousek in
    `ApJ 374, 344 (1991) <https://ui.adsabs.harvard.edu/abs/1991ApJ...374..344K>`_.
    The XMM-Newton upper limit server used the same formalism.

    Parameters
    ----------
    N : int or np.int32/np.int64
        Total observed count number
    B : float or np.float32/np.float64
        Background count rate (assumed to be known with negligible error
        from a large background area).
    CL : float or np.float32/np.float64
       Confidence level (number between 0 and 1)

    Returns
    -------
    S : source count limit

    Notes
    -----
    This functions has an optional dependency: Either :mod:`scipy` or `mpmath
    <http://mpmath.org/>`_  need to be available. (Scipy only works for
    N < 100).
    r   )	HAS_SCIPY
HAS_MPMATHrh   z1Need mpmath package for input numbers this large.z$Either scipy or mpmath are required.N)Z"astropy.utils.compat.optional_depsr   r   r   OverflowErrorr7   r   ImportError)rx   ry   r   r   r   r   r   r   r[     s    r[   c             C   s  yddl m}m} W n$ tk
r8   ddlm}m} Y nX | dk sJ| dkrRtd| d| k r~d||| d|  |d    S | d| k r"||  d  d }t|d	 ||  d d	 d  }| | | |  }}d||d ||d  d|  ||d  d|    ||d	   ||   S | d
kr:|d	 dks^| |d d|  kr$|d	 dkr$tt	|d|   d }| ||  }	|	|d  |	d | |	d	 | dd	|    |	| |d  dd	|   |  ||d  |d	  |d	    }
|
||| d|  ||  || d   }|
 S | t| }tdd| }d	d|d	  |d	  d  td|d	  |d	   
 }|d	 d|d	  |d	  d  td|d	  |d	   
 }|d|  d |  S dS )a  Compute the false positive probability for the Kuiper statistic.

    Uses the set of four formulas described in Paltani 2004; they report
    the resulting function never underestimates the false positive
    probability but can be a bit high in the N=40..50 range.
    (They quote a factor 1.5 at the 1e-7 level.)

    Parameters
    ----------
    D : float
        The Kuiper test score.
    N : float
        The effective sample size.

    Returns
    -------
    fpp : float
        The probability of a score this large arising from the null hypothesis.

    Notes
    -----
    Eq 7 of Paltani 2004 appears to incorrectly quote the original formula
    (Stephens 1965). This function implements the original formula, as it
    produces a result closer to Monte Carlo simulations.

    References
    ----------

    .. [1] Paltani, S., "Searching for periods in X-ray observations using
           Kuiper's test. Application to the ROSAT PSPC archive",
           Astronomy and Astrophysics, v.240, p.789-790, 2004.

    .. [2] Stephens, M. A., "The goodness-of-fit statistic VN: distribution
           and significance points", Biometrika, v.52, p.309, 1965.

    r   )rs   combg        g       @z2Must have 0<=D<=2 by definition of the Kuiper testg      ?r   g      @r3   g      ?   gR2@r4      N)r;   rs   r   r   Z
scipy.miscr7   r   r<   ru   floorrv   rt   )Drx   rs   r   r@   rabrd   yZTttermzmsZS1ZS2r   r   r   r     s4    % "J<
b*:>c             C   s   | S )Nr   )rL   r   r   r   r   Z  r   r   r   c             C   sj   t | } || f| }t| }t |t |t|  t t |d t| |  }|t||fS )a  Compute the Kuiper statistic.

    Use the Kuiper statistic version of the Kolmogorov-Smirnov test to
    find the probability that a sample like ``data`` was drawn from the
    distribution whose CDF is given as ``cdf``.

    .. warning::
        This will not work correctly for distributions that are actually
        discrete (Poisson, for example).

    Parameters
    ----------
    data : array-like
        The data values.
    cdf : callable
        A callable to evaluate the CDF of the distribution being tested
        against. Will be called with a vector of all values at once.
        The default is a uniform distribution.
    args : list-like, optional
        Additional arguments to be supplied to cdf.

    Returns
    -------
    D : float
        The raw statistic.
    fpp : float
        The probability of a D this large arising with a sample drawn from
        the distribution whose CDF is cdf.

    Notes
    -----
    The Kuiper statistic resembles the Kolmogorov-Smirnov test in that
    it is nonparametric and invariant under reparameterizations of the data.
    The Kuiper statistic, in addition, is equally sensitive throughout
    the domain, and it is also invariant under cyclic permutations (making
    it particularly appropriate for analyzing circular data).

    Returns (D, fpp), where D is the Kuiper D number and fpp is the
    probability that a value as large as D would occur if data was
    drawn from cdf.

    .. warning::
        The fpp is calculated only approximately, and it can be
        as much as 1.5 times the true value.

    Stephens 1970 claims this is more effective than the KS at detecting
    changes in the variance of a distribution; the KS is (he claims) more
    sensitive at detecting changes in the mean.

    If cdf was obtained from data by fitting, then fpp is not correct and
    it will be necessary to do Monte Carlo simulations to interpret D.
    D should normally be independent of the shape of CDF.

    References
    ----------

    .. [1] Stephens, M. A., "Use of the Kolmogorov-Smirnov, Cramer-Von Mises
           and Related Statistics Without Extensive Tables", Journal of the
           Royal Statistical Society. Series B (Methodological), Vol. 32,
           No. 1. (1970), pp. 115-122.


    r   )r   sortr%   Zamaxru   r>   r   )r,   Zcdfr   Zcdfvrx   r   r   r   r   r   Z  s    A
"c             C   s   t | } t |}| j\}|j\}t g | j|jg}t |t jrTt |t jr\tdt 	| d sxt 	|d rtdt
t | |t ||}t| t| tt| t|  }|t||fS )a  Compute the Kuiper statistic to compare two samples.

    Parameters
    ----------
    data1 : array-like
        The first set of data values.
    data2 : array-like
        The second set of data values.

    Returns
    -------
    D : float
        The raw test statistic.
    fpp : float
        The probability of obtaining two samples this different from
        the same distribution.

    .. warning::
        The fpp is quite approximate, especially for small samples.

    z#kuiper_two only accepts real inputsrI   z&kuiper_two only accepts non-nan inputs)r   r   r*   Zfind_common_typerW   Z
issubdtypenumberZcomplexfloatingr7   r^   r   Zks_2sampr!   r%   r>   r   )Zdata1Zdata2Zn1Zn2Zcommon_typer   ZNer   r   r   r     s    

$c             C   s  g }t  }d}xv| D ]n\}}}|t|t| | 7 }|d }|| |d|| f |d }|| ||d| f qW |d |d t|}tdd t|D }	t	t
|d }
|
|7 }
x.|D ]&\}}}|
|	| |	|   |7  < qW t||
fS )a  Fold the weighted intervals to the interval (0,1).

    Convert a list of intervals (ai, bi, wi) to a list of non-overlapping
    intervals covering (0,1). Each output interval has a weight equal
    to the sum of the wis of all the intervals that include it. All intervals
    are interpreted modulo 1, and weights are accumulated counting
    multiplicity. This is appropriate, for example, if you have one or more
    blocks of observation and you want to determine how much observation
    time was spent on different parts of a system's orbit (the blocks
    should be converted to units of the orbital period first).

    Parameters
    ----------
    intervals : list of (3,) tuple
        For each tuple (ai,bi,wi); ai and bi are the limits of the interval,
        and wi is the weight to apply to the interval.

    Returns
    -------
    breaks : (N,) array of float
        The endpoints of a set of intervals covering [0,1]; breaks[0]=0 and
        breaks[-1] = 1
    weights : (N-1,) array of float
        The ith element is the sum of number of times the interval
        breaks[i],breaks[i+1] is included in each interval times the weight
        associated with that interval.

    r   r   g        g      ?c             S   s   g | ]\}}||fqS r   r   )r   rp   fr   r   r   r     s    z"fold_intervals.<locals>.<listcomp>)setr   ceilr   addappendsorteddict	enumeratezerosr%   r?   )Z	intervalsr   breaksZtotr   r   wtfaZfbZ
breaks_maptotalsr   r   r   r     s(    



 c          	      s   | d dks| d dkr t dtt| dkr<t dt|dk rRt dt|dkrht d|   tdt|t  fd   fd	d
S )aN  Construct a callable piecewise-linear CDF from a pair of arrays.

    Take a pair of arrays in the format returned by fold_intervals and
    make a callable cumulative distribution function on the interval
    (0,1).

    Parameters
    ----------
    breaks : (N,) array of float
        The boundaries of successive intervals.
    totals : (N-1,) array of float
        The weight for each interval.

    Returns
    -------
    f : callable
        A cumulative distribution function corresponding to the
        piecewise-constant probability distribution given by breaks, weights

    r   rI   r   z%Intervals must be restricted to [0,1]z"Breaks must be strictly increasingz5Total weights in each subinterval must be nonnegativez1At least one interval must have positive exposure)r   c                s   t |  ddS )Nr   r   )r   Zinterp)rL   )r   cr   r   r      r   z$cdf_from_intervals.<locals>.<lambda>)r7   r   r:   diffallr\   ZconcatenateZcumsum)r   r   r   )r   r   r   r     s    c             C   sh   | \}}|\}}||k r>||k r$dS ||k r4|| S || S n&||k r`||k rV|| S || S ndS dS )a	  Compute the length of overlap of two intervals.

    Parameters
    ----------
    i1, i2 : (float, float)
        The two intervals, (interval 1, interval 2).

    Returns
    -------
    l : float
        The length of the overlap between the two intervals.

    g        r   Nr   )i1i2r   r   r   dr   r   r   r   #  s    

c       	      C   s   t | }|d }x|tt|D ]l}||d  }xVt| D ]J}tt||  t|d |  f||f}||  |d|   ||  7  < q:W |}q W |S )a  Histogram of a piecewise-constant weight function.

    This function takes a piecewise-constant weight function and
    computes the average weight in each histogram bin.

    Parameters
    ----------
    n : int
        The number of bins
    breaks : (N,) array of float
        Endpoints of the intervals in the PDF
    totals : (N-1,) array of float
        Probability densities in each bin

    Returns
    -------
    h : array of float
        The average weight for each bin

    r   r   g      ?)r   r   r+   r%   r   r>   )	rA   r   r   hstartrp   endjolr   r   r   r   C  s    
$)r0   r1   )rG   Nr0   r1   )rR   r   r   N)NNF)NNF)r   )rh   NN)"__doc__r   numpyr   Zastropy.utils.decoratorsr    r   __all__Z__doctest_skip__Z__doctest_requires__r<   logr   r   r-   r   r   rQ   r	   r
   r   r   r   r   r   r[   r   r   r   r   r   r   r   r   r   r   r   <module>
   sV   

-
 d
  4	
   
X
=
*
qO^*LJ'4$ 