自助法：计算机时代的统计革命
自助重抽样是一种强大的统计技术，通过反复从观测数据中重抽样来估计任何统计量的抽样分布。它由 Bradley Efron 在 1979 年提出，无需依赖数学公式或分布假设就能进行复杂统计量的推断，彻底革新了统计推断的方式。
自助法背后的核心思想优雅而简洁：你的样本就是对总体的最佳估计。通过有放回地从样本中重抽样，你模拟了反复从总体中抽样会发生什么。这种方法对标准差尤为有价值，因为传统的置信区间公式假设正态性——而这个假设在实际中经常不成立。
自助法已成为现代数据科学不可或缺的工具，因为它适用于任何统计量（中位数、相关系数、回归系数、神经网络权重），且不对数据的底层分布做任何假设。
为什么用自助法计算标准差？
传统的标准差置信区间假设数据来自正态分布。当这个假设不成立时（这很常见），这些区间可能严重偏离实际。自助法提供了一种不依赖分布假设的替代方案。
传统方法失效时
自助法计算标准差的主要优势：
- 无分布假设：在正态、偏态或厚尾数据上表现同样好
- 小样本表现：当 n < 30 时通常比参数方法更准确
- 处理复杂统计量：同样的方法适用于截断标准差、MAD 或自定义的变异性指标
- 直观的可视化：自助分布让你看到正在发生什么，而不仅仅是最终数字
自助法步骤
自助法的算法非常直观。从你原始的 n 个观测值开始：
抽取自助样本
计算统计量
重复多次
分析分布
为什么要有放回？
需要多少自助样本？B = 1,000 通常足以进行粗略估计和假设检验。对于置信区间，B = 10,000 能提供稳定的百分位数。对于发表级别的 BCa 区间，建议 B = 15,000 以上。
自助置信区间方法
从自助样本构建置信区间有几种方法，各有优劣：
1. 百分位法（最简单）
最直观的方法：直接取自助分布的百分位数。
百分位法置信区间
对于 10,000 个自助样本，就是排序后的第 250 个和第 9,750 个值。简单，但当自助分布偏态时可能有偏。
2. 基本（枢轴）自助法
利用样本统计量与自助统计量之间的关系：
基本自助法置信区间
其中 θ̂ 是原始样本标准差。这种方法将百分位区间关于样本估计“镜像翻转”。
3. BCa 法（偏差校正加速法）
准确度的金标准。BCa 同时校正了自助分布的偏差和加速度（标准误差如何随参数值变化）。计算更复杂，但能提供二阶精确的区间。
| 方法 | 优点 | 缺点 |
|---|---|---|
| 百分位法 | 简单、直观 | 偏态数据时可能有偏 |
| 基本法 | 对称区间 | 可能产生负值 |
| BCa 法 | 最准确，尊重变换 | 计算量大 |
计算示例：非正态数据
考虑 15 次反应时间测量（毫秒）：245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312。这组数据呈右偏（部分反应非常慢）。
计算样本标准差
生成自助样本
计算自助标准差
求百分位数
构建 95% 置信区间
自助法的置信区间是不对称的（高端更宽），反映了数据的右偏特性。卡方分布的置信区间无法捕捉到这种不对称性。
Python 实现
包含多种置信区间方法的完整自助法实现：
import numpy as np
from scipy import stats
def bootstrap_sd_ci(data, n_bootstrap=10000, ci=0.95, method='percentile'):
"""
Bootstrap confidence interval for standard deviation.
Parameters:
-----------
data : array-like - Original sample
n_bootstrap : int - Number of bootstrap samples
ci : float - Confidence level (e.g., 0.95)
method : str - 'percentile', 'basic', or 'bca'
Returns:
--------
tuple : (lower_bound, upper_bound, bootstrap_sds)
"""
data = np.array(data)
n = len(data)
original_sd = np.std(data, ddof=1)
# Generate bootstrap samples and calculate SDs
bootstrap_sds = np.array([
np.std(np.random.choice(data, size=n, replace=True), ddof=1)
for _ in range(n_bootstrap)
])
alpha = 1 - ci
if method == 'percentile':
lower = np.percentile(bootstrap_sds, 100 * alpha/2)
upper = np.percentile(bootstrap_sds, 100 * (1 - alpha/2))
elif method == 'basic':
lower = 2*original_sd - np.percentile(bootstrap_sds, 100*(1-alpha/2))
upper = 2*original_sd - np.percentile(bootstrap_sds, 100*alpha/2)
elif method == 'bca':
# Bias correction
prop_less = np.mean(bootstrap_sds < original_sd)
z0 = stats.norm.ppf(prop_less)
# Acceleration (jackknife estimate)
jackknife_sds = np.array([
np.std(np.delete(data, i), ddof=1) for i in range(n)
])
jack_mean = jackknife_sds.mean()
a = np.sum((jack_mean - jackknife_sds)**3) / \
(6 * np.sum((jack_mean - jackknife_sds)**2)**1.5)
# Adjusted percentiles
z_alpha = stats.norm.ppf([alpha/2, 1-alpha/2])
adj_percentiles = stats.norm.cdf(
z0 + (z0 + z_alpha) / (1 - a*(z0 + z_alpha))
) * 100
lower = np.percentile(bootstrap_sds, adj_percentiles[0])
upper = np.percentile(bootstrap_sds, adj_percentiles[1])
return lower, upper, bootstrap_sds
# Example usage
response_times = [245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312]
for method in ['percentile', 'basic', 'bca']:
lower, upper, _ = bootstrap_sd_ci(response_times, method=method)
print(f"{method.upper():12s} 95% CI: [{lower:.1f}, {upper:.1f}]")