# stats.morestats¶

## Module Contents¶

### Functions¶

 bayes_mvs(data,alpha=0.9) r mvsdist(data) ‘Frozen’ distributions for mean, variance, and standard deviation of data. kstat(data,n=2) r kstatvar(data,n=2) r _calc_uniform_order_statistic_medians(n) Approximations of uniform order statistic medians. _parse_dist_kw(dist,enforce_subclass=True) Parse dist keyword. _add_axis_labels_title(plot,xlabel,ylabel,title) Helper function to add axes labels and a title to stats plots probplot(x,sparams=tuple,dist=”norm”,fit=True,plot=None,rvalue=False) Calculate quantiles for a probability plot, and optionally show the plot. ppcc_max(x,brack=tuple,dist=”tukeylambda”) Calculate the shape parameter that maximizes the PPCC ppcc_plot(x,a,b,dist=”tukeylambda”,plot=None,N=80) Calculate and optionally plot probability plot correlation coefficient. boxcox_llf(lmb,data) rThe boxcox log-likelihood function. _boxcox_conf_interval(x,lmax,alpha) boxcox(x,lmbda=None,alpha=None) r boxcox_normmax(x,brack=tuple,method=”pearsonr”) Compute optimal Box-Cox transform parameter for input data. boxcox_normplot(x,la,lb,plot=None,N=80) Compute parameters for a Box-Cox normality plot, optionally show it. shapiro(x) Perform the Shapiro-Wilk test for normality. anderson(x,dist=”norm”) Anderson-Darling test for data coming from a particular distribution _anderson_ksamp_midrank(samples,Z,Zstar,k,n,N) Compute A2akN equation 7 of Scholz and Stephens. _anderson_ksamp_right(samples,Z,Zstar,k,n,N) Compute A2akN equation 6 of Scholz & Stephens. anderson_ksamp(samples,midrank=True) The Anderson-Darling test for k-samples. ansari(x,y) Perform the Ansari-Bradley test for equal scale parameters bartlett(*args) Perform Bartlett’s test for equal variances levene(*args,**kwds) Perform Levene test for equal variances. binom_test(x,n=None,p=0.5,alternative=”two-sided”) Perform a test that the probability of success is p. _apply_func(x,g,func) fligner(*args,**kwds) Perform Fligner-Killeen test for equality of variance. mood(x,y,axis=0) Perform Mood’s test for equal scale parameters. wilcoxon(x,y=None,zero_method=”wilcox”,correction=False) Calculate the Wilcoxon signed-rank test. median_test(*args,**kwds) Mood’s median test. _circfuncs_common(samples,high,low) circmean(samples,high=None,low=0,axis=None) Compute the circular mean for samples in a range. circvar(samples,high=None,low=0,axis=None) Compute the circular variance for samples assumed to be in a range circstd(samples,high=None,low=0,axis=None) Compute the circular standard deviation for samples assumed to be in the
bayes_mvs(data, alpha=0.9)

r Bayesian confidence intervals for the mean, var, and std.

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.
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.

mvsdist

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))

T.E. Oliphant, “A Bayesian perspective on estimating mean, variance, and standard-deviation from data”, http://scholarsarchive.byu.edu/facpub/278, 2006.

First a basic example to demonstrate the outputs:

>>> from scipy import stats
>>> data = [6, 9, 12, 7, 8, 8, 13]
>>> mean, var, std = stats.bayes_mvs(data)
>>> mean
Mean(statistic=9.0, minmax=(7.1036502226125329, 10.896349777387467))
>>> var
Variance(statistic=10.0, minmax=(3.176724206..., 24.45910382...))
>>> std
Std_dev(statistic=2.9724954732045084, minmax=(1.7823367265645143, 4.9456146050146295))


Now we generate some normally distributed random data, and get estimates of mean and standard deviation with 95% confidence intervals for those estimates:

>>> n_samples = 100000
>>> data = stats.norm.rvs(size=n_samples)
>>> res_mean, res_var, res_std = stats.bayes_mvs(data, alpha=0.95)

>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.hist(data, bins=100, normed=True, label='Histogram of data')
>>> ax.vlines(res_mean.statistic, 0, 0.5, colors='r', label='Estimated mean')
>>> ax.axvspan(res_mean.minmax,res_mean.minmax, facecolor='r',
...            alpha=0.2, label=r'Estimated mean (95% limits)')
>>> ax.vlines(res_std.statistic, 0, 0.5, colors='g', label='Estimated scale')
>>> ax.axvspan(res_std.minmax,res_std.minmax, facecolor='g', alpha=0.2,
...            label=r'Estimated scale (95% limits)')

>>> ax.legend(fontsize=10)
>>> ax.set_xlim([-4, 4])
>>> ax.set_ylim([0, 0.5])
>>> plt.show()

mvsdist(data)

‘Frozen’ distributions for mean, variance, and standard deviation of data.

data : array_like
Input array. Converted to 1-D using ravel. Requires 2 or more data-points.
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

bayes_mvs

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.

T.E. Oliphant, “A Bayesian perspective on estimating mean, variance, and standard-deviation from data”, http://scholarsarchive.byu.edu/facpub/278, 2006.

>>> from scipy import stats
>>> data = [6, 9, 12, 7, 8, 8, 13]
>>> mean, var, std = stats.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

kstat(data, n=2)

r Return the nth k-statistic (1<=n<=4 so far).

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

data : array_like
Input array. Note that n-D input gets flattened.
n : int, {1, 2, 3, 4}, optional
Default is equal to 2.
kstat : float
The nth k-statistic.

kstatvar: Returns an unbiased estimator of the variance of the k-statistic. moment: Returns the n-th central moment about the mean for a sample.

For a sample size n, the first few k-statistics are given by:

where is the sample mean, is the sample variance, and is the i-th sample central moment.

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

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

>>> from scipy import stats
>>> rndm = np.random.RandomState(1234)


As sample size increases, n-th moment and n-th k-statistic converge to the same number (although they aren’t identical). In the case of the normal distribution, they converge to zero.

