Coverage for /usr/lib/python3/dist-packages/scipy/stats/_morestats.py: 10%
1068 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from __future__ import annotations
2import math
3import warnings
4from collections import namedtuple
6import numpy as np
7from numpy import (isscalar, r_, log, around, unique, asarray, zeros,
8 arange, sort, amin, amax, atleast_1d, sqrt, array,
9 compress, pi, exp, ravel, count_nonzero, sin, cos,
10 arctan2, hypot)
12from scipy import optimize, special, interpolate, stats
13from scipy._lib._bunch import _make_tuple_bunch
14from scipy._lib._util import _rename_parameter, _contains_nan, _get_nan
16from . import _statlib
17from . import _stats_py
18from ._fit import FitResult
19from ._stats_py import find_repeats, _normtest_finish, SignificanceResult
20from .contingency import chi2_contingency
21from . import distributions
22from ._distn_infrastructure import rv_generic
23from ._hypotests import _get_wilcoxon_distr
24from ._axis_nan_policy import _axis_nan_policy_factory
25from .._lib.deprecation import _deprecated
28__all__ = ['mvsdist',
29 'bayes_mvs', 'kstat', 'kstatvar', 'probplot', 'ppcc_max', 'ppcc_plot',
30 'boxcox_llf', 'boxcox', 'boxcox_normmax', 'boxcox_normplot',
31 'shapiro', 'anderson', 'ansari', 'bartlett', 'levene', 'binom_test',
32 'fligner', 'mood', 'wilcoxon', 'median_test',
33 'circmean', 'circvar', 'circstd', 'anderson_ksamp',
34 'yeojohnson_llf', 'yeojohnson', 'yeojohnson_normmax',
35 'yeojohnson_normplot', 'directional_stats',
36 'false_discovery_control'
37 ]
40Mean = namedtuple('Mean', ('statistic', 'minmax'))
41Variance = namedtuple('Variance', ('statistic', 'minmax'))
42Std_dev = namedtuple('Std_dev', ('statistic', 'minmax'))
45def bayes_mvs(data, alpha=0.90):
46 r"""
47 Bayesian confidence intervals for the mean, var, and std.
49 Parameters
50 ----------
51 data : array_like
52 Input data, if multi-dimensional it is flattened to 1-D by `bayes_mvs`.
53 Requires 2 or more data points.
54 alpha : float, optional
55 Probability that the returned confidence interval contains
56 the true parameter.
58 Returns
59 -------
60 mean_cntr, var_cntr, std_cntr : tuple
61 The three results are for the mean, variance and standard deviation,
62 respectively. Each result is a tuple of the form::
64 (center, (lower, upper))
66 with `center` the mean of the conditional pdf of the value given the
67 data, and `(lower, upper)` a confidence interval, centered on the
68 median, containing the estimate to a probability ``alpha``.
70 See Also
71 --------
72 mvsdist
74 Notes
75 -----
76 Each tuple of mean, variance, and standard deviation estimates represent
77 the (center, (lower, upper)) with center the mean of the conditional pdf
78 of the value given the data and (lower, upper) is a confidence interval
79 centered on the median, containing the estimate to a probability
80 ``alpha``.
82 Converts data to 1-D and assumes all data has the same mean and variance.
83 Uses Jeffrey's prior for variance and std.
85 Equivalent to ``tuple((x.mean(), x.interval(alpha)) for x in mvsdist(dat))``
87 References
88 ----------
89 T.E. Oliphant, "A Bayesian perspective on estimating mean, variance, and
90 standard-deviation from data", https://scholarsarchive.byu.edu/facpub/278,
91 2006.
93 Examples
94 --------
95 First a basic example to demonstrate the outputs:
97 >>> from scipy import stats
98 >>> data = [6, 9, 12, 7, 8, 8, 13]
99 >>> mean, var, std = stats.bayes_mvs(data)
100 >>> mean
101 Mean(statistic=9.0, minmax=(7.103650222612533, 10.896349777387467))
102 >>> var
103 Variance(statistic=10.0, minmax=(3.176724206..., 24.45910382...))
104 >>> std
105 Std_dev(statistic=2.9724954732045084, minmax=(1.7823367265645143, 4.945614605014631))
107 Now we generate some normally distributed random data, and get estimates of
108 mean and standard deviation with 95% confidence intervals for those
109 estimates:
111 >>> n_samples = 100000
112 >>> data = stats.norm.rvs(size=n_samples)
113 >>> res_mean, res_var, res_std = stats.bayes_mvs(data, alpha=0.95)
115 >>> import matplotlib.pyplot as plt
116 >>> fig = plt.figure()
117 >>> ax = fig.add_subplot(111)
118 >>> ax.hist(data, bins=100, density=True, label='Histogram of data')
119 >>> ax.vlines(res_mean.statistic, 0, 0.5, colors='r', label='Estimated mean')
120 >>> ax.axvspan(res_mean.minmax[0],res_mean.minmax[1], facecolor='r',
121 ... alpha=0.2, label=r'Estimated mean (95% limits)')
122 >>> ax.vlines(res_std.statistic, 0, 0.5, colors='g', label='Estimated scale')
123 >>> ax.axvspan(res_std.minmax[0],res_std.minmax[1], facecolor='g', alpha=0.2,
124 ... label=r'Estimated scale (95% limits)')
126 >>> ax.legend(fontsize=10)
127 >>> ax.set_xlim([-4, 4])
128 >>> ax.set_ylim([0, 0.5])
129 >>> plt.show()
131 """
132 m, v, s = mvsdist(data)
133 if alpha >= 1 or alpha <= 0:
134 raise ValueError("0 < alpha < 1 is required, but alpha=%s was given."
135 % alpha)
137 m_res = Mean(m.mean(), m.interval(alpha))
138 v_res = Variance(v.mean(), v.interval(alpha))
139 s_res = Std_dev(s.mean(), s.interval(alpha))
141 return m_res, v_res, s_res
144def mvsdist(data):
145 """
146 'Frozen' distributions for mean, variance, and standard deviation of data.
148 Parameters
149 ----------
150 data : array_like
151 Input array. Converted to 1-D using ravel.
152 Requires 2 or more data-points.
154 Returns
155 -------
156 mdist : "frozen" distribution object
157 Distribution object representing the mean of the data.
158 vdist : "frozen" distribution object
159 Distribution object representing the variance of the data.
160 sdist : "frozen" distribution object
161 Distribution object representing the standard deviation of the data.
163 See Also
164 --------
165 bayes_mvs
167 Notes
168 -----
169 The return values from ``bayes_mvs(data)`` is equivalent to
170 ``tuple((x.mean(), x.interval(0.90)) for x in mvsdist(data))``.
172 In other words, calling ``<dist>.mean()`` and ``<dist>.interval(0.90)``
173 on the three distribution objects returned from this function will give
174 the same results that are returned from `bayes_mvs`.
176 References
177 ----------
178 T.E. Oliphant, "A Bayesian perspective on estimating mean, variance, and
179 standard-deviation from data", https://scholarsarchive.byu.edu/facpub/278,
180 2006.
182 Examples
183 --------
184 >>> from scipy import stats
185 >>> data = [6, 9, 12, 7, 8, 8, 13]
186 >>> mean, var, std = stats.mvsdist(data)
188 We now have frozen distribution objects "mean", "var" and "std" that we can
189 examine:
191 >>> mean.mean()
192 9.0
193 >>> mean.interval(0.95)
194 (6.6120585482655692, 11.387941451734431)
195 >>> mean.std()
196 1.1952286093343936
198 """
199 x = ravel(data)
200 n = len(x)
201 if n < 2:
202 raise ValueError("Need at least 2 data-points.")
203 xbar = x.mean()
204 C = x.var()
205 if n > 1000: # gaussian approximations for large n
206 mdist = distributions.norm(loc=xbar, scale=math.sqrt(C / n))
207 sdist = distributions.norm(loc=math.sqrt(C), scale=math.sqrt(C / (2. * n)))
208 vdist = distributions.norm(loc=C, scale=math.sqrt(2.0 / n) * C)
209 else:
210 nm1 = n - 1
211 fac = n * C / 2.
212 val = nm1 / 2.
213 mdist = distributions.t(nm1, loc=xbar, scale=math.sqrt(C / nm1))
214 sdist = distributions.gengamma(val, -2, scale=math.sqrt(fac))
215 vdist = distributions.invgamma(val, scale=fac)
216 return mdist, vdist, sdist
219@_axis_nan_policy_factory(
220 lambda x: x, result_to_tuple=lambda x: (x,), n_outputs=1, default_axis=None
221)
222def kstat(data, n=2):
223 r"""
224 Return the nth k-statistic (1<=n<=4 so far).
226 The nth k-statistic k_n is the unique symmetric unbiased estimator of the
227 nth cumulant kappa_n.
229 Parameters
230 ----------
231 data : array_like
232 Input array. Note that n-D input gets flattened.
233 n : int, {1, 2, 3, 4}, optional
234 Default is equal to 2.
236 Returns
237 -------
238 kstat : float
239 The nth k-statistic.
241 See Also
242 --------
243 kstatvar : Returns an unbiased estimator of the variance of the k-statistic
244 moment : Returns the n-th central moment about the mean for a sample.
246 Notes
247 -----
248 For a sample size n, the first few k-statistics are given by:
250 .. math::
252 k_{1} = \mu
253 k_{2} = \frac{n}{n-1} m_{2}
254 k_{3} = \frac{ n^{2} } {(n-1) (n-2)} m_{3}
255 k_{4} = \frac{ n^{2} [(n + 1)m_{4} - 3(n - 1) m^2_{2}]} {(n-1) (n-2) (n-3)}
257 where :math:`\mu` is the sample mean, :math:`m_2` is the sample
258 variance, and :math:`m_i` is the i-th sample central moment.
260 References
261 ----------
262 http://mathworld.wolfram.com/k-Statistic.html
264 http://mathworld.wolfram.com/Cumulant.html
266 Examples
267 --------
268 >>> from scipy import stats
269 >>> from numpy.random import default_rng
270 >>> rng = default_rng()
272 As sample size increases, n-th moment and n-th k-statistic converge to the
273 same number (although they aren't identical). In the case of the normal
274 distribution, they converge to zero.
276 >>> for n in [2, 3, 4, 5, 6, 7]:
277 ... x = rng.normal(size=10**n)
278 ... m, k = stats.moment(x, 3), stats.kstat(x, 3)
279 ... print("%.3g %.3g %.3g" % (m, k, m-k))
280 -0.631 -0.651 0.0194 # random
281 0.0282 0.0283 -8.49e-05
282 -0.0454 -0.0454 1.36e-05
283 7.53e-05 7.53e-05 -2.26e-09
284 0.00166 0.00166 -4.99e-09
285 -2.88e-06 -2.88e-06 8.63e-13
286 """
287 if n > 4 or n < 1:
288 raise ValueError("k-statistics only supported for 1<=n<=4")
289 n = int(n)
290 S = np.zeros(n + 1, np.float64)
291 data = ravel(data)
292 N = data.size
294 # raise ValueError on empty input
295 if N == 0:
296 raise ValueError("Data input must not be empty")
298 # on nan input, return nan without warning
299 if np.isnan(np.sum(data)):
300 return np.nan
302 for k in range(1, n + 1):
303 S[k] = np.sum(data**k, axis=0)
304 if n == 1:
305 return S[1] * 1.0/N
306 elif n == 2:
307 return (N*S[2] - S[1]**2.0) / (N*(N - 1.0))
308 elif n == 3:
309 return (2*S[1]**3 - 3*N*S[1]*S[2] + N*N*S[3]) / (N*(N - 1.0)*(N - 2.0))
310 elif n == 4:
311 return ((-6*S[1]**4 + 12*N*S[1]**2 * S[2] - 3*N*(N-1.0)*S[2]**2 -
312 4*N*(N+1)*S[1]*S[3] + N*N*(N+1)*S[4]) /
313 (N*(N-1.0)*(N-2.0)*(N-3.0)))
314 else:
315 raise ValueError("Should not be here.")
318@_axis_nan_policy_factory(
319 lambda x: x, result_to_tuple=lambda x: (x,), n_outputs=1, default_axis=None
320)
321def kstatvar(data, n=2):
322 r"""Return an unbiased estimator of the variance of the k-statistic.
324 See `kstat` for more details of the k-statistic.
326 Parameters
327 ----------
328 data : array_like
329 Input array. Note that n-D input gets flattened.
330 n : int, {1, 2}, optional
331 Default is equal to 2.
333 Returns
334 -------
335 kstatvar : float
336 The nth k-statistic variance.
338 See Also
339 --------
340 kstat : Returns the n-th k-statistic.
341 moment : Returns the n-th central moment about the mean for a sample.
343 Notes
344 -----
345 The variances of the first few k-statistics are given by:
347 .. math::
349 var(k_{1}) = \frac{\kappa^2}{n}
350 var(k_{2}) = \frac{\kappa^4}{n} + \frac{2\kappa^2_{2}}{n - 1}
351 var(k_{3}) = \frac{\kappa^6}{n} + \frac{9 \kappa_2 \kappa_4}{n - 1} +
352 \frac{9 \kappa^2_{3}}{n - 1} +
353 \frac{6 n \kappa^3_{2}}{(n-1) (n-2)}
354 var(k_{4}) = \frac{\kappa^8}{n} + \frac{16 \kappa_2 \kappa_6}{n - 1} +
355 \frac{48 \kappa_{3} \kappa_5}{n - 1} +
356 \frac{34 \kappa^2_{4}}{n-1} + \frac{72 n \kappa^2_{2} \kappa_4}{(n - 1) (n - 2)} +
357 \frac{144 n \kappa_{2} \kappa^2_{3}}{(n - 1) (n - 2)} +
358 \frac{24 (n + 1) n \kappa^4_{2}}{(n - 1) (n - 2) (n - 3)}
359 """
360 data = ravel(data)
361 N = len(data)
362 if n == 1:
363 return kstat(data, n=2) * 1.0/N
364 elif n == 2:
365 k2 = kstat(data, n=2)
366 k4 = kstat(data, n=4)
367 return (2*N*k2**2 + (N-1)*k4) / (N*(N+1))
368 else:
369 raise ValueError("Only n=1 or n=2 supported.")
372def _calc_uniform_order_statistic_medians(n):
373 """Approximations of uniform order statistic medians.
375 Parameters
376 ----------
377 n : int
378 Sample size.
380 Returns
381 -------
382 v : 1d float array
383 Approximations of the order statistic medians.
385 References
386 ----------
387 .. [1] James J. Filliben, "The Probability Plot Correlation Coefficient
388 Test for Normality", Technometrics, Vol. 17, pp. 111-117, 1975.
390 Examples
391 --------
392 Order statistics of the uniform distribution on the unit interval
393 are marginally distributed according to beta distributions.
394 The expectations of these order statistic are evenly spaced across
395 the interval, but the distributions are skewed in a way that
396 pushes the medians slightly towards the endpoints of the unit interval:
398 >>> import numpy as np
399 >>> n = 4
400 >>> k = np.arange(1, n+1)
401 >>> from scipy.stats import beta
402 >>> a = k
403 >>> b = n-k+1
404 >>> beta.mean(a, b)
405 array([0.2, 0.4, 0.6, 0.8])
406 >>> beta.median(a, b)
407 array([0.15910358, 0.38572757, 0.61427243, 0.84089642])
409 The Filliben approximation uses the exact medians of the smallest
410 and greatest order statistics, and the remaining medians are approximated
411 by points spread evenly across a sub-interval of the unit interval:
413 >>> from scipy.stats._morestats import _calc_uniform_order_statistic_medians
414 >>> _calc_uniform_order_statistic_medians(n)
415 array([0.15910358, 0.38545246, 0.61454754, 0.84089642])
417 This plot shows the skewed distributions of the order statistics
418 of a sample of size four from a uniform distribution on the unit interval:
420 >>> import matplotlib.pyplot as plt
421 >>> x = np.linspace(0.0, 1.0, num=50, endpoint=True)
422 >>> pdfs = [beta.pdf(x, a[i], b[i]) for i in range(n)]
423 >>> plt.figure()
424 >>> plt.plot(x, pdfs[0], x, pdfs[1], x, pdfs[2], x, pdfs[3])
426 """
427 v = np.empty(n, dtype=np.float64)
428 v[-1] = 0.5**(1.0 / n)
429 v[0] = 1 - v[-1]
430 i = np.arange(2, n)
431 v[1:-1] = (i - 0.3175) / (n + 0.365)
432 return v
435def _parse_dist_kw(dist, enforce_subclass=True):
436 """Parse `dist` keyword.
438 Parameters
439 ----------
440 dist : str or stats.distributions instance.
441 Several functions take `dist` as a keyword, hence this utility
442 function.
443 enforce_subclass : bool, optional
444 If True (default), `dist` needs to be a
445 `_distn_infrastructure.rv_generic` instance.
446 It can sometimes be useful to set this keyword to False, if a function
447 wants to accept objects that just look somewhat like such an instance
448 (for example, they have a ``ppf`` method).
450 """
451 if isinstance(dist, rv_generic):
452 pass
453 elif isinstance(dist, str):
454 try:
455 dist = getattr(distributions, dist)
456 except AttributeError as e:
457 raise ValueError("%s is not a valid distribution name" % dist) from e
458 elif enforce_subclass:
459 msg = ("`dist` should be a stats.distributions instance or a string "
460 "with the name of such a distribution.")
461 raise ValueError(msg)
463 return dist
466def _add_axis_labels_title(plot, xlabel, ylabel, title):
467 """Helper function to add axes labels and a title to stats plots."""
468 try:
469 if hasattr(plot, 'set_title'):
470 # Matplotlib Axes instance or something that looks like it
471 plot.set_title(title)
472 plot.set_xlabel(xlabel)
473 plot.set_ylabel(ylabel)
474 else:
475 # matplotlib.pyplot module
476 plot.title(title)
477 plot.xlabel(xlabel)
478 plot.ylabel(ylabel)
479 except Exception:
480 # Not an MPL object or something that looks (enough) like it.
481 # Don't crash on adding labels or title
482 pass
485def probplot(x, sparams=(), dist='norm', fit=True, plot=None, rvalue=False):
486 """
487 Calculate quantiles for a probability plot, and optionally show the plot.
489 Generates a probability plot of sample data against the quantiles of a
490 specified theoretical distribution (the normal distribution by default).
491 `probplot` optionally calculates a best-fit line for the data and plots the
492 results using Matplotlib or a given plot function.
494 Parameters
495 ----------
496 x : array_like
497 Sample/response data from which `probplot` creates the plot.
498 sparams : tuple, optional
499 Distribution-specific shape parameters (shape parameters plus location
500 and scale).
501 dist : str or stats.distributions instance, optional
502 Distribution or distribution function name. The default is 'norm' for a
503 normal probability plot. Objects that look enough like a
504 stats.distributions instance (i.e. they have a ``ppf`` method) are also
505 accepted.
506 fit : bool, optional
507 Fit a least-squares regression (best-fit) line to the sample data if
508 True (default).
509 plot : object, optional
510 If given, plots the quantiles.
511 If given and `fit` is True, also plots the least squares fit.
512 `plot` is an object that has to have methods "plot" and "text".
513 The `matplotlib.pyplot` module or a Matplotlib Axes object can be used,
514 or a custom object with the same methods.
515 Default is None, which means that no plot is created.
516 rvalue : bool, optional
517 If `plot` is provided and `fit` is True, setting `rvalue` to True
518 includes the coefficient of determination on the plot.
519 Default is False.
521 Returns
522 -------
523 (osm, osr) : tuple of ndarrays
524 Tuple of theoretical quantiles (osm, or order statistic medians) and
525 ordered responses (osr). `osr` is simply sorted input `x`.
526 For details on how `osm` is calculated see the Notes section.
527 (slope, intercept, r) : tuple of floats, optional
528 Tuple containing the result of the least-squares fit, if that is
529 performed by `probplot`. `r` is the square root of the coefficient of
530 determination. If ``fit=False`` and ``plot=None``, this tuple is not
531 returned.
533 Notes
534 -----
535 Even if `plot` is given, the figure is not shown or saved by `probplot`;
536 ``plt.show()`` or ``plt.savefig('figname.png')`` should be used after
537 calling `probplot`.
539 `probplot` generates a probability plot, which should not be confused with
540 a Q-Q or a P-P plot. Statsmodels has more extensive functionality of this
541 type, see ``statsmodels.api.ProbPlot``.
543 The formula used for the theoretical quantiles (horizontal axis of the
544 probability plot) is Filliben's estimate::
546 quantiles = dist.ppf(val), for
548 0.5**(1/n), for i = n
549 val = (i - 0.3175) / (n + 0.365), for i = 2, ..., n-1
550 1 - 0.5**(1/n), for i = 1
552 where ``i`` indicates the i-th ordered value and ``n`` is the total number
553 of values.
555 Examples
556 --------
557 >>> import numpy as np
558 >>> from scipy import stats
559 >>> import matplotlib.pyplot as plt
560 >>> nsample = 100
561 >>> rng = np.random.default_rng()
563 A t distribution with small degrees of freedom:
565 >>> ax1 = plt.subplot(221)
566 >>> x = stats.t.rvs(3, size=nsample, random_state=rng)
567 >>> res = stats.probplot(x, plot=plt)
569 A t distribution with larger degrees of freedom:
571 >>> ax2 = plt.subplot(222)
572 >>> x = stats.t.rvs(25, size=nsample, random_state=rng)
573 >>> res = stats.probplot(x, plot=plt)
575 A mixture of two normal distributions with broadcasting:
577 >>> ax3 = plt.subplot(223)
578 >>> x = stats.norm.rvs(loc=[0,5], scale=[1,1.5],
579 ... size=(nsample//2,2), random_state=rng).ravel()
580 >>> res = stats.probplot(x, plot=plt)
582 A standard normal distribution:
584 >>> ax4 = plt.subplot(224)
585 >>> x = stats.norm.rvs(loc=0, scale=1, size=nsample, random_state=rng)
586 >>> res = stats.probplot(x, plot=plt)
588 Produce a new figure with a loggamma distribution, using the ``dist`` and
589 ``sparams`` keywords:
591 >>> fig = plt.figure()
592 >>> ax = fig.add_subplot(111)
593 >>> x = stats.loggamma.rvs(c=2.5, size=500, random_state=rng)
594 >>> res = stats.probplot(x, dist=stats.loggamma, sparams=(2.5,), plot=ax)
595 >>> ax.set_title("Probplot for loggamma dist with shape parameter 2.5")
597 Show the results with Matplotlib:
599 >>> plt.show()
601 """
602 x = np.asarray(x)
603 if x.size == 0:
604 if fit:
605 return (x, x), (np.nan, np.nan, 0.0)
606 else:
607 return x, x
609 osm_uniform = _calc_uniform_order_statistic_medians(len(x))
610 dist = _parse_dist_kw(dist, enforce_subclass=False)
611 if sparams is None:
612 sparams = ()
613 if isscalar(sparams):
614 sparams = (sparams,)
615 if not isinstance(sparams, tuple):
616 sparams = tuple(sparams)
618 osm = dist.ppf(osm_uniform, *sparams)
619 osr = sort(x)
620 if fit:
621 # perform a linear least squares fit.
622 slope, intercept, r, prob, _ = _stats_py.linregress(osm, osr)
624 if plot is not None:
625 plot.plot(osm, osr, 'bo')
626 if fit:
627 plot.plot(osm, slope*osm + intercept, 'r-')
628 _add_axis_labels_title(plot, xlabel='Theoretical quantiles',
629 ylabel='Ordered Values',
630 title='Probability Plot')
632 # Add R^2 value to the plot as text
633 if fit and rvalue:
634 xmin = amin(osm)
635 xmax = amax(osm)
636 ymin = amin(x)
637 ymax = amax(x)
638 posx = xmin + 0.70 * (xmax - xmin)
639 posy = ymin + 0.01 * (ymax - ymin)
640 plot.text(posx, posy, "$R^2=%1.4f$" % r**2)
642 if fit:
643 return (osm, osr), (slope, intercept, r)
644 else:
645 return osm, osr
648def ppcc_max(x, brack=(0.0, 1.0), dist='tukeylambda'):
649 """Calculate the shape parameter that maximizes the PPCC.
651 The probability plot correlation coefficient (PPCC) plot can be used
652 to determine the optimal shape parameter for a one-parameter family
653 of distributions. ``ppcc_max`` returns the shape parameter that would
654 maximize the probability plot correlation coefficient for the given
655 data to a one-parameter family of distributions.
657 Parameters
658 ----------
659 x : array_like
660 Input array.
661 brack : tuple, optional
662 Triple (a,b,c) where (a<b<c). If bracket consists of two numbers (a, c)
663 then they are assumed to be a starting interval for a downhill bracket
664 search (see `scipy.optimize.brent`).
665 dist : str or stats.distributions instance, optional
666 Distribution or distribution function name. Objects that look enough
667 like a stats.distributions instance (i.e. they have a ``ppf`` method)
668 are also accepted. The default is ``'tukeylambda'``.
670 Returns
671 -------
672 shape_value : float
673 The shape parameter at which the probability plot correlation
674 coefficient reaches its max value.
676 See Also
677 --------
678 ppcc_plot, probplot, boxcox
680 Notes
681 -----
682 The brack keyword serves as a starting point which is useful in corner
683 cases. One can use a plot to obtain a rough visual estimate of the location
684 for the maximum to start the search near it.
686 References
687 ----------
688 .. [1] J.J. Filliben, "The Probability Plot Correlation Coefficient Test
689 for Normality", Technometrics, Vol. 17, pp. 111-117, 1975.
690 .. [2] Engineering Statistics Handbook, NIST/SEMATEC,
691 https://www.itl.nist.gov/div898/handbook/eda/section3/ppccplot.htm
693 Examples
694 --------
695 First we generate some random data from a Weibull distribution
696 with shape parameter 2.5:
698 >>> import numpy as np
699 >>> from scipy import stats
700 >>> import matplotlib.pyplot as plt
701 >>> rng = np.random.default_rng()
702 >>> c = 2.5
703 >>> x = stats.weibull_min.rvs(c, scale=4, size=2000, random_state=rng)
705 Generate the PPCC plot for this data with the Weibull distribution.
707 >>> fig, ax = plt.subplots(figsize=(8, 6))
708 >>> res = stats.ppcc_plot(x, c/2, 2*c, dist='weibull_min', plot=ax)
710 We calculate the value where the shape should reach its maximum and a
711 red line is drawn there. The line should coincide with the highest
712 point in the PPCC graph.
714 >>> cmax = stats.ppcc_max(x, brack=(c/2, 2*c), dist='weibull_min')
715 >>> ax.axvline(cmax, color='r')
716 >>> plt.show()
718 """
719 dist = _parse_dist_kw(dist)
720 osm_uniform = _calc_uniform_order_statistic_medians(len(x))
721 osr = sort(x)
723 # this function computes the x-axis values of the probability plot
724 # and computes a linear regression (including the correlation)
725 # and returns 1-r so that a minimization function maximizes the
726 # correlation
727 def tempfunc(shape, mi, yvals, func):
728 xvals = func(mi, shape)
729 r, prob = _stats_py.pearsonr(xvals, yvals)
730 return 1 - r
732 return optimize.brent(tempfunc, brack=brack,
733 args=(osm_uniform, osr, dist.ppf))
736def ppcc_plot(x, a, b, dist='tukeylambda', plot=None, N=80):
737 """Calculate and optionally plot probability plot correlation coefficient.
739 The probability plot correlation coefficient (PPCC) plot can be used to
740 determine the optimal shape parameter for a one-parameter family of
741 distributions. It cannot be used for distributions without shape
742 parameters
743 (like the normal distribution) or with multiple shape parameters.
745 By default a Tukey-Lambda distribution (`stats.tukeylambda`) is used. A
746 Tukey-Lambda PPCC plot interpolates from long-tailed to short-tailed
747 distributions via an approximately normal one, and is therefore
748 particularly useful in practice.
750 Parameters
751 ----------
752 x : array_like
753 Input array.
754 a, b : scalar
755 Lower and upper bounds of the shape parameter to use.
756 dist : str or stats.distributions instance, optional
757 Distribution or distribution function name. Objects that look enough
758 like a stats.distributions instance (i.e. they have a ``ppf`` method)
759 are also accepted. The default is ``'tukeylambda'``.
760 plot : object, optional
761 If given, plots PPCC against the shape parameter.
762 `plot` is an object that has to have methods "plot" and "text".
763 The `matplotlib.pyplot` module or a Matplotlib Axes object can be used,
764 or a custom object with the same methods.
765 Default is None, which means that no plot is created.
766 N : int, optional
767 Number of points on the horizontal axis (equally distributed from
768 `a` to `b`).
770 Returns
771 -------
772 svals : ndarray
773 The shape values for which `ppcc` was calculated.
774 ppcc : ndarray
775 The calculated probability plot correlation coefficient values.
777 See Also
778 --------
779 ppcc_max, probplot, boxcox_normplot, tukeylambda
781 References
782 ----------
783 J.J. Filliben, "The Probability Plot Correlation Coefficient Test for
784 Normality", Technometrics, Vol. 17, pp. 111-117, 1975.
786 Examples
787 --------
788 First we generate some random data from a Weibull distribution
789 with shape parameter 2.5, and plot the histogram of the data:
791 >>> import numpy as np
792 >>> from scipy import stats
793 >>> import matplotlib.pyplot as plt
794 >>> rng = np.random.default_rng()
795 >>> c = 2.5
796 >>> x = stats.weibull_min.rvs(c, scale=4, size=2000, random_state=rng)
798 Take a look at the histogram of the data.
800 >>> fig1, ax = plt.subplots(figsize=(9, 4))
801 >>> ax.hist(x, bins=50)
802 >>> ax.set_title('Histogram of x')
803 >>> plt.show()
805 Now we explore this data with a PPCC plot as well as the related
806 probability plot and Box-Cox normplot. A red line is drawn where we
807 expect the PPCC value to be maximal (at the shape parameter ``c``
808 used above):
810 >>> fig2 = plt.figure(figsize=(12, 4))
811 >>> ax1 = fig2.add_subplot(1, 3, 1)
812 >>> ax2 = fig2.add_subplot(1, 3, 2)
813 >>> ax3 = fig2.add_subplot(1, 3, 3)
814 >>> res = stats.probplot(x, plot=ax1)
815 >>> res = stats.boxcox_normplot(x, -4, 4, plot=ax2)
816 >>> res = stats.ppcc_plot(x, c/2, 2*c, dist='weibull_min', plot=ax3)
817 >>> ax3.axvline(c, color='r')
818 >>> plt.show()
820 """
821 if b <= a:
822 raise ValueError("`b` has to be larger than `a`.")
824 svals = np.linspace(a, b, num=N)
825 ppcc = np.empty_like(svals)
826 for k, sval in enumerate(svals):
827 _, r2 = probplot(x, sval, dist=dist, fit=True)
828 ppcc[k] = r2[-1]
830 if plot is not None:
831 plot.plot(svals, ppcc, 'x')
832 _add_axis_labels_title(plot, xlabel='Shape Values',
833 ylabel='Prob Plot Corr. Coef.',
834 title='(%s) PPCC Plot' % dist)
836 return svals, ppcc
839def boxcox_llf(lmb, data):
840 r"""The boxcox log-likelihood function.
842 Parameters
843 ----------
844 lmb : scalar
845 Parameter for Box-Cox transformation. See `boxcox` for details.
846 data : array_like
847 Data to calculate Box-Cox log-likelihood for. If `data` is
848 multi-dimensional, the log-likelihood is calculated along the first
849 axis.
851 Returns
852 -------
853 llf : float or ndarray
854 Box-Cox log-likelihood of `data` given `lmb`. A float for 1-D `data`,
855 an array otherwise.
857 See Also
858 --------
859 boxcox, probplot, boxcox_normplot, boxcox_normmax
861 Notes
862 -----
863 The Box-Cox log-likelihood function is defined here as
865 .. math::
867 llf = (\lambda - 1) \sum_i(\log(x_i)) -
868 N/2 \log(\sum_i (y_i - \bar{y})^2 / N),
870 where ``y`` is the Box-Cox transformed input data ``x``.
872 Examples
873 --------
874 >>> import numpy as np
875 >>> from scipy import stats
876 >>> import matplotlib.pyplot as plt
877 >>> from mpl_toolkits.axes_grid1.inset_locator import inset_axes
879 Generate some random variates and calculate Box-Cox log-likelihood values
880 for them for a range of ``lmbda`` values:
882 >>> rng = np.random.default_rng()
883 >>> x = stats.loggamma.rvs(5, loc=10, size=1000, random_state=rng)
884 >>> lmbdas = np.linspace(-2, 10)
885 >>> llf = np.zeros(lmbdas.shape, dtype=float)
886 >>> for ii, lmbda in enumerate(lmbdas):
887 ... llf[ii] = stats.boxcox_llf(lmbda, x)
889 Also find the optimal lmbda value with `boxcox`:
891 >>> x_most_normal, lmbda_optimal = stats.boxcox(x)
893 Plot the log-likelihood as function of lmbda. Add the optimal lmbda as a
894 horizontal line to check that that's really the optimum:
896 >>> fig = plt.figure()
897 >>> ax = fig.add_subplot(111)
898 >>> ax.plot(lmbdas, llf, 'b.-')
899 >>> ax.axhline(stats.boxcox_llf(lmbda_optimal, x), color='r')
900 >>> ax.set_xlabel('lmbda parameter')
901 >>> ax.set_ylabel('Box-Cox log-likelihood')
903 Now add some probability plots to show that where the log-likelihood is
904 maximized the data transformed with `boxcox` looks closest to normal:
906 >>> locs = [3, 10, 4] # 'lower left', 'center', 'lower right'
907 >>> for lmbda, loc in zip([-1, lmbda_optimal, 9], locs):
908 ... xt = stats.boxcox(x, lmbda=lmbda)
909 ... (osm, osr), (slope, intercept, r_sq) = stats.probplot(xt)
910 ... ax_inset = inset_axes(ax, width="20%", height="20%", loc=loc)
911 ... ax_inset.plot(osm, osr, 'c.', osm, slope*osm + intercept, 'k-')
912 ... ax_inset.set_xticklabels([])
913 ... ax_inset.set_yticklabels([])
914 ... ax_inset.set_title(r'$\lambda=%1.2f$' % lmbda)
916 >>> plt.show()
918 """
919 data = np.asarray(data)
920 N = data.shape[0]
921 if N == 0:
922 return np.nan
924 logdata = np.log(data)
926 # Compute the variance of the transformed data.
927 if lmb == 0:
928 variance = np.var(logdata, axis=0)
929 else:
930 # Transform without the constant offset 1/lmb. The offset does
931 # not effect the variance, and the subtraction of the offset can
932 # lead to loss of precision.
933 variance = np.var(data**lmb / lmb, axis=0)
935 return (lmb - 1) * np.sum(logdata, axis=0) - N/2 * np.log(variance)
938def _boxcox_conf_interval(x, lmax, alpha):
939 # Need to find the lambda for which
940 # f(x,lmbda) >= f(x,lmax) - 0.5*chi^2_alpha;1
941 fac = 0.5 * distributions.chi2.ppf(1 - alpha, 1)
942 target = boxcox_llf(lmax, x) - fac
944 def rootfunc(lmbda, data, target):
945 return boxcox_llf(lmbda, data) - target
947 # Find positive endpoint of interval in which answer is to be found
948 newlm = lmax + 0.5
949 N = 0
950 while (rootfunc(newlm, x, target) > 0.0) and (N < 500):
951 newlm += 0.1
952 N += 1
954 if N == 500:
955 raise RuntimeError("Could not find endpoint.")
957 lmplus = optimize.brentq(rootfunc, lmax, newlm, args=(x, target))
959 # Now find negative interval in the same way
960 newlm = lmax - 0.5
961 N = 0
962 while (rootfunc(newlm, x, target) > 0.0) and (N < 500):
963 newlm -= 0.1
964 N += 1
966 if N == 500:
967 raise RuntimeError("Could not find endpoint.")
969 lmminus = optimize.brentq(rootfunc, newlm, lmax, args=(x, target))
970 return lmminus, lmplus
973def boxcox(x, lmbda=None, alpha=None, optimizer=None):
974 r"""Return a dataset transformed by a Box-Cox power transformation.
976 Parameters
977 ----------
978 x : ndarray
979 Input array to be transformed.
981 If `lmbda` is not None, this is an alias of
982 `scipy.special.boxcox`.
983 Returns nan if ``x < 0``; returns -inf if ``x == 0 and lmbda < 0``.
985 If `lmbda` is None, array must be positive, 1-dimensional, and
986 non-constant.
988 lmbda : scalar, optional
989 If `lmbda` is None (default), find the value of `lmbda` that maximizes
990 the log-likelihood function and return it as the second output
991 argument.
993 If `lmbda` is not None, do the transformation for that value.
995 alpha : float, optional
996 If `lmbda` is None and `alpha` is not None (default), return the
997 ``100 * (1-alpha)%`` confidence interval for `lmbda` as the third
998 output argument. Must be between 0.0 and 1.0.
1000 If `lmbda` is not None, `alpha` is ignored.
1001 optimizer : callable, optional
1002 If `lmbda` is None, `optimizer` is the scalar optimizer used to find
1003 the value of `lmbda` that minimizes the negative log-likelihood
1004 function. `optimizer` is a callable that accepts one argument:
1006 fun : callable
1007 The objective function, which evaluates the negative
1008 log-likelihood function at a provided value of `lmbda`
1010 and returns an object, such as an instance of
1011 `scipy.optimize.OptimizeResult`, which holds the optimal value of
1012 `lmbda` in an attribute `x`.
1014 See the example in `boxcox_normmax` or the documentation of
1015 `scipy.optimize.minimize_scalar` for more information.
1017 If `lmbda` is not None, `optimizer` is ignored.
1019 Returns
1020 -------
1021 boxcox : ndarray
1022 Box-Cox power transformed array.
1023 maxlog : float, optional
1024 If the `lmbda` parameter is None, the second returned argument is
1025 the `lmbda` that maximizes the log-likelihood function.
1026 (min_ci, max_ci) : tuple of float, optional
1027 If `lmbda` parameter is None and `alpha` is not None, this returned
1028 tuple of floats represents the minimum and maximum confidence limits
1029 given `alpha`.
1031 See Also
1032 --------
1033 probplot, boxcox_normplot, boxcox_normmax, boxcox_llf
1035 Notes
1036 -----
1037 The Box-Cox transform is given by::
1039 y = (x**lmbda - 1) / lmbda, for lmbda != 0
1040 log(x), for lmbda = 0
1042 `boxcox` requires the input data to be positive. Sometimes a Box-Cox
1043 transformation provides a shift parameter to achieve this; `boxcox` does
1044 not. Such a shift parameter is equivalent to adding a positive constant to
1045 `x` before calling `boxcox`.
1047 The confidence limits returned when `alpha` is provided give the interval
1048 where:
1050 .. math::
1052 llf(\hat{\lambda}) - llf(\lambda) < \frac{1}{2}\chi^2(1 - \alpha, 1),
1054 with ``llf`` the log-likelihood function and :math:`\chi^2` the chi-squared
1055 function.
1057 References
1058 ----------
1059 G.E.P. Box and D.R. Cox, "An Analysis of Transformations", Journal of the
1060 Royal Statistical Society B, 26, 211-252 (1964).
1062 Examples
1063 --------
1064 >>> from scipy import stats
1065 >>> import matplotlib.pyplot as plt
1067 We generate some random variates from a non-normal distribution and make a
1068 probability plot for it, to show it is non-normal in the tails:
1070 >>> fig = plt.figure()
1071 >>> ax1 = fig.add_subplot(211)
1072 >>> x = stats.loggamma.rvs(5, size=500) + 5
1073 >>> prob = stats.probplot(x, dist=stats.norm, plot=ax1)
1074 >>> ax1.set_xlabel('')
1075 >>> ax1.set_title('Probplot against normal distribution')
1077 We now use `boxcox` to transform the data so it's closest to normal:
1079 >>> ax2 = fig.add_subplot(212)
1080 >>> xt, _ = stats.boxcox(x)
1081 >>> prob = stats.probplot(xt, dist=stats.norm, plot=ax2)
1082 >>> ax2.set_title('Probplot after Box-Cox transformation')
1084 >>> plt.show()
1086 """
1087 x = np.asarray(x)
1089 if lmbda is not None: # single transformation
1090 return special.boxcox(x, lmbda)
1092 if x.ndim != 1:
1093 raise ValueError("Data must be 1-dimensional.")
1095 if x.size == 0:
1096 return x
1098 if np.all(x == x[0]):
1099 raise ValueError("Data must not be constant.")
1101 if np.any(x <= 0):
1102 raise ValueError("Data must be positive.")
1104 # If lmbda=None, find the lmbda that maximizes the log-likelihood function.
1105 lmax = boxcox_normmax(x, method='mle', optimizer=optimizer)
1106 y = boxcox(x, lmax)
1108 if alpha is None:
1109 return y, lmax
1110 else:
1111 # Find confidence interval
1112 interval = _boxcox_conf_interval(x, lmax, alpha)
1113 return y, lmax, interval
1116def boxcox_normmax(x, brack=None, method='pearsonr', optimizer=None):
1117 """Compute optimal Box-Cox transform parameter for input data.
1119 Parameters
1120 ----------
1121 x : array_like
1122 Input array.
1123 brack : 2-tuple, optional, default (-2.0, 2.0)
1124 The starting interval for a downhill bracket search for the default
1125 `optimize.brent` solver. Note that this is in most cases not
1126 critical; the final result is allowed to be outside this bracket.
1127 If `optimizer` is passed, `brack` must be None.
1128 method : str, optional
1129 The method to determine the optimal transform parameter (`boxcox`
1130 ``lmbda`` parameter). Options are:
1132 'pearsonr' (default)
1133 Maximizes the Pearson correlation coefficient between
1134 ``y = boxcox(x)`` and the expected values for ``y`` if `x` would be
1135 normally-distributed.
1137 'mle'
1138 Minimizes the log-likelihood `boxcox_llf`. This is the method used
1139 in `boxcox`.
1141 'all'
1142 Use all optimization methods available, and return all results.
1143 Useful to compare different methods.
1144 optimizer : callable, optional
1145 `optimizer` is a callable that accepts one argument:
1147 fun : callable
1148 The objective function to be optimized. `fun` accepts one argument,
1149 the Box-Cox transform parameter `lmbda`, and returns the negative
1150 log-likelihood function at the provided value. The job of `optimizer`
1151 is to find the value of `lmbda` that minimizes `fun`.
1153 and returns an object, such as an instance of
1154 `scipy.optimize.OptimizeResult`, which holds the optimal value of
1155 `lmbda` in an attribute `x`.
1157 See the example below or the documentation of
1158 `scipy.optimize.minimize_scalar` for more information.
1160 Returns
1161 -------
1162 maxlog : float or ndarray
1163 The optimal transform parameter found. An array instead of a scalar
1164 for ``method='all'``.
1166 See Also
1167 --------
1168 boxcox, boxcox_llf, boxcox_normplot, scipy.optimize.minimize_scalar
1170 Examples
1171 --------
1172 >>> import numpy as np
1173 >>> from scipy import stats
1174 >>> import matplotlib.pyplot as plt
1176 We can generate some data and determine the optimal ``lmbda`` in various
1177 ways:
1179 >>> rng = np.random.default_rng()
1180 >>> x = stats.loggamma.rvs(5, size=30, random_state=rng) + 5
1181 >>> y, lmax_mle = stats.boxcox(x)
1182 >>> lmax_pearsonr = stats.boxcox_normmax(x)
1184 >>> lmax_mle
1185 2.217563431465757
1186 >>> lmax_pearsonr
1187 2.238318660200961
1188 >>> stats.boxcox_normmax(x, method='all')
1189 array([2.23831866, 2.21756343])
1191 >>> fig = plt.figure()
1192 >>> ax = fig.add_subplot(111)
1193 >>> prob = stats.boxcox_normplot(x, -10, 10, plot=ax)
1194 >>> ax.axvline(lmax_mle, color='r')
1195 >>> ax.axvline(lmax_pearsonr, color='g', ls='--')
1197 >>> plt.show()
1199 Alternatively, we can define our own `optimizer` function. Suppose we
1200 are only interested in values of `lmbda` on the interval [6, 7], we
1201 want to use `scipy.optimize.minimize_scalar` with ``method='bounded'``,
1202 and we want to use tighter tolerances when optimizing the log-likelihood
1203 function. To do this, we define a function that accepts positional argument
1204 `fun` and uses `scipy.optimize.minimize_scalar` to minimize `fun` subject
1205 to the provided bounds and tolerances:
1207 >>> from scipy import optimize
1208 >>> options = {'xatol': 1e-12} # absolute tolerance on `x`
1209 >>> def optimizer(fun):
1210 ... return optimize.minimize_scalar(fun, bounds=(6, 7),
1211 ... method="bounded", options=options)
1212 >>> stats.boxcox_normmax(x, optimizer=optimizer)
1213 6.000...
1214 """
1215 # If optimizer is not given, define default 'brent' optimizer.
1216 if optimizer is None:
1218 # Set default value for `brack`.
1219 if brack is None:
1220 brack = (-2.0, 2.0)
1222 def _optimizer(func, args):
1223 return optimize.brent(func, args=args, brack=brack)
1225 # Otherwise check optimizer.
1226 else:
1227 if not callable(optimizer):
1228 raise ValueError("`optimizer` must be a callable")
1230 if brack is not None:
1231 raise ValueError("`brack` must be None if `optimizer` is given")
1233 # `optimizer` is expected to return a `OptimizeResult` object, we here
1234 # get the solution to the optimization problem.
1235 def _optimizer(func, args):
1236 def func_wrapped(x):
1237 return func(x, *args)
1238 return getattr(optimizer(func_wrapped), 'x', None)
1240 def _pearsonr(x):
1241 osm_uniform = _calc_uniform_order_statistic_medians(len(x))
1242 xvals = distributions.norm.ppf(osm_uniform)
1244 def _eval_pearsonr(lmbda, xvals, samps):
1245 # This function computes the x-axis values of the probability plot
1246 # and computes a linear regression (including the correlation) and
1247 # returns ``1 - r`` so that a minimization function maximizes the
1248 # correlation.
1249 y = boxcox(samps, lmbda)
1250 yvals = np.sort(y)
1251 r, prob = _stats_py.pearsonr(xvals, yvals)
1252 return 1 - r
1254 return _optimizer(_eval_pearsonr, args=(xvals, x))
1256 def _mle(x):
1257 def _eval_mle(lmb, data):
1258 # function to minimize
1259 return -boxcox_llf(lmb, data)
1261 return _optimizer(_eval_mle, args=(x,))
1263 def _all(x):
1264 maxlog = np.empty(2, dtype=float)
1265 maxlog[0] = _pearsonr(x)
1266 maxlog[1] = _mle(x)
1267 return maxlog
1269 methods = {'pearsonr': _pearsonr,
1270 'mle': _mle,
1271 'all': _all}
1272 if method not in methods.keys():
1273 raise ValueError("Method %s not recognized." % method)
1275 optimfunc = methods[method]
1276 res = optimfunc(x)
1277 if res is None:
1278 message = ("`optimizer` must return an object containing the optimal "
1279 "`lmbda` in attribute `x`")
1280 raise ValueError(message)
1281 return res
1284def _normplot(method, x, la, lb, plot=None, N=80):
1285 """Compute parameters for a Box-Cox or Yeo-Johnson normality plot,
1286 optionally show it.
1288 See `boxcox_normplot` or `yeojohnson_normplot` for details.
1289 """
1291 if method == 'boxcox':
1292 title = 'Box-Cox Normality Plot'
1293 transform_func = boxcox
1294 else:
1295 title = 'Yeo-Johnson Normality Plot'
1296 transform_func = yeojohnson
1298 x = np.asarray(x)
1299 if x.size == 0:
1300 return x
1302 if lb <= la:
1303 raise ValueError("`lb` has to be larger than `la`.")
1305 if method == 'boxcox' and np.any(x <= 0):
1306 raise ValueError("Data must be positive.")
1308 lmbdas = np.linspace(la, lb, num=N)
1309 ppcc = lmbdas * 0.0
1310 for i, val in enumerate(lmbdas):
1311 # Determine for each lmbda the square root of correlation coefficient
1312 # of transformed x
1313 z = transform_func(x, lmbda=val)
1314 _, (_, _, r) = probplot(z, dist='norm', fit=True)
1315 ppcc[i] = r
1317 if plot is not None:
1318 plot.plot(lmbdas, ppcc, 'x')
1319 _add_axis_labels_title(plot, xlabel='$\\lambda$',
1320 ylabel='Prob Plot Corr. Coef.',
1321 title=title)
1323 return lmbdas, ppcc
1326def boxcox_normplot(x, la, lb, plot=None, N=80):
1327 """Compute parameters for a Box-Cox normality plot, optionally show it.
1329 A Box-Cox normality plot shows graphically what the best transformation
1330 parameter is to use in `boxcox` to obtain a distribution that is close
1331 to normal.
1333 Parameters
1334 ----------
1335 x : array_like
1336 Input array.
1337 la, lb : scalar
1338 The lower and upper bounds for the ``lmbda`` values to pass to `boxcox`
1339 for Box-Cox transformations. These are also the limits of the
1340 horizontal axis of the plot if that is generated.
1341 plot : object, optional
1342 If given, plots the quantiles and least squares fit.
1343 `plot` is an object that has to have methods "plot" and "text".
1344 The `matplotlib.pyplot` module or a Matplotlib Axes object can be used,
1345 or a custom object with the same methods.
1346 Default is None, which means that no plot is created.
1347 N : int, optional
1348 Number of points on the horizontal axis (equally distributed from
1349 `la` to `lb`).
1351 Returns
1352 -------
1353 lmbdas : ndarray
1354 The ``lmbda`` values for which a Box-Cox transform was done.
1355 ppcc : ndarray
1356 Probability Plot Correlelation Coefficient, as obtained from `probplot`
1357 when fitting the Box-Cox transformed input `x` against a normal
1358 distribution.
1360 See Also
1361 --------
1362 probplot, boxcox, boxcox_normmax, boxcox_llf, ppcc_max
1364 Notes
1365 -----
1366 Even if `plot` is given, the figure is not shown or saved by
1367 `boxcox_normplot`; ``plt.show()`` or ``plt.savefig('figname.png')``
1368 should be used after calling `probplot`.
1370 Examples
1371 --------
1372 >>> from scipy import stats
1373 >>> import matplotlib.pyplot as plt
1375 Generate some non-normally distributed data, and create a Box-Cox plot:
1377 >>> x = stats.loggamma.rvs(5, size=500) + 5
1378 >>> fig = plt.figure()
1379 >>> ax = fig.add_subplot(111)
1380 >>> prob = stats.boxcox_normplot(x, -20, 20, plot=ax)
1382 Determine and plot the optimal ``lmbda`` to transform ``x`` and plot it in
1383 the same plot:
1385 >>> _, maxlog = stats.boxcox(x)
1386 >>> ax.axvline(maxlog, color='r')
1388 >>> plt.show()
1390 """
1391 return _normplot('boxcox', x, la, lb, plot, N)
1394def yeojohnson(x, lmbda=None):
1395 r"""Return a dataset transformed by a Yeo-Johnson power transformation.
1397 Parameters
1398 ----------
1399 x : ndarray
1400 Input array. Should be 1-dimensional.
1401 lmbda : float, optional
1402 If ``lmbda`` is ``None``, find the lambda that maximizes the
1403 log-likelihood function and return it as the second output argument.
1404 Otherwise the transformation is done for the given value.
1406 Returns
1407 -------
1408 yeojohnson: ndarray
1409 Yeo-Johnson power transformed array.
1410 maxlog : float, optional
1411 If the `lmbda` parameter is None, the second returned argument is
1412 the lambda that maximizes the log-likelihood function.
1414 See Also
1415 --------
1416 probplot, yeojohnson_normplot, yeojohnson_normmax, yeojohnson_llf, boxcox
1418 Notes
1419 -----
1420 The Yeo-Johnson transform is given by::
1422 y = ((x + 1)**lmbda - 1) / lmbda, for x >= 0, lmbda != 0
1423 log(x + 1), for x >= 0, lmbda = 0
1424 -((-x + 1)**(2 - lmbda) - 1) / (2 - lmbda), for x < 0, lmbda != 2
1425 -log(-x + 1), for x < 0, lmbda = 2
1427 Unlike `boxcox`, `yeojohnson` does not require the input data to be
1428 positive.
1430 .. versionadded:: 1.2.0
1433 References
1434 ----------
1435 I. Yeo and R.A. Johnson, "A New Family of Power Transformations to
1436 Improve Normality or Symmetry", Biometrika 87.4 (2000):
1439 Examples
1440 --------
1441 >>> from scipy import stats
1442 >>> import matplotlib.pyplot as plt
1444 We generate some random variates from a non-normal distribution and make a
1445 probability plot for it, to show it is non-normal in the tails:
1447 >>> fig = plt.figure()
1448 >>> ax1 = fig.add_subplot(211)
1449 >>> x = stats.loggamma.rvs(5, size=500) + 5
1450 >>> prob = stats.probplot(x, dist=stats.norm, plot=ax1)
1451 >>> ax1.set_xlabel('')
1452 >>> ax1.set_title('Probplot against normal distribution')
1454 We now use `yeojohnson` to transform the data so it's closest to normal:
1456 >>> ax2 = fig.add_subplot(212)
1457 >>> xt, lmbda = stats.yeojohnson(x)
1458 >>> prob = stats.probplot(xt, dist=stats.norm, plot=ax2)
1459 >>> ax2.set_title('Probplot after Yeo-Johnson transformation')
1461 >>> plt.show()
1463 """
1464 x = np.asarray(x)
1465 if x.size == 0:
1466 return x
1468 if np.issubdtype(x.dtype, np.complexfloating):
1469 raise ValueError('Yeo-Johnson transformation is not defined for '
1470 'complex numbers.')
1472 if np.issubdtype(x.dtype, np.integer):
1473 x = x.astype(np.float64, copy=False)
1475 if lmbda is not None:
1476 return _yeojohnson_transform(x, lmbda)
1478 # if lmbda=None, find the lmbda that maximizes the log-likelihood function.
1479 lmax = yeojohnson_normmax(x)
1480 y = _yeojohnson_transform(x, lmax)
1482 return y, lmax
1485def _yeojohnson_transform(x, lmbda):
1486 """Returns `x` transformed by the Yeo-Johnson power transform with given
1487 parameter `lmbda`.
1488 """
1489 out = np.zeros_like(x)
1490 pos = x >= 0 # binary mask
1492 # when x >= 0
1493 if abs(lmbda) < np.spacing(1.):
1494 out[pos] = np.log1p(x[pos])
1495 else: # lmbda != 0
1496 out[pos] = (np.power(x[pos] + 1, lmbda) - 1) / lmbda
1498 # when x < 0
1499 if abs(lmbda - 2) > np.spacing(1.):
1500 out[~pos] = -(np.power(-x[~pos] + 1, 2 - lmbda) - 1) / (2 - lmbda)
1501 else: # lmbda == 2
1502 out[~pos] = -np.log1p(-x[~pos])
1504 return out
1507def yeojohnson_llf(lmb, data):
1508 r"""The yeojohnson log-likelihood function.
1510 Parameters
1511 ----------
1512 lmb : scalar
1513 Parameter for Yeo-Johnson transformation. See `yeojohnson` for
1514 details.
1515 data : array_like
1516 Data to calculate Yeo-Johnson log-likelihood for. If `data` is
1517 multi-dimensional, the log-likelihood is calculated along the first
1518 axis.
1520 Returns
1521 -------
1522 llf : float
1523 Yeo-Johnson log-likelihood of `data` given `lmb`.
1525 See Also
1526 --------
1527 yeojohnson, probplot, yeojohnson_normplot, yeojohnson_normmax
1529 Notes
1530 -----
1531 The Yeo-Johnson log-likelihood function is defined here as
1533 .. math::
1535 llf = -N/2 \log(\hat{\sigma}^2) + (\lambda - 1)
1536 \sum_i \text{ sign }(x_i)\log(|x_i| + 1)
1538 where :math:`\hat{\sigma}^2` is estimated variance of the Yeo-Johnson
1539 transformed input data ``x``.
1541 .. versionadded:: 1.2.0
1543 Examples
1544 --------
1545 >>> import numpy as np
1546 >>> from scipy import stats
1547 >>> import matplotlib.pyplot as plt
1548 >>> from mpl_toolkits.axes_grid1.inset_locator import inset_axes
1550 Generate some random variates and calculate Yeo-Johnson log-likelihood
1551 values for them for a range of ``lmbda`` values:
1553 >>> x = stats.loggamma.rvs(5, loc=10, size=1000)
1554 >>> lmbdas = np.linspace(-2, 10)
1555 >>> llf = np.zeros(lmbdas.shape, dtype=float)
1556 >>> for ii, lmbda in enumerate(lmbdas):
1557 ... llf[ii] = stats.yeojohnson_llf(lmbda, x)
1559 Also find the optimal lmbda value with `yeojohnson`:
1561 >>> x_most_normal, lmbda_optimal = stats.yeojohnson(x)
1563 Plot the log-likelihood as function of lmbda. Add the optimal lmbda as a
1564 horizontal line to check that that's really the optimum:
1566 >>> fig = plt.figure()
1567 >>> ax = fig.add_subplot(111)
1568 >>> ax.plot(lmbdas, llf, 'b.-')
1569 >>> ax.axhline(stats.yeojohnson_llf(lmbda_optimal, x), color='r')
1570 >>> ax.set_xlabel('lmbda parameter')
1571 >>> ax.set_ylabel('Yeo-Johnson log-likelihood')
1573 Now add some probability plots to show that where the log-likelihood is
1574 maximized the data transformed with `yeojohnson` looks closest to normal:
1576 >>> locs = [3, 10, 4] # 'lower left', 'center', 'lower right'
1577 >>> for lmbda, loc in zip([-1, lmbda_optimal, 9], locs):
1578 ... xt = stats.yeojohnson(x, lmbda=lmbda)
1579 ... (osm, osr), (slope, intercept, r_sq) = stats.probplot(xt)
1580 ... ax_inset = inset_axes(ax, width="20%", height="20%", loc=loc)
1581 ... ax_inset.plot(osm, osr, 'c.', osm, slope*osm + intercept, 'k-')
1582 ... ax_inset.set_xticklabels([])
1583 ... ax_inset.set_yticklabels([])
1584 ... ax_inset.set_title(r'$\lambda=%1.2f$' % lmbda)
1586 >>> plt.show()
1588 """
1589 data = np.asarray(data)
1590 n_samples = data.shape[0]
1592 if n_samples == 0:
1593 return np.nan
1595 trans = _yeojohnson_transform(data, lmb)
1596 trans_var = trans.var(axis=0)
1597 loglike = np.empty_like(trans_var)
1599 # Avoid RuntimeWarning raised by np.log when the variance is too low
1600 tiny_variance = trans_var < np.finfo(trans_var.dtype).tiny
1601 loglike[tiny_variance] = np.inf
1603 loglike[~tiny_variance] = (
1604 -n_samples / 2 * np.log(trans_var[~tiny_variance]))
1605 loglike[~tiny_variance] += (
1606 (lmb - 1) * (np.sign(data) * np.log(np.abs(data) + 1)).sum(axis=0))
1607 return loglike
1610def yeojohnson_normmax(x, brack=(-2, 2)):
1611 """Compute optimal Yeo-Johnson transform parameter.
1613 Compute optimal Yeo-Johnson transform parameter for input data, using
1614 maximum likelihood estimation.
1616 Parameters
1617 ----------
1618 x : array_like
1619 Input array.
1620 brack : 2-tuple, optional
1621 The starting interval for a downhill bracket search with
1622 `optimize.brent`. Note that this is in most cases not critical; the
1623 final result is allowed to be outside this bracket.
1625 Returns
1626 -------
1627 maxlog : float
1628 The optimal transform parameter found.
1630 See Also
1631 --------
1632 yeojohnson, yeojohnson_llf, yeojohnson_normplot
1634 Notes
1635 -----
1636 .. versionadded:: 1.2.0
1638 Examples
1639 --------
1640 >>> import numpy as np
1641 >>> from scipy import stats
1642 >>> import matplotlib.pyplot as plt
1644 Generate some data and determine optimal ``lmbda``
1646 >>> rng = np.random.default_rng()
1647 >>> x = stats.loggamma.rvs(5, size=30, random_state=rng) + 5
1648 >>> lmax = stats.yeojohnson_normmax(x)
1650 >>> fig = plt.figure()
1651 >>> ax = fig.add_subplot(111)
1652 >>> prob = stats.yeojohnson_normplot(x, -10, 10, plot=ax)
1653 >>> ax.axvline(lmax, color='r')
1655 >>> plt.show()
1657 """
1658 def _neg_llf(lmbda, data):
1659 llf = yeojohnson_llf(lmbda, data)
1660 # reject likelihoods that are inf which are likely due to small
1661 # variance in the transformed space
1662 llf[np.isinf(llf)] = -np.inf
1663 return -llf
1665 with np.errstate(invalid='ignore'):
1666 return optimize.brent(_neg_llf, brack=brack, args=(x,))
1669def yeojohnson_normplot(x, la, lb, plot=None, N=80):
1670 """Compute parameters for a Yeo-Johnson normality plot, optionally show it.
1672 A Yeo-Johnson normality plot shows graphically what the best
1673 transformation parameter is to use in `yeojohnson` to obtain a
1674 distribution that is close to normal.
1676 Parameters
1677 ----------
1678 x : array_like
1679 Input array.
1680 la, lb : scalar
1681 The lower and upper bounds for the ``lmbda`` values to pass to
1682 `yeojohnson` for Yeo-Johnson transformations. These are also the
1683 limits of the horizontal axis of the plot if that is generated.
1684 plot : object, optional
1685 If given, plots the quantiles and least squares fit.
1686 `plot` is an object that has to have methods "plot" and "text".
1687 The `matplotlib.pyplot` module or a Matplotlib Axes object can be used,
1688 or a custom object with the same methods.
1689 Default is None, which means that no plot is created.
1690 N : int, optional
1691 Number of points on the horizontal axis (equally distributed from
1692 `la` to `lb`).
1694 Returns
1695 -------
1696 lmbdas : ndarray
1697 The ``lmbda`` values for which a Yeo-Johnson transform was done.
1698 ppcc : ndarray
1699 Probability Plot Correlelation Coefficient, as obtained from `probplot`
1700 when fitting the Box-Cox transformed input `x` against a normal
1701 distribution.
1703 See Also
1704 --------
1705 probplot, yeojohnson, yeojohnson_normmax, yeojohnson_llf, ppcc_max
1707 Notes
1708 -----
1709 Even if `plot` is given, the figure is not shown or saved by
1710 `boxcox_normplot`; ``plt.show()`` or ``plt.savefig('figname.png')``
1711 should be used after calling `probplot`.
1713 .. versionadded:: 1.2.0
1715 Examples
1716 --------
1717 >>> from scipy import stats
1718 >>> import matplotlib.pyplot as plt
1720 Generate some non-normally distributed data, and create a Yeo-Johnson plot:
1722 >>> x = stats.loggamma.rvs(5, size=500) + 5
1723 >>> fig = plt.figure()
1724 >>> ax = fig.add_subplot(111)
1725 >>> prob = stats.yeojohnson_normplot(x, -20, 20, plot=ax)
1727 Determine and plot the optimal ``lmbda`` to transform ``x`` and plot it in
1728 the same plot:
1730 >>> _, maxlog = stats.yeojohnson(x)
1731 >>> ax.axvline(maxlog, color='r')
1733 >>> plt.show()
1735 """
1736 return _normplot('yeojohnson', x, la, lb, plot, N)
1739ShapiroResult = namedtuple('ShapiroResult', ('statistic', 'pvalue'))
1742def shapiro(x):
1743 r"""Perform the Shapiro-Wilk test for normality.
1745 The Shapiro-Wilk test tests the null hypothesis that the
1746 data was drawn from a normal distribution.
1748 Parameters
1749 ----------
1750 x : array_like
1751 Array of sample data.
1753 Returns
1754 -------
1755 statistic : float
1756 The test statistic.
1757 p-value : float
1758 The p-value for the hypothesis test.
1760 See Also
1761 --------
1762 anderson : The Anderson-Darling test for normality
1763 kstest : The Kolmogorov-Smirnov test for goodness of fit.
1765 Notes
1766 -----
1767 The algorithm used is described in [4]_ but censoring parameters as
1768 described are not implemented. For N > 5000 the W test statistic is
1769 accurate, but the p-value may not be.
1771 References
1772 ----------
1773 .. [1] https://www.itl.nist.gov/div898/handbook/prc/section2/prc213.htm
1774 .. [2] Shapiro, S. S. & Wilk, M.B (1965). An analysis of variance test for
1775 normality (complete samples), Biometrika, Vol. 52, pp. 591-611.
1776 .. [3] Razali, N. M. & Wah, Y. B. (2011) Power comparisons of Shapiro-Wilk,
1777 Kolmogorov-Smirnov, Lilliefors and Anderson-Darling tests, Journal of
1778 Statistical Modeling and Analytics, Vol. 2, pp. 21-33.
1779 .. [4] ALGORITHM AS R94 APPL. STATIST. (1995) VOL. 44, NO. 4.
1780 .. [5] B. Phipson and G. K. Smyth. "Permutation P-values Should Never Be
1781 Zero: Calculating Exact P-values When Permutations Are Randomly
1782 Drawn." Statistical Applications in Genetics and Molecular Biology
1783 9.1 (2010).
1784 .. [6] Panagiotakos, D. B. (2008). The value of p-value in biomedical
1785 research. The open cardiovascular medicine journal, 2, 97.
1787 Examples
1788 --------
1789 Suppose we wish to infer from measurements whether the weights of adult
1790 human males in a medical study are not normally distributed [2]_.
1791 The weights (lbs) are recorded in the array ``x`` below.
1793 >>> import numpy as np
1794 >>> x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236])
1796 The normality test of [1]_ and [2]_ begins by computing a statistic based
1797 on the relationship between the observations and the expected order
1798 statistics of a normal distribution.
1800 >>> from scipy import stats
1801 >>> res = stats.shapiro(x)
1802 >>> res.statistic
1803 0.7888147830963135
1805 The value of this statistic tends to be high (close to 1) for samples drawn
1806 from a normal distribution.
1808 The test is performed by comparing the observed value of the statistic
1809 against the null distribution: the distribution of statistic values formed
1810 under the null hypothesis that the weights were drawn from a normal
1811 distribution. For this normality test, the null distribution is not easy to
1812 calculate exactly, so it is usually approximated by Monte Carlo methods,
1813 that is, drawing many samples of the same size as ``x`` from a normal
1814 distribution and computing the values of the statistic for each.
1816 >>> def statistic(x):
1817 ... # Get only the `shapiro` statistic; ignore its p-value
1818 ... return stats.shapiro(x).statistic
1819 >>> ref = stats.monte_carlo_test(x, stats.norm.rvs, statistic,
1820 ... alternative='less')
1821 >>> import matplotlib.pyplot as plt
1822 >>> fig, ax = plt.subplots(figsize=(8, 5))
1823 >>> bins = np.linspace(0.65, 1, 50)
1824 >>> def plot(ax): # we'll re-use this
1825 ... ax.hist(ref.null_distribution, density=True, bins=bins)
1826 ... ax.set_title("Shapiro-Wilk Test Null Distribution \n"
1827 ... "(Monte Carlo Approximation, 11 Observations)")
1828 ... ax.set_xlabel("statistic")
1829 ... ax.set_ylabel("probability density")
1830 >>> plot(ax)
1831 >>> plt.show()
1833 The comparison is quantified by the p-value: the proportion of values in
1834 the null distribution less than or equal to the observed value of the
1835 statistic.
1837 >>> fig, ax = plt.subplots(figsize=(8, 5))
1838 >>> plot(ax)
1839 >>> annotation = (f'p-value={res.pvalue:.6f}\n(highlighted area)')
1840 >>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
1841 >>> _ = ax.annotate(annotation, (0.75, 0.1), (0.68, 0.7), arrowprops=props)
1842 >>> i_extreme = np.where(bins <= res.statistic)[0]
1843 >>> for i in i_extreme:
1844 ... ax.patches[i].set_color('C1')
1845 >>> plt.xlim(0.65, 0.9)
1846 >>> plt.ylim(0, 4)
1847 >>> plt.show
1848 >>> res.pvalue
1849 0.006703833118081093
1851 If the p-value is "small" - that is, if there is a low probability of
1852 sampling data from a normally distributed population that produces such an
1853 extreme value of the statistic - this may be taken as evidence against
1854 the null hypothesis in favor of the alternative: the weights were not
1855 drawn from a normal distribution. Note that:
1857 - The inverse is not true; that is, the test is not used to provide
1858 evidence *for* the null hypothesis.
1859 - The threshold for values that will be considered "small" is a choice that
1860 should be made before the data is analyzed [5]_ with consideration of the
1861 risks of both false positives (incorrectly rejecting the null hypothesis)
1862 and false negatives (failure to reject a false null hypothesis).
1864 """
1865 x = np.ravel(x)
1867 N = len(x)
1868 if N < 3:
1869 raise ValueError("Data must be at least length 3.")
1871 x = x - np.median(x)
1873 a = zeros(N, 'f')
1874 init = 0
1876 y = sort(x)
1877 a, w, pw, ifault = _statlib.swilk(y, a[:N//2], init)
1878 if ifault not in [0, 2]:
1879 warnings.warn("Input data for shapiro has range zero. The results "
1880 "may not be accurate.")
1881 if N > 5000:
1882 warnings.warn("p-value may not be accurate for N > 5000.")
1884 # `swilk` can return negative p-values for N==3; see gh-18322.
1885 if N == 3:
1886 # Potential improvement: precision for small p-values
1887 pw = 1 - 6/np.pi*np.arccos(np.sqrt(w))
1888 return ShapiroResult(w, pw)
1891# Values from Stephens, M A, "EDF Statistics for Goodness of Fit and
1892# Some Comparisons", Journal of the American Statistical
1893# Association, Vol. 69, Issue 347, Sept. 1974, pp 730-737
1894_Avals_norm = array([0.576, 0.656, 0.787, 0.918, 1.092])
1895_Avals_expon = array([0.922, 1.078, 1.341, 1.606, 1.957])
1896# From Stephens, M A, "Goodness of Fit for the Extreme Value Distribution",
1897# Biometrika, Vol. 64, Issue 3, Dec. 1977, pp 583-588.
1898_Avals_gumbel = array([0.474, 0.637, 0.757, 0.877, 1.038])
1899# From Stephens, M A, "Tests of Fit for the Logistic Distribution Based
1900# on the Empirical Distribution Function.", Biometrika,
1901# Vol. 66, Issue 3, Dec. 1979, pp 591-595.
1902_Avals_logistic = array([0.426, 0.563, 0.660, 0.769, 0.906, 1.010])
1903# From Richard A. Lockhart and Michael A. Stephens "Estimation and Tests of
1904# Fit for the Three-Parameter Weibull Distribution"
1905# Journal of the Royal Statistical Society.Series B(Methodological)
1906# Vol. 56, No. 3 (1994), pp. 491-500, table 1. Keys are c*100
1907_Avals_weibull = [[0.292, 0.395, 0.467, 0.522, 0.617, 0.711, 0.836, 0.931],
1908 [0.295, 0.399, 0.471, 0.527, 0.623, 0.719, 0.845, 0.941],
1909 [0.298, 0.403, 0.476, 0.534, 0.631, 0.728, 0.856, 0.954],
1910 [0.301, 0.408, 0.483, 0.541, 0.640, 0.738, 0.869, 0.969],
1911 [0.305, 0.414, 0.490, 0.549, 0.650, 0.751, 0.885, 0.986],
1912 [0.309, 0.421, 0.498, 0.559, 0.662, 0.765, 0.902, 1.007],
1913 [0.314, 0.429, 0.508, 0.570, 0.676, 0.782, 0.923, 1.030],
1914 [0.320, 0.438, 0.519, 0.583, 0.692, 0.802, 0.947, 1.057],
1915 [0.327, 0.448, 0.532, 0.598, 0.711, 0.824, 0.974, 1.089],
1916 [0.334, 0.469, 0.547, 0.615, 0.732, 0.850, 1.006, 1.125],
1917 [0.342, 0.472, 0.563, 0.636, 0.757, 0.879, 1.043, 1.167]]
1918_Avals_weibull = np.array(_Avals_weibull)
1919_cvals_weibull = np.linspace(0, 0.5, 11)
1920_get_As_weibull = interpolate.interp1d(_cvals_weibull, _Avals_weibull.T,
1921 kind='linear', bounds_error=False,
1922 fill_value=_Avals_weibull[-1])
1925def _weibull_fit_check(params, x):
1926 # Refine the fit returned by `weibull_min.fit` to ensure that the first
1927 # order necessary conditions are satisfied. If not, raise an error.
1928 # Here, use `m` for the shape parameter to be consistent with [7]
1929 # and avoid confusion with `c` as defined in [7].
1930 n = len(x)
1931 m, u, s = params
1933 def dnllf_dm(m, u):
1934 # Partial w.r.t. shape w/ optimal scale. See [7] Equation 5.
1935 xu = x-u
1936 return (1/m - (xu**m*np.log(xu)).sum()/(xu**m).sum()
1937 + np.log(xu).sum()/n)
1939 def dnllf_du(m, u):
1940 # Partial w.r.t. loc w/ optimal scale. See [7] Equation 6.
1941 xu = x-u
1942 return (m-1)/m*(xu**-1).sum() - n*(xu**(m-1)).sum()/(xu**m).sum()
1944 def get_scale(m, u):
1945 # Partial w.r.t. scale solved in terms of shape and location.
1946 # See [7] Equation 7.
1947 return ((x-u)**m/n).sum()**(1/m)
1949 def dnllf(params):
1950 # Partial derivatives of the NLLF w.r.t. parameters, i.e.
1951 # first order necessary conditions for MLE fit.
1952 return [dnllf_dm(*params), dnllf_du(*params)]
1954 suggestion = ("Maximum likelihood estimation is known to be challenging "
1955 "for the three-parameter Weibull distribution. Consider "
1956 "performing a custom goodness-of-fit test using "
1957 "`scipy.stats.monte_carlo_test`.")
1959 if np.allclose(u, np.min(x)) or m < 1:
1960 # The critical values provided by [7] don't seem to control the
1961 # Type I error rate in this case. Error out.
1962 message = ("Maximum likelihood estimation has converged to "
1963 "a solution in which the location is equal to the minimum "
1964 "of the data, the shape parameter is less than 2, or both. "
1965 "The table of critical values in [7] does not "
1966 "include this case. " + suggestion)
1967 raise ValueError(message)
1969 try:
1970 # Refine the MLE / verify that first-order necessary conditions are
1971 # satisfied. If so, the critical values provided in [7] seem reliable.
1972 with np.errstate(over='raise', invalid='raise'):
1973 res = optimize.root(dnllf, params[:-1])
1975 message = ("Solution of MLE first-order conditions failed: "
1976 f"{res.message}. `anderson` cannot continue. " + suggestion)
1977 if not res.success:
1978 raise ValueError(message)
1980 except (FloatingPointError, ValueError) as e:
1981 message = ("An error occurred while fitting the Weibull distribution "
1982 "to the data, so `anderson` cannot continue. " + suggestion)
1983 raise ValueError(message) from e
1985 m, u = res.x
1986 s = get_scale(m, u)
1987 return m, u, s
1990AndersonResult = _make_tuple_bunch('AndersonResult',
1991 ['statistic', 'critical_values',
1992 'significance_level'], ['fit_result'])
1995def anderson(x, dist='norm'):
1996 """Anderson-Darling test for data coming from a particular distribution.
1998 The Anderson-Darling test tests the null hypothesis that a sample is
1999 drawn from a population that follows a particular distribution.
2000 For the Anderson-Darling test, the critical values depend on
2001 which distribution is being tested against. This function works
2002 for normal, exponential, logistic, weibull_min, or Gumbel (Extreme Value
2003 Type I) distributions.
2005 Parameters
2006 ----------
2007 x : array_like
2008 Array of sample data.
2009 dist : {'norm', 'expon', 'logistic', 'gumbel', 'gumbel_l', 'gumbel_r', 'extreme1', 'weibull_min'}, optional
2010 The type of distribution to test against. The default is 'norm'.
2011 The names 'extreme1', 'gumbel_l' and 'gumbel' are synonyms for the
2012 same distribution.
2014 Returns
2015 -------
2016 result : AndersonResult
2017 An object with the following attributes:
2019 statistic : float
2020 The Anderson-Darling test statistic.
2021 critical_values : list
2022 The critical values for this distribution.
2023 significance_level : list
2024 The significance levels for the corresponding critical values
2025 in percents. The function returns critical values for a
2026 differing set of significance levels depending on the
2027 distribution that is being tested against.
2028 fit_result : `~scipy.stats._result_classes.FitResult`
2029 An object containing the results of fitting the distribution to
2030 the data.
2032 See Also
2033 --------
2034 kstest : The Kolmogorov-Smirnov test for goodness-of-fit.
2036 Notes
2037 -----
2038 Critical values provided are for the following significance levels:
2040 normal/exponential
2041 15%, 10%, 5%, 2.5%, 1%
2042 logistic
2043 25%, 10%, 5%, 2.5%, 1%, 0.5%
2044 gumbel_l / gumbel_r
2045 25%, 10%, 5%, 2.5%, 1%
2046 weibull_min
2047 50%, 25%, 15%, 10%, 5%, 2.5%, 1%, 0.5%
2049 If the returned statistic is larger than these critical values then
2050 for the corresponding significance level, the null hypothesis that
2051 the data come from the chosen distribution can be rejected.
2052 The returned statistic is referred to as 'A2' in the references.
2054 For `weibull_min`, maximum likelihood estimation is known to be
2055 challenging. If the test returns successfully, then the first order
2056 conditions for a maximum likehood estimate have been verified and
2057 the critical values correspond relatively well to the significance levels,
2058 provided that the sample is sufficiently large (>10 observations [7]).
2059 However, for some data - especially data with no left tail - `anderson`
2060 is likely to result in an error message. In this case, consider
2061 performing a custom goodness of fit test using
2062 `scipy.stats.monte_carlo_test`.
2064 References
2065 ----------
2066 .. [1] https://www.itl.nist.gov/div898/handbook/prc/section2/prc213.htm
2067 .. [2] Stephens, M. A. (1974). EDF Statistics for Goodness of Fit and
2068 Some Comparisons, Journal of the American Statistical Association,
2069 Vol. 69, pp. 730-737.
2070 .. [3] Stephens, M. A. (1976). Asymptotic Results for Goodness-of-Fit
2071 Statistics with Unknown Parameters, Annals of Statistics, Vol. 4,
2072 pp. 357-369.
2073 .. [4] Stephens, M. A. (1977). Goodness of Fit for the Extreme Value
2074 Distribution, Biometrika, Vol. 64, pp. 583-588.
2075 .. [5] Stephens, M. A. (1977). Goodness of Fit with Special Reference
2076 to Tests for Exponentiality , Technical Report No. 262,
2077 Department of Statistics, Stanford University, Stanford, CA.
2078 .. [6] Stephens, M. A. (1979). Tests of Fit for the Logistic Distribution
2079 Based on the Empirical Distribution Function, Biometrika, Vol. 66,
2080 pp. 591-595.
2081 .. [7] Richard A. Lockhart and Michael A. Stephens "Estimation and Tests of
2082 Fit for the Three-Parameter Weibull Distribution"
2083 Journal of the Royal Statistical Society.Series B(Methodological)
2084 Vol. 56, No. 3 (1994), pp. 491-500, Table 0.
2086 Examples
2087 --------
2088 Test the null hypothesis that a random sample was drawn from a normal
2089 distribution (with unspecified mean and standard deviation).
2091 >>> import numpy as np
2092 >>> from scipy.stats import anderson
2093 >>> rng = np.random.default_rng()
2094 >>> data = rng.random(size=35)
2095 >>> res = anderson(data)
2096 >>> res.statistic
2097 0.8398018749744764
2098 >>> res.critical_values
2099 array([0.527, 0.6 , 0.719, 0.839, 0.998])
2100 >>> res.significance_level
2101 array([15. , 10. , 5. , 2.5, 1. ])
2103 The value of the statistic (barely) exceeds the critical value associated
2104 with a significance level of 2.5%, so the null hypothesis may be rejected
2105 at a significance level of 2.5%, but not at a significance level of 1%.
2107 """ # noqa
2108 dist = dist.lower()
2109 if dist in {'extreme1', 'gumbel'}:
2110 dist = 'gumbel_l'
2111 dists = {'norm', 'expon', 'gumbel_l',
2112 'gumbel_r', 'logistic', 'weibull_min'}
2114 if dist not in dists:
2115 raise ValueError(f"Invalid distribution; dist must be in {dists}.")
2116 y = sort(x)
2117 xbar = np.mean(x, axis=0)
2118 N = len(y)
2119 if dist == 'norm':
2120 s = np.std(x, ddof=1, axis=0)
2121 w = (y - xbar) / s
2122 fit_params = xbar, s
2123 logcdf = distributions.norm.logcdf(w)
2124 logsf = distributions.norm.logsf(w)
2125 sig = array([15, 10, 5, 2.5, 1])
2126 critical = around(_Avals_norm / (1.0 + 4.0/N - 25.0/N/N), 3)
2127 elif dist == 'expon':
2128 w = y / xbar
2129 fit_params = 0, xbar
2130 logcdf = distributions.expon.logcdf(w)
2131 logsf = distributions.expon.logsf(w)
2132 sig = array([15, 10, 5, 2.5, 1])
2133 critical = around(_Avals_expon / (1.0 + 0.6/N), 3)
2134 elif dist == 'logistic':
2135 def rootfunc(ab, xj, N):
2136 a, b = ab
2137 tmp = (xj - a) / b
2138 tmp2 = exp(tmp)
2139 val = [np.sum(1.0/(1+tmp2), axis=0) - 0.5*N,
2140 np.sum(tmp*(1.0-tmp2)/(1+tmp2), axis=0) + N]
2141 return array(val)
2143 sol0 = array([xbar, np.std(x, ddof=1, axis=0)])
2144 sol = optimize.fsolve(rootfunc, sol0, args=(x, N), xtol=1e-5)
2145 w = (y - sol[0]) / sol[1]
2146 fit_params = sol
2147 logcdf = distributions.logistic.logcdf(w)
2148 logsf = distributions.logistic.logsf(w)
2149 sig = array([25, 10, 5, 2.5, 1, 0.5])
2150 critical = around(_Avals_logistic / (1.0 + 0.25/N), 3)
2151 elif dist == 'gumbel_r':
2152 xbar, s = distributions.gumbel_r.fit(x)
2153 w = (y - xbar) / s
2154 fit_params = xbar, s
2155 logcdf = distributions.gumbel_r.logcdf(w)
2156 logsf = distributions.gumbel_r.logsf(w)
2157 sig = array([25, 10, 5, 2.5, 1])
2158 critical = around(_Avals_gumbel / (1.0 + 0.2/sqrt(N)), 3)
2159 elif dist == 'gumbel_l':
2160 xbar, s = distributions.gumbel_l.fit(x)
2161 w = (y - xbar) / s
2162 fit_params = xbar, s
2163 logcdf = distributions.gumbel_l.logcdf(w)
2164 logsf = distributions.gumbel_l.logsf(w)
2165 sig = array([25, 10, 5, 2.5, 1])
2166 critical = around(_Avals_gumbel / (1.0 + 0.2/sqrt(N)), 3)
2167 elif dist == 'weibull_min':
2168 message = ("Critical values of the test statistic are given for the "
2169 "asymptotic distribution. These may not be accurate for "
2170 "samples with fewer than 10 observations. Consider using "
2171 "`scipy.stats.monte_carlo_test`.")
2172 if N < 10:
2173 warnings.warn(message, stacklevel=2)
2174 # [7] writes our 'c' as 'm', and they write `c = 1/m`. Use their names.
2175 m, loc, scale = distributions.weibull_min.fit(y)
2176 m, loc, scale = _weibull_fit_check((m, loc, scale), y)
2177 fit_params = m, loc, scale
2178 logcdf = stats.weibull_min(*fit_params).logcdf(y)
2179 logsf = stats.weibull_min(*fit_params).logsf(y)
2180 c = 1 / m # m and c are as used in [7]
2181 sig = array([0.5, 0.75, 0.85, 0.9, 0.95, 0.975, 0.99, 0.995])
2182 critical = _get_As_weibull(c)
2183 # Goodness-of-fit tests should only be used to provide evidence
2184 # _against_ the null hypothesis. Be conservative and round up.
2185 critical = np.round(critical + 0.0005, decimals=3)
2187 i = arange(1, N + 1)
2188 A2 = -N - np.sum((2*i - 1.0) / N * (logcdf + logsf[::-1]), axis=0)
2190 # FitResult initializer expects an optimize result, so let's work with it
2191 message = '`anderson` successfully fit the distribution to the data.'
2192 res = optimize.OptimizeResult(success=True, message=message)
2193 res.x = np.array(fit_params)
2194 fit_result = FitResult(getattr(distributions, dist), y,
2195 discrete=False, res=res)
2197 return AndersonResult(A2, critical, sig, fit_result=fit_result)
2200def _anderson_ksamp_midrank(samples, Z, Zstar, k, n, N):
2201 """Compute A2akN equation 7 of Scholz and Stephens.
2203 Parameters
2204 ----------
2205 samples : sequence of 1-D array_like
2206 Array of sample arrays.
2207 Z : array_like
2208 Sorted array of all observations.
2209 Zstar : array_like
2210 Sorted array of unique observations.
2211 k : int
2212 Number of samples.
2213 n : array_like
2214 Number of observations in each sample.
2215 N : int
2216 Total number of observations.
2218 Returns
2219 -------
2220 A2aKN : float
2221 The A2aKN statistics of Scholz and Stephens 1987.
2223 """
2224 A2akN = 0.
2225 Z_ssorted_left = Z.searchsorted(Zstar, 'left')
2226 if N == Zstar.size:
2227 lj = 1.
2228 else:
2229 lj = Z.searchsorted(Zstar, 'right') - Z_ssorted_left
2230 Bj = Z_ssorted_left + lj / 2.
2231 for i in arange(0, k):
2232 s = np.sort(samples[i])
2233 s_ssorted_right = s.searchsorted(Zstar, side='right')
2234 Mij = s_ssorted_right.astype(float)
2235 fij = s_ssorted_right - s.searchsorted(Zstar, 'left')
2236 Mij -= fij / 2.
2237 inner = lj / float(N) * (N*Mij - Bj*n[i])**2 / (Bj*(N - Bj) - N*lj/4.)
2238 A2akN += inner.sum() / n[i]
2239 A2akN *= (N - 1.) / N
2240 return A2akN
2243def _anderson_ksamp_right(samples, Z, Zstar, k, n, N):
2244 """Compute A2akN equation 6 of Scholz & Stephens.
2246 Parameters
2247 ----------
2248 samples : sequence of 1-D array_like
2249 Array of sample arrays.
2250 Z : array_like
2251 Sorted array of all observations.
2252 Zstar : array_like
2253 Sorted array of unique observations.
2254 k : int
2255 Number of samples.
2256 n : array_like
2257 Number of observations in each sample.
2258 N : int
2259 Total number of observations.
2261 Returns
2262 -------
2263 A2KN : float
2264 The A2KN statistics of Scholz and Stephens 1987.
2266 """
2267 A2kN = 0.
2268 lj = Z.searchsorted(Zstar[:-1], 'right') - Z.searchsorted(Zstar[:-1],
2269 'left')
2270 Bj = lj.cumsum()
2271 for i in arange(0, k):
2272 s = np.sort(samples[i])
2273 Mij = s.searchsorted(Zstar[:-1], side='right')
2274 inner = lj / float(N) * (N * Mij - Bj * n[i])**2 / (Bj * (N - Bj))
2275 A2kN += inner.sum() / n[i]
2276 return A2kN
2279Anderson_ksampResult = _make_tuple_bunch(
2280 'Anderson_ksampResult',
2281 ['statistic', 'critical_values', 'pvalue'], []
2282)
2285def anderson_ksamp(samples, midrank=True):
2286 """The Anderson-Darling test for k-samples.
2288 The k-sample Anderson-Darling test is a modification of the
2289 one-sample Anderson-Darling test. It tests the null hypothesis
2290 that k-samples are drawn from the same population without having
2291 to specify the distribution function of that population. The
2292 critical values depend on the number of samples.
2294 Parameters
2295 ----------
2296 samples : sequence of 1-D array_like
2297 Array of sample data in arrays.
2298 midrank : bool, optional
2299 Type of Anderson-Darling test which is computed. Default
2300 (True) is the midrank test applicable to continuous and
2301 discrete populations. If False, the right side empirical
2302 distribution is used.
2304 Returns
2305 -------
2306 res : Anderson_ksampResult
2307 An object containing attributes:
2309 statistic : float
2310 Normalized k-sample Anderson-Darling test statistic.
2311 critical_values : array
2312 The critical values for significance levels 25%, 10%, 5%, 2.5%, 1%,
2313 0.5%, 0.1%.
2314 pvalue : float
2315 The approximate p-value of the test. The value is floored / capped
2316 at 0.1% / 25%.
2318 Raises
2319 ------
2320 ValueError
2321 If less than 2 samples are provided, a sample is empty, or no
2322 distinct observations are in the samples.
2324 See Also
2325 --------
2326 ks_2samp : 2 sample Kolmogorov-Smirnov test
2327 anderson : 1 sample Anderson-Darling test
2329 Notes
2330 -----
2331 [1]_ defines three versions of the k-sample Anderson-Darling test:
2332 one for continuous distributions and two for discrete
2333 distributions, in which ties between samples may occur. The
2334 default of this routine is to compute the version based on the
2335 midrank empirical distribution function. This test is applicable
2336 to continuous and discrete data. If midrank is set to False, the
2337 right side empirical distribution is used for a test for discrete
2338 data. According to [1]_, the two discrete test statistics differ
2339 only slightly if a few collisions due to round-off errors occur in
2340 the test not adjusted for ties between samples.
2342 The critical values corresponding to the significance levels from 0.01
2343 to 0.25 are taken from [1]_. p-values are floored / capped
2344 at 0.1% / 25%. Since the range of critical values might be extended in
2345 future releases, it is recommended not to test ``p == 0.25``, but rather
2346 ``p >= 0.25`` (analogously for the lower bound).
2348 .. versionadded:: 0.14.0
2350 References
2351 ----------
2352 .. [1] Scholz, F. W and Stephens, M. A. (1987), K-Sample
2353 Anderson-Darling Tests, Journal of the American Statistical
2354 Association, Vol. 82, pp. 918-924.
2356 Examples
2357 --------
2358 >>> import numpy as np
2359 >>> from scipy import stats
2360 >>> rng = np.random.default_rng()
2361 >>> res = stats.anderson_ksamp([rng.normal(size=50),
2362 ... rng.normal(loc=0.5, size=30)])
2363 >>> res.statistic, res.pvalue
2364 (1.974403288713695, 0.04991293614572478)
2365 >>> res.critical_values
2366 array([0.325, 1.226, 1.961, 2.718, 3.752, 4.592, 6.546])
2368 The null hypothesis that the two random samples come from the same
2369 distribution can be rejected at the 5% level because the returned
2370 test value is greater than the critical value for 5% (1.961) but
2371 not at the 2.5% level. The interpolation gives an approximate
2372 p-value of 4.99%.
2374 >>> res = stats.anderson_ksamp([rng.normal(size=50),
2375 ... rng.normal(size=30), rng.normal(size=20)])
2376 >>> res.statistic, res.pvalue
2377 (-0.29103725200789504, 0.25)
2378 >>> res.critical_values
2379 array([ 0.44925884, 1.3052767 , 1.9434184 , 2.57696569, 3.41634856,
2380 4.07210043, 5.56419101])
2382 The null hypothesis cannot be rejected for three samples from an
2383 identical distribution. The reported p-value (25%) has been capped and
2384 may not be very accurate (since it corresponds to the value 0.449
2385 whereas the statistic is -0.291).
2387 """
2388 k = len(samples)
2389 if (k < 2):
2390 raise ValueError("anderson_ksamp needs at least two samples")
2392 samples = list(map(np.asarray, samples))
2393 Z = np.sort(np.hstack(samples))
2394 N = Z.size
2395 Zstar = np.unique(Z)
2396 if Zstar.size < 2:
2397 raise ValueError("anderson_ksamp needs more than one distinct "
2398 "observation")
2400 n = np.array([sample.size for sample in samples])
2401 if np.any(n == 0):
2402 raise ValueError("anderson_ksamp encountered sample without "
2403 "observations")
2405 if midrank:
2406 A2kN = _anderson_ksamp_midrank(samples, Z, Zstar, k, n, N)
2407 else:
2408 A2kN = _anderson_ksamp_right(samples, Z, Zstar, k, n, N)
2410 H = (1. / n).sum()
2411 hs_cs = (1. / arange(N - 1, 1, -1)).cumsum()
2412 h = hs_cs[-1] + 1
2413 g = (hs_cs / arange(2, N)).sum()
2415 a = (4*g - 6) * (k - 1) + (10 - 6*g)*H
2416 b = (2*g - 4)*k**2 + 8*h*k + (2*g - 14*h - 4)*H - 8*h + 4*g - 6
2417 c = (6*h + 2*g - 2)*k**2 + (4*h - 4*g + 6)*k + (2*h - 6)*H + 4*h
2418 d = (2*h + 6)*k**2 - 4*h*k
2419 sigmasq = (a*N**3 + b*N**2 + c*N + d) / ((N - 1.) * (N - 2.) * (N - 3.))
2420 m = k - 1
2421 A2 = (A2kN - m) / math.sqrt(sigmasq)
2423 # The b_i values are the interpolation coefficients from Table 2
2424 # of Scholz and Stephens 1987
2425 b0 = np.array([0.675, 1.281, 1.645, 1.96, 2.326, 2.573, 3.085])
2426 b1 = np.array([-0.245, 0.25, 0.678, 1.149, 1.822, 2.364, 3.615])
2427 b2 = np.array([-0.105, -0.305, -0.362, -0.391, -0.396, -0.345, -0.154])
2428 critical = b0 + b1 / math.sqrt(m) + b2 / m
2430 sig = np.array([0.25, 0.1, 0.05, 0.025, 0.01, 0.005, 0.001])
2431 if A2 < critical.min():
2432 p = sig.max()
2433 warnings.warn("p-value capped: true value larger than {}".format(p),
2434 stacklevel=2)
2435 elif A2 > critical.max():
2436 p = sig.min()
2437 warnings.warn("p-value floored: true value smaller than {}".format(p),
2438 stacklevel=2)
2439 else:
2440 # interpolation of probit of significance level
2441 pf = np.polyfit(critical, log(sig), 2)
2442 p = math.exp(np.polyval(pf, A2))
2444 # create result object with alias for backward compatibility
2445 res = Anderson_ksampResult(A2, critical, p)
2446 res.significance_level = p
2447 return res
2450AnsariResult = namedtuple('AnsariResult', ('statistic', 'pvalue'))
2453class _ABW:
2454 """Distribution of Ansari-Bradley W-statistic under the null hypothesis."""
2455 # TODO: calculate exact distribution considering ties
2456 # We could avoid summing over more than half the frequencies,
2457 # but inititally it doesn't seem worth the extra complexity
2459 def __init__(self):
2460 """Minimal initializer."""
2461 self.m = None
2462 self.n = None
2463 self.astart = None
2464 self.total = None
2465 self.freqs = None
2467 def _recalc(self, n, m):
2468 """When necessary, recalculate exact distribution."""
2469 if n != self.n or m != self.m:
2470 self.n, self.m = n, m
2471 # distribution is NOT symmetric when m + n is odd
2472 # n is len(x), m is len(y), and ratio of scales is defined x/y
2473 astart, a1, _ = _statlib.gscale(n, m)
2474 self.astart = astart # minimum value of statistic
2475 # Exact distribution of test statistic under null hypothesis
2476 # expressed as frequencies/counts/integers to maintain precision.
2477 # Stored as floats to avoid overflow of sums.
2478 self.freqs = a1.astype(np.float64)
2479 self.total = self.freqs.sum() # could calculate from m and n
2480 # probability mass is self.freqs / self.total;
2482 def pmf(self, k, n, m):
2483 """Probability mass function."""
2484 self._recalc(n, m)
2485 # The convention here is that PMF at k = 12.5 is the same as at k = 12,
2486 # -> use `floor` in case of ties.
2487 ind = np.floor(k - self.astart).astype(int)
2488 return self.freqs[ind] / self.total
2490 def cdf(self, k, n, m):
2491 """Cumulative distribution function."""
2492 self._recalc(n, m)
2493 # Null distribution derived without considering ties is
2494 # approximate. Round down to avoid Type I error.
2495 ind = np.ceil(k - self.astart).astype(int)
2496 return self.freqs[:ind+1].sum() / self.total
2498 def sf(self, k, n, m):
2499 """Survival function."""
2500 self._recalc(n, m)
2501 # Null distribution derived without considering ties is
2502 # approximate. Round down to avoid Type I error.
2503 ind = np.floor(k - self.astart).astype(int)
2504 return self.freqs[ind:].sum() / self.total
2507# Maintain state for faster repeat calls to ansari w/ method='exact'
2508_abw_state = _ABW()
2511def ansari(x, y, alternative='two-sided'):
2512 """Perform the Ansari-Bradley test for equal scale parameters.
2514 The Ansari-Bradley test ([1]_, [2]_) is a non-parametric test
2515 for the equality of the scale parameter of the distributions
2516 from which two samples were drawn. The null hypothesis states that
2517 the ratio of the scale of the distribution underlying `x` to the scale
2518 of the distribution underlying `y` is 1.
2520 Parameters
2521 ----------
2522 x, y : array_like
2523 Arrays of sample data.
2524 alternative : {'two-sided', 'less', 'greater'}, optional
2525 Defines the alternative hypothesis. Default is 'two-sided'.
2526 The following options are available:
2528 * 'two-sided': the ratio of scales is not equal to 1.
2529 * 'less': the ratio of scales is less than 1.
2530 * 'greater': the ratio of scales is greater than 1.
2532 .. versionadded:: 1.7.0
2534 Returns
2535 -------
2536 statistic : float
2537 The Ansari-Bradley test statistic.
2538 pvalue : float
2539 The p-value of the hypothesis test.
2541 See Also
2542 --------
2543 fligner : A non-parametric test for the equality of k variances
2544 mood : A non-parametric test for the equality of two scale parameters
2546 Notes
2547 -----
2548 The p-value given is exact when the sample sizes are both less than
2549 55 and there are no ties, otherwise a normal approximation for the
2550 p-value is used.
2552 References
2553 ----------
2554 .. [1] Ansari, A. R. and Bradley, R. A. (1960) Rank-sum tests for
2555 dispersions, Annals of Mathematical Statistics, 31, 1174-1189.
2556 .. [2] Sprent, Peter and N.C. Smeeton. Applied nonparametric
2557 statistical methods. 3rd ed. Chapman and Hall/CRC. 2001.
2558 Section 5.8.2.
2559 .. [3] Nathaniel E. Helwig "Nonparametric Dispersion and Equality
2560 Tests" at http://users.stat.umn.edu/~helwig/notes/npde-Notes.pdf
2562 Examples
2563 --------
2564 >>> import numpy as np
2565 >>> from scipy.stats import ansari
2566 >>> rng = np.random.default_rng()
2568 For these examples, we'll create three random data sets. The first
2569 two, with sizes 35 and 25, are drawn from a normal distribution with
2570 mean 0 and standard deviation 2. The third data set has size 25 and
2571 is drawn from a normal distribution with standard deviation 1.25.
2573 >>> x1 = rng.normal(loc=0, scale=2, size=35)
2574 >>> x2 = rng.normal(loc=0, scale=2, size=25)
2575 >>> x3 = rng.normal(loc=0, scale=1.25, size=25)
2577 First we apply `ansari` to `x1` and `x2`. These samples are drawn
2578 from the same distribution, so we expect the Ansari-Bradley test
2579 should not lead us to conclude that the scales of the distributions
2580 are different.
2582 >>> ansari(x1, x2)
2583 AnsariResult(statistic=541.0, pvalue=0.9762532927399098)
2585 With a p-value close to 1, we cannot conclude that there is a
2586 significant difference in the scales (as expected).
2588 Now apply the test to `x1` and `x3`:
2590 >>> ansari(x1, x3)
2591 AnsariResult(statistic=425.0, pvalue=0.0003087020407974518)
2593 The probability of observing such an extreme value of the statistic
2594 under the null hypothesis of equal scales is only 0.03087%. We take this
2595 as evidence against the null hypothesis in favor of the alternative:
2596 the scales of the distributions from which the samples were drawn
2597 are not equal.
2599 We can use the `alternative` parameter to perform a one-tailed test.
2600 In the above example, the scale of `x1` is greater than `x3` and so
2601 the ratio of scales of `x1` and `x3` is greater than 1. This means
2602 that the p-value when ``alternative='greater'`` should be near 0 and
2603 hence we should be able to reject the null hypothesis:
2605 >>> ansari(x1, x3, alternative='greater')
2606 AnsariResult(statistic=425.0, pvalue=0.0001543510203987259)
2608 As we can see, the p-value is indeed quite low. Use of
2609 ``alternative='less'`` should thus yield a large p-value:
2611 >>> ansari(x1, x3, alternative='less')
2612 AnsariResult(statistic=425.0, pvalue=0.9998643258449039)
2614 """
2615 if alternative not in {'two-sided', 'greater', 'less'}:
2616 raise ValueError("'alternative' must be 'two-sided',"
2617 " 'greater', or 'less'.")
2618 x, y = asarray(x), asarray(y)
2619 n = len(x)
2620 m = len(y)
2621 if m < 1:
2622 raise ValueError("Not enough other observations.")
2623 if n < 1:
2624 raise ValueError("Not enough test observations.")
2626 N = m + n
2627 xy = r_[x, y] # combine
2628 rank = _stats_py.rankdata(xy)
2629 symrank = amin(array((rank, N - rank + 1)), 0)
2630 AB = np.sum(symrank[:n], axis=0)
2631 uxy = unique(xy)
2632 repeats = (len(uxy) != len(xy))
2633 exact = ((m < 55) and (n < 55) and not repeats)
2634 if repeats and (m < 55 or n < 55):
2635 warnings.warn("Ties preclude use of exact statistic.")
2636 if exact:
2637 if alternative == 'two-sided':
2638 pval = 2.0 * np.minimum(_abw_state.cdf(AB, n, m),
2639 _abw_state.sf(AB, n, m))
2640 elif alternative == 'greater':
2641 # AB statistic is _smaller_ when ratio of scales is larger,
2642 # so this is the opposite of the usual calculation
2643 pval = _abw_state.cdf(AB, n, m)
2644 else:
2645 pval = _abw_state.sf(AB, n, m)
2646 return AnsariResult(AB, min(1.0, pval))
2648 # otherwise compute normal approximation
2649 if N % 2: # N odd
2650 mnAB = n * (N+1.0)**2 / 4.0 / N
2651 varAB = n * m * (N+1.0) * (3+N**2) / (48.0 * N**2)
2652 else:
2653 mnAB = n * (N+2.0) / 4.0
2654 varAB = m * n * (N+2) * (N-2.0) / 48 / (N-1.0)
2655 if repeats: # adjust variance estimates
2656 # compute np.sum(tj * rj**2,axis=0)
2657 fac = np.sum(symrank**2, axis=0)
2658 if N % 2: # N odd
2659 varAB = m * n * (16*N*fac - (N+1)**4) / (16.0 * N**2 * (N-1))
2660 else: # N even
2661 varAB = m * n * (16*fac - N*(N+2)**2) / (16.0 * N * (N-1))
2663 # Small values of AB indicate larger dispersion for the x sample.
2664 # Large values of AB indicate larger dispersion for the y sample.
2665 # This is opposite to the way we define the ratio of scales. see [1]_.
2666 z = (mnAB - AB) / sqrt(varAB)
2667 z, pval = _normtest_finish(z, alternative)
2668 return AnsariResult(AB, pval)
2671BartlettResult = namedtuple('BartlettResult', ('statistic', 'pvalue'))
2674def bartlett(*samples):
2675 r"""Perform Bartlett's test for equal variances.
2677 Bartlett's test tests the null hypothesis that all input samples
2678 are from populations with equal variances. For samples
2679 from significantly non-normal populations, Levene's test
2680 `levene` is more robust.
2682 Parameters
2683 ----------
2684 sample1, sample2, ... : array_like
2685 arrays of sample data. Only 1d arrays are accepted, they may have
2686 different lengths.
2688 Returns
2689 -------
2690 statistic : float
2691 The test statistic.
2692 pvalue : float
2693 The p-value of the test.
2695 See Also
2696 --------
2697 fligner : A non-parametric test for the equality of k variances
2698 levene : A robust parametric test for equality of k variances
2700 Notes
2701 -----
2702 Conover et al. (1981) examine many of the existing parametric and
2703 nonparametric tests by extensive simulations and they conclude that the
2704 tests proposed by Fligner and Killeen (1976) and Levene (1960) appear to be
2705 superior in terms of robustness of departures from normality and power
2706 ([3]_).
2708 References
2709 ----------
2710 .. [1] https://www.itl.nist.gov/div898/handbook/eda/section3/eda357.htm
2711 .. [2] Snedecor, George W. and Cochran, William G. (1989), Statistical
2712 Methods, Eighth Edition, Iowa State University Press.
2713 .. [3] Park, C. and Lindsay, B. G. (1999). Robust Scale Estimation and
2714 Hypothesis Testing based on Quadratic Inference Function. Technical
2715 Report #99-03, Center for Likelihood Studies, Pennsylvania State
2716 University.
2717 .. [4] Bartlett, M. S. (1937). Properties of Sufficiency and Statistical
2718 Tests. Proceedings of the Royal Society of London. Series A,
2719 Mathematical and Physical Sciences, Vol. 160, No.901, pp. 268-282.
2720 .. [5] C.I. BLISS (1952), The Statistics of Bioassay: With Special
2721 Reference to the Vitamins, pp 499-503,
2722 :doi:`10.1016/C2013-0-12584-6`.
2723 .. [6] B. Phipson and G. K. Smyth. "Permutation P-values Should Never Be
2724 Zero: Calculating Exact P-values When Permutations Are Randomly
2725 Drawn." Statistical Applications in Genetics and Molecular Biology
2726 9.1 (2010).
2727 .. [7] Ludbrook, J., & Dudley, H. (1998). Why permutation tests are
2728 superior to t and F tests in biomedical research. The American
2729 Statistician, 52(2), 127-132.
2731 Examples
2732 --------
2733 In [5]_, the influence of vitamin C on the tooth growth of guinea pigs
2734 was investigated. In a control study, 60 subjects were divided into
2735 small dose, medium dose, and large dose groups that received
2736 daily doses of 0.5, 1.0 and 2.0 mg of vitamin C, respectively.
2737 After 42 days, the tooth growth was measured.
2739 The ``small_dose``, ``medium_dose``, and ``large_dose`` arrays below record
2740 tooth growth measurements of the three groups in microns.
2742 >>> import numpy as np
2743 >>> small_dose = np.array([
2744 ... 4.2, 11.5, 7.3, 5.8, 6.4, 10, 11.2, 11.2, 5.2, 7,
2745 ... 15.2, 21.5, 17.6, 9.7, 14.5, 10, 8.2, 9.4, 16.5, 9.7
2746 ... ])
2747 >>> medium_dose = np.array([
2748 ... 16.5, 16.5, 15.2, 17.3, 22.5, 17.3, 13.6, 14.5, 18.8, 15.5,
2749 ... 19.7, 23.3, 23.6, 26.4, 20, 25.2, 25.8, 21.2, 14.5, 27.3
2750 ... ])
2751 >>> large_dose = np.array([
2752 ... 23.6, 18.5, 33.9, 25.5, 26.4, 32.5, 26.7, 21.5, 23.3, 29.5,
2753 ... 25.5, 26.4, 22.4, 24.5, 24.8, 30.9, 26.4, 27.3, 29.4, 23
2754 ... ])
2756 The `bartlett` statistic is sensitive to differences in variances
2757 between the samples.
2759 >>> from scipy import stats
2760 >>> res = stats.bartlett(small_dose, medium_dose, large_dose)
2761 >>> res.statistic
2762 0.6654670663030519
2764 The value of the statistic tends to be high when there is a large
2765 difference in variances.
2767 We can test for inequality of variance among the groups by comparing the
2768 observed value of the statistic against the null distribution: the
2769 distribution of statistic values derived under the null hypothesis that
2770 the population variances of the three groups are equal.
2772 For this test, the null distribution follows the chi-square distribution
2773 as shown below.
2775 >>> import matplotlib.pyplot as plt
2776 >>> k = 3 # number of samples
2777 >>> dist = stats.chi2(df=k-1)
2778 >>> val = np.linspace(0, 5, 100)
2779 >>> pdf = dist.pdf(val)
2780 >>> fig, ax = plt.subplots(figsize=(8, 5))
2781 >>> def plot(ax): # we'll re-use this
2782 ... ax.plot(val, pdf, color='C0')
2783 ... ax.set_title("Bartlett Test Null Distribution")
2784 ... ax.set_xlabel("statistic")
2785 ... ax.set_ylabel("probability density")
2786 ... ax.set_xlim(0, 5)
2787 ... ax.set_ylim(0, 1)
2788 >>> plot(ax)
2789 >>> plt.show()
2791 The comparison is quantified by the p-value: the proportion of values in
2792 the null distribution greater than or equal to the observed value of the
2793 statistic.
2795 >>> fig, ax = plt.subplots(figsize=(8, 5))
2796 >>> plot(ax)
2797 >>> pvalue = dist.sf(res.statistic)
2798 >>> annotation = (f'p-value={pvalue:.3f}\n(shaded area)')
2799 >>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
2800 >>> _ = ax.annotate(annotation, (1.5, 0.22), (2.25, 0.3), arrowprops=props)
2801 >>> i = val >= res.statistic
2802 >>> ax.fill_between(val[i], y1=0, y2=pdf[i], color='C0')
2803 >>> plt.show()
2805 >>> res.pvalue
2806 0.71696121509966
2808 If the p-value is "small" - that is, if there is a low probability of
2809 sampling data from distributions with identical variances that produces
2810 such an extreme value of the statistic - this may be taken as evidence
2811 against the null hypothesis in favor of the alternative: the variances of
2812 the groups are not equal. Note that:
2814 - The inverse is not true; that is, the test is not used to provide
2815 evidence for the null hypothesis.
2816 - The threshold for values that will be considered "small" is a choice that
2817 should be made before the data is analyzed [6]_ with consideration of the
2818 risks of both false positives (incorrectly rejecting the null hypothesis)
2819 and false negatives (failure to reject a false null hypothesis).
2820 - Small p-values are not evidence for a *large* effect; rather, they can
2821 only provide evidence for a "significant" effect, meaning that they are
2822 unlikely to have occurred under the null hypothesis.
2824 Note that the chi-square distribution provides the null distribution
2825 when the observations are normally distributed. For small samples
2826 drawn from non-normal populations, it may be more appropriate to
2827 perform a
2828 permutation test: Under the null hypothesis that all three samples were
2829 drawn from the same population, each of the measurements is equally likely
2830 to have been observed in any of the three samples. Therefore, we can form
2831 a randomized null distribution by calculating the statistic under many
2832 randomly-generated partitionings of the observations into the three
2833 samples.
2835 >>> def statistic(*samples):
2836 ... return stats.bartlett(*samples).statistic
2837 >>> ref = stats.permutation_test(
2838 ... (small_dose, medium_dose, large_dose), statistic,
2839 ... permutation_type='independent', alternative='greater'
2840 ... )
2841 >>> fig, ax = plt.subplots(figsize=(8, 5))
2842 >>> plot(ax)
2843 >>> bins = np.linspace(0, 5, 25)
2844 >>> ax.hist(
2845 ... ref.null_distribution, bins=bins, density=True, facecolor="C1"
2846 ... )
2847 >>> ax.legend(['aymptotic approximation\n(many observations)',
2848 ... 'randomized null distribution'])
2849 >>> plot(ax)
2850 >>> plt.show()
2852 >>> ref.pvalue # randomized test p-value
2853 0.5387 # may vary
2855 Note that there is significant disagreement between the p-value calculated
2856 here and the asymptotic approximation returned by `bartlett` above.
2857 The statistical inferences that can be drawn rigorously from a permutation
2858 test are limited; nonetheless, they may be the preferred approach in many
2859 circumstances [7]_.
2861 Following is another generic example where the null hypothesis would be
2862 rejected.
2864 Test whether the lists `a`, `b` and `c` come from populations
2865 with equal variances.
2867 >>> a = [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
2868 >>> b = [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
2869 >>> c = [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
2870 >>> stat, p = stats.bartlett(a, b, c)
2871 >>> p
2872 1.1254782518834628e-05
2874 The very small p-value suggests that the populations do not have equal
2875 variances.
2877 This is not surprising, given that the sample variance of `b` is much
2878 larger than that of `a` and `c`:
2880 >>> [np.var(x, ddof=1) for x in [a, b, c]]
2881 [0.007054444444444413, 0.13073888888888888, 0.008890000000000002]
2883 """
2884 # Handle empty input and input that is not 1d
2885 for sample in samples:
2886 if np.asanyarray(sample).size == 0:
2887 return BartlettResult(np.nan, np.nan)
2888 if np.asanyarray(sample).ndim > 1:
2889 raise ValueError('Samples must be one-dimensional.')
2891 k = len(samples)
2892 if k < 2:
2893 raise ValueError("Must enter at least two input sample vectors.")
2894 Ni = np.empty(k)
2895 ssq = np.empty(k, 'd')
2896 for j in range(k):
2897 Ni[j] = len(samples[j])
2898 ssq[j] = np.var(samples[j], ddof=1)
2899 Ntot = np.sum(Ni, axis=0)
2900 spsq = np.sum((Ni - 1)*ssq, axis=0) / (1.0*(Ntot - k))
2901 numer = (Ntot*1.0 - k) * log(spsq) - np.sum((Ni - 1.0)*log(ssq), axis=0)
2902 denom = 1.0 + 1.0/(3*(k - 1)) * ((np.sum(1.0/(Ni - 1.0), axis=0)) -
2903 1.0/(Ntot - k))
2904 T = numer / denom
2905 pval = distributions.chi2.sf(T, k - 1) # 1 - cdf
2907 return BartlettResult(T, pval)
2910LeveneResult = namedtuple('LeveneResult', ('statistic', 'pvalue'))
2913def levene(*samples, center='median', proportiontocut=0.05):
2914 r"""Perform Levene test for equal variances.
2916 The Levene test tests the null hypothesis that all input samples
2917 are from populations with equal variances. Levene's test is an
2918 alternative to Bartlett's test `bartlett` in the case where
2919 there are significant deviations from normality.
2921 Parameters
2922 ----------
2923 sample1, sample2, ... : array_like
2924 The sample data, possibly with different lengths. Only one-dimensional
2925 samples are accepted.
2926 center : {'mean', 'median', 'trimmed'}, optional
2927 Which function of the data to use in the test. The default
2928 is 'median'.
2929 proportiontocut : float, optional
2930 When `center` is 'trimmed', this gives the proportion of data points
2931 to cut from each end. (See `scipy.stats.trim_mean`.)
2932 Default is 0.05.
2934 Returns
2935 -------
2936 statistic : float
2937 The test statistic.
2938 pvalue : float
2939 The p-value for the test.
2941 See Also
2942 --------
2943 fligner : A non-parametric test for the equality of k variances
2944 bartlett : A parametric test for equality of k variances in normal samples
2946 Notes
2947 -----
2948 Three variations of Levene's test are possible. The possibilities
2949 and their recommended usages are:
2951 * 'median' : Recommended for skewed (non-normal) distributions>
2952 * 'mean' : Recommended for symmetric, moderate-tailed distributions.
2953 * 'trimmed' : Recommended for heavy-tailed distributions.
2955 The test version using the mean was proposed in the original article
2956 of Levene ([2]_) while the median and trimmed mean have been studied by
2957 Brown and Forsythe ([3]_), sometimes also referred to as Brown-Forsythe
2958 test.
2960 References
2961 ----------
2962 .. [1] https://www.itl.nist.gov/div898/handbook/eda/section3/eda35a.htm
2963 .. [2] Levene, H. (1960). In Contributions to Probability and Statistics:
2964 Essays in Honor of Harold Hotelling, I. Olkin et al. eds.,
2965 Stanford University Press, pp. 278-292.
2966 .. [3] Brown, M. B. and Forsythe, A. B. (1974), Journal of the American
2967 Statistical Association, 69, 364-367
2968 .. [4] C.I. BLISS (1952), The Statistics of Bioassay: With Special
2969 Reference to the Vitamins, pp 499-503,
2970 :doi:`10.1016/C2013-0-12584-6`.
2971 .. [5] B. Phipson and G. K. Smyth. "Permutation P-values Should Never Be
2972 Zero: Calculating Exact P-values When Permutations Are Randomly
2973 Drawn." Statistical Applications in Genetics and Molecular Biology
2974 9.1 (2010).
2975 .. [6] Ludbrook, J., & Dudley, H. (1998). Why permutation tests are
2976 superior to t and F tests in biomedical research. The American
2977 Statistician, 52(2), 127-132.
2979 Examples
2980 --------
2981 In [4]_, the influence of vitamin C on the tooth growth of guinea pigs
2982 was investigated. In a control study, 60 subjects were divided into
2983 small dose, medium dose, and large dose groups that received
2984 daily doses of 0.5, 1.0 and 2.0 mg of vitamin C, respectively.
2985 After 42 days, the tooth growth was measured.
2987 The ``small_dose``, ``medium_dose``, and ``large_dose`` arrays below record
2988 tooth growth measurements of the three groups in microns.
2990 >>> import numpy as np
2991 >>> small_dose = np.array([
2992 ... 4.2, 11.5, 7.3, 5.8, 6.4, 10, 11.2, 11.2, 5.2, 7,
2993 ... 15.2, 21.5, 17.6, 9.7, 14.5, 10, 8.2, 9.4, 16.5, 9.7
2994 ... ])
2995 >>> medium_dose = np.array([
2996 ... 16.5, 16.5, 15.2, 17.3, 22.5, 17.3, 13.6, 14.5, 18.8, 15.5,
2997 ... 19.7, 23.3, 23.6, 26.4, 20, 25.2, 25.8, 21.2, 14.5, 27.3
2998 ... ])
2999 >>> large_dose = np.array([
3000 ... 23.6, 18.5, 33.9, 25.5, 26.4, 32.5, 26.7, 21.5, 23.3, 29.5,
3001 ... 25.5, 26.4, 22.4, 24.5, 24.8, 30.9, 26.4, 27.3, 29.4, 23
3002 ... ])
3004 The `levene` statistic is sensitive to differences in variances
3005 between the samples.
3007 >>> from scipy import stats
3008 >>> res = stats.levene(small_dose, medium_dose, large_dose)
3009 >>> res.statistic
3010 0.6457341109631506
3012 The value of the statistic tends to be high when there is a large
3013 difference in variances.
3015 We can test for inequality of variance among the groups by comparing the
3016 observed value of the statistic against the null distribution: the
3017 distribution of statistic values derived under the null hypothesis that
3018 the population variances of the three groups are equal.
3020 For this test, the null distribution follows the F distribution as shown
3021 below.
3023 >>> import matplotlib.pyplot as plt
3024 >>> k, n = 3, 60 # number of samples, total number of observations
3025 >>> dist = stats.f(dfn=k-1, dfd=n-k)
3026 >>> val = np.linspace(0, 5, 100)
3027 >>> pdf = dist.pdf(val)
3028 >>> fig, ax = plt.subplots(figsize=(8, 5))
3029 >>> def plot(ax): # we'll re-use this
3030 ... ax.plot(val, pdf, color='C0')
3031 ... ax.set_title("Levene Test Null Distribution")
3032 ... ax.set_xlabel("statistic")
3033 ... ax.set_ylabel("probability density")
3034 ... ax.set_xlim(0, 5)
3035 ... ax.set_ylim(0, 1)
3036 >>> plot(ax)
3037 >>> plt.show()
3039 The comparison is quantified by the p-value: the proportion of values in
3040 the null distribution greater than or equal to the observed value of the
3041 statistic.
3043 >>> fig, ax = plt.subplots(figsize=(8, 5))
3044 >>> plot(ax)
3045 >>> pvalue = dist.sf(res.statistic)
3046 >>> annotation = (f'p-value={pvalue:.3f}\n(shaded area)')
3047 >>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
3048 >>> _ = ax.annotate(annotation, (1.5, 0.22), (2.25, 0.3), arrowprops=props)
3049 >>> i = val >= res.statistic
3050 >>> ax.fill_between(val[i], y1=0, y2=pdf[i], color='C0')
3051 >>> plt.show()
3053 >>> res.pvalue
3054 0.5280694573759905
3056 If the p-value is "small" - that is, if there is a low probability of
3057 sampling data from distributions with identical variances that produces
3058 such an extreme value of the statistic - this may be taken as evidence
3059 against the null hypothesis in favor of the alternative: the variances of
3060 the groups are not equal. Note that:
3062 - The inverse is not true; that is, the test is not used to provide
3063 evidence for the null hypothesis.
3064 - The threshold for values that will be considered "small" is a choice that
3065 should be made before the data is analyzed [5]_ with consideration of the
3066 risks of both false positives (incorrectly rejecting the null hypothesis)
3067 and false negatives (failure to reject a false null hypothesis).
3068 - Small p-values are not evidence for a *large* effect; rather, they can
3069 only provide evidence for a "significant" effect, meaning that they are
3070 unlikely to have occurred under the null hypothesis.
3072 Note that the F distribution provides an asymptotic approximation of the
3073 null distribution.
3074 For small samples, it may be more appropriate to perform a permutation
3075 test: Under the null hypothesis that all three samples were drawn from
3076 the same population, each of the measurements is equally likely to have
3077 been observed in any of the three samples. Therefore, we can form a
3078 randomized null distribution by calculating the statistic under many
3079 randomly-generated partitionings of the observations into the three
3080 samples.
3082 >>> def statistic(*samples):
3083 ... return stats.levene(*samples).statistic
3084 >>> ref = stats.permutation_test(
3085 ... (small_dose, medium_dose, large_dose), statistic,
3086 ... permutation_type='independent', alternative='greater'
3087 ... )
3088 >>> fig, ax = plt.subplots(figsize=(8, 5))
3089 >>> plot(ax)
3090 >>> bins = np.linspace(0, 5, 25)
3091 >>> ax.hist(
3092 ... ref.null_distribution, bins=bins, density=True, facecolor="C1"
3093 ... )
3094 >>> ax.legend(['aymptotic approximation\n(many observations)',
3095 ... 'randomized null distribution'])
3096 >>> plot(ax)
3097 >>> plt.show()
3099 >>> ref.pvalue # randomized test p-value
3100 0.4559 # may vary
3102 Note that there is significant disagreement between the p-value calculated
3103 here and the asymptotic approximation returned by `levene` above.
3104 The statistical inferences that can be drawn rigorously from a permutation
3105 test are limited; nonetheless, they may be the preferred approach in many
3106 circumstances [6]_.
3108 Following is another generic example where the null hypothesis would be
3109 rejected.
3111 Test whether the lists `a`, `b` and `c` come from populations
3112 with equal variances.
3114 >>> a = [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
3115 >>> b = [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
3116 >>> c = [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
3117 >>> stat, p = stats.levene(a, b, c)
3118 >>> p
3119 0.002431505967249681
3121 The small p-value suggests that the populations do not have equal
3122 variances.
3124 This is not surprising, given that the sample variance of `b` is much
3125 larger than that of `a` and `c`:
3127 >>> [np.var(x, ddof=1) for x in [a, b, c]]
3128 [0.007054444444444413, 0.13073888888888888, 0.008890000000000002]
3130 """
3131 if center not in ['mean', 'median', 'trimmed']:
3132 raise ValueError("center must be 'mean', 'median' or 'trimmed'.")
3134 k = len(samples)
3135 if k < 2:
3136 raise ValueError("Must enter at least two input sample vectors.")
3137 # check for 1d input
3138 for j in range(k):
3139 if np.asanyarray(samples[j]).ndim > 1:
3140 raise ValueError('Samples must be one-dimensional.')
3142 Ni = np.empty(k)
3143 Yci = np.empty(k, 'd')
3145 if center == 'median':
3147 def func(x):
3148 return np.median(x, axis=0)
3150 elif center == 'mean':
3152 def func(x):
3153 return np.mean(x, axis=0)
3155 else: # center == 'trimmed'
3156 samples = tuple(_stats_py.trimboth(np.sort(sample), proportiontocut)
3157 for sample in samples)
3159 def func(x):
3160 return np.mean(x, axis=0)
3162 for j in range(k):
3163 Ni[j] = len(samples[j])
3164 Yci[j] = func(samples[j])
3165 Ntot = np.sum(Ni, axis=0)
3167 # compute Zij's
3168 Zij = [None] * k
3169 for i in range(k):
3170 Zij[i] = abs(asarray(samples[i]) - Yci[i])
3172 # compute Zbari
3173 Zbari = np.empty(k, 'd')
3174 Zbar = 0.0
3175 for i in range(k):
3176 Zbari[i] = np.mean(Zij[i], axis=0)
3177 Zbar += Zbari[i] * Ni[i]
3179 Zbar /= Ntot
3180 numer = (Ntot - k) * np.sum(Ni * (Zbari - Zbar)**2, axis=0)
3182 # compute denom_variance
3183 dvar = 0.0
3184 for i in range(k):
3185 dvar += np.sum((Zij[i] - Zbari[i])**2, axis=0)
3187 denom = (k - 1.0) * dvar
3189 W = numer / denom
3190 pval = distributions.f.sf(W, k-1, Ntot-k) # 1 - cdf
3191 return LeveneResult(W, pval)
3194@_deprecated("'binom_test' is deprecated in favour of"
3195 " 'binomtest' from version 1.7.0 and will"
3196 " be removed in Scipy 1.12.0.")
3197def binom_test(x, n=None, p=0.5, alternative='two-sided'):
3198 """Perform a test that the probability of success is p.
3200 This is an exact, two-sided test of the null hypothesis
3201 that the probability of success in a Bernoulli experiment
3202 is `p`.
3204 .. deprecated:: 1.10.0
3205 `binom_test` is deprecated in favour of `binomtest` and will
3206 be removed in Scipy 1.12.0.
3208 Parameters
3209 ----------
3210 x : int or array_like
3211 The number of successes, or if x has length 2, it is the
3212 number of successes and the number of failures.
3213 n : int
3214 The number of trials. This is ignored if x gives both the
3215 number of successes and failures.
3216 p : float, optional
3217 The hypothesized probability of success. ``0 <= p <= 1``. The
3218 default value is ``p = 0.5``.
3219 alternative : {'two-sided', 'greater', 'less'}, optional
3220 Indicates the alternative hypothesis. The default value is
3221 'two-sided'.
3223 Returns
3224 -------
3225 p-value : float
3226 The p-value of the hypothesis test.
3228 References
3229 ----------
3230 .. [1] https://en.wikipedia.org/wiki/Binomial_test
3232 Examples
3233 --------
3234 >>> from scipy import stats
3236 A car manufacturer claims that no more than 10% of their cars are unsafe.
3237 15 cars are inspected for safety, 3 were found to be unsafe. Test the
3238 manufacturer's claim:
3240 >>> stats.binom_test(3, n=15, p=0.1, alternative='greater')
3241 0.18406106910639114
3243 The null hypothesis cannot be rejected at the 5% level of significance
3244 because the returned p-value is greater than the critical value of 5%.
3246 """
3247 x = atleast_1d(x).astype(np.int_)
3248 if len(x) == 2:
3249 n = x[1] + x[0]
3250 x = x[0]
3251 elif len(x) == 1:
3252 x = x[0]
3253 if n is None or n < x:
3254 raise ValueError("n must be >= x")
3255 n = np.int_(n)
3256 else:
3257 raise ValueError("Incorrect length for x.")
3259 if (p > 1.0) or (p < 0.0):
3260 raise ValueError("p must be in range [0,1]")
3262 if alternative not in ('two-sided', 'less', 'greater'):
3263 raise ValueError("alternative not recognized\n"
3264 "should be 'two-sided', 'less' or 'greater'")
3266 if alternative == 'less':
3267 pval = distributions.binom.cdf(x, n, p)
3268 return pval
3270 if alternative == 'greater':
3271 pval = distributions.binom.sf(x-1, n, p)
3272 return pval
3274 # if alternative was neither 'less' nor 'greater', then it's 'two-sided'
3275 d = distributions.binom.pmf(x, n, p)
3276 rerr = 1 + 1e-7
3277 if x == p * n:
3278 # special case as shortcut, would also be handled by `else` below
3279 pval = 1.
3280 elif x < p * n:
3281 i = np.arange(np.ceil(p * n), n+1)
3282 y = np.sum(distributions.binom.pmf(i, n, p) <= d*rerr, axis=0)
3283 pval = (distributions.binom.cdf(x, n, p) +
3284 distributions.binom.sf(n - y, n, p))
3285 else:
3286 i = np.arange(np.floor(p*n) + 1)
3287 y = np.sum(distributions.binom.pmf(i, n, p) <= d*rerr, axis=0)
3288 pval = (distributions.binom.cdf(y-1, n, p) +
3289 distributions.binom.sf(x-1, n, p))
3291 return min(1.0, pval)
3294def _apply_func(x, g, func):
3295 # g is list of indices into x
3296 # separating x into different groups
3297 # func should be applied over the groups
3298 g = unique(r_[0, g, len(x)])
3299 output = [func(x[g[k]:g[k+1]]) for k in range(len(g) - 1)]
3301 return asarray(output)
3304FlignerResult = namedtuple('FlignerResult', ('statistic', 'pvalue'))
3307def fligner(*samples, center='median', proportiontocut=0.05):
3308 r"""Perform Fligner-Killeen test for equality of variance.
3310 Fligner's test tests the null hypothesis that all input samples
3311 are from populations with equal variances. Fligner-Killeen's test is
3312 distribution free when populations are identical [2]_.
3314 Parameters
3315 ----------
3316 sample1, sample2, ... : array_like
3317 Arrays of sample data. Need not be the same length.
3318 center : {'mean', 'median', 'trimmed'}, optional
3319 Keyword argument controlling which function of the data is used in
3320 computing the test statistic. The default is 'median'.
3321 proportiontocut : float, optional
3322 When `center` is 'trimmed', this gives the proportion of data points
3323 to cut from each end. (See `scipy.stats.trim_mean`.)
3324 Default is 0.05.
3326 Returns
3327 -------
3328 statistic : float
3329 The test statistic.
3330 pvalue : float
3331 The p-value for the hypothesis test.
3333 See Also
3334 --------
3335 bartlett : A parametric test for equality of k variances in normal samples
3336 levene : A robust parametric test for equality of k variances
3338 Notes
3339 -----
3340 As with Levene's test there are three variants of Fligner's test that
3341 differ by the measure of central tendency used in the test. See `levene`
3342 for more information.
3344 Conover et al. (1981) examine many of the existing parametric and
3345 nonparametric tests by extensive simulations and they conclude that the
3346 tests proposed by Fligner and Killeen (1976) and Levene (1960) appear to be
3347 superior in terms of robustness of departures from normality and power
3348 [3]_.
3350 References
3351 ----------
3352 .. [1] Park, C. and Lindsay, B. G. (1999). Robust Scale Estimation and
3353 Hypothesis Testing based on Quadratic Inference Function. Technical
3354 Report #99-03, Center for Likelihood Studies, Pennsylvania State
3355 University.
3356 https://cecas.clemson.edu/~cspark/cv/paper/qif/draftqif2.pdf
3357 .. [2] Fligner, M.A. and Killeen, T.J. (1976). Distribution-free two-sample
3358 tests for scale. 'Journal of the American Statistical Association.'
3359 71(353), 210-213.
3360 .. [3] Park, C. and Lindsay, B. G. (1999). Robust Scale Estimation and
3361 Hypothesis Testing based on Quadratic Inference Function. Technical
3362 Report #99-03, Center for Likelihood Studies, Pennsylvania State
3363 University.
3364 .. [4] Conover, W. J., Johnson, M. E. and Johnson M. M. (1981). A
3365 comparative study of tests for homogeneity of variances, with
3366 applications to the outer continental shelf biding data.
3367 Technometrics, 23(4), 351-361.
3368 .. [5] C.I. BLISS (1952), The Statistics of Bioassay: With Special
3369 Reference to the Vitamins, pp 499-503,
3370 :doi:`10.1016/C2013-0-12584-6`.
3371 .. [6] B. Phipson and G. K. Smyth. "Permutation P-values Should Never Be
3372 Zero: Calculating Exact P-values When Permutations Are Randomly
3373 Drawn." Statistical Applications in Genetics and Molecular Biology
3374 9.1 (2010).
3375 .. [7] Ludbrook, J., & Dudley, H. (1998). Why permutation tests are
3376 superior to t and F tests in biomedical research. The American
3377 Statistician, 52(2), 127-132.
3379 Examples
3380 --------
3381 In [5]_, the influence of vitamin C on the tooth growth of guinea pigs
3382 was investigated. In a control study, 60 subjects were divided into
3383 small dose, medium dose, and large dose groups that received
3384 daily doses of 0.5, 1.0 and 2.0 mg of vitamin C, respectively.
3385 After 42 days, the tooth growth was measured.
3387 The ``small_dose``, ``medium_dose``, and ``large_dose`` arrays below record
3388 tooth growth measurements of the three groups in microns.
3390 >>> import numpy as np
3391 >>> small_dose = np.array([
3392 ... 4.2, 11.5, 7.3, 5.8, 6.4, 10, 11.2, 11.2, 5.2, 7,
3393 ... 15.2, 21.5, 17.6, 9.7, 14.5, 10, 8.2, 9.4, 16.5, 9.7
3394 ... ])
3395 >>> medium_dose = np.array([
3396 ... 16.5, 16.5, 15.2, 17.3, 22.5, 17.3, 13.6, 14.5, 18.8, 15.5,
3397 ... 19.7, 23.3, 23.6, 26.4, 20, 25.2, 25.8, 21.2, 14.5, 27.3
3398 ... ])
3399 >>> large_dose = np.array([
3400 ... 23.6, 18.5, 33.9, 25.5, 26.4, 32.5, 26.7, 21.5, 23.3, 29.5,
3401 ... 25.5, 26.4, 22.4, 24.5, 24.8, 30.9, 26.4, 27.3, 29.4, 23
3402 ... ])
3404 The `fligner` statistic is sensitive to differences in variances
3405 between the samples.
3407 >>> from scipy import stats
3408 >>> res = stats.fligner(small_dose, medium_dose, large_dose)
3409 >>> res.statistic
3410 1.3878943408857916
3412 The value of the statistic tends to be high when there is a large
3413 difference in variances.
3415 We can test for inequality of variance among the groups by comparing the
3416 observed value of the statistic against the null distribution: the
3417 distribution of statistic values derived under the null hypothesis that
3418 the population variances of the three groups are equal.
3420 For this test, the null distribution follows the chi-square distribution
3421 as shown below.
3423 >>> import matplotlib.pyplot as plt
3424 >>> k = 3 # number of samples
3425 >>> dist = stats.chi2(df=k-1)
3426 >>> val = np.linspace(0, 8, 100)
3427 >>> pdf = dist.pdf(val)
3428 >>> fig, ax = plt.subplots(figsize=(8, 5))
3429 >>> def plot(ax): # we'll re-use this
3430 ... ax.plot(val, pdf, color='C0')
3431 ... ax.set_title("Fligner Test Null Distribution")
3432 ... ax.set_xlabel("statistic")
3433 ... ax.set_ylabel("probability density")
3434 ... ax.set_xlim(0, 8)
3435 ... ax.set_ylim(0, 0.5)
3436 >>> plot(ax)
3437 >>> plt.show()
3439 The comparison is quantified by the p-value: the proportion of values in
3440 the null distribution greater than or equal to the observed value of the
3441 statistic.
3443 >>> fig, ax = plt.subplots(figsize=(8, 5))
3444 >>> plot(ax)
3445 >>> pvalue = dist.sf(res.statistic)
3446 >>> annotation = (f'p-value={pvalue:.4f}\n(shaded area)')
3447 >>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
3448 >>> _ = ax.annotate(annotation, (1.5, 0.22), (2.25, 0.3), arrowprops=props)
3449 >>> i = val >= res.statistic
3450 >>> ax.fill_between(val[i], y1=0, y2=pdf[i], color='C0')
3451 >>> plt.show()
3453 >>> res.pvalue
3454 0.49960016501182125
3456 If the p-value is "small" - that is, if there is a low probability of
3457 sampling data from distributions with identical variances that produces
3458 such an extreme value of the statistic - this may be taken as evidence
3459 against the null hypothesis in favor of the alternative: the variances of
3460 the groups are not equal. Note that:
3462 - The inverse is not true; that is, the test is not used to provide
3463 evidence for the null hypothesis.
3464 - The threshold for values that will be considered "small" is a choice that
3465 should be made before the data is analyzed [6]_ with consideration of the
3466 risks of both false positives (incorrectly rejecting the null hypothesis)
3467 and false negatives (failure to reject a false null hypothesis).
3468 - Small p-values are not evidence for a *large* effect; rather, they can
3469 only provide evidence for a "significant" effect, meaning that they are
3470 unlikely to have occurred under the null hypothesis.
3472 Note that the chi-square distribution provides an asymptotic approximation
3473 of the null distribution.
3474 For small samples, it may be more appropriate to perform a
3475 permutation test: Under the null hypothesis that all three samples were
3476 drawn from the same population, each of the measurements is equally likely
3477 to have been observed in any of the three samples. Therefore, we can form
3478 a randomized null distribution by calculating the statistic under many
3479 randomly-generated partitionings of the observations into the three
3480 samples.
3482 >>> def statistic(*samples):
3483 ... return stats.fligner(*samples).statistic
3484 >>> ref = stats.permutation_test(
3485 ... (small_dose, medium_dose, large_dose), statistic,
3486 ... permutation_type='independent', alternative='greater'
3487 ... )
3488 >>> fig, ax = plt.subplots(figsize=(8, 5))
3489 >>> plot(ax)
3490 >>> bins = np.linspace(0, 8, 25)
3491 >>> ax.hist(
3492 ... ref.null_distribution, bins=bins, density=True, facecolor="C1"
3493 ... )
3494 >>> ax.legend(['aymptotic approximation\n(many observations)',
3495 ... 'randomized null distribution'])
3496 >>> plot(ax)
3497 >>> plt.show()
3499 >>> ref.pvalue # randomized test p-value
3500 0.4332 # may vary
3502 Note that there is significant disagreement between the p-value calculated
3503 here and the asymptotic approximation returned by `fligner` above.
3504 The statistical inferences that can be drawn rigorously from a permutation
3505 test are limited; nonetheless, they may be the preferred approach in many
3506 circumstances [7]_.
3508 Following is another generic example where the null hypothesis would be
3509 rejected.
3511 Test whether the lists `a`, `b` and `c` come from populations
3512 with equal variances.
3514 >>> a = [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
3515 >>> b = [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
3516 >>> c = [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
3517 >>> stat, p = stats.fligner(a, b, c)
3518 >>> p
3519 0.00450826080004775
3521 The small p-value suggests that the populations do not have equal
3522 variances.
3524 This is not surprising, given that the sample variance of `b` is much
3525 larger than that of `a` and `c`:
3527 >>> [np.var(x, ddof=1) for x in [a, b, c]]
3528 [0.007054444444444413, 0.13073888888888888, 0.008890000000000002]
3530 """
3531 if center not in ['mean', 'median', 'trimmed']:
3532 raise ValueError("center must be 'mean', 'median' or 'trimmed'.")
3534 # Handle empty input
3535 for sample in samples:
3536 if np.asanyarray(sample).size == 0:
3537 return FlignerResult(np.nan, np.nan)
3539 k = len(samples)
3540 if k < 2:
3541 raise ValueError("Must enter at least two input sample vectors.")
3543 if center == 'median':
3545 def func(x):
3546 return np.median(x, axis=0)
3548 elif center == 'mean':
3550 def func(x):
3551 return np.mean(x, axis=0)
3553 else: # center == 'trimmed'
3554 samples = tuple(_stats_py.trimboth(sample, proportiontocut)
3555 for sample in samples)
3557 def func(x):
3558 return np.mean(x, axis=0)
3560 Ni = asarray([len(samples[j]) for j in range(k)])
3561 Yci = asarray([func(samples[j]) for j in range(k)])
3562 Ntot = np.sum(Ni, axis=0)
3563 # compute Zij's
3564 Zij = [abs(asarray(samples[i]) - Yci[i]) for i in range(k)]
3565 allZij = []
3566 g = [0]
3567 for i in range(k):
3568 allZij.extend(list(Zij[i]))
3569 g.append(len(allZij))
3571 ranks = _stats_py.rankdata(allZij)
3572 sample = distributions.norm.ppf(ranks / (2*(Ntot + 1.0)) + 0.5)
3574 # compute Aibar
3575 Aibar = _apply_func(sample, g, np.sum) / Ni
3576 anbar = np.mean(sample, axis=0)
3577 varsq = np.var(sample, axis=0, ddof=1)
3578 Xsq = np.sum(Ni * (asarray(Aibar) - anbar)**2.0, axis=0) / varsq
3579 pval = distributions.chi2.sf(Xsq, k - 1) # 1 - cdf
3580 return FlignerResult(Xsq, pval)
3583@_axis_nan_policy_factory(lambda x1: (x1,), n_samples=4, n_outputs=1)
3584def _mood_inner_lc(xy, x, diffs, sorted_xy, n, m, N) -> float:
3585 # Obtain the unique values and their frequencies from the pooled samples.
3586 # "a_j, + b_j, = t_j, for j = 1, ... k" where `k` is the number of unique
3587 # classes, and "[t]he number of values associated with the x's and y's in
3588 # the jth class will be denoted by a_j, and b_j respectively."
3589 # (Mielke, 312)
3590 # Reuse previously computed sorted array and `diff` arrays to obtain the
3591 # unique values and counts. Prepend `diffs` with a non-zero to indicate
3592 # that the first element should be marked as not matching what preceded it.
3593 diffs_prep = np.concatenate(([1], diffs))
3594 # Unique elements are where the was a difference between elements in the
3595 # sorted array
3596 uniques = sorted_xy[diffs_prep != 0]
3597 # The count of each element is the bin size for each set of consecutive
3598 # differences where the difference is zero. Replace nonzero differences
3599 # with 1 and then use the cumulative sum to count the indices.
3600 t = np.bincount(np.cumsum(np.asarray(diffs_prep != 0, dtype=int)))[1:]
3601 k = len(uniques)
3602 js = np.arange(1, k + 1, dtype=int)
3603 # the `b` array mentioned in the paper is not used, outside of the
3604 # calculation of `t`, so we do not need to calculate it separately. Here
3605 # we calculate `a`. In plain language, `a[j]` is the number of values in
3606 # `x` that equal `uniques[j]`.
3607 sorted_xyx = np.sort(np.concatenate((xy, x)))
3608 diffs = np.diff(sorted_xyx)
3609 diffs_prep = np.concatenate(([1], diffs))
3610 diff_is_zero = np.asarray(diffs_prep != 0, dtype=int)
3611 xyx_counts = np.bincount(np.cumsum(diff_is_zero))[1:]
3612 a = xyx_counts - t
3613 # "Define .. a_0 = b_0 = t_0 = S_0 = 0" (Mielke 312) so we shift `a`
3614 # and `t` arrays over 1 to allow a first element of 0 to accommodate this
3615 # indexing.
3616 t = np.concatenate(([0], t))
3617 a = np.concatenate(([0], a))
3618 # S is built from `t`, so it does not need a preceding zero added on.
3619 S = np.cumsum(t)
3620 # define a copy of `S` with a prepending zero for later use to avoid
3621 # the need for indexing.
3622 S_i_m1 = np.concatenate(([0], S[:-1]))
3624 # Psi, as defined by the 6th unnumbered equation on page 313 (Mielke).
3625 # Note that in the paper there is an error where the denominator `2` is
3626 # squared when it should be the entire equation.
3627 def psi(indicator):
3628 return (indicator - (N + 1)/2)**2
3630 # define summation range for use in calculation of phi, as seen in sum
3631 # in the unnumbered equation on the bottom of page 312 (Mielke).
3632 s_lower = S[js - 1] + 1
3633 s_upper = S[js] + 1
3634 phi_J = [np.arange(s_lower[idx], s_upper[idx]) for idx in range(k)]
3636 # for every range in the above array, determine the sum of psi(I) for
3637 # every element in the range. Divide all the sums by `t`. Following the
3638 # last unnumbered equation on page 312.
3639 phis = [np.sum(psi(I_j)) for I_j in phi_J] / t[js]
3641 # `T` is equal to a[j] * phi[j], per the first unnumbered equation on
3642 # page 312. `phis` is already in the order based on `js`, so we index
3643 # into `a` with `js` as well.
3644 T = sum(phis * a[js])
3646 # The approximate statistic
3647 E_0_T = n * (N * N - 1) / 12
3649 varM = (m * n * (N + 1.0) * (N ** 2 - 4) / 180 -
3650 m * n / (180 * N * (N - 1)) * np.sum(
3651 t * (t**2 - 1) * (t**2 - 4 + (15 * (N - S - S_i_m1) ** 2))
3652 ))
3654 return ((T - E_0_T) / np.sqrt(varM),)
3657def mood(x, y, axis=0, alternative="two-sided"):
3658 """Perform Mood's test for equal scale parameters.
3660 Mood's two-sample test for scale parameters is a non-parametric
3661 test for the null hypothesis that two samples are drawn from the
3662 same distribution with the same scale parameter.
3664 Parameters
3665 ----------
3666 x, y : array_like
3667 Arrays of sample data.
3668 axis : int, optional
3669 The axis along which the samples are tested. `x` and `y` can be of
3670 different length along `axis`.
3671 If `axis` is None, `x` and `y` are flattened and the test is done on
3672 all values in the flattened arrays.
3673 alternative : {'two-sided', 'less', 'greater'}, optional
3674 Defines the alternative hypothesis. Default is 'two-sided'.
3675 The following options are available:
3677 * 'two-sided': the scales of the distributions underlying `x` and `y`
3678 are different.
3679 * 'less': the scale of the distribution underlying `x` is less than
3680 the scale of the distribution underlying `y`.
3681 * 'greater': the scale of the distribution underlying `x` is greater
3682 than the scale of the distribution underlying `y`.
3684 .. versionadded:: 1.7.0
3686 Returns
3687 -------
3688 res : SignificanceResult
3689 An object containing attributes:
3691 statistic : scalar or ndarray
3692 The z-score for the hypothesis test. For 1-D inputs a scalar is
3693 returned.
3694 pvalue : scalar ndarray
3695 The p-value for the hypothesis test.
3697 See Also
3698 --------
3699 fligner : A non-parametric test for the equality of k variances
3700 ansari : A non-parametric test for the equality of 2 variances
3701 bartlett : A parametric test for equality of k variances in normal samples
3702 levene : A parametric test for equality of k variances
3704 Notes
3705 -----
3706 The data are assumed to be drawn from probability distributions ``f(x)``
3707 and ``f(x/s) / s`` respectively, for some probability density function f.
3708 The null hypothesis is that ``s == 1``.
3710 For multi-dimensional arrays, if the inputs are of shapes
3711 ``(n0, n1, n2, n3)`` and ``(n0, m1, n2, n3)``, then if ``axis=1``, the
3712 resulting z and p values will have shape ``(n0, n2, n3)``. Note that
3713 ``n1`` and ``m1`` don't have to be equal, but the other dimensions do.
3715 References
3716 ----------
3717 [1] Mielke, Paul W. "Note on Some Squared Rank Tests with Existing Ties."
3718 Technometrics, vol. 9, no. 2, 1967, pp. 312-14. JSTOR,
3719 https://doi.org/10.2307/1266427. Accessed 18 May 2022.
3721 Examples
3722 --------
3723 >>> import numpy as np
3724 >>> from scipy import stats
3725 >>> rng = np.random.default_rng()
3726 >>> x2 = rng.standard_normal((2, 45, 6, 7))
3727 >>> x1 = rng.standard_normal((2, 30, 6, 7))
3728 >>> res = stats.mood(x1, x2, axis=1)
3729 >>> res.pvalue.shape
3730 (2, 6, 7)
3732 Find the number of points where the difference in scale is not significant:
3734 >>> (res.pvalue > 0.1).sum()
3735 78
3737 Perform the test with different scales:
3739 >>> x1 = rng.standard_normal((2, 30))
3740 >>> x2 = rng.standard_normal((2, 35)) * 10.0
3741 >>> stats.mood(x1, x2, axis=1)
3742 SignificanceResult(statistic=array([-5.76174136, -6.12650783]),
3743 pvalue=array([8.32505043e-09, 8.98287869e-10]))
3745 """
3746 x = np.asarray(x, dtype=float)
3747 y = np.asarray(y, dtype=float)
3749 if axis is None:
3750 x = x.flatten()
3751 y = y.flatten()
3752 axis = 0
3754 if axis < 0:
3755 axis = x.ndim + axis
3757 # Determine shape of the result arrays
3758 res_shape = tuple([x.shape[ax] for ax in range(len(x.shape)) if ax != axis])
3759 if not (res_shape == tuple([y.shape[ax] for ax in range(len(y.shape)) if
3760 ax != axis])):
3761 raise ValueError("Dimensions of x and y on all axes except `axis` "
3762 "should match")
3764 n = x.shape[axis]
3765 m = y.shape[axis]
3766 N = m + n
3767 if N < 3:
3768 raise ValueError("Not enough observations.")
3770 xy = np.concatenate((x, y), axis=axis)
3771 # determine if any of the samples contain ties
3772 sorted_xy = np.sort(xy, axis=axis)
3773 diffs = np.diff(sorted_xy, axis=axis)
3774 if 0 in diffs:
3775 z = np.asarray(_mood_inner_lc(xy, x, diffs, sorted_xy, n, m, N,
3776 axis=axis))
3777 else:
3778 if axis != 0:
3779 xy = np.moveaxis(xy, axis, 0)
3781 xy = xy.reshape(xy.shape[0], -1)
3782 # Generalized to the n-dimensional case by adding the axis argument,
3783 # and using for loops, since rankdata is not vectorized. For improving
3784 # performance consider vectorizing rankdata function.
3785 all_ranks = np.empty_like(xy)
3786 for j in range(xy.shape[1]):
3787 all_ranks[:, j] = _stats_py.rankdata(xy[:, j])
3789 Ri = all_ranks[:n]
3790 M = np.sum((Ri - (N + 1.0) / 2) ** 2, axis=0)
3791 # Approx stat.
3792 mnM = n * (N * N - 1.0) / 12
3793 varM = m * n * (N + 1.0) * (N + 2) * (N - 2) / 180
3794 z = (M - mnM) / sqrt(varM)
3795 z, pval = _normtest_finish(z, alternative)
3797 if res_shape == ():
3798 # Return scalars, not 0-D arrays
3799 z = z[0]
3800 pval = pval[0]
3801 else:
3802 z.shape = res_shape
3803 pval.shape = res_shape
3804 return SignificanceResult(z, pval)
3807WilcoxonResult = _make_tuple_bunch('WilcoxonResult', ['statistic', 'pvalue'])
3810def wilcoxon_result_unpacker(res):
3811 if hasattr(res, 'zstatistic'):
3812 return res.statistic, res.pvalue, res.zstatistic
3813 else:
3814 return res.statistic, res.pvalue
3817def wilcoxon_result_object(statistic, pvalue, zstatistic=None):
3818 res = WilcoxonResult(statistic, pvalue)
3819 if zstatistic is not None:
3820 res.zstatistic = zstatistic
3821 return res
3824def wilcoxon_outputs(kwds):
3825 method = kwds.get('method', 'auto')
3826 if method == 'approx':
3827 return 3
3828 return 2
3831@_rename_parameter("mode", "method")
3832@_axis_nan_policy_factory(
3833 wilcoxon_result_object, paired=True,
3834 n_samples=lambda kwds: 2 if kwds.get('y', None) is not None else 1,
3835 result_to_tuple=wilcoxon_result_unpacker, n_outputs=wilcoxon_outputs,
3836)
3837def wilcoxon(x, y=None, zero_method="wilcox", correction=False,
3838 alternative="two-sided", method='auto'):
3839 """Calculate the Wilcoxon signed-rank test.
3841 The Wilcoxon signed-rank test tests the null hypothesis that two
3842 related paired samples come from the same distribution. In particular,
3843 it tests whether the distribution of the differences ``x - y`` is symmetric
3844 about zero. It is a non-parametric version of the paired T-test.
3846 Parameters
3847 ----------
3848 x : array_like
3849 Either the first set of measurements (in which case ``y`` is the second
3850 set of measurements), or the differences between two sets of
3851 measurements (in which case ``y`` is not to be specified.) Must be
3852 one-dimensional.
3853 y : array_like, optional
3854 Either the second set of measurements (if ``x`` is the first set of
3855 measurements), or not specified (if ``x`` is the differences between
3856 two sets of measurements.) Must be one-dimensional.
3858 .. warning::
3859 When `y` is provided, `wilcoxon` calculates the test statistic
3860 based on the ranks of the absolute values of ``d = x - y``.
3861 Roundoff error in the subtraction can result in elements of ``d``
3862 being assigned different ranks even when they would be tied with
3863 exact arithmetic. Rather than passing `x` and `y` separately,
3864 consider computing the difference ``x - y``, rounding as needed to
3865 ensure that only truly unique elements are numerically distinct,
3866 and passing the result as `x`, leaving `y` at the default (None).
3868 zero_method : {"wilcox", "pratt", "zsplit"}, optional
3869 There are different conventions for handling pairs of observations
3870 with equal values ("zero-differences", or "zeros").
3872 * "wilcox": Discards all zero-differences (default); see [4]_.
3873 * "pratt": Includes zero-differences in the ranking process,
3874 but drops the ranks of the zeros (more conservative); see [3]_.
3875 In this case, the normal approximation is adjusted as in [5]_.
3876 * "zsplit": Includes zero-differences in the ranking process and
3877 splits the zero rank between positive and negative ones.
3879 correction : bool, optional
3880 If True, apply continuity correction by adjusting the Wilcoxon rank
3881 statistic by 0.5 towards the mean value when computing the
3882 z-statistic if a normal approximation is used. Default is False.
3883 alternative : {"two-sided", "greater", "less"}, optional
3884 Defines the alternative hypothesis. Default is 'two-sided'.
3885 In the following, let ``d`` represent the difference between the paired
3886 samples: ``d = x - y`` if both ``x`` and ``y`` are provided, or
3887 ``d = x`` otherwise.
3889 * 'two-sided': the distribution underlying ``d`` is not symmetric
3890 about zero.
3891 * 'less': the distribution underlying ``d`` is stochastically less
3892 than a distribution symmetric about zero.
3893 * 'greater': the distribution underlying ``d`` is stochastically
3894 greater than a distribution symmetric about zero.
3896 method : {"auto", "exact", "approx"}, optional
3897 Method to calculate the p-value, see Notes. Default is "auto".
3899 Returns
3900 -------
3901 An object with the following attributes.
3903 statistic : array_like
3904 If `alternative` is "two-sided", the sum of the ranks of the
3905 differences above or below zero, whichever is smaller.
3906 Otherwise the sum of the ranks of the differences above zero.
3907 pvalue : array_like
3908 The p-value for the test depending on `alternative` and `method`.
3909 zstatistic : array_like
3910 When ``method = 'approx'``, this is the normalized z-statistic::
3912 z = (T - mn - d) / se
3914 where ``T`` is `statistic` as defined above, ``mn`` is the mean of the
3915 distribution under the null hypothesis, ``d`` is a continuity
3916 correction, and ``se`` is the standard error.
3917 When ``method != 'approx'``, this attribute is not available.
3919 See Also
3920 --------
3921 kruskal, mannwhitneyu
3923 Notes
3924 -----
3925 In the following, let ``d`` represent the difference between the paired
3926 samples: ``d = x - y`` if both ``x`` and ``y`` are provided, or ``d = x``
3927 otherwise. Assume that all elements of ``d`` are independent and
3928 identically distributed observations, and all are distinct and nonzero.
3930 - When ``len(d)`` is sufficiently large, the null distribution of the
3931 normalized test statistic (`zstatistic` above) is approximately normal,
3932 and ``method = 'approx'`` can be used to compute the p-value.
3934 - When ``len(d)`` is small, the normal approximation may not be accurate,
3935 and ``method='exact'`` is preferred (at the cost of additional
3936 execution time).
3938 - The default, ``method='auto'``, selects between the two: when
3939 ``len(d) <= 50``, the exact method is used; otherwise, the approximate
3940 method is used.
3942 The presence of "ties" (i.e. not all elements of ``d`` are unique) and
3943 "zeros" (i.e. elements of ``d`` are zero) changes the null distribution
3944 of the test statistic, and ``method='exact'`` no longer calculates
3945 the exact p-value. If ``method='approx'``, the z-statistic is adjusted
3946 for more accurate comparison against the standard normal, but still,
3947 for finite sample sizes, the standard normal is only an approximation of
3948 the true null distribution of the z-statistic. There is no clear
3949 consensus among references on which method most accurately approximates
3950 the p-value for small samples in the presence of zeros and/or ties. In any
3951 case, this is the behavior of `wilcoxon` when ``method='auto':
3952 ``method='exact'`` is used when ``len(d) <= 50`` *and there are no zeros*;
3953 otherwise, ``method='approx'`` is used.
3955 References
3956 ----------
3957 .. [1] https://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test
3958 .. [2] Conover, W.J., Practical Nonparametric Statistics, 1971.
3959 .. [3] Pratt, J.W., Remarks on Zeros and Ties in the Wilcoxon Signed
3960 Rank Procedures, Journal of the American Statistical Association,
3961 Vol. 54, 1959, pp. 655-667. :doi:`10.1080/01621459.1959.10501526`
3962 .. [4] Wilcoxon, F., Individual Comparisons by Ranking Methods,
3963 Biometrics Bulletin, Vol. 1, 1945, pp. 80-83. :doi:`10.2307/3001968`
3964 .. [5] Cureton, E.E., The Normal Approximation to the Signed-Rank
3965 Sampling Distribution When Zero Differences are Present,
3966 Journal of the American Statistical Association, Vol. 62, 1967,
3967 pp. 1068-1069. :doi:`10.1080/01621459.1967.10500917`
3969 Examples
3970 --------
3971 In [4]_, the differences in height between cross- and self-fertilized
3972 corn plants is given as follows:
3974 >>> d = [6, 8, 14, 16, 23, 24, 28, 29, 41, -48, 49, 56, 60, -67, 75]
3976 Cross-fertilized plants appear to be higher. To test the null
3977 hypothesis that there is no height difference, we can apply the
3978 two-sided test:
3980 >>> from scipy.stats import wilcoxon
3981 >>> res = wilcoxon(d)
3982 >>> res.statistic, res.pvalue
3983 (24.0, 0.041259765625)
3985 Hence, we would reject the null hypothesis at a confidence level of 5%,
3986 concluding that there is a difference in height between the groups.
3987 To confirm that the median of the differences can be assumed to be
3988 positive, we use:
3990 >>> res = wilcoxon(d, alternative='greater')
3991 >>> res.statistic, res.pvalue
3992 (96.0, 0.0206298828125)
3994 This shows that the null hypothesis that the median is negative can be
3995 rejected at a confidence level of 5% in favor of the alternative that
3996 the median is greater than zero. The p-values above are exact. Using the
3997 normal approximation gives very similar values:
3999 >>> res = wilcoxon(d, method='approx')
4000 >>> res.statistic, res.pvalue
4001 (24.0, 0.04088813291185591)
4003 Note that the statistic changed to 96 in the one-sided case (the sum
4004 of ranks of positive differences) whereas it is 24 in the two-sided
4005 case (the minimum of sum of ranks above and below zero).
4007 In the example above, the differences in height between paired plants are
4008 provided to `wilcoxon` directly. Alternatively, `wilcoxon` accepts two
4009 samples of equal length, calculates the differences between paired
4010 elements, then performs the test. Consider the samples ``x`` and ``y``:
4012 >>> import numpy as np
4013 >>> x = np.array([0.5, 0.825, 0.375, 0.5])
4014 >>> y = np.array([0.525, 0.775, 0.325, 0.55])
4015 >>> res = wilcoxon(x, y, alternative='greater')
4016 >>> res
4017 WilcoxonResult(statistic=5.0, pvalue=0.5625)
4019 Note that had we calculated the differences by hand, the test would have
4020 produced different results:
4022 >>> d = [-0.025, 0.05, 0.05, -0.05]
4023 >>> ref = wilcoxon(d, alternative='greater')
4024 >>> ref
4025 WilcoxonResult(statistic=6.0, pvalue=0.4375)
4027 The substantial difference is due to roundoff error in the results of
4028 ``x-y``:
4030 >>> d - (x-y)
4031 array([2.08166817e-17, 6.93889390e-17, 1.38777878e-17, 4.16333634e-17])
4033 Even though we expected all the elements of ``(x-y)[1:]`` to have the same
4034 magnitude ``0.05``, they have slightly different magnitudes in practice,
4035 and therefore are assigned different ranks in the test. Before performing
4036 the test, consider calculating ``d`` and adjusting it as necessary to
4037 ensure that theoretically identically values are not numerically distinct.
4038 For example:
4040 >>> d2 = np.around(x - y, decimals=3)
4041 >>> wilcoxon(d2, alternative='greater')
4042 WilcoxonResult(statistic=6.0, pvalue=0.4375)
4044 """
4045 mode = method
4047 if mode not in ["auto", "approx", "exact"]:
4048 raise ValueError("mode must be either 'auto', 'approx' or 'exact'")
4050 if zero_method not in ["wilcox", "pratt", "zsplit"]:
4051 raise ValueError("Zero method must be either 'wilcox' "
4052 "or 'pratt' or 'zsplit'")
4054 if alternative not in ["two-sided", "less", "greater"]:
4055 raise ValueError("Alternative must be either 'two-sided', "
4056 "'greater' or 'less'")
4058 if y is None:
4059 d = asarray(x)
4060 if d.ndim > 1:
4061 raise ValueError('Sample x must be one-dimensional.')
4062 else:
4063 x, y = map(asarray, (x, y))
4064 if x.ndim > 1 or y.ndim > 1:
4065 raise ValueError('Samples x and y must be one-dimensional.')
4066 if len(x) != len(y):
4067 raise ValueError('The samples x and y must have the same length.')
4068 # Future enhancement: consider warning when elements of `d` appear to
4069 # be tied but are numerically distinct.
4070 d = x - y
4072 if len(d) == 0:
4073 NaN = _get_nan(d)
4074 res = WilcoxonResult(NaN, NaN)
4075 if method == 'approx':
4076 res.zstatistic = NaN
4077 return res
4079 if mode == "auto":
4080 if len(d) <= 50:
4081 mode = "exact"
4082 else:
4083 mode = "approx"
4085 n_zero = np.sum(d == 0)
4086 if n_zero > 0 and mode == "exact":
4087 mode = "approx"
4088 warnings.warn("Exact p-value calculation does not work if there are "
4089 "zeros. Switching to normal approximation.")
4091 if mode == "approx":
4092 if zero_method in ["wilcox", "pratt"]:
4093 if n_zero == len(d):
4094 raise ValueError("zero_method 'wilcox' and 'pratt' do not "
4095 "work if x - y is zero for all elements.")
4096 if zero_method == "wilcox":
4097 # Keep all non-zero differences
4098 d = compress(np.not_equal(d, 0), d)
4100 count = len(d)
4101 if count < 10 and mode == "approx":
4102 warnings.warn("Sample size too small for normal approximation.")
4104 r = _stats_py.rankdata(abs(d))
4105 r_plus = np.sum((d > 0) * r)
4106 r_minus = np.sum((d < 0) * r)
4108 if zero_method == "zsplit":
4109 r_zero = np.sum((d == 0) * r)
4110 r_plus += r_zero / 2.
4111 r_minus += r_zero / 2.
4113 # return min for two-sided test, but r_plus for one-sided test
4114 # the literature is not consistent here
4115 # r_plus is more informative since r_plus + r_minus = count*(count+1)/2,
4116 # i.e. the sum of the ranks, so r_minus and the min can be inferred
4117 # (If alternative='pratt', r_plus + r_minus = count*(count+1)/2 - r_zero.)
4118 # [3] uses the r_plus for the one-sided test, keep min for two-sided test
4119 # to keep backwards compatibility
4120 if alternative == "two-sided":
4121 T = min(r_plus, r_minus)
4122 else:
4123 T = r_plus
4125 if mode == "approx":
4126 mn = count * (count + 1.) * 0.25
4127 se = count * (count + 1.) * (2. * count + 1.)
4129 if zero_method == "pratt":
4130 r = r[d != 0]
4131 # normal approximation needs to be adjusted, see Cureton (1967)
4132 mn -= n_zero * (n_zero + 1.) * 0.25
4133 se -= n_zero * (n_zero + 1.) * (2. * n_zero + 1.)
4135 replist, repnum = find_repeats(r)
4136 if repnum.size != 0:
4137 # Correction for repeated elements.
4138 se -= 0.5 * (repnum * (repnum * repnum - 1)).sum()
4140 se = sqrt(se / 24)
4142 # apply continuity correction if applicable
4143 d = 0
4144 if correction:
4145 if alternative == "two-sided":
4146 d = 0.5 * np.sign(T - mn)
4147 elif alternative == "less":
4148 d = -0.5
4149 else:
4150 d = 0.5
4152 # compute statistic and p-value using normal approximation
4153 z = (T - mn - d) / se
4154 if alternative == "two-sided":
4155 prob = 2. * distributions.norm.sf(abs(z))
4156 elif alternative == "greater":
4157 # large T = r_plus indicates x is greater than y; i.e.
4158 # accept alternative in that case and return small p-value (sf)
4159 prob = distributions.norm.sf(z)
4160 else:
4161 prob = distributions.norm.cdf(z)
4162 elif mode == "exact":
4163 # get pmf of the possible positive ranksums r_plus
4164 pmf = _get_wilcoxon_distr(count)
4165 # note: r_plus is int (ties not allowed), need int for slices below
4166 r_plus = int(r_plus)
4167 if alternative == "two-sided":
4168 if r_plus == (len(pmf) - 1) // 2:
4169 # r_plus is the center of the distribution.
4170 prob = 1.0
4171 else:
4172 p_less = np.sum(pmf[:r_plus + 1])
4173 p_greater = np.sum(pmf[r_plus:])
4174 prob = 2*min(p_greater, p_less)
4175 elif alternative == "greater":
4176 prob = np.sum(pmf[r_plus:])
4177 else:
4178 prob = np.sum(pmf[:r_plus + 1])
4179 prob = np.clip(prob, 0, 1)
4181 res = WilcoxonResult(T, prob)
4182 if method == 'approx':
4183 res.zstatistic = z
4184 return res
4187MedianTestResult = _make_tuple_bunch(
4188 'MedianTestResult',
4189 ['statistic', 'pvalue', 'median', 'table'], []
4190)
4193def median_test(*samples, ties='below', correction=True, lambda_=1,
4194 nan_policy='propagate'):
4195 """Perform a Mood's median test.
4197 Test that two or more samples come from populations with the same median.
4199 Let ``n = len(samples)`` be the number of samples. The "grand median" of
4200 all the data is computed, and a contingency table is formed by
4201 classifying the values in each sample as being above or below the grand
4202 median. The contingency table, along with `correction` and `lambda_`,
4203 are passed to `scipy.stats.chi2_contingency` to compute the test statistic
4204 and p-value.
4206 Parameters
4207 ----------
4208 sample1, sample2, ... : array_like
4209 The set of samples. There must be at least two samples.
4210 Each sample must be a one-dimensional sequence containing at least
4211 one value. The samples are not required to have the same length.
4212 ties : str, optional
4213 Determines how values equal to the grand median are classified in
4214 the contingency table. The string must be one of::
4216 "below":
4217 Values equal to the grand median are counted as "below".
4218 "above":
4219 Values equal to the grand median are counted as "above".
4220 "ignore":
4221 Values equal to the grand median are not counted.
4223 The default is "below".
4224 correction : bool, optional
4225 If True, *and* there are just two samples, apply Yates' correction
4226 for continuity when computing the test statistic associated with
4227 the contingency table. Default is True.
4228 lambda_ : float or str, optional
4229 By default, the statistic computed in this test is Pearson's
4230 chi-squared statistic. `lambda_` allows a statistic from the
4231 Cressie-Read power divergence family to be used instead. See
4232 `power_divergence` for details.
4233 Default is 1 (Pearson's chi-squared statistic).
4234 nan_policy : {'propagate', 'raise', 'omit'}, optional
4235 Defines how to handle when input contains nan. 'propagate' returns nan,
4236 'raise' throws an error, 'omit' performs the calculations ignoring nan
4237 values. Default is 'propagate'.
4239 Returns
4240 -------
4241 res : MedianTestResult
4242 An object containing attributes:
4244 statistic : float
4245 The test statistic. The statistic that is returned is determined
4246 by `lambda_`. The default is Pearson's chi-squared statistic.
4247 pvalue : float
4248 The p-value of the test.
4249 median : float
4250 The grand median.
4251 table : ndarray
4252 The contingency table. The shape of the table is (2, n), where
4253 n is the number of samples. The first row holds the counts of the
4254 values above the grand median, and the second row holds the counts
4255 of the values below the grand median. The table allows further
4256 analysis with, for example, `scipy.stats.chi2_contingency`, or with
4257 `scipy.stats.fisher_exact` if there are two samples, without having
4258 to recompute the table. If ``nan_policy`` is "propagate" and there
4259 are nans in the input, the return value for ``table`` is ``None``.
4261 See Also
4262 --------
4263 kruskal : Compute the Kruskal-Wallis H-test for independent samples.
4264 mannwhitneyu : Computes the Mann-Whitney rank test on samples x and y.
4266 Notes
4267 -----
4268 .. versionadded:: 0.15.0
4270 References
4271 ----------
4272 .. [1] Mood, A. M., Introduction to the Theory of Statistics. McGraw-Hill
4273 (1950), pp. 394-399.
4274 .. [2] Zar, J. H., Biostatistical Analysis, 5th ed. Prentice Hall (2010).
4275 See Sections 8.12 and 10.15.
4277 Examples
4278 --------
4279 A biologist runs an experiment in which there are three groups of plants.
4280 Group 1 has 16 plants, group 2 has 15 plants, and group 3 has 17 plants.
4281 Each plant produces a number of seeds. The seed counts for each group
4282 are::
4284 Group 1: 10 14 14 18 20 22 24 25 31 31 32 39 43 43 48 49
4285 Group 2: 28 30 31 33 34 35 36 40 44 55 57 61 91 92 99
4286 Group 3: 0 3 9 22 23 25 25 33 34 34 40 45 46 48 62 67 84
4288 The following code applies Mood's median test to these samples.
4290 >>> g1 = [10, 14, 14, 18, 20, 22, 24, 25, 31, 31, 32, 39, 43, 43, 48, 49]
4291 >>> g2 = [28, 30, 31, 33, 34, 35, 36, 40, 44, 55, 57, 61, 91, 92, 99]
4292 >>> g3 = [0, 3, 9, 22, 23, 25, 25, 33, 34, 34, 40, 45, 46, 48, 62, 67, 84]
4293 >>> from scipy.stats import median_test
4294 >>> res = median_test(g1, g2, g3)
4296 The median is
4298 >>> res.median
4299 34.0
4301 and the contingency table is
4303 >>> res.table
4304 array([[ 5, 10, 7],
4305 [11, 5, 10]])
4307 `p` is too large to conclude that the medians are not the same:
4309 >>> res.pvalue
4310 0.12609082774093244
4312 The "G-test" can be performed by passing ``lambda_="log-likelihood"`` to
4313 `median_test`.
4315 >>> res = median_test(g1, g2, g3, lambda_="log-likelihood")
4316 >>> res.pvalue
4317 0.12224779737117837
4319 The median occurs several times in the data, so we'll get a different
4320 result if, for example, ``ties="above"`` is used:
4322 >>> res = median_test(g1, g2, g3, ties="above")
4323 >>> res.pvalue
4324 0.063873276069553273
4326 >>> res.table
4327 array([[ 5, 11, 9],
4328 [11, 4, 8]])
4330 This example demonstrates that if the data set is not large and there
4331 are values equal to the median, the p-value can be sensitive to the
4332 choice of `ties`.
4334 """
4335 if len(samples) < 2:
4336 raise ValueError('median_test requires two or more samples.')
4338 ties_options = ['below', 'above', 'ignore']
4339 if ties not in ties_options:
4340 raise ValueError("invalid 'ties' option '{}'; 'ties' must be one "
4341 "of: {}".format(ties, str(ties_options)[1:-1]))
4343 data = [np.asarray(sample) for sample in samples]
4345 # Validate the sizes and shapes of the arguments.
4346 for k, d in enumerate(data):
4347 if d.size == 0:
4348 raise ValueError("Sample %d is empty. All samples must "
4349 "contain at least one value." % (k + 1))
4350 if d.ndim != 1:
4351 raise ValueError("Sample %d has %d dimensions. All "
4352 "samples must be one-dimensional sequences." %
4353 (k + 1, d.ndim))
4355 cdata = np.concatenate(data)
4356 contains_nan, nan_policy = _contains_nan(cdata, nan_policy)
4357 if contains_nan and nan_policy == 'propagate':
4358 return MedianTestResult(np.nan, np.nan, np.nan, None)
4360 if contains_nan:
4361 grand_median = np.median(cdata[~np.isnan(cdata)])
4362 else:
4363 grand_median = np.median(cdata)
4364 # When the minimum version of numpy supported by scipy is 1.9.0,
4365 # the above if/else statement can be replaced by the single line:
4366 # grand_median = np.nanmedian(cdata)
4368 # Create the contingency table.
4369 table = np.zeros((2, len(data)), dtype=np.int64)
4370 for k, sample in enumerate(data):
4371 sample = sample[~np.isnan(sample)]
4373 nabove = count_nonzero(sample > grand_median)
4374 nbelow = count_nonzero(sample < grand_median)
4375 nequal = sample.size - (nabove + nbelow)
4376 table[0, k] += nabove
4377 table[1, k] += nbelow
4378 if ties == "below":
4379 table[1, k] += nequal
4380 elif ties == "above":
4381 table[0, k] += nequal
4383 # Check that no row or column of the table is all zero.
4384 # Such a table can not be given to chi2_contingency, because it would have
4385 # a zero in the table of expected frequencies.
4386 rowsums = table.sum(axis=1)
4387 if rowsums[0] == 0:
4388 raise ValueError("All values are below the grand median (%r)." %
4389 grand_median)
4390 if rowsums[1] == 0:
4391 raise ValueError("All values are above the grand median (%r)." %
4392 grand_median)
4393 if ties == "ignore":
4394 # We already checked that each sample has at least one value, but it
4395 # is possible that all those values equal the grand median. If `ties`
4396 # is "ignore", that would result in a column of zeros in `table`. We
4397 # check for that case here.
4398 zero_cols = np.nonzero((table == 0).all(axis=0))[0]
4399 if len(zero_cols) > 0:
4400 msg = ("All values in sample %d are equal to the grand "
4401 "median (%r), so they are ignored, resulting in an "
4402 "empty sample." % (zero_cols[0] + 1, grand_median))
4403 raise ValueError(msg)
4405 stat, p, dof, expected = chi2_contingency(table, lambda_=lambda_,
4406 correction=correction)
4407 return MedianTestResult(stat, p, grand_median, table)
4410def _circfuncs_common(samples, high, low, nan_policy='propagate'):
4411 # Ensure samples are array-like and size is not zero
4412 samples = np.asarray(samples)
4413 if samples.size == 0:
4414 return np.nan, np.asarray(np.nan), np.asarray(np.nan), None
4416 # Recast samples as radians that range between 0 and 2 pi and calculate
4417 # the sine and cosine
4418 sin_samp = sin((samples - low)*2.*pi / (high - low))
4419 cos_samp = cos((samples - low)*2.*pi / (high - low))
4421 # Apply the NaN policy
4422 contains_nan, nan_policy = _contains_nan(samples, nan_policy)
4423 if contains_nan and nan_policy == 'omit':
4424 mask = np.isnan(samples)
4425 # Set the sines and cosines that are NaN to zero
4426 sin_samp[mask] = 0.0
4427 cos_samp[mask] = 0.0
4428 else:
4429 mask = None
4431 return samples, sin_samp, cos_samp, mask
4434def circmean(samples, high=2*pi, low=0, axis=None, nan_policy='propagate'):
4435 """Compute the circular mean for samples in a range.
4437 Parameters
4438 ----------
4439 samples : array_like
4440 Input array.
4441 high : float or int, optional
4442 High boundary for the sample range. Default is ``2*pi``.
4443 low : float or int, optional
4444 Low boundary for the sample range. Default is 0.
4445 axis : int, optional
4446 Axis along which means are computed. The default is to compute
4447 the mean of the flattened array.
4448 nan_policy : {'propagate', 'raise', 'omit'}, optional
4449 Defines how to handle when input contains nan. 'propagate' returns nan,
4450 'raise' throws an error, 'omit' performs the calculations ignoring nan
4451 values. Default is 'propagate'.
4453 Returns
4454 -------
4455 circmean : float
4456 Circular mean.
4458 See Also
4459 --------
4460 circstd : Circular standard deviation.
4461 circvar : Circular variance.
4463 Examples
4464 --------
4465 For simplicity, all angles are printed out in degrees.
4467 >>> import numpy as np
4468 >>> from scipy.stats import circmean
4469 >>> import matplotlib.pyplot as plt
4470 >>> angles = np.deg2rad(np.array([20, 30, 330]))
4471 >>> circmean = circmean(angles)
4472 >>> np.rad2deg(circmean)
4473 7.294976657784009
4475 >>> mean = angles.mean()
4476 >>> np.rad2deg(mean)
4477 126.66666666666666
4479 Plot and compare the circular mean against the arithmetic mean.
4481 >>> plt.plot(np.cos(np.linspace(0, 2*np.pi, 500)),
4482 ... np.sin(np.linspace(0, 2*np.pi, 500)),
4483 ... c='k')
4484 >>> plt.scatter(np.cos(angles), np.sin(angles), c='k')
4485 >>> plt.scatter(np.cos(circmean), np.sin(circmean), c='b',
4486 ... label='circmean')
4487 >>> plt.scatter(np.cos(mean), np.sin(mean), c='r', label='mean')
4488 >>> plt.legend()
4489 >>> plt.axis('equal')
4490 >>> plt.show()
4492 """
4493 samples, sin_samp, cos_samp, nmask = _circfuncs_common(samples, high, low,
4494 nan_policy=nan_policy)
4495 sin_sum = sin_samp.sum(axis=axis)
4496 cos_sum = cos_samp.sum(axis=axis)
4497 res = arctan2(sin_sum, cos_sum)
4499 mask_nan = ~np.isnan(res)
4500 if mask_nan.ndim > 0:
4501 mask = res[mask_nan] < 0
4502 else:
4503 mask = res < 0
4505 if mask.ndim > 0:
4506 mask_nan[mask_nan] = mask
4507 res[mask_nan] += 2*pi
4508 elif mask:
4509 res += 2*pi
4511 # Set output to NaN if no samples went into the mean
4512 if nmask is not None:
4513 if nmask.all():
4514 res = np.full(shape=res.shape, fill_value=np.nan)
4515 else:
4516 # Find out if any of the axis that are being averaged consist
4517 # entirely of NaN. If one exists, set the result (res) to NaN
4518 nshape = 0 if axis is None else axis
4519 smask = nmask.shape[nshape] == nmask.sum(axis=axis)
4520 if smask.any():
4521 res[smask] = np.nan
4523 return res*(high - low)/2.0/pi + low
4526def circvar(samples, high=2*pi, low=0, axis=None, nan_policy='propagate'):
4527 """Compute the circular variance for samples assumed to be in a range.
4529 Parameters
4530 ----------
4531 samples : array_like
4532 Input array.
4533 high : float or int, optional
4534 High boundary for the sample range. Default is ``2*pi``.
4535 low : float or int, optional
4536 Low boundary for the sample range. Default is 0.
4537 axis : int, optional
4538 Axis along which variances are computed. The default is to compute
4539 the variance of the flattened array.
4540 nan_policy : {'propagate', 'raise', 'omit'}, optional
4541 Defines how to handle when input contains nan. 'propagate' returns nan,
4542 'raise' throws an error, 'omit' performs the calculations ignoring nan
4543 values. Default is 'propagate'.
4545 Returns
4546 -------
4547 circvar : float
4548 Circular variance.
4550 See Also
4551 --------
4552 circmean : Circular mean.
4553 circstd : Circular standard deviation.
4555 Notes
4556 -----
4557 This uses the following definition of circular variance: ``1-R``, where
4558 ``R`` is the mean resultant vector. The
4559 returned value is in the range [0, 1], 0 standing for no variance, and 1
4560 for a large variance. In the limit of small angles, this value is similar
4561 to half the 'linear' variance.
4563 References
4564 ----------
4565 .. [1] Fisher, N.I. *Statistical analysis of circular data*. Cambridge
4566 University Press, 1993.
4568 Examples
4569 --------
4570 >>> import numpy as np
4571 >>> from scipy.stats import circvar
4572 >>> import matplotlib.pyplot as plt
4573 >>> samples_1 = np.array([0.072, -0.158, 0.077, 0.108, 0.286,
4574 ... 0.133, -0.473, -0.001, -0.348, 0.131])
4575 >>> samples_2 = np.array([0.111, -0.879, 0.078, 0.733, 0.421,
4576 ... 0.104, -0.136, -0.867, 0.012, 0.105])
4577 >>> circvar_1 = circvar(samples_1)
4578 >>> circvar_2 = circvar(samples_2)
4580 Plot the samples.
4582 >>> fig, (left, right) = plt.subplots(ncols=2)
4583 >>> for image in (left, right):
4584 ... image.plot(np.cos(np.linspace(0, 2*np.pi, 500)),
4585 ... np.sin(np.linspace(0, 2*np.pi, 500)),
4586 ... c='k')
4587 ... image.axis('equal')
4588 ... image.axis('off')
4589 >>> left.scatter(np.cos(samples_1), np.sin(samples_1), c='k', s=15)
4590 >>> left.set_title(f"circular variance: {np.round(circvar_1, 2)!r}")
4591 >>> right.scatter(np.cos(samples_2), np.sin(samples_2), c='k', s=15)
4592 >>> right.set_title(f"circular variance: {np.round(circvar_2, 2)!r}")
4593 >>> plt.show()
4595 """
4596 samples, sin_samp, cos_samp, mask = _circfuncs_common(samples, high, low,
4597 nan_policy=nan_policy)
4598 if mask is None:
4599 sin_mean = sin_samp.mean(axis=axis)
4600 cos_mean = cos_samp.mean(axis=axis)
4601 else:
4602 nsum = np.asarray(np.sum(~mask, axis=axis).astype(float))
4603 nsum[nsum == 0] = np.nan
4604 sin_mean = sin_samp.sum(axis=axis) / nsum
4605 cos_mean = cos_samp.sum(axis=axis) / nsum
4606 # hypot can go slightly above 1 due to rounding errors
4607 with np.errstate(invalid='ignore'):
4608 R = np.minimum(1, hypot(sin_mean, cos_mean))
4610 res = 1. - R
4611 return res
4614def circstd(samples, high=2*pi, low=0, axis=None, nan_policy='propagate', *,
4615 normalize=False):
4616 """
4617 Compute the circular standard deviation for samples assumed to be in the
4618 range [low to high].
4620 Parameters
4621 ----------
4622 samples : array_like
4623 Input array.
4624 high : float or int, optional
4625 High boundary for the sample range. Default is ``2*pi``.
4626 low : float or int, optional
4627 Low boundary for the sample range. Default is 0.
4628 axis : int, optional
4629 Axis along which standard deviations are computed. The default is
4630 to compute the standard deviation of the flattened array.
4631 nan_policy : {'propagate', 'raise', 'omit'}, optional
4632 Defines how to handle when input contains nan. 'propagate' returns nan,
4633 'raise' throws an error, 'omit' performs the calculations ignoring nan
4634 values. Default is 'propagate'.
4635 normalize : boolean, optional
4636 If True, the returned value is equal to ``sqrt(-2*log(R))`` and does
4637 not depend on the variable units. If False (default), the returned
4638 value is scaled by ``((high-low)/(2*pi))``.
4640 Returns
4641 -------
4642 circstd : float
4643 Circular standard deviation.
4645 See Also
4646 --------
4647 circmean : Circular mean.
4648 circvar : Circular variance.
4650 Notes
4651 -----
4652 This uses a definition of circular standard deviation from [1]_.
4653 Essentially, the calculation is as follows.
4655 .. code-block:: python
4657 import numpy as np
4658 C = np.cos(samples).mean()
4659 S = np.sin(samples).mean()
4660 R = np.sqrt(C**2 + S**2)
4661 l = 2*np.pi / (high-low)
4662 circstd = np.sqrt(-2*np.log(R)) / l
4664 In the limit of small angles, it returns a number close to the 'linear'
4665 standard deviation.
4667 References
4668 ----------
4669 .. [1] Mardia, K. V. (1972). 2. In *Statistics of Directional Data*
4670 (pp. 18-24). Academic Press. :doi:`10.1016/C2013-0-07425-7`.
4672 Examples
4673 --------
4674 >>> import numpy as np
4675 >>> from scipy.stats import circstd
4676 >>> import matplotlib.pyplot as plt
4677 >>> samples_1 = np.array([0.072, -0.158, 0.077, 0.108, 0.286,
4678 ... 0.133, -0.473, -0.001, -0.348, 0.131])
4679 >>> samples_2 = np.array([0.111, -0.879, 0.078, 0.733, 0.421,
4680 ... 0.104, -0.136, -0.867, 0.012, 0.105])
4681 >>> circstd_1 = circstd(samples_1)
4682 >>> circstd_2 = circstd(samples_2)
4684 Plot the samples.
4686 >>> fig, (left, right) = plt.subplots(ncols=2)
4687 >>> for image in (left, right):
4688 ... image.plot(np.cos(np.linspace(0, 2*np.pi, 500)),
4689 ... np.sin(np.linspace(0, 2*np.pi, 500)),
4690 ... c='k')
4691 ... image.axis('equal')
4692 ... image.axis('off')
4693 >>> left.scatter(np.cos(samples_1), np.sin(samples_1), c='k', s=15)
4694 >>> left.set_title(f"circular std: {np.round(circstd_1, 2)!r}")
4695 >>> right.plot(np.cos(np.linspace(0, 2*np.pi, 500)),
4696 ... np.sin(np.linspace(0, 2*np.pi, 500)),
4697 ... c='k')
4698 >>> right.scatter(np.cos(samples_2), np.sin(samples_2), c='k', s=15)
4699 >>> right.set_title(f"circular std: {np.round(circstd_2, 2)!r}")
4700 >>> plt.show()
4702 """
4703 samples, sin_samp, cos_samp, mask = _circfuncs_common(samples, high, low,
4704 nan_policy=nan_policy)
4705 if mask is None:
4706 sin_mean = sin_samp.mean(axis=axis) # [1] (2.2.3)
4707 cos_mean = cos_samp.mean(axis=axis) # [1] (2.2.3)
4708 else:
4709 nsum = np.asarray(np.sum(~mask, axis=axis).astype(float))
4710 nsum[nsum == 0] = np.nan
4711 sin_mean = sin_samp.sum(axis=axis) / nsum
4712 cos_mean = cos_samp.sum(axis=axis) / nsum
4713 # hypot can go slightly above 1 due to rounding errors
4714 with np.errstate(invalid='ignore'):
4715 R = np.minimum(1, hypot(sin_mean, cos_mean)) # [1] (2.2.4)
4717 res = sqrt(-2*log(R))
4718 if not normalize:
4719 res *= (high-low)/(2.*pi) # [1] (2.3.14) w/ (2.3.7)
4720 return res
4723class DirectionalStats:
4724 def __init__(self, mean_direction, mean_resultant_length):
4725 self.mean_direction = mean_direction
4726 self.mean_resultant_length = mean_resultant_length
4728 def __repr__(self):
4729 return (f"DirectionalStats(mean_direction={self.mean_direction},"
4730 f" mean_resultant_length={self.mean_resultant_length})")
4733def directional_stats(samples, *, axis=0, normalize=True):
4734 """
4735 Computes sample statistics for directional data.
4737 Computes the directional mean (also called the mean direction vector) and
4738 mean resultant length of a sample of vectors.
4740 The directional mean is a measure of "preferred direction" of vector data.
4741 It is analogous to the sample mean, but it is for use when the length of
4742 the data is irrelevant (e.g. unit vectors).
4744 The mean resultant length is a value between 0 and 1 used to quantify the
4745 dispersion of directional data: the smaller the mean resultant length, the
4746 greater the dispersion. Several definitions of directional variance
4747 involving the mean resultant length are given in [1]_ and [2]_.
4749 Parameters
4750 ----------
4751 samples : array_like
4752 Input array. Must be at least two-dimensional, and the last axis of the
4753 input must correspond with the dimensionality of the vector space.
4754 When the input is exactly two dimensional, this means that each row
4755 of the data is a vector observation.
4756 axis : int, default: 0
4757 Axis along which the directional mean is computed.
4758 normalize: boolean, default: True
4759 If True, normalize the input to ensure that each observation is a
4760 unit vector. It the observations are already unit vectors, consider
4761 setting this to False to avoid unnecessary computation.
4763 Returns
4764 -------
4765 res : DirectionalStats
4766 An object containing attributes:
4768 mean_direction : ndarray
4769 Directional mean.
4770 mean_resultant_length : ndarray
4771 The mean resultant length [1]_.
4773 See Also
4774 --------
4775 circmean: circular mean; i.e. directional mean for 2D *angles*
4776 circvar: circular variance; i.e. directional variance for 2D *angles*
4778 Notes
4779 -----
4780 This uses a definition of directional mean from [1]_.
4781 Assuming the observations are unit vectors, the calculation is as follows.
4783 .. code-block:: python
4785 mean = samples.mean(axis=0)
4786 mean_resultant_length = np.linalg.norm(mean)
4787 mean_direction = mean / mean_resultant_length
4789 This definition is appropriate for *directional* data (i.e. vector data
4790 for which the magnitude of each observation is irrelevant) but not
4791 for *axial* data (i.e. vector data for which the magnitude and *sign* of
4792 each observation is irrelevant).
4794 Several definitions of directional variance involving the mean resultant
4795 length ``R`` have been proposed, including ``1 - R`` [1]_, ``1 - R**2``
4796 [2]_, and ``2 * (1 - R)`` [2]_. Rather than choosing one, this function
4797 returns ``R`` as attribute `mean_resultant_length` so the user can compute
4798 their preferred measure of dispersion.
4800 References
4801 ----------
4802 .. [1] Mardia, Jupp. (2000). *Directional Statistics*
4803 (p. 163). Wiley.
4805 .. [2] https://en.wikipedia.org/wiki/Directional_statistics
4807 Examples
4808 --------
4809 >>> import numpy as np
4810 >>> from scipy.stats import directional_stats
4811 >>> data = np.array([[3, 4], # first observation, 2D vector space
4812 ... [6, -8]]) # second observation
4813 >>> dirstats = directional_stats(data)
4814 >>> dirstats.mean_direction
4815 array([1., 0.])
4817 In contrast, the regular sample mean of the vectors would be influenced
4818 by the magnitude of each observation. Furthermore, the result would not be
4819 a unit vector.
4821 >>> data.mean(axis=0)
4822 array([4.5, -2.])
4824 An exemplary use case for `directional_stats` is to find a *meaningful*
4825 center for a set of observations on a sphere, e.g. geographical locations.
4827 >>> data = np.array([[0.8660254, 0.5, 0.],
4828 ... [0.8660254, -0.5, 0.]])
4829 >>> dirstats = directional_stats(data)
4830 >>> dirstats.mean_direction
4831 array([1., 0., 0.])
4833 The regular sample mean on the other hand yields a result which does not
4834 lie on the surface of the sphere.
4836 >>> data.mean(axis=0)
4837 array([0.8660254, 0., 0.])
4839 The function also returns the mean resultant length, which
4840 can be used to calculate a directional variance. For example, using the
4841 definition ``Var(z) = 1 - R`` from [2]_ where ``R`` is the
4842 mean resultant length, we can calculate the directional variance of the
4843 vectors in the above example as:
4845 >>> 1 - dirstats.mean_resultant_length
4846 0.13397459716167093
4847 """
4848 samples = np.asarray(samples)
4849 if samples.ndim < 2:
4850 raise ValueError("samples must at least be two-dimensional. "
4851 f"Instead samples has shape: {samples.shape!r}")
4852 samples = np.moveaxis(samples, axis, 0)
4853 if normalize:
4854 vectornorms = np.linalg.norm(samples, axis=-1, keepdims=True)
4855 samples = samples/vectornorms
4856 mean = np.mean(samples, axis=0)
4857 mean_resultant_length = np.linalg.norm(mean, axis=-1, keepdims=True)
4858 mean_direction = mean / mean_resultant_length
4859 return DirectionalStats(mean_direction,
4860 mean_resultant_length.squeeze(-1)[()])
4863def false_discovery_control(ps, *, axis=0, method='bh'):
4864 """Adjust p-values to control the false discovery rate.
4866 The false discovery rate (FDR) is the expected proportion of rejected null
4867 hypotheses that are actually true.
4868 If the null hypothesis is rejected when the *adjusted* p-value falls below
4869 a specified level, the false discovery rate is controlled at that level.
4871 Parameters
4872 ----------
4873 ps : 1D array_like
4874 The p-values to adjust. Elements must be real numbers between 0 and 1.
4875 axis : int
4876 The axis along which to perform the adjustment. The adjustment is
4877 performed independently along each axis-slice. If `axis` is None, `ps`
4878 is raveled before performing the adjustment.
4879 method : {'bh', 'by'}
4880 The false discovery rate control procedure to apply: ``'bh'`` is for
4881 Benjamini-Hochberg [1]_ (Eq. 1), ``'by'`` is for Benjaminini-Yekutieli
4882 [2]_ (Theorem 1.3). The latter is more conservative, but it is
4883 guaranteed to control the FDR even when the p-values are not from
4884 independent tests.
4886 Returns
4887 -------
4888 ps_adusted : array_like
4889 The adjusted p-values. If the null hypothesis is rejected where these
4890 fall below a specified level, the false discovery rate is controlled
4891 at that level.
4893 See Also
4894 --------
4895 combine_pvalues
4896 statsmodels.stats.multitest.multipletests
4898 Notes
4899 -----
4900 In multiple hypothesis testing, false discovery control procedures tend to
4901 offer higher power than familywise error rate control procedures (e.g.
4902 Bonferroni correction [1]_).
4904 If the p-values correspond with independent tests (or tests with
4905 "positive regression dependencies" [2]_), rejecting null hypotheses
4906 corresponding with Benjamini-Hochberg-adjusted p-values below :math:`q`
4907 controls the false discovery rate at a level less than or equal to
4908 :math:`q m_0 / m`, where :math:`m_0` is the number of true null hypotheses
4909 and :math:`m` is the total number of null hypotheses tested. The same is
4910 true even for dependent tests when the p-values are adjusted accorded to
4911 the more conservative Benjaminini-Yekutieli procedure.
4913 The adjusted p-values produced by this function are comparable to those
4914 produced by the R function ``p.adjust`` and the statsmodels function
4915 `statsmodels.stats.multitest.multipletests`. Please consider the latter
4916 for more advanced methods of multiple comparison correction.
4918 References
4919 ----------
4920 .. [1] Benjamini, Yoav, and Yosef Hochberg. "Controlling the false
4921 discovery rate: a practical and powerful approach to multiple
4922 testing." Journal of the Royal statistical society: series B
4923 (Methodological) 57.1 (1995): 289-300.
4925 .. [2] Benjamini, Yoav, and Daniel Yekutieli. "The control of the false
4926 discovery rate in multiple testing under dependency." Annals of
4927 statistics (2001): 1165-1188.
4929 .. [3] TileStats. FDR - Benjamini-Hochberg explained - Youtube.
4930 https://www.youtube.com/watch?v=rZKa4tW2NKs.
4932 .. [4] Neuhaus, Karl-Ludwig, et al. "Improved thrombolysis in acute
4933 myocardial infarction with front-loaded administration of alteplase:
4934 results of the rt-PA-APSAC patency study (TAPS)." Journal of the
4935 American College of Cardiology 19.5 (1992): 885-891.
4937 Examples
4938 --------
4939 We follow the example from [1]_.
4941 Thrombolysis with recombinant tissue-type plasminogen activator (rt-PA)
4942 and anisoylated plasminogen streptokinase activator (APSAC) in
4943 myocardial infarction has been proved to reduce mortality. [4]_
4944 investigated the effects of a new front-loaded administration of rt-PA
4945 versus those obtained with a standard regimen of APSAC, in a randomized
4946 multicentre trial in 421 patients with acute myocardial infarction.
4948 There were four families of hypotheses tested in the study, the last of
4949 which was "cardiac and other events after the start of thrombolitic
4950 treatment". FDR control may be desired in this family of hypotheses
4951 because it would not be appropriate to conclude that the front-loaded
4952 treatment is better if it is merely equivalent to the previous treatment.
4954 The p-values corresponding with the 15 hypotheses in this family were
4956 >>> ps = [0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344,
4957 ... 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000]
4959 If the chosen significance level is 0.05, we may be tempted to reject the
4960 null hypotheses for the tests corresponding with the first nine p-values,
4961 as the first nine p-values fall below the chosen significance level.
4962 However, this would ignore the problem of "multiplicity": if we fail to
4963 correct for the fact that multiple comparisons are being performed, we
4964 are more likely to incorrectly reject true null hypotheses.
4966 One approach to the multiplicity problem is to control the family-wise
4967 error rate (FWER), that is, the rate at which the null hypothesis is
4968 rejected when it is actually true. A common procedure of this kind is the
4969 Bonferroni correction [1]_. We begin by multiplying the p-values by the
4970 number of hypotheses tested.
4972 >>> import numpy as np
4973 >>> np.array(ps) * len(ps)
4974 array([1.5000e-03, 6.0000e-03, 2.8500e-02, 1.4250e-01, 3.0150e-01,
4975 4.1700e-01, 4.4700e-01, 5.1600e-01, 6.8850e-01, 4.8600e+00,
4976 6.3930e+00, 8.5785e+00, 9.7920e+00, 1.1385e+01, 1.5000e+01])
4978 To control the FWER at 5%, we reject only the hypotheses corresponding
4979 with adjusted p-values less than 0.05. In this case, only the hypotheses
4980 corresponding with the first three p-values can be rejected. According to
4981 [1]_, these three hypotheses concerned "allergic reaction" and "two
4982 different aspects of bleeding."
4984 An alternative approach is to control the false discovery rate: the
4985 expected fraction of rejected null hypotheses that are actually true. The
4986 advantage of this approach is that it typically affords greater power: an
4987 increased rate of rejecting the null hypothesis when it is indeed false. To
4988 control the false discovery rate at 5%, we apply the Benjamini-Hochberg
4989 p-value adjustment.
4991 >>> from scipy import stats
4992 >>> stats.false_discovery_control(ps)
4993 array([0.0015 , 0.003 , 0.0095 , 0.035625 , 0.0603 ,
4994 0.06385714, 0.06385714, 0.0645 , 0.0765 , 0.486 ,
4995 0.58118182, 0.714875 , 0.75323077, 0.81321429, 1. ])
4997 Now, the first *four* adjusted p-values fall below 0.05, so we would reject
4998 the null hypotheses corresponding with these *four* p-values. Rejection
4999 of the fourth null hypothesis was particularly important to the original
5000 study as it led to the conclusion that the new treatment had a
5001 "substantially lower in-hospital mortality rate."
5003 """
5004 # Input Validation and Special Cases
5005 ps = np.asarray(ps)
5007 ps_in_range = (np.issubdtype(ps.dtype, np.number)
5008 and np.all(ps == np.clip(ps, 0, 1)))
5009 if not ps_in_range:
5010 raise ValueError("`ps` must include only numbers between 0 and 1.")
5012 methods = {'bh', 'by'}
5013 if method.lower() not in methods:
5014 raise ValueError(f"Unrecognized `method` '{method}'."
5015 f"Method must be one of {methods}.")
5016 method = method.lower()
5018 if axis is None:
5019 axis = 0
5020 ps = ps.ravel()
5022 axis = np.asarray(axis)[()]
5023 if not np.issubdtype(axis.dtype, np.integer) or axis.size != 1:
5024 raise ValueError("`axis` must be an integer or `None`")
5026 if ps.size <= 1 or ps.shape[axis] <= 1:
5027 return ps[()]
5029 ps = np.moveaxis(ps, axis, -1)
5030 m = ps.shape[-1]
5032 # Main Algorithm
5033 # Equivalent to the ideas of [1] and [2], except that this adjusts the
5034 # p-values as described in [3]. The results are similar to those produced
5035 # by R's p.adjust.
5037 # "Let [ps] be the ordered observed p-values..."
5038 order = np.argsort(ps, axis=-1)
5039 ps = np.take_along_axis(ps, order, axis=-1) # this copies ps
5041 # Equation 1 of [1] rearranged to reject when p is less than specified q
5042 i = np.arange(1, m+1)
5043 ps *= m / i
5045 # Theorem 1.3 of [2]
5046 if method == 'by':
5047 ps *= np.sum(1 / i)
5049 # accounts for rejecting all null hypotheses i for i < k, where k is
5050 # defined in Eq. 1 of either [1] or [2]. See [3]. Starting with the index j
5051 # of the second to last element, we replace element j with element j+1 if
5052 # the latter is smaller.
5053 np.minimum.accumulate(ps[..., ::-1], out=ps[..., ::-1], axis=-1)
5055 # Restore original order of axes and data
5056 np.put_along_axis(ps, order, values=ps.copy(), axis=-1)
5057 ps = np.moveaxis(ps, -1, axis)
5059 return np.clip(ps, 0, 1)