import warnings import pytest import numpy as np from numpy.testing import assert_allclose, assert_equal from scipy._lib._util import rng_integers from scipy._lib._array_api import (is_numpy, make_xp_test_case, xp_default_dtype, xp_size, array_namespace, _xp_copy_to_numpy) from scipy._lib._array_api_no_0d import xp_assert_close, xp_assert_equal from scipy._lib import array_api_extra as xpx from scipy import stats, special from scipy.fft.tests.test_fftlog import skip_xp_backends from scipy.optimize import root from scipy.stats import bootstrap, monte_carlo_test, permutation_test, power import scipy.stats._resampling as _resampling @make_xp_test_case(bootstrap) class TestBootstrap: @skip_xp_backends('numpy', reason='NumPy does not raise') def test_bootstrap_iv_other(self, xp): message = f"When using array library {xp.__name__}" with pytest.raises(TypeError, match=message): bootstrap((xp.asarray([1, 2, 3]),), lambda x: xp.mean(x)) def test_bootstrap_iv(self, xp): sample = xp.asarray([1, 2, 3]) message = "`data` must contain at least one sample." with pytest.raises(ValueError, match=message): bootstrap(tuple(), xp.mean) message = "each sample in `data` must contain two or more observations..." with pytest.raises(ValueError, match=message): bootstrap((sample, xp.asarray([1])), xp.mean) message = ("When `paired is True`, all samples must have the same length ") with pytest.raises(ValueError, match=message): bootstrap((sample, xp.asarray([1, 2, 3, 4])), xp.mean, paired=True) message = "`vectorized` must be `True`, `False`, or `None`." with pytest.raises(ValueError, match=message): bootstrap(sample, xp.mean, vectorized='ekki') message = "`axis` must be an integer." with pytest.raises(ValueError, match=message): bootstrap((sample,), xp.mean, axis=1.5) message = "could not convert string to float" with pytest.raises(ValueError, match=message): bootstrap((sample,), xp.mean, confidence_level='ni') message = "`n_resamples` must be a non-negative integer." with pytest.raises(ValueError, match=message): bootstrap((sample,), xp.mean, n_resamples=-1000) message = "`n_resamples` must be a non-negative integer." with pytest.raises(ValueError, match=message): bootstrap((sample,), xp.mean, n_resamples=1000.5) message = "`batch` must be a positive integer or None." with pytest.raises(ValueError, match=message): bootstrap((sample,), xp.mean, batch=-1000) message = "`batch` must be a positive integer or None." with pytest.raises(ValueError, match=message): bootstrap((sample,), xp.mean, batch=1000.5) message = "`method` must be in" with pytest.raises(ValueError, match=message): bootstrap((sample,), xp.mean, method='ekki') message = "`bootstrap_result` must have attribute `bootstrap_distribution'" with pytest.raises(ValueError, match=message): bootstrap((sample,), xp.mean, bootstrap_result=10) message = "Either `bootstrap_result.bootstrap_distribution.size`" with pytest.raises(ValueError, match=message): bootstrap((sample,), xp.mean, n_resamples=0) message = "SeedSequence expects int or sequence of ints" with pytest.raises(TypeError, match=message): bootstrap((sample,), xp.mean, rng='herring') @pytest.mark.parametrize("method", ['basic', 'percentile', 'BCa']) @pytest.mark.parametrize("axis", [0, 1, 2]) def test_bootstrap_batch(self, method, axis, xp): # for one-sample statistics, batch size shouldn't affect the result rng = np.random.RandomState(0) x = rng.rand(10, 11, 12) # SPEC-007 leave one call with random_state to ensure it still works res1 = bootstrap((xp.asarray(x),), xp.mean, batch=None, method=method, random_state=0, axis=axis, n_resamples=100) rng = np.random.RandomState(0) res2 = bootstrap((xp.asarray(x),), xp.mean, batch=10, method=method, axis=axis, n_resamples=100, random_state=rng) xp_assert_equal(res2.confidence_interval.low, res1.confidence_interval.low) xp_assert_equal(res2.confidence_interval.high, res1.confidence_interval.high) xp_assert_equal(res2.standard_error, res1.standard_error) @pytest.mark.parametrize("method", ['basic', 'percentile', 'BCa']) def test_bootstrap_paired(self, method, xp): # test that `paired` works as expected rng = np.random.RandomState(0) n = 100 x = xp.asarray(rng.rand(n)) y = xp.asarray(rng.rand(n)) def my_statistic(x, y, axis=-1): return xp.mean((x-y)**2, axis=axis) def my_paired_statistic(i, axis=-1): a = x[i] b = y[i] res = my_statistic(a, b) return res i = xp.arange(x.shape[0]) res1 = bootstrap((i,), my_paired_statistic, rng=0) res2 = bootstrap((x, y), my_statistic, paired=True, rng=0) xp_assert_close(res1.confidence_interval.low, res2.confidence_interval.low) xp_assert_close(res1.confidence_interval.high, res2.confidence_interval.high) xp_assert_close(res1.standard_error, res2.standard_error) @pytest.mark.parametrize("method", ['basic', 'percentile', 'BCa']) @pytest.mark.parametrize("axis", [0, 1, 2]) @pytest.mark.parametrize("paired", [True, False]) def test_bootstrap_vectorized(self, method, axis, paired, xp): # test that paired is vectorized as expected: when samples are tiled, # CI and standard_error of each axis-slice is the same as those of the # original 1d sample rng = np.random.RandomState(0) def my_statistic(x, y, z, axis=-1): return xp.mean(x, axis=axis) + xp.mean(y, axis=axis) + xp.mean(z, axis=axis) shape = 10, 11, 12 n_samples = shape[axis] x = rng.rand(n_samples) y = rng.rand(n_samples) z = rng.rand(n_samples) x, y, z = xp.asarray(x), xp.asarray(y), xp.asarray(z) res1 = bootstrap((x, y, z), my_statistic, paired=paired, method=method, rng=0, axis=0, n_resamples=100) assert (res1.bootstrap_distribution.shape == res1.standard_error.shape + (100,)) reshape = [1, 1, 1] reshape[axis] = n_samples reshape = tuple(reshape) x = xp.broadcast_to(xp.reshape(x, reshape), shape) y = xp.broadcast_to(xp.reshape(y, reshape), shape) z = xp.broadcast_to(xp.reshape(z, reshape), shape) res2 = bootstrap((x, y, z), my_statistic, paired=paired, method=method, rng=0, axis=axis, n_resamples=100) result_shape = list(shape) result_shape.pop(axis) ref_ci_low = xp.broadcast_to(res2.confidence_interval.low, result_shape) ref_ci_high = xp.broadcast_to(res2.confidence_interval.high, result_shape) ref_ci_standard_error = xp.broadcast_to(res2.standard_error, result_shape) xp_assert_close(res2.confidence_interval.low, ref_ci_low) xp_assert_close(res2.confidence_interval.high, ref_ci_high) xp_assert_close(res2.standard_error, ref_ci_standard_error) @pytest.mark.slow @pytest.mark.xfail_on_32bit("MemoryError with BCa observed in CI") @pytest.mark.parametrize("method", ['basic', 'percentile', 'BCa']) def test_bootstrap_against_theory(self, method, xp): # based on https://www.statology.org/confidence-intervals-python/ rng = np.random.default_rng(2442101192988600726) data = stats.norm.rvs(loc=5, scale=2, size=5000, random_state=rng) alpha = 0.95 dist = stats.t(df=len(data)-1, loc=np.mean(data), scale=stats.sem(data)) ref_low, ref_high = dist.interval(confidence=alpha) ref_se = dist.std() config = dict(data=(xp.asarray(data),), statistic=xp.mean, n_resamples=5000, method=method, rng=rng) res = bootstrap(**config, confidence_level=alpha) xp_assert_close(res.confidence_interval.low, xp.asarray(ref_low), rtol=5e-4) xp_assert_close(res.confidence_interval.high, xp.asarray(ref_high), rtol=5e-4) xp_assert_close(res.standard_error, xp.asarray(ref_se), atol=3e-4) config.update(dict(n_resamples=0, bootstrap_result=res)) res = bootstrap(**config, confidence_level=alpha, alternative='less') xp_assert_close(res.confidence_interval.high, xp.asarray(dist.ppf(alpha)), rtol=5e-4) config.update(dict(n_resamples=0, bootstrap_result=res)) res = bootstrap(**config, confidence_level=alpha, alternative='greater') xp_assert_close(res.confidence_interval.low, xp.asarray(dist.ppf(1-alpha)), rtol=5e-4) tests_R = [("basic", 23.77, 79.12), ("percentile", 28.86, 84.21), ("BCa", 32.31, 91.43)] @pytest.mark.parametrize("method, ref_low, ref_high", tests_R) def test_bootstrap_against_R(self, method, ref_low, ref_high, xp): # Compare against R's "boot" library # library(boot) # stat <- function (x, a) { # mean(x[a]) # } # x <- c(10, 12, 12.5, 12.5, 13.9, 15, 21, 22, # 23, 34, 50, 81, 89, 121, 134, 213) # # Use a large value so we get a few significant digits for the CI. # n = 1000000 # bootresult = boot(x, stat, n) # result <- boot.ci(bootresult) # print(result) x = xp.asarray([10, 12, 12.5, 12.5, 13.9, 15, 21, 22, 23, 34, 50, 81, 89, 121, 134, 213]) res = bootstrap((x,), xp.mean, n_resamples=1000000, method=method, rng=0) xp_assert_close(res.confidence_interval.low, xp.asarray(ref_low), rtol=0.005) xp_assert_close(res.confidence_interval.high, xp.asarray(ref_high), rtol=0.005) def test_multisample_BCa_against_R(self, xp): # Because bootstrap is stochastic, it's tricky to test against reference # behavior. Here, we show that SciPy's BCa CI matches R wboot's BCa CI # much more closely than the other SciPy CIs do. # arbitrary skewed data x = xp.asarray([0.75859206, 0.5910282, -0.4419409, -0.36654601, 0.34955357, -1.38835871, 0.76735821]) y = xp.asarray([1.41186073, 0.49775975, 0.08275588, 0.24086388, 0.03567057, 0.52024419, 0.31966611, 1.32067634]) # a multi-sample statistic for which the BCa CI tends to be different # from the other CIs def statistic(x, y, axis): s1 = stats.skew(x, axis=axis) s2 = stats.skew(y, axis=axis) return s1 - s2 # compute confidence intervals using each method rng = np.random.default_rng(468865032284792692) res_basic = stats.bootstrap((x, y), statistic, method='basic', batch=100, rng=rng) res_percent = stats.bootstrap((x, y), statistic, method='percentile', batch=100, rng=rng) res_bca = stats.bootstrap((x, y), statistic, method='bca', batch=100, rng=rng) # compute midpoints so we can compare just one number for each mid_basic = xp.mean(xp.stack(res_basic.confidence_interval)) mid_percent = xp.mean(xp.stack(res_percent.confidence_interval)) mid_bca = xp.mean(xp.stack(res_bca.confidence_interval)) # reference for BCA CI computed using R wboot package: # library(wBoot) # library(moments) # x = c(0.75859206, 0.5910282, -0.4419409, -0.36654601, # 0.34955357, -1.38835871, 0.76735821) # y = c(1.41186073, 0.49775975, 0.08275588, 0.24086388, # 0.03567057, 0.52024419, 0.31966611, 1.32067634) # twoskew <- function(x1, y1) {skewness(x1) - skewness(y1)} # boot.two.bca(x, y, skewness, conf.level = 0.95, # R = 9999, stacked = FALSE) mid_wboot = -1.5519 # compute percent difference relative to wboot BCA method diff_basic = (mid_basic - mid_wboot)/abs(mid_wboot) diff_percent = (mid_percent - mid_wboot)/abs(mid_wboot) diff_bca = (mid_bca - mid_wboot)/abs(mid_wboot) # SciPy's BCa CI midpoint is much closer than that of the other methods assert diff_basic < -0.15 assert diff_percent > 0.15 assert abs(diff_bca) < 0.03 def test_BCa_acceleration_against_reference(self, xp): # Compare the (deterministic) acceleration parameter for a multi-sample # problem against a reference value. The example is from [1], but Efron's # value seems inaccurate. Straightforward code for computing the # reference acceleration (0.011008228344026734) is available at: # https://github.com/scipy/scipy/pull/16455#issuecomment-1193400981 y = xp.asarray([10., 27., 31., 40., 46., 50., 52., 104., 146.]) z = xp.asarray([16., 23., 38., 94., 99., 141., 197.]) def statistic(z, y, axis=0): return xp.mean(z, axis=axis) - xp.mean(y, axis=axis) data = [z, y] res = stats.bootstrap(data, statistic) axis = -1 alpha = 0.95 theta_hat_b = res.bootstrap_distribution batch = 100 _, _, a_hat = _resampling._bca_interval(data, statistic, axis, alpha, theta_hat_b, batch, xp) xp_assert_close(a_hat, xp.asarray(0.011008228344026734)) tests_against_itself_1samp = {"basic": 1789, "percentile": 1790, "BCa": 1789} @pytest.mark.slow @pytest.mark.parametrize("method, expected", tests_against_itself_1samp.items()) def test_bootstrap_against_itself_1samp(self, method, expected, xp): # The expected values in this test were generated using bootstrap # to check for unintended changes in behavior. The test also makes sure # that bootstrap works with multi-sample statistics and that the # `axis` argument works as expected / function is vectorized. rng = np.random.default_rng(9123847) n = 100 # size of sample n_resamples = 999 # number of bootstrap resamples used to form each CI confidence_level = 0.9 # The true mean is 5 dist = stats.norm(loc=5, scale=1) stat_true = float(dist.mean()) # Do the same thing 2000 times. (The code is fully vectorized.) n_replications = 2000 data = dist.rvs(size=(n_replications, n), random_state=rng) res = bootstrap((xp.asarray(data),), statistic=xp.mean, confidence_level=confidence_level, n_resamples=n_resamples, batch=50, method=method, axis=-1, rng=rng) ci = res.confidence_interval # ci contains vectors of lower and upper confidence interval bounds ci_contains_true = xp.count_nonzero((ci[0] < stat_true) & (stat_true < ci[1])) assert ci_contains_true == expected # ci_contains_true is not inconsistent with confidence_level pvalue = stats.binomtest(int(ci_contains_true), n_replications, confidence_level).pvalue assert pvalue > 0.1 tests_against_itself_2samp = {"basic": 892, "percentile": 890} @pytest.mark.slow @pytest.mark.parametrize("method, expected", tests_against_itself_2samp.items()) def test_bootstrap_against_itself_2samp(self, method, expected, xp): # The expected values in this test were generated using bootstrap # to check for unintended changes in behavior. The test also makes sure # that bootstrap works with multi-sample statistics and that the # `axis` argument works as expected / function is vectorized. rng = np.random.RandomState(0) n1 = 100 # size of sample 1 n2 = 120 # size of sample 2 n_resamples = 999 # number of bootstrap resamples used to form each CI confidence_level = 0.9 # The statistic we're interested in is the difference in means def my_stat(data1, data2, axis=-1): mean1 = xp.mean(data1, axis=axis) mean2 = xp.mean(data2, axis=axis) return mean1 - mean2 # The true difference in the means is -0.1 dist1 = stats.norm(loc=0, scale=1) dist2 = stats.norm(loc=0.1, scale=1) stat_true = float(dist1.mean() - dist2.mean()) # Do the same thing 1000 times. (The code is fully vectorized.) n_replications = 1000 data1 = dist1.rvs(size=(n_replications, n1), random_state=rng) data2 = dist2.rvs(size=(n_replications, n2), random_state=rng) res = bootstrap((xp.asarray(data1), xp.asarray(data2)), statistic=my_stat, confidence_level=confidence_level, n_resamples=n_resamples, batch=50, method=method, axis=-1, random_state=rng) ci = res.confidence_interval # ci contains vectors of lower and upper confidence interval bounds ci_contains_true = xp.count_nonzero((ci[0] < stat_true) & (stat_true < ci[1])) assert ci_contains_true == expected # ci_contains_true is not inconsistent with confidence_level pvalue = stats.binomtest(int(ci_contains_true), n_replications, confidence_level).pvalue assert pvalue > 0.1 @pytest.mark.parametrize("method", ["basic", "percentile", "BCa"]) @pytest.mark.parametrize("axis", [0, 1]) @pytest.mark.parametrize("dtype", ['float32', 'float64']) def test_bootstrap_vectorized_3samp(self, method, axis, dtype, xp): def statistic(*data, axis=0): # an arbitrary, vectorized statistic return sum(xp.mean(sample, axis=axis) for sample in data) def statistic_1d(*data): # the same statistic, not vectorized for sample in data: assert sample.ndim == 1 return np.asarray(sum(sample.mean() for sample in data), dtype=dtype) rng = np.random.RandomState(0) x = rng.rand(4, 5).astype(dtype) y = rng.rand(4, 5).astype(dtype) z = rng.rand(4, 5).astype(dtype) res1 = bootstrap((xp.asarray(x), xp.asarray(y), xp.asarray(z)), statistic, vectorized=True, axis=axis, n_resamples=100, method=method, rng=0) res2 = bootstrap((x, y, z), statistic_1d, vectorized=False, axis=axis, n_resamples=100, method=method, rng=0) rtol = 1e-6 if dtype == 'float32' else 1e-14 xp_assert_close(res1.confidence_interval.low, xp.asarray(res2.confidence_interval.low), rtol=rtol) xp_assert_close(res1.confidence_interval.high, xp.asarray(res2.confidence_interval.high), rtol=rtol) xp_assert_close(res1.standard_error, xp.asarray(res2.standard_error), rtol=rtol) @pytest.mark.xfail_on_32bit("Failure is not concerning; see gh-14107") @pytest.mark.parametrize("method", ["basic", "percentile", "BCa"]) @pytest.mark.parametrize("axis", [0, 1]) def test_bootstrap_vectorized_1samp(self, method, axis, xp): def statistic(x, axis=0): # an arbitrary, vectorized statistic return xp.mean(x, axis=axis) def statistic_1d(x): # the same statistic, not vectorized assert x.ndim == 1 return x.mean(axis=0) rng = np.random.default_rng(7939952824) x = rng.random((2, 3)) res1 = bootstrap((xp.asarray(x),), statistic, vectorized=True, axis=axis, n_resamples=100, batch=None, method=method, rng=0) res2 = bootstrap((x,), statistic_1d, vectorized=False, axis=axis, n_resamples=100, batch=10, method=method, rng=0) xp_assert_close(res1.confidence_interval.low, xp.asarray(res2.confidence_interval.low)) xp_assert_close(res1.confidence_interval.high, xp.asarray(res2.confidence_interval.high)) xp_assert_close(res1.standard_error, xp.asarray(res2.standard_error)) @pytest.mark.parametrize("method", ["basic", "percentile", "BCa"]) def test_bootstrap_degenerate(self, method, xp): data = xp.full((35,), 10000.) if method == "BCa": with np.errstate(invalid='ignore'): msg = "The BCa confidence interval cannot be calculated" with pytest.warns(stats.DegenerateDataWarning, match=msg): res = bootstrap([data, ], xp.mean, method=method) xp_assert_equal(res.confidence_interval.low, xp.asarray(xp.nan)) xp_assert_equal(res.confidence_interval.high, xp.asarray(xp.nan)) else: res = bootstrap([data, ], xp.mean, method=method) xp_assert_equal(res.confidence_interval.low, xp.asarray(10000.)) xp_assert_equal(res.confidence_interval.high, xp.asarray(10000.)) xp_assert_equal(res.standard_error, xp.asarray(0.)) @pytest.mark.parametrize("method", ["BCa", "basic", "percentile"]) def test_bootstrap_gh15678(self, method, xp): # Check that gh-15678 is fixed: when statistic function returned a Python # float, method="BCa" failed when trying to add a dimension to the float rng = np.random.default_rng(354645618886684) dist = stats.norm(loc=2, scale=4) data = dist.rvs(size=100, random_state=rng) res = bootstrap((xp.asarray(data),), stats.skew, method=method, n_resamples=100, rng=np.random.default_rng(9563)) # this always worked because np.apply_along_axis returns NumPy data type ref = bootstrap((data,), stats.skew, method=method, n_resamples=100, rng=np.random.default_rng(9563), vectorized=False) xp_assert_close(res.confidence_interval.low, xp.asarray(ref.confidence_interval.low)) xp_assert_close(res.confidence_interval.high, xp.asarray(ref.confidence_interval.high)) xp_assert_close(res.standard_error, xp.asarray(ref.standard_error)) def test_bootstrap_min(self, xp): # Check that gh-15883 is fixed: percentileofscore should # behave according to the 'mean' behavior and not trigger nan for BCa rng = np.random.default_rng(1891289180021102) dist = stats.norm(loc=2, scale=4) data = dist.rvs(size=100, random_state=rng) true_min = np.min(data) data = xp.asarray(data) res = bootstrap((data,), xp.min, method="BCa", n_resamples=100, rng=np.random.default_rng(3942)) xp_assert_equal(res.confidence_interval.low, xp.asarray(true_min)) res2 = bootstrap((-data,), xp.max, method="BCa", n_resamples=100, rng=np.random.default_rng(3942)) xp_assert_close(-res.confidence_interval.low, res2.confidence_interval.high) xp_assert_close(-res.confidence_interval.high, res2.confidence_interval.low) @pytest.mark.parametrize("additional_resamples", [0, 1000]) def test_re_bootstrap(self, additional_resamples, xp): # Test behavior of parameter `bootstrap_result` rng = np.random.default_rng(8958153316228384) x = rng.random(size=100) n1 = 1000 n2 = additional_resamples n3 = n1 + additional_resamples rng = np.random.default_rng(296689032789913033) res = stats.bootstrap((xp.asarray(x),), xp.mean, n_resamples=n1, rng=rng, confidence_level=0.95, method='percentile') res = stats.bootstrap((xp.asarray(x),), xp.mean, n_resamples=n2, rng=rng, confidence_level=0.90, method='BCa', bootstrap_result=res) rng = np.random.default_rng(296689032789913033) ref = stats.bootstrap((xp.asarray(x),), xp.mean, n_resamples=n3, rng=rng, confidence_level=0.90, method='BCa') xp_assert_close(res.confidence_interval.low, xp.asarray(ref.confidence_interval.low), rtol=1e-14) xp_assert_close(res.confidence_interval.high, xp.asarray(ref.confidence_interval.high), rtol=1e-14) xp_assert_close(res.standard_error, xp.asarray(ref.standard_error), rtol=1e-14) @pytest.mark.xfail_on_32bit("Sensitive to machine precision") @pytest.mark.parametrize("method", ['basic', 'percentile', 'BCa']) def test_bootstrap_alternative(self, method, xp): rng = np.random.default_rng(5894822712842015040) dist = stats.norm(loc=2, scale=4) data = (xp.asarray(dist.rvs(size=(100), random_state=rng)),) config = dict(data=data, statistic=xp.std, rng=rng, axis=-1) t = stats.bootstrap(**config, confidence_level=0.9) config.update(dict(n_resamples=0, bootstrap_result=t)) l = stats.bootstrap(**config, confidence_level=0.95, alternative='less') g = stats.bootstrap(**config, confidence_level=0.95, alternative='greater') xp_assert_close(l.confidence_interval.high, t.confidence_interval.high, rtol=1e-14) xp_assert_close(g.confidence_interval.low, t.confidence_interval.low, rtol=1e-14) assert xp.isinf(l.confidence_interval.low) and l.confidence_interval.low < 0 assert xp.isinf(g.confidence_interval.high) and g.confidence_interval.high > 0 with pytest.raises(ValueError, match='`alternative` must be one of'): stats.bootstrap(**config, alternative='ekki-ekki') def test_jackknife_resample(self, xp): shape = 3, 4, 5, 6 rng = np.random.default_rng(5274950392) x = rng.random(size=shape) y = next(_resampling._jackknife_resample(xp.asarray(x), xp=xp)) for i in range(shape[-1]): # each resample is indexed along second to last axis # (last axis is the one the statistic will be taken over / consumed) slc = y[..., i, :] expected = np.delete(x, i, axis=-1) xp_assert_equal(slc, xp.asarray(expected)) y2 = list(_resampling._jackknife_resample(xp.asarray(x), batch=2, xp=xp)) xp_assert_equal(xp.concat(y2, axis=-2), y) @pytest.mark.skip_xp_backends("array_api_strict", reason="Test uses ... + fancy indexing") @pytest.mark.parametrize("rng_name", ["RandomState", "default_rng"]) def test_bootstrap_resample(self, rng_name, xp): rng_gen = getattr(np.random, rng_name) rng1 = rng_gen(3949441460) rng2 = rng_gen(3949441460) n_resamples = 10 shape = 3, 4, 5, 6 rng = np.random.default_rng(5274950392) x = xp.asarray(rng.random(shape)) y = _resampling._bootstrap_resample(x, n_resamples, rng=rng1, xp=xp) for i in range(n_resamples): # each resample is indexed along second to last axis # (last axis is the one the statistic will be taken over / consumed) slc = y[..., i, :] js = xp.asarray(rng_integers(rng2, 0, shape[-1], shape[-1])) expected = x[..., js] xp_assert_equal(slc, expected) @pytest.mark.parametrize("score", [0, 0.5, 1]) @pytest.mark.parametrize("axis", [0, 1, 2]) def test_percentile_of_score(self, score, axis, xp): shape = 10, 20, 30 rng = np.random.default_rng(5903363153) x = rng.random(shape) dtype = xp_default_dtype(xp) p = _resampling._percentile_of_score(xp.asarray(x, dtype=dtype), xp.asarray(score, dtype=dtype), axis=-1, xp=xp) def vectorized_pos(a, score, axis): return np.apply_along_axis(stats.percentileofscore, axis, a, score) p2 = vectorized_pos(x, score, axis=-1)/100 xp_assert_close(p, xp.asarray(p2, dtype=dtype), rtol=1e-15) @pytest.mark.parametrize("axis", [0, 1, 2]) def test_vectorize_statistic(self, axis): # test that _vectorize_statistic vectorizes a statistic along `axis` # only used with NumPy arrays def statistic(*data, axis): # an arbitrary, vectorized statistic return sum(sample.mean(axis) for sample in data) def statistic_1d(*data): # the same statistic, not vectorized for sample in data: assert sample.ndim == 1 return statistic(*data, axis=0) # vectorize the non-vectorized statistic statistic2 = _resampling._vectorize_statistic(statistic_1d) rng = np.random.RandomState(0) x = rng.rand(4, 5, 6) y = rng.rand(4, 1, 6) z = rng.rand(1, 5, 6) res1 = statistic(x, y, z, axis=axis) res2 = statistic2(x, y, z, axis=axis) assert_allclose(res1, res2) @pytest.mark.slow @pytest.mark.parametrize("method", ["basic", "percentile", "BCa"]) def test_vector_valued_statistic(self, method, xp): # Generate 95% confidence interval around MLE of normal distribution # parameters. Repeat 100 times, each time on sample of size 100. # Check that confidence interval contains true parameters ~95 times. # Confidence intervals are estimated and stochastic; a test failure # does not necessarily indicate that something is wrong. More important # than values of `counts` below is that the shapes of the outputs are # correct. rng = np.random.default_rng(2196847219) params = 1, 0.5 sample = xp.asarray(rng.normal(*params, size=(100, 100))) def statistic(data, axis): return xp.stack([xp.mean(data, axis=axis), xp.std(data, axis=axis, correction=1)]) res = bootstrap((sample,), statistic, method=method, axis=-1, n_resamples=9999, batch=200, random_state=rng) params = xp.asarray([1, 0.5]) counts = xp.count_nonzero((res.confidence_interval.low.T < params) & (res.confidence_interval.high.T > params), axis=0) assert xp.all(counts >= 90) assert xp.all(counts <= 100) assert res.confidence_interval.low.shape == (2, 100) assert res.confidence_interval.high.shape == (2, 100) assert res.standard_error.shape == (2, 100) assert res.bootstrap_distribution.shape == (2, 100, 9999) @pytest.mark.slow @pytest.mark.filterwarnings('ignore::RuntimeWarning') def test_vector_valued_statistic_gh17715(self): # gh-17715 reported a mistake introduced in the extension of BCa to # multi-sample statistics; a `len` should have been `.shape[-1]`. Check # that this is resolved. # If this is resolved for NumPy, it's resolved for the rest, # let's only run this for NumPy. rng = np.random.default_rng(141921000979291141) def concordance(x, y, axis): xm = x.mean(axis) ym = y.mean(axis) cov = ((x - xm[..., None]) * (y - ym[..., None])).mean(axis) return (2 * cov) / (x.var(axis) + y.var(axis) + (xm - ym) ** 2) def statistic(tp, tn, fp, fn, axis): actual = tp + fp expected = tp + fn return np.nan_to_num(concordance(actual, expected, axis)) def statistic_extradim(*args, axis): return statistic(*args, axis)[np.newaxis, ...] data = [[4, 0, 0, 2], # (tp, tn, fp, fn) [2, 1, 2, 1], [0, 6, 0, 0], [0, 6, 3, 0], [0, 8, 1, 0]] data = np.array(data).T res = bootstrap(data, statistic_extradim, rng=rng, paired=True) ref = bootstrap(data, statistic, rng=rng, paired=True) assert_allclose(res.confidence_interval.low[0], ref.confidence_interval.low, atol=1e-15) assert_allclose(res.confidence_interval.high[0], ref.confidence_interval.high, atol=1e-15) # torch doesn't have T distribution CDF and can't fall back to NumPy on GPU @pytest.mark.skip_xp_backends(np_only=True, exceptions=['cupy']) def test_gh_20850(self, xp): rng = np.random.default_rng(2085020850) x = xp.asarray(rng.random((10, 2))) y = xp.asarray(rng.random((11, 2))) def statistic(x, y, axis): return stats.ttest_ind(x, y, axis=axis).statistic # The shapes do *not* need to be the same along axis stats.bootstrap((x, y), statistic) stats.bootstrap((x.T, y.T), statistic, axis=1) # But even when the shapes *are* the same along axis, the lengths # along other dimensions have to be the same (or `bootstrap` warns). message = "Array shapes are incompatible for broadcasting." with pytest.raises(ValueError, match=message): stats.bootstrap((x, y[:10, 0]), statistic) # this won't work after 1.16 stats.bootstrap((x, y[:10, 0:1]), statistic) # this will stats.bootstrap((x.T, y.T[0:1, :10]), statistic, axis=1) # this will # --- Test Monte Carlo Hypothesis Test --- # @make_xp_test_case(monte_carlo_test) class TestMonteCarloHypothesisTest: atol = 2.5e-2 # for comparing p-value def get_rvs(self, rvs_in, rs, dtype=None, xp=np): return lambda *args, **kwds: xp.asarray(rvs_in(*args, random_state=rs, **kwds), dtype=dtype) def get_statistic(self, xp): def statistic(x, axis): m = xp.mean(x, axis=axis) v = xp.var(x, axis=axis, correction=1) n = x.shape[axis] return m / (v/n)**0.5 # return stats.ttest_1samp(x, popmean=0., axis=axis).statistic) return statistic def test_input_validation(self, xp): # test that the appropriate error messages are raised for invalid input data = xp.asarray([1., 2., 3.]) def stat(x, axis=None): return xp.mean(x, axis=axis) message = "Array shapes are incompatible for broadcasting." temp = (xp.zeros((2, 5)), xp.zeros((3, 5))) rvs = (stats.norm.rvs, stats.norm.rvs) with pytest.raises(ValueError, match=message): monte_carlo_test(temp, rvs, lambda x, y, axis: 1, axis=-1) message = "`axis` must be an integer." with pytest.raises(ValueError, match=message): monte_carlo_test(data, stats.norm.rvs, stat, axis=1.5) message = "`vectorized` must be `True`, `False`, or `None`." with pytest.raises(ValueError, match=message): monte_carlo_test(data, stats.norm.rvs, stat, vectorized=1.5) message = "`rvs` must be callable or sequence of callables." with pytest.raises(TypeError, match=message): monte_carlo_test(data, None, stat) with pytest.raises(TypeError, match=message): temp = xp.asarray([[1., 2.], [3., 4.]]) monte_carlo_test(temp, [lambda x: x, None], stat) message = "If `rvs` is a sequence..." with pytest.raises(ValueError, match=message): temp = xp.asarray([[1., 2., 3.]]) monte_carlo_test(temp, [lambda x: x, lambda x: x], stat) message = "`statistic` must be callable." with pytest.raises(TypeError, match=message): monte_carlo_test(data, stats.norm.rvs, None) message = "`n_resamples` must be a positive integer." with pytest.raises(ValueError, match=message): monte_carlo_test(data, stats.norm.rvs, stat, n_resamples=-1000) message = "`n_resamples` must be a positive integer." with pytest.raises(ValueError, match=message): monte_carlo_test(data, stats.norm.rvs, stat, n_resamples=1000.5) message = "`batch` must be a positive integer or None." with pytest.raises(ValueError, match=message): monte_carlo_test(data, stats.norm.rvs, stat, batch=-1000) message = "`batch` must be a positive integer or None." with pytest.raises(ValueError, match=message): monte_carlo_test(data, stats.norm.rvs, stat, batch=1000.5) message = "`alternative` must be in..." with pytest.raises(ValueError, match=message): monte_carlo_test(data, stats.norm.rvs, stat, alternative='ekki') # *If* this raises a value error, make sure it has the intended message message = "Signature inspection of statistic" def rvs(size): return xp.asarray(stats.norm.rvs(size=size)) try: monte_carlo_test(data, rvs, xp.mean) except ValueError as e: assert str(e).startswith(message) def test_input_validation_xp(self, xp): def non_vectorized_statistic(x): return xp.mean(x) message = "`statistic` must be vectorized..." sample = xp.asarray([1., 2., 3.]) if is_numpy(xp): monte_carlo_test(sample, stats.norm.rvs, non_vectorized_statistic) return with pytest.raises(ValueError, match=message): monte_carlo_test(sample, stats.norm.rvs, non_vectorized_statistic) with pytest.raises(ValueError, match=message): monte_carlo_test(sample, stats.norm.rvs, xp.mean, vectorized=False) @pytest.mark.xslow def test_batch(self, xp): # make sure that the `batch` parameter is respected by checking the # maximum batch size provided in calls to `statistic` rng = np.random.default_rng(23492340193) x = xp.asarray(rng.standard_normal(size=10)) def statistic(x, axis): batch_size = 1 if x.ndim == 1 else x.shape[0] statistic.batch_size = max(batch_size, statistic.batch_size) statistic.counter += 1 return self.get_statistic(xp)(x, axis=axis) statistic.counter = 0 statistic.batch_size = 0 kwds = {'sample': x, 'statistic': statistic, 'n_resamples': 1000, 'vectorized': True} kwds['rvs'] = self.get_rvs(stats.norm.rvs, np.random.default_rng(328423), xp=xp) res1 = monte_carlo_test(batch=1, **kwds) assert_equal(statistic.counter, 1001) assert_equal(statistic.batch_size, 1) kwds['rvs'] = self.get_rvs(stats.norm.rvs, np.random.default_rng(328423), xp=xp) statistic.counter = 0 res2 = monte_carlo_test(batch=50, **kwds) assert_equal(statistic.counter, 21) assert_equal(statistic.batch_size, 50) kwds['rvs'] = self.get_rvs(stats.norm.rvs, np.random.default_rng(328423), xp=xp) statistic.counter = 0 res3 = monte_carlo_test(**kwds) assert_equal(statistic.counter, 2) assert_equal(statistic.batch_size, 1000) xp_assert_equal(res1.pvalue, res3.pvalue) xp_assert_equal(res2.pvalue, res3.pvalue) @pytest.mark.parametrize('axis', range(-3, 3)) def test_axis_dtype(self, axis, xp): # test that Nd-array samples are handled correctly for valid values # of the `axis` parameter; also make sure non-default dtype is maintained rng = np.random.default_rng(2389234) size = [2, 3, 4] size[axis] = 100 # Determine non-default dtype dtype_default = xp.asarray(1.).dtype dtype_str = 'float32'if ("64" in str(dtype_default)) else 'float64' dtype_np = getattr(np, dtype_str) dtype = getattr(xp, dtype_str) # ttest_1samp is CPU array-API compatible, but it would be good to # include CuPy in this test. We'll perform ttest_1samp with a # NumPy array, but all the rest with be done with fully array-API # compatible code. x = rng.standard_normal(size=size, dtype=dtype_np) expected = stats.ttest_1samp(x, popmean=0., axis=axis) x = xp.asarray(x, dtype=dtype) statistic = self.get_statistic(xp) rvs = self.get_rvs(stats.norm.rvs, rng, dtype=dtype, xp=xp) res = monte_carlo_test(x, rvs, statistic, vectorized=True, n_resamples=20000, axis=axis) ref_statistic = xp.asarray(expected.statistic, dtype=dtype) ref_pvalue = xp.asarray(expected.pvalue, dtype=dtype) xp_assert_close(res.statistic, ref_statistic) xp_assert_close(res.pvalue, ref_pvalue, atol=self.atol) @pytest.mark.parametrize('alternative', ("two-sided", "less", "greater")) def test_alternative(self, alternative, xp): # test that `alternative` is working as expected rng = np.random.default_rng(65723433) x = rng.standard_normal(size=30) ref = stats.ttest_1samp(x, 0., alternative=alternative) x = xp.asarray(x) statistic = self.get_statistic(xp) rvs = self.get_rvs(stats.norm.rvs, rng, xp=xp) res = monte_carlo_test(x, rvs, statistic, alternative=alternative) xp_assert_close(res.statistic, xp.asarray(ref.statistic)) xp_assert_close(res.pvalue, xp.asarray(ref.pvalue), atol=self.atol) # Tests below involve statistics that are not yet array-API compatible. # They can be converted when the statistics are converted. @pytest.mark.slow @pytest.mark.parametrize('alternative', ("less", "greater")) @pytest.mark.parametrize('a', np.linspace(-0.5, 0.5, 5)) # skewness def test_against_ks_1samp(self, alternative, a): # test that monte_carlo_test can reproduce pvalue of ks_1samp rng = np.random.default_rng(65723433) x = stats.skewnorm.rvs(a=a, size=30, random_state=rng) expected = stats.ks_1samp(x, stats.norm.cdf, alternative=alternative) def statistic1d(x): return stats.ks_1samp(x, stats.norm.cdf, mode='asymp', alternative=alternative).statistic norm_rvs = self.get_rvs(stats.norm.rvs, rng) res = monte_carlo_test(x, norm_rvs, statistic1d, n_resamples=1000, vectorized=False, alternative=alternative) assert_allclose(res.statistic, expected.statistic) if alternative == 'greater': assert_allclose(res.pvalue, expected.pvalue, atol=self.atol) elif alternative == 'less': assert_allclose(1-res.pvalue, expected.pvalue, atol=self.atol) @pytest.mark.parametrize('hypotest', (stats.skewtest, stats.kurtosistest)) @pytest.mark.parametrize('alternative', ("less", "greater", "two-sided")) @pytest.mark.parametrize('a', np.linspace(-2, 2, 5)) # skewness def test_against_normality_tests(self, hypotest, alternative, a): # test that monte_carlo_test can reproduce pvalue of normality tests rng = np.random.default_rng(85723405) x = stats.skewnorm.rvs(a=a, size=150, random_state=rng) expected = hypotest(x, alternative=alternative) def statistic(x, axis): return hypotest(x, axis=axis).statistic norm_rvs = self.get_rvs(stats.norm.rvs, rng) res = monte_carlo_test(x, norm_rvs, statistic, vectorized=True, alternative=alternative) assert_allclose(res.statistic, expected.statistic) assert_allclose(res.pvalue, expected.pvalue, atol=self.atol) @pytest.mark.parametrize('a', np.arange(-2, 3)) # skewness parameter def test_against_normaltest(self, a): # test that monte_carlo_test can reproduce pvalue of normaltest rng = np.random.default_rng(12340513) x = stats.skewnorm.rvs(a=a, size=150, random_state=rng) expected = stats.normaltest(x) def statistic(x, axis): return stats.normaltest(x, axis=axis).statistic norm_rvs = self.get_rvs(stats.norm.rvs, rng) res = monte_carlo_test(x, norm_rvs, statistic, vectorized=True, alternative='greater') assert_allclose(res.statistic, expected.statistic) assert_allclose(res.pvalue, expected.pvalue, atol=self.atol) @pytest.mark.xslow @pytest.mark.parametrize('a', np.linspace(-0.5, 0.5, 5)) # skewness def test_against_cramervonmises(self, a): # test that monte_carlo_test can reproduce pvalue of cramervonmises rng = np.random.default_rng(234874135) x = stats.skewnorm.rvs(a=a, size=30, random_state=rng) expected = stats.cramervonmises(x, stats.norm.cdf) def statistic1d(x): return stats.cramervonmises(x, stats.norm.cdf).statistic norm_rvs = self.get_rvs(stats.norm.rvs, rng) res = monte_carlo_test(x, norm_rvs, statistic1d, n_resamples=1000, vectorized=False, alternative='greater') assert_allclose(res.statistic, expected.statistic) assert_allclose(res.pvalue, expected.pvalue, atol=self.atol) @pytest.mark.slow @pytest.mark.parametrize('dist_name', ['norm', 'logistic']) @pytest.mark.parametrize('target_statistic', [0.6, 0.7, 0.8]) def test_against_anderson(self, dist_name, target_statistic): # test that monte_carlo_test can reproduce results of `anderson`. # find the skewness for which the sample statistic is within the range of # values tubulated by `anderson` def fun(a): rng = np.random.default_rng(394295467) x = stats.tukeylambda.rvs(a, size=100, random_state=rng) expected = stats.anderson(x, dist_name, method='interpolate') return expected.statistic - target_statistic with warnings.catch_warnings(): warnings.simplefilter("ignore", RuntimeWarning) sol = root(fun, x0=0) assert sol.success # get the significance level (p-value) associated with that critical # value a = sol.x[0] rng = np.random.default_rng(394295467) x = stats.tukeylambda.rvs(a, size=100, random_state=rng) expected = stats.anderson(x, dist_name, method='interpolate') expected_stat = expected.statistic expected_p = expected.pvalue # perform equivalent Monte Carlo test and compare results def statistic1d(x): return stats.anderson(x, dist_name, method='interpolate').statistic dist_rvs = self.get_rvs(getattr(stats, dist_name).rvs, rng) with warnings.catch_warnings(): warnings.simplefilter("ignore", RuntimeWarning) res = monte_carlo_test(x, dist_rvs, statistic1d, n_resamples=1000, vectorized=False, alternative='greater') assert_allclose(res.statistic, expected_stat) assert_allclose(res.pvalue, expected_p, atol=2*self.atol) def test_p_never_zero(self): # Use biased estimate of p-value to ensure that p-value is never zero # per monte_carlo_test reference [1] rng = np.random.default_rng(2190176673029737545) x = np.zeros(100) res = monte_carlo_test(x, rng.random, np.mean, vectorized=True, alternative='less') assert res.pvalue == 0.0001 def test_against_ttest_ind(self): # test that `monte_carlo_test` can reproduce results of `ttest_ind`. rng = np.random.default_rng(219017667302737545) data = rng.random(size=(2, 5)), rng.random(size=7) # broadcastable rvs = rng.normal, rng.normal def statistic(x, y, axis): return stats.ttest_ind(x, y, axis=axis).statistic res = stats.monte_carlo_test(data, rvs, statistic, axis=-1) ref = stats.ttest_ind(data[0], [data[1]], axis=-1) assert_allclose(res.statistic, ref.statistic) assert_allclose(res.pvalue, ref.pvalue, rtol=2e-2) def test_against_f_oneway(self): # test that `monte_carlo_test` can reproduce results of `f_oneway`. rng = np.random.default_rng(219017667302737545) data = (rng.random(size=(2, 100)), rng.random(size=(2, 101)), rng.random(size=(2, 102)), rng.random(size=(2, 103))) rvs = rng.normal, rng.normal, rng.normal, rng.normal def statistic(*args, axis): return stats.f_oneway(*args, axis=axis).statistic res = stats.monte_carlo_test(data, rvs, statistic, axis=-1, alternative='greater') ref = stats.f_oneway(*data, axis=-1) assert_allclose(res.statistic, ref.statistic) assert_allclose(res.pvalue, ref.pvalue, atol=1e-2) @pytest.mark.fail_slow(2) @pytest.mark.xfail_on_32bit("Statistic may not depend on sample order on 32-bit") def test_finite_precision_statistic(self): # Some statistics return numerically distinct values when the values # should be equal in theory. Test that `monte_carlo_test` accounts # for this in some way. rng = np.random.default_rng(2549824598234528) n_resamples = 9999 def rvs(size): return 1. * stats.bernoulli(p=0.333).rvs(size=size, random_state=rng) x = rvs(100) res = stats.monte_carlo_test(x, rvs, np.var, alternative='less', n_resamples=n_resamples) # show that having a tolerance matters c0 = np.sum(res.null_distribution <= res.statistic) c1 = np.sum(res.null_distribution <= res.statistic*(1+1e-15)) assert c0 != c1 assert res.pvalue == (c1 + 1)/(n_resamples + 1) @make_xp_test_case(power) class TestPower: def xp_normal(self, rng, *, xp, dtype=None): dtype = xp_default_dtype(xp) if dtype is None else dtype return lambda size: xp.asarray(rng.normal(size=size), dtype=dtype) def test_input_validation(self, xp): # test that the appropriate error messages are raised for invalid input rng = np.random.default_rng(8519895914314711673) test = stats.ttest_ind rvs = (self.xp_normal(rng, xp=xp), self.xp_normal(rng, xp=xp)) n_observations = (xp.asarray(10), xp.asarray(12)) message = "`vectorized` must be `True`, `False`, or `None`." with pytest.raises(ValueError, match=message): power(test, rvs, n_observations, vectorized=1.5) message = "`rvs` must be callable or sequence of callables." with pytest.raises(TypeError, match=message): power(test, None, n_observations) with pytest.raises(TypeError, match=message): power(test, (self.xp_normal(rng, xp=xp), 'ekki'), n_observations) message = "If `rvs` is a sequence..." with pytest.raises(ValueError, match=message): power(test, (self.xp_normal(rng, xp=xp),), n_observations) with pytest.raises(ValueError, match=message): power(test, rvs, (10,)) message = "`significance` must contain floats between 0 and 1." with pytest.raises(ValueError, match=message): power(test, rvs, n_observations, significance=2) with pytest.raises(ValueError, match=message): power(test, rvs, n_observations, significance=xp.linspace(-1, 1, 50)) message = "`kwargs` must be a dictionary" with pytest.raises(TypeError, match=message): power(test, rvs, n_observations, kwargs=(1, 2, 3)) message = "not be broadcast|Chunks do not add|Incompatible shapes" with pytest.raises((ValueError, RuntimeError), match=message): power(test, rvs, (xp.asarray([10, 11]), xp.asarray([12, 13, 14]))) with pytest.raises((ValueError, RuntimeError), match=message): power(test, rvs, (xp.asarray([10, 11]), xp.asarray([12, 13])), kwargs={'x': xp.asarray([1, 2, 3])}) message = "`test` must be callable" with pytest.raises(TypeError, match=message): power(None, rvs, n_observations) message = "`n_resamples` must be a positive integer" with pytest.raises(ValueError, match=message): power(test, rvs, n_observations, n_resamples=-10) with pytest.raises(ValueError, match=message): power(test, rvs, n_observations, n_resamples=10.5) message = "`batch` must be a positive integer" with pytest.raises(ValueError, match=message): power(test, rvs, n_observations, batch=-10) with pytest.raises(ValueError, match=message): power(test, rvs, n_observations, batch=10.5) @pytest.mark.slow def test_batch(self, xp): # make sure that the `batch` parameter is respected by checking the # maximum batch size provided in calls to `test` rng = np.random.default_rng(23492340193) def test(x, axis): batch_size = 1 if x.ndim == 1 else x.shape[0] test.batch_size = max(batch_size, test.batch_size) test.counter += 1 return stats.ttest_1samp(x, xp.asarray(0.), axis=axis).pvalue test.counter = 0 test.batch_size = 0 kwds = dict(test=test, n_observations=xp.asarray(10), n_resamples=1000) rng = np.random.default_rng(23492340193) res1 = power(**kwds, rvs=self.xp_normal(rng, xp=xp), batch=1) assert test.counter == 1000 assert test.batch_size == 1 rng = np.random.default_rng(23492340193) test.counter = 0 res2 = power(**kwds, rvs=self.xp_normal(rng, xp=xp), batch=50) assert test.counter == 20 assert test.batch_size == 50 rng = np.random.default_rng(23492340193) test.counter = 0 res3 = power(**kwds, rvs=self.xp_normal(rng, xp=xp), batch=1000) assert test.counter == 1 assert test.batch_size == 1000 xp_assert_equal(res1.power, res3.power) xp_assert_equal(res2.power, res3.power) @pytest.mark.slow def test_vectorization(self, xp): # Test that `power` is vectorized as expected rng = np.random.default_rng(25495254834552) alternatives = {-1: 'less', 0:'two-sided', 1: 'greater'} # Single vectorized call popmeans = xp.asarray([0, 0.2]) def test(x, alternative, axis=-1): # ensure that popmeans axis is zeroth and orthogonal to the rest popmeans_expanded = xpx.expand_dims(popmeans, axis=tuple(range(1, x.ndim + 1))) alternative = alternatives[int(alternative)] return stats.ttest_1samp(x, popmeans_expanded, alternative=alternative, axis=axis) # nx and kwargs broadcast against one another nx = xp.asarray([10, 15, 20, 50, 100])[:, xp.newaxis] kwargs = {'alternative': xp.asarray([-1, 0, 1])} # This dimension is added to the beginning significance = xp.asarray([0.01, 0.025, 0.05, 0.1]) res = stats.power(test, self.xp_normal(rng, xp=xp), nx, significance=significance, kwargs=kwargs) # Looping over all combinations ref = [] for significance_i in significance: for i in range(nx.shape[0]): nx_i = nx[i, ...] for alternative_i in kwargs['alternative']: for j in range(popmeans.shape[0]): popmean_j = popmeans[j] def test2(x, axis=-1): return stats.ttest_1samp(x, popmean_j, axis=axis, alternative=alternatives[int(alternative_i)]) tmp = stats.power(test2, self.xp_normal(rng, xp=xp), nx_i, significance=significance_i) ref.append(tmp.power) ref = xp.reshape(xp.stack(ref), res.power.shape) # Show that results are similar xp_assert_close(res.power, ref, rtol=2e-2, atol=1e-2) @pytest.mark.skip_xp_backends(cpu_only=True, exceptions=['cupy', 'jax.numpy', 'dask.array']) def test_ttest_ind_null(self, xp): # Check that the p-values of `ttest_ind` are uniformly distributed under # the null hypothesis rng = np.random.default_rng(254952548345528) test = stats.ttest_ind n_observations = (xp.asarray(rng.integers(10, 100, size=(10))), xp.asarray(rng.integers(10, 100, size=(10)))) rvs = self.xp_normal(rng, xp=xp), self.xp_normal(rng, xp=xp) significance = xp.asarray([0.01, 0.05, 0.1]) res = stats.power(test, rvs, n_observations, significance=significance) significance = xp.broadcast_to(significance[:, xp.newaxis], res.power.shape) xp_assert_close(res.power, significance, atol=1e-2) @pytest.mark.skip_xp_backends('array_api_strict', reason='currently combines integer and float arrays') @pytest.mark.skip_xp_backends(cpu_only=True, exceptions=['cupy', 'jax.numpy', 'dask.array']) def test_ttest_1samp_power(self, xp): # Check simulated ttest_1samp power against reference rng = np.random.default_rng(254952548345528) # Reference values computed with statmodels # import numpy as np # from statsmodels.stats.power import tt_solve_power # tt_solve_power = np.vectorize(tt_solve_power) # tt_solve_power([0.1, 0.5, 0.9], [[10], [20]], [[[0.01]], [[0.05]]]) ref = [[[0.0126515 , 0.10269751, 0.40415802], [0.01657775, 0.29734608, 0.86228288]], [[0.0592903 , 0.29317561, 0.71718121], [0.07094116, 0.56450441, 0.96815163]]] kwargs = {'popmean': xp.asarray([0.1, 0.5, 0.9])} n_observations = xp.asarray([[10], [20]]) significance = xp.asarray([0.01, 0.05]) res = stats.power(stats.ttest_1samp, self.xp_normal(rng, xp=xp), n_observations, significance=significance, kwargs=kwargs) xp_assert_close(res.power, xp.asarray(ref), atol=1e-2) @make_xp_test_case(permutation_test) class TestPermutationTest: rtol = 1e-14 def setup_method(self): self.rng = np.random.default_rng(7170559330470561044) # -- Input validation -- # def test_permutation_test_iv(self, xp): def stat(x, y, axis): return stats.ttest_ind((x, y), axis).statistic data = (xp.asarray([1, 2, 3]), xp.asarray([1, 2, 3])) message = "each sample in `data` must contain two or more ..." with pytest.raises(ValueError, match=message): permutation_test((data[0], xp.asarray([0])), stat) message = "`data` must be a tuple containing at least two samples" with pytest.raises(ValueError, match=message): permutation_test((1,), stat) with pytest.raises(TypeError, match=message): permutation_test(1, stat) message = "`axis` must be an integer." with pytest.raises(ValueError, match=message): permutation_test(data, stat, axis=1.5) message = "`permutation_type` must be in..." with pytest.raises(ValueError, match=message): permutation_test(data, stat, permutation_type="ekki") message = "`vectorized` must be `True`, `False`, or `None`." with pytest.raises(ValueError, match=message): permutation_test(data, stat, vectorized=1.5) message = "`n_resamples` must be a positive integer." with pytest.raises(ValueError, match=message): permutation_test(data, stat, n_resamples=-1000) message = "`n_resamples` must be a positive integer." with pytest.raises(ValueError, match=message): permutation_test(data, stat, n_resamples=1000.5) message = "`batch` must be a positive integer or None." with pytest.raises(ValueError, match=message): permutation_test(data, stat, batch=-1000) message = "`batch` must be a positive integer or None." with pytest.raises(ValueError, match=message): permutation_test(data, stat, batch=1000.5) message = "`alternative` must be in..." with pytest.raises(ValueError, match=message): permutation_test(data, stat, alternative='ekki') message = "SeedSequence expects int or sequence of ints" with pytest.raises(TypeError, match=message): permutation_test(data, stat, rng='herring') # -- Test Parameters -- # # SPEC-007 leave one call with seed to check it still works @pytest.mark.parametrize('random_state', [np.random.RandomState, np.random.default_rng]) @pytest.mark.parametrize('permutation_type', ['pairings', 'samples', 'independent']) def test_batch(self, permutation_type, random_state, xp): # make sure that the `batch` parameter is respected by checking the # maximum batch size provided in calls to `statistic` x = xp.asarray(self.rng.random(10)) y = xp.asarray(self.rng.random(10)) def statistic(x, y, axis): batch_size = 1 if x.ndim == 1 else x.shape[0] statistic.batch_size = max(batch_size, statistic.batch_size) statistic.counter += 1 return xp.mean(x, axis=axis) - xp.mean(y, axis=axis) statistic.counter = 0 statistic.batch_size = 0 kwds = {'n_resamples': 1000, 'permutation_type': permutation_type, 'vectorized': True} res1 = stats.permutation_test((x, y), statistic, batch=1, random_state=random_state(0), **kwds) assert statistic.counter == 1001 assert statistic.batch_size == 1 statistic.counter = 0 res2 = stats.permutation_test((x, y), statistic, batch=50, random_state=random_state(0), **kwds) assert statistic.counter == 21 assert statistic.batch_size == 50 statistic.counter = 0 res3 = stats.permutation_test((x, y), statistic, batch=1000, random_state=random_state(0), **kwds) assert statistic.counter == 2 assert statistic.batch_size == 1000 xp_assert_equal(res1.pvalue, res3.pvalue) xp_assert_equal(res2.pvalue, res3.pvalue) # SPEC-007 leave at least one call with seed to check it still works @pytest.mark.parametrize('random_state', [np.random.RandomState, np.random.default_rng]) @pytest.mark.parametrize('permutation_type, exact_size', [('pairings', special.factorial(3)**2), ('samples', 2**3), ('independent', special.binom(6, 3))]) def test_permutations(self, permutation_type, exact_size, random_state, xp): # make sure that the `permutations` parameter is respected by checking # the size of the null distribution x = xp.asarray(self.rng.random(3)) y = xp.asarray(self.rng.random(3)) def statistic(x, y, axis): return xp.mean(x, axis=axis) - xp.mean(y, axis=axis) kwds = {'permutation_type': permutation_type, 'vectorized': True} res = stats.permutation_test((x, y), statistic, n_resamples=3, random_state=random_state(0), **kwds) assert xp_size(res.null_distribution) == 3 res = stats.permutation_test((x, y), statistic, **kwds) assert xp_size(res.null_distribution) == exact_size # -- Randomized Permutation Tests -- # # To get reasonable accuracy, these next three tests are somewhat slow. # Originally, I had them passing for all combinations of permutation type, # alternative, and RNG, but that takes too long for CI. Instead, split # into three tests, each testing a particular combination of the three # parameters. def test_randomized_test_against_exact_both(self, xp): # check that the randomized and exact tests agree to reasonable # precision for permutation_type='both alternative, rng = 'less', 0 nx, ny, permutations = 8, 9, 24000 assert special.binom(nx + ny, nx) > permutations rng = np.random.default_rng(8235259808) x = xp.asarray(rng.standard_normal(size=nx)) y = xp.asarray(rng.standard_normal(size=ny)) data = x, y def statistic(x, y, axis): return xp.mean(x, axis=axis) - xp.mean(y, axis=axis) kwds = {'vectorized': True, 'permutation_type': 'independent', 'batch': 100, 'alternative': alternative, 'rng': rng} res = permutation_test(data, statistic, n_resamples=permutations, **kwds) res2 = permutation_test(data, statistic, n_resamples=xp.inf, **kwds) assert res.statistic == res2.statistic xp_assert_close(res.pvalue, res2.pvalue, atol=1e-2) @pytest.mark.slow() def test_randomized_test_against_exact_samples(self, xp): # check that the randomized and exact tests agree to reasonable # precision for permutation_type='samples' alternative, rng = 'greater', None nx, ny, permutations = 15, 15, 32000 assert 2**nx > permutations rng = np.random.default_rng(8235259808) x = xp.asarray(rng.standard_normal(size=nx)) y = xp.asarray(rng.standard_normal(size=ny)) data = x, y def statistic(x, y, axis): return xp.mean(x - y, axis=axis) kwds = {'vectorized': True, 'permutation_type': 'samples', 'batch': 100, 'alternative': alternative, 'rng': rng} res = permutation_test(data, statistic, n_resamples=permutations, **kwds) res2 = permutation_test(data, statistic, n_resamples=xp.inf, **kwds) assert res.statistic == res2.statistic xp_assert_close(res.pvalue, res2.pvalue, atol=1e-2) # I only need to skip torch on GPU because it doesn't have betaincc for pearsonr @pytest.mark.skip_xp_backends(cpu_only=True, exceptions=['cupy', 'jax.numpy']) def test_randomized_test_against_exact_pairings(self, xp): # check that the randomized and exact tests agree to reasonable # precision for permutation_type='pairings' alternative, rng = 'two-sided', self.rng nx, ny, permutations = 8, 8, 40000 assert special.factorial(nx) > permutations rng = np.random.default_rng(8235259808) x = xp.asarray(rng.standard_normal(size=nx)) y = xp.asarray(rng.standard_normal(size=ny)) data = [x] def statistic(x, axis): return stats.pearsonr(x, y, axis=axis).statistic kwds = {'vectorized': True, 'permutation_type': 'samples', 'batch': 100, 'alternative': alternative, 'rng': rng} res = permutation_test(data, statistic, n_resamples=permutations, **kwds) res2 = permutation_test(data, statistic, n_resamples=xp.inf, **kwds) assert res.statistic == res2.statistic xp_assert_close(res.pvalue, res2.pvalue, atol=1e-2) # -- Independent (Unpaired) Sample Tests -- # @pytest.mark.skip_xp_backends(eager_only=True) # TODO: change to jax_jit=False @pytest.mark.parametrize('alternative', ("less", "greater", "two-sided")) def test_against_ks_2samp(self, alternative, xp): x = self.rng.normal(size=4, scale=1) y = self.rng.normal(size=5, loc=3, scale=3) expected = stats.ks_2samp(x, y, alternative=alternative, mode='exact') def statistic(x, y, axis): # todo: use `xp` as backend when `ks_2samp` is translated to array API x, y = _xp_copy_to_numpy(x), _xp_copy_to_numpy(y) res = stats.ks_2samp(x, y, axis=axis, mode='asymp', alternative=alternative) res = xp.asarray(res.statistic) return res[()] if res.ndim == 0 else res # ks_2samp is always a one-tailed 'greater' test # it's the statistic that changes (D+ vs D- vs max(D+, D-)) x, y = xp.asarray(x), xp.asarray(y) res = permutation_test((x, y), statistic, n_resamples=np.inf, alternative='greater', rng=self.rng) xp_assert_close(res.statistic, xp.asarray(expected.statistic), rtol=self.rtol) xp_assert_close(res.pvalue, xp.asarray(expected.pvalue), rtol=self.rtol) @pytest.mark.skip_xp_backends(eager_only=True) # TODO: change to jax_jit=False @pytest.mark.parametrize('alternative', ("less", "greater", "two-sided")) def test_against_ansari(self, alternative, xp): x = self.rng.normal(size=4, scale=1) y = self.rng.normal(size=5, scale=3) # ansari has a different convention for 'alternative' alternative_correspondence = {"less": "greater", "greater": "less", "two-sided": "two-sided"} alternative_scipy = alternative_correspondence[alternative] expected = stats.ansari(x, y, alternative=alternative_scipy) def statistic(x, y, axis): # todo: use `xp` as backend when `ansari` is translated to array API x, y = _xp_copy_to_numpy(x), _xp_copy_to_numpy(y) res = stats.ansari(x, y, axis=axis) res = xp.asarray(res.statistic) return res[()] if res.ndim == 0 else res x, y = xp.asarray(x), xp.asarray(y) res = permutation_test((x, y), statistic, n_resamples=np.inf, alternative=alternative, rng=self.rng) xp_assert_close(res.statistic, xp.asarray(expected.statistic), rtol=self.rtol) xp_assert_close(res.pvalue, xp.asarray(expected.pvalue), rtol=self.rtol) @skip_xp_backends('cupy', reason='needs mannwhitneyu') @skip_xp_backends(eager_only=True) # mannwhitneyu does input validation @pytest.mark.parametrize('alternative', ("less", "greater", "two-sided")) def test_against_mannwhitneyu(self, alternative, xp): x = stats.uniform.rvs(size=(3, 5, 2), loc=0, random_state=self.rng) y = stats.uniform.rvs(size=(3, 5, 2), loc=0.05, random_state=self.rng) expected = stats.mannwhitneyu(x, y, axis=1, alternative=alternative) def statistic(x, y, axis): return stats.mannwhitneyu(x, y, axis=axis, method='asymptotic').statistic x, y = xp.asarray(x), xp.asarray(y) res = permutation_test((x, y), statistic, vectorized=True, n_resamples=xp.inf, alternative=alternative, axis=1, rng=self.rng) xp_assert_close(res.statistic, xp.asarray(expected.statistic), rtol=self.rtol) xp_assert_close(res.pvalue, xp.asarray(expected.pvalue), rtol=self.rtol) @skip_xp_backends('cupy', reason='needs cramervonmises_2samp') @skip_xp_backends(eager_only=True) # cramervonmises_2samp does input validation @skip_xp_backends(cpu_only=True) # torch doesn't have `kv` def test_against_cvm(self, xp): x = stats.norm.rvs(size=4, scale=1, random_state=self.rng) y = stats.norm.rvs(size=5, loc=3, scale=3, random_state=self.rng) expected = stats.cramervonmises_2samp(x, y, method='exact') def statistic(x, y, axis): res = stats.cramervonmises_2samp(x, y, axis=axis, method='asymptotic') return res.statistic # cramervonmises_2samp has only one alternative, greater x, y = xp.asarray(x), xp.asarray(y) res = permutation_test((x, y), statistic, n_resamples=np.inf, alternative='greater', rng=self.rng) xp_assert_close(res.statistic, xp.asarray(expected.statistic), rtol=self.rtol) xp_assert_close(res.pvalue, xp.asarray(expected.pvalue), rtol=self.rtol) @skip_xp_backends('cupy', reason='needs kruskal') @skip_xp_backends(eager_only=True) # kruskal does input validation @pytest.mark.parametrize('axis', (-1, 2)) def test_vectorized_nsamp_ptype_both(self, axis, xp): # statistic only available for NumPy # Test that permutation_test with permutation_type='independent' works # properly for a 3-sample statistic with nd array samples of different # (but compatible) shapes and ndims. Show that exact permutation test # and random permutation tests approximate SciPy's asymptotic pvalues # and that exact and random permutation test results are even closer # to one another (than they are to the asymptotic results). # Three samples, different (but compatible) shapes with different ndims rng = np.random.default_rng(6709265303529651545) x = rng.random(size=(3)) y = rng.random(size=(1, 3, 2)) z = rng.random(size=(2, 1, 4)) data = (x, y, z) expected = stats.kruskal(*data, axis=axis) # Define the statistic (and pvalue for comparison) def statistic(*data, axis): return stats.kruskal(*data, axis=axis).statistic # Calculate exact and randomized permutation results kwds = {'axis': axis, 'alternative': 'greater', 'permutation_type': 'independent', 'rng': rng} data = [xp.asarray(data_) for data_ in data] res = permutation_test(data, statistic, n_resamples=xp.inf, **kwds) res2 = permutation_test(data, statistic, n_resamples=1000, **kwds) # Check results xp_assert_close(res.statistic, xp.asarray(expected.statistic), rtol=self.rtol*5) xp_assert_close(res.statistic, res2.statistic, rtol=self.rtol*5) xp_assert_close(res.pvalue, xp.asarray(expected.pvalue), atol=6e-2) xp_assert_close(res.pvalue, res2.pvalue, atol=3e-2) # -- Paired-Sample Tests -- # @pytest.mark.skip_xp_backends(eager_only=True) # TODO: change to jax_jit=False @pytest.mark.parametrize('alternative', ("less", "greater", "two-sided")) def test_against_wilcoxon(self, alternative, xp): x = stats.uniform.rvs(size=(3, 6, 2), loc=0, random_state=self.rng) y = stats.uniform.rvs(size=(3, 6, 2), loc=0.05, random_state=self.rng) expected = stats.wilcoxon(x, y, alternative=alternative, axis=1) # We'll check both 1- and 2-sample versions of the same test; # we expect identical results to wilcoxon in all cases. def statistic_1samp_1d(z, axis): # todo: use `xp` as backend when `wilcoxon` is translated to array API # 'less' ensures we get the same of two statistics every time z = _xp_copy_to_numpy(z) res = stats.wilcoxon(z, alternative='less', axis=axis) res = xp.asarray(res.statistic) return res[()] if res.ndim == 0 else res def statistic_2samp_1d(x, y, axis): # todo: use `xp` as backend when `wilcoxon` is translated to array API x, y = _xp_copy_to_numpy(x), _xp_copy_to_numpy(y) res = stats.wilcoxon(x, y, alternative='less', axis=axis) res = xp.asarray(res.statistic) return res[()] if res.ndim == 0 else res x, y = xp.asarray(x), xp.asarray(y) kwds = {'axis': 1, 'alternative': alternative, 'permutation_type': 'samples', 'rng': self.rng, 'n_resamples': np.inf} res1 = permutation_test((x-y,), statistic_1samp_1d, **kwds) res2 = permutation_test((x, y), statistic_2samp_1d, **kwds) # `wilcoxon` returns a different statistic with 'two-sided' xp_assert_close(res1.statistic, res2.statistic, rtol=self.rtol) if alternative != 'two-sided': xp_assert_close(res2.statistic, xp.asarray(expected.statistic), rtol=self.rtol) xp_assert_close(res2.pvalue, xp.asarray(expected.pvalue), rtol=self.rtol) xp_assert_close(res1.pvalue, res2.pvalue, rtol=self.rtol) @pytest.mark.parametrize('alternative', ("less", "greater", "two-sided")) def test_against_binomtest(self, alternative, xp): x = self.rng.integers(0, 2, size=10) x[x == 0] = -1 # More naturally, the test would flip elements between 0 and one. # However, permutation_test will flip the _signs_ of the elements. # So we have to work with +1/-1 instead of 1/0. def statistic(x, axis=0): xp_ = array_namespace(x) return xp_.count_nonzero(x > 0, axis=axis) k, n, p = statistic(x), 10, 0.5 expected = stats.binomtest(k, n, p, alternative=alternative) res = stats.permutation_test((xp.asarray(x, dtype=xp.float64),), statistic, vectorized=True, permutation_type='samples', n_resamples=xp.inf, rng=self.rng, alternative=alternative) xp_assert_close(res.pvalue, xp.asarray(expected.pvalue), rtol=self.rtol) # -- Exact Association Tests -- # @pytest.mark.skip_xp_backends(eager_only=True) # TODO: change to jax_jit=False def test_against_kendalltau(self, xp): x = self.rng.normal(size=6) y = x + self.rng.normal(size=6) expected = stats.kendalltau(x, y, method='exact') def statistic(x, axis): # todo: use `xp` as backend when `kendalltau` is translated to array API x = _xp_copy_to_numpy(x) res = stats.kendalltau(x, y, method='asymptotic', axis=axis) res = xp.asarray(res.statistic) return res[()] if res.ndim == 0 else res # kendalltau currently has only one alternative, two-sided x = xp.asarray(x) res = permutation_test((x,), statistic, permutation_type='pairings', n_resamples=np.inf, rng=self.rng) xp_assert_close(res.statistic, xp.asarray(expected.statistic), rtol=self.rtol) xp_assert_close(res.pvalue, xp.asarray(expected.pvalue), rtol=self.rtol) @pytest.mark.parametrize('alternative', ('less', 'greater', 'two-sided')) def test_against_fisher_exact(self, alternative, xp): # x and y are binary random variables with some dependence rng = np.random.default_rng(6235696159000529929) x = (rng.random(7) > 0.6).astype(float) y = (rng.random(7) + 0.25*x > 0.6).astype(float) tab = stats.contingency.crosstab(x, y)[1] x, y = xp.asarray(x), xp.asarray(y) def statistic(x, axis): return xp.count_nonzero((x == 1) & (y == 1), axis=axis) res = permutation_test((x,), statistic, permutation_type='pairings', n_resamples=xp.inf, alternative=alternative, rng=rng) res2 = stats.fisher_exact(tab, alternative=alternative) xp_assert_close(res.pvalue, xp.asarray(res2.pvalue, dtype=x.dtype)) @pytest.mark.xslow() @pytest.mark.parametrize('axis', (-2, 1)) def test_vectorized_nsamp_ptype_samples(self, axis): # statistic only available for NumPy, and it's a pain to vectorize # Test that permutation_test with permutation_type='samples' works # properly for a 3-sample statistic with nd array samples of different # (but compatible) shapes and ndims. Show that exact permutation test # reproduces SciPy's exact pvalue and that random permutation test # approximates it. x = self.rng.random(size=(2, 4, 3)) y = self.rng.random(size=(1, 4, 3)) z = self.rng.random(size=(2, 4, 1)) x = stats.rankdata(x, axis=axis) y = stats.rankdata(y, axis=axis) z = stats.rankdata(z, axis=axis) y = y[0] # to check broadcast with different ndim data = (x, y, z) def statistic1d(*data): return stats.page_trend_test(data, ranked=True, method='asymptotic').statistic def pvalue1d(*data): return stats.page_trend_test(data, ranked=True, method='exact').pvalue statistic = _resampling._vectorize_statistic(statistic1d) pvalue = _resampling._vectorize_statistic(pvalue1d) expected_statistic = statistic(*np.broadcast_arrays(*data), axis=axis) expected_pvalue = pvalue(*np.broadcast_arrays(*data), axis=axis) # Let's forgive this use of an integer seed, please. kwds = {'vectorized': False, 'axis': axis, 'alternative': 'greater', 'permutation_type': 'pairings', 'rng': 0} res = permutation_test(data, statistic1d, n_resamples=np.inf, **kwds) res2 = permutation_test(data, statistic1d, n_resamples=5000, **kwds) assert_allclose(res.statistic, expected_statistic, rtol=self.rtol) assert_allclose(res.statistic, res2.statistic, rtol=self.rtol) assert_allclose(res.pvalue, expected_pvalue, rtol=self.rtol) assert_allclose(res.pvalue, res2.pvalue, atol=3e-2) # -- Test Against External References -- # tie_case_1 = {'x': [1, 2, 3, 4], 'y': [1.5, 2, 2.5], 'expected_less': 0.2000000000, 'expected_2sided': 0.4, # 2*expected_less 'expected_Pr_gte_S_mean': 0.3428571429, # see note below 'expected_statistic': 7.5, 'expected_avg': 9.142857, 'expected_std': 1.40698} tie_case_2 = {'x': [111, 107, 100, 99, 102, 106, 109, 108], 'y': [107, 108, 106, 98, 105, 103, 110, 105, 104], 'expected_less': 0.1555738379, 'expected_2sided': 0.3111476758, 'expected_Pr_gte_S_mean': 0.2969971205, # see note below 'expected_statistic': 32.5, 'expected_avg': 38.117647, 'expected_std': 5.172124} @pytest.mark.skip_xp_backends(eager_only=True) # TODO: change to jax_jit=False @pytest.mark.xslow() # only the second case is slow, really @pytest.mark.parametrize('case', (tie_case_1, tie_case_2)) def test_with_ties(self, case, xp): """ Results above from SAS PROC NPAR1WAY, e.g. DATA myData; INPUT X Y; CARDS; 1 1 1 2 1 3 1 4 2 1.5 2 2 2 2.5 ods graphics on; proc npar1way AB data=myData; class X; EXACT; run; ods graphics off; Note: SAS provides Pr >= |S-Mean|, which is different from our definition of a two-sided p-value. """ x = case['x'] y = case['y'] expected_statistic = xp.asarray(case['expected_statistic']) expected_less = xp.asarray(case['expected_less']) expected_2sided = xp.asarray(case['expected_2sided']) expected_Pr_gte_S_mean = xp.asarray(case['expected_Pr_gte_S_mean']) expected_avg = xp.asarray(case['expected_avg']) expected_std = xp.asarray(case['expected_std']) def statistic(x, y, axis): # todo: use `xp` as backend when `ansari` is translated to array API x, y = _xp_copy_to_numpy(x), _xp_copy_to_numpy(y) res = stats.ansari(x, y, axis=axis) res = xp.asarray(res.statistic) return res[()] if res.ndim == 0 else res dtype = xp_default_dtype(xp) x, y = xp.asarray(x, dtype=dtype), xp.asarray(y, dtype=dtype) with warnings.catch_warnings(): warnings.filterwarnings( "ignore", "Ties preclude use of exact statistic", UserWarning) res = permutation_test((x, y), statistic, n_resamples=np.inf, alternative='less') res2 = permutation_test((x, y), statistic, n_resamples=np.inf, alternative='two-sided') xp_assert_close(res.statistic, expected_statistic, rtol=self.rtol) xp_assert_close(res.pvalue, expected_less, atol=1e-10) xp_assert_close(res2.pvalue, expected_2sided, atol=1e-10) xp_assert_close(xp.mean(res2.null_distribution), expected_avg, rtol=1e-6) xp_assert_close(xp.std(res2.null_distribution), expected_std, rtol=1e-6) # SAS provides Pr >= |S-Mean|; might as well check against that, too S = res.statistic mean = xp.mean(res.null_distribution) n = res.null_distribution.shape[0] Pr_gte_S_mean = xp.astype(xp.count_nonzero( xp.abs(res.null_distribution-mean) >= xp.abs(S-mean)), S.dtype) / n xp_assert_close(Pr_gte_S_mean, expected_Pr_gte_S_mean) @pytest.mark.slow @pytest.mark.parametrize('alternative, expected_pvalue', (('less', 0.9708333333333), ('greater', 0.05138888888889), ('two-sided', 0.1027777777778))) # I only need to skip torch on GPU because it doesn't have betaincc for pearsonr @pytest.mark.skip_xp_backends(cpu_only=True, exceptions=['cupy', 'jax.numpy']) @pytest.mark.skip_xp_backends(eager_only=True) # TODO: change to jax_jit=False def test_against_spearmanr_in_R(self, alternative, expected_pvalue, xp): """ Results above from R cor.test, e.g. options(digits=16) x <- c(1.76405235, 0.40015721, 0.97873798, 2.2408932, 1.86755799, -0.97727788) y <- c(2.71414076, 0.2488, 0.87551913, 2.6514917, 2.01160156, 0.47699563) cor.test(x, y, method = "spearm", alternative = "t") """ # data comes from # np.random.seed(0) # x = stats.norm.rvs(size=6) # y = x + stats.norm.rvs(size=6) x = xp.asarray([1.76405235, 0.40015721, 0.97873798, 2.2408932, 1.86755799, -0.97727788]) y = xp.asarray([2.71414076, 0.2488, 0.87551913, 2.6514917, 2.01160156, 0.47699563]) expected_statistic = 0.7714285714285715 y = xp.asarray(stats.rankdata(_xp_copy_to_numpy(y))) def statistic(x, axis): # `spearmanr` is not array api compatible, but `pearsonr` is. So for now # use _xp_copy_to_numpy just for ranking so we can run this test w/ CuPy. # TODO: use `xp` as backend when cupy works with `rankdata` x = xp.asarray(stats.rankdata(_xp_copy_to_numpy(x), axis=axis)) return stats.pearsonr(x, y, axis=axis).statistic res = permutation_test((x,), statistic, permutation_type='pairings', n_resamples=xp.inf, alternative=alternative) xp_assert_close(res.statistic, xp.asarray(expected_statistic), rtol=self.rtol) xp_assert_close(res.pvalue, xp.asarray(expected_pvalue), atol=1e-13) @pytest.mark.parametrize("batch", (-1, 0)) def test_batch_generator_iv(self, batch): with pytest.raises(ValueError, match="`batch` must be positive."): list(_resampling._batch_generator([1, 2, 3], batch)) batch_generator_cases = [(range(0), 3, []), (range(6), 3, [[0, 1, 2], [3, 4, 5]]), (range(8), 3, [[0, 1, 2], [3, 4, 5], [6, 7]])] @pytest.mark.parametrize("iterable, batch, expected", batch_generator_cases) def test_batch_generator(self, iterable, batch, expected): got = list(_resampling._batch_generator(iterable, batch)) assert got == expected @pytest.mark.fail_slow(2) # I only need to skip torch on GPU because it doesn't have betaincc for pearsonr @pytest.mark.skip_xp_backends(cpu_only=True, exceptions=['cupy', 'jax.numpy']) def test_finite_precision_statistic(self, xp): # Some statistics return numerically distinct values when the values # should be equal in theory. Test that `permutation_test` accounts # for this in some way. x = xp.asarray([1., 2., 4., 3.], dtype=xp.float64) y = xp.asarray([2., 4., 6., 8.], dtype=xp.float64) def statistic(x, y, axis): return stats.pearsonr(x, y, axis=axis)[0] res = stats.permutation_test((x, y), statistic, permutation_type='pairings') r, pvalue, null = res.statistic, res.pvalue, res.null_distribution correct_p = 2 * float(xp.count_nonzero(null >= r - 1e-14)) / null.shape[0] assert pvalue == correct_p == 1/3 # Compare against other exact correlation tests using R corr.test # options(digits=16) # x = c(1, 2, 4, 3) # y = c(2, 4, 6, 8) # cor.test(x, y, alternative = "t", method = "spearman") # 0.333333333 # cor.test(x, y, alternative = "t", method = "kendall") # 0.333333333 def test_all_partitions_concatenated(): # make sure that _all_paritions_concatenated produces the correct number # of partitions of the data into samples of the given sizes and that # all are unique n = np.array([3, 2, 4], dtype=int) nc = np.cumsum(n) all_partitions = set() counter = 0 for partition_concatenated in _resampling._all_partitions_concatenated(n): counter += 1 partitioning = np.split(partition_concatenated, nc[:-1]) all_partitions.add(tuple([frozenset(i) for i in partitioning])) expected = np.prod([special.binom(sum(n[i:]), sum(n[i+1:])) for i in range(len(n)-1)]) assert_equal(counter, expected) assert_equal(len(all_partitions), expected) @pytest.mark.parametrize('fun_name', ['bootstrap', 'permutation_test', 'monte_carlo_test']) def test_parameter_vectorized(fun_name): # Check that parameter `vectorized` is working as desired for all # resampling functions. Results don't matter; just don't fail asserts. rng = np.random.default_rng(75245098234592) sample = rng.random(size=10) def rvs(size): # needed by `monte_carlo_test` return stats.norm.rvs(size=size, random_state=rng) fun_options = {'bootstrap': {'data': (sample,), 'rng': rng, 'method': 'percentile'}, 'permutation_test': {'data': (sample,), 'rng': rng, 'permutation_type': 'samples'}, 'monte_carlo_test': {'sample': sample, 'rvs': rvs}} common_options = {'n_resamples': 100} fun = getattr(stats, fun_name) options = fun_options[fun_name] options.update(common_options) def statistic(x, axis): assert x.ndim > 1 or np.array_equal(x, sample) return np.mean(x, axis=axis) fun(statistic=statistic, vectorized=None, **options) fun(statistic=statistic, vectorized=True, **options) def statistic(x): assert x.ndim == 1 return np.mean(x) fun(statistic=statistic, vectorized=None, **options) fun(statistic=statistic, vectorized=False, **options) class TestMonteCarloMethod: def test_rvs_and_random_state(self): message = "Use of `rvs` and `rng` are mutually exclusive." rng = np.random.default_rng(34982345) with pytest.raises(ValueError, match=message): stats.MonteCarloMethod(rvs=rng.random, rng=rng)