非参数Bootstrap方法

• 设总体的分布$F$未知，但按放回抽样的方法抽取了一个容量为$n$的样本，称为Bootstrap样本或称为自助样本。独立地取多个Bootstrap样本，利用这些样本信息对总体$F$进行推断，这种方法称为非参数Bootstrap方法，又称为自助法。这一方法可用于对总体知之甚少地情况。
• 优点：不需要对总体分布有任何假设，而且可以使用于小样本，且能用于各种统计量。

估计量的标准误差的Bootstrap估计

• 在估计总体位置参数$\theta$时，不仅要给出$\theta$的估计$\hat{\theta }$，还需要这一估计的标准误差，常用${\sigma }_{\hat{\theta }}=\sqrt{D(\hat{\theta })}$度量。

步骤

1. 相继、独立地从已知的容量为$n$的样本中抽出$B(B\ge 1000)$个容量为$n$的Bootstrap样本

$x_1^{\ast i},x_2^{\ast i},...,x_n^{\ast i},i=1,2,...,B$

2. 计算

$${\hat{\sigma }}_{\hat{\theta }}=\sqrt{\frac{1}{B-1}\sum _{i=1}^B({\hat{\theta }}_i^{\ast }-{\bar{\theta }}^{\ast }{)}^2}$$
式中： ${\hat{\theta }}_i^{\ast }=\hat{\theta }(x_1^{\ast i},x_2^{\ast i},...,x_n^{\ast i}),i=1,2,...,B{\bar{\theta }}^{\ast }=\frac{1}{B}{\sum }_{i=1}^B{\hat{\theta }}_i^{\ast }$

估计量的均方误差的Bootstrap估计

例1

• 有30窝仔猪出生时各窝猪的存活只数为

9 8 10 12 11 12 7 9 11 8 9 7 7 8 9 7 9 9 10 9 9 9 12 10 10 9 13 11 13 9
（1）试求中位数估计$M$的标准误差的Bootstrap估计。
（2）求均方误差$MSE=E[(M-\theta {)}^2]$的估计。

data=[9  8  10  12  11  12  7  9  11  8  9  7  7  8  9  7  9  9  10  9  9  9  12  10  10  9  13  11  13  9];
b=bootstrp(1000,@(x)quantile(x,0.5),data);%1000组样本，@(x)quantile(x,0.5)求中位数的函数，data原始数据
b_std=std(b);%计算b的标准差
b_var=mean((b-quantile(data,0.5)).^2);%求均方误差


b1=bootstrp(1000,@mean,data);%平均数
b1_mean=std(b1);

Bootstrap置信区间

• 设$X=(X_1,X_2,....,X_n)$是来自总体$F$容量为$n$的样本，$x=(x_1,x_2,...,x_n)$是一个已知的样本值。$F$中含有位置参数$\theta$，$\hat{\theta }=\hat{\theta }(X_1,X_2,...,X_n)$是$\theta$的估计量。现在求$\theta$的置信水平为$1-\alpha$的置信区间。

步骤及原理

1. 独立从样本$x$中抽出$B$个容量为$n$的Bootstrap样本，对于每个Bootstrap样本求出$\theta$的Bootstrap估计：${\hat{\theta }}_1^{\ast },{\hat{\theta }}_2^{\ast },...,{\hat{\theta }}_B^{\ast }$。
2. 将估计从小到大排序：${\hat{\theta }}_{(1)}^{\ast },{\hat{\theta }}_{(2)}^{\ast },...,{\hat{\theta }}_{(B)}^{\ast }$。
3. 求出近似分位数${\hat{\theta }}_{\alpha /2}^{\ast },{\hat{\theta }}_{1-\alpha /2}^{\ast }$使得
   P { θ ^ α / 2 ∗ < θ ^ ∗ < θ ^ 1 − α / 2 ∗ } = 1 − α P\{\hat\theta_{\alpha/2}^*<\hat\theta^*<\hat\theta_{1-\alpha/2}^*\}=1-\alpha P{θ^α/2∗​<θ^∗<θ^1−α/2∗​}=1−α
   于是近似的有
   P { θ ^ α / 2 ∗ < θ ^ < θ ^ 1 − α / 2 ∗ } = 1 − α P\{\hat\theta_{\alpha/2}^*<\hat\theta<\hat\theta_{1-\alpha/2}^*\}=1-\alpha P{θ^α/2∗​<θ^<θ^1−α/2∗​}=1−α
