Coverage for /usr/lib/python3/dist-packages/scipy/stats/_mstats_basic.py: 11%
991 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2An extension of scipy.stats._stats_py to support masked arrays
4"""
5# Original author (2007): Pierre GF Gerard-Marchant
8__all__ = ['argstoarray',
9 'count_tied_groups',
10 'describe',
11 'f_oneway', 'find_repeats','friedmanchisquare',
12 'kendalltau','kendalltau_seasonal','kruskal','kruskalwallis',
13 'ks_twosamp', 'ks_2samp', 'kurtosis', 'kurtosistest',
14 'ks_1samp', 'kstest',
15 'linregress',
16 'mannwhitneyu', 'meppf','mode','moment','mquantiles','msign',
17 'normaltest',
18 'obrientransform',
19 'pearsonr','plotting_positions','pointbiserialr',
20 'rankdata',
21 'scoreatpercentile','sem',
22 'sen_seasonal_slopes','skew','skewtest','spearmanr',
23 'siegelslopes', 'theilslopes',
24 'tmax','tmean','tmin','trim','trimboth',
25 'trimtail','trima','trimr','trimmed_mean','trimmed_std',
26 'trimmed_stde','trimmed_var','tsem','ttest_1samp','ttest_onesamp',
27 'ttest_ind','ttest_rel','tvar',
28 'variation',
29 'winsorize',
30 'brunnermunzel',
31 ]
33import numpy as np
34from numpy import ndarray
35import numpy.ma as ma
36from numpy.ma import masked, nomask
37import math
39import itertools
40import warnings
41from collections import namedtuple
43from . import distributions
44from scipy._lib._util import _rename_parameter, _contains_nan
45from scipy._lib._bunch import _make_tuple_bunch
46import scipy.special as special
47import scipy.stats._stats_py
49from ._stats_mstats_common import (
50 _find_repeats,
51 linregress as stats_linregress,
52 LinregressResult as stats_LinregressResult,
53 theilslopes as stats_theilslopes,
54 siegelslopes as stats_siegelslopes
55 )
58def _chk_asarray(a, axis):
59 # Always returns a masked array, raveled for axis=None
60 a = ma.asanyarray(a)
61 if axis is None:
62 a = ma.ravel(a)
63 outaxis = 0
64 else:
65 outaxis = axis
66 return a, outaxis
69def _chk2_asarray(a, b, axis):
70 a = ma.asanyarray(a)
71 b = ma.asanyarray(b)
72 if axis is None:
73 a = ma.ravel(a)
74 b = ma.ravel(b)
75 outaxis = 0
76 else:
77 outaxis = axis
78 return a, b, outaxis
81def _chk_size(a, b):
82 a = ma.asanyarray(a)
83 b = ma.asanyarray(b)
84 (na, nb) = (a.size, b.size)
85 if na != nb:
86 raise ValueError("The size of the input array should match!"
87 " ({} <> {})".format(na, nb))
88 return (a, b, na)
91def argstoarray(*args):
92 """
93 Constructs a 2D array from a group of sequences.
95 Sequences are filled with missing values to match the length of the longest
96 sequence.
98 Parameters
99 ----------
100 *args : sequences
101 Group of sequences.
103 Returns
104 -------
105 argstoarray : MaskedArray
106 A ( `m` x `n` ) masked array, where `m` is the number of arguments and
107 `n` the length of the longest argument.
109 Notes
110 -----
111 `numpy.ma.vstack` has identical behavior, but is called with a sequence
112 of sequences.
114 Examples
115 --------
116 A 2D masked array constructed from a group of sequences is returned.
118 >>> from scipy.stats.mstats import argstoarray
119 >>> argstoarray([1, 2, 3], [4, 5, 6])
120 masked_array(
121 data=[[1.0, 2.0, 3.0],
122 [4.0, 5.0, 6.0]],
123 mask=[[False, False, False],
124 [False, False, False]],
125 fill_value=1e+20)
127 The returned masked array filled with missing values when the lengths of
128 sequences are different.
130 >>> argstoarray([1, 3], [4, 5, 6])
131 masked_array(
132 data=[[1.0, 3.0, --],
133 [4.0, 5.0, 6.0]],
134 mask=[[False, False, True],
135 [False, False, False]],
136 fill_value=1e+20)
138 """
139 if len(args) == 1 and not isinstance(args[0], ndarray):
140 output = ma.asarray(args[0])
141 if output.ndim != 2:
142 raise ValueError("The input should be 2D")
143 else:
144 n = len(args)
145 m = max([len(k) for k in args])
146 output = ma.array(np.empty((n,m), dtype=float), mask=True)
147 for (k,v) in enumerate(args):
148 output[k,:len(v)] = v
150 output[np.logical_not(np.isfinite(output._data))] = masked
151 return output
154def find_repeats(arr):
155 """Find repeats in arr and return a tuple (repeats, repeat_count).
157 The input is cast to float64. Masked values are discarded.
159 Parameters
160 ----------
161 arr : sequence
162 Input array. The array is flattened if it is not 1D.
164 Returns
165 -------
166 repeats : ndarray
167 Array of repeated values.
168 counts : ndarray
169 Array of counts.
171 Examples
172 --------
173 >>> from scipy.stats import mstats
174 >>> mstats.find_repeats([2, 1, 2, 3, 2, 2, 5])
175 (array([2.]), array([4]))
177 In the above example, 2 repeats 4 times.
179 >>> mstats.find_repeats([[10, 20, 1, 2], [5, 5, 4, 4]])
180 (array([4., 5.]), array([2, 2]))
182 In the above example, both 4 and 5 repeat 2 times.
184 """
185 # Make sure we get a copy. ma.compressed promises a "new array", but can
186 # actually return a reference.
187 compr = np.asarray(ma.compressed(arr), dtype=np.float64)
188 try:
189 need_copy = np.may_share_memory(compr, arr)
190 except AttributeError:
191 # numpy < 1.8.2 bug: np.may_share_memory([], []) raises,
192 # while in numpy 1.8.2 and above it just (correctly) returns False.
193 need_copy = False
194 if need_copy:
195 compr = compr.copy()
196 return _find_repeats(compr)
199def count_tied_groups(x, use_missing=False):
200 """
201 Counts the number of tied values.
203 Parameters
204 ----------
205 x : sequence
206 Sequence of data on which to counts the ties
207 use_missing : bool, optional
208 Whether to consider missing values as tied.
210 Returns
211 -------
212 count_tied_groups : dict
213 Returns a dictionary (nb of ties: nb of groups).
215 Examples
216 --------
217 >>> from scipy.stats import mstats
218 >>> import numpy as np
219 >>> z = [0, 0, 0, 2, 2, 2, 3, 3, 4, 5, 6]
220 >>> mstats.count_tied_groups(z)
221 {2: 1, 3: 2}
223 In the above example, the ties were 0 (3x), 2 (3x) and 3 (2x).
225 >>> z = np.ma.array([0, 0, 1, 2, 2, 2, 3, 3, 4, 5, 6])
226 >>> mstats.count_tied_groups(z)
227 {2: 2, 3: 1}
228 >>> z[[1,-1]] = np.ma.masked
229 >>> mstats.count_tied_groups(z, use_missing=True)
230 {2: 2, 3: 1}
232 """
233 nmasked = ma.getmask(x).sum()
234 # We need the copy as find_repeats will overwrite the initial data
235 data = ma.compressed(x).copy()
236 (ties, counts) = find_repeats(data)
237 nties = {}
238 if len(ties):
239 nties = dict(zip(np.unique(counts), itertools.repeat(1)))
240 nties.update(dict(zip(*find_repeats(counts))))
242 if nmasked and use_missing:
243 try:
244 nties[nmasked] += 1
245 except KeyError:
246 nties[nmasked] = 1
248 return nties
251def rankdata(data, axis=None, use_missing=False):
252 """Returns the rank (also known as order statistics) of each data point
253 along the given axis.
255 If some values are tied, their rank is averaged.
256 If some values are masked, their rank is set to 0 if use_missing is False,
257 or set to the average rank of the unmasked values if use_missing is True.
259 Parameters
260 ----------
261 data : sequence
262 Input data. The data is transformed to a masked array
263 axis : {None,int}, optional
264 Axis along which to perform the ranking.
265 If None, the array is first flattened. An exception is raised if
266 the axis is specified for arrays with a dimension larger than 2
267 use_missing : bool, optional
268 Whether the masked values have a rank of 0 (False) or equal to the
269 average rank of the unmasked values (True).
271 """
272 def _rank1d(data, use_missing=False):
273 n = data.count()
274 rk = np.empty(data.size, dtype=float)
275 idx = data.argsort()
276 rk[idx[:n]] = np.arange(1,n+1)
278 if use_missing:
279 rk[idx[n:]] = (n+1)/2.
280 else:
281 rk[idx[n:]] = 0
283 repeats = find_repeats(data.copy())
284 for r in repeats[0]:
285 condition = (data == r).filled(False)
286 rk[condition] = rk[condition].mean()
287 return rk
289 data = ma.array(data, copy=False)
290 if axis is None:
291 if data.ndim > 1:
292 return _rank1d(data.ravel(), use_missing).reshape(data.shape)
293 else:
294 return _rank1d(data, use_missing)
295 else:
296 return ma.apply_along_axis(_rank1d,axis,data,use_missing).view(ndarray)
299ModeResult = namedtuple('ModeResult', ('mode', 'count'))
302def mode(a, axis=0):
303 """
304 Returns an array of the modal (most common) value in the passed array.
306 Parameters
307 ----------
308 a : array_like
309 n-dimensional array of which to find mode(s).
310 axis : int or None, optional
311 Axis along which to operate. Default is 0. If None, compute over
312 the whole array `a`.
314 Returns
315 -------
316 mode : ndarray
317 Array of modal values.
318 count : ndarray
319 Array of counts for each mode.
321 Notes
322 -----
323 For more details, see `scipy.stats.mode`.
325 Examples
326 --------
327 >>> import numpy as np
328 >>> from scipy import stats
329 >>> from scipy.stats import mstats
330 >>> m_arr = np.ma.array([1, 1, 0, 0, 0, 0], mask=[0, 0, 1, 1, 1, 0])
331 >>> mstats.mode(m_arr) # note that most zeros are masked
332 ModeResult(mode=array([1.]), count=array([2.]))
334 """
335 return _mode(a, axis=axis, keepdims=True)
338def _mode(a, axis=0, keepdims=True):
339 # Don't want to expose `keepdims` from the public `mstats.mode`
340 a, axis = _chk_asarray(a, axis)
342 def _mode1D(a):
343 (rep,cnt) = find_repeats(a)
344 if not cnt.ndim:
345 return (0, 0)
346 elif cnt.size:
347 return (rep[cnt.argmax()], cnt.max())
348 else:
349 return (a.min(), 1)
351 if axis is None:
352 output = _mode1D(ma.ravel(a))
353 output = (ma.array(output[0]), ma.array(output[1]))
354 else:
355 output = ma.apply_along_axis(_mode1D, axis, a)
356 if keepdims is None or keepdims:
357 newshape = list(a.shape)
358 newshape[axis] = 1
359 slices = [slice(None)] * output.ndim
360 slices[axis] = 0
361 modes = output[tuple(slices)].reshape(newshape)
362 slices[axis] = 1
363 counts = output[tuple(slices)].reshape(newshape)
364 output = (modes, counts)
365 else:
366 output = np.moveaxis(output, axis, 0)
368 return ModeResult(*output)
371def _betai(a, b, x):
372 x = np.asanyarray(x)
373 x = ma.where(x < 1.0, x, 1.0) # if x > 1 then return 1.0
374 return special.betainc(a, b, x)
377def msign(x):
378 """Returns the sign of x, or 0 if x is masked."""
379 return ma.filled(np.sign(x), 0)
382def pearsonr(x, y):
383 r"""
384 Pearson correlation coefficient and p-value for testing non-correlation.
386 The Pearson correlation coefficient [1]_ measures the linear relationship
387 between two datasets. The calculation of the p-value relies on the
388 assumption that each dataset is normally distributed. (See Kowalski [3]_
389 for a discussion of the effects of non-normality of the input on the
390 distribution of the correlation coefficient.) Like other correlation
391 coefficients, this one varies between -1 and +1 with 0 implying no
392 correlation. Correlations of -1 or +1 imply an exact linear relationship.
394 Parameters
395 ----------
396 x : (N,) array_like
397 Input array.
398 y : (N,) array_like
399 Input array.
401 Returns
402 -------
403 r : float
404 Pearson's correlation coefficient.
405 p-value : float
406 Two-tailed p-value.
408 Warns
409 -----
410 PearsonRConstantInputWarning
411 Raised if an input is a constant array. The correlation coefficient
412 is not defined in this case, so ``np.nan`` is returned.
414 PearsonRNearConstantInputWarning
415 Raised if an input is "nearly" constant. The array ``x`` is considered
416 nearly constant if ``norm(x - mean(x)) < 1e-13 * abs(mean(x))``.
417 Numerical errors in the calculation ``x - mean(x)`` in this case might
418 result in an inaccurate calculation of r.
420 See Also
421 --------
422 spearmanr : Spearman rank-order correlation coefficient.
423 kendalltau : Kendall's tau, a correlation measure for ordinal data.
425 Notes
426 -----
427 The correlation coefficient is calculated as follows:
429 .. math::
431 r = \frac{\sum (x - m_x) (y - m_y)}
432 {\sqrt{\sum (x - m_x)^2 \sum (y - m_y)^2}}
434 where :math:`m_x` is the mean of the vector x and :math:`m_y` is
435 the mean of the vector y.
437 Under the assumption that x and y are drawn from
438 independent normal distributions (so the population correlation coefficient
439 is 0), the probability density function of the sample correlation
440 coefficient r is ([1]_, [2]_):
442 .. math::
444 f(r) = \frac{{(1-r^2)}^{n/2-2}}{\mathrm{B}(\frac{1}{2},\frac{n}{2}-1)}
446 where n is the number of samples, and B is the beta function. This
447 is sometimes referred to as the exact distribution of r. This is
448 the distribution that is used in `pearsonr` to compute the p-value.
449 The distribution is a beta distribution on the interval [-1, 1],
450 with equal shape parameters a = b = n/2 - 1. In terms of SciPy's
451 implementation of the beta distribution, the distribution of r is::
453 dist = scipy.stats.beta(n/2 - 1, n/2 - 1, loc=-1, scale=2)
455 The p-value returned by `pearsonr` is a two-sided p-value. The p-value
456 roughly indicates the probability of an uncorrelated system
457 producing datasets that have a Pearson correlation at least as extreme
458 as the one computed from these datasets. More precisely, for a
459 given sample with correlation coefficient r, the p-value is
460 the probability that abs(r') of a random sample x' and y' drawn from
461 the population with zero correlation would be greater than or equal
462 to abs(r). In terms of the object ``dist`` shown above, the p-value
463 for a given r and length n can be computed as::
465 p = 2*dist.cdf(-abs(r))
467 When n is 2, the above continuous distribution is not well-defined.
468 One can interpret the limit of the beta distribution as the shape
469 parameters a and b approach a = b = 0 as a discrete distribution with
470 equal probability masses at r = 1 and r = -1. More directly, one
471 can observe that, given the data x = [x1, x2] and y = [y1, y2], and
472 assuming x1 != x2 and y1 != y2, the only possible values for r are 1
473 and -1. Because abs(r') for any sample x' and y' with length 2 will
474 be 1, the two-sided p-value for a sample of length 2 is always 1.
476 References
477 ----------
478 .. [1] "Pearson correlation coefficient", Wikipedia,
479 https://en.wikipedia.org/wiki/Pearson_correlation_coefficient
480 .. [2] Student, "Probable error of a correlation coefficient",
481 Biometrika, Volume 6, Issue 2-3, 1 September 1908, pp. 302-310.
482 .. [3] C. J. Kowalski, "On the Effects of Non-Normality on the Distribution
483 of the Sample Product-Moment Correlation Coefficient"
484 Journal of the Royal Statistical Society. Series C (Applied
485 Statistics), Vol. 21, No. 1 (1972), pp. 1-12.
487 Examples
488 --------
489 >>> import numpy as np
490 >>> from scipy import stats
491 >>> from scipy.stats import mstats
492 >>> mstats.pearsonr([1, 2, 3, 4, 5], [10, 9, 2.5, 6, 4])
493 (-0.7426106572325057, 0.1505558088534455)
495 There is a linear dependence between x and y if y = a + b*x + e, where
496 a,b are constants and e is a random error term, assumed to be independent
497 of x. For simplicity, assume that x is standard normal, a=0, b=1 and let
498 e follow a normal distribution with mean zero and standard deviation s>0.
500 >>> s = 0.5
501 >>> x = stats.norm.rvs(size=500)
502 >>> e = stats.norm.rvs(scale=s, size=500)
503 >>> y = x + e
504 >>> mstats.pearsonr(x, y)
505 (0.9029601878969703, 8.428978827629898e-185) # may vary
507 This should be close to the exact value given by
509 >>> 1/np.sqrt(1 + s**2)
510 0.8944271909999159
512 For s=0.5, we observe a high level of correlation. In general, a large
513 variance of the noise reduces the correlation, while the correlation
514 approaches one as the variance of the error goes to zero.
516 It is important to keep in mind that no correlation does not imply
517 independence unless (x, y) is jointly normal. Correlation can even be zero
518 when there is a very simple dependence structure: if X follows a
519 standard normal distribution, let y = abs(x). Note that the correlation
520 between x and y is zero. Indeed, since the expectation of x is zero,
521 cov(x, y) = E[x*y]. By definition, this equals E[x*abs(x)] which is zero
522 by symmetry. The following lines of code illustrate this observation:
524 >>> y = np.abs(x)
525 >>> mstats.pearsonr(x, y)
526 (-0.016172891856853524, 0.7182823678751942) # may vary
528 A non-zero correlation coefficient can be misleading. For example, if X has
529 a standard normal distribution, define y = x if x < 0 and y = 0 otherwise.
530 A simple calculation shows that corr(x, y) = sqrt(2/Pi) = 0.797...,
531 implying a high level of correlation:
533 >>> y = np.where(x < 0, x, 0)
534 >>> mstats.pearsonr(x, y)
535 (0.8537091583771509, 3.183461621422181e-143) # may vary
537 This is unintuitive since there is no dependence of x and y if x is larger
538 than zero which happens in about half of the cases if we sample x and y.
539 """
540 (x, y, n) = _chk_size(x, y)
541 (x, y) = (x.ravel(), y.ravel())
542 # Get the common mask and the total nb of unmasked elements
543 m = ma.mask_or(ma.getmask(x), ma.getmask(y))
544 n -= m.sum()
545 df = n-2
546 if df < 0:
547 return (masked, masked)
549 return scipy.stats._stats_py.pearsonr(
550 ma.masked_array(x, mask=m).compressed(),
551 ma.masked_array(y, mask=m).compressed())
554def spearmanr(x, y=None, use_ties=True, axis=None, nan_policy='propagate',
555 alternative='two-sided'):
556 """
557 Calculates a Spearman rank-order correlation coefficient and the p-value
558 to test for non-correlation.
560 The Spearman correlation is a nonparametric measure of the linear
561 relationship between two datasets. Unlike the Pearson correlation, the
562 Spearman correlation does not assume that both datasets are normally
563 distributed. Like other correlation coefficients, this one varies
564 between -1 and +1 with 0 implying no correlation. Correlations of -1 or
565 +1 imply a monotonic relationship. Positive correlations imply that
566 as `x` increases, so does `y`. Negative correlations imply that as `x`
567 increases, `y` decreases.
569 Missing values are discarded pair-wise: if a value is missing in `x`, the
570 corresponding value in `y` is masked.
572 The p-value roughly indicates the probability of an uncorrelated system
573 producing datasets that have a Spearman correlation at least as extreme
574 as the one computed from these datasets. The p-values are not entirely
575 reliable but are probably reasonable for datasets larger than 500 or so.
577 Parameters
578 ----------
579 x, y : 1D or 2D array_like, y is optional
580 One or two 1-D or 2-D arrays containing multiple variables and
581 observations. When these are 1-D, each represents a vector of
582 observations of a single variable. For the behavior in the 2-D case,
583 see under ``axis``, below.
584 use_ties : bool, optional
585 DO NOT USE. Does not do anything, keyword is only left in place for
586 backwards compatibility reasons.
587 axis : int or None, optional
588 If axis=0 (default), then each column represents a variable, with
589 observations in the rows. If axis=1, the relationship is transposed:
590 each row represents a variable, while the columns contain observations.
591 If axis=None, then both arrays will be raveled.
592 nan_policy : {'propagate', 'raise', 'omit'}, optional
593 Defines how to handle when input contains nan. 'propagate' returns nan,
594 'raise' throws an error, 'omit' performs the calculations ignoring nan
595 values. Default is 'propagate'.
596 alternative : {'two-sided', 'less', 'greater'}, optional
597 Defines the alternative hypothesis. Default is 'two-sided'.
598 The following options are available:
600 * 'two-sided': the correlation is nonzero
601 * 'less': the correlation is negative (less than zero)
602 * 'greater': the correlation is positive (greater than zero)
604 .. versionadded:: 1.7.0
606 Returns
607 -------
608 res : SignificanceResult
609 An object containing attributes:
611 statistic : float or ndarray (2-D square)
612 Spearman correlation matrix or correlation coefficient (if only 2
613 variables are given as parameters). Correlation matrix is square
614 with length equal to total number of variables (columns or rows) in
615 ``a`` and ``b`` combined.
616 pvalue : float
617 The p-value for a hypothesis test whose null hypothesis
618 is that two sets of data are linearly uncorrelated. See
619 `alternative` above for alternative hypotheses. `pvalue` has the
620 same shape as `statistic`.
622 References
623 ----------
624 [CRCProbStat2000] section 14.7
626 """
627 if not use_ties:
628 raise ValueError("`use_ties=False` is not supported in SciPy >= 1.2.0")
630 # Always returns a masked array, raveled if axis=None
631 x, axisout = _chk_asarray(x, axis)
632 if y is not None:
633 # Deal only with 2-D `x` case.
634 y, _ = _chk_asarray(y, axis)
635 if axisout == 0:
636 x = ma.column_stack((x, y))
637 else:
638 x = ma.row_stack((x, y))
640 if axisout == 1:
641 # To simplify the code that follow (always use `n_obs, n_vars` shape)
642 x = x.T
644 if nan_policy == 'omit':
645 x = ma.masked_invalid(x)
647 def _spearmanr_2cols(x):
648 # Mask the same observations for all variables, and then drop those
649 # observations (can't leave them masked, rankdata is weird).
650 x = ma.mask_rowcols(x, axis=0)
651 x = x[~x.mask.any(axis=1), :]
653 # If either column is entirely NaN or Inf
654 if not np.any(x.data):
655 res = scipy.stats._stats_py.SignificanceResult(np.nan, np.nan)
656 res.correlation = np.nan
657 return res
659 m = ma.getmask(x)
660 n_obs = x.shape[0]
661 dof = n_obs - 2 - int(m.sum(axis=0)[0])
662 if dof < 0:
663 raise ValueError("The input must have at least 3 entries!")
665 # Gets the ranks and rank differences
666 x_ranked = rankdata(x, axis=0)
667 rs = ma.corrcoef(x_ranked, rowvar=False).data
669 # rs can have elements equal to 1, so avoid zero division warnings
670 with np.errstate(divide='ignore'):
671 # clip the small negative values possibly caused by rounding
672 # errors before taking the square root
673 t = rs * np.sqrt((dof / ((rs+1.0) * (1.0-rs))).clip(0))
675 t, prob = scipy.stats._stats_py._ttest_finish(dof, t, alternative)
677 # For backwards compatibility, return scalars when comparing 2 columns
678 if rs.shape == (2, 2):
679 res = scipy.stats._stats_py.SignificanceResult(rs[1, 0],
680 prob[1, 0])
681 res.correlation = rs[1, 0]
682 return res
683 else:
684 res = scipy.stats._stats_py.SignificanceResult(rs, prob)
685 res.correlation = rs
686 return res
688 # Need to do this per pair of variables, otherwise the dropped observations
689 # in a third column mess up the result for a pair.
690 n_vars = x.shape[1]
691 if n_vars == 2:
692 return _spearmanr_2cols(x)
693 else:
694 rs = np.ones((n_vars, n_vars), dtype=float)
695 prob = np.zeros((n_vars, n_vars), dtype=float)
696 for var1 in range(n_vars - 1):
697 for var2 in range(var1+1, n_vars):
698 result = _spearmanr_2cols(x[:, [var1, var2]])
699 rs[var1, var2] = result.correlation
700 rs[var2, var1] = result.correlation
701 prob[var1, var2] = result.pvalue
702 prob[var2, var1] = result.pvalue
704 res = scipy.stats._stats_py.SignificanceResult(rs, prob)
705 res.correlation = rs
706 return res
709def _kendall_p_exact(n, c, alternative='two-sided'):
711 # Use the fact that distribution is symmetric: always calculate a CDF in
712 # the left tail.
713 # This will be the one-sided p-value if `c` is on the side of
714 # the null distribution predicted by the alternative hypothesis.
715 # The two-sided p-value will be twice this value.
716 # If `c` is on the other side of the null distribution, we'll need to
717 # take the complement and add back the probability mass at `c`.
718 in_right_tail = (c >= (n*(n-1))//2 - c)
719 alternative_greater = (alternative == 'greater')
720 c = int(min(c, (n*(n-1))//2 - c))
722 # Exact p-value, see Maurice G. Kendall, "Rank Correlation Methods"
723 # (4th Edition), Charles Griffin & Co., 1970.
724 if n <= 0:
725 raise ValueError(f'n ({n}) must be positive')
726 elif c < 0 or 4*c > n*(n-1):
727 raise ValueError(f'c ({c}) must satisfy 0 <= 4c <= n(n-1) = {n*(n-1)}.')
728 elif n == 1:
729 prob = 1.0
730 p_mass_at_c = 1
731 elif n == 2:
732 prob = 1.0
733 p_mass_at_c = 0.5
734 elif c == 0:
735 prob = 2.0/math.factorial(n) if n < 171 else 0.0
736 p_mass_at_c = prob/2
737 elif c == 1:
738 prob = 2.0/math.factorial(n-1) if n < 172 else 0.0
739 p_mass_at_c = (n-1)/math.factorial(n)
740 elif 4*c == n*(n-1) and alternative == 'two-sided':
741 # I'm sure there's a simple formula for p_mass_at_c in this
742 # case, but I don't know it. Use generic formula for one-sided p-value.
743 prob = 1.0
744 elif n < 171:
745 new = np.zeros(c+1)
746 new[0:2] = 1.0
747 for j in range(3,n+1):
748 new = np.cumsum(new)
749 if j <= c:
750 new[j:] -= new[:c+1-j]
751 prob = 2.0*np.sum(new)/math.factorial(n)
752 p_mass_at_c = new[-1]/math.factorial(n)
753 else:
754 new = np.zeros(c+1)
755 new[0:2] = 1.0
756 for j in range(3, n+1):
757 new = np.cumsum(new)/j
758 if j <= c:
759 new[j:] -= new[:c+1-j]
760 prob = np.sum(new)
761 p_mass_at_c = new[-1]/2
763 if alternative != 'two-sided':
764 # if the alternative hypothesis and alternative agree,
765 # one-sided p-value is half the two-sided p-value
766 if in_right_tail == alternative_greater:
767 prob /= 2
768 else:
769 prob = 1 - prob/2 + p_mass_at_c
771 prob = np.clip(prob, 0, 1)
773 return prob
776def kendalltau(x, y, use_ties=True, use_missing=False, method='auto',
777 alternative='two-sided'):
778 """
779 Computes Kendall's rank correlation tau on two variables *x* and *y*.
781 Parameters
782 ----------
783 x : sequence
784 First data list (for example, time).
785 y : sequence
786 Second data list.
787 use_ties : {True, False}, optional
788 Whether ties correction should be performed.
789 use_missing : {False, True}, optional
790 Whether missing data should be allocated a rank of 0 (False) or the
791 average rank (True)
792 method : {'auto', 'asymptotic', 'exact'}, optional
793 Defines which method is used to calculate the p-value [1]_.
794 'asymptotic' uses a normal approximation valid for large samples.
795 'exact' computes the exact p-value, but can only be used if no ties
796 are present. As the sample size increases, the 'exact' computation
797 time may grow and the result may lose some precision.
798 'auto' is the default and selects the appropriate
799 method based on a trade-off between speed and accuracy.
800 alternative : {'two-sided', 'less', 'greater'}, optional
801 Defines the alternative hypothesis. Default is 'two-sided'.
802 The following options are available:
804 * 'two-sided': the rank correlation is nonzero
805 * 'less': the rank correlation is negative (less than zero)
806 * 'greater': the rank correlation is positive (greater than zero)
808 Returns
809 -------
810 res : SignificanceResult
811 An object containing attributes:
813 statistic : float
814 The tau statistic.
815 pvalue : float
816 The p-value for a hypothesis test whose null hypothesis is
817 an absence of association, tau = 0.
819 References
820 ----------
821 .. [1] Maurice G. Kendall, "Rank Correlation Methods" (4th Edition),
822 Charles Griffin & Co., 1970.
824 """
825 (x, y, n) = _chk_size(x, y)
826 (x, y) = (x.flatten(), y.flatten())
827 m = ma.mask_or(ma.getmask(x), ma.getmask(y))
828 if m is not nomask:
829 x = ma.array(x, mask=m, copy=True)
830 y = ma.array(y, mask=m, copy=True)
831 # need int() here, otherwise numpy defaults to 32 bit
832 # integer on all Windows architectures, causing overflow.
833 # int() will keep it infinite precision.
834 n -= int(m.sum())
836 if n < 2:
837 res = scipy.stats._stats_py.SignificanceResult(np.nan, np.nan)
838 res.correlation = np.nan
839 return res
841 rx = ma.masked_equal(rankdata(x, use_missing=use_missing), 0)
842 ry = ma.masked_equal(rankdata(y, use_missing=use_missing), 0)
843 idx = rx.argsort()
844 (rx, ry) = (rx[idx], ry[idx])
845 C = np.sum([((ry[i+1:] > ry[i]) * (rx[i+1:] > rx[i])).filled(0).sum()
846 for i in range(len(ry)-1)], dtype=float)
847 D = np.sum([((ry[i+1:] < ry[i])*(rx[i+1:] > rx[i])).filled(0).sum()
848 for i in range(len(ry)-1)], dtype=float)
849 xties = count_tied_groups(x)
850 yties = count_tied_groups(y)
851 if use_ties:
852 corr_x = np.sum([v*k*(k-1) for (k,v) in xties.items()], dtype=float)
853 corr_y = np.sum([v*k*(k-1) for (k,v) in yties.items()], dtype=float)
854 denom = ma.sqrt((n*(n-1)-corr_x)/2. * (n*(n-1)-corr_y)/2.)
855 else:
856 denom = n*(n-1)/2.
857 tau = (C-D) / denom
859 if method == 'exact' and (xties or yties):
860 raise ValueError("Ties found, exact method cannot be used.")
862 if method == 'auto':
863 if (not xties and not yties) and (n <= 33 or min(C, n*(n-1)/2.0-C) <= 1):
864 method = 'exact'
865 else:
866 method = 'asymptotic'
868 if not xties and not yties and method == 'exact':
869 prob = _kendall_p_exact(n, C, alternative)
871 elif method == 'asymptotic':
872 var_s = n*(n-1)*(2*n+5)
873 if use_ties:
874 var_s -= np.sum([v*k*(k-1)*(2*k+5)*1. for (k,v) in xties.items()])
875 var_s -= np.sum([v*k*(k-1)*(2*k+5)*1. for (k,v) in yties.items()])
876 v1 = (np.sum([v*k*(k-1) for (k, v) in xties.items()], dtype=float) *
877 np.sum([v*k*(k-1) for (k, v) in yties.items()], dtype=float))
878 v1 /= 2.*n*(n-1)
879 if n > 2:
880 v2 = np.sum([v*k*(k-1)*(k-2) for (k,v) in xties.items()],
881 dtype=float) * \
882 np.sum([v*k*(k-1)*(k-2) for (k,v) in yties.items()],
883 dtype=float)
884 v2 /= 9.*n*(n-1)*(n-2)
885 else:
886 v2 = 0
887 else:
888 v1 = v2 = 0
890 var_s /= 18.
891 var_s += (v1 + v2)
892 z = (C-D)/np.sqrt(var_s)
893 _, prob = scipy.stats._stats_py._normtest_finish(z, alternative)
894 else:
895 raise ValueError("Unknown method "+str(method)+" specified, please "
896 "use auto, exact or asymptotic.")
898 res = scipy.stats._stats_py.SignificanceResult(tau, prob)
899 res.correlation = tau
900 return res
903def kendalltau_seasonal(x):
904 """
905 Computes a multivariate Kendall's rank correlation tau, for seasonal data.
907 Parameters
908 ----------
909 x : 2-D ndarray
910 Array of seasonal data, with seasons in columns.
912 """
913 x = ma.array(x, subok=True, copy=False, ndmin=2)
914 (n,m) = x.shape
915 n_p = x.count(0)
917 S_szn = sum(msign(x[i:]-x[i]).sum(0) for i in range(n))
918 S_tot = S_szn.sum()
920 n_tot = x.count()
921 ties = count_tied_groups(x.compressed())
922 corr_ties = sum(v*k*(k-1) for (k,v) in ties.items())
923 denom_tot = ma.sqrt(1.*n_tot*(n_tot-1)*(n_tot*(n_tot-1)-corr_ties))/2.
925 R = rankdata(x, axis=0, use_missing=True)
926 K = ma.empty((m,m), dtype=int)
927 covmat = ma.empty((m,m), dtype=float)
928 denom_szn = ma.empty(m, dtype=float)
929 for j in range(m):
930 ties_j = count_tied_groups(x[:,j].compressed())
931 corr_j = sum(v*k*(k-1) for (k,v) in ties_j.items())
932 cmb = n_p[j]*(n_p[j]-1)
933 for k in range(j,m,1):
934 K[j,k] = sum(msign((x[i:,j]-x[i,j])*(x[i:,k]-x[i,k])).sum()
935 for i in range(n))
936 covmat[j,k] = (K[j,k] + 4*(R[:,j]*R[:,k]).sum() -
937 n*(n_p[j]+1)*(n_p[k]+1))/3.
938 K[k,j] = K[j,k]
939 covmat[k,j] = covmat[j,k]
941 denom_szn[j] = ma.sqrt(cmb*(cmb-corr_j)) / 2.
943 var_szn = covmat.diagonal()
945 z_szn = msign(S_szn) * (abs(S_szn)-1) / ma.sqrt(var_szn)
946 z_tot_ind = msign(S_tot) * (abs(S_tot)-1) / ma.sqrt(var_szn.sum())
947 z_tot_dep = msign(S_tot) * (abs(S_tot)-1) / ma.sqrt(covmat.sum())
949 prob_szn = special.erfc(abs(z_szn)/np.sqrt(2))
950 prob_tot_ind = special.erfc(abs(z_tot_ind)/np.sqrt(2))
951 prob_tot_dep = special.erfc(abs(z_tot_dep)/np.sqrt(2))
953 chi2_tot = (z_szn*z_szn).sum()
954 chi2_trd = m * z_szn.mean()**2
955 output = {'seasonal tau': S_szn/denom_szn,
956 'global tau': S_tot/denom_tot,
957 'global tau (alt)': S_tot/denom_szn.sum(),
958 'seasonal p-value': prob_szn,
959 'global p-value (indep)': prob_tot_ind,
960 'global p-value (dep)': prob_tot_dep,
961 'chi2 total': chi2_tot,
962 'chi2 trend': chi2_trd,
963 }
964 return output
967PointbiserialrResult = namedtuple('PointbiserialrResult', ('correlation',
968 'pvalue'))
971def pointbiserialr(x, y):
972 """Calculates a point biserial correlation coefficient and its p-value.
974 Parameters
975 ----------
976 x : array_like of bools
977 Input array.
978 y : array_like
979 Input array.
981 Returns
982 -------
983 correlation : float
984 R value
985 pvalue : float
986 2-tailed p-value
988 Notes
989 -----
990 Missing values are considered pair-wise: if a value is missing in x,
991 the corresponding value in y is masked.
993 For more details on `pointbiserialr`, see `scipy.stats.pointbiserialr`.
995 """
996 x = ma.fix_invalid(x, copy=True).astype(bool)
997 y = ma.fix_invalid(y, copy=True).astype(float)
998 # Get rid of the missing data
999 m = ma.mask_or(ma.getmask(x), ma.getmask(y))
1000 if m is not nomask:
1001 unmask = np.logical_not(m)
1002 x = x[unmask]
1003 y = y[unmask]
1005 n = len(x)
1006 # phat is the fraction of x values that are True
1007 phat = x.sum() / float(n)
1008 y0 = y[~x] # y-values where x is False
1009 y1 = y[x] # y-values where x is True
1010 y0m = y0.mean()
1011 y1m = y1.mean()
1013 rpb = (y1m - y0m)*np.sqrt(phat * (1-phat)) / y.std()
1015 df = n-2
1016 t = rpb*ma.sqrt(df/(1.0-rpb**2))
1017 prob = _betai(0.5*df, 0.5, df/(df+t*t))
1019 return PointbiserialrResult(rpb, prob)
1022def linregress(x, y=None):
1023 r"""
1024 Linear regression calculation
1026 Note that the non-masked version is used, and that this docstring is
1027 replaced by the non-masked docstring + some info on missing data.
1029 """
1030 if y is None:
1031 x = ma.array(x)
1032 if x.shape[0] == 2:
1033 x, y = x
1034 elif x.shape[1] == 2:
1035 x, y = x.T
1036 else:
1037 raise ValueError("If only `x` is given as input, "
1038 "it has to be of shape (2, N) or (N, 2), "
1039 f"provided shape was {x.shape}")
1040 else:
1041 x = ma.array(x)
1042 y = ma.array(y)
1044 x = x.flatten()
1045 y = y.flatten()
1047 if np.amax(x) == np.amin(x) and len(x) > 1:
1048 raise ValueError("Cannot calculate a linear regression "
1049 "if all x values are identical")
1051 m = ma.mask_or(ma.getmask(x), ma.getmask(y), shrink=False)
1052 if m is not nomask:
1053 x = ma.array(x, mask=m)
1054 y = ma.array(y, mask=m)
1055 if np.any(~m):
1056 result = stats_linregress(x.data[~m], y.data[~m])
1057 else:
1058 # All data is masked
1059 result = stats_LinregressResult(slope=None, intercept=None,
1060 rvalue=None, pvalue=None,
1061 stderr=None,
1062 intercept_stderr=None)
1063 else:
1064 result = stats_linregress(x.data, y.data)
1066 return result
1069def theilslopes(y, x=None, alpha=0.95, method='separate'):
1070 r"""
1071 Computes the Theil-Sen estimator for a set of points (x, y).
1073 `theilslopes` implements a method for robust linear regression. It
1074 computes the slope as the median of all slopes between paired values.
1076 Parameters
1077 ----------
1078 y : array_like
1079 Dependent variable.
1080 x : array_like or None, optional
1081 Independent variable. If None, use ``arange(len(y))`` instead.
1082 alpha : float, optional
1083 Confidence degree between 0 and 1. Default is 95% confidence.
1084 Note that `alpha` is symmetric around 0.5, i.e. both 0.1 and 0.9 are
1085 interpreted as "find the 90% confidence interval".
1086 method : {'joint', 'separate'}, optional
1087 Method to be used for computing estimate for intercept.
1088 Following methods are supported,
1090 * 'joint': Uses np.median(y - slope * x) as intercept.
1091 * 'separate': Uses np.median(y) - slope * np.median(x)
1092 as intercept.
1094 The default is 'separate'.
1096 .. versionadded:: 1.8.0
1098 Returns
1099 -------
1100 result : ``TheilslopesResult`` instance
1101 The return value is an object with the following attributes:
1103 slope : float
1104 Theil slope.
1105 intercept : float
1106 Intercept of the Theil line.
1107 low_slope : float
1108 Lower bound of the confidence interval on `slope`.
1109 high_slope : float
1110 Upper bound of the confidence interval on `slope`.
1112 See Also
1113 --------
1114 siegelslopes : a similar technique using repeated medians
1117 Notes
1118 -----
1119 For more details on `theilslopes`, see `scipy.stats.theilslopes`.
1121 """
1122 y = ma.asarray(y).flatten()
1123 if x is None:
1124 x = ma.arange(len(y), dtype=float)
1125 else:
1126 x = ma.asarray(x).flatten()
1127 if len(x) != len(y):
1128 raise ValueError(f"Incompatible lengths ! ({len(y)}<>{len(x)})")
1130 m = ma.mask_or(ma.getmask(x), ma.getmask(y))
1131 y._mask = x._mask = m
1132 # Disregard any masked elements of x or y
1133 y = y.compressed()
1134 x = x.compressed().astype(float)
1135 # We now have unmasked arrays so can use `scipy.stats.theilslopes`
1136 return stats_theilslopes(y, x, alpha=alpha, method=method)
1139def siegelslopes(y, x=None, method="hierarchical"):
1140 r"""
1141 Computes the Siegel estimator for a set of points (x, y).
1143 `siegelslopes` implements a method for robust linear regression
1144 using repeated medians to fit a line to the points (x, y).
1145 The method is robust to outliers with an asymptotic breakdown point
1146 of 50%.
1148 Parameters
1149 ----------
1150 y : array_like
1151 Dependent variable.
1152 x : array_like or None, optional
1153 Independent variable. If None, use ``arange(len(y))`` instead.
1154 method : {'hierarchical', 'separate'}
1155 If 'hierarchical', estimate the intercept using the estimated
1156 slope ``slope`` (default option).
1157 If 'separate', estimate the intercept independent of the estimated
1158 slope. See Notes for details.
1160 Returns
1161 -------
1162 result : ``SiegelslopesResult`` instance
1163 The return value is an object with the following attributes:
1165 slope : float
1166 Estimate of the slope of the regression line.
1167 intercept : float
1168 Estimate of the intercept of the regression line.
1170 See Also
1171 --------
1172 theilslopes : a similar technique without repeated medians
1174 Notes
1175 -----
1176 For more details on `siegelslopes`, see `scipy.stats.siegelslopes`.
1178 """
1179 y = ma.asarray(y).ravel()
1180 if x is None:
1181 x = ma.arange(len(y), dtype=float)
1182 else:
1183 x = ma.asarray(x).ravel()
1184 if len(x) != len(y):
1185 raise ValueError(f"Incompatible lengths ! ({len(y)}<>{len(x)})")
1187 m = ma.mask_or(ma.getmask(x), ma.getmask(y))
1188 y._mask = x._mask = m
1189 # Disregard any masked elements of x or y
1190 y = y.compressed()
1191 x = x.compressed().astype(float)
1192 # We now have unmasked arrays so can use `scipy.stats.siegelslopes`
1193 return stats_siegelslopes(y, x, method=method)
1196SenSeasonalSlopesResult = _make_tuple_bunch('SenSeasonalSlopesResult',
1197 ['intra_slope', 'inter_slope'])
1200def sen_seasonal_slopes(x):
1201 r"""
1202 Computes seasonal Theil-Sen and Kendall slope estimators.
1204 The seasonal generalization of Sen's slope computes the slopes between all
1205 pairs of values within a "season" (column) of a 2D array. It returns an
1206 array containing the median of these "within-season" slopes for each
1207 season (the Theil-Sen slope estimator of each season), and it returns the
1208 median of the within-season slopes across all seasons (the seasonal Kendall
1209 slope estimator).
1211 Parameters
1212 ----------
1213 x : 2D array_like
1214 Each column of `x` contains measurements of the dependent variable
1215 within a season. The independent variable (usually time) of each season
1216 is assumed to be ``np.arange(x.shape[0])``.
1218 Returns
1219 -------
1220 result : ``SenSeasonalSlopesResult`` instance
1221 The return value is an object with the following attributes:
1223 intra_slope : ndarray
1224 For each season, the Theil-Sen slope estimator: the median of
1225 within-season slopes.
1226 inter_slope : float
1227 The seasonal Kendall slope estimateor: the median of within-season
1228 slopes *across all* seasons.
1230 See Also
1231 --------
1232 theilslopes : the analogous function for non-seasonal data
1233 scipy.stats.theilslopes : non-seasonal slopes for non-masked arrays
1235 Notes
1236 -----
1237 The slopes :math:`d_{ijk}` within season :math:`i` are:
1239 .. math::
1241 d_{ijk} = \frac{x_{ij} - x_{ik}}
1242 {j - k}
1244 for pairs of distinct integer indices :math:`j, k` of :math:`x`.
1246 Element :math:`i` of the returned `intra_slope` array is the median of the
1247 :math:`d_{ijk}` over all :math:`j < k`; this is the Theil-Sen slope
1248 estimator of season :math:`i`. The returned `inter_slope` value, better
1249 known as the seasonal Kendall slope estimator, is the median of the
1250 :math:`d_{ijk}` over all :math:`i, j, k`.
1252 References
1253 ----------
1254 .. [1] Hirsch, Robert M., James R. Slack, and Richard A. Smith.
1255 "Techniques of trend analysis for monthly water quality data."
1256 *Water Resources Research* 18.1 (1982): 107-121.
1258 Examples
1259 --------
1260 Suppose we have 100 observations of a dependent variable for each of four
1261 seasons:
1263 >>> import numpy as np
1264 >>> rng = np.random.default_rng()
1265 >>> x = rng.random(size=(100, 4))
1267 We compute the seasonal slopes as:
1269 >>> from scipy import stats
1270 >>> intra_slope, inter_slope = stats.mstats.sen_seasonal_slopes(x)
1272 If we define a function to compute all slopes between observations within
1273 a season:
1275 >>> def dijk(yi):
1276 ... n = len(yi)
1277 ... x = np.arange(n)
1278 ... dy = yi - yi[:, np.newaxis]
1279 ... dx = x - x[:, np.newaxis]
1280 ... # we only want unique pairs of distinct indices
1281 ... mask = np.triu(np.ones((n, n), dtype=bool), k=1)
1282 ... return dy[mask]/dx[mask]
1284 then element ``i`` of ``intra_slope`` is the median of ``dijk[x[:, i]]``:
1286 >>> i = 2
1287 >>> np.allclose(np.median(dijk(x[:, i])), intra_slope[i])
1288 True
1290 and ``inter_slope`` is the median of the values returned by ``dijk`` for
1291 all seasons:
1293 >>> all_slopes = np.concatenate([dijk(x[:, i]) for i in range(x.shape[1])])
1294 >>> np.allclose(np.median(all_slopes), inter_slope)
1295 True
1297 Because the data are randomly generated, we would expect the median slopes
1298 to be nearly zero both within and across all seasons, and indeed they are:
1300 >>> intra_slope.data
1301 array([ 0.00124504, -0.00277761, -0.00221245, -0.00036338])
1302 >>> inter_slope
1303 -0.0010511779872922058
1305 """
1306 x = ma.array(x, subok=True, copy=False, ndmin=2)
1307 (n,_) = x.shape
1308 # Get list of slopes per season
1309 szn_slopes = ma.vstack([(x[i+1:]-x[i])/np.arange(1,n-i)[:,None]
1310 for i in range(n)])
1311 szn_medslopes = ma.median(szn_slopes, axis=0)
1312 medslope = ma.median(szn_slopes, axis=None)
1313 return SenSeasonalSlopesResult(szn_medslopes, medslope)
1316Ttest_1sampResult = namedtuple('Ttest_1sampResult', ('statistic', 'pvalue'))
1319def ttest_1samp(a, popmean, axis=0, alternative='two-sided'):
1320 """
1321 Calculates the T-test for the mean of ONE group of scores.
1323 Parameters
1324 ----------
1325 a : array_like
1326 sample observation
1327 popmean : float or array_like
1328 expected value in null hypothesis, if array_like than it must have the
1329 same shape as `a` excluding the axis dimension
1330 axis : int or None, optional
1331 Axis along which to compute test. If None, compute over the whole
1332 array `a`.
1333 alternative : {'two-sided', 'less', 'greater'}, optional
1334 Defines the alternative hypothesis.
1335 The following options are available (default is 'two-sided'):
1337 * 'two-sided': the mean of the underlying distribution of the sample
1338 is different than the given population mean (`popmean`)
1339 * 'less': the mean of the underlying distribution of the sample is
1340 less than the given population mean (`popmean`)
1341 * 'greater': the mean of the underlying distribution of the sample is
1342 greater than the given population mean (`popmean`)
1344 .. versionadded:: 1.7.0
1346 Returns
1347 -------
1348 statistic : float or array
1349 t-statistic
1350 pvalue : float or array
1351 The p-value
1353 Notes
1354 -----
1355 For more details on `ttest_1samp`, see `scipy.stats.ttest_1samp`.
1357 """
1358 a, axis = _chk_asarray(a, axis)
1359 if a.size == 0:
1360 return (np.nan, np.nan)
1362 x = a.mean(axis=axis)
1363 v = a.var(axis=axis, ddof=1)
1364 n = a.count(axis=axis)
1365 # force df to be an array for masked division not to throw a warning
1366 df = ma.asanyarray(n - 1.0)
1367 svar = ((n - 1.0) * v) / df
1368 with np.errstate(divide='ignore', invalid='ignore'):
1369 t = (x - popmean) / ma.sqrt(svar / n)
1371 t, prob = scipy.stats._stats_py._ttest_finish(df, t, alternative)
1372 return Ttest_1sampResult(t, prob)
1375ttest_onesamp = ttest_1samp
1378Ttest_indResult = namedtuple('Ttest_indResult', ('statistic', 'pvalue'))
1381def ttest_ind(a, b, axis=0, equal_var=True, alternative='two-sided'):
1382 """
1383 Calculates the T-test for the means of TWO INDEPENDENT samples of scores.
1385 Parameters
1386 ----------
1387 a, b : array_like
1388 The arrays must have the same shape, except in the dimension
1389 corresponding to `axis` (the first, by default).
1390 axis : int or None, optional
1391 Axis along which to compute test. If None, compute over the whole
1392 arrays, `a`, and `b`.
1393 equal_var : bool, optional
1394 If True, perform a standard independent 2 sample test that assumes equal
1395 population variances.
1396 If False, perform Welch's t-test, which does not assume equal population
1397 variance.
1399 .. versionadded:: 0.17.0
1400 alternative : {'two-sided', 'less', 'greater'}, optional
1401 Defines the alternative hypothesis.
1402 The following options are available (default is 'two-sided'):
1404 * 'two-sided': the means of the distributions underlying the samples
1405 are unequal.
1406 * 'less': the mean of the distribution underlying the first sample
1407 is less than the mean of the distribution underlying the second
1408 sample.
1409 * 'greater': the mean of the distribution underlying the first
1410 sample is greater than the mean of the distribution underlying
1411 the second sample.
1413 .. versionadded:: 1.7.0
1415 Returns
1416 -------
1417 statistic : float or array
1418 The calculated t-statistic.
1419 pvalue : float or array
1420 The p-value.
1422 Notes
1423 -----
1424 For more details on `ttest_ind`, see `scipy.stats.ttest_ind`.
1426 """
1427 a, b, axis = _chk2_asarray(a, b, axis)
1429 if a.size == 0 or b.size == 0:
1430 return Ttest_indResult(np.nan, np.nan)
1432 (x1, x2) = (a.mean(axis), b.mean(axis))
1433 (v1, v2) = (a.var(axis=axis, ddof=1), b.var(axis=axis, ddof=1))
1434 (n1, n2) = (a.count(axis), b.count(axis))
1436 if equal_var:
1437 # force df to be an array for masked division not to throw a warning
1438 df = ma.asanyarray(n1 + n2 - 2.0)
1439 svar = ((n1-1)*v1+(n2-1)*v2) / df
1440 denom = ma.sqrt(svar*(1.0/n1 + 1.0/n2)) # n-D computation here!
1441 else:
1442 vn1 = v1/n1
1443 vn2 = v2/n2
1444 with np.errstate(divide='ignore', invalid='ignore'):
1445 df = (vn1 + vn2)**2 / (vn1**2 / (n1 - 1) + vn2**2 / (n2 - 1))
1447 # If df is undefined, variances are zero.
1448 # It doesn't matter what df is as long as it is not NaN.
1449 df = np.where(np.isnan(df), 1, df)
1450 denom = ma.sqrt(vn1 + vn2)
1452 with np.errstate(divide='ignore', invalid='ignore'):
1453 t = (x1-x2) / denom
1455 t, prob = scipy.stats._stats_py._ttest_finish(df, t, alternative)
1456 return Ttest_indResult(t, prob)
1459Ttest_relResult = namedtuple('Ttest_relResult', ('statistic', 'pvalue'))
1462def ttest_rel(a, b, axis=0, alternative='two-sided'):
1463 """
1464 Calculates the T-test on TWO RELATED samples of scores, a and b.
1466 Parameters
1467 ----------
1468 a, b : array_like
1469 The arrays must have the same shape.
1470 axis : int or None, optional
1471 Axis along which to compute test. If None, compute over the whole
1472 arrays, `a`, and `b`.
1473 alternative : {'two-sided', 'less', 'greater'}, optional
1474 Defines the alternative hypothesis.
1475 The following options are available (default is 'two-sided'):
1477 * 'two-sided': the means of the distributions underlying the samples
1478 are unequal.
1479 * 'less': the mean of the distribution underlying the first sample
1480 is less than the mean of the distribution underlying the second
1481 sample.
1482 * 'greater': the mean of the distribution underlying the first
1483 sample is greater than the mean of the distribution underlying
1484 the second sample.
1486 .. versionadded:: 1.7.0
1488 Returns
1489 -------
1490 statistic : float or array
1491 t-statistic
1492 pvalue : float or array
1493 two-tailed p-value
1495 Notes
1496 -----
1497 For more details on `ttest_rel`, see `scipy.stats.ttest_rel`.
1499 """
1500 a, b, axis = _chk2_asarray(a, b, axis)
1501 if len(a) != len(b):
1502 raise ValueError('unequal length arrays')
1504 if a.size == 0 or b.size == 0:
1505 return Ttest_relResult(np.nan, np.nan)
1507 n = a.count(axis)
1508 df = ma.asanyarray(n-1.0)
1509 d = (a-b).astype('d')
1510 dm = d.mean(axis)
1511 v = d.var(axis=axis, ddof=1)
1512 denom = ma.sqrt(v / n)
1513 with np.errstate(divide='ignore', invalid='ignore'):
1514 t = dm / denom
1516 t, prob = scipy.stats._stats_py._ttest_finish(df, t, alternative)
1517 return Ttest_relResult(t, prob)
1520MannwhitneyuResult = namedtuple('MannwhitneyuResult', ('statistic',
1521 'pvalue'))
1524def mannwhitneyu(x,y, use_continuity=True):
1525 """
1526 Computes the Mann-Whitney statistic
1528 Missing values in `x` and/or `y` are discarded.
1530 Parameters
1531 ----------
1532 x : sequence
1533 Input
1534 y : sequence
1535 Input
1536 use_continuity : {True, False}, optional
1537 Whether a continuity correction (1/2.) should be taken into account.
1539 Returns
1540 -------
1541 statistic : float
1542 The minimum of the Mann-Whitney statistics
1543 pvalue : float
1544 Approximate two-sided p-value assuming a normal distribution.
1546 """
1547 x = ma.asarray(x).compressed().view(ndarray)
1548 y = ma.asarray(y).compressed().view(ndarray)
1549 ranks = rankdata(np.concatenate([x,y]))
1550 (nx, ny) = (len(x), len(y))
1551 nt = nx + ny
1552 U = ranks[:nx].sum() - nx*(nx+1)/2.
1553 U = max(U, nx*ny - U)
1554 u = nx*ny - U
1556 mu = (nx*ny)/2.
1557 sigsq = (nt**3 - nt)/12.
1558 ties = count_tied_groups(ranks)
1559 sigsq -= sum(v*(k**3-k) for (k,v) in ties.items())/12.
1560 sigsq *= nx*ny/float(nt*(nt-1))
1562 if use_continuity:
1563 z = (U - 1/2. - mu) / ma.sqrt(sigsq)
1564 else:
1565 z = (U - mu) / ma.sqrt(sigsq)
1567 prob = special.erfc(abs(z)/np.sqrt(2))
1568 return MannwhitneyuResult(u, prob)
1571KruskalResult = namedtuple('KruskalResult', ('statistic', 'pvalue'))
1574def kruskal(*args):
1575 """
1576 Compute the Kruskal-Wallis H-test for independent samples
1578 Parameters
1579 ----------
1580 sample1, sample2, ... : array_like
1581 Two or more arrays with the sample measurements can be given as
1582 arguments.
1584 Returns
1585 -------
1586 statistic : float
1587 The Kruskal-Wallis H statistic, corrected for ties
1588 pvalue : float
1589 The p-value for the test using the assumption that H has a chi
1590 square distribution
1592 Notes
1593 -----
1594 For more details on `kruskal`, see `scipy.stats.kruskal`.
1596 Examples
1597 --------
1598 >>> from scipy.stats.mstats import kruskal
1600 Random samples from three different brands of batteries were tested
1601 to see how long the charge lasted. Results were as follows:
1603 >>> a = [6.3, 5.4, 5.7, 5.2, 5.0]
1604 >>> b = [6.9, 7.0, 6.1, 7.9]
1605 >>> c = [7.2, 6.9, 6.1, 6.5]
1607 Test the hypotesis that the distribution functions for all of the brands'
1608 durations are identical. Use 5% level of significance.
1610 >>> kruskal(a, b, c)
1611 KruskalResult(statistic=7.113812154696133, pvalue=0.028526948491942164)
1613 The null hypothesis is rejected at the 5% level of significance
1614 because the returned p-value is less than the critical value of 5%.
1616 """
1617 output = argstoarray(*args)
1618 ranks = ma.masked_equal(rankdata(output, use_missing=False), 0)
1619 sumrk = ranks.sum(-1)
1620 ngrp = ranks.count(-1)
1621 ntot = ranks.count()
1622 H = 12./(ntot*(ntot+1)) * (sumrk**2/ngrp).sum() - 3*(ntot+1)
1623 # Tie correction
1624 ties = count_tied_groups(ranks)
1625 T = 1. - sum(v*(k**3-k) for (k,v) in ties.items())/float(ntot**3-ntot)
1626 if T == 0:
1627 raise ValueError('All numbers are identical in kruskal')
1629 H /= T
1630 df = len(output) - 1
1631 prob = distributions.chi2.sf(H, df)
1632 return KruskalResult(H, prob)
1635kruskalwallis = kruskal
1638@_rename_parameter("mode", "method")
1639def ks_1samp(x, cdf, args=(), alternative="two-sided", method='auto'):
1640 """
1641 Computes the Kolmogorov-Smirnov test on one sample of masked values.
1643 Missing values in `x` are discarded.
1645 Parameters
1646 ----------
1647 x : array_like
1648 a 1-D array of observations of random variables.
1649 cdf : str or callable
1650 If a string, it should be the name of a distribution in `scipy.stats`.
1651 If a callable, that callable is used to calculate the cdf.
1652 args : tuple, sequence, optional
1653 Distribution parameters, used if `cdf` is a string.
1654 alternative : {'two-sided', 'less', 'greater'}, optional
1655 Indicates the alternative hypothesis. Default is 'two-sided'.
1656 method : {'auto', 'exact', 'asymp'}, optional
1657 Defines the method used for calculating the p-value.
1658 The following options are available (default is 'auto'):
1660 * 'auto' : use 'exact' for small size arrays, 'asymp' for large
1661 * 'exact' : use approximation to exact distribution of test statistic
1662 * 'asymp' : use asymptotic distribution of test statistic
1664 Returns
1665 -------
1666 d : float
1667 Value of the Kolmogorov Smirnov test
1668 p : float
1669 Corresponding p-value.
1671 """
1672 alternative = {'t': 'two-sided', 'g': 'greater', 'l': 'less'}.get(
1673 alternative.lower()[0], alternative)
1674 return scipy.stats._stats_py.ks_1samp(
1675 x, cdf, args=args, alternative=alternative, method=method)
1678@_rename_parameter("mode", "method")
1679def ks_2samp(data1, data2, alternative="two-sided", method='auto'):
1680 """
1681 Computes the Kolmogorov-Smirnov test on two samples.
1683 Missing values in `x` and/or `y` are discarded.
1685 Parameters
1686 ----------
1687 data1 : array_like
1688 First data set
1689 data2 : array_like
1690 Second data set
1691 alternative : {'two-sided', 'less', 'greater'}, optional
1692 Indicates the alternative hypothesis. Default is 'two-sided'.
1693 method : {'auto', 'exact', 'asymp'}, optional
1694 Defines the method used for calculating the p-value.
1695 The following options are available (default is 'auto'):
1697 * 'auto' : use 'exact' for small size arrays, 'asymp' for large
1698 * 'exact' : use approximation to exact distribution of test statistic
1699 * 'asymp' : use asymptotic distribution of test statistic
1701 Returns
1702 -------
1703 d : float
1704 Value of the Kolmogorov Smirnov test
1705 p : float
1706 Corresponding p-value.
1708 """
1709 # Ideally this would be accomplished by
1710 # ks_2samp = scipy.stats._stats_py.ks_2samp
1711 # but the circular dependencies between _mstats_basic and stats prevent that.
1712 alternative = {'t': 'two-sided', 'g': 'greater', 'l': 'less'}.get(
1713 alternative.lower()[0], alternative)
1714 return scipy.stats._stats_py.ks_2samp(data1, data2,
1715 alternative=alternative,
1716 method=method)
1719ks_twosamp = ks_2samp
1722@_rename_parameter("mode", "method")
1723def kstest(data1, data2, args=(), alternative='two-sided', method='auto'):
1724 """
1726 Parameters
1727 ----------
1728 data1 : array_like
1729 data2 : str, callable or array_like
1730 args : tuple, sequence, optional
1731 Distribution parameters, used if `data1` or `data2` are strings.
1732 alternative : str, as documented in stats.kstest
1733 method : str, as documented in stats.kstest
1735 Returns
1736 -------
1737 tuple of (K-S statistic, probability)
1739 """
1740 return scipy.stats._stats_py.kstest(data1, data2, args,
1741 alternative=alternative, method=method)
1744def trima(a, limits=None, inclusive=(True,True)):
1745 """
1746 Trims an array by masking the data outside some given limits.
1748 Returns a masked version of the input array.
1750 Parameters
1751 ----------
1752 a : array_like
1753 Input array.
1754 limits : {None, tuple}, optional
1755 Tuple of (lower limit, upper limit) in absolute values.
1756 Values of the input array lower (greater) than the lower (upper) limit
1757 will be masked. A limit is None indicates an open interval.
1758 inclusive : (bool, bool) tuple, optional
1759 Tuple of (lower flag, upper flag), indicating whether values exactly
1760 equal to the lower (upper) limit are allowed.
1762 Examples
1763 --------
1764 >>> from scipy.stats.mstats import trima
1765 >>> import numpy as np
1767 >>> a = np.arange(10)
1769 The interval is left-closed and right-open, i.e., `[2, 8)`.
1770 Trim the array by keeping only values in the interval.
1772 >>> trima(a, limits=(2, 8), inclusive=(True, False))
1773 masked_array(data=[--, --, 2, 3, 4, 5, 6, 7, --, --],
1774 mask=[ True, True, False, False, False, False, False, False,
1775 True, True],
1776 fill_value=999999)
1778 """
1779 a = ma.asarray(a)
1780 a.unshare_mask()
1781 if (limits is None) or (limits == (None, None)):
1782 return a
1784 (lower_lim, upper_lim) = limits
1785 (lower_in, upper_in) = inclusive
1786 condition = False
1787 if lower_lim is not None:
1788 if lower_in:
1789 condition |= (a < lower_lim)
1790 else:
1791 condition |= (a <= lower_lim)
1793 if upper_lim is not None:
1794 if upper_in:
1795 condition |= (a > upper_lim)
1796 else:
1797 condition |= (a >= upper_lim)
1799 a[condition.filled(True)] = masked
1800 return a
1803def trimr(a, limits=None, inclusive=(True, True), axis=None):
1804 """
1805 Trims an array by masking some proportion of the data on each end.
1806 Returns a masked version of the input array.
1808 Parameters
1809 ----------
1810 a : sequence
1811 Input array.
1812 limits : {None, tuple}, optional
1813 Tuple of the percentages to cut on each side of the array, with respect
1814 to the number of unmasked data, as floats between 0. and 1.
1815 Noting n the number of unmasked data before trimming, the
1816 (n*limits[0])th smallest data and the (n*limits[1])th largest data are
1817 masked, and the total number of unmasked data after trimming is
1818 n*(1.-sum(limits)). The value of one limit can be set to None to
1819 indicate an open interval.
1820 inclusive : {(True,True) tuple}, optional
1821 Tuple of flags indicating whether the number of data being masked on
1822 the left (right) end should be truncated (True) or rounded (False) to
1823 integers.
1824 axis : {None,int}, optional
1825 Axis along which to trim. If None, the whole array is trimmed, but its
1826 shape is maintained.
1828 """
1829 def _trimr1D(a, low_limit, up_limit, low_inclusive, up_inclusive):
1830 n = a.count()
1831 idx = a.argsort()
1832 if low_limit:
1833 if low_inclusive:
1834 lowidx = int(low_limit*n)
1835 else:
1836 lowidx = int(np.round(low_limit*n))
1837 a[idx[:lowidx]] = masked
1838 if up_limit is not None:
1839 if up_inclusive:
1840 upidx = n - int(n*up_limit)
1841 else:
1842 upidx = n - int(np.round(n*up_limit))
1843 a[idx[upidx:]] = masked
1844 return a
1846 a = ma.asarray(a)
1847 a.unshare_mask()
1848 if limits is None:
1849 return a
1851 # Check the limits
1852 (lolim, uplim) = limits
1853 errmsg = "The proportion to cut from the %s should be between 0. and 1."
1854 if lolim is not None:
1855 if lolim > 1. or lolim < 0:
1856 raise ValueError(errmsg % 'beginning' + "(got %s)" % lolim)
1857 if uplim is not None:
1858 if uplim > 1. or uplim < 0:
1859 raise ValueError(errmsg % 'end' + "(got %s)" % uplim)
1861 (loinc, upinc) = inclusive
1863 if axis is None:
1864 shp = a.shape
1865 return _trimr1D(a.ravel(),lolim,uplim,loinc,upinc).reshape(shp)
1866 else:
1867 return ma.apply_along_axis(_trimr1D, axis, a, lolim,uplim,loinc,upinc)
1870trimdoc = """
1871 Parameters
1872 ----------
1873 a : sequence
1874 Input array
1875 limits : {None, tuple}, optional
1876 If `relative` is False, tuple (lower limit, upper limit) in absolute values.
1877 Values of the input array lower (greater) than the lower (upper) limit are
1878 masked.
1880 If `relative` is True, tuple (lower percentage, upper percentage) to cut
1881 on each side of the array, with respect to the number of unmasked data.
1883 Noting n the number of unmasked data before trimming, the (n*limits[0])th
1884 smallest data and the (n*limits[1])th largest data are masked, and the
1885 total number of unmasked data after trimming is n*(1.-sum(limits))
1886 In each case, the value of one limit can be set to None to indicate an
1887 open interval.
1889 If limits is None, no trimming is performed
1890 inclusive : {(bool, bool) tuple}, optional
1891 If `relative` is False, tuple indicating whether values exactly equal
1892 to the absolute limits are allowed.
1893 If `relative` is True, tuple indicating whether the number of data
1894 being masked on each side should be rounded (True) or truncated
1895 (False).
1896 relative : bool, optional
1897 Whether to consider the limits as absolute values (False) or proportions
1898 to cut (True).
1899 axis : int, optional
1900 Axis along which to trim.
1901"""
1904def trim(a, limits=None, inclusive=(True,True), relative=False, axis=None):
1905 """
1906 Trims an array by masking the data outside some given limits.
1908 Returns a masked version of the input array.
1910 %s
1912 Examples
1913 --------
1914 >>> from scipy.stats.mstats import trim
1915 >>> z = [ 1, 2, 3, 4, 5, 6, 7, 8, 9,10]
1916 >>> print(trim(z,(3,8)))
1917 [-- -- 3 4 5 6 7 8 -- --]
1918 >>> print(trim(z,(0.1,0.2),relative=True))
1919 [-- 2 3 4 5 6 7 8 -- --]
1921 """
1922 if relative:
1923 return trimr(a, limits=limits, inclusive=inclusive, axis=axis)
1924 else:
1925 return trima(a, limits=limits, inclusive=inclusive)
1928if trim.__doc__:
1929 trim.__doc__ = trim.__doc__ % trimdoc
1932def trimboth(data, proportiontocut=0.2, inclusive=(True,True), axis=None):
1933 """
1934 Trims the smallest and largest data values.
1936 Trims the `data` by masking the ``int(proportiontocut * n)`` smallest and
1937 ``int(proportiontocut * n)`` largest values of data along the given axis,
1938 where n is the number of unmasked values before trimming.
1940 Parameters
1941 ----------
1942 data : ndarray
1943 Data to trim.
1944 proportiontocut : float, optional
1945 Percentage of trimming (as a float between 0 and 1).
1946 If n is the number of unmasked values before trimming, the number of
1947 values after trimming is ``(1 - 2*proportiontocut) * n``.
1948 Default is 0.2.
1949 inclusive : {(bool, bool) tuple}, optional
1950 Tuple indicating whether the number of data being masked on each side
1951 should be rounded (True) or truncated (False).
1952 axis : int, optional
1953 Axis along which to perform the trimming.
1954 If None, the input array is first flattened.
1956 """
1957 return trimr(data, limits=(proportiontocut,proportiontocut),
1958 inclusive=inclusive, axis=axis)
1961def trimtail(data, proportiontocut=0.2, tail='left', inclusive=(True,True),
1962 axis=None):
1963 """
1964 Trims the data by masking values from one tail.
1966 Parameters
1967 ----------
1968 data : array_like
1969 Data to trim.
1970 proportiontocut : float, optional
1971 Percentage of trimming. If n is the number of unmasked values
1972 before trimming, the number of values after trimming is
1973 ``(1 - proportiontocut) * n``. Default is 0.2.
1974 tail : {'left','right'}, optional
1975 If 'left' the `proportiontocut` lowest values will be masked.
1976 If 'right' the `proportiontocut` highest values will be masked.
1977 Default is 'left'.
1978 inclusive : {(bool, bool) tuple}, optional
1979 Tuple indicating whether the number of data being masked on each side
1980 should be rounded (True) or truncated (False). Default is
1981 (True, True).
1982 axis : int, optional
1983 Axis along which to perform the trimming.
1984 If None, the input array is first flattened. Default is None.
1986 Returns
1987 -------
1988 trimtail : ndarray
1989 Returned array of same shape as `data` with masked tail values.
1991 """
1992 tail = str(tail).lower()[0]
1993 if tail == 'l':
1994 limits = (proportiontocut,None)
1995 elif tail == 'r':
1996 limits = (None, proportiontocut)
1997 else:
1998 raise TypeError("The tail argument should be in ('left','right')")
2000 return trimr(data, limits=limits, axis=axis, inclusive=inclusive)
2003trim1 = trimtail
2006def trimmed_mean(a, limits=(0.1,0.1), inclusive=(1,1), relative=True,
2007 axis=None):
2008 """Returns the trimmed mean of the data along the given axis.
2010 %s
2012 """
2013 if (not isinstance(limits,tuple)) and isinstance(limits,float):
2014 limits = (limits, limits)
2015 if relative:
2016 return trimr(a,limits=limits,inclusive=inclusive,axis=axis).mean(axis=axis)
2017 else:
2018 return trima(a,limits=limits,inclusive=inclusive).mean(axis=axis)
2021if trimmed_mean.__doc__:
2022 trimmed_mean.__doc__ = trimmed_mean.__doc__ % trimdoc
2025def trimmed_var(a, limits=(0.1,0.1), inclusive=(1,1), relative=True,
2026 axis=None, ddof=0):
2027 """Returns the trimmed variance of the data along the given axis.
2029 %s
2030 ddof : {0,integer}, optional
2031 Means Delta Degrees of Freedom. The denominator used during computations
2032 is (n-ddof). DDOF=0 corresponds to a biased estimate, DDOF=1 to an un-
2033 biased estimate of the variance.
2035 """
2036 if (not isinstance(limits,tuple)) and isinstance(limits,float):
2037 limits = (limits, limits)
2038 if relative:
2039 out = trimr(a,limits=limits, inclusive=inclusive,axis=axis)
2040 else:
2041 out = trima(a,limits=limits,inclusive=inclusive)
2043 return out.var(axis=axis, ddof=ddof)
2046if trimmed_var.__doc__:
2047 trimmed_var.__doc__ = trimmed_var.__doc__ % trimdoc
2050def trimmed_std(a, limits=(0.1,0.1), inclusive=(1,1), relative=True,
2051 axis=None, ddof=0):
2052 """Returns the trimmed standard deviation of the data along the given axis.
2054 %s
2055 ddof : {0,integer}, optional
2056 Means Delta Degrees of Freedom. The denominator used during computations
2057 is (n-ddof). DDOF=0 corresponds to a biased estimate, DDOF=1 to an un-
2058 biased estimate of the variance.
2060 """
2061 if (not isinstance(limits,tuple)) and isinstance(limits,float):
2062 limits = (limits, limits)
2063 if relative:
2064 out = trimr(a,limits=limits,inclusive=inclusive,axis=axis)
2065 else:
2066 out = trima(a,limits=limits,inclusive=inclusive)
2067 return out.std(axis=axis,ddof=ddof)
2070if trimmed_std.__doc__:
2071 trimmed_std.__doc__ = trimmed_std.__doc__ % trimdoc
2074def trimmed_stde(a, limits=(0.1,0.1), inclusive=(1,1), axis=None):
2075 """
2076 Returns the standard error of the trimmed mean along the given axis.
2078 Parameters
2079 ----------
2080 a : sequence
2081 Input array
2082 limits : {(0.1,0.1), tuple of float}, optional
2083 tuple (lower percentage, upper percentage) to cut on each side of the
2084 array, with respect to the number of unmasked data.
2086 If n is the number of unmasked data before trimming, the values
2087 smaller than ``n * limits[0]`` and the values larger than
2088 ``n * `limits[1]`` are masked, and the total number of unmasked
2089 data after trimming is ``n * (1.-sum(limits))``. In each case,
2090 the value of one limit can be set to None to indicate an open interval.
2091 If `limits` is None, no trimming is performed.
2092 inclusive : {(bool, bool) tuple} optional
2093 Tuple indicating whether the number of data being masked on each side
2094 should be rounded (True) or truncated (False).
2095 axis : int, optional
2096 Axis along which to trim.
2098 Returns
2099 -------
2100 trimmed_stde : scalar or ndarray
2102 """
2103 def _trimmed_stde_1D(a, low_limit, up_limit, low_inclusive, up_inclusive):
2104 "Returns the standard error of the trimmed mean for a 1D input data."
2105 n = a.count()
2106 idx = a.argsort()
2107 if low_limit:
2108 if low_inclusive:
2109 lowidx = int(low_limit*n)
2110 else:
2111 lowidx = np.round(low_limit*n)
2112 a[idx[:lowidx]] = masked
2113 if up_limit is not None:
2114 if up_inclusive:
2115 upidx = n - int(n*up_limit)
2116 else:
2117 upidx = n - np.round(n*up_limit)
2118 a[idx[upidx:]] = masked
2119 a[idx[:lowidx]] = a[idx[lowidx]]
2120 a[idx[upidx:]] = a[idx[upidx-1]]
2121 winstd = a.std(ddof=1)
2122 return winstd / ((1-low_limit-up_limit)*np.sqrt(len(a)))
2124 a = ma.array(a, copy=True, subok=True)
2125 a.unshare_mask()
2126 if limits is None:
2127 return a.std(axis=axis,ddof=1)/ma.sqrt(a.count(axis))
2128 if (not isinstance(limits,tuple)) and isinstance(limits,float):
2129 limits = (limits, limits)
2131 # Check the limits
2132 (lolim, uplim) = limits
2133 errmsg = "The proportion to cut from the %s should be between 0. and 1."
2134 if lolim is not None:
2135 if lolim > 1. or lolim < 0:
2136 raise ValueError(errmsg % 'beginning' + "(got %s)" % lolim)
2137 if uplim is not None:
2138 if uplim > 1. or uplim < 0:
2139 raise ValueError(errmsg % 'end' + "(got %s)" % uplim)
2141 (loinc, upinc) = inclusive
2142 if (axis is None):
2143 return _trimmed_stde_1D(a.ravel(),lolim,uplim,loinc,upinc)
2144 else:
2145 if a.ndim > 2:
2146 raise ValueError("Array 'a' must be at most two dimensional, "
2147 "but got a.ndim = %d" % a.ndim)
2148 return ma.apply_along_axis(_trimmed_stde_1D, axis, a,
2149 lolim,uplim,loinc,upinc)
2152def _mask_to_limits(a, limits, inclusive):
2153 """Mask an array for values outside of given limits.
2155 This is primarily a utility function.
2157 Parameters
2158 ----------
2159 a : array
2160 limits : (float or None, float or None)
2161 A tuple consisting of the (lower limit, upper limit). Values in the
2162 input array less than the lower limit or greater than the upper limit
2163 will be masked out. None implies no limit.
2164 inclusive : (bool, bool)
2165 A tuple consisting of the (lower flag, upper flag). These flags
2166 determine whether values exactly equal to lower or upper are allowed.
2168 Returns
2169 -------
2170 A MaskedArray.
2172 Raises
2173 ------
2174 A ValueError if there are no values within the given limits.
2175 """
2176 lower_limit, upper_limit = limits
2177 lower_include, upper_include = inclusive
2178 am = ma.MaskedArray(a)
2179 if lower_limit is not None:
2180 if lower_include:
2181 am = ma.masked_less(am, lower_limit)
2182 else:
2183 am = ma.masked_less_equal(am, lower_limit)
2185 if upper_limit is not None:
2186 if upper_include:
2187 am = ma.masked_greater(am, upper_limit)
2188 else:
2189 am = ma.masked_greater_equal(am, upper_limit)
2191 if am.count() == 0:
2192 raise ValueError("No array values within given limits")
2194 return am
2197def tmean(a, limits=None, inclusive=(True, True), axis=None):
2198 """
2199 Compute the trimmed mean.
2201 Parameters
2202 ----------
2203 a : array_like
2204 Array of values.
2205 limits : None or (lower limit, upper limit), optional
2206 Values in the input array less than the lower limit or greater than the
2207 upper limit will be ignored. When limits is None (default), then all
2208 values are used. Either of the limit values in the tuple can also be
2209 None representing a half-open interval.
2210 inclusive : (bool, bool), optional
2211 A tuple consisting of the (lower flag, upper flag). These flags
2212 determine whether values exactly equal to the lower or upper limits
2213 are included. The default value is (True, True).
2214 axis : int or None, optional
2215 Axis along which to operate. If None, compute over the
2216 whole array. Default is None.
2218 Returns
2219 -------
2220 tmean : float
2222 Notes
2223 -----
2224 For more details on `tmean`, see `scipy.stats.tmean`.
2226 Examples
2227 --------
2228 >>> import numpy as np
2229 >>> from scipy.stats import mstats
2230 >>> a = np.array([[6, 8, 3, 0],
2231 ... [3, 9, 1, 2],
2232 ... [8, 7, 8, 2],
2233 ... [5, 6, 0, 2],
2234 ... [4, 5, 5, 2]])
2235 ...
2236 ...
2237 >>> mstats.tmean(a, (2,5))
2238 3.3
2239 >>> mstats.tmean(a, (2,5), axis=0)
2240 masked_array(data=[4.0, 5.0, 4.0, 2.0],
2241 mask=[False, False, False, False],
2242 fill_value=1e+20)
2244 """
2245 return trima(a, limits=limits, inclusive=inclusive).mean(axis=axis)
2248def tvar(a, limits=None, inclusive=(True, True), axis=0, ddof=1):
2249 """
2250 Compute the trimmed variance
2252 This function computes the sample variance of an array of values,
2253 while ignoring values which are outside of given `limits`.
2255 Parameters
2256 ----------
2257 a : array_like
2258 Array of values.
2259 limits : None or (lower limit, upper limit), optional
2260 Values in the input array less than the lower limit or greater than the
2261 upper limit will be ignored. When limits is None, then all values are
2262 used. Either of the limit values in the tuple can also be None
2263 representing a half-open interval. The default value is None.
2264 inclusive : (bool, bool), optional
2265 A tuple consisting of the (lower flag, upper flag). These flags
2266 determine whether values exactly equal to the lower or upper limits
2267 are included. The default value is (True, True).
2268 axis : int or None, optional
2269 Axis along which to operate. If None, compute over the
2270 whole array. Default is zero.
2271 ddof : int, optional
2272 Delta degrees of freedom. Default is 1.
2274 Returns
2275 -------
2276 tvar : float
2277 Trimmed variance.
2279 Notes
2280 -----
2281 For more details on `tvar`, see `scipy.stats.tvar`.
2283 """
2284 a = a.astype(float).ravel()
2285 if limits is None:
2286 n = (~a.mask).sum() # todo: better way to do that?
2287 return np.ma.var(a) * n/(n-1.)
2288 am = _mask_to_limits(a, limits=limits, inclusive=inclusive)
2290 return np.ma.var(am, axis=axis, ddof=ddof)
2293def tmin(a, lowerlimit=None, axis=0, inclusive=True):
2294 """
2295 Compute the trimmed minimum
2297 Parameters
2298 ----------
2299 a : array_like
2300 array of values
2301 lowerlimit : None or float, optional
2302 Values in the input array less than the given limit will be ignored.
2303 When lowerlimit is None, then all values are used. The default value
2304 is None.
2305 axis : int or None, optional
2306 Axis along which to operate. Default is 0. If None, compute over the
2307 whole array `a`.
2308 inclusive : {True, False}, optional
2309 This flag determines whether values exactly equal to the lower limit
2310 are included. The default value is True.
2312 Returns
2313 -------
2314 tmin : float, int or ndarray
2316 Notes
2317 -----
2318 For more details on `tmin`, see `scipy.stats.tmin`.
2320 Examples
2321 --------
2322 >>> import numpy as np
2323 >>> from scipy.stats import mstats
2324 >>> a = np.array([[6, 8, 3, 0],
2325 ... [3, 2, 1, 2],
2326 ... [8, 1, 8, 2],
2327 ... [5, 3, 0, 2],
2328 ... [4, 7, 5, 2]])
2329 ...
2330 >>> mstats.tmin(a, 5)
2331 masked_array(data=[5, 7, 5, --],
2332 mask=[False, False, False, True],
2333 fill_value=999999)
2335 """
2336 a, axis = _chk_asarray(a, axis)
2337 am = trima(a, (lowerlimit, None), (inclusive, False))
2338 return ma.minimum.reduce(am, axis)
2341def tmax(a, upperlimit=None, axis=0, inclusive=True):
2342 """
2343 Compute the trimmed maximum
2345 This function computes the maximum value of an array along a given axis,
2346 while ignoring values larger than a specified upper limit.
2348 Parameters
2349 ----------
2350 a : array_like
2351 array of values
2352 upperlimit : None or float, optional
2353 Values in the input array greater than the given limit will be ignored.
2354 When upperlimit is None, then all values are used. The default value
2355 is None.
2356 axis : int or None, optional
2357 Axis along which to operate. Default is 0. If None, compute over the
2358 whole array `a`.
2359 inclusive : {True, False}, optional
2360 This flag determines whether values exactly equal to the upper limit
2361 are included. The default value is True.
2363 Returns
2364 -------
2365 tmax : float, int or ndarray
2367 Notes
2368 -----
2369 For more details on `tmax`, see `scipy.stats.tmax`.
2371 Examples
2372 --------
2373 >>> import numpy as np
2374 >>> from scipy.stats import mstats
2375 >>> a = np.array([[6, 8, 3, 0],
2376 ... [3, 9, 1, 2],
2377 ... [8, 7, 8, 2],
2378 ... [5, 6, 0, 2],
2379 ... [4, 5, 5, 2]])
2380 ...
2381 ...
2382 >>> mstats.tmax(a, 4)
2383 masked_array(data=[4, --, 3, 2],
2384 mask=[False, True, False, False],
2385 fill_value=999999)
2387 """
2388 a, axis = _chk_asarray(a, axis)
2389 am = trima(a, (None, upperlimit), (False, inclusive))
2390 return ma.maximum.reduce(am, axis)
2393def tsem(a, limits=None, inclusive=(True, True), axis=0, ddof=1):
2394 """
2395 Compute the trimmed standard error of the mean.
2397 This function finds the standard error of the mean for given
2398 values, ignoring values outside the given `limits`.
2400 Parameters
2401 ----------
2402 a : array_like
2403 array of values
2404 limits : None or (lower limit, upper limit), optional
2405 Values in the input array less than the lower limit or greater than the
2406 upper limit will be ignored. When limits is None, then all values are
2407 used. Either of the limit values in the tuple can also be None
2408 representing a half-open interval. The default value is None.
2409 inclusive : (bool, bool), optional
2410 A tuple consisting of the (lower flag, upper flag). These flags
2411 determine whether values exactly equal to the lower or upper limits
2412 are included. The default value is (True, True).
2413 axis : int or None, optional
2414 Axis along which to operate. If None, compute over the
2415 whole array. Default is zero.
2416 ddof : int, optional
2417 Delta degrees of freedom. Default is 1.
2419 Returns
2420 -------
2421 tsem : float
2423 Notes
2424 -----
2425 For more details on `tsem`, see `scipy.stats.tsem`.
2427 """
2428 a = ma.asarray(a).ravel()
2429 if limits is None:
2430 n = float(a.count())
2431 return a.std(axis=axis, ddof=ddof)/ma.sqrt(n)
2433 am = trima(a.ravel(), limits, inclusive)
2434 sd = np.sqrt(am.var(axis=axis, ddof=ddof))
2435 return sd / np.sqrt(am.count())
2438def winsorize(a, limits=None, inclusive=(True, True), inplace=False,
2439 axis=None, nan_policy='propagate'):
2440 """Returns a Winsorized version of the input array.
2442 The (limits[0])th lowest values are set to the (limits[0])th percentile,
2443 and the (limits[1])th highest values are set to the (1 - limits[1])th
2444 percentile.
2445 Masked values are skipped.
2448 Parameters
2449 ----------
2450 a : sequence
2451 Input array.
2452 limits : {None, tuple of float}, optional
2453 Tuple of the percentages to cut on each side of the array, with respect
2454 to the number of unmasked data, as floats between 0. and 1.
2455 Noting n the number of unmasked data before trimming, the
2456 (n*limits[0])th smallest data and the (n*limits[1])th largest data are
2457 masked, and the total number of unmasked data after trimming
2458 is n*(1.-sum(limits)) The value of one limit can be set to None to
2459 indicate an open interval.
2460 inclusive : {(True, True) tuple}, optional
2461 Tuple indicating whether the number of data being masked on each side
2462 should be truncated (True) or rounded (False).
2463 inplace : {False, True}, optional
2464 Whether to winsorize in place (True) or to use a copy (False)
2465 axis : {None, int}, optional
2466 Axis along which to trim. If None, the whole array is trimmed, but its
2467 shape is maintained.
2468 nan_policy : {'propagate', 'raise', 'omit'}, optional
2469 Defines how to handle when input contains nan.
2470 The following options are available (default is 'propagate'):
2472 * 'propagate': allows nan values and may overwrite or propagate them
2473 * 'raise': throws an error
2474 * 'omit': performs the calculations ignoring nan values
2476 Notes
2477 -----
2478 This function is applied to reduce the effect of possibly spurious outliers
2479 by limiting the extreme values.
2481 Examples
2482 --------
2483 >>> import numpy as np
2484 >>> from scipy.stats.mstats import winsorize
2486 A shuffled array contains integers from 1 to 10.
2488 >>> a = np.array([10, 4, 9, 8, 5, 3, 7, 2, 1, 6])
2490 The 10% of the lowest value (i.e., `1`) and the 20% of the highest
2491 values (i.e., `9` and `10`) are replaced.
2493 >>> winsorize(a, limits=[0.1, 0.2])
2494 masked_array(data=[8, 4, 8, 8, 5, 3, 7, 2, 2, 6],
2495 mask=False,
2496 fill_value=999999)
2498 """
2499 def _winsorize1D(a, low_limit, up_limit, low_include, up_include,
2500 contains_nan, nan_policy):
2501 n = a.count()
2502 idx = a.argsort()
2503 if contains_nan:
2504 nan_count = np.count_nonzero(np.isnan(a))
2505 if low_limit:
2506 if low_include:
2507 lowidx = int(low_limit * n)
2508 else:
2509 lowidx = np.round(low_limit * n).astype(int)
2510 if contains_nan and nan_policy == 'omit':
2511 lowidx = min(lowidx, n-nan_count-1)
2512 a[idx[:lowidx]] = a[idx[lowidx]]
2513 if up_limit is not None:
2514 if up_include:
2515 upidx = n - int(n * up_limit)
2516 else:
2517 upidx = n - np.round(n * up_limit).astype(int)
2518 if contains_nan and nan_policy == 'omit':
2519 a[idx[upidx:-nan_count]] = a[idx[upidx - 1]]
2520 else:
2521 a[idx[upidx:]] = a[idx[upidx - 1]]
2522 return a
2524 contains_nan, nan_policy = _contains_nan(a, nan_policy)
2525 # We are going to modify a: better make a copy
2526 a = ma.array(a, copy=np.logical_not(inplace))
2528 if limits is None:
2529 return a
2530 if (not isinstance(limits, tuple)) and isinstance(limits, float):
2531 limits = (limits, limits)
2533 # Check the limits
2534 (lolim, uplim) = limits
2535 errmsg = "The proportion to cut from the %s should be between 0. and 1."
2536 if lolim is not None:
2537 if lolim > 1. or lolim < 0:
2538 raise ValueError(errmsg % 'beginning' + "(got %s)" % lolim)
2539 if uplim is not None:
2540 if uplim > 1. or uplim < 0:
2541 raise ValueError(errmsg % 'end' + "(got %s)" % uplim)
2543 (loinc, upinc) = inclusive
2545 if axis is None:
2546 shp = a.shape
2547 return _winsorize1D(a.ravel(), lolim, uplim, loinc, upinc,
2548 contains_nan, nan_policy).reshape(shp)
2549 else:
2550 return ma.apply_along_axis(_winsorize1D, axis, a, lolim, uplim, loinc,
2551 upinc, contains_nan, nan_policy)
2554def moment(a, moment=1, axis=0):
2555 """
2556 Calculates the nth moment about the mean for a sample.
2558 Parameters
2559 ----------
2560 a : array_like
2561 data
2562 moment : int, optional
2563 order of central moment that is returned
2564 axis : int or None, optional
2565 Axis along which the central moment is computed. Default is 0.
2566 If None, compute over the whole array `a`.
2568 Returns
2569 -------
2570 n-th central moment : ndarray or float
2571 The appropriate moment along the given axis or over all values if axis
2572 is None. The denominator for the moment calculation is the number of
2573 observations, no degrees of freedom correction is done.
2575 Notes
2576 -----
2577 For more details about `moment`, see `scipy.stats.moment`.
2579 """
2580 a, axis = _chk_asarray(a, axis)
2581 if a.size == 0:
2582 moment_shape = list(a.shape)
2583 del moment_shape[axis]
2584 dtype = a.dtype.type if a.dtype.kind in 'fc' else np.float64
2585 # empty array, return nan(s) with shape matching `moment`
2586 out_shape = (moment_shape if np.isscalar(moment)
2587 else [len(moment)] + moment_shape)
2588 if len(out_shape) == 0:
2589 return dtype(np.nan)
2590 else:
2591 return ma.array(np.full(out_shape, np.nan, dtype=dtype))
2593 # for array_like moment input, return a value for each.
2594 if not np.isscalar(moment):
2595 mean = a.mean(axis, keepdims=True)
2596 mmnt = [_moment(a, i, axis, mean=mean) for i in moment]
2597 return ma.array(mmnt)
2598 else:
2599 return _moment(a, moment, axis)
2602# Moment with optional pre-computed mean, equal to a.mean(axis, keepdims=True)
2603def _moment(a, moment, axis, *, mean=None):
2604 if np.abs(moment - np.round(moment)) > 0:
2605 raise ValueError("All moment parameters must be integers")
2607 if moment == 0 or moment == 1:
2608 # By definition the zeroth moment about the mean is 1, and the first
2609 # moment is 0.
2610 shape = list(a.shape)
2611 del shape[axis]
2612 dtype = a.dtype.type if a.dtype.kind in 'fc' else np.float64
2614 if len(shape) == 0:
2615 return dtype(1.0 if moment == 0 else 0.0)
2616 else:
2617 return (ma.ones(shape, dtype=dtype) if moment == 0
2618 else ma.zeros(shape, dtype=dtype))
2619 else:
2620 # Exponentiation by squares: form exponent sequence
2621 n_list = [moment]
2622 current_n = moment
2623 while current_n > 2:
2624 if current_n % 2:
2625 current_n = (current_n-1)/2
2626 else:
2627 current_n /= 2
2628 n_list.append(current_n)
2630 # Starting point for exponentiation by squares
2631 mean = a.mean(axis, keepdims=True) if mean is None else mean
2632 a_zero_mean = a - mean
2633 if n_list[-1] == 1:
2634 s = a_zero_mean.copy()
2635 else:
2636 s = a_zero_mean**2
2638 # Perform multiplications
2639 for n in n_list[-2::-1]:
2640 s = s**2
2641 if n % 2:
2642 s *= a_zero_mean
2643 return s.mean(axis)
2646def variation(a, axis=0, ddof=0):
2647 """
2648 Compute the coefficient of variation.
2650 The coefficient of variation is the standard deviation divided by the
2651 mean. This function is equivalent to::
2653 np.std(x, axis=axis, ddof=ddof) / np.mean(x)
2655 The default for ``ddof`` is 0, but many definitions of the coefficient
2656 of variation use the square root of the unbiased sample variance
2657 for the sample standard deviation, which corresponds to ``ddof=1``.
2659 Parameters
2660 ----------
2661 a : array_like
2662 Input array.
2663 axis : int or None, optional
2664 Axis along which to calculate the coefficient of variation. Default
2665 is 0. If None, compute over the whole array `a`.
2666 ddof : int, optional
2667 Delta degrees of freedom. Default is 0.
2669 Returns
2670 -------
2671 variation : ndarray
2672 The calculated variation along the requested axis.
2674 Notes
2675 -----
2676 For more details about `variation`, see `scipy.stats.variation`.
2678 Examples
2679 --------
2680 >>> import numpy as np
2681 >>> from scipy.stats.mstats import variation
2682 >>> a = np.array([2,8,4])
2683 >>> variation(a)
2684 0.5345224838248487
2685 >>> b = np.array([2,8,3,4])
2686 >>> c = np.ma.masked_array(b, mask=[0,0,1,0])
2687 >>> variation(c)
2688 0.5345224838248487
2690 In the example above, it can be seen that this works the same as
2691 `scipy.stats.variation` except 'stats.mstats.variation' ignores masked
2692 array elements.
2694 """
2695 a, axis = _chk_asarray(a, axis)
2696 return a.std(axis, ddof=ddof)/a.mean(axis)
2699def skew(a, axis=0, bias=True):
2700 """
2701 Computes the skewness of a data set.
2703 Parameters
2704 ----------
2705 a : ndarray
2706 data
2707 axis : int or None, optional
2708 Axis along which skewness is calculated. Default is 0.
2709 If None, compute over the whole array `a`.
2710 bias : bool, optional
2711 If False, then the calculations are corrected for statistical bias.
2713 Returns
2714 -------
2715 skewness : ndarray
2716 The skewness of values along an axis, returning 0 where all values are
2717 equal.
2719 Notes
2720 -----
2721 For more details about `skew`, see `scipy.stats.skew`.
2723 """
2724 a, axis = _chk_asarray(a,axis)
2725 mean = a.mean(axis, keepdims=True)
2726 m2 = _moment(a, 2, axis, mean=mean)
2727 m3 = _moment(a, 3, axis, mean=mean)
2728 zero = (m2 <= (np.finfo(m2.dtype).resolution * mean.squeeze(axis))**2)
2729 with np.errstate(all='ignore'):
2730 vals = ma.where(zero, 0, m3 / m2**1.5)
2732 if not bias and zero is not ma.masked and m2 is not ma.masked:
2733 n = a.count(axis)
2734 can_correct = ~zero & (n > 2)
2735 if can_correct.any():
2736 n = np.extract(can_correct, n)
2737 m2 = np.extract(can_correct, m2)
2738 m3 = np.extract(can_correct, m3)
2739 nval = ma.sqrt((n-1.0)*n)/(n-2.0)*m3/m2**1.5
2740 np.place(vals, can_correct, nval)
2741 return vals
2744def kurtosis(a, axis=0, fisher=True, bias=True):
2745 """
2746 Computes the kurtosis (Fisher or Pearson) of a dataset.
2748 Kurtosis is the fourth central moment divided by the square of the
2749 variance. If Fisher's definition is used, then 3.0 is subtracted from
2750 the result to give 0.0 for a normal distribution.
2752 If bias is False then the kurtosis is calculated using k statistics to
2753 eliminate bias coming from biased moment estimators
2755 Use `kurtosistest` to see if result is close enough to normal.
2757 Parameters
2758 ----------
2759 a : array
2760 data for which the kurtosis is calculated
2761 axis : int or None, optional
2762 Axis along which the kurtosis is calculated. Default is 0.
2763 If None, compute over the whole array `a`.
2764 fisher : bool, optional
2765 If True, Fisher's definition is used (normal ==> 0.0). If False,
2766 Pearson's definition is used (normal ==> 3.0).
2767 bias : bool, optional
2768 If False, then the calculations are corrected for statistical bias.
2770 Returns
2771 -------
2772 kurtosis : array
2773 The kurtosis of values along an axis. If all values are equal,
2774 return -3 for Fisher's definition and 0 for Pearson's definition.
2776 Notes
2777 -----
2778 For more details about `kurtosis`, see `scipy.stats.kurtosis`.
2780 """
2781 a, axis = _chk_asarray(a, axis)
2782 mean = a.mean(axis, keepdims=True)
2783 m2 = _moment(a, 2, axis, mean=mean)
2784 m4 = _moment(a, 4, axis, mean=mean)
2785 zero = (m2 <= (np.finfo(m2.dtype).resolution * mean.squeeze(axis))**2)
2786 with np.errstate(all='ignore'):
2787 vals = ma.where(zero, 0, m4 / m2**2.0)
2789 if not bias and zero is not ma.masked and m2 is not ma.masked:
2790 n = a.count(axis)
2791 can_correct = ~zero & (n > 3)
2792 if can_correct.any():
2793 n = np.extract(can_correct, n)
2794 m2 = np.extract(can_correct, m2)
2795 m4 = np.extract(can_correct, m4)
2796 nval = 1.0/(n-2)/(n-3)*((n*n-1.0)*m4/m2**2.0-3*(n-1)**2.0)
2797 np.place(vals, can_correct, nval+3.0)
2798 if fisher:
2799 return vals - 3
2800 else:
2801 return vals
2804DescribeResult = namedtuple('DescribeResult', ('nobs', 'minmax', 'mean',
2805 'variance', 'skewness',
2806 'kurtosis'))
2809def describe(a, axis=0, ddof=0, bias=True):
2810 """
2811 Computes several descriptive statistics of the passed array.
2813 Parameters
2814 ----------
2815 a : array_like
2816 Data array
2817 axis : int or None, optional
2818 Axis along which to calculate statistics. Default 0. If None,
2819 compute over the whole array `a`.
2820 ddof : int, optional
2821 degree of freedom (default 0); note that default ddof is different
2822 from the same routine in stats.describe
2823 bias : bool, optional
2824 If False, then the skewness and kurtosis calculations are corrected for
2825 statistical bias.
2827 Returns
2828 -------
2829 nobs : int
2830 (size of the data (discarding missing values)
2832 minmax : (int, int)
2833 min, max
2835 mean : float
2836 arithmetic mean
2838 variance : float
2839 unbiased variance
2841 skewness : float
2842 biased skewness
2844 kurtosis : float
2845 biased kurtosis
2847 Examples
2848 --------
2849 >>> import numpy as np
2850 >>> from scipy.stats.mstats import describe
2851 >>> ma = np.ma.array(range(6), mask=[0, 0, 0, 1, 1, 1])
2852 >>> describe(ma)
2853 DescribeResult(nobs=3, minmax=(masked_array(data=0,
2854 mask=False,
2855 fill_value=999999), masked_array(data=2,
2856 mask=False,
2857 fill_value=999999)), mean=1.0, variance=0.6666666666666666,
2858 skewness=masked_array(data=0., mask=False, fill_value=1e+20),
2859 kurtosis=-1.5)
2861 """
2862 a, axis = _chk_asarray(a, axis)
2863 n = a.count(axis)
2864 mm = (ma.minimum.reduce(a, axis=axis), ma.maximum.reduce(a, axis=axis))
2865 m = a.mean(axis)
2866 v = a.var(axis, ddof=ddof)
2867 sk = skew(a, axis, bias=bias)
2868 kurt = kurtosis(a, axis, bias=bias)
2870 return DescribeResult(n, mm, m, v, sk, kurt)
2873def stde_median(data, axis=None):
2874 """Returns the McKean-Schrader estimate of the standard error of the sample
2875 median along the given axis. masked values are discarded.
2877 Parameters
2878 ----------
2879 data : ndarray
2880 Data to trim.
2881 axis : {None,int}, optional
2882 Axis along which to perform the trimming.
2883 If None, the input array is first flattened.
2885 """
2886 def _stdemed_1D(data):
2887 data = np.sort(data.compressed())
2888 n = len(data)
2889 z = 2.5758293035489004
2890 k = int(np.round((n+1)/2. - z * np.sqrt(n/4.),0))
2891 return ((data[n-k] - data[k-1])/(2.*z))
2893 data = ma.array(data, copy=False, subok=True)
2894 if (axis is None):
2895 return _stdemed_1D(data)
2896 else:
2897 if data.ndim > 2:
2898 raise ValueError("Array 'data' must be at most two dimensional, "
2899 "but got data.ndim = %d" % data.ndim)
2900 return ma.apply_along_axis(_stdemed_1D, axis, data)
2903SkewtestResult = namedtuple('SkewtestResult', ('statistic', 'pvalue'))
2906def skewtest(a, axis=0, alternative='two-sided'):
2907 """
2908 Tests whether the skew is different from the normal distribution.
2910 Parameters
2911 ----------
2912 a : array_like
2913 The data to be tested
2914 axis : int or None, optional
2915 Axis along which statistics are calculated. Default is 0.
2916 If None, compute over the whole array `a`.
2917 alternative : {'two-sided', 'less', 'greater'}, optional
2918 Defines the alternative hypothesis. Default is 'two-sided'.
2919 The following options are available:
2921 * 'two-sided': the skewness of the distribution underlying the sample
2922 is different from that of the normal distribution (i.e. 0)
2923 * 'less': the skewness of the distribution underlying the sample
2924 is less than that of the normal distribution
2925 * 'greater': the skewness of the distribution underlying the sample
2926 is greater than that of the normal distribution
2928 .. versionadded:: 1.7.0
2930 Returns
2931 -------
2932 statistic : array_like
2933 The computed z-score for this test.
2934 pvalue : array_like
2935 A p-value for the hypothesis test
2937 Notes
2938 -----
2939 For more details about `skewtest`, see `scipy.stats.skewtest`.
2941 """
2942 a, axis = _chk_asarray(a, axis)
2943 if axis is None:
2944 a = a.ravel()
2945 axis = 0
2946 b2 = skew(a,axis)
2947 n = a.count(axis)
2948 if np.min(n) < 8:
2949 raise ValueError(
2950 "skewtest is not valid with less than 8 samples; %i samples"
2951 " were given." % np.min(n))
2953 y = b2 * ma.sqrt(((n+1)*(n+3)) / (6.0*(n-2)))
2954 beta2 = (3.0*(n*n+27*n-70)*(n+1)*(n+3)) / ((n-2.0)*(n+5)*(n+7)*(n+9))
2955 W2 = -1 + ma.sqrt(2*(beta2-1))
2956 delta = 1/ma.sqrt(0.5*ma.log(W2))
2957 alpha = ma.sqrt(2.0/(W2-1))
2958 y = ma.where(y == 0, 1, y)
2959 Z = delta*ma.log(y/alpha + ma.sqrt((y/alpha)**2+1))
2961 return SkewtestResult(*scipy.stats._stats_py._normtest_finish(Z, alternative))
2964KurtosistestResult = namedtuple('KurtosistestResult', ('statistic', 'pvalue'))
2967def kurtosistest(a, axis=0, alternative='two-sided'):
2968 """
2969 Tests whether a dataset has normal kurtosis
2971 Parameters
2972 ----------
2973 a : array_like
2974 array of the sample data
2975 axis : int or None, optional
2976 Axis along which to compute test. Default is 0. If None,
2977 compute over the whole array `a`.
2978 alternative : {'two-sided', 'less', 'greater'}, optional
2979 Defines the alternative hypothesis.
2980 The following options are available (default is 'two-sided'):
2982 * 'two-sided': the kurtosis of the distribution underlying the sample
2983 is different from that of the normal distribution
2984 * 'less': the kurtosis of the distribution underlying the sample
2985 is less than that of the normal distribution
2986 * 'greater': the kurtosis of the distribution underlying the sample
2987 is greater than that of the normal distribution
2989 .. versionadded:: 1.7.0
2991 Returns
2992 -------
2993 statistic : array_like
2994 The computed z-score for this test.
2995 pvalue : array_like
2996 The p-value for the hypothesis test
2998 Notes
2999 -----
3000 For more details about `kurtosistest`, see `scipy.stats.kurtosistest`.
3002 """
3003 a, axis = _chk_asarray(a, axis)
3004 n = a.count(axis=axis)
3005 if np.min(n) < 5:
3006 raise ValueError(
3007 "kurtosistest requires at least 5 observations; %i observations"
3008 " were given." % np.min(n))
3009 if np.min(n) < 20:
3010 warnings.warn(
3011 "kurtosistest only valid for n>=20 ... continuing anyway, n=%i" %
3012 np.min(n))
3014 b2 = kurtosis(a, axis, fisher=False)
3015 E = 3.0*(n-1) / (n+1)
3016 varb2 = 24.0*n*(n-2.)*(n-3) / ((n+1)*(n+1.)*(n+3)*(n+5))
3017 x = (b2-E)/ma.sqrt(varb2)
3018 sqrtbeta1 = 6.0*(n*n-5*n+2)/((n+7)*(n+9)) * np.sqrt((6.0*(n+3)*(n+5)) /
3019 (n*(n-2)*(n-3)))
3020 A = 6.0 + 8.0/sqrtbeta1 * (2.0/sqrtbeta1 + np.sqrt(1+4.0/(sqrtbeta1**2)))
3021 term1 = 1 - 2./(9.0*A)
3022 denom = 1 + x*ma.sqrt(2/(A-4.0))
3023 if np.ma.isMaskedArray(denom):
3024 # For multi-dimensional array input
3025 denom[denom == 0.0] = masked
3026 elif denom == 0.0:
3027 denom = masked
3029 term2 = np.ma.where(denom > 0, ma.power((1-2.0/A)/denom, 1/3.0),
3030 -ma.power(-(1-2.0/A)/denom, 1/3.0))
3031 Z = (term1 - term2) / np.sqrt(2/(9.0*A))
3033 return KurtosistestResult(
3034 *scipy.stats._stats_py._normtest_finish(Z, alternative)
3035 )
3038NormaltestResult = namedtuple('NormaltestResult', ('statistic', 'pvalue'))
3041def normaltest(a, axis=0):
3042 """
3043 Tests whether a sample differs from a normal distribution.
3045 Parameters
3046 ----------
3047 a : array_like
3048 The array containing the data to be tested.
3049 axis : int or None, optional
3050 Axis along which to compute test. Default is 0. If None,
3051 compute over the whole array `a`.
3053 Returns
3054 -------
3055 statistic : float or array
3056 ``s^2 + k^2``, where ``s`` is the z-score returned by `skewtest` and
3057 ``k`` is the z-score returned by `kurtosistest`.
3058 pvalue : float or array
3059 A 2-sided chi squared probability for the hypothesis test.
3061 Notes
3062 -----
3063 For more details about `normaltest`, see `scipy.stats.normaltest`.
3065 """
3066 a, axis = _chk_asarray(a, axis)
3067 s, _ = skewtest(a, axis)
3068 k, _ = kurtosistest(a, axis)
3069 k2 = s*s + k*k
3071 return NormaltestResult(k2, distributions.chi2.sf(k2, 2))
3074def mquantiles(a, prob=list([.25,.5,.75]), alphap=.4, betap=.4, axis=None,
3075 limit=()):
3076 """
3077 Computes empirical quantiles for a data array.
3079 Samples quantile are defined by ``Q(p) = (1-gamma)*x[j] + gamma*x[j+1]``,
3080 where ``x[j]`` is the j-th order statistic, and gamma is a function of
3081 ``j = floor(n*p + m)``, ``m = alphap + p*(1 - alphap - betap)`` and
3082 ``g = n*p + m - j``.
3084 Reinterpreting the above equations to compare to **R** lead to the
3085 equation: ``p(k) = (k - alphap)/(n + 1 - alphap - betap)``
3087 Typical values of (alphap,betap) are:
3088 - (0,1) : ``p(k) = k/n`` : linear interpolation of cdf
3089 (**R** type 4)
3090 - (.5,.5) : ``p(k) = (k - 1/2.)/n`` : piecewise linear function
3091 (**R** type 5)
3092 - (0,0) : ``p(k) = k/(n+1)`` :
3093 (**R** type 6)
3094 - (1,1) : ``p(k) = (k-1)/(n-1)``: p(k) = mode[F(x[k])].
3095 (**R** type 7, **R** default)
3096 - (1/3,1/3): ``p(k) = (k-1/3)/(n+1/3)``: Then p(k) ~ median[F(x[k])].
3097 The resulting quantile estimates are approximately median-unbiased
3098 regardless of the distribution of x.
3099 (**R** type 8)
3100 - (3/8,3/8): ``p(k) = (k-3/8)/(n+1/4)``: Blom.
3101 The resulting quantile estimates are approximately unbiased
3102 if x is normally distributed
3103 (**R** type 9)
3104 - (.4,.4) : approximately quantile unbiased (Cunnane)
3105 - (.35,.35): APL, used with PWM
3107 Parameters
3108 ----------
3109 a : array_like
3110 Input data, as a sequence or array of dimension at most 2.
3111 prob : array_like, optional
3112 List of quantiles to compute.
3113 alphap : float, optional
3114 Plotting positions parameter, default is 0.4.
3115 betap : float, optional
3116 Plotting positions parameter, default is 0.4.
3117 axis : int, optional
3118 Axis along which to perform the trimming.
3119 If None (default), the input array is first flattened.
3120 limit : tuple, optional
3121 Tuple of (lower, upper) values.
3122 Values of `a` outside this open interval are ignored.
3124 Returns
3125 -------
3126 mquantiles : MaskedArray
3127 An array containing the calculated quantiles.
3129 Notes
3130 -----
3131 This formulation is very similar to **R** except the calculation of
3132 ``m`` from ``alphap`` and ``betap``, where in **R** ``m`` is defined
3133 with each type.
3135 References
3136 ----------
3137 .. [1] *R* statistical software: https://www.r-project.org/
3138 .. [2] *R* ``quantile`` function:
3139 http://stat.ethz.ch/R-manual/R-devel/library/stats/html/quantile.html
3141 Examples
3142 --------
3143 >>> import numpy as np
3144 >>> from scipy.stats.mstats import mquantiles
3145 >>> a = np.array([6., 47., 49., 15., 42., 41., 7., 39., 43., 40., 36.])
3146 >>> mquantiles(a)
3147 array([ 19.2, 40. , 42.8])
3149 Using a 2D array, specifying axis and limit.
3151 >>> data = np.array([[ 6., 7., 1.],
3152 ... [ 47., 15., 2.],
3153 ... [ 49., 36., 3.],
3154 ... [ 15., 39., 4.],
3155 ... [ 42., 40., -999.],
3156 ... [ 41., 41., -999.],
3157 ... [ 7., -999., -999.],
3158 ... [ 39., -999., -999.],
3159 ... [ 43., -999., -999.],
3160 ... [ 40., -999., -999.],
3161 ... [ 36., -999., -999.]])
3162 >>> print(mquantiles(data, axis=0, limit=(0, 50)))
3163 [[19.2 14.6 1.45]
3164 [40. 37.5 2.5 ]
3165 [42.8 40.05 3.55]]
3167 >>> data[:, 2] = -999.
3168 >>> print(mquantiles(data, axis=0, limit=(0, 50)))
3169 [[19.200000000000003 14.6 --]
3170 [40.0 37.5 --]
3171 [42.800000000000004 40.05 --]]
3173 """
3174 def _quantiles1D(data,m,p):
3175 x = np.sort(data.compressed())
3176 n = len(x)
3177 if n == 0:
3178 return ma.array(np.empty(len(p), dtype=float), mask=True)
3179 elif n == 1:
3180 return ma.array(np.resize(x, p.shape), mask=nomask)
3181 aleph = (n*p + m)
3182 k = np.floor(aleph.clip(1, n-1)).astype(int)
3183 gamma = (aleph-k).clip(0,1)
3184 return (1.-gamma)*x[(k-1).tolist()] + gamma*x[k.tolist()]
3186 data = ma.array(a, copy=False)
3187 if data.ndim > 2:
3188 raise TypeError("Array should be 2D at most !")
3190 if limit:
3191 condition = (limit[0] < data) & (data < limit[1])
3192 data[~condition.filled(True)] = masked
3194 p = np.array(prob, copy=False, ndmin=1)
3195 m = alphap + p*(1.-alphap-betap)
3196 # Computes quantiles along axis (or globally)
3197 if (axis is None):
3198 return _quantiles1D(data, m, p)
3200 return ma.apply_along_axis(_quantiles1D, axis, data, m, p)
3203def scoreatpercentile(data, per, limit=(), alphap=.4, betap=.4):
3204 """Calculate the score at the given 'per' percentile of the
3205 sequence a. For example, the score at per=50 is the median.
3207 This function is a shortcut to mquantile
3209 """
3210 if (per < 0) or (per > 100.):
3211 raise ValueError("The percentile should be between 0. and 100. !"
3212 " (got %s)" % per)
3214 return mquantiles(data, prob=[per/100.], alphap=alphap, betap=betap,
3215 limit=limit, axis=0).squeeze()
3218def plotting_positions(data, alpha=0.4, beta=0.4):
3219 """
3220 Returns plotting positions (or empirical percentile points) for the data.
3222 Plotting positions are defined as ``(i-alpha)/(n+1-alpha-beta)``, where:
3223 - i is the rank order statistics
3224 - n is the number of unmasked values along the given axis
3225 - `alpha` and `beta` are two parameters.
3227 Typical values for `alpha` and `beta` are:
3228 - (0,1) : ``p(k) = k/n``, linear interpolation of cdf (R, type 4)
3229 - (.5,.5) : ``p(k) = (k-1/2.)/n``, piecewise linear function
3230 (R, type 5)
3231 - (0,0) : ``p(k) = k/(n+1)``, Weibull (R type 6)
3232 - (1,1) : ``p(k) = (k-1)/(n-1)``, in this case,
3233 ``p(k) = mode[F(x[k])]``. That's R default (R type 7)
3234 - (1/3,1/3): ``p(k) = (k-1/3)/(n+1/3)``, then
3235 ``p(k) ~ median[F(x[k])]``.
3236 The resulting quantile estimates are approximately median-unbiased
3237 regardless of the distribution of x. (R type 8)
3238 - (3/8,3/8): ``p(k) = (k-3/8)/(n+1/4)``, Blom.
3239 The resulting quantile estimates are approximately unbiased
3240 if x is normally distributed (R type 9)
3241 - (.4,.4) : approximately quantile unbiased (Cunnane)
3242 - (.35,.35): APL, used with PWM
3243 - (.3175, .3175): used in scipy.stats.probplot
3245 Parameters
3246 ----------
3247 data : array_like
3248 Input data, as a sequence or array of dimension at most 2.
3249 alpha : float, optional
3250 Plotting positions parameter. Default is 0.4.
3251 beta : float, optional
3252 Plotting positions parameter. Default is 0.4.
3254 Returns
3255 -------
3256 positions : MaskedArray
3257 The calculated plotting positions.
3259 """
3260 data = ma.array(data, copy=False).reshape(1,-1)
3261 n = data.count()
3262 plpos = np.empty(data.size, dtype=float)
3263 plpos[n:] = 0
3264 plpos[data.argsort(axis=None)[:n]] = ((np.arange(1, n+1) - alpha) /
3265 (n + 1.0 - alpha - beta))
3266 return ma.array(plpos, mask=data._mask)
3269meppf = plotting_positions
3272def obrientransform(*args):
3273 """
3274 Computes a transform on input data (any number of columns). Used to
3275 test for homogeneity of variance prior to running one-way stats. Each
3276 array in ``*args`` is one level of a factor. If an `f_oneway()` run on
3277 the transformed data and found significant, variances are unequal. From
3278 Maxwell and Delaney, p.112.
3280 Returns: transformed data for use in an ANOVA
3281 """
3282 data = argstoarray(*args).T
3283 v = data.var(axis=0,ddof=1)
3284 m = data.mean(0)
3285 n = data.count(0).astype(float)
3286 # result = ((N-1.5)*N*(a-m)**2 - 0.5*v*(n-1))/((n-1)*(n-2))
3287 data -= m
3288 data **= 2
3289 data *= (n-1.5)*n
3290 data -= 0.5*v*(n-1)
3291 data /= (n-1.)*(n-2.)
3292 if not ma.allclose(v,data.mean(0)):
3293 raise ValueError("Lack of convergence in obrientransform.")
3295 return data
3298def sem(a, axis=0, ddof=1):
3299 """
3300 Calculates the standard error of the mean of the input array.
3302 Also sometimes called standard error of measurement.
3304 Parameters
3305 ----------
3306 a : array_like
3307 An array containing the values for which the standard error is
3308 returned.
3309 axis : int or None, optional
3310 If axis is None, ravel `a` first. If axis is an integer, this will be
3311 the axis over which to operate. Defaults to 0.
3312 ddof : int, optional
3313 Delta degrees-of-freedom. How many degrees of freedom to adjust
3314 for bias in limited samples relative to the population estimate
3315 of variance. Defaults to 1.
3317 Returns
3318 -------
3319 s : ndarray or float
3320 The standard error of the mean in the sample(s), along the input axis.
3322 Notes
3323 -----
3324 The default value for `ddof` changed in scipy 0.15.0 to be consistent with
3325 `scipy.stats.sem` as well as with the most common definition used (like in
3326 the R documentation).
3328 Examples
3329 --------
3330 Find standard error along the first axis:
3332 >>> import numpy as np
3333 >>> from scipy import stats
3334 >>> a = np.arange(20).reshape(5,4)
3335 >>> print(stats.mstats.sem(a))
3336 [2.8284271247461903 2.8284271247461903 2.8284271247461903
3337 2.8284271247461903]
3339 Find standard error across the whole array, using n degrees of freedom:
3341 >>> print(stats.mstats.sem(a, axis=None, ddof=0))
3342 1.2893796958227628
3344 """
3345 a, axis = _chk_asarray(a, axis)
3346 n = a.count(axis=axis)
3347 s = a.std(axis=axis, ddof=ddof) / ma.sqrt(n)
3348 return s
3351F_onewayResult = namedtuple('F_onewayResult', ('statistic', 'pvalue'))
3354def f_oneway(*args):
3355 """
3356 Performs a 1-way ANOVA, returning an F-value and probability given
3357 any number of groups. From Heiman, pp.394-7.
3359 Usage: ``f_oneway(*args)``, where ``*args`` is 2 or more arrays,
3360 one per treatment group.
3362 Returns
3363 -------
3364 statistic : float
3365 The computed F-value of the test.
3366 pvalue : float
3367 The associated p-value from the F-distribution.
3369 """
3370 # Construct a single array of arguments: each row is a group
3371 data = argstoarray(*args)
3372 ngroups = len(data)
3373 ntot = data.count()
3374 sstot = (data**2).sum() - (data.sum())**2/float(ntot)
3375 ssbg = (data.count(-1) * (data.mean(-1)-data.mean())**2).sum()
3376 sswg = sstot-ssbg
3377 dfbg = ngroups-1
3378 dfwg = ntot - ngroups
3379 msb = ssbg/float(dfbg)
3380 msw = sswg/float(dfwg)
3381 f = msb/msw
3382 prob = special.fdtrc(dfbg, dfwg, f) # equivalent to stats.f.sf
3384 return F_onewayResult(f, prob)
3387FriedmanchisquareResult = namedtuple('FriedmanchisquareResult',
3388 ('statistic', 'pvalue'))
3391def friedmanchisquare(*args):
3392 """Friedman Chi-Square is a non-parametric, one-way within-subjects ANOVA.
3393 This function calculates the Friedman Chi-square test for repeated measures
3394 and returns the result, along with the associated probability value.
3396 Each input is considered a given group. Ideally, the number of treatments
3397 among each group should be equal. If this is not the case, only the first
3398 n treatments are taken into account, where n is the number of treatments
3399 of the smallest group.
3400 If a group has some missing values, the corresponding treatments are masked
3401 in the other groups.
3402 The test statistic is corrected for ties.
3404 Masked values in one group are propagated to the other groups.
3406 Returns
3407 -------
3408 statistic : float
3409 the test statistic.
3410 pvalue : float
3411 the associated p-value.
3413 """
3414 data = argstoarray(*args).astype(float)
3415 k = len(data)
3416 if k < 3:
3417 raise ValueError("Less than 3 groups (%i): " % k +
3418 "the Friedman test is NOT appropriate.")
3420 ranked = ma.masked_values(rankdata(data, axis=0), 0)
3421 if ranked._mask is not nomask:
3422 ranked = ma.mask_cols(ranked)
3423 ranked = ranked.compressed().reshape(k,-1).view(ndarray)
3424 else:
3425 ranked = ranked._data
3426 (k,n) = ranked.shape
3427 # Ties correction
3428 repeats = [find_repeats(row) for row in ranked.T]
3429 ties = np.array([y for x, y in repeats if x.size > 0])
3430 tie_correction = 1 - (ties**3-ties).sum()/float(n*(k**3-k))
3432 ssbg = np.sum((ranked.sum(-1) - n*(k+1)/2.)**2)
3433 chisq = ssbg * 12./(n*k*(k+1)) * 1./tie_correction
3435 return FriedmanchisquareResult(chisq,
3436 distributions.chi2.sf(chisq, k-1))
3439BrunnerMunzelResult = namedtuple('BrunnerMunzelResult', ('statistic', 'pvalue'))
3442def brunnermunzel(x, y, alternative="two-sided", distribution="t"):
3443 """
3444 Computes the Brunner-Munzel test on samples x and y
3446 Missing values in `x` and/or `y` are discarded.
3448 Parameters
3449 ----------
3450 x, y : array_like
3451 Array of samples, should be one-dimensional.
3452 alternative : 'less', 'two-sided', or 'greater', optional
3453 Whether to get the p-value for the one-sided hypothesis ('less'
3454 or 'greater') or for the two-sided hypothesis ('two-sided').
3455 Defaults value is 'two-sided' .
3456 distribution : 't' or 'normal', optional
3457 Whether to get the p-value by t-distribution or by standard normal
3458 distribution.
3459 Defaults value is 't' .
3461 Returns
3462 -------
3463 statistic : float
3464 The Brunner-Munzer W statistic.
3465 pvalue : float
3466 p-value assuming an t distribution. One-sided or
3467 two-sided, depending on the choice of `alternative` and `distribution`.
3469 See Also
3470 --------
3471 mannwhitneyu : Mann-Whitney rank test on two samples.
3473 Notes
3474 -----
3475 For more details on `brunnermunzel`, see `scipy.stats.brunnermunzel`.
3477 """
3478 x = ma.asarray(x).compressed().view(ndarray)
3479 y = ma.asarray(y).compressed().view(ndarray)
3480 nx = len(x)
3481 ny = len(y)
3482 if nx == 0 or ny == 0:
3483 return BrunnerMunzelResult(np.nan, np.nan)
3484 rankc = rankdata(np.concatenate((x,y)))
3485 rankcx = rankc[0:nx]
3486 rankcy = rankc[nx:nx+ny]
3487 rankcx_mean = np.mean(rankcx)
3488 rankcy_mean = np.mean(rankcy)
3489 rankx = rankdata(x)
3490 ranky = rankdata(y)
3491 rankx_mean = np.mean(rankx)
3492 ranky_mean = np.mean(ranky)
3494 Sx = np.sum(np.power(rankcx - rankx - rankcx_mean + rankx_mean, 2.0))
3495 Sx /= nx - 1
3496 Sy = np.sum(np.power(rankcy - ranky - rankcy_mean + ranky_mean, 2.0))
3497 Sy /= ny - 1
3499 wbfn = nx * ny * (rankcy_mean - rankcx_mean)
3500 wbfn /= (nx + ny) * np.sqrt(nx * Sx + ny * Sy)
3502 if distribution == "t":
3503 df_numer = np.power(nx * Sx + ny * Sy, 2.0)
3504 df_denom = np.power(nx * Sx, 2.0) / (nx - 1)
3505 df_denom += np.power(ny * Sy, 2.0) / (ny - 1)
3506 df = df_numer / df_denom
3507 p = distributions.t.cdf(wbfn, df)
3508 elif distribution == "normal":
3509 p = distributions.norm.cdf(wbfn)
3510 else:
3511 raise ValueError(
3512 "distribution should be 't' or 'normal'")
3514 if alternative == "greater":
3515 pass
3516 elif alternative == "less":
3517 p = 1 - p
3518 elif alternative == "two-sided":
3519 p = 2 * np.min([p, 1-p])
3520 else:
3521 raise ValueError(
3522 "alternative should be 'less', 'greater' or 'two-sided'")
3524 return BrunnerMunzelResult(wbfn, p)