r语言Bootstrap自助法重采样构建统计量T抽样分布近似值可视化|代码分享
r语言Bootstrap自助法重采样构建统计量T抽样分布近似值可视化|代码分享
全文链接：https://tecdat.cn/?p=33939
统计量T是数据的一个函数，不依赖于任何未知参数（即我们可以根据数据计算得到它）（点击文末“阅读原文”获取完整代码数据）。
相关视频
这意味着给定数据值x1,x2,⋯,xn，统计量T就是一个"数字"。然而，在观察到数据之前，"数据"是随机变量X1,X2,⋯,Xn，而我们的统计量T作为随机变量的函数，也是一个随机变量。T的分布被称为"抽样分布"。
例如，如果我们有以下数据：
感兴趣的统计量是X¯=1/n∑ni=1Xi，我们知道
这就是X¯的抽样分布。统计量的抽样分布并不总是容易找到。让我们考虑两种抽样分布更难以通过解析方法找到的情况。
情况1
假设我们有来自一个倾斜分布的40个数据点。下面给出了数据的直方图。
我们首先计算样本均值和样本标准差。
#数据的平均值
mean(x)
#数据的方差
var(x)
中心极限定理告诉我们，当n很大时，样本均值将服从正态分布。但是这里有一个重要问题：我们怎么知道n是否足够大呢？尽管数据倾斜严重，我们应该相信CLT的近似吗？
情况2
考虑一组新的200个数据点（我们将这些数据称为yi）。直方图如下所示
点击标题查阅往期内容
01
02
03
04
这个分布是右偏的还是对称的？很难说。回想一下，分布的总体偏度定义为
这个参数的一个简单估计量（统计量）是下面给出的"样本偏度"
其中y¯和s是数据的样本均值和标准差。那么问题来了，γ^的抽样分布是什么？这个分布肯定是非常难以解析计算的。
自助法
自助法最初由Bradley Efron在1979年发表，被称为一种重新采样技术。Bootstrap（至少这个版本）被称为一种非参数方法，因为它不需要我们对数据做出任何特定分布的假设，这是一个有用的特性。
基本思想是，如果样本数据准确反映了总体，我们可以"重新采样"数据并构建统计量T的抽样分布的近似值。这个近似值有时被称为T的"Bootstrap分布"。需要记住的是，像大多数统计方法一样，当样本量非常小时，Bootstrap可能会失败。
算法其实相当简单，步骤如下：
- 通过从原始数据中（有放回地）抽样，创建一个“新”数据集，直到你有一个大小为 n 的新数据集。
- 计算这个新数据集的检验统计量，并将其称为 T1。
- 重复步骤 1 和 2 多次（比如说 B 次），这样你就得到了一系列的估计值 T1，T2，⋯，TB。这是对 T 的抽样分布的数值近似。
情况1 - 使用自助法
在这个例子中，我们可以使用自助法来近似样本均值 X¯ 的抽样分布。如果自助法的分布看起来近似正态分布，那么我们可以合理地认为中心极限定理（CLT）会给出一个不错的近似结果。如果不是，我们应该怀疑是否能够信任CLT对于这个数据的适用性。
B <- 1000 #设置 B 为一个较大的数值
boot......) #创建一个向量来存储自助法的估计值
for(i in 1:B){
x_new <......ce=T) #创建新数据集
boot_......(x_new) #存储自助法的估计值
}
现在，我们已经构建了自助法的分布，我们可以绘制它并检查其是否服从正态分布。
par(m......1,2)) #将图形放置在一行的两个子区域中
#绘制带有叠加正态密度曲线的自助法分布直方图
hist(boo......)), add=T, col='red', lwd=2)
#创建自助法分布的 QQ 图
qqnorm(......
从这些图中可以明显看出，样本均值 X¯ 的抽样分布稍微右倾。严格来说，在我们完全相信CLT之前，可能需要更多的样本。不过，自助法的分布近似正态分布，因此CLT可能会给出合理的答案。
情况2 - 使用自助法
我们可以首先计算原始数据的样本偏度。
python
#计算样本偏度
n = len(y)
......
我们可以观察到，偏度是正的，表明数据略微向右倾斜。但这个结果有多显著呢？由于样本大小相当大，这是一个很好的自助法（bootstrap）的应用场景。让我们使用以下方法来近似估计 γ̂。
n = len(y) # 获取样本大小
B = 1000 # 设置一个较大的B
boot_sample......
NA, B) # 创建一个向量来存储自助法估计值
for i in 1:B:
y_new = sam......
ace=T) # 创建新的数据集
boot_sam......
) / sd(y_new)^3 # 存储自助法估计值
hist(boot_s......
s=20) # 显示自助法分布
现在，我们已经得到了对抽样分布的近似，我们可以找到一个自助置信区间来表示 γ̂。对于给定的置信水平 C ∈ (0,1)，我们可以找到包含中间 C×100% 的自助法估计值的区间 (a,b)。在R中，可以通过以下方式轻松完成。
# 将置信水平设置为0.95
C = 0.95
alpha = 1 - C
# 获取置信区间
CI = quantile(boot_s......
2))
CI
# 绘制自助法分布并显示置信区间
hist(boot_sampl......
ty=3)
对于这个区间的解释大致如下：我们有95%的置信度，真实的总体偏度在 0.132 和 0.618 之间。因此我们在某种程度上可以相信这个分布的偏度是正的。