参数 Bootstrap 方法

假设所研究的总体的分布函数 $F(x;\beta)$ 的形式已知，但其中包含未知参数 $\beta$ （ $\beta$ 可以是向量）。现在已知有一个来自 $F(x;\beta)$ 的样本

利用这一样本求出 $\beta$ 的最大似然估计 $\hat{\beta}$ 。在 $F(x;\beta)$ 中以 $\hat{\beta}$ 代替 $\beta$ 得到 $F(x;\hat{\beta})$ ，接着在 $F(x;\hat{\beta})$ 中产生容量为 $n$ 的样本

X_1^*,X_2^*,\cdots,X_n^* \sim F(x;\hat{\beta})

这种样本可以产生很多个，如产生 $B(B\geq1000)$ 个，就可以利用这些样本对总体进行统计推断。

例1 已知某种电子元件的寿命（单位：h）服从威布尔分布，其分布函数为

F(x) = \left\{ {\begin{array}{{20}}{1 - {e^{ - {{\left( {x/\eta } \right)}^\beta }}},}&{x > 0,}\\{0,}&{\text{其他,}}\end{array}} \beta>0 ,\eta>0 。\right.

概率密度为

f\left( x \right) = \left\{ {\begin{array}{*{20}{c}}
{\frac{\beta }{{{\eta ^\beta }}}{x^{\beta - 1}}{e^{ - {{\left( {x/\eta } \right)}^\beta }}},}&{x > 0,}\\
{0,}&{\text{其他,}}
\end{array}} \right.

已知参数 $\beta=2$ 。今有样本 $142.84$ $97.04$ $32.46$ $69.14$ $85.67$ $114.43$ $41.76$ $163.07$ $108.22$ $63.28$

（1）确定参数 $\eta$ 的最大似然估计。

（2）对于时刻 $t_0=50$ ，求可靠性 $R(50)=1-F(50)=e^{-(50/{\eta})^2}$ 的置信水平为0.95的Bootstrap置信区间。

解（1）设有样本 $x_1,x_2,\cdots,x_n$ ，似然函数为（已将 $\beta=2$ 代入）

L=\prod_{i=1}^n \frac{2}{{\eta}^2} x_i e^{-(x_i/{\eta})^2}
=\frac{2^n}{\eta^{2n}}(\prod_{i=1}^n x_i) e^{-(\sum_{i=1}^n x_i^2)/\eta^2},

\ln{L}=n\ln{2} - 2n\ln{\eta} + \sum_{i=1}^n {\ln{x_i}} - \frac{1}{\eta^2} \sum_{i=1}^n {x_i^n},

令 $\frac{d}{d\eta} \ln{L} = 0$ ，得

\frac{-2n}{\eta} + \frac{2}{\eta^3} \sum_{i=1}^n {x_i^2} = 0

\hat{\eta} = \sqrt{\frac{\sum_{i=1}^n {x_i^2}}{n}} \text{。}

以数据代入，得 $\eta$ 的最大似然估计为 $\hat{\eta}=100.0696$ 。

（2）对于参数 $\beta = 2,\eta = \hat{\eta} = 100.0696$ ，产生服从对应威布尔分布的10000个容量为10的Bootstrap样本。

R代码

eta_est <- function(x) sqrt(mean(x^2))

bootstrap <- function(B, FUN, alpha = 0.05)
{
  rnd <- replicate(B, rweibull(10, 2, 100.0696)) # 产生威布尔随机数
  est <- sort(apply(rnd, 2, FUN))# 计算Bootstrap估计并排序
  
  upper <- est[floor(R * (1 - alpha/2))] # 置信上限
  lower <- est[floor(R * alpha/2)] # 置信下限
  c("CL" = lower, "CU" = upper)
}
bootstrap(1e4, eta_est)
#      CL       CU 
# 69.9376 130.7169

于是， $R(50)$ 的置信水平为0.95的Bootstrap置信区间为 $[0.5998,0.8639]$ 。

非参数 Bootstrap 方法

在很多情况下，总体的分布函数 $F(x;\beta)$ 的形式通常是未知的，这时就无法产生模拟样本，需要另外的方法。

现在设分布 $F$ 未知， $x_1,x_2,\cdots,x_n$ 是来自 $F$ 的样本值， $F_n$ 是相应的经验分布函数。当 $n$ 很大时， $F_n$ 接近 $F$ 。在 $F_n$ 中抽样，就是在原始样本 $x_1,x_2,\cdots,x_n$ 中每次随机地取一个个体作放回抽样。如此得到一个容量为 $n$ 的样本称为Bootstrap样本。相应地、独立地抽得 $B$ 个Bootstrap样本，以这些样本分别求出 $\beta$ 的相应地Bootstrap估计如下： Bootstrap样本1 $x_1^{*1},x_2^{*1},\cdots,x_n^{*1}$ ，Bootstrap估计 $\hat{\beta_1^*}$ ； Bootstrap样本2 $x_1^{*2},x_2^{*2},\cdots,x_n^{*2}$ ，Bootstrap估计 $\hat{\beta_2^*}$ ；

Bootstrap样本B $x_1^{*B},x_2^{*B},\cdots,x_n^{*B}$ ，Bootstrap估计 $\hat{\beta_B^*}$ 。则 $\hat{\beta}$ 的标准误差 $\sqrt{D(\hat{\beta})}$ ，就以

\hat{\sigma_{\hat{\beta}}} = \sqrt{\frac{1}{B-1} \sum_{i=1}^B {(\hat{\beta_i^*} - \bar{\beta^*})^2}}

来估计，其中 $\bar{\beta^*} = \frac{1}{B} \sum_{i=1}^B {\hat{\beta_i^*}}$ 。

例2 设金属元素铂的升华热是具有分布函数 $F$ 的连续型随机变量， $F$ 的中位数 $\theta$ 是未知参数，现测得以下数据（以 kcal/mol 计）： 136.3 136.6 135.8 135.4 134.7 135.0 134.1 143.3 147.8 148.8 134.8 135.2 134.9 149.5 141.2 135.4 134.8 135.8 135.0 133.7 134.4 134.9 134.8 134.5 134.3 135.2 以样本中位数 $M=M(X)$ 作为总体中位数 $\theta$ 的估计，试求均方误差 $MSE=E[(M-\theta)^2]$ 的Bootstrap估计。

R代码

bootstrap <- function(B, x)
{
  m <- median(x) # 样本中位数
  bs <- replicate(B, sample(x, length(x), replace = T)) # 抽取Bootstrap样本
  bs <- apply(bs, 2, median) # 计算Bootstrap样本中位数
  sum((bs - m)^2)/(B - 1) # 计算均方误差
  
}
x <- c(136.3, 136.6, 135.8, 135.4, 134.7,
       135.0, 134.1, 143.3, 147.8, 148.8,
       134.8, 135.2, 134.9, 149.5, 141.2,
       135.4, 134.8, 135.8, 135.0, 133.7,
       134.4, 134.9, 134.8, 134.5, 134.3, 135.2)
bootstrap(1e4, x)
# [1] 0.07338434

Bootstrap方法估计参数

参数 Bootstrap 方法

这种样本可以产生很多个，如产生个，就可以利用这些样本对总体进行统计推断。

非参数 Bootstrap 方法

这种样本可以产生很多个，如产生 $B(B\geq1000)$ 个，就可以利用这些样本对总体进行统计推断。