>>> for n in [2, 3, 4, 5, 6, 7]:
...     x = rndm.normal(size=10**n)
...     m, k = stats.moment(x, 3), stats.kstat(x, 3)
...     print("%.3g %.3g %.3g" % (m, k, m-k))
-0.631 -0.651 0.0194
0.0282 0.0283 -8.49e-05
-0.0454 -0.0454 1.36e-05
7.53e-05 7.53e-05 -2.26e-09
0.00166 0.00166 -4.99e-09
-2.88e-06 -2.88e-06 8.63e-13

kstatvar(data, n=2)

r Returns an unbiased estimator of the variance of the k-statistic.

See kstat for more details of the k-statistic.

data : array_like
Input array. Note that n-D input gets flattened.
n : int, {1, 2}, optional
Default is equal to 2.
kstatvar : float
The nth k-statistic variance.

kstat: Returns the n-th k-statistic. moment: Returns the n-th central moment about the mean for a sample.

The variances of the first few k-statistics are given by:

_calc_uniform_order_statistic_medians(n)

Approximations of uniform order statistic medians.

n : int
Sample size.
v : 1d float array
Approximations of the order statistic medians.
  James J. Filliben, “The Probability Plot Correlation Coefficient Test for Normality”, Technometrics, Vol. 17, pp. 111-117, 1975.

Order statistics of the uniform distribution on the unit interval are marginally distributed according to beta distributions. The expectations of these order statistic are evenly spaced across the interval, but the distributions are skewed in a way that pushes the medians slightly towards the endpoints of the unit interval:

>>> n = 4
>>> k = np.arange(1, n+1)
>>> from scipy.stats import beta
>>> a = k
>>> b = n-k+1
>>> beta.mean(a, b)
array([ 0.2,  0.4,  0.6,  0.8])
>>> beta.median(a, b)
array([ 0.15910358,  0.38572757,  0.61427243,  0.84089642])


The Filliben approximation uses the exact medians of the smallest and greatest order statistics, and the remaining medians are approximated by points spread evenly across a sub-interval of the unit interval:

>>> from scipy.morestats import _calc_uniform_order_statistic_medians
>>> _calc_uniform_order_statistic_medians(n)
array([ 0.15910358,  0.38545246,  0.61454754,  0.84089642])


This plot shows the skewed distributions of the order statistics of a sample of size four from a uniform distribution on the unit interval:

>>> import matplotlib.pyplot as plt
>>> x = np.linspace(0.0, 1.0, num=50, endpoint=True)
>>> pdfs = [beta.pdf(x, a[i], b[i]) for i in range(n)]
>>> plt.figure()
>>> plt.plot(x, pdfs, x, pdfs, x, pdfs, x, pdfs)

_parse_dist_kw(dist, enforce_subclass=True)

Parse dist keyword.

dist : str or stats.distributions instance.
Several functions take dist as a keyword, hence this utility function.
enforce_subclass : bool, optional
If True (default), dist needs to be a _distn_infrastructure.rv_generic instance. It can sometimes be useful to set this keyword to False, if a function wants to accept objects that just look somewhat like such an instance (for example, they have a ppf method).
_add_axis_labels_title(plot, xlabel, ylabel, title)

Helper function to add axes labels and a title to stats plots

probplot(x, sparams=tuple, dist="norm", fit=True, plot=None, rvalue=False)

Calculate quantiles for a probability plot, and optionally show the plot.

Generates a probability plot of sample data against the quantiles of a specified theoretical distribution (the normal distribution by default). probplot optionally calculates a best-fit line for the data and plots the results using Matplotlib or a given plot function.

x : array_like
Sample/response data from which probplot creates the plot.
sparams : tuple, optional
Distribution-specific shape parameters (shape parameters plus location and scale).
dist : str or stats.distributions instance, optional
Distribution or distribution function name. The default is ‘norm’ for a normal probability plot. Objects that look enough like a stats.distributions instance (i.e. they have a ppf method) are also accepted.
fit : bool, optional
Fit a least-squares regression (best-fit) line to the sample data if True (default).
plot : object, optional
If given, plots the quantiles and least squares fit. plot is an object that has to have methods “plot” and “text”. The matplotlib.pyplot module or a Matplotlib Axes object can be used, or a custom object with the same methods. Default is None, which means that no plot is created.
(osm, osr) : tuple of ndarrays
Tuple of theoretical quantiles (osm, or order statistic medians) and ordered responses (osr). osr is simply sorted input x. For details on how osm is calculated see the Notes section.
(slope, intercept, r) : tuple of floats, optional
Tuple containing the result of the least-squares fit, if that is performed by probplot. r is the square root of the coefficient of determination. If fit=False and plot=None, this tuple is not returned.

Even if plot is given, the figure is not shown or saved by probplot; plt.show() or plt.savefig('figname.png') should be used after calling probplot.

probplot generates a probability plot, which should not be confused with a Q-Q or a P-P plot. Statsmodels has more extensive functionality of this type, see statsmodels.api.ProbPlot.

The formula used for the theoretical quantiles (horizontal axis of the probability plot) is Filliben’s estimate:

quantiles = dist.ppf(val), for

0.5**(1/n),                  for i = n
val = (i - 0.3175) / (n + 0.365),  for i = 2, ..., n-1
1 - 0.5**(1/n),              for i = 1


where i indicates the i-th ordered value and n is the total number of values.

>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> nsample = 100
>>> np.random.seed(7654321)


A t distribution with small degrees of freedom:

>>> ax1 = plt.subplot(221)
>>> x = stats.t.rvs(3, size=nsample)
>>> res = stats.probplot(x, plot=plt)


A t distribution with larger degrees of freedom:

>>> ax2 = plt.subplot(222)
>>> x = stats.t.rvs(25, size=nsample)
>>> res = stats.probplot(x, plot=plt)


A mixture of two normal distributions with broadcasting:

