_bws_test.py 7.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179
  1. import numpy as np
  2. from functools import partial
  3. from scipy import stats
  4. from scipy._lib._array_api import xp_capabilities
  5. def _bws_input_validation(x, y, alternative, method):
  6. ''' Input validation and standardization for bws test'''
  7. x, y = np.atleast_1d(x, y)
  8. if x.ndim > 1 or y.ndim > 1:
  9. raise ValueError('`x` and `y` must be exactly one-dimensional.')
  10. if np.isnan(x).any() or np.isnan(y).any():
  11. raise ValueError('`x` and `y` must not contain NaNs.')
  12. if np.size(x) == 0 or np.size(y) == 0:
  13. raise ValueError('`x` and `y` must be of nonzero size.')
  14. z = stats.rankdata(np.concatenate((x, y)))
  15. x, y = z[:len(x)], z[len(x):]
  16. alternatives = {'two-sided', 'less', 'greater'}
  17. alternative = alternative.lower()
  18. if alternative not in alternatives:
  19. raise ValueError(f'`alternative` must be one of {alternatives}.')
  20. method = stats.PermutationMethod() if method is None else method
  21. if not isinstance(method, stats.PermutationMethod):
  22. raise ValueError('`method` must be an instance of '
  23. '`scipy.stats.PermutationMethod`')
  24. return x, y, alternative, method
  25. def _bws_statistic(x, y, alternative, axis):
  26. '''Compute the BWS test statistic for two independent samples'''
  27. # Public function currently does not accept `axis`, but `permutation_test`
  28. # uses `axis` to make vectorized call.
  29. Ri, Hj = np.sort(x, axis=axis), np.sort(y, axis=axis)
  30. n, m = Ri.shape[axis], Hj.shape[axis]
  31. i, j = np.arange(1, n+1), np.arange(1, m+1)
  32. Bx_num = Ri - (m + n)/n * i
  33. By_num = Hj - (m + n)/m * j
  34. if alternative == 'two-sided':
  35. Bx_num *= Bx_num
  36. By_num *= By_num
  37. else:
  38. Bx_num *= np.abs(Bx_num)
  39. By_num *= np.abs(By_num)
  40. Bx_den = i/(n+1) * (1 - i/(n+1)) * m*(m+n)/n
  41. By_den = j/(m+1) * (1 - j/(m+1)) * n*(m+n)/m
  42. Bx = 1/n * np.sum(Bx_num/Bx_den, axis=axis)
  43. By = 1/m * np.sum(By_num/By_den, axis=axis)
  44. B = (Bx + By) / 2 if alternative == 'two-sided' else (Bx - By) / 2
  45. return B
  46. @xp_capabilities(np_only=True)
  47. def bws_test(x, y, *, alternative="two-sided", method=None):
  48. r'''Perform the Baumgartner-Weiss-Schindler test on two independent samples.
  49. The Baumgartner-Weiss-Schindler (BWS) test is a nonparametric test of
  50. the null hypothesis that the distribution underlying sample `x`
  51. is the same as the distribution underlying sample `y`. Unlike
  52. the Kolmogorov-Smirnov, Wilcoxon, and Cramer-Von Mises tests,
  53. the BWS test weights the integral by the variance of the difference
  54. in cumulative distribution functions (CDFs), emphasizing the tails of the
  55. distributions, which increases the power of the test in many applications.
  56. Parameters
  57. ----------
  58. x, y : array-like
  59. 1-d arrays of samples.
  60. alternative : {'two-sided', 'less', 'greater'}, optional
  61. Defines the alternative hypothesis. Default is 'two-sided'.
  62. Let *F(u)* and *G(u)* be the cumulative distribution functions of the
  63. distributions underlying `x` and `y`, respectively. Then the following
  64. alternative hypotheses are available:
  65. * 'two-sided': the distributions are not equal, i.e. *F(u) ≠ G(u)* for
  66. at least one *u*.
  67. * 'less': the distribution underlying `x` is stochastically less than
  68. the distribution underlying `y`, i.e. *F(u) >= G(u)* for all *u*.
  69. * 'greater': the distribution underlying `x` is stochastically greater
  70. than the distribution underlying `y`, i.e. *F(u) <= G(u)* for all
  71. *u*.
  72. Under a more restrictive set of assumptions, the alternative hypotheses
  73. can be expressed in terms of the locations of the distributions;
  74. see [2] section 5.1.
  75. method : PermutationMethod, optional
  76. Configures the method used to compute the p-value. The default is
  77. the default `PermutationMethod` object.
  78. Returns
  79. -------
  80. res : PermutationTestResult
  81. An object with attributes:
  82. statistic : float
  83. The observed test statistic of the data.
  84. pvalue : float
  85. The p-value for the given alternative.
  86. null_distribution : ndarray
  87. The values of the test statistic generated under the null hypothesis.
  88. See also
  89. --------
  90. scipy.stats.wilcoxon, scipy.stats.mannwhitneyu, scipy.stats.ttest_ind
  91. Notes
  92. -----
  93. When ``alternative=='two-sided'``, the statistic is defined by the
  94. equations given in [1]_ Section 2. This statistic is not appropriate for
  95. one-sided alternatives; in that case, the statistic is the *negative* of
  96. that given by the equations in [1]_ Section 2. Consequently, when the
  97. distribution of the first sample is stochastically greater than that of the
  98. second sample, the statistic will tend to be positive.
  99. References
  100. ----------
  101. .. [1] Neuhäuser, M. (2005). Exact Tests Based on the
  102. Baumgartner-Weiss-Schindler Statistic: A Survey. Statistical Papers,
  103. 46(1), 1-29.
  104. .. [2] Fay, M. P., & Proschan, M. A. (2010). Wilcoxon-Mann-Whitney or t-test?
  105. On assumptions for hypothesis tests and multiple interpretations of
  106. decision rules. Statistics surveys, 4, 1.
  107. Examples
  108. --------
  109. We follow the example of table 3 in [1]_: Fourteen children were divided
  110. randomly into two groups. Their ranks at performing a specific tests are
  111. as follows.
  112. >>> import numpy as np
  113. >>> x = [1, 2, 3, 4, 6, 7, 8]
  114. >>> y = [5, 9, 10, 11, 12, 13, 14]
  115. We use the BWS test to assess whether there is a statistically significant
  116. difference between the two groups.
  117. The null hypothesis is that there is no difference in the distributions of
  118. performance between the two groups. We decide that a significance level of
  119. 1% is required to reject the null hypothesis in favor of the alternative
  120. that the distributions are different.
  121. Since the number of samples is very small, we can compare the observed test
  122. statistic against the *exact* distribution of the test statistic under the
  123. null hypothesis.
  124. >>> from scipy.stats import bws_test
  125. >>> res = bws_test(x, y)
  126. >>> print(res.statistic)
  127. 5.132167152575315
  128. This agrees with :math:`B = 5.132` reported in [1]_. The *p*-value produced
  129. by `bws_test` also agrees with :math:`p = 0.0029` reported in [1]_.
  130. >>> print(res.pvalue)
  131. 0.002913752913752914
  132. Because the p-value is below our threshold of 1%, we take this as evidence
  133. against the null hypothesis in favor of the alternative that there is a
  134. difference in performance between the two groups.
  135. '''
  136. x, y, alternative, method = _bws_input_validation(x, y, alternative,
  137. method)
  138. bws_statistic = partial(_bws_statistic, alternative=alternative)
  139. permutation_alternative = 'less' if alternative == 'less' else 'greater'
  140. res = stats.permutation_test((x, y), bws_statistic,
  141. alternative=permutation_alternative,
  142. **method._asdict())
  143. return res