共计 6335 个字符，预计需要花费 16 分钟才能阅读完成。
“自举”（翻译自 bootstrap）这个词汇在多个畛域可能见到，它字面意思是提着靴子上的带子把本人提起来，这当然是不可能的，在机器学习畛域能够了解为原样本本身的数据再抽样得出新的样本及统计量，也有被翻译为自助法的。
Bayesian Bootstrap 是一个弱小的办法，它比其余的自举法更快，并且能够给出更严密的置信区间，并防止许多极其状况。在本文中咱们将具体地探讨这个简略但功能强大的过程。
自举
自举是通过对数据进行随机重采样和替换来计算估计量属性的过程，它首先由 Efron(1979)提出。这个过程非常简单，包含以下步骤：
假如一个 i.i.d. 样本 {Xᵢ}ᵢⁿ，并且咱们想用估计量 θ̂(X) 计算一个统计 θ。能够近似 θ̂的散布如下：
- 从样本{Xᵢ}ᵢⁿ中替换{X̃ᵢ}ᵢⁿ的 n 个察看样本。
- 计算估计量 θ̂-bootstrap(X̃)。
- 大量反复步骤 1 和步骤 2。
这样 θ̂-bootstrap 的散布很好地迫近了 θ̂的散布。
为什么 bootstrap 是无效的呢?
首先，它很容易实现。因为咱们只有反复做一件事件: 估算 θ，并且反复屡次就能够了。这其实也是自举的一个次要毛病：如果评估过程很慢，那么自举法的计算成本就会变得很高。
第二，自举不做散布假如。它只假如你的样本是总体的代表，察看后果是互相独立的。当察看后果彼此紧密联系时，比方在钻研社交网络或市场互动时，可能会违反此假如。
bootstrap 仅仅是加权吗?
当咱们从新抽样时，咱们所做的其实就是给咱们的察看值调配整数权重，这样它们的和就等于样本容量 n。这样的散布就是多项式散布。
咱们绘制大小为 10.000 的样本来看看多项式散布是什么样子的。从 src.utils 导入一些函数。
from src.utils import *
N = 10_000
np.random.seed(1)
bootstrap_weights = np.random.multinomial(N, np.ones(N)/N)
np.sum(bootstrap_weights)
#后果：10000首先，咱们确认权重之和是否的确等于 1000，或者说，咱们重采样生成了的是一个雷同大小的数据。
咱们当初能够画出权重的散布。
sns.countplot(bootstrap_weights, color='C0').set(title='Bootstrap Weights');正如咱们所看到的，大概 3600 次察看失去的权重为零，而一些察看失去的权重为 6。或者说大概 3600 个察看后果没有被从新采样，而一些察看后果被重采样多达 6 次。
这里可能就有一个问题：为什么不必间断权值来代替离散权值呢?
贝叶斯自举就是这个问题的答案。
贝叶斯自举
Bayesian bootstrap 是由 Rubin(1981)提出的，它基于一个非常简单的想法: 为什么不画一个更平滑的权重散布? 多项式散布的间断等价是狄利克雷散布。上面的图绘制了一次观测的多项和狄利克雷权重的概率分布(它们别离是泊松散布和伽马散布)。
贝叶斯自举的长处
第一个也是最直观的是，因为其间断的加权计划，它提供的估计值比一般的自举法更润滑。
此外间断加权计划阻止了极其状况的呈现（没有察看到的 0 权重）。例如在线性回归中，如果原始样本中没有共线性，则不会呈现共线性问题。
最初作为一种贝叶斯办法：估计量的预计散布能够解释为具备非信息先验的后验散布。
当初，让咱们画一个狄利克雷权重
bayesian_weights = np.random.dirichlet(alpha=np.ones(N), size=1)[0] * N
np.sum(bayesian_weights)
# 10000.000000000005权重的总和 (大概) 为 1，所以咱们必须将它们乘以 N。咱们能够画出权重的散布，当初为了失去了间断的权重，咱们必须近似散布。
sns.histplot(bayesian_weights, color='C1', kde=True).set(title='Dirichlet Weights');狄利克雷散布有一个参数 α，咱们对所有观测值都设置为 1。它是做什么的?
α 参数实质上决定被抽样的相对概率和绝对概率。减少所有观测值的 α 值能够缩小散布的偏斜，使所有观测值具备更类似的权重。对于 α→∞，所有的观测值得到雷同的权重。
那么咱们应该如何抉择 α 的值?Shao 和 Tu(1995)提出以下倡议。
The distribution of the random weight vector does not have to be restricted to the Diri(l, … , 1). Later investigations found that the weights having a scaled Diri(4, … ,4) distribution give better approximations (Tu and Zheng, 1987)
依据下面的倡议，让咱们来看看 α = 4 的狄利克雷散布和之前 α = 1 的狄利克雷散布是如何比拟的。
bayesian_weights2 = np.random.dirichlet(np.ones(N) * 4, 1)[0] * N
sns.histplot(bayesian_weights, color='C1')
sns.histplot(bayesian_weights2, color='C2').set(title='Comparing Dirichlet Weights');
plt.legend([r'$\alpha = 1$', r'$\alpha = 4$']);新的散布不那么歪斜了，并且更集中在平均值 1 左近。
例子
让咱们看几个例子，其中咱们比拟两个推断过程。
1、偏态散布的均值
首先，让咱们看一看最简略、最常见的估计量: 样本均值。首先咱们从帕累托散布中得出 100 个察看值。
np.random.seed(2)
X = pd.Series(np.random.pareto(2, 100))
sns.histplot(X);这种散布是十分歪斜的，几个观测值的值比平均值要高得多。
让咱们应用经典自举进行重采样，而后进行评估。
def classic_boot(df, estimator, seed=1):
df_boot = df.sample(n=len(df), replace=True, random_state=seed)
estimate = estimator(df_boot)
return estimate而后，让咱们应用一组随机权重的贝叶斯自举过程。
def bayes_boot(df, estimator, seed=1):
np.random.seed(seed)
w = np.random.dirichlet(np.ones(len(df)), 1)[0]
result = estimator(df, weights=w)
return result咱们这里应用 joblib 库并行化计算。
from joblib import Parallel, delayed
def bootstrap(boot_method, df, estimator, K):
r = Parallel(n_jobs=8)(delayed(boot_method)(df, estimator, seed=i) for i in range(K))
return r最初，让咱们写一个比拟后果的函数。
def compare_boot(df, boot1, boot2, estimator, title, K=1000):
s1 = bootstrap(boot1, df, estimator, K)
s2 = bootstrap(boot2, df, estimator, K)
df = pd.DataFrame({'Estimate': s1 + s2,
'Estimator': ['Classic']*K + ['Bayes']*K})
sns.histplot(data=df, x='Estimate', hue='Estimator')
plt.legend([f'Bayes: {np.mean(s2):.2f} ({np.std(s2):.2f})',
f'Classic: {np.mean(s1):.2f} ({np.std(s1):.2f})'])
plt.title(f'Bootstrap Estimates of {title}')当初开始比拟：
compare_boot(X, classic_boot, bayes_boot, np.average, 'Sample Mean')在这种状况下，这两个过程都给出了十分类似的后果。这两个散布十分靠近，而且估计量的预计平均值和标准偏差简直雷同，与咱们抉择的自举无关。
那么哪个过程更快呢？
import time
def compare_time(df, boot1, boot2, estimator, K=1000):
t1, t2 = np.zeros(K), np.zeros(K)
for k in range(K):
# Classic bootstrap
start = time.time()
boot1(df, estimator)
t1[k] = time.time() - start
# Bayesian bootstrap
start = time.time()
boot2(df, estimator)
t2[k] = time.time() - start
print(f"Bayes wins {np.mean(t1 > t2)*100}% of the time (by {np.mean((t1 - t2)/t1*100):.2f}%)")后果如下：
Bayes wins 97.8% of the time (by 66.46%)贝叶斯自举快了很多。
2、没有权重怎么办？也没问题
如果咱们有一个不承受权重的估计量，例如中位数？咱们能够进行两级抽样：咱们采样权重，而后依据权重采样观测值。
def twolv_boot(df, estimator, seed=1):
np.random.seed(seed)
w = np.random.dirichlet(np.ones(len(df))*4, 1)[0]
df_boot = df.sample(n=len(df)*10, replace=True, weights=w, random_state=seed)
result = estimator(df_boot)
return result看看后果：
np.random.seed(1)
X = pd.Series(np.random.normal(0, 10, 1000))
compare_boot(X, classic_boot, twolv_boot, np.median, "Sample Median")能够看到，贝叶斯自举也比个别的自举更准确，因为 α = 4 时的权重散布更密集。
3、逻辑回归和离群值
当初，看看个别的自举过程可能产生问题的案例。假如咱们察看到一个正态分布的特色 X 和二元后果 y。咱们对两个变量之间的关系进行钻研。
N = 100
np.random.seed(1)
x = np.random.normal(0, 1, N)
y = np.rint(np.random.normal(x, 1, N) > 2)
df = pd.DataFrame({'x': x, 'y': y})
df.head()在该样本中，咱们仅在 100 个察看后果中的到真值后果。因为后果是二进制的，因而能够应用逻辑回归模型。
smf.logit('y ~ x', data=df).fit(disp=False).summary().tables[1]咱们失去一个 -23 点的估计值，置信区间十分严密的。
咱们能自举估计量的散布吗? 上面计算 1000 个自举样本的逻辑回归系数。
estimate_logit = lambda df: smf.logit('y ~ x', data=df).fit(disp=False).params[1]
for i in range(1000):
try:
classic_boot(df, estimate_logit, seed=i)
except Exception as e:
print(f'Error for bootstrap number {i}: {e}')会失去以下谬误：
Error for bootstrap number 92: Perfect separation detected, results not available
Error for bootstrap number 521: Perfect separation detected, results not available
Error for bootstrap number 545: Perfect separation detected, results not available
Error for bootstrap number 721: Perfect separation detected, results not available
Error for bootstrap number 835: Perfect separation detected, results not available对于 1000 个样本中的 5 个，咱们无奈计算估计值。然而这种状况是不会产生在贝叶斯自举过程中的。
因为对于贝叶斯自举能够疏忽这些察看后果。
4、应用 Treated Units 进行回归
假如咱们察看到二元特色 X 和间断的后果 y。咱们再次钻研两个变量之间的关系。
N = 100
np.random.seed(1)
x = np.random.binomial(1, 5/N, N)
y = np.random.normal(1 + 2*x, 1, N)
df = pd.DataFrame({'x': x, 'y': y})
df.head()让咱们比拟 X 上 Y 回归系数的两个估计量。
estimate_beta = lambda df, **kwargs: smf.wls('y ~ x', data=df, **kwargs).fit().params[1]
compare_boot(df, classic_boot, bayes_boot, estimate_beta, "beta")个别的自举过程预计的差别比贝叶斯自举大 50%。为什么？如果咱们更认真地查看就会发现在将近 20 个从新采样的样本中，会失去一个十分不寻常的预计！
这个问题的起因是在某些样本中，咱们可能没有任何察看后果 x = 1。因而在这些重采样的样本中，预计的系数为零。贝叶斯自举程序不会产生这种状况，因为它不会放弃任何察看（所有察看后果总是会蕴含所有的后果）。
并且这里的一个重要的问题是咱们没有收到任何谬误音讯或正告，这样咱们很容易会漠视这个问题！
总结
在本文中咱们介绍了贝叶斯自举法，它的要害的想法是，每当咱们的估计量以加权估计量示意时，自举过程就等于用多项式权重随机加权。贝叶斯自举等同于用狄利克雷权重加权，这是多项式散布的间断等效物。具备间断的权重防止了极其的样本，并且能够生成估计量的平滑散布。
本文参考
[1] B. Efron Bootstrap Methods: Another Look at the Jackknife (1979), The Annals of Statistics.
[2] D. Rubin, The Bayesian Bootstrap (1981), The Annals of Statistics.
[3] A. Lo, A Large Sample Study of the Bayesian Bootstrap (1987), The Annals of Statistics.
[4] J. Shao, D. Tu, Jacknife and Bootstrap (1995), Springer.
代码：
https://avoid.overfit.cn/post/8a0dab8d9ae34bb3926966915f639163
作者：Matteo Courthoud