>>> ax3 = plt.subplot(223)
>>> x = stats.norm.rvs(loc=[0,5], scale=[1,1.5],
...                    size=(nsample//2,2)).ravel()
>>> res = stats.probplot(x, plot=plt)


A standard normal distribution:

>>> ax4 = plt.subplot(224)
>>> x = stats.norm.rvs(loc=0, scale=1, size=nsample)
>>> res = stats.probplot(x, plot=plt)


Produce a new figure with a loggamma distribution, using the dist and sparams keywords:

>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> x = stats.loggamma.rvs(c=2.5, size=500)
>>> res = stats.probplot(x, dist=stats.loggamma, sparams=(2.5,), plot=ax)
>>> ax.set_title("Probplot for loggamma dist with shape parameter 2.5")


Show the results with Matplotlib:

>>> plt.show()

ppcc_max(x, brack=tuple, dist="tukeylambda")

Calculate the shape parameter that maximizes the PPCC

The probability plot correlation coefficient (PPCC) plot can be used to determine the optimal shape parameter for a one-parameter family of distributions. ppcc_max returns the shape parameter that would maximize the probability plot correlation coefficient for the given data to a one-parameter family of distributions.

x : array_like
Input array.
brack : tuple, optional
Triple (a,b,c) where (a<b<c). If bracket consists of two numbers (a, c) then they are assumed to be a starting interval for a downhill bracket search (see scipy.optimize.brent).
dist : str or stats.distributions instance, optional
Distribution or distribution function name. Objects that look enough like a stats.distributions instance (i.e. they have a ppf method) are also accepted. The default is 'tukeylambda'.
shape_value : float
The shape parameter at which the probability plot correlation coefficient reaches its max value.

ppcc_plot, probplot, boxcox

The brack keyword serves as a starting point which is useful in corner cases. One can use a plot to obtain a rough visual estimate of the location for the maximum to start the search near it.

  J.J. Filliben, “The Probability Plot Correlation Coefficient Test for Normality”, Technometrics, Vol. 17, pp. 111-117, 1975.

First we generate some random data from a Tukey-Lambda distribution, with shape parameter -0.7:

>>> from scipy import stats
>>> x = stats.tukeylambda.rvs(-0.7, loc=2, scale=0.5, size=10000,
...                           random_state=1234567) + 1e4


Now we explore this data with a PPCC plot as well as the related probability plot and Box-Cox normplot. A red line is drawn where we expect the PPCC value to be maximal (at the shape parameter -0.7 used above):

>>> import matplotlib.pyplot as plt
>>> fig = plt.figure(figsize=(8, 6))
>>> ax = fig.add_subplot(111)
>>> res = stats.ppcc_plot(x, -5, 5, plot=ax)


We calculate the value where the shape should reach its maximum and a red line is drawn there. The line should coincide with the highest point in the ppcc_plot.

>>> max = stats.ppcc_max(x)
>>> ax.vlines(max, 0, 1, colors='r', label='Expected shape value')

>>> plt.show()

ppcc_plot(x, a, b, dist="tukeylambda", plot=None, N=80)

Calculate and optionally plot probability plot correlation coefficient.

The probability plot correlation coefficient (PPCC) plot can be used to determine the optimal shape parameter for a one-parameter family of distributions. It cannot be used for distributions without shape parameters (like the normal distribution) or with multiple shape parameters.

By default a Tukey-Lambda distribution (stats.tukeylambda) is used. A Tukey-Lambda PPCC plot interpolates from long-tailed to short-tailed distributions via an approximately normal one, and is therefore particularly useful in practice.

x : array_like
Input array.
a, b: scalar
Lower and upper bounds of the shape parameter to use.
dist : str or stats.distributions instance, optional
Distribution or distribution function name. Objects that look enough like a stats.distributions instance (i.e. they have a ppf method) are also accepted. The default is 'tukeylambda'.
plot : object, optional
If given, plots PPCC against the shape parameter. plot is an object that has to have methods “plot” and “text”. The matplotlib.pyplot module or a Matplotlib Axes object can be used, or a custom object with the same methods. Default is None, which means that no plot is created.
N : int, optional
Number of points on the horizontal axis (equally distributed from a to b).
svals : ndarray
The shape values for which ppcc was calculated.
ppcc : ndarray
The calculated probability plot correlation coefficient values.

ppcc_max, probplot, boxcox_normplot, tukeylambda

J.J. Filliben, “The Probability Plot Correlation Coefficient Test for Normality”, Technometrics, Vol. 17, pp. 111-117, 1975.

First we generate some random data from a Tukey-Lambda distribution, with shape parameter -0.7:

>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> np.random.seed(1234567)
>>> x = stats.tukeylambda.rvs(-0.7, loc=2, scale=0.5, size=10000) + 1e4


Now we explore this data with a PPCC plot as well as the related probability plot and Box-Cox normplot. A red line is drawn where we expect the PPCC value to be maximal (at the shape parameter -0.7 used above):

>>> fig = plt.figure(figsize=(12, 4))
>>> ax1 = fig.add_subplot(131)
>>> ax2 = fig.add_subplot(132)
>>> ax3 = fig.add_subplot(133)
>>> res = stats.probplot(x, plot=ax1)
>>> res = stats.boxcox_normplot(x, -5, 5, plot=ax2)
>>> res = stats.ppcc_plot(x, -5, 5, plot=ax3)
>>> ax3.vlines(-0.7, 0, 1, colors='r', label='Expected shape value')
>>> plt.show()

boxcox_llf(lmb, data)

rThe boxcox log-likelihood function.

lmb : scalar
Parameter for Box-Cox transformation. See boxcox for details.
data : array_like
Data to calculate Box-Cox log-likelihood for. If data is multi-dimensional, the log-likelihood is calculated along the first axis.
llf : float or ndarray
Box-Cox log-likelihood of data given lmb. A float for 1-D data, an array otherwise.

boxcox, probplot, boxcox_normplot, boxcox_normmax

The Box-Cox log-likelihood function is defined here as

where y is the Box-Cox transformed input data x.

>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.axes_grid1.inset_locator import inset_axes
>>> np.random.seed(1245)


Generate some random variates and calculate Box-Cox log-likelihood values for them for a range of lmbda values:

>>> x = stats.loggamma.rvs(5, loc=10, size=1000)
>>> lmbdas = np.linspace(-2, 10)
>>> llf = np.zeros(lmbdas.shape, dtype=float)
>>> for ii, lmbda in enumerate(lmbdas):
...     llf[ii] = stats.boxcox_llf(lmbda, x)


Also find the optimal lmbda value with boxcox:

>>> x_most_normal, lmbda_optimal = stats.boxcox(x)


Plot the log-likelihood as function of lmbda. Add the optimal lmbda as a horizontal line to check that that’s really the optimum:

>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(lmbdas, llf, 'b.-')
>>> ax.axhline(stats.boxcox_llf(lmbda_optimal, x), color='r')
>>> ax.set_xlabel('lmbda parameter')
>>> ax.set_ylabel('Box-Cox log-likelihood')


Now add some probability plots to show that where the log-likelihood is maximized the data transformed with boxcox looks closest to normal:

>>> locs = [3, 10, 4]  # 'lower left', 'center', 'lower right'
>>> for lmbda, loc in zip([-1, lmbda_optimal, 9], locs):
...     xt = stats.boxcox(x, lmbda=lmbda)
...     (osm, osr), (slope, intercept, r_sq) = stats.probplot(xt)
...     ax_inset = inset_axes(ax, width="20%", height="20%", loc=loc)
...     ax_inset.plot(osm, osr, 'c.', osm, slope*osm + intercept, 'k-')
...     ax_inset.set_xticklabels([])
...     ax_inset.set_yticklabels([])
...     ax_inset.set_title('$\lambda=%1.2f$' % lmbda)

>>> plt.show()

_boxcox_conf_interval(x, lmax, alpha)
boxcox(x, lmbda=None, alpha=None)

r Return a positive dataset transformed by a Box-Cox power transformation.

x : ndarray
Input array. Should be 1-dimensional.
lmbda : {None, scalar}, optional

If lmbda is not None, do the transformation for that value.

If lmbda is None, find the lambda that maximizes the log-likelihood function and return it as the second output argument.

alpha : {None, float}, optional
If alpha is not None, return the 100 * (1-alpha)% confidence interval for lmbda as the third output argument. Must be between 0.0 and 1.0.
boxcox : ndarray
Box-Cox power transformed array.
maxlog : float, optional
If the lmbda parameter is None, the second returned argument is the lambda that maximizes the log-likelihood function.
(min_ci, max_ci) : tuple of float, optional
If lmbda parameter is None and alpha is not None, this returned tuple of floats represents the minimum and maximum confidence limits given alpha.

probplot, boxcox_normplot, boxcox_normmax, boxcox_llf

The Box-Cox transform is given by:

y = (x**lmbda - 1) / lmbda,  for lmbda > 0
log(x),                  for lmbda = 0


boxcox requires the input data to be positive. Sometimes a Box-Cox transformation provides a shift parameter to achieve this; boxcox does not. Such a shift parameter is equivalent to adding a positive constant to x before calling boxcox.

The confidence limits returned when alpha is provided give the interval where:

with llf the log-likelihood function and the chi-squared function.

G.E.P. Box and D.R. Cox, “An Analysis of Transformations”, Journal of the Royal Statistical Society B, 26, 211-252 (1964).

>>> from scipy import stats
>>> import matplotlib.pyplot as plt


We generate some random variates from a non-normal distribution and make a probability plot for it, to show it is non-normal in the tails:

>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(211)
>>> x = stats.loggamma.rvs(5, size=500) + 5
>>> prob = stats.probplot(x, dist=stats.norm, plot=ax1)
>>> ax1.set_xlabel('')
>>> ax1.set_title('Probplot against normal distribution')


We now use boxcox to transform the data so it’s closest to normal:

>>> ax2 = fig.add_subplot(212)
>>> xt, _ = stats.boxcox(x)
>>> prob = stats.probplot(xt, dist=stats.norm, plot=ax2)
>>> ax2.set_title('Probplot after Box-Cox transformation')

>>> plt.show()

boxcox_normmax(x, brack=tuple, method="pearsonr")

Compute optimal Box-Cox transform parameter for input data.

x : array_like
Input array.
brack : 2-tuple, optional
The starting interval for a downhill bracket search with optimize.brent. Note that this is in most cases not critical; the final result is allowed to be outside this bracket.
method : str, optional

The method to determine the optimal transform parameter (boxcox lmbda parameter). Options are:

‘pearsonr’ (default)
Maximizes the Pearson correlation coefficient between y = boxcox(x) and the expected values for y if x would be normally-distributed.
‘mle’
Minimizes the log-likelihood boxcox_llf. This is the method used in boxcox.
‘all’
Use all optimization methods available, and return all results. Useful to compare different methods.
maxlog : float or ndarray
The optimal transform parameter found. An array instead of a scalar for method='all'.

boxcox, boxcox_llf, boxcox_normplot

>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> np.random.seed(1234)  # make this example reproducible


Generate some data and determine optimal lmbda in various ways:

>>> x = stats.loggamma.rvs(5, size=30) + 5
>>> y, lmax_mle = stats.boxcox(x)
>>> lmax_pearsonr = stats.boxcox_normmax(x)

>>> lmax_mle
7.177...
>>> lmax_pearsonr
7.916...
>>> stats.boxcox_normmax(x, method='all')
array([ 7.91667384,  7.17718692])

>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> prob = stats.boxcox_normplot(x, -10, 10, plot=ax)
>>> ax.axvline(lmax_mle, color='r')
>>> ax.axvline(lmax_pearsonr, color='g', ls='--')

>>> plt.show()

boxcox_normplot(x, la, lb, plot=None, N=80)

Compute parameters for a Box-Cox normality plot, optionally show it.

A Box-Cox normality plot shows graphically what the best transformation parameter is to use in boxcox to obtain a distribution that is close to normal.

x : array_like
Input array.
la, lb : scalar
The lower and upper bounds for the lmbda values to pass to boxcox for Box-Cox transformations. These are also the limits of the horizontal axis of the plot if that is generated.
plot : object, optional
If given, plots the quantiles and least squares fit. plot is an object that has to have methods “plot” and “text”. The matplotlib.pyplot module or a Matplotlib Axes object can be used, or a custom object with the same methods. Default is None, which means that no plot is created.
N : int, optional
Number of points on the horizontal axis (equally distributed from la to lb).
lmbdas : ndarray
The lmbda values for which a Box-Cox transform was done.
ppcc : ndarray
Probability Plot Correlelation Coefficient, as obtained from probplot when fitting the Box-Cox transformed input x against a normal distribution.

probplot, boxcox, boxcox_normmax, boxcox_llf, ppcc_max

Even if plot is given, the figure is not shown or saved by boxcox_normplot; plt.show() or plt.savefig('figname.png') should be used after calling probplot.

>>> from scipy import stats
>>> import matplotlib.pyplot as plt


Generate some non-normally distributed data, and create a Box-Cox plot:

>>> x = stats.loggamma.rvs(5, size=500) + 5
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> prob = stats.boxcox_normplot(x, -20, 20, plot=ax)


Determine and plot the optimal lmbda to transform x and plot it in the same plot:

>>> _, maxlog = stats.boxcox(x)
>>> ax.axvline(maxlog, color='r')

>>> plt.show()

shapiro(x)

Perform the Shapiro-Wilk test for normality.

The Shapiro-Wilk test tests the null hypothesis that the data was drawn from a normal distribution.

x : array_like
Array of sample data.
W : float
The test statistic.
p-value : float
The p-value for the hypothesis test.

anderson : The Anderson-Darling test for normality kstest : The Kolmogorov-Smirnov test for goodness of fit.

The algorithm used is described in _ but censoring parameters as described are not implemented. For N > 5000 the W test statistic is accurate but the p-value may not be.

The chance of rejecting the null hypothesis when it is true is close to 5% regardless of sample size.

  Shapiro, S. S. & Wilk, M.B (1965). An analysis of variance test for normality (complete samples), Biometrika, Vol. 52, pp. 591-611.
  Razali, N. M. & Wah, Y. B. (2011) Power comparisons of Shapiro-Wilk, Kolmogorov-Smirnov, Lilliefors and Anderson-Darling tests, Journal of Statistical Modeling and Analytics, Vol. 2, pp. 21-33.
  ALGORITHM AS R94 APPL. STATIST. (1995) VOL. 44, NO. 4.
>>> from scipy import stats
>>> np.random.seed(12345678)
>>> x = stats.norm.rvs(loc=5, scale=3, size=100)
>>> stats.shapiro(x)
(0.9772805571556091, 0.08144091814756393)

anderson(x, dist="norm")

Anderson-Darling test for data coming from a particular distribution

The Anderson-Darling test is a modification of the Kolmogorov- Smirnov test kstest for the null hypothesis that a sample is drawn from a population that follows a particular distribution. For the Anderson-Darling test, the critical values depend on which distribution is being tested against. This function works for normal, exponential, logistic, or Gumbel (Extreme Value Type I) distributions.

x : array_like
array of sample data
dist : {‘norm’,’expon’,’logistic’,’gumbel’,’gumbel_l’, gumbel_r’,
‘extreme1’}, optional the type of distribution to test against. The default is ‘norm’ and ‘extreme1’, ‘gumbel_l’ and ‘gumbel’ are synonyms.
statistic : float
The Anderson-Darling test statistic
critical_values : list
The critical values for this distribution
significance_level : list
The significance levels for the corresponding critical values in percents. The function returns critical values for a differing set of significance levels depending on the distribution that is being tested against.

Critical values provided are for the following significance levels:

normal/exponenential
15%, 10%, 5%, 2.5%, 1%
logistic
25%, 10%, 5%, 2.5%, 1%, 0.5%
Gumbel
25%, 10%, 5%, 2.5%, 1%

If A2 is larger than these critical values then for the corresponding significance level, the null hypothesis that the data come from the chosen distribution can be rejected.

  Stephens, M. A. (1974). EDF Statistics for Goodness of Fit and Some Comparisons, Journal of the American Statistical Association, Vol. 69, pp. 730-737.
  Stephens, M. A. (1976). Asymptotic Results for Goodness-of-Fit Statistics with Unknown Parameters, Annals of Statistics, Vol. 4, pp. 357-369.
  Stephens, M. A. (1977). Goodness of Fit for the Extreme Value Distribution, Biometrika, Vol. 64, pp. 583-588.
  Stephens, M. A. (1977). Goodness of Fit with Special Reference to Tests for Exponentiality , Technical Report No. 262, Department of Statistics, Stanford University, Stanford, CA.
  Stephens, M. A. (1979). Tests of Fit for the Logistic Distribution Based on the Empirical Distribution Function, Biometrika, Vol. 66, pp. 591-595.
_anderson_ksamp_midrank(samples, Z, Zstar, k, n, N)

Compute A2akN equation 7 of Scholz and Stephens.

samples : sequence of 1-D array_like
Array of sample arrays.
Z : array_like
Sorted array of all observations.
Zstar : array_like
Sorted array of unique observations.
k : int
Number of samples.
n : array_like
Number of observations in each sample.
N : int
Total number of observations.
A2aKN : float
The A2aKN statistics of Scholz and Stephens 1987.
_anderson_ksamp_right(samples, Z, Zstar, k, n, N)

Compute A2akN equation 6 of Scholz & Stephens.

samples : sequence of 1-D array_like
Array of sample arrays.
Z : array_like
Sorted array of all observations.
Zstar : array_like
Sorted array of unique observations.
k : int
Number of samples.
n : array_like
Number of observations in each sample.
N : int
Total number of observations.
A2KN : float
The A2KN statistics of Scholz and Stephens 1987.
anderson_ksamp(samples, midrank=True)

The Anderson-Darling test for k-samples.

The k-sample Anderson-Darling test is a modification of the one-sample Anderson-Darling test. It tests the null hypothesis that k-samples are drawn from the same population without having to specify the distribution function of that population. The critical values depend on the number of samples.

samples : sequence of 1-D array_like
Array of sample data in arrays.
midrank : bool, optional
Type of Anderson-Darling test which is computed. Default (True) is the midrank test applicable to continuous and discrete populations. If False, the right side empirical distribution is used.
statistic : float
Normalized k-sample Anderson-Darling test statistic.
critical_values : array
The critical values for significance levels 25%, 10%, 5%, 2.5%, 1%.
significance_level : float
An approximate significance level at which the null hypothesis for the provided samples can be rejected.
ValueError
If less than 2 samples are provided, a sample is empty, or no distinct observations are in the samples.

ks_2samp : 2 sample Kolmogorov-Smirnov test anderson : 1 sample Anderson-Darling test

_ Defines three versions of the k-sample Anderson-Darling test: one for continuous distributions and two for discrete distributions, in which ties between samples may occur. The default of this routine is to compute the version based on the midrank empirical distribution function. This test is applicable to continuous and discrete data. If midrank is set to False, the right side empirical distribution is used for a test for discrete data. According to _, the two discrete test statistics differ only slightly if a few collisions due to round-off errors occur in the test not adjusted for ties between samples.

New in version 0.14.0.

  Scholz, F. W and Stephens, M. A. (1987), K-Sample Anderson-Darling Tests, Journal of the American Statistical Association, Vol. 82, pp. 918-924.
>>> from scipy import stats
>>> np.random.seed(314159)


The null hypothesis that the two random samples come from the same distribution can be rejected at the 5% level because the returned test value is greater than the critical value for 5% (1.961) but not at the 2.5% level. The interpolation gives an approximate significance level of 3.1%:

>>> stats.anderson_ksamp([np.random.normal(size=50),
... np.random.normal(loc=0.5, size=30)])
(2.4615796189876105,
array([ 0.325,  1.226,  1.961,  2.718,  3.752]),
0.03134990135800783)


The null hypothesis cannot be rejected for three samples from an identical distribution. The approximate p-value (87%) has to be computed by extrapolation and may not be very accurate:

>>> stats.anderson_ksamp([np.random.normal(size=50),
... np.random.normal(size=30), np.random.normal(size=20)])
(-0.73091722665244196,
array([ 0.44925884,  1.3052767 ,  1.9434184 ,  2.57696569,  3.41634856]),
0.8789283903979661)

ansari(x, y)

Perform the Ansari-Bradley test for equal scale parameters

The Ansari-Bradley test is a non-parametric test for the equality of the scale parameter of the distributions from which two samples were drawn.

x, y : array_like
arrays of sample data
statistic : float
The Ansari-Bradley test statistic
pvalue : float
The p-value of the hypothesis test

fligner : A non-parametric test for the equality of k variances mood : A non-parametric test for the equality of two scale parameters

The p-value given is exact when the sample sizes are both less than 55 and there are no ties, otherwise a normal approximation for the p-value is used.

  Sprent, Peter and N.C. Smeeton. Applied nonparametric statistical methods. 3rd ed. Chapman and Hall/CRC. 2001. Section 5.8.2.
bartlett(*args)

Perform Bartlett’s test for equal variances

Bartlett’s test tests the null hypothesis that all input samples are from populations with equal variances. For samples from significantly non-normal populations, Levene’s test levene is more robust.

sample1, sample2,… : array_like
arrays of sample data. May be different lengths.
statistic : float
The test statistic.
pvalue : float
The p-value of the test.

fligner : A non-parametric test for the equality of k variances levene : A robust parametric test for equality of k variances

Conover et al. (1981) examine many of the existing parametric and nonparametric tests by extensive simulations and they conclude that the tests proposed by Fligner and Killeen (1976) and Levene (1960) appear to be superior in terms of robustness of departures from normality and power _.

  Snedecor, George W. and Cochran, William G. (1989), Statistical Methods, Eighth Edition, Iowa State University Press.
  Park, C. and Lindsay, B. G. (1999). Robust Scale Estimation and Hypothesis Testing based on Quadratic Inference Function. Technical Report #99-03, Center for Likelihood Studies, Pennsylvania State University.
  Bartlett, M. S. (1937). Properties of Sufficiency and Statistical Tests. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, Vol. 160, No.901, pp. 268-282.
levene(*args, **kwds)

Perform Levene test for equal variances.

The Levene test tests the null hypothesis that all input samples are from populations with equal variances. Levene’s test is an alternative to Bartlett’s test bartlett in the case where there are significant deviations from normality.

sample1, sample2, … : array_like
The sample data, possibly with different lengths
center : {‘mean’, ‘median’, ‘trimmed’}, optional
Which function of the data to use in the test. The default is ‘median’.
proportiontocut : float, optional
When center is ‘trimmed’, this gives the proportion of data points to cut from each end. (See scipy.stats.trim_mean.) Default is 0.05.
statistic : float
The test statistic.
pvalue : float
The p-value for the test.

Three variations of Levene’s test are possible. The possibilities and their recommended usages are:

• ‘median’ : Recommended for skewed (non-normal) distributions>
• ‘mean’ : Recommended for symmetric, moderate-tailed distributions.
• ‘trimmed’ : Recommended for heavy-tailed distributions.
  Levene, H. (1960). In Contributions to Probability and Statistics: Essays in Honor of Harold Hotelling, I. Olkin et al. eds., Stanford University Press, pp. 278-292.
  Brown, M. B. and Forsythe, A. B. (1974), Journal of the American Statistical Association, 69, 364-367
binom_test(x, n=None, p=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.

x : integer or array_like
the number of successes, or if x has length 2, it is the number of successes and the number of failures.
n : integer
the number of trials. This is ignored if x gives both the number of successes and failures
p : float, optional
The hypothesized probability of success. 0 <= p <= 1. The default value is p = 0.5
alternative : {‘two-sided’, ‘greater’, ‘less’}, optional
Indicates the alternative hypothesis. The default value is ‘two-sided’.
p-value : float
The p-value of the hypothesis test
_apply_func(x, g, func)
fligner(*args, **kwds)

Perform Fligner-Killeen test for equality of variance.

Fligner’s test tests the null hypothesis that all input samples are from populations with equal variances. Fligner-Killeen’s test is distribution free when populations are identical _.

sample1, sample2, … : array_like
Arrays of sample data. Need not be the same length.
center : {‘mean’, ‘median’, ‘trimmed’}, optional
Keyword argument controlling which function of the data is used in computing the test statistic. The default is ‘median’.
proportiontocut : float, optional
When center is ‘trimmed’, this gives the proportion of data points to cut from each end. (See scipy.stats.trim_mean.) Default is 0.05.
statistic : float
The test statistic.
pvalue : float
The p-value for the hypothesis test.

bartlett : A parametric test for equality of k variances in normal samples levene : A robust parametric test for equality of k variances

As with Levene’s test there are three variants of Fligner’s test that differ by the measure of central tendency used in the test. See levene for more information.

Conover et al. (1981) examine many of the existing parametric and nonparametric tests by extensive simulations and they conclude that the tests proposed by Fligner and Killeen (1976) and Levene (1960) appear to be superior in terms of robustness of departures from normality and power _.

  Park, C. and Lindsay, B. G. (1999). Robust Scale Estimation and Hypothesis Testing based on Quadratic Inference Function. Technical Report #99-03, Center for Likelihood Studies, Pennsylvania State University. http://cecas.clemson.edu/~cspark/cv/paper/qif/draftqif2.pdf
  Fligner, M.A. and Killeen, T.J. (1976). Distribution-free two-sample tests for scale. ‘Journal of the American Statistical Association.’ 71(353), 210-213.
  Park, C. and Lindsay, B. G. (1999). Robust Scale Estimation and Hypothesis Testing based on Quadratic Inference Function. Technical Report #99-03, Center for Likelihood Studies, Pennsylvania State University.
  Conover, W. J., Johnson, M. E. and Johnson M. M. (1981). A comparative study of tests for homogeneity of variances, with applications to the outer continental shelf biding data. Technometrics, 23(4), 351-361.
mood(x, y, axis=0)

Perform Mood’s test for equal scale parameters.

Mood’s two-sample test for scale parameters is a non-parametric test for the null hypothesis that two samples are drawn from the same distribution with the same scale parameter.

x, y : array_like
Arrays of sample data.
axis : int, optional
The axis along which the samples are tested. x and y can be of different length along axis. If axis is None, x and y are flattened and the test is done on all values in the flattened arrays.
z : scalar or ndarray
The z-score for the hypothesis test. For 1-D inputs a scalar is returned.
p-value : scalar ndarray
The p-value for the hypothesis test.

fligner : A non-parametric test for the equality of k variances ansari : A non-parametric test for the equality of 2 variances bartlett : A parametric test for equality of k variances in normal samples levene : A parametric test for equality of k variances

The data are assumed to be drawn from probability distributions f(x) and f(x/s) / s respectively, for some probability density function f. The null hypothesis is that s == 1.

For multi-dimensional arrays, if the inputs are of shapes (n0, n1, n2, n3) and (n0, m1, n2, n3), then if axis=1, the resulting z and p values will have shape (n0, n2, n3). Note that n1 and m1 don’t have to be equal, but the other dimensions do.

>>> from scipy import stats
>>> np.random.seed(1234)
>>> x2 = np.random.randn(2, 45, 6, 7)
>>> x1 = np.random.randn(2, 30, 6, 7)
>>> z, p = stats.mood(x1, x2, axis=1)
>>> p.shape
(2, 6, 7)


Find the number of points where the difference in scale is not significant:

>>> (p > 0.1).sum()
74


Perform the test with different scales:

>>> x1 = np.random.randn(2, 30)
>>> x2 = np.random.randn(2, 35) * 10.0
>>> stats.mood(x1, x2, axis=1)
(array([-5.7178125 , -5.25342163]), array([  1.07904114e-08,   1.49299218e-07]))

wilcoxon(x, y=None, zero_method="wilcox", correction=False)

Calculate the Wilcoxon signed-rank test.

The Wilcoxon signed-rank test tests the null hypothesis that two related paired samples come from the same distribution. In particular, it tests whether the distribution of the differences x - y is symmetric about zero. It is a non-parametric version of the paired T-test.

x : array_like
The first set of measurements.
y : array_like, optional
The second set of measurements. If y is not given, then the x array is considered to be the differences between the two sets of measurements.
zero_method : string, {“pratt”, “wilcox”, “zsplit”}, optional
“pratt”:
Pratt treatment: includes zero-differences in the ranking process (more conservative)
“wilcox”:
Wilcox treatment: discards all zero-differences
“zsplit”:
Zero rank split: just like Pratt, but spliting the zero rank between positive and negative ones
correction : bool, optional
If True, apply continuity correction by adjusting the Wilcoxon rank statistic by 0.5 towards the mean value when computing the z-statistic. Default is False.
statistic : float
The sum of the ranks of the differences above or below zero, whichever is smaller.
pvalue : float
The two-sided p-value for the test.

Because the normal approximation is used for the calculations, the samples used should be large. A typical rule is to require that n > 20.

median_test(*args, **kwds)

Mood’s median test.

Test that two or more samples come from populations with the same median.

Let n = len(args) be the number of samples. The “grand median” of all the data is computed, and a contingency table is formed by classifying the values in each sample as being above or below the grand median. The contingency table, along with correction and lambda_, are passed to scipy.stats.chi2_contingency to compute the test statistic and p-value.

sample1, sample2, … : array_like
The set of samples. There must be at least two samples. Each sample must be a one-dimensional sequence containing at least one value. The samples are not required to have the same length.
ties : str, optional

Determines how values equal to the grand median are classified in the contingency table. The string must be one of:

"below":
Values equal to the grand median are counted as "below".
"above":
Values equal to the grand median are counted as "above".
"ignore":
Values equal to the grand median are not counted.


The default is “below”.

correction : bool, optional
If True, and there are just two samples, apply Yates’ correction for continuity when computing the test statistic associated with the contingency table. Default is True.
lambda_ : float or str, optional.
By default, the statistic computed in this test is Pearson’s chi-squared statistic. lambda_ allows a statistic from the Cressie-Read power divergence family to be used instead. See power_divergence for details. Default is 1 (Pearson’s chi-squared statistic).
nan_policy : {‘propagate’, ‘raise’, ‘omit’}, optional
Defines how to handle when input contains nan. ‘propagate’ returns nan, ‘raise’ throws an error, ‘omit’ performs the calculations ignoring nan values. Default is ‘propagate’.
stat : float
The test statistic. The statistic that is returned is determined by lambda_. The default is Pearson’s chi-squared statistic.
p : float
The p-value of the test.
m : float
The grand median.
table : ndarray
The contingency table. The shape of the table is (2, n), where n is the number of samples. The first row holds the counts of the values above the grand median, and the second row holds the counts of the values below the grand median. The table allows further analysis with, for example, scipy.stats.chi2_contingency, or with scipy.stats.fisher_exact if there are two samples, without having to recompute the table. If nan_policy is “propagate” and there are nans in the input, the return value for table is None.

kruskal : Compute the Kruskal-Wallis H-test for independent samples. mannwhitneyu : Computes the Mann-Whitney rank test on samples x and y.

New in version 0.15.0.

  Mood, A. M., Introduction to the Theory of Statistics. McGraw-Hill (1950), pp. 394-399.
  Zar, J. H., Biostatistical Analysis, 5th ed. Prentice Hall (2010). See Sections 8.12 and 10.15.

A biologist runs an experiment in which there are three groups of plants. Group 1 has 16 plants, group 2 has 15 plants, and group 3 has 17 plants. Each plant produces a number of seeds. The seed counts for each group are:

Group 1: 10 14 14 18 20 22 24 25 31 31 32 39 43 43 48 49
Group 2: 28 30 31 33 34 35 36 40 44 55 57 61 91 92 99
Group 3:  0  3  9 22 23 25 25 33 34 34 40 45 46 48 62 67 84


The following code applies Mood’s median test to these samples.

>>> g1 = [10, 14, 14, 18, 20, 22, 24, 25, 31, 31, 32, 39, 43, 43, 48, 49]
>>> g2 = [28, 30, 31, 33, 34, 35, 36, 40, 44, 55, 57, 61, 91, 92, 99]
>>> g3 = [0, 3, 9, 22, 23, 25, 25, 33, 34, 34, 40, 45, 46, 48, 62, 67, 84]
>>> from scipy.stats import median_test
>>> stat, p, med, tbl = median_test(g1, g2, g3)


The median is

>>> med
34.0


and the contingency table is

>>> tbl
array([[ 5, 10,  7],
[11,  5, 10]])


p is too large to conclude that the medians are not the same:

>>> p
0.12609082774093244


The “G-test” can be performed by passing lambda_="log-likelihood" to median_test.

>>> g, p, med, tbl = median_test(g1, g2, g3, lambda_="log-likelihood")
>>> p
0.12224779737117837


The median occurs several times in the data, so we’ll get a different result if, for example, ties="above" is used:

>>> stat, p, med, tbl = median_test(g1, g2, g3, ties="above")
>>> p
0.063873276069553273

>>> tbl
array([[ 5, 11,  9],
[11,  4,  8]])


This example demonstrates that if the data set is not large and there are values equal to the median, the p-value can be sensitive to the choice of ties.

_circfuncs_common(samples, high, low)
circmean(samples, high=None, low=0, axis=None)

Compute the circular mean for samples in a range.

samples : array_like
Input array.
high : float or int, optional
High boundary for circular mean range. Default is 2*pi.
low : float or int, optional
Low boundary for circular mean range. Default is 0.
axis : int, optional
Axis along which means are computed. The default is to compute the mean of the flattened array.
circmean : float
Circular mean.
>>> from scipy.stats import circmean
>>> circmean([0.1, 2*np.pi+0.2, 6*np.pi+0.3])
0.2

>>> from scipy.stats import circmean
>>> circmean([0.2, 1.4, 2.6], high = 1, low = 0)
0.4

circvar(samples, high=None, low=0, axis=None)

Compute the circular variance for samples assumed to be in a range

samples : array_like
Input array.
low : float or int, optional
Low boundary for circular variance range. Default is 0.
high : float or int, optional
High boundary for circular variance range. Default is 2*pi.
axis : int, optional
Axis along which variances are computed. The default is to compute the variance of the flattened array.
circvar : float
Circular variance.

This uses a definition of circular variance that in the limit of small angles returns a number close to the ‘linear’ variance.

>>> from scipy.stats import circvar
>>> circvar([0, 2*np.pi/3, 5*np.pi/3])
2.19722457734

circstd(samples, high=None, low=0, axis=None)

Compute the circular standard deviation for samples assumed to be in the range [low to high].

samples : array_like
Input array.
low : float or int, optional
Low boundary for circular standard deviation range. Default is 0.
high : float or int, optional
High boundary for circular standard deviation range. Default is 2*pi.
axis : int, optional
Axis along which standard deviations are computed. The default is to compute the standard deviation of the flattened array.
circstd : float
Circular standard deviation.

This uses a definition of circular standard deviation that in the limit of small angles returns a number close to the ‘linear’ standard deviation.

>>> from scipy.stats import circstd
>>> circstd([0, 0.1*np.pi/2, 0.001*np.pi, 0.03*np.pi/2])
0.063564063306