r语言Bootstrap自助法重采样构建统计量T抽样分布近似值可视化|代码分享

全文链接：https://tecdat.cn/?p=33939

相关视频

这意味着给定数据值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）。直方图如下所示

这个分布是右偏的还是对称的？很难说。回想一下，分布的总体偏度定义为

这个参数的一个简单估计量（统计量）是下面给出的"样本偏度"

其中y¯和s是数据的样本均值和标准差。那么问题来了，γ^的抽样分布是什么？这个分布肯定是非常难以解析计算的。

自助法

自助法最初由Bradley Efron在1979年发表，被称为一种重新采样技术。Bootstrap（至少这个版本）被称为一种非参数方法，因为它不需要我们对数据做出任何特定分布的假设，这是一个有用的特性。

基本思想是，如果样本数据准确反映了总体，我们可以"重新采样"数据并构建统计量T的抽样分布的近似值。这个近似值有时被称为T的"Bootstrap分布"。需要记住的是，像大多数统计方法一样，当样本量非常小时，Bootstrap可能会失败。

算法其实相当简单，步骤如下：

通过从原始数据中（有放回地）抽样，创建一个“新”数据集，直到你有一个大小为 n 的新数据集。

计算这个新数据集的检验统计量，并将其称为 T1。

重复步骤 1 和 2 多次（比如说 B 次），这样你就得到了一系列的估计值 T1，T2，⋯，TB。这是对 T 的抽样分布的数值近似。

情况1 - 使用自助法

在这个例子中，我们可以使用自助法来近似样本均值 X¯ 的抽样分布。如果自助法的分布看起来近似正态分布，那么我们可以合理地认为中心极限定理（CLT）会给出一个不错的近似结果。如果不是，我们应该怀疑是否能够信任CLT对于这个数据的适用性。

B

boot......) #创建一个向量来存储自助法的估计值

for(i in 1:B){

x_new

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 之间。因此我们在某种程度上可以相信这个分布的偏度是正的。