bootstrap#
- scipy.stats.bootstrap(data, statistic, *, n_resamples=9999, batch=None, vectorized=None, paired=False, axis=0, confidence_level=0.95, alternative='two-sided', method='BCa', bootstrap_result=None, random_state=None)[源代码][源代码]#
计算统计量的双侧自助置信区间。
当 method 为
'percentile'且 alternative 为'two-sided'时，将根据以下步骤计算自举置信区间。重采样数据：对于 data 中的每个样本和 n_resamples 中的每个样本，从原始样本中随机抽取一个与原始样本大小相同的样本（允许重复）。
计算统计量的bootstrap分布：对于每一组重采样，计算检验统计量。
确定置信区间：找到自举分布的区间，即
关于中位数对称且
包含重采样统计值的 置信水平。
虽然
'percentile'方法是最直观的，但在实践中很少使用。还有两种更常见的方法，``’basic’``（’反百分位数’）和 ``’BCa’``（’偏差校正和加速’）；它们在步骤 3 的执行方式上有所不同。如果 data 中的样本是从各自的分布中随机抽取的 \(n\) 次，
bootstrap返回的置信区间将包含这些分布的统计量的真实值，大约 confidence_level\(\, \times \, n\) 次。- 参数:
- 数据类数组的序列
data 的每个元素都是一个包含从基础分布中得到的标量观测值的样本。data 的元素必须能够广播到相同的形状（axis 指定的维度可能除外）。
- 统计可调用
要计算置信区间的统计量。statistic 必须是一个可调用对象，它接受
len(data)个样本作为单独的参数并返回结果统计量。如果 vectorized 设置为True，statistic 还必须接受一个关键字参数 axis，并且能够向量化以沿着提供的 axis 计算统计量。- n_resamples : int, 默认值:
9999int, 默认值: 用于形成统计量的bootstrap分布的重采样次数。
- 批处理int, 可选
在每次向量化调用 statistic 时处理的重新采样次数。内存使用量为 O( batch *
n)，其中n是样本大小。默认值为None，在这种情况下，batch = n_resamples``（或 ``batch = max(n_resamples, n)用于method='BCa'）。- 矢量化bool, 可选
如果 vectorized 设置为
False，statistic 将不会传递关键字参数 axis，并且预期仅计算 1D 样本的统计量。如果为True，statistic 将传递关键字参数 axis，并且预期在传递 ND 样本数组时沿 axis 计算统计量。如果为None``（默认），如果 `statistic` 的参数中包含 `axis`，则 `vectorized` 将设置为 ``True。使用向量化统计量通常会减少计算时间。- paired : bool, 默认:
Falsebool, 默认值: 统计量是否将 data 中样本的对应元素视为配对。
- 轴 : int, 默认值:
0int, 默认值: 在 data 中样本的轴，沿着该轴计算 statistic。
- confidence_level : float, 默认值:
0.95浮点数，默认值： 置信区间的置信水平。
- alternative : {‘two-sided’, ‘less’, ‘greater’}, 默认:
'two-sided'{‘双侧’, ‘小于’, ‘大于’}，默认： 选择
'two-sided'``（默认）以获得双侧置信区间，选择 ``'less'以获得下限为-np.inf的单侧置信区间，选择'greater'以获得上限为np.inf的单侧置信区间。单侧置信区间的另一边界与双侧置信区间相同，但 confidence_level 距离1.0的距离是双侧置信区间的一半；例如，95%'less'置信区间的上限与90%'two-sided'置信区间的上限相同。- 方法 : {‘百分位数’, ‘基本’, ‘bca’}, 默认:
'BCa'{‘百分位数’, ‘基本’, ‘BCa’}, 默认: 是否返回 ‘百分位数’ 自举置信区间 (
'percentile')，’基本’ (又名 ‘反向’) 自举置信区间 ('basic')，或偏差校正和加速自举置信区间 ('BCa')。- bootstrap_resultBootstrapResult, 可选
提供先前调用
bootstrap返回的结果对象，以在新引导分布中包含先前的引导分布。例如，这可以用于更改 confidence_level、更改 method，或查看在不重复计算的情况下执行额外重采样的效果。- random_state{None, int,}
用于生成重采样的伪随机数生成器状态。
如果 random_state 是
None（或 np.random），则使用numpy.random.RandomState单例。如果 random_state 是整数，则使用一个新的RandomState实例，并以 random_state 为种子。如果 random_state 已经是Generator或RandomState实例，则使用该实例。
- 返回:
- resBootstrapResult
一个带有属性的对象：
- 置信区间置信区间
bootstrap 置信区间作为
collections.namedtuple的一个实例，具有 low 和 high 属性。- bootstrap_distributionndarray
bootstrap 分布，即每个重采样的 statistic 值。最后一个维度对应于重采样（例如
res.bootstrap_distribution.shape[-1] == n_resamples）。- 标准误差浮点数或ndarray
bootstrap 标准误差，即 bootstrap 分布的样本标准差。
- 警告:
DegenerateDataWarning当
method='BCa'且 bootstrap 分布是退化的（例如所有元素都相同）时生成。
注释
对于
method='BCa'，如果 bootstrap 分布是退化的（例如所有元素都相同），置信区间的元素可能是 NaN。在这种情况下，考虑使用另一种 method 或检查 data 以确定其他分析可能更合适（例如所有观测值都相同）。参考文献
[1]B. Efron and R. J. Tibshirani, An Introduction to the Bootstrap, Chapman & Hall/CRC, Boca Raton, FL, USA (1993)
[2]Nathaniel E. Helwig, “Bootstrap置信区间”, http://users.stat.umn.edu/~helwig/notes/bootci-Notes.pdf
[3]Bootstrapping (统计学), 维基百科, https://en.wikipedia.org/wiki/Bootstrapping_%28statistics%29
示例
假设我们从某个未知分布中采样了数据。
>>> import numpy as np >>> rng = np.random.default_rng() >>> from scipy.stats import norm >>> dist = norm(loc=2, scale=4) # our "unknown" distribution >>> data = dist.rvs(size=100, random_state=rng)
我们对分布的标准差感兴趣。
>>> std_true = dist.std() # the true value of the statistic >>> print(std_true) 4.0 >>> std_sample = np.std(data) # the sample statistic >>> print(std_sample) 3.9460644295563863
引导法用于近似我们在从未知分布中反复抽样并每次计算样本统计量时所期望的变异性。它通过反复从*原始样本*中进行有放回的抽样并计算每个重抽样样本的统计量来实现这一点。这导致了一个统计量的“引导分布”。
>>> import matplotlib.pyplot as plt >>> from scipy.stats import bootstrap >>> data = (data,) # samples must be in a sequence >>> res = bootstrap(data, np.std, confidence_level=0.9, ... random_state=rng) >>> fig, ax = plt.subplots() >>> ax.hist(res.bootstrap_distribution, bins=25) >>> ax.set_title('Bootstrap Distribution') >>> ax.set_xlabel('statistic value') >>> ax.set_ylabel('frequency') >>> plt.show()
标准误差量化了这种变异性。它是通过计算自举分布的标准偏差来计算的。
>>> res.standard_error 0.24427002125829136 >>> res.standard_error == np.std(res.bootstrap_distribution, ddof=1) True
统计量的bootstrap分布通常近似于标准误差等于其标准差的正态分布。
>>> x = np.linspace(3, 5) >>> pdf = norm.pdf(x, loc=std_sample, scale=res.standard_error) >>> fig, ax = plt.subplots() >>> ax.hist(res.bootstrap_distribution, bins=25, density=True) >>> ax.plot(x, pdf) >>> ax.set_title('Normal Approximation of the Bootstrap Distribution') >>> ax.set_xlabel('statistic value') >>> ax.set_ylabel('pdf') >>> plt.show()
这表明我们可以基于该正态分布的分位数，为统计量构建一个90%的置信区间。
>>> norm.interval(0.9, loc=std_sample, scale=res.standard_error) (3.5442759991341726, 4.3478528599786)
由于中心极限定理，这种正态近似对于各种统计量和样本基础分布是准确的；然而，这种近似并非在所有情况下都可靠。因为
bootstrap旨在处理任意基础分布和统计量，它使用更高级的技术来生成准确的置信区间。>>> print(res.confidence_interval) ConfidenceInterval(low=3.57655333533867, high=4.382043696342881)
如果我们从原始分布中抽样100次，并为每次抽样形成一个bootstrap置信区间，置信区间大约有90%的时间包含统计量的真实值。
>>> n_trials = 100 >>> ci_contains_true_std = 0 >>> for i in range(n_trials): ... data = (dist.rvs(size=100, random_state=rng),) ... res = bootstrap(data, np.std, confidence_level=0.9, ... n_resamples=999, random_state=rng) ... ci = res.confidence_interval ... if ci[0] < std_true < ci[1]: ... ci_contains_true_std += 1 >>> print(ci_contains_true_std) 88
我们也可以一次性确定所有100个样本的置信区间，而不是编写一个循环。
>>> data = (dist.rvs(size=(n_trials, 100), random_state=rng),) >>> res = bootstrap(data, np.std, axis=-1, confidence_level=0.9, ... n_resamples=999, random_state=rng) >>> ci_l, ci_u = res.confidence_interval
在这里，ci_l 和 ci_u 包含了每个
n_trials = 100样本的置信区间。>>> print(ci_l[:5]) [3.86401283 3.33304394 3.52474647 3.54160981 3.80569252] >>> print(ci_u[:5]) [4.80217409 4.18143252 4.39734707 4.37549713 4.72843584]
再次，大约90% 包含真实值，
std_true = 4。>>> print(np.sum((ci_l < std_true) & (std_true < ci_u))) 93
bootstrap也可以用于估计多样本统计的置信区间。例如，要获取均值差异的置信区间，我们编写一个接受两个样本参数并仅返回统计量的函数。使用axis参数确保所有均值计算都在一次向量化调用中执行，这比在 Python 中循环遍历重采样对更快。>>> def my_statistic(sample1, sample2, axis=-1): ... mean1 = np.mean(sample1, axis=axis) ... mean2 = np.mean(sample2, axis=axis) ... return mean1 - mean2
在这里，我们使用默认95%置信水平的’百分位数’方法。
>>> sample1 = norm.rvs(scale=1, size=100, random_state=rng) >>> sample2 = norm.rvs(scale=2, size=100, random_state=rng) >>> data = (sample1, sample2) >>> res = bootstrap(data, my_statistic, method='basic', random_state=rng) >>> print(my_statistic(sample1, sample2)) 0.16661030792089523 >>> print(res.confidence_interval) ConfidenceInterval(low=-0.29087973240818693, high=0.6371338699912273)
标准误差的bootstrap估计也是可用的。
>>> print(res.standard_error) 0.238323948262459
配对样本统计也适用。例如，考虑皮尔逊相关系数。
>>> from scipy.stats import pearsonr >>> n = 100 >>> x = np.linspace(0, 10, n) >>> y = x + rng.uniform(size=n) >>> print(pearsonr(x, y)[0]) # element 0 is the statistic 0.9954306665125647
我们封装
pearsonr以便它只返回统计量，确保我们使用 axis 参数，因为它可用。>>> def my_statistic(x, y, axis=-1): ... return pearsonr(x, y, axis=axis)[0]
我们使用
paired=True调用bootstrap。>>> res = bootstrap((x, y), my_statistic, paired=True, random_state=rng) >>> print(res.confidence_interval) ConfidenceInterval(low=0.9941504301315878, high=0.996377412215445)
结果对象可以传递回
bootstrap以执行额外的重采样：>>> len(res.bootstrap_distribution) 9999 >>> res = bootstrap((x, y), my_statistic, paired=True, ... n_resamples=1000, random_state=rng, ... bootstrap_result=res) >>> len(res.bootstrap_distribution) 10999
或更改置信区间选项：
>>> res2 = bootstrap((x, y), my_statistic, paired=True, ... n_resamples=0, random_state=rng, bootstrap_result=res, ... method='percentile', confidence_level=0.9) >>> np.testing.assert_equal(res2.bootstrap_distribution, ... res.bootstrap_distribution) >>> res.confidence_interval ConfidenceInterval(low=0.9941574828235082, high=0.9963781698210212)
无需重复计算原始的 bootstrap 分布。