4. 记$k_1=[B\times \frac{\alpha }{2}],k_2=[B\times (1-\frac{\alpha }{2})]$，以${\hat{\theta }}_{(k_1)}^{\ast }$和${\hat{\theta }}_{(k_2)}^{\ast }$分别作为分位数${\hat{\theta }}_{\alpha /2}^{\ast },{\hat{\theta }}_{1-\alpha /2}^{\ast }$的估计，得到近似等式
   P { θ ^ ( k 1 ) ∗ < θ ^ < θ ^ ( k 2 ) ∗ } = 1 − α P\{\hat\theta_{(k_1)}^*<\hat\theta<\hat\theta_{(k_2)}^*\}=1-\alpha P{θ^(k1​)∗​<θ^<θ^(k2​)∗​}=1−α
   得到$\theta$的置信水平为$1-\alpha$的Bootstrap置信区间为$({\hat{\theta }}_{(k_1)}^{\ast },{\hat{\theta }}_{(k_1)}^{\ast })$，这种方法称为分位数法。

续例1

• 以样本均值$\bar{x}$作为总体均值$\mu$的估计，以标准差$s$作为总体标准差$\sigma$的估计，按分位数法求$\mu ，\sigma$的置信水平为0.90的Bootstrap置信区间。

b=bootci(1000,{@(x)[mean(x),std(x)],data},'alpha',0.1)%bootci得到Bootstrap置信区间

结果：第一列是均值的置信区间，第二列是标准差的置信区间
b =
    9.0667    1.4368
   10.0667    2.0609

参数Bootstrap方法

• 假设所研究的总体的分布函数$F=(x;\beta )$的形式已知，但其中高喊未知参数$\beta$。现在已知有一个来自总体的样本$X_1,X_2,...,X_n$。

步骤

1. 由该样本得到$\beta$的极大似然估计$\hat{\beta }$。
2. 在总体分布为$F(x;\hat{\beta })$中产生$B(B\ge 1000)$个容量为$n$的样本
   $$X_1^{\ast },X_2^{\ast },...,X_n^{\ast }\sim F(x;\hat{\beta })$$
3. 用非参数Bootstrap置信区间的方法得到$\beta$的Bootstrap置信区间。

例2：

• 已知某电子原件的寿命服从威布尔分布，其分布函数如下
  $$F(x)=\{\begin{array}{ll}1-e^{-(x/\eta {)}^{\beta }}& x>0\\ 0& 其他\end{array}\beta >0,\eta >0$$
  已知参数$\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）求似然估计是概率论中的基础，不在此详细阐述，结果为$\hat{\eta }=\sqrt{\frac{{\sum }_{i=1}^n{x}_i^2}{n}}=100.0696$。
  （2）对于参数$\beta =2,\eta =\hat{\eta }=100.0696$，产生服从对应韦布尔分布的5000个容量为10的Bootstrap样本。
  对于每个样本计算$\eta$的Bootstrap估计${\hat{\eta }}_i^{\ast }=\sqrt{\frac{{\sum }_{j=1}^{10}(x_j^{\ast i}{)}^2}{10}}$。
  将5000个${\eta }_i^{\ast }$自小到大排列，取坐起第250（[5000×0.05]=250）位，得${\eta }_{(250)}^{\ast }=73.3758$。
  所以置信水平位0.95得Bootstrap单侧置信下限为$e^{-(50/{\hat{\eta }}_{(250)}^{\ast }{)}^2}=0.6286$。

clc,clear
a=[142.84  97.04  32.46  69.14  85.67  114.43  41.76  163.07  108.22  63.28];
eta=sqrt(mean(a.^2));
beta=2;B=5000;alpha=0.05;
b=wblrnd(eta,beta,[B,10]);%产生shape=(B,10)的韦布尔分布的随机数
etahat=sqrt(mean(b.^2,2));%计算每个样本对应的最大似然估计
seteta=sort(etahat);%把etahat从小到大排序
k=floor(B*alpha);%取整数部分
se=seteta(k);%提取相应位置的估计量
Rt0=exp(-(50/se)^2);