Coverage for /usr/lib/python3/dist-packages/scipy/stats/_resampling.py: 14%
453 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from __future__ import annotations
3import warnings
4import numpy as np
5from itertools import combinations, permutations, product
6from collections.abc import Sequence
7import inspect
9from scipy._lib._util import check_random_state, _rename_parameter
10from scipy.special import ndtr, ndtri, comb, factorial
11from scipy._lib._util import rng_integers
12from dataclasses import dataclass
13from ._common import ConfidenceInterval
14from ._axis_nan_policy import _broadcast_concatenate, _broadcast_arrays
15from ._warnings_errors import DegenerateDataWarning
17__all__ = ['bootstrap', 'monte_carlo_test', 'permutation_test']
20def _vectorize_statistic(statistic):
21 """Vectorize an n-sample statistic"""
22 # This is a little cleaner than np.nditer at the expense of some data
23 # copying: concatenate samples together, then use np.apply_along_axis
24 def stat_nd(*data, axis=0):
25 lengths = [sample.shape[axis] for sample in data]
26 split_indices = np.cumsum(lengths)[:-1]
27 z = _broadcast_concatenate(data, axis)
29 # move working axis to position 0 so that new dimensions in the output
30 # of `statistic` are _prepended_. ("This axis is removed, and replaced
31 # with new dimensions...")
32 z = np.moveaxis(z, axis, 0)
34 def stat_1d(z):
35 data = np.split(z, split_indices)
36 return statistic(*data)
38 return np.apply_along_axis(stat_1d, 0, z)[()]
39 return stat_nd
42def _jackknife_resample(sample, batch=None):
43 """Jackknife resample the sample. Only one-sample stats for now."""
44 n = sample.shape[-1]
45 batch_nominal = batch or n
47 for k in range(0, n, batch_nominal):
48 # col_start:col_end are the observations to remove
49 batch_actual = min(batch_nominal, n-k)
51 # jackknife - each row leaves out one observation
52 j = np.ones((batch_actual, n), dtype=bool)
53 np.fill_diagonal(j[:, k:k+batch_actual], False)
54 i = np.arange(n)
55 i = np.broadcast_to(i, (batch_actual, n))
56 i = i[j].reshape((batch_actual, n-1))
58 resamples = sample[..., i]
59 yield resamples
62def _bootstrap_resample(sample, n_resamples=None, random_state=None):
63 """Bootstrap resample the sample."""
64 n = sample.shape[-1]
66 # bootstrap - each row is a random resample of original observations
67 i = rng_integers(random_state, 0, n, (n_resamples, n))
69 resamples = sample[..., i]
70 return resamples
73def _percentile_of_score(a, score, axis):
74 """Vectorized, simplified `scipy.stats.percentileofscore`.
75 Uses logic of the 'mean' value of percentileofscore's kind parameter.
77 Unlike `stats.percentileofscore`, the percentile returned is a fraction
78 in [0, 1].
79 """
80 B = a.shape[axis]
81 return ((a < score).sum(axis=axis) + (a <= score).sum(axis=axis)) / (2 * B)
84def _percentile_along_axis(theta_hat_b, alpha):
85 """`np.percentile` with different percentile for each slice."""
86 # the difference between _percentile_along_axis and np.percentile is that
87 # np.percentile gets _all_ the qs for each axis slice, whereas
88 # _percentile_along_axis gets the q corresponding with each axis slice
89 shape = theta_hat_b.shape[:-1]
90 alpha = np.broadcast_to(alpha, shape)
91 percentiles = np.zeros_like(alpha, dtype=np.float64)
92 for indices, alpha_i in np.ndenumerate(alpha):
93 if np.isnan(alpha_i):
94 # e.g. when bootstrap distribution has only one unique element
95 msg = (
96 "The BCa confidence interval cannot be calculated."
97 " This problem is known to occur when the distribution"
98 " is degenerate or the statistic is np.min."
99 )
100 warnings.warn(DegenerateDataWarning(msg))
101 percentiles[indices] = np.nan
102 else:
103 theta_hat_b_i = theta_hat_b[indices]
104 percentiles[indices] = np.percentile(theta_hat_b_i, alpha_i)
105 return percentiles[()] # return scalar instead of 0d array
108def _bca_interval(data, statistic, axis, alpha, theta_hat_b, batch):
109 """Bias-corrected and accelerated interval."""
110 # closely follows [1] 14.3 and 15.4 (Eq. 15.36)
112 # calculate z0_hat
113 theta_hat = np.asarray(statistic(*data, axis=axis))[..., None]
114 percentile = _percentile_of_score(theta_hat_b, theta_hat, axis=-1)
115 z0_hat = ndtri(percentile)
117 # calculate a_hat
118 theta_hat_ji = [] # j is for sample of data, i is for jackknife resample
119 for j, sample in enumerate(data):
120 # _jackknife_resample will add an axis prior to the last axis that
121 # corresponds with the different jackknife resamples. Do the same for
122 # each sample of the data to ensure broadcastability. We need to
123 # create a copy of the list containing the samples anyway, so do this
124 # in the loop to simplify the code. This is not the bottleneck...
125 samples = [np.expand_dims(sample, -2) for sample in data]
126 theta_hat_i = []
127 for jackknife_sample in _jackknife_resample(sample, batch):
128 samples[j] = jackknife_sample
129 broadcasted = _broadcast_arrays(samples, axis=-1)
130 theta_hat_i.append(statistic(*broadcasted, axis=-1))
131 theta_hat_ji.append(theta_hat_i)
133 theta_hat_ji = [np.concatenate(theta_hat_i, axis=-1)
134 for theta_hat_i in theta_hat_ji]
136 n_j = [theta_hat_i.shape[-1] for theta_hat_i in theta_hat_ji]
138 theta_hat_j_dot = [theta_hat_i.mean(axis=-1, keepdims=True)
139 for theta_hat_i in theta_hat_ji]
141 U_ji = [(n - 1) * (theta_hat_dot - theta_hat_i)
142 for theta_hat_dot, theta_hat_i, n
143 in zip(theta_hat_j_dot, theta_hat_ji, n_j)]
145 nums = [(U_i**3).sum(axis=-1)/n**3 for U_i, n in zip(U_ji, n_j)]
146 dens = [(U_i**2).sum(axis=-1)/n**2 for U_i, n in zip(U_ji, n_j)]
147 a_hat = 1/6 * sum(nums) / sum(dens)**(3/2)
149 # calculate alpha_1, alpha_2
150 z_alpha = ndtri(alpha)
151 z_1alpha = -z_alpha
152 num1 = z0_hat + z_alpha
153 alpha_1 = ndtr(z0_hat + num1/(1 - a_hat*num1))
154 num2 = z0_hat + z_1alpha
155 alpha_2 = ndtr(z0_hat + num2/(1 - a_hat*num2))
156 return alpha_1, alpha_2, a_hat # return a_hat for testing
159def _bootstrap_iv(data, statistic, vectorized, paired, axis, confidence_level,
160 alternative, n_resamples, batch, method, bootstrap_result,
161 random_state):
162 """Input validation and standardization for `bootstrap`."""
164 if vectorized not in {True, False, None}:
165 raise ValueError("`vectorized` must be `True`, `False`, or `None`.")
167 if vectorized is None:
168 vectorized = 'axis' in inspect.signature(statistic).parameters
170 if not vectorized:
171 statistic = _vectorize_statistic(statistic)
173 axis_int = int(axis)
174 if axis != axis_int:
175 raise ValueError("`axis` must be an integer.")
177 n_samples = 0
178 try:
179 n_samples = len(data)
180 except TypeError:
181 raise ValueError("`data` must be a sequence of samples.")
183 if n_samples == 0:
184 raise ValueError("`data` must contain at least one sample.")
186 data_iv = []
187 for sample in data:
188 sample = np.atleast_1d(sample)
189 if sample.shape[axis_int] <= 1:
190 raise ValueError("each sample in `data` must contain two or more "
191 "observations along `axis`.")
192 sample = np.moveaxis(sample, axis_int, -1)
193 data_iv.append(sample)
195 if paired not in {True, False}:
196 raise ValueError("`paired` must be `True` or `False`.")
198 if paired:
199 n = data_iv[0].shape[-1]
200 for sample in data_iv[1:]:
201 if sample.shape[-1] != n:
202 message = ("When `paired is True`, all samples must have the "
203 "same length along `axis`")
204 raise ValueError(message)
206 # to generate the bootstrap distribution for paired-sample statistics,
207 # resample the indices of the observations
208 def statistic(i, axis=-1, data=data_iv, unpaired_statistic=statistic):
209 data = [sample[..., i] for sample in data]
210 return unpaired_statistic(*data, axis=axis)
212 data_iv = [np.arange(n)]
214 confidence_level_float = float(confidence_level)
216 alternative = alternative.lower()
217 alternatives = {'two-sided', 'less', 'greater'}
218 if alternative not in alternatives:
219 raise ValueError(f"`alternative` must be one of {alternatives}")
221 n_resamples_int = int(n_resamples)
222 if n_resamples != n_resamples_int or n_resamples_int < 0:
223 raise ValueError("`n_resamples` must be a non-negative integer.")
225 if batch is None:
226 batch_iv = batch
227 else:
228 batch_iv = int(batch)
229 if batch != batch_iv or batch_iv <= 0:
230 raise ValueError("`batch` must be a positive integer or None.")
232 methods = {'percentile', 'basic', 'bca'}
233 method = method.lower()
234 if method not in methods:
235 raise ValueError(f"`method` must be in {methods}")
237 message = "`bootstrap_result` must have attribute `bootstrap_distribution'"
238 if (bootstrap_result is not None
239 and not hasattr(bootstrap_result, "bootstrap_distribution")):
240 raise ValueError(message)
242 message = ("Either `bootstrap_result.bootstrap_distribution.size` or "
243 "`n_resamples` must be positive.")
244 if ((not bootstrap_result or
245 not bootstrap_result.bootstrap_distribution.size)
246 and n_resamples_int == 0):
247 raise ValueError(message)
249 random_state = check_random_state(random_state)
251 return (data_iv, statistic, vectorized, paired, axis_int,
252 confidence_level_float, alternative, n_resamples_int, batch_iv,
253 method, bootstrap_result, random_state)
256@dataclass
257class BootstrapResult:
258 """Result object returned by `scipy.stats.bootstrap`.
260 Attributes
261 ----------
262 confidence_interval : ConfidenceInterval
263 The bootstrap confidence interval as an instance of
264 `collections.namedtuple` with attributes `low` and `high`.
265 bootstrap_distribution : ndarray
266 The bootstrap distribution, that is, the value of `statistic` for
267 each resample. The last dimension corresponds with the resamples
268 (e.g. ``res.bootstrap_distribution.shape[-1] == n_resamples``).
269 standard_error : float or ndarray
270 The bootstrap standard error, that is, the sample standard
271 deviation of the bootstrap distribution.
273 """
274 confidence_interval: ConfidenceInterval
275 bootstrap_distribution: np.ndarray
276 standard_error: float | np.ndarray
279def bootstrap(data, statistic, *, n_resamples=9999, batch=None,
280 vectorized=None, paired=False, axis=0, confidence_level=0.95,
281 alternative='two-sided', method='BCa', bootstrap_result=None,
282 random_state=None):
283 r"""
284 Compute a two-sided bootstrap confidence interval of a statistic.
286 When `method` is ``'percentile'`` and `alternative` is ``'two-sided'``,
287 a bootstrap confidence interval is computed according to the following
288 procedure.
290 1. Resample the data: for each sample in `data` and for each of
291 `n_resamples`, take a random sample of the original sample
292 (with replacement) of the same size as the original sample.
294 2. Compute the bootstrap distribution of the statistic: for each set of
295 resamples, compute the test statistic.
297 3. Determine the confidence interval: find the interval of the bootstrap
298 distribution that is
300 - symmetric about the median and
301 - contains `confidence_level` of the resampled statistic values.
303 While the ``'percentile'`` method is the most intuitive, it is rarely
304 used in practice. Two more common methods are available, ``'basic'``
305 ('reverse percentile') and ``'BCa'`` ('bias-corrected and accelerated');
306 they differ in how step 3 is performed.
308 If the samples in `data` are taken at random from their respective
309 distributions :math:`n` times, the confidence interval returned by
310 `bootstrap` will contain the true value of the statistic for those
311 distributions approximately `confidence_level`:math:`\, \times \, n` times.
313 Parameters
314 ----------
315 data : sequence of array-like
316 Each element of data is a sample from an underlying distribution.
317 statistic : callable
318 Statistic for which the confidence interval is to be calculated.
319 `statistic` must be a callable that accepts ``len(data)`` samples
320 as separate arguments and returns the resulting statistic.
321 If `vectorized` is set ``True``,
322 `statistic` must also accept a keyword argument `axis` and be
323 vectorized to compute the statistic along the provided `axis`.
324 n_resamples : int, default: ``9999``
325 The number of resamples performed to form the bootstrap distribution
326 of the statistic.
327 batch : int, optional
328 The number of resamples to process in each vectorized call to
329 `statistic`. Memory usage is O(`batch`*``n``), where ``n`` is the
330 sample size. Default is ``None``, in which case ``batch = n_resamples``
331 (or ``batch = max(n_resamples, n)`` for ``method='BCa'``).
332 vectorized : bool, optional
333 If `vectorized` is set ``False``, `statistic` will not be passed
334 keyword argument `axis` and is expected to calculate the statistic
335 only for 1D samples. If ``True``, `statistic` will be passed keyword
336 argument `axis` and is expected to calculate the statistic along `axis`
337 when passed an ND sample array. If ``None`` (default), `vectorized`
338 will be set ``True`` if ``axis`` is a parameter of `statistic`. Use of
339 a vectorized statistic typically reduces computation time.
340 paired : bool, default: ``False``
341 Whether the statistic treats corresponding elements of the samples
342 in `data` as paired.
343 axis : int, default: ``0``
344 The axis of the samples in `data` along which the `statistic` is
345 calculated.
346 confidence_level : float, default: ``0.95``
347 The confidence level of the confidence interval.
348 alternative : {'two-sided', 'less', 'greater'}, default: ``'two-sided'``
349 Choose ``'two-sided'`` (default) for a two-sided confidence interval,
350 ``'less'`` for a one-sided confidence interval with the lower bound
351 at ``-np.inf``, and ``'greater'`` for a one-sided confidence interval
352 with the upper bound at ``np.inf``. The other bound of the one-sided
353 confidence intervals is the same as that of a two-sided confidence
354 interval with `confidence_level` twice as far from 1.0; e.g. the upper
355 bound of a 95% ``'less'`` confidence interval is the same as the upper
356 bound of a 90% ``'two-sided'`` confidence interval.
357 method : {'percentile', 'basic', 'bca'}, default: ``'BCa'``
358 Whether to return the 'percentile' bootstrap confidence interval
359 (``'percentile'``), the 'basic' (AKA 'reverse') bootstrap confidence
360 interval (``'basic'``), or the bias-corrected and accelerated bootstrap
361 confidence interval (``'BCa'``).
362 bootstrap_result : BootstrapResult, optional
363 Provide the result object returned by a previous call to `bootstrap`
364 to include the previous bootstrap distribution in the new bootstrap
365 distribution. This can be used, for example, to change
366 `confidence_level`, change `method`, or see the effect of performing
367 additional resampling without repeating computations.
368 random_state : {None, int, `numpy.random.Generator`,
369 `numpy.random.RandomState`}, optional
371 Pseudorandom number generator state used to generate resamples.
373 If `random_state` is ``None`` (or `np.random`), the
374 `numpy.random.RandomState` singleton is used.
375 If `random_state` is an int, a new ``RandomState`` instance is used,
376 seeded with `random_state`.
377 If `random_state` is already a ``Generator`` or ``RandomState``
378 instance then that instance is used.
380 Returns
381 -------
382 res : BootstrapResult
383 An object with attributes:
385 confidence_interval : ConfidenceInterval
386 The bootstrap confidence interval as an instance of
387 `collections.namedtuple` with attributes `low` and `high`.
388 bootstrap_distribution : ndarray
389 The bootstrap distribution, that is, the value of `statistic` for
390 each resample. The last dimension corresponds with the resamples
391 (e.g. ``res.bootstrap_distribution.shape[-1] == n_resamples``).
392 standard_error : float or ndarray
393 The bootstrap standard error, that is, the sample standard
394 deviation of the bootstrap distribution.
396 Warns
397 -----
398 `~scipy.stats.DegenerateDataWarning`
399 Generated when ``method='BCa'`` and the bootstrap distribution is
400 degenerate (e.g. all elements are identical).
402 Notes
403 -----
404 Elements of the confidence interval may be NaN for ``method='BCa'`` if
405 the bootstrap distribution is degenerate (e.g. all elements are identical).
406 In this case, consider using another `method` or inspecting `data` for
407 indications that other analysis may be more appropriate (e.g. all
408 observations are identical).
410 References
411 ----------
412 .. [1] B. Efron and R. J. Tibshirani, An Introduction to the Bootstrap,
413 Chapman & Hall/CRC, Boca Raton, FL, USA (1993)
414 .. [2] Nathaniel E. Helwig, "Bootstrap Confidence Intervals",
415 http://users.stat.umn.edu/~helwig/notes/bootci-Notes.pdf
416 .. [3] Bootstrapping (statistics), Wikipedia,
417 https://en.wikipedia.org/wiki/Bootstrapping_%28statistics%29
419 Examples
420 --------
421 Suppose we have sampled data from an unknown distribution.
423 >>> import numpy as np
424 >>> rng = np.random.default_rng()
425 >>> from scipy.stats import norm
426 >>> dist = norm(loc=2, scale=4) # our "unknown" distribution
427 >>> data = dist.rvs(size=100, random_state=rng)
429 We are interested in the standard deviation of the distribution.
431 >>> std_true = dist.std() # the true value of the statistic
432 >>> print(std_true)
433 4.0
434 >>> std_sample = np.std(data) # the sample statistic
435 >>> print(std_sample)
436 3.9460644295563863
438 The bootstrap is used to approximate the variability we would expect if we
439 were to repeatedly sample from the unknown distribution and calculate the
440 statistic of the sample each time. It does this by repeatedly resampling
441 values *from the original sample* with replacement and calculating the
442 statistic of each resample. This results in a "bootstrap distribution" of
443 the statistic.
445 >>> import matplotlib.pyplot as plt
446 >>> from scipy.stats import bootstrap
447 >>> data = (data,) # samples must be in a sequence
448 >>> res = bootstrap(data, np.std, confidence_level=0.9,
449 ... random_state=rng)
450 >>> fig, ax = plt.subplots()
451 >>> ax.hist(res.bootstrap_distribution, bins=25)
452 >>> ax.set_title('Bootstrap Distribution')
453 >>> ax.set_xlabel('statistic value')
454 >>> ax.set_ylabel('frequency')
455 >>> plt.show()
457 The standard error quantifies this variability. It is calculated as the
458 standard deviation of the bootstrap distribution.
460 >>> res.standard_error
461 0.24427002125829136
462 >>> res.standard_error == np.std(res.bootstrap_distribution, ddof=1)
463 True
465 The bootstrap distribution of the statistic is often approximately normal
466 with scale equal to the standard error.
468 >>> x = np.linspace(3, 5)
469 >>> pdf = norm.pdf(x, loc=std_sample, scale=res.standard_error)
470 >>> fig, ax = plt.subplots()
471 >>> ax.hist(res.bootstrap_distribution, bins=25, density=True)
472 >>> ax.plot(x, pdf)
473 >>> ax.set_title('Normal Approximation of the Bootstrap Distribution')
474 >>> ax.set_xlabel('statistic value')
475 >>> ax.set_ylabel('pdf')
476 >>> plt.show()
478 This suggests that we could construct a 90% confidence interval on the
479 statistic based on quantiles of this normal distribution.
481 >>> norm.interval(0.9, loc=std_sample, scale=res.standard_error)
482 (3.5442759991341726, 4.3478528599786)
484 Due to central limit theorem, this normal approximation is accurate for a
485 variety of statistics and distributions underlying the samples; however,
486 the approximation is not reliable in all cases. Because `bootstrap` is
487 designed to work with arbitrary underlying distributions and statistics,
488 it uses more advanced techniques to generate an accurate confidence
489 interval.
491 >>> print(res.confidence_interval)
492 ConfidenceInterval(low=3.57655333533867, high=4.382043696342881)
494 If we sample from the original distribution 1000 times and form a bootstrap
495 confidence interval for each sample, the confidence interval
496 contains the true value of the statistic approximately 90% of the time.
498 >>> n_trials = 1000
499 >>> ci_contains_true_std = 0
500 >>> for i in range(n_trials):
501 ... data = (dist.rvs(size=100, random_state=rng),)
502 ... ci = bootstrap(data, np.std, confidence_level=0.9, n_resamples=1000,
503 ... random_state=rng).confidence_interval
504 ... if ci[0] < std_true < ci[1]:
505 ... ci_contains_true_std += 1
506 >>> print(ci_contains_true_std)
507 875
509 Rather than writing a loop, we can also determine the confidence intervals
510 for all 1000 samples at once.
512 >>> data = (dist.rvs(size=(n_trials, 100), random_state=rng),)
513 >>> res = bootstrap(data, np.std, axis=-1, confidence_level=0.9,
514 ... n_resamples=1000, random_state=rng)
515 >>> ci_l, ci_u = res.confidence_interval
517 Here, `ci_l` and `ci_u` contain the confidence interval for each of the
518 ``n_trials = 1000`` samples.
520 >>> print(ci_l[995:])
521 [3.77729695 3.75090233 3.45829131 3.34078217 3.48072829]
522 >>> print(ci_u[995:])
523 [4.88316666 4.86924034 4.32032996 4.2822427 4.59360598]
525 And again, approximately 90% contain the true value, ``std_true = 4``.
527 >>> print(np.sum((ci_l < std_true) & (std_true < ci_u)))
528 900
530 `bootstrap` can also be used to estimate confidence intervals of
531 multi-sample statistics, including those calculated by hypothesis
532 tests. `scipy.stats.mood` perform's Mood's test for equal scale parameters,
533 and it returns two outputs: a statistic, and a p-value. To get a
534 confidence interval for the test statistic, we first wrap
535 `scipy.stats.mood` in a function that accepts two sample arguments,
536 accepts an `axis` keyword argument, and returns only the statistic.
538 >>> from scipy.stats import mood
539 >>> def my_statistic(sample1, sample2, axis):
540 ... statistic, _ = mood(sample1, sample2, axis=-1)
541 ... return statistic
543 Here, we use the 'percentile' method with the default 95% confidence level.
545 >>> sample1 = norm.rvs(scale=1, size=100, random_state=rng)
546 >>> sample2 = norm.rvs(scale=2, size=100, random_state=rng)
547 >>> data = (sample1, sample2)
548 >>> res = bootstrap(data, my_statistic, method='basic', random_state=rng)
549 >>> print(mood(sample1, sample2)[0]) # element 0 is the statistic
550 -5.521109549096542
551 >>> print(res.confidence_interval)
552 ConfidenceInterval(low=-7.255994487314675, high=-4.016202624747605)
554 The bootstrap estimate of the standard error is also available.
556 >>> print(res.standard_error)
557 0.8344963846318795
559 Paired-sample statistics work, too. For example, consider the Pearson
560 correlation coefficient.
562 >>> from scipy.stats import pearsonr
563 >>> n = 100
564 >>> x = np.linspace(0, 10, n)
565 >>> y = x + rng.uniform(size=n)
566 >>> print(pearsonr(x, y)[0]) # element 0 is the statistic
567 0.9962357936065914
569 We wrap `pearsonr` so that it returns only the statistic.
571 >>> def my_statistic(x, y):
572 ... return pearsonr(x, y)[0]
574 We call `bootstrap` using ``paired=True``.
575 Also, since ``my_statistic`` isn't vectorized to calculate the statistic
576 along a given axis, we pass in ``vectorized=False``.
578 >>> res = bootstrap((x, y), my_statistic, vectorized=False, paired=True,
579 ... random_state=rng)
580 >>> print(res.confidence_interval)
581 ConfidenceInterval(low=0.9950085825848624, high=0.9971212407917498)
583 The result object can be passed back into `bootstrap` to perform additional
584 resampling:
586 >>> len(res.bootstrap_distribution)
587 9999
588 >>> res = bootstrap((x, y), my_statistic, vectorized=False, paired=True,
589 ... n_resamples=1001, random_state=rng,
590 ... bootstrap_result=res)
591 >>> len(res.bootstrap_distribution)
592 11000
594 or to change the confidence interval options:
596 >>> res2 = bootstrap((x, y), my_statistic, vectorized=False, paired=True,
597 ... n_resamples=0, random_state=rng, bootstrap_result=res,
598 ... method='percentile', confidence_level=0.9)
599 >>> np.testing.assert_equal(res2.bootstrap_distribution,
600 ... res.bootstrap_distribution)
601 >>> res.confidence_interval
602 ConfidenceInterval(low=0.9950035351407804, high=0.9971170323404578)
604 without repeating computation of the original bootstrap distribution.
606 """
607 # Input validation
608 args = _bootstrap_iv(data, statistic, vectorized, paired, axis,
609 confidence_level, alternative, n_resamples, batch,
610 method, bootstrap_result, random_state)
611 (data, statistic, vectorized, paired, axis, confidence_level,
612 alternative, n_resamples, batch, method, bootstrap_result,
613 random_state) = args
615 theta_hat_b = ([] if bootstrap_result is None
616 else [bootstrap_result.bootstrap_distribution])
618 batch_nominal = batch or n_resamples or 1
620 for k in range(0, n_resamples, batch_nominal):
621 batch_actual = min(batch_nominal, n_resamples-k)
622 # Generate resamples
623 resampled_data = []
624 for sample in data:
625 resample = _bootstrap_resample(sample, n_resamples=batch_actual,
626 random_state=random_state)
627 resampled_data.append(resample)
629 # Compute bootstrap distribution of statistic
630 theta_hat_b.append(statistic(*resampled_data, axis=-1))
631 theta_hat_b = np.concatenate(theta_hat_b, axis=-1)
633 # Calculate percentile interval
634 alpha = ((1 - confidence_level)/2 if alternative == 'two-sided'
635 else (1 - confidence_level))
636 if method == 'bca':
637 interval = _bca_interval(data, statistic, axis=-1, alpha=alpha,
638 theta_hat_b=theta_hat_b, batch=batch)[:2]
639 percentile_fun = _percentile_along_axis
640 else:
641 interval = alpha, 1-alpha
643 def percentile_fun(a, q):
644 return np.percentile(a=a, q=q, axis=-1)
646 # Calculate confidence interval of statistic
647 ci_l = percentile_fun(theta_hat_b, interval[0]*100)
648 ci_u = percentile_fun(theta_hat_b, interval[1]*100)
649 if method == 'basic': # see [3]
650 theta_hat = statistic(*data, axis=-1)
651 ci_l, ci_u = 2*theta_hat - ci_u, 2*theta_hat - ci_l
653 if alternative == 'less':
654 ci_l = np.full_like(ci_l, -np.inf)
655 elif alternative == 'greater':
656 ci_u = np.full_like(ci_u, np.inf)
658 return BootstrapResult(confidence_interval=ConfidenceInterval(ci_l, ci_u),
659 bootstrap_distribution=theta_hat_b,
660 standard_error=np.std(theta_hat_b, ddof=1, axis=-1))
663def _monte_carlo_test_iv(data, rvs, statistic, vectorized, n_resamples,
664 batch, alternative, axis):
665 """Input validation for `monte_carlo_test`."""
667 axis_int = int(axis)
668 if axis != axis_int:
669 raise ValueError("`axis` must be an integer.")
671 if vectorized not in {True, False, None}:
672 raise ValueError("`vectorized` must be `True`, `False`, or `None`.")
674 if not isinstance(rvs, Sequence):
675 rvs = (rvs,)
676 data = (data,)
677 for rvs_i in rvs:
678 if not callable(rvs_i):
679 raise TypeError("`rvs` must be callable or sequence of callables.")
681 if not len(rvs) == len(data):
682 message = "If `rvs` is a sequence, `len(rvs)` must equal `len(data)`."
683 raise ValueError(message)
685 if not callable(statistic):
686 raise TypeError("`statistic` must be callable.")
688 if vectorized is None:
689 vectorized = 'axis' in inspect.signature(statistic).parameters
691 if not vectorized:
692 statistic_vectorized = _vectorize_statistic(statistic)
693 else:
694 statistic_vectorized = statistic
696 data = _broadcast_arrays(data, axis)
697 data_iv = []
698 for sample in data:
699 sample = np.atleast_1d(sample)
700 sample = np.moveaxis(sample, axis_int, -1)
701 data_iv.append(sample)
703 n_resamples_int = int(n_resamples)
704 if n_resamples != n_resamples_int or n_resamples_int <= 0:
705 raise ValueError("`n_resamples` must be a positive integer.")
707 if batch is None:
708 batch_iv = batch
709 else:
710 batch_iv = int(batch)
711 if batch != batch_iv or batch_iv <= 0:
712 raise ValueError("`batch` must be a positive integer or None.")
714 alternatives = {'two-sided', 'greater', 'less'}
715 alternative = alternative.lower()
716 if alternative not in alternatives:
717 raise ValueError(f"`alternative` must be in {alternatives}")
719 return (data_iv, rvs, statistic_vectorized, vectorized, n_resamples_int,
720 batch_iv, alternative, axis_int)
723@dataclass
724class MonteCarloTestResult:
725 """Result object returned by `scipy.stats.monte_carlo_test`.
727 Attributes
728 ----------
729 statistic : float or ndarray
730 The observed test statistic of the sample.
731 pvalue : float or ndarray
732 The p-value for the given alternative.
733 null_distribution : ndarray
734 The values of the test statistic generated under the null
735 hypothesis.
736 """
737 statistic: float | np.ndarray
738 pvalue: float | np.ndarray
739 null_distribution: np.ndarray
742@_rename_parameter('sample', 'data')
743def monte_carlo_test(data, rvs, statistic, *, vectorized=None,
744 n_resamples=9999, batch=None, alternative="two-sided",
745 axis=0):
746 r"""Perform a Monte Carlo hypothesis test.
748 `data` contains a sample or a sequence of one or more samples. `rvs`
749 specifies the distribution(s) of the sample(s) in `data` under the null
750 hypothesis. The value of `statistic` for the given `data` is compared
751 against a Monte Carlo null distribution: the value of the statistic for
752 each of `n_resamples` sets of samples generated using `rvs`. This gives
753 the p-value, the probability of observing such an extreme value of the
754 test statistic under the null hypothesis.
756 Parameters
757 ----------
758 data : array-like or sequence of array-like
759 An array or sequence of arrays of observations.
760 rvs : callable or tuple of callables
761 A callable or sequence of callables that generates random variates
762 under the null hypothesis. Each element of `rvs` must be a callable
763 that accepts keyword argument ``size`` (e.g. ``rvs(size=(m, n))``) and
764 returns an N-d array sample of that shape. If `rvs` is a sequence, the
765 number of callables in `rvs` must match the number of samples in
766 `data`, i.e. ``len(rvs) == len(data)``. If `rvs` is a single callable,
767 `data` is treated as a single sample.
768 statistic : callable
769 Statistic for which the p-value of the hypothesis test is to be
770 calculated. `statistic` must be a callable that accepts a sample
771 (e.g. ``statistic(sample)``) or ``len(rvs)`` separate samples (e.g.
772 ``statistic(samples1, sample2)`` if `rvs` contains two callables and
773 `data` contains two samples) and returns the resulting statistic.
774 If `vectorized` is set ``True``, `statistic` must also accept a keyword
775 argument `axis` and be vectorized to compute the statistic along the
776 provided `axis` of the samples in `data`.
777 vectorized : bool, optional
778 If `vectorized` is set ``False``, `statistic` will not be passed
779 keyword argument `axis` and is expected to calculate the statistic
780 only for 1D samples. If ``True``, `statistic` will be passed keyword
781 argument `axis` and is expected to calculate the statistic along `axis`
782 when passed ND sample arrays. If ``None`` (default), `vectorized`
783 will be set ``True`` if ``axis`` is a parameter of `statistic`. Use of
784 a vectorized statistic typically reduces computation time.
785 n_resamples : int, default: 9999
786 Number of samples drawn from each of the callables of `rvs`.
787 Equivalently, the number statistic values under the null hypothesis
788 used as the Monte Carlo null distribution.
789 batch : int, optional
790 The number of Monte Carlo samples to process in each call to
791 `statistic`. Memory usage is O(`batch`*``sample.size[axis]``). Default
792 is ``None``, in which case `batch` equals `n_resamples`.
793 alternative : {'two-sided', 'less', 'greater'}
794 The alternative hypothesis for which the p-value is calculated.
795 For each alternative, the p-value is defined as follows.
797 - ``'greater'`` : the percentage of the null distribution that is
798 greater than or equal to the observed value of the test statistic.
799 - ``'less'`` : the percentage of the null distribution that is
800 less than or equal to the observed value of the test statistic.
801 - ``'two-sided'`` : twice the smaller of the p-values above.
803 axis : int, default: 0
804 The axis of `data` (or each sample within `data`) over which to
805 calculate the statistic.
807 Returns
808 -------
809 res : MonteCarloTestResult
810 An object with attributes:
812 statistic : float or ndarray
813 The test statistic of the observed `data`.
814 pvalue : float or ndarray
815 The p-value for the given alternative.
816 null_distribution : ndarray
817 The values of the test statistic generated under the null
818 hypothesis.
820 References
821 ----------
823 .. [1] B. Phipson and G. K. Smyth. "Permutation P-values Should Never Be
824 Zero: Calculating Exact P-values When Permutations Are Randomly Drawn."
825 Statistical Applications in Genetics and Molecular Biology 9.1 (2010).
827 Examples
828 --------
830 Suppose we wish to test whether a small sample has been drawn from a normal
831 distribution. We decide that we will use the skew of the sample as a
832 test statistic, and we will consider a p-value of 0.05 to be statistically
833 significant.
835 >>> import numpy as np
836 >>> from scipy import stats
837 >>> def statistic(x, axis):
838 ... return stats.skew(x, axis)
840 After collecting our data, we calculate the observed value of the test
841 statistic.
843 >>> rng = np.random.default_rng()
844 >>> x = stats.skewnorm.rvs(a=1, size=50, random_state=rng)
845 >>> statistic(x, axis=0)
846 0.12457412450240658
848 To determine the probability of observing such an extreme value of the
849 skewness by chance if the sample were drawn from the normal distribution,
850 we can perform a Monte Carlo hypothesis test. The test will draw many
851 samples at random from their normal distribution, calculate the skewness
852 of each sample, and compare our original skewness against this
853 distribution to determine an approximate p-value.
855 >>> from scipy.stats import monte_carlo_test
856 >>> # because our statistic is vectorized, we pass `vectorized=True`
857 >>> rvs = lambda size: stats.norm.rvs(size=size, random_state=rng)
858 >>> res = monte_carlo_test(x, rvs, statistic, vectorized=True)
859 >>> print(res.statistic)
860 0.12457412450240658
861 >>> print(res.pvalue)
862 0.7012
864 The probability of obtaining a test statistic less than or equal to the
865 observed value under the null hypothesis is ~70%. This is greater than
866 our chosen threshold of 5%, so we cannot consider this to be significant
867 evidence against the null hypothesis.
869 Note that this p-value essentially matches that of
870 `scipy.stats.skewtest`, which relies on an asymptotic distribution of a
871 test statistic based on the sample skewness.
873 >>> stats.skewtest(x).pvalue
874 0.6892046027110614
876 This asymptotic approximation is not valid for small sample sizes, but
877 `monte_carlo_test` can be used with samples of any size.
879 >>> x = stats.skewnorm.rvs(a=1, size=7, random_state=rng)
880 >>> # stats.skewtest(x) would produce an error due to small sample
881 >>> res = monte_carlo_test(x, rvs, statistic, vectorized=True)
883 The Monte Carlo distribution of the test statistic is provided for
884 further investigation.
886 >>> import matplotlib.pyplot as plt
887 >>> fig, ax = plt.subplots()
888 >>> ax.hist(res.null_distribution, bins=50)
889 >>> ax.set_title("Monte Carlo distribution of test statistic")
890 >>> ax.set_xlabel("Value of Statistic")
891 >>> ax.set_ylabel("Frequency")
892 >>> plt.show()
894 """
895 args = _monte_carlo_test_iv(data, rvs, statistic, vectorized,
896 n_resamples, batch, alternative, axis)
897 (data, rvs, statistic, vectorized,
898 n_resamples, batch, alternative, axis) = args
900 # Some statistics return plain floats; ensure they're at least np.float64
901 observed = np.asarray(statistic(*data, axis=-1))[()]
903 n_observations = [sample.shape[-1] for sample in data]
904 batch_nominal = batch or n_resamples
905 null_distribution = []
906 for k in range(0, n_resamples, batch_nominal):
907 batch_actual = min(batch_nominal, n_resamples - k)
908 resamples = [rvs_i(size=(batch_actual, n_observations_i))
909 for rvs_i, n_observations_i in zip(rvs, n_observations)]
910 null_distribution.append(statistic(*resamples, axis=-1))
911 null_distribution = np.concatenate(null_distribution)
912 null_distribution = null_distribution.reshape([-1] + [1]*observed.ndim)
914 def less(null_distribution, observed):
915 cmps = null_distribution <= observed
916 pvalues = (cmps.sum(axis=0) + 1) / (n_resamples + 1) # see [1]
917 return pvalues
919 def greater(null_distribution, observed):
920 cmps = null_distribution >= observed
921 pvalues = (cmps.sum(axis=0) + 1) / (n_resamples + 1) # see [1]
922 return pvalues
924 def two_sided(null_distribution, observed):
925 pvalues_less = less(null_distribution, observed)
926 pvalues_greater = greater(null_distribution, observed)
927 pvalues = np.minimum(pvalues_less, pvalues_greater) * 2
928 return pvalues
930 compare = {"less": less,
931 "greater": greater,
932 "two-sided": two_sided}
934 pvalues = compare[alternative](null_distribution, observed)
935 pvalues = np.clip(pvalues, 0, 1)
937 return MonteCarloTestResult(observed, pvalues, null_distribution)
940@dataclass
941class PermutationTestResult:
942 """Result object returned by `scipy.stats.permutation_test`.
944 Attributes
945 ----------
946 statistic : float or ndarray
947 The observed test statistic of the data.
948 pvalue : float or ndarray
949 The p-value for the given alternative.
950 null_distribution : ndarray
951 The values of the test statistic generated under the null
952 hypothesis.
953 """
954 statistic: float | np.ndarray
955 pvalue: float | np.ndarray
956 null_distribution: np.ndarray
959def _all_partitions_concatenated(ns):
960 """
961 Generate all partitions of indices of groups of given sizes, concatenated
963 `ns` is an iterable of ints.
964 """
965 def all_partitions(z, n):
966 for c in combinations(z, n):
967 x0 = set(c)
968 x1 = z - x0
969 yield [x0, x1]
971 def all_partitions_n(z, ns):
972 if len(ns) == 0:
973 yield [z]
974 return
975 for c in all_partitions(z, ns[0]):
976 for d in all_partitions_n(c[1], ns[1:]):
977 yield c[0:1] + d
979 z = set(range(np.sum(ns)))
980 for partitioning in all_partitions_n(z, ns[:]):
981 x = np.concatenate([list(partition)
982 for partition in partitioning]).astype(int)
983 yield x
986def _batch_generator(iterable, batch):
987 """A generator that yields batches of elements from an iterable"""
988 iterator = iter(iterable)
989 if batch <= 0:
990 raise ValueError("`batch` must be positive.")
991 z = [item for i, item in zip(range(batch), iterator)]
992 while z: # we don't want StopIteration without yielding an empty list
993 yield z
994 z = [item for i, item in zip(range(batch), iterator)]
997def _pairings_permutations_gen(n_permutations, n_samples, n_obs_sample, batch,
998 random_state):
999 # Returns a generator that yields arrays of size
1000 # `(batch, n_samples, n_obs_sample)`.
1001 # Each row is an independent permutation of indices 0 to `n_obs_sample`.
1002 batch = min(batch, n_permutations)
1004 if hasattr(random_state, 'permuted'):
1005 def batched_perm_generator():
1006 indices = np.arange(n_obs_sample)
1007 indices = np.tile(indices, (batch, n_samples, 1))
1008 for k in range(0, n_permutations, batch):
1009 batch_actual = min(batch, n_permutations-k)
1010 # Don't permute in place, otherwise results depend on `batch`
1011 permuted_indices = random_state.permuted(indices, axis=-1)
1012 yield permuted_indices[:batch_actual]
1013 else: # RandomState and early Generators don't have `permuted`
1014 def batched_perm_generator():
1015 for k in range(0, n_permutations, batch):
1016 batch_actual = min(batch, n_permutations-k)
1017 size = (batch_actual, n_samples, n_obs_sample)
1018 x = random_state.random(size=size)
1019 yield np.argsort(x, axis=-1)[:batch_actual]
1021 return batched_perm_generator()
1024def _calculate_null_both(data, statistic, n_permutations, batch,
1025 random_state=None):
1026 """
1027 Calculate null distribution for independent sample tests.
1028 """
1029 n_samples = len(data)
1031 # compute number of permutations
1032 # (distinct partitions of data into samples of these sizes)
1033 n_obs_i = [sample.shape[-1] for sample in data] # observations per sample
1034 n_obs_ic = np.cumsum(n_obs_i)
1035 n_obs = n_obs_ic[-1] # total number of observations
1036 n_max = np.prod([comb(n_obs_ic[i], n_obs_ic[i-1])
1037 for i in range(n_samples-1, 0, -1)])
1039 # perm_generator is an iterator that produces permutations of indices
1040 # from 0 to n_obs. We'll concatenate the samples, use these indices to
1041 # permute the data, then split the samples apart again.
1042 if n_permutations >= n_max:
1043 exact_test = True
1044 n_permutations = n_max
1045 perm_generator = _all_partitions_concatenated(n_obs_i)
1046 else:
1047 exact_test = False
1048 # Neither RandomState.permutation nor Generator.permutation
1049 # can permute axis-slices independently. If this feature is
1050 # added in the future, batches of the desired size should be
1051 # generated in a single call.
1052 perm_generator = (random_state.permutation(n_obs)
1053 for i in range(n_permutations))
1055 batch = batch or int(n_permutations)
1056 null_distribution = []
1058 # First, concatenate all the samples. In batches, permute samples with
1059 # indices produced by the `perm_generator`, split them into new samples of
1060 # the original sizes, compute the statistic for each batch, and add these
1061 # statistic values to the null distribution.
1062 data = np.concatenate(data, axis=-1)
1063 for indices in _batch_generator(perm_generator, batch=batch):
1064 indices = np.array(indices)
1066 # `indices` is 2D: each row is a permutation of the indices.
1067 # We use it to index `data` along its last axis, which corresponds
1068 # with observations.
1069 # After indexing, the second to last axis of `data_batch` corresponds
1070 # with permutations, and the last axis corresponds with observations.
1071 data_batch = data[..., indices]
1073 # Move the permutation axis to the front: we'll concatenate a list
1074 # of batched statistic values along this zeroth axis to form the
1075 # null distribution.
1076 data_batch = np.moveaxis(data_batch, -2, 0)
1077 data_batch = np.split(data_batch, n_obs_ic[:-1], axis=-1)
1078 null_distribution.append(statistic(*data_batch, axis=-1))
1079 null_distribution = np.concatenate(null_distribution, axis=0)
1081 return null_distribution, n_permutations, exact_test
1084def _calculate_null_pairings(data, statistic, n_permutations, batch,
1085 random_state=None):
1086 """
1087 Calculate null distribution for association tests.
1088 """
1089 n_samples = len(data)
1091 # compute number of permutations (factorial(n) permutations of each sample)
1092 n_obs_sample = data[0].shape[-1] # observations per sample; same for each
1093 n_max = factorial(n_obs_sample)**n_samples
1095 # `perm_generator` is an iterator that produces a list of permutations of
1096 # indices from 0 to n_obs_sample, one for each sample.
1097 if n_permutations >= n_max:
1098 exact_test = True
1099 n_permutations = n_max
1100 batch = batch or int(n_permutations)
1101 # cartesian product of the sets of all permutations of indices
1102 perm_generator = product(*(permutations(range(n_obs_sample))
1103 for i in range(n_samples)))
1104 batched_perm_generator = _batch_generator(perm_generator, batch=batch)
1105 else:
1106 exact_test = False
1107 batch = batch or int(n_permutations)
1108 # Separate random permutations of indices for each sample.
1109 # Again, it would be nice if RandomState/Generator.permutation
1110 # could permute each axis-slice separately.
1111 args = n_permutations, n_samples, n_obs_sample, batch, random_state
1112 batched_perm_generator = _pairings_permutations_gen(*args)
1114 null_distribution = []
1116 for indices in batched_perm_generator:
1117 indices = np.array(indices)
1119 # `indices` is 3D: the zeroth axis is for permutations, the next is
1120 # for samples, and the last is for observations. Swap the first two
1121 # to make the zeroth axis correspond with samples, as it does for
1122 # `data`.
1123 indices = np.swapaxes(indices, 0, 1)
1125 # When we're done, `data_batch` will be a list of length `n_samples`.
1126 # Each element will be a batch of random permutations of one sample.
1127 # The zeroth axis of each batch will correspond with permutations,
1128 # and the last will correspond with observations. (This makes it
1129 # easy to pass into `statistic`.)
1130 data_batch = [None]*n_samples
1131 for i in range(n_samples):
1132 data_batch[i] = data[i][..., indices[i]]
1133 data_batch[i] = np.moveaxis(data_batch[i], -2, 0)
1135 null_distribution.append(statistic(*data_batch, axis=-1))
1136 null_distribution = np.concatenate(null_distribution, axis=0)
1138 return null_distribution, n_permutations, exact_test
1141def _calculate_null_samples(data, statistic, n_permutations, batch,
1142 random_state=None):
1143 """
1144 Calculate null distribution for paired-sample tests.
1145 """
1146 n_samples = len(data)
1148 # By convention, the meaning of the "samples" permutations type for
1149 # data with only one sample is to flip the sign of the observations.
1150 # Achieve this by adding a second sample - the negative of the original.
1151 if n_samples == 1:
1152 data = [data[0], -data[0]]
1154 # The "samples" permutation strategy is the same as the "pairings"
1155 # strategy except the roles of samples and observations are flipped.
1156 # So swap these axes, then we'll use the function for the "pairings"
1157 # strategy to do all the work!
1158 data = np.swapaxes(data, 0, -1)
1160 # (Of course, the user's statistic doesn't know what we've done here,
1161 # so we need to pass it what it's expecting.)
1162 def statistic_wrapped(*data, axis):
1163 data = np.swapaxes(data, 0, -1)
1164 if n_samples == 1:
1165 data = data[0:1]
1166 return statistic(*data, axis=axis)
1168 return _calculate_null_pairings(data, statistic_wrapped, n_permutations,
1169 batch, random_state)
1172def _permutation_test_iv(data, statistic, permutation_type, vectorized,
1173 n_resamples, batch, alternative, axis, random_state):
1174 """Input validation for `permutation_test`."""
1176 axis_int = int(axis)
1177 if axis != axis_int:
1178 raise ValueError("`axis` must be an integer.")
1180 permutation_types = {'samples', 'pairings', 'independent'}
1181 permutation_type = permutation_type.lower()
1182 if permutation_type not in permutation_types:
1183 raise ValueError(f"`permutation_type` must be in {permutation_types}.")
1185 if vectorized not in {True, False, None}:
1186 raise ValueError("`vectorized` must be `True`, `False`, or `None`.")
1188 if vectorized is None:
1189 vectorized = 'axis' in inspect.signature(statistic).parameters
1191 if not vectorized:
1192 statistic = _vectorize_statistic(statistic)
1194 message = "`data` must be a tuple containing at least two samples"
1195 try:
1196 if len(data) < 2 and permutation_type == 'independent':
1197 raise ValueError(message)
1198 except TypeError:
1199 raise TypeError(message)
1201 data = _broadcast_arrays(data, axis)
1202 data_iv = []
1203 for sample in data:
1204 sample = np.atleast_1d(sample)
1205 if sample.shape[axis] <= 1:
1206 raise ValueError("each sample in `data` must contain two or more "
1207 "observations along `axis`.")
1208 sample = np.moveaxis(sample, axis_int, -1)
1209 data_iv.append(sample)
1211 n_resamples_int = (int(n_resamples) if not np.isinf(n_resamples)
1212 else np.inf)
1213 if n_resamples != n_resamples_int or n_resamples_int <= 0:
1214 raise ValueError("`n_resamples` must be a positive integer.")
1216 if batch is None:
1217 batch_iv = batch
1218 else:
1219 batch_iv = int(batch)
1220 if batch != batch_iv or batch_iv <= 0:
1221 raise ValueError("`batch` must be a positive integer or None.")
1223 alternatives = {'two-sided', 'greater', 'less'}
1224 alternative = alternative.lower()
1225 if alternative not in alternatives:
1226 raise ValueError(f"`alternative` must be in {alternatives}")
1228 random_state = check_random_state(random_state)
1230 return (data_iv, statistic, permutation_type, vectorized, n_resamples_int,
1231 batch_iv, alternative, axis_int, random_state)
1234def permutation_test(data, statistic, *, permutation_type='independent',
1235 vectorized=None, n_resamples=9999, batch=None,
1236 alternative="two-sided", axis=0, random_state=None):
1237 r"""
1238 Performs a permutation test of a given statistic on provided data.
1240 For independent sample statistics, the null hypothesis is that the data are
1241 randomly sampled from the same distribution.
1242 For paired sample statistics, two null hypothesis can be tested:
1243 that the data are paired at random or that the data are assigned to samples
1244 at random.
1246 Parameters
1247 ----------
1248 data : iterable of array-like
1249 Contains the samples, each of which is an array of observations.
1250 Dimensions of sample arrays must be compatible for broadcasting except
1251 along `axis`.
1252 statistic : callable
1253 Statistic for which the p-value of the hypothesis test is to be
1254 calculated. `statistic` must be a callable that accepts samples
1255 as separate arguments (e.g. ``statistic(*data)``) and returns the
1256 resulting statistic.
1257 If `vectorized` is set ``True``, `statistic` must also accept a keyword
1258 argument `axis` and be vectorized to compute the statistic along the
1259 provided `axis` of the sample arrays.
1260 permutation_type : {'independent', 'samples', 'pairings'}, optional
1261 The type of permutations to be performed, in accordance with the
1262 null hypothesis. The first two permutation types are for paired sample
1263 statistics, in which all samples contain the same number of
1264 observations and observations with corresponding indices along `axis`
1265 are considered to be paired; the third is for independent sample
1266 statistics.
1268 - ``'samples'`` : observations are assigned to different samples
1269 but remain paired with the same observations from other samples.
1270 This permutation type is appropriate for paired sample hypothesis
1271 tests such as the Wilcoxon signed-rank test and the paired t-test.
1272 - ``'pairings'`` : observations are paired with different observations,
1273 but they remain within the same sample. This permutation type is
1274 appropriate for association/correlation tests with statistics such
1275 as Spearman's :math:`\rho`, Kendall's :math:`\tau`, and Pearson's
1276 :math:`r`.
1277 - ``'independent'`` (default) : observations are assigned to different
1278 samples. Samples may contain different numbers of observations. This
1279 permutation type is appropriate for independent sample hypothesis
1280 tests such as the Mann-Whitney :math:`U` test and the independent
1281 sample t-test.
1283 Please see the Notes section below for more detailed descriptions
1284 of the permutation types.
1286 vectorized : bool, optional
1287 If `vectorized` is set ``False``, `statistic` will not be passed
1288 keyword argument `axis` and is expected to calculate the statistic
1289 only for 1D samples. If ``True``, `statistic` will be passed keyword
1290 argument `axis` and is expected to calculate the statistic along `axis`
1291 when passed an ND sample array. If ``None`` (default), `vectorized`
1292 will be set ``True`` if ``axis`` is a parameter of `statistic`. Use
1293 of a vectorized statistic typically reduces computation time.
1294 n_resamples : int or np.inf, default: 9999
1295 Number of random permutations (resamples) used to approximate the null
1296 distribution. If greater than or equal to the number of distinct
1297 permutations, the exact null distribution will be computed.
1298 Note that the number of distinct permutations grows very rapidly with
1299 the sizes of samples, so exact tests are feasible only for very small
1300 data sets.
1301 batch : int, optional
1302 The number of permutations to process in each call to `statistic`.
1303 Memory usage is O(`batch`*``n``), where ``n`` is the total size
1304 of all samples, regardless of the value of `vectorized`. Default is
1305 ``None``, in which case ``batch`` is the number of permutations.
1306 alternative : {'two-sided', 'less', 'greater'}, optional
1307 The alternative hypothesis for which the p-value is calculated.
1308 For each alternative, the p-value is defined for exact tests as
1309 follows.
1311 - ``'greater'`` : the percentage of the null distribution that is
1312 greater than or equal to the observed value of the test statistic.
1313 - ``'less'`` : the percentage of the null distribution that is
1314 less than or equal to the observed value of the test statistic.
1315 - ``'two-sided'`` (default) : twice the smaller of the p-values above.
1317 Note that p-values for randomized tests are calculated according to the
1318 conservative (over-estimated) approximation suggested in [2]_ and [3]_
1319 rather than the unbiased estimator suggested in [4]_. That is, when
1320 calculating the proportion of the randomized null distribution that is
1321 as extreme as the observed value of the test statistic, the values in
1322 the numerator and denominator are both increased by one. An
1323 interpretation of this adjustment is that the observed value of the
1324 test statistic is always included as an element of the randomized
1325 null distribution.
1326 The convention used for two-sided p-values is not universal;
1327 the observed test statistic and null distribution are returned in
1328 case a different definition is preferred.
1330 axis : int, default: 0
1331 The axis of the (broadcasted) samples over which to calculate the
1332 statistic. If samples have a different number of dimensions,
1333 singleton dimensions are prepended to samples with fewer dimensions
1334 before `axis` is considered.
1335 random_state : {None, int, `numpy.random.Generator`,
1336 `numpy.random.RandomState`}, optional
1338 Pseudorandom number generator state used to generate permutations.
1340 If `random_state` is ``None`` (default), the
1341 `numpy.random.RandomState` singleton is used.
1342 If `random_state` is an int, a new ``RandomState`` instance is used,
1343 seeded with `random_state`.
1344 If `random_state` is already a ``Generator`` or ``RandomState``
1345 instance then that instance is used.
1347 Returns
1348 -------
1349 res : PermutationTestResult
1350 An object with attributes:
1352 statistic : float or ndarray
1353 The observed test statistic of the data.
1354 pvalue : float or ndarray
1355 The p-value for the given alternative.
1356 null_distribution : ndarray
1357 The values of the test statistic generated under the null
1358 hypothesis.
1360 Notes
1361 -----
1363 The three types of permutation tests supported by this function are
1364 described below.
1366 **Unpaired statistics** (``permutation_type='independent'``):
1368 The null hypothesis associated with this permutation type is that all
1369 observations are sampled from the same underlying distribution and that
1370 they have been assigned to one of the samples at random.
1372 Suppose ``data`` contains two samples; e.g. ``a, b = data``.
1373 When ``1 < n_resamples < binom(n, k)``, where
1375 * ``k`` is the number of observations in ``a``,
1376 * ``n`` is the total number of observations in ``a`` and ``b``, and
1377 * ``binom(n, k)`` is the binomial coefficient (``n`` choose ``k``),
1379 the data are pooled (concatenated), randomly assigned to either the first
1380 or second sample, and the statistic is calculated. This process is
1381 performed repeatedly, `permutation` times, generating a distribution of the
1382 statistic under the null hypothesis. The statistic of the original
1383 data is compared to this distribution to determine the p-value.
1385 When ``n_resamples >= binom(n, k)``, an exact test is performed: the data
1386 are *partitioned* between the samples in each distinct way exactly once,
1387 and the exact null distribution is formed.
1388 Note that for a given partitioning of the data between the samples,
1389 only one ordering/permutation of the data *within* each sample is
1390 considered. For statistics that do not depend on the order of the data
1391 within samples, this dramatically reduces computational cost without
1392 affecting the shape of the null distribution (because the frequency/count
1393 of each value is affected by the same factor).
1395 For ``a = [a1, a2, a3, a4]`` and ``b = [b1, b2, b3]``, an example of this
1396 permutation type is ``x = [b3, a1, a2, b2]`` and ``y = [a4, b1, a3]``.
1397 Because only one ordering/permutation of the data *within* each sample
1398 is considered in an exact test, a resampling like ``x = [b3, a1, b2, a2]``
1399 and ``y = [a4, a3, b1]`` would *not* be considered distinct from the
1400 example above.
1402 ``permutation_type='independent'`` does not support one-sample statistics,
1403 but it can be applied to statistics with more than two samples. In this
1404 case, if ``n`` is an array of the number of observations within each
1405 sample, the number of distinct partitions is::
1407 np.prod([binom(sum(n[i:]), sum(n[i+1:])) for i in range(len(n)-1)])
1409 **Paired statistics, permute pairings** (``permutation_type='pairings'``):
1411 The null hypothesis associated with this permutation type is that
1412 observations within each sample are drawn from the same underlying
1413 distribution and that pairings with elements of other samples are
1414 assigned at random.
1416 Suppose ``data`` contains only one sample; e.g. ``a, = data``, and we
1417 wish to consider all possible pairings of elements of ``a`` with elements
1418 of a second sample, ``b``. Let ``n`` be the number of observations in
1419 ``a``, which must also equal the number of observations in ``b``.
1421 When ``1 < n_resamples < factorial(n)``, the elements of ``a`` are
1422 randomly permuted. The user-supplied statistic accepts one data argument,
1423 say ``a_perm``, and calculates the statistic considering ``a_perm`` and
1424 ``b``. This process is performed repeatedly, `permutation` times,
1425 generating a distribution of the statistic under the null hypothesis.
1426 The statistic of the original data is compared to this distribution to
1427 determine the p-value.
1429 When ``n_resamples >= factorial(n)``, an exact test is performed:
1430 ``a`` is permuted in each distinct way exactly once. Therefore, the
1431 `statistic` is computed for each unique pairing of samples between ``a``
1432 and ``b`` exactly once.
1434 For ``a = [a1, a2, a3]`` and ``b = [b1, b2, b3]``, an example of this
1435 permutation type is ``a_perm = [a3, a1, a2]`` while ``b`` is left
1436 in its original order.
1438 ``permutation_type='pairings'`` supports ``data`` containing any number
1439 of samples, each of which must contain the same number of observations.
1440 All samples provided in ``data`` are permuted *independently*. Therefore,
1441 if ``m`` is the number of samples and ``n`` is the number of observations
1442 within each sample, then the number of permutations in an exact test is::
1444 factorial(n)**m
1446 Note that if a two-sample statistic, for example, does not inherently
1447 depend on the order in which observations are provided - only on the
1448 *pairings* of observations - then only one of the two samples should be
1449 provided in ``data``. This dramatically reduces computational cost without
1450 affecting the shape of the null distribution (because the frequency/count
1451 of each value is affected by the same factor).
1453 **Paired statistics, permute samples** (``permutation_type='samples'``):
1455 The null hypothesis associated with this permutation type is that
1456 observations within each pair are drawn from the same underlying
1457 distribution and that the sample to which they are assigned is random.
1459 Suppose ``data`` contains two samples; e.g. ``a, b = data``.
1460 Let ``n`` be the number of observations in ``a``, which must also equal
1461 the number of observations in ``b``.
1463 When ``1 < n_resamples < 2**n``, the elements of ``a`` are ``b`` are
1464 randomly swapped between samples (maintaining their pairings) and the
1465 statistic is calculated. This process is performed repeatedly,
1466 `permutation` times, generating a distribution of the statistic under the
1467 null hypothesis. The statistic of the original data is compared to this
1468 distribution to determine the p-value.
1470 When ``n_resamples >= 2**n``, an exact test is performed: the observations
1471 are assigned to the two samples in each distinct way (while maintaining
1472 pairings) exactly once.
1474 For ``a = [a1, a2, a3]`` and ``b = [b1, b2, b3]``, an example of this
1475 permutation type is ``x = [b1, a2, b3]`` and ``y = [a1, b2, a3]``.
1477 ``permutation_type='samples'`` supports ``data`` containing any number
1478 of samples, each of which must contain the same number of observations.
1479 If ``data`` contains more than one sample, paired observations within
1480 ``data`` are exchanged between samples *independently*. Therefore, if ``m``
1481 is the number of samples and ``n`` is the number of observations within
1482 each sample, then the number of permutations in an exact test is::
1484 factorial(m)**n
1486 Several paired-sample statistical tests, such as the Wilcoxon signed rank
1487 test and paired-sample t-test, can be performed considering only the
1488 *difference* between two paired elements. Accordingly, if ``data`` contains
1489 only one sample, then the null distribution is formed by independently
1490 changing the *sign* of each observation.
1492 .. warning::
1493 The p-value is calculated by counting the elements of the null
1494 distribution that are as extreme or more extreme than the observed
1495 value of the statistic. Due to the use of finite precision arithmetic,
1496 some statistic functions return numerically distinct values when the
1497 theoretical values would be exactly equal. In some cases, this could
1498 lead to a large error in the calculated p-value. `permutation_test`
1499 guards against this by considering elements in the null distribution
1500 that are "close" (within a factor of ``1+1e-14``) to the observed
1501 value of the test statistic as equal to the observed value of the
1502 test statistic. However, the user is advised to inspect the null
1503 distribution to assess whether this method of comparison is
1504 appropriate, and if not, calculate the p-value manually. See example
1505 below.
1507 References
1508 ----------
1510 .. [1] R. A. Fisher. The Design of Experiments, 6th Ed (1951).
1511 .. [2] B. Phipson and G. K. Smyth. "Permutation P-values Should Never Be
1512 Zero: Calculating Exact P-values When Permutations Are Randomly Drawn."
1513 Statistical Applications in Genetics and Molecular Biology 9.1 (2010).
1514 .. [3] M. D. Ernst. "Permutation Methods: A Basis for Exact Inference".
1515 Statistical Science (2004).
1516 .. [4] B. Efron and R. J. Tibshirani. An Introduction to the Bootstrap
1517 (1993).
1519 Examples
1520 --------
1522 Suppose we wish to test whether two samples are drawn from the same
1523 distribution. Assume that the underlying distributions are unknown to us,
1524 and that before observing the data, we hypothesized that the mean of the
1525 first sample would be less than that of the second sample. We decide that
1526 we will use the difference between the sample means as a test statistic,
1527 and we will consider a p-value of 0.05 to be statistically significant.
1529 For efficiency, we write the function defining the test statistic in a
1530 vectorized fashion: the samples ``x`` and ``y`` can be ND arrays, and the
1531 statistic will be calculated for each axis-slice along `axis`.
1533 >>> import numpy as np
1534 >>> def statistic(x, y, axis):
1535 ... return np.mean(x, axis=axis) - np.mean(y, axis=axis)
1537 After collecting our data, we calculate the observed value of the test
1538 statistic.
1540 >>> from scipy.stats import norm
1541 >>> rng = np.random.default_rng()
1542 >>> x = norm.rvs(size=5, random_state=rng)
1543 >>> y = norm.rvs(size=6, loc = 3, random_state=rng)
1544 >>> statistic(x, y, 0)
1545 -3.5411688580987266
1547 Indeed, the test statistic is negative, suggesting that the true mean of
1548 the distribution underlying ``x`` is less than that of the distribution
1549 underlying ``y``. To determine the probability of this occuring by chance
1550 if the two samples were drawn from the same distribution, we perform
1551 a permutation test.
1553 >>> from scipy.stats import permutation_test
1554 >>> # because our statistic is vectorized, we pass `vectorized=True`
1555 >>> # `n_resamples=np.inf` indicates that an exact test is to be performed
1556 >>> res = permutation_test((x, y), statistic, vectorized=True,
1557 ... n_resamples=np.inf, alternative='less')
1558 >>> print(res.statistic)
1559 -3.5411688580987266
1560 >>> print(res.pvalue)
1561 0.004329004329004329
1563 The probability of obtaining a test statistic less than or equal to the
1564 observed value under the null hypothesis is 0.4329%. This is less than our
1565 chosen threshold of 5%, so we consider this to be significant evidence
1566 against the null hypothesis in favor of the alternative.
1568 Because the size of the samples above was small, `permutation_test` could
1569 perform an exact test. For larger samples, we resort to a randomized
1570 permutation test.
1572 >>> x = norm.rvs(size=100, random_state=rng)
1573 >>> y = norm.rvs(size=120, loc=0.3, random_state=rng)
1574 >>> res = permutation_test((x, y), statistic, n_resamples=100000,
1575 ... vectorized=True, alternative='less',
1576 ... random_state=rng)
1577 >>> print(res.statistic)
1578 -0.5230459671240913
1579 >>> print(res.pvalue)
1580 0.00016999830001699983
1582 The approximate probability of obtaining a test statistic less than or
1583 equal to the observed value under the null hypothesis is 0.0225%. This is
1584 again less than our chosen threshold of 5%, so again we have significant
1585 evidence to reject the null hypothesis in favor of the alternative.
1587 For large samples and number of permutations, the result is comparable to
1588 that of the corresponding asymptotic test, the independent sample t-test.
1590 >>> from scipy.stats import ttest_ind
1591 >>> res_asymptotic = ttest_ind(x, y, alternative='less')
1592 >>> print(res_asymptotic.pvalue)
1593 0.00012688101537979522
1595 The permutation distribution of the test statistic is provided for
1596 further investigation.
1598 >>> import matplotlib.pyplot as plt
1599 >>> plt.hist(res.null_distribution, bins=50)
1600 >>> plt.title("Permutation distribution of test statistic")
1601 >>> plt.xlabel("Value of Statistic")
1602 >>> plt.ylabel("Frequency")
1603 >>> plt.show()
1605 Inspection of the null distribution is essential if the statistic suffers
1606 from inaccuracy due to limited machine precision. Consider the following
1607 case:
1609 >>> from scipy.stats import pearsonr
1610 >>> x = [1, 2, 4, 3]
1611 >>> y = [2, 4, 6, 8]
1612 >>> def statistic(x, y):
1613 ... return pearsonr(x, y).statistic
1614 >>> res = permutation_test((x, y), statistic, vectorized=False,
1615 ... permutation_type='pairings',
1616 ... alternative='greater')
1617 >>> r, pvalue, null = res.statistic, res.pvalue, res.null_distribution
1619 In this case, some elements of the null distribution differ from the
1620 observed value of the correlation coefficient ``r`` due to numerical noise.
1621 We manually inspect the elements of the null distribution that are nearly
1622 the same as the observed value of the test statistic.
1624 >>> r
1625 0.8
1626 >>> unique = np.unique(null)
1627 >>> unique
1628 array([-1. , -0.8, -0.8, -0.6, -0.4, -0.2, -0.2, 0. , 0.2, 0.2, 0.4,
1629 0.6, 0.8, 0.8, 1. ]) # may vary
1630 >>> unique[np.isclose(r, unique)].tolist()
1631 [0.7999999999999999, 0.8]
1633 If `permutation_test` were to perform the comparison naively, the
1634 elements of the null distribution with value ``0.7999999999999999`` would
1635 not be considered as extreme or more extreme as the observed value of the
1636 statistic, so the calculated p-value would be too small.
1638 >>> incorrect_pvalue = np.count_nonzero(null >= r) / len(null)
1639 >>> incorrect_pvalue
1640 0.1111111111111111 # may vary
1642 Instead, `permutation_test` treats elements of the null distribution that
1643 are within ``max(1e-14, abs(r)*1e-14)`` of the observed value of the
1644 statistic ``r`` to be equal to ``r``.
1646 >>> correct_pvalue = np.count_nonzero(null >= r - 1e-14) / len(null)
1647 >>> correct_pvalue
1648 0.16666666666666666
1649 >>> res.pvalue == correct_pvalue
1650 True
1652 This method of comparison is expected to be accurate in most practical
1653 situations, but the user is advised to assess this by inspecting the
1654 elements of the null distribution that are close to the observed value
1655 of the statistic. Also, consider the use of statistics that can be
1656 calculated using exact arithmetic (e.g. integer statistics).
1658 """
1659 args = _permutation_test_iv(data, statistic, permutation_type, vectorized,
1660 n_resamples, batch, alternative, axis,
1661 random_state)
1662 (data, statistic, permutation_type, vectorized, n_resamples, batch,
1663 alternative, axis, random_state) = args
1665 observed = statistic(*data, axis=-1)
1667 null_calculators = {"pairings": _calculate_null_pairings,
1668 "samples": _calculate_null_samples,
1669 "independent": _calculate_null_both}
1670 null_calculator_args = (data, statistic, n_resamples,
1671 batch, random_state)
1672 calculate_null = null_calculators[permutation_type]
1673 null_distribution, n_resamples, exact_test = (
1674 calculate_null(*null_calculator_args))
1676 # See References [2] and [3]
1677 adjustment = 0 if exact_test else 1
1679 # relative tolerance for detecting numerically distinct but
1680 # theoretically equal values in the null distribution
1681 eps = 1e-14
1682 gamma = np.maximum(eps, np.abs(eps * observed))
1684 def less(null_distribution, observed):
1685 cmps = null_distribution <= observed + gamma
1686 pvalues = (cmps.sum(axis=0) + adjustment) / (n_resamples + adjustment)
1687 return pvalues
1689 def greater(null_distribution, observed):
1690 cmps = null_distribution >= observed - gamma
1691 pvalues = (cmps.sum(axis=0) + adjustment) / (n_resamples + adjustment)
1692 return pvalues
1694 def two_sided(null_distribution, observed):
1695 pvalues_less = less(null_distribution, observed)
1696 pvalues_greater = greater(null_distribution, observed)
1697 pvalues = np.minimum(pvalues_less, pvalues_greater) * 2
1698 return pvalues
1700 compare = {"less": less,
1701 "greater": greater,
1702 "two-sided": two_sided}
1704 pvalues = compare[alternative](null_distribution, observed)
1705 pvalues = np.clip(pvalues, 0, 1)
1707 return PermutationTestResult(observed, pvalues, null_distribution)
1710@dataclass
1711class ResamplingMethod:
1712 """Configuration information for a statistical resampling method.
1714 Instances of this class can be passed into the `method` parameter of some
1715 hypothesis test functions to perform a resampling or Monte Carlo version
1716 of the hypothesis test.
1718 Attributes
1719 ----------
1720 n_resamples : int
1721 The number of resamples to perform or Monte Carlo samples to draw.
1722 batch : int, optional
1723 The number of resamples to process in each vectorized call to
1724 the statistic. Batch sizes >>1 tend to be faster when the statistic
1725 is vectorized, but memory usage scales linearly with the batch size.
1726 Default is ``None``, which processes all resamples in a single batch.
1727 """
1728 n_resamples: int = 9999
1729 batch: int = None # type: ignore[assignment]
1732@dataclass
1733class MonteCarloMethod(ResamplingMethod):
1734 """Configuration information for a Monte Carlo hypothesis test.
1736 Instances of this class can be passed into the `method` parameter of some
1737 hypothesis test functions to perform a Monte Carlo version of the
1738 hypothesis tests.
1740 Attributes
1741 ----------
1742 n_resamples : int, optional
1743 The number of Monte Carlo samples to draw. Default is 9999.
1744 batch : int, optional
1745 The number of Monte Carlo samples to process in each vectorized call to
1746 the statistic. Batch sizes >>1 tend to be faster when the statistic
1747 is vectorized, but memory usage scales linearly with the batch size.
1748 Default is ``None``, which processes all samples in a single batch.
1749 rvs : callable or tuple of callables, optional
1750 A callable or sequence of callables that generates random variates
1751 under the null hypothesis. Each element of `rvs` must be a callable
1752 that accepts keyword argument ``size`` (e.g. ``rvs(size=(m, n))``) and
1753 returns an N-d array sample of that shape. If `rvs` is a sequence, the
1754 number of callables in `rvs` must match the number of samples passed
1755 to the hypothesis test in which the `MonteCarloMethod` is used. Default
1756 is ``None``, in which case the hypothesis test function chooses values
1757 to match the standard version of the hypothesis test. For example,
1758 the null hypothesis of `scipy.stats.pearsonr` is typically that the
1759 samples are drawn from the standard normal distribution, so
1760 ``rvs = (rng.normal, rng.normal)`` where
1761 ``rng = np.random.default_rng()``.
1762 """
1763 rvs: object = None
1765 def _asdict(self):
1766 # `dataclasses.asdict` deepcopies; we don't want that.
1767 return dict(n_resamples=self.n_resamples, batch=self.batch,
1768 rvs=self.rvs)
1771@dataclass
1772class PermutationMethod(ResamplingMethod):
1773 """Configuration information for a permutation hypothesis test.
1775 Instances of this class can be passed into the `method` parameter of some
1776 hypothesis test functions to perform a permutation version of the
1777 hypothesis tests.
1779 Attributes
1780 ----------
1781 n_resamples : int, optional
1782 The number of resamples to perform. Default is 9999.
1783 batch : int, optional
1784 The number of resamples to process in each vectorized call to
1785 the statistic. Batch sizes >>1 tend to be faster when the statistic
1786 is vectorized, but memory usage scales linearly with the batch size.
1787 Default is ``None``, which processes all resamples in a single batch.
1788 random_state : {None, int, `numpy.random.Generator`,
1789 `numpy.random.RandomState`}, optional
1791 Pseudorandom number generator state used to generate resamples.
1793 If `random_state` is already a ``Generator`` or ``RandomState``
1794 instance, then that instance is used.
1795 If `random_state` is an int, a new ``RandomState`` instance is used,
1796 seeded with `random_state`.
1797 If `random_state` is ``None`` (default), the
1798 `numpy.random.RandomState` singleton is used.
1799 """
1800 random_state: object = None
1802 def _asdict(self):
1803 # `dataclasses.asdict` deepcopies; we don't want that.
1804 return dict(n_resamples=self.n_resamples, batch=self.batch,
1805 random_state=self.random_state)
1808@dataclass
1809class BootstrapMethod(ResamplingMethod):
1810 """Configuration information for a bootstrap confidence interval.
1812 Instances of this class can be passed into the `method` parameter of some
1813 confidence interval methods to generate a bootstrap confidence interval.
1815 Attributes
1816 ----------
1817 n_resamples : int, optional
1818 The number of resamples to perform. Default is 9999.
1819 batch : int, optional
1820 The number of resamples to process in each vectorized call to
1821 the statistic. Batch sizes >>1 tend to be faster when the statistic
1822 is vectorized, but memory usage scales linearly with the batch size.
1823 Default is ``None``, which processes all resamples in a single batch.
1824 random_state : {None, int, `numpy.random.Generator`,
1825 `numpy.random.RandomState`}, optional
1827 Pseudorandom number generator state used to generate resamples.
1829 If `random_state` is already a ``Generator`` or ``RandomState``
1830 instance, then that instance is used.
1831 If `random_state` is an int, a new ``RandomState`` instance is used,
1832 seeded with `random_state`.
1833 If `random_state` is ``None`` (default), the
1834 `numpy.random.RandomState` singleton is used.
1836 method : {'bca', 'percentile', 'basic'}
1837 Whether to use the 'percentile' bootstrap ('percentile'), the 'basic'
1838 (AKA 'reverse') bootstrap ('basic'), or the bias-corrected and
1839 accelerated bootstrap ('BCa', default).
1840 """
1841 random_state: object = None
1842 method: str = 'BCa'
1844 def _asdict(self):
1845 # `dataclasses.asdict` deepcopies; we don't want that.
1846 return dict(n_resamples=self.n_resamples, batch=self.batch,
1847 random_state=self.random_state, method=self.method)