scipy.stats.monte_carlo_test — SciPy v1.13.0 Manual (2024)

Perform a Monte Carlo hypothesis test.

data contains a sample or a sequence of one or more samples. rvsspecifies the distribution(s) of the sample(s) in data under the nullhypothesis. The value of statistic for the given data is comparedagainst a Monte Carlo null distribution: the value of the statistic foreach of n_resamples sets of samples generated using rvs. This givesthe p-value, the probability of observing such an extreme value of thetest statistic under the null hypothesis.

Parameters:

dataarray-like or sequence of array-like

An array or sequence of arrays of observations.

rvscallable or tuple of callables

A callable or sequence of callables that generates random variatesunder the null hypothesis. Each element of rvs must be a callablethat accepts keyword argument size (e.g. rvs(size=(m, n))) andreturns an N-d array sample of that shape. If rvs is a sequence, thenumber of callables in rvs must match the number of samples indata, i.e. len(rvs) == len(data). If rvs is a single callable,data is treated as a single sample.

statisticcallable

Statistic for which the p-value of the hypothesis test is to becalculated. statistic must be a callable that accepts a sample(e.g. statistic(sample)) or len(rvs) separate samples (e.g.statistic(samples1, sample2) if rvs contains two callables anddata contains two samples) and returns the resulting statistic.If vectorized is set True, statistic must also accept a keywordargument axis and be vectorized to compute the statistic along theprovided axis of the samples in data.

vectorizedbool, optional

If vectorized is set False, statistic will not be passedkeyword argument axis and is expected to calculate the statisticonly for 1D samples. If True, statistic will be passed keywordargument axis and is expected to calculate the statistic along axiswhen passed ND sample arrays. If None (default), vectorizedwill be set True if axis is a parameter of statistic. Use ofa vectorized statistic typically reduces computation time.

n_resamplesint, default: 9999

Number of samples drawn from each of the callables of rvs.Equivalently, the number statistic values under the null hypothesisused as the Monte Carlo null distribution.

batchint, optional

The number of Monte Carlo samples to process in each call tostatistic. Memory usage is O( batch * sample.size[axis] ). Defaultis None, in which case batch equals n_resamples.

alternative{‘two-sided’, ‘less’, ‘greater’}

The alternative hypothesis for which the p-value is calculated.For each alternative, the p-value is defined as follows.

• 'greater' : the percentage of the null distribution that isgreater than or equal to the observed value of the test statistic.

• 'less' : the percentage of the null distribution that isless than or equal to the observed value of the test statistic.

• 'two-sided' : twice the smaller of the p-values above.

  See Also

  Monte Carlo Methods (Monographs on Applied Probability & Statistics)Tennis Best Bets: Monte Carlo Masters picks and predictions for Monday, April 8Making Sense of Monte Carlo

axisint, default: 0

The axis of data (or each sample within data) over which tocalculate the statistic.

Returns:

resMonteCarloTestResult

An object with attributes:

statisticfloat or ndarray

The test statistic of the observed data.

pvaluefloat or ndarray

The p-value for the given alternative.

null_distributionndarray

The values of the test statistic generated under the nullhypothesis.

Warning

The p-value is calculated by counting the elements of the nulldistribution that are as extreme or more extreme than the observedvalue of the statistic. Due to the use of finite precision arithmetic,some statistic functions return numerically distinct values when thetheoretical values would be exactly equal. In some cases, this couldlead to a large error in the calculated p-value. monte_carlo_testguards against this by considering elements in the null distributionthat are “close” (within a relative tolerance of 100 times thefloating point epsilon of inexact dtypes) to the observedvalue of the test statistic as equal to the observed value of thetest statistic. However, the user is advised to inspect the nulldistribution to assess whether this method of comparison isappropriate, and if not, calculate the p-value manually.

References

Examples

Suppose we wish to test whether a small sample has been drawn from a normaldistribution. We decide that we will use the skew of the sample as atest statistic, and we will consider a p-value of 0.05 to be statisticallysignificant.

>>> import numpy as np>>> from scipy import stats>>> def statistic(x, axis):...  return stats.skew(x, axis)

After collecting our data, we calculate the observed value of the teststatistic.

>>> rng = np.random.default_rng()>>> x = stats.skewnorm.rvs(a=1, size=50, random_state=rng)>>> statistic(x, axis=0)0.12457412450240658

To determine the probability of observing such an extreme value of theskewness by chance if the sample were drawn from the normal distribution,we can perform a Monte Carlo hypothesis test. The test will draw manysamples at random from their normal distribution, calculate the skewnessof each sample, and compare our original skewness against thisdistribution to determine an approximate p-value.

>>> from scipy.stats import monte_carlo_test>>> # because our statistic is vectorized, we pass `vectorized=True`>>> rvs = lambda size: stats.norm.rvs(size=size, random_state=rng)>>> res = monte_carlo_test(x, rvs, statistic, vectorized=True)>>> print(res.statistic)0.12457412450240658>>> print(res.pvalue)0.7012

The probability of obtaining a test statistic less than or equal to theobserved value under the null hypothesis is ~70%. This is greater thanour chosen threshold of 5%, so we cannot consider this to be significantevidence against the null hypothesis.

Note that this p-value essentially matches that ofscipy.stats.skewtest, which relies on an asymptotic distribution of atest statistic based on the sample skewness.

>>> stats.skewtest(x).pvalue0.6892046027110614

This asymptotic approximation is not valid for small sample sizes, butmonte_carlo_test can be used with samples of any size.

>>> x = stats.skewnorm.rvs(a=1, size=7, random_state=rng)>>> # stats.skewtest(x) would produce an error due to small sample>>> res = monte_carlo_test(x, rvs, statistic, vectorized=True)

The Monte Carlo distribution of the test statistic is provided forfurther investigation.

>>> import matplotlib.pyplot as plt>>> fig, ax = plt.subplots()>>> ax.hist(res.null_distribution, bins=50)>>> ax.set_title("Monte Carlo distribution of test statistic")>>> ax.set_xlabel("Value of Statistic")>>> ax.set_ylabel("Frequency")>>> plt.show()