# Author:  Travis Oliphant, 2002
#
# Further updates and enhancements by many SciPy developers.
#
from __future__ import division, print_function, absolute_import

import math
import types
import warnings

from . import statlib
from . import stats
from .stats import find_repeats
from . import distributions
from numpy import isscalar, r_, log, sum, around, unique, asarray
from numpy import zeros, arange, sort, amin, amax, any, where, \
     atleast_1d, sqrt, ceil, floor, array, poly1d, compress, not_equal, \
     pi, exp, ravel, angle
import numpy as np
import scipy.optimize as optimize
from numpy.testing.decorators import setastest


__all__ = ['mvsdist',
           'bayes_mvs', 'kstat', 'kstatvar', 'probplot', 'ppcc_max', 'ppcc_plot',
           'boxcox_llf', 'boxcox', 'boxcox_normmax', 'boxcox_normplot',
           'shapiro', 'anderson', 'ansari', 'bartlett', 'levene', 'binom_test',
           'fligner', 'mood', 'oneway', 'wilcoxon',
           'pdf_fromgamma', 'circmean', 'circvar', 'circstd',
          ]


def bayes_mvs(data, alpha=0.90):
    """
    Bayesian confidence intervals for the mean, var, and std.

    Parameters
    ----------
    data : array_like
        Input data, if multi-dimensional it is flattened to 1-D by `bayes_mvs`.
        Requires 2 or more data points.
    alpha : float, optional
        Probability that the returned confidence interval contains
        the true parameter.

    Returns
    -------
    mean_cntr, var_cntr, std_cntr : tuple
        The three results are for the mean, variance and standard deviation,
        respectively.  Each result is a tuple of the form::

            (center, (lower, upper))

        with `center` the mean of the conditional pdf of the value given the
        data, and `(lower, upper)` a confidence interval, centered on the
        median, containing the estimate to a probability `alpha`.

    Notes
    -----
    Each tuple of mean, variance, and standard deviation estimates represent
    the (center, (lower, upper)) with center the mean of the conditional pdf
    of the value given the data and (lower, upper) is a confidence interval
    centered on the median, containing the estimate to a probability
    `alpha`.

    Converts data to 1-D and assumes all data has the same mean and variance.
    Uses Jeffrey's prior for variance and std.

    Equivalent to tuple((x.mean(), x.interval(alpha)) for x in mvsdist(dat))

    References
    ----------
    T.E. Oliphant, "A Bayesian perspective on estimating mean, variance, and
    standard-deviation from data", http://hdl.handle.net/1877/438, 2006.

    """
    res = mvsdist(data)
    if alpha >= 1 or alpha <= 0:
        raise ValueError("0 < alpha < 1 is required, but alpha=%s was given." % alpha)
    return tuple((x.mean(), x.interval(alpha)) for x in res)

def mvsdist(data):
    """
    'Frozen' distributions for mean, variance, and standard deviation of data.

    Parameters
    ----------
    data : array_like
        Input array. Converted to 1-D using ravel.
        Requires 2 or more data-points.

    Returns
    -------
    mdist : "frozen" distribution object
        Distribution object representing the mean of the data
    vdist : "frozen" distribution object
        Distribution object representing the variance of the data
    sdist : "frozen" distribution object
        Distribution object representing the standard deviation of the data

    Notes
    -----
    The return values from bayes_mvs(data) is equivalent to
    ``tuple((x.mean(), x.interval(0.90)) for x in mvsdist(data))``.

    In other words, calling ``<dist>.mean()`` and ``<dist>.interval(0.90)``
    on the three distribution objects returned from this function will give
    the same results that are returned from `bayes_mvs`.

    Examples
    --------
    >>> from scipy.stats import mvsdist
    >>> data = [6, 9, 12, 7, 8, 8, 13]
    >>> mean, var, std = mvsdist(data)

    We now have frozen distribution objects "mean", "var" and "std" that we can
    examine:

    >>> mean.mean()
    9.0
    >>> mean.interval(0.95)
    (6.6120585482655692, 11.387941451734431)
    >>> mean.std()
    1.1952286093343936

    """
    x = ravel(data)
    n = len(x)
    if (n < 2):
        raise ValueError("Need at least 2 data-points.")
    xbar = x.mean()
    C = x.var()
    if (n > 1000): # gaussian approximations for large n
        mdist = distributions.norm(loc=xbar, scale=math.sqrt(C/n))
        sdist = distributions.norm(loc=math.sqrt(C), scale=math.sqrt(C/(2.*n)))
        vdist = distributions.norm(loc=C, scale=math.sqrt(2.0/n)*C)
    else:
        nm1 = n-1
        fac = n*C/2.
        val = nm1/2.
        mdist = distributions.t(nm1,loc=xbar,scale=math.sqrt(C/nm1))
        sdist = distributions.gengamma(val,-2,scale=math.sqrt(fac))
        vdist = distributions.invgamma(val,scale=fac)
    return mdist, vdist, sdist


