Bootstrap的三种置信区间：正态、枢轴量、分位数

最新推荐文章于 2024-06-03 09:55:46 发布

klinfelter

最新推荐文章于 2024-06-03 09:55:46 发布

阅读量3.2k

点赞数 4

文章标签： bootstrap python

本文链接：https://blog.csdn.net/weixin_43851599/article/details/128221615

版权

Bootstrap置信区间

Bootstrap的原理在很多教材和博客上都有详细的解释，这里就不赘述了。在构造Bootstrap置信区间时我们有三种方法：正态近似法、枢轴量法、分位数法。

给定样本 $X_1,\cdots ,X_n$ ，假设统计量 $T_n = g(X1,\cdots ,X_n)$ 是参数 $\theta = T(F)$ 的一个估计，Bootstrap抽样得到 $T_{n,1}^{\star}, \cdots ,T_{n,B}^{\star}$ ，用 $\hat{\theta}_{n,1}^{\star}, \cdots \hat{\theta}_{n,B}^{\star}$ 表示 $\hat{\theta}_n$ 的Bootstrap复本， $\theta_{\alpha}^{\star}$ 表示 $\hat{\theta}_{n,1}^{\star}, \cdots \hat{\theta}_{n,B}^{\star}$ 的下 $\alpha$ 分位数。求 $\theta$ 的置信水平为 $1-\alpha$ 置信区间：

· 正态置信区间

正态近似法是这三种方法中最常用、最简单的。

$(\hat{\theta} - z_{\frac{\alpha}{2}}\sqrt{v_{boot}}, \hat{\theta} + z_{\frac{\alpha}{2}}\sqrt{v_{boot}})$

· 枢轴量置信区间

$(2\hat{\theta} - \theta_{1-\tfrac{\alpha}{2}}^{\star}, 2\hat{\theta} - \theta_{\tfrac{\alpha}{2}}^{\star})$

· 分位数置信区间

$(\hat{\theta}_{\tfrac{\alpha}{2}}^{\star}, \hat{\theta}_{1-\tfrac{\alpha}{2}}^{\star})$

代码实现

以下面这个题目为例实现Bootstrap和三种方法的置信区间计算：

（计算机实验）用随机模拟比较不同的Bootstrap置信区间方法. 令n=50，并令 $T(F) = \int \tfrac{(x-\mu)^3}{\sigma^3}dF(x)$ 为偏度. 抽取随机样本使得 $Y_1, \cdots , Y_n\sim N(0,1)$ ，令 $X_i = e^{Y_i}, i=1,\cdots ,n$ . 根据数据 $X_1, \cdots ,X_n$ ，为 $T(F)$ 构造三种类型的95%的Bootstrap置信区间. 重复整个过程若干次，估计这三个区间的真实值.

三种方法的Bootstrap过程是一样的，仅在获得Bootstrap复本后的计算式不同。下面的代码可以直接运行，可以得到三个置信区间的输出。

T_F()这个函数实现了题目中偏度的计算。对于任何其他的估计函数，修改T_F()这个函数即可。

import numpy as np
from scipy.stats import norm

y = np.random.normal(loc = 0.0, scale = 1.0, size = 50)
x = np.exp(y)

def T_F(data):
    mu = np.mean(data)
    sigma = np.std(data)
    temp = ((data - mu) / sigma) ** 3
    return (sum(temp) / len(temp))

def bootstrap_replicate_1d(data, func):
    bs_sample=np.random.choice(data, len(data))
    return func(bs_sample)

def bootstrap_res(data, func, size=1):
    bs_replicates = np.empty(size)
    for i in range(size):
        bs_replicates[i] = bootstrap_replicate_1d(data, func)
    return bs_replicates

T_hat = T_F(x)
B_hat = bootstrap_res(x, T_F, size = 10000)
B_sort = sorted(B_hat)    # 这一步其实是不必要的，但如果你想用其他方法求分位数可能会用到这一排序结果
std = np.std(B_hat)
alpha = 0.05
z = norm.ppf(q = 1-alpha/2, loc = 0, scale = 1)

print("Bootstrap normal interval: (", T_hat - z * std, ", ", T_hat + z * std, ")")
print("Bootstrap pivotal interval: (", 2 * T_hat - np.percentile(B_sort, (1-alpha/2)*100), ", ", 2 * T_hat - np.percentile(B_sort, alpha/2*100), ")")
print("Bootstrap percentile interval: (", np.percentile(B_sort, alpha/2*100), ", ", np.percentile(B_sort, (1-alpha/2)*100), ")")