def kstat(data,n=2):
    """
    Return the nth k-statistic (1<=n<=4 so far).

    The nth k-statistic is the unique symmetric unbiased estimator of the nth
    cumulant kappa_n.

    Parameters
    ----------
    data : array_like
        Input array.
    n : int, {1, 2, 3, 4}, optional
        Default is equal to 2.

    Returns
    -------
    kstat : float
        The nth k-statistic.

    See Also
    --------
    kstatvar: Returns an unbiased estimator of the variance of the k-statistic.

    Notes
    -----
    The cumulants are related to central moments but are specifically defined
    using a power series expansion of the logarithm of the characteristic
    function (which is the Fourier transform of the PDF).
    In particular let phi(t) be the characteristic function, then::

        ln phi(t) = > kappa_n (it)^n / n!    (sum from n=0 to inf)

    The first few cumulants (kappa_n)  in terms of central moments (mu_n) are::

        kappa_1 = mu_1
        kappa_2 = mu_2
        kappa_3 = mu_3
        kappa_4 = mu_4 - 3*mu_2**2
        kappa_5 = mu_5 - 10*mu_2 * mu_3

    References
    -------
---
    http://mathworld.wolfram.com/k-Statistic.html


    http://mathworld.wolfram.com/Cumulant.html

    """
    if n > 4 or n < 1:
        raise ValueError("k-statistics only supported for 1<=n<=4")
    n = int(n)
    S = zeros(n+1,'d')
    data = ravel(data)
    N = len(data)
    for k in range(1,n+1):
        S[k] = sum(data**k,axis=0)
    if n==1:
        return S[1]*1.0/N
    elif n==2:
        return (N*S[2]-S[1]**2.0)/(N*(N-1.0))
    elif n==3:
        return (2*S[1]**3 - 3*N*S[1]*S[2]+N*N*S[3]) / (N*(N-1.0)*(N-2.0))
    elif n==4:
        return (-6*S[1]**4 + 12*N*S[1]**2 * S[2] - 3*N*(N-1.0)*S[2]**2 - \
                4*N*(N+1)*S[1]*S[3] + N*N*(N+1)*S[4]) / \
                (N*(N-1.0)*(N-2.0)*(N-3.0))
    else:
        raise ValueError("Should not be here.")

def kstatvar(data,n=2):
    """
    Returns an unbiased estimator of the variance of the k-statistic.

    See `kstat` for more details of the k-statistic.

    Parameters
    ----------
    data : array_like
        Input array.
    n : int, {1, 2}, optional
        Default is equal to 2.

    Returns
    -------
    kstatvar : float
        The nth k-statistic variance.

    See Also
    --------
    kstat

    """
    data = ravel(data)
    N = len(data)
    if n == 1:
        return kstat(data,n=2)*1.0/N
    elif n == 2:
        k2 = kstat(data,n=2)
        k4 = kstat(data,n=4)
        return (2*k2*k2*N + (N-1)*k4)/(N*(N+1))
    else:
        raise ValueError("Only n=1 or n=2 supported.")


def probplot(x, sparams=(), dist='norm', fit=True, plot=None):
    """
    Calculate quantiles for a probability plot of sample data against a
    specified theoretical distribution.

    `probplot` optionally calculates a best-fit line for the data and plots the
    results using Matplotlib or a given plot function.

    Parameters
    ----------
    x : array_like
        Sample/response data from which `probplot` creates the plot.
    sparams : tuple, optional
        Distribution-specific shape parameters (location(s) and scale(s)).
    dist : str, optional
        Distribution function name. The default is 'norm' for a normal
        probability plot.
    fit : bool, optional
   