数理统计中常用的方差估计方法

背景

随着科技的发展, 我们正进入大数据时代,对数据进行统计推断的需求也越来越多。统计推断中的一个重要问题便是参数估计,衡量一个估计量好坏的标准之一便是这个估计量的方差。通常，方差是未知的，只能利用本身现有样本数据进行估计。这篇文章主要介绍一下最近学习的一些方差估计方法，主要参考Wolter的《Introduction to variance estimation》一书。

随机组方法(Random Groups)

随机组是一种实验设计方法，其中样本通过随机抽样被分配到不同的组别中，在不同组里，我们可以利用多个独立样本构建参数的估计量，并通过这些估计量之间的方差来估计方差。这种方法可以简化方差估计的过程，该方法已经比较成熟,也是我们目前最常用的方差估计方法之一.
现在考虑有来自某总体的$k$组样本，需要估计的参数为$\theta$, 记${\hat{\theta }}_{\alpha }$是第$\alpha$组对$\theta$的估计量$(\alpha =1,2,\cdots ,k)$,下面的定理给出了随机组方差估计的形式：

定理1.1.设${\hat{\theta }}_1,\dots ,{\hat{\theta }}_k$ 不相关且期望均为 E { θ ^ 1 } = μ \mathrm{E}\{\hat{\theta}_1\}=\mu E{θ^1​}=μ.记 $\hat{\bar{\theta }}$ 定义为
$$\hat{\bar{\theta }}=\sum _{\alpha =1}^k{\hat{\theta }}_{\alpha }/k.$$
则 E { θ ˉ ^ } = μ \mathrm{E}\{\hat{\bar{\theta}}\}=\mu E{θˉ^}=μ ,且其方差 Var ⁡ { θ ˉ ^ } \operatorname{Var}\{\hat{\bar{\theta}}\} Var{θˉ^} 的无偏估计为
$$v(\hat{\bar{\theta }})=\sum _{\alpha =1}^k{\left({\hat{\theta }}_{\alpha }-\hat{\bar{\theta }}\right)}^2/k(k-1).$$
注1：$\hat{\bar{\theta }}$可以作为全样本下参数$\mu$的估计,$v(\hat{\bar{\theta }})$是$\hat{\bar{\theta }}$的随机组方差估计量(random groups).注2：一般都取$\theta =\mu$, 可以发现，当每一个组别的样本量都是$1$时, $\hat{\bar{\theta }}$和$v(\hat{\bar{\theta }})$分别是我们平时最熟悉的平均值和标准差.对估计量$\hat{\bar{\theta }}$而言,有了方差估计便可以尝试算置信区间或者做假设检验了,此外，一般会有正态分布和t分布假设. 设${\hat{\theta }}_1,\dots ,{\hat{\theta }}_k$独立且分布均为$N\left(\theta ,{\sigma }^2\right)$,则有:
(i)统计量$z=(\hat{\bar{\theta }}-\theta )/\sqrt{{\sigma }^2/k}\sim N(0,1).$
(ii)统计量$t=(\hat{\bar{\theta }}-\theta )/\sqrt{v(\hat{\bar{\theta }})}\sim t(k-1).$
注：当k足够大时候$\theta$的 $(1-\alpha )100\%$ 置信区间为$\left(\hat{\bar{\theta }}-z_{\alpha /2}\sqrt{v(\hat{\bar{\theta }})},\hat{\theta }+z_{\alpha /2}\sqrt{v(\hat{\bar{\theta }}}\right)),$
其中 $z_{\alpha /2}$ 标准正态分布$N(0,1)$ 的上侧 $\alpha /2$ 分位点. 可以发现，这个定理其实是我们在学数理统计时的抽样分布.
前面我们介绍了随机组方差估计, 定理也表明其是无偏的, 但应该如何衡量方差估计的稳定性呢?还有分组数$k$应该怎么选择呢?这里作者给出了一个常用的准则,便是变异系数:
CV ⁡ { v ( θ ˉ ^ ) } = [ Var ⁡ { v ( θ ˉ ^ ) } ] 1 / 2 / Var ⁡ { θ ˉ ^ } . \operatorname{CV}\{v(\hat{\bar{\theta}})\}=[\operatorname{Var}\{v(\hat{\bar{\theta}})\}]^{1 / 2} / \operatorname{Var}\{\hat{\bar{\theta}}\} . CV{v(θˉ^)}=[Var{v(θˉ^)}]1/2/Var{θˉ^}.
其中分子是该估计量的方差开根号,分母是其期望.
定理1.2. 设${\hat{\theta }}_1,\dots ,{\hat{\theta }}_k$ 独立同分布, $v(\hat{\bar{\theta }})$ 为之前定义.则
CV ⁡ { v ( θ ˉ ^ ) } = { β 4 ( θ ^ 1 ) − ( k − 3 ) / ( k − 1 ) k } 1 / 2 , \operatorname{CV}\{v(\hat{\bar{\theta}})\}=\left\{\frac{\beta_4\left(\hat{\theta}_1\right)-(k-3) /(k-1)}{k}\right\}^{1 / 2}, CV{v(θˉ^)}=⎩ ⎨ ⎧​kβ4​(θ^1​)−(k−3)/(k−1)​⎭ ⎬ ⎫​1/2,
其中
$$\begin{array}{rl}{\beta }_4\left({\hat{\theta }}_1\right)& =\frac{\mathrm{E}\left\{{\left({\hat{\theta }}_1-\mu \right)}^4\right\}}{{\left[\mathrm{E}\left\{{\left({\hat{\theta }}_1-\mu \right)}^2\right\}\right]}^2},\\ \mu & =\mathrm{E}\left\{{\hat{\theta }}_1\right\}.\end{array}$$
可以看出, CV 依赖于${\beta }_4\left({\hat{\theta }}_1\right)$和$k$, 当$k$较小或者峰度 ${\beta }_4\left({\hat{\theta }}_1\right)$较大时,CV 偏大. 因此,在选择组别数$k$时,应当越多越好.

Jackknife

Jackknife方法最早由Quenouille（1949年）提出, 并用于估计量的纠偏(debias).
Tukey（1958年）提出,可以将各个子样本估计量合理地视为独立且同分布的随机变量,从而得到一个简单却实用的方差估计量.
首先介绍Jackknife方法的基本定义,设 $Y_1,\dots ,Y_n$是来自总体$F(y)$的独立同分布样本. 对于感兴趣的参数$\theta$, 假设$\hat{\theta }$ 是基于全部样本的估计量. 现在考虑将完整样本分成$k$组，每组有$m$个观测值，其中$n=mk$, 对于每个$\alpha =1,\dots ,k$，计算去掉第$\alpha$组后剩余 $m(k-1)$个观测样本的估计量${\hat{\theta }}_{(-\alpha )}$, 然后定义
$${\hat{\theta }}_{\alpha }=k\hat{\theta }-(k-1){\hat{\theta }}_{(-\alpha )}.$$
Quenouille提出的估计量便是便是这些${\hat{\theta }}_{\alpha }$的均值,
$${\hat{\theta }}_{jack}=\sum _{\alpha =1}^k{\hat{\theta }}_{\alpha }/k,$$
注：当$k=n,m=1$时, ${\hat{\theta }}_{(-\alpha )}$便是去掉单个样本的情况.
相信大家看到这个${\hat{\theta }}_{\alpha }$形式的时候有点懵,这里我尝试用一个特殊的例子解释一下${\hat{\theta }}_{(-\alpha )}$是如何达到纠偏效果的. 假设$k=n,m=1$, $\hat{\theta }$估计$\theta$具有一阶偏差,即
E { θ ^ } = θ + a 1 / n , E\{\hat{\theta}\}=\theta+a_1 / n , E{θ^}=θ+a1​/n,
则有 E { θ ^ ( − α ) } = θ + a 1 / ( n − 1 ) , E\{\hat{\theta}_{(-\alpha)}\}=\theta+a_1 / (n - 1), E{θ^(−α)​}=θ+a1​/(n−1), 为了消去偏差带有$a_1$项的偏差, 自然可以通过作差,简单计算可以得到 E { n θ ^ − ( n − 1 ) θ ^ ( − α ) } = θ E\{n\hat{\theta} - (n - 1)\hat{\theta}_{(-\alpha)}\} = \theta E{nθ^−(n−1)θ^(−α)​}=θ, 即$n\hat{\theta }-(n-1){\hat{\theta }}_{(-\alpha )}$是$\theta$的无偏估计, 此时纠偏的目的就达到了,这也可以简单解释Jackknife估计量是如何得到的.
按照Tukey的建议，将${\hat{\theta }}_{\alpha }$视为近似独立且同分布的随机变量, 此时Jackknife方差估计为:
$$v_1({\hat{\theta }}_{jack})=\frac{{\sum }_{\alpha =1}^k{\left({\hat{\theta }}_{\alpha }-{\hat{\theta }}_{jack}\right)}^2}{k(k-1)}.$$
这里${\hat{\theta }}_{jack}$也可以直接用$\hat{\theta }$代替,
即$$v_2({\hat{\theta }}_{jack})=\frac{{\sum }_{\alpha =1}^k{\left({\hat{\theta }}_{\alpha }-\hat{\theta }\right)}^2}{k(k-1)}.$$
这两个方差估计的关系是
$$v_2({\hat{\theta }}_{jack})=v_1({\hat{\theta }}_{jack})+(\hat{\theta }-{\hat{\theta }}_{jack}{)}^2/(k-1)$$
接下来简单介绍一些Jackknife方法的性质.
许多重要参数可以表示为 $\theta =g(\mu )$的形式，其中$\mu$表示期望 E { Y i } = μ \mathrm{E}\left\{Y_i\right\}=\mu E{Yi​}=μ。尽管
$$\bar{Y}=n^{-1}\sum _{j=1}^n{Y}_j$$
是$\mu$的无偏估计量，但$\hat{\theta }=g(\bar{Y})$通常是$\theta =g(\mu )$的有偏估计量. Quenouille给出的估计量形式为
$${\hat{\theta }}_{jack}=kg(\bar{Y})-(k-1)k^{-1}\sum _{\alpha =1}^{k}g\left({\bar{Y}}_{(-\alpha )}\right),$$
其中 ${\bar{Y}}_{(-\alpha )}$表示去掉第$\alpha$组观测值后$m(k-1)$个观测值的样本均值.
下面的定理保证了其渐近性结果.
定理2.1. 令 { Y j } \left\{Y_j\right\} {Yj​}是独立同分布的随机变量序列，其均值为$\mu$，方差为$0<{\sigma }^2<\infty$.令$g(\cdot )$是定义在实数域上的函数，在 $\mu$的邻域内具有有界二阶导数。那么，$k\to \infty ,k^{1/2}({\hat{\theta }}_{jack}-\theta )$在分布上收敛于均值为零，方差为 σ 2 { g ′ ( μ ) } 2 \sigma^2\left\{g^{\prime}(\mu)\right\}^2 σ2{g′(μ)}2的正态分布随机变量，其中$g^{\prime }(\mu )$是 $g(\cdot )$在$\mu$处的一阶导数。
该定理保证了参数一般形式下的渐近性结果,极限分布有点像我们后面要介绍的Delta方法.
最后又是关于分组数$k$的选择了,主要有两个考量因素：计算成本和估计量的准确性,$k$越大,计算成本越大,精度更高;$k$越小,计算成本越小,但准确性可能较差.在大数据集上,在选择时希望在计算成本和精确度之间找到平衡.
最后这里仅仅介绍了Jackknife在方差估计方面的结果,有兴趣的同学可以了解他的纠偏理论.

Bootstrap

接下来介绍大名鼎鼎的Bootstrap方法.在前面小节中，我们讨论了几种基于重复的方差估计方法, Bootstrap核心思想是有放回抽样,也是所谓的重复方法,那么Bootstrap与其他复制方法有何不同？在最简单的情况下，随机组是基于样本量为$n/k$的重复；Jackknife方法则使用样本量为$n-1$的重复。相比之下，Bootstrap使用潜在样本量为$n^{\star }$的任何重复.
现考虑常规的Bootstrap方法,设$Y_1,Y_2,\cdots ,Y_n$是来自分布函数为$F$总体的独立同分布样本,设$\theta$是感兴趣的未知参数, $\hat{\theta }$是基于样本对$\theta$的估计,我们希望能估计$\hat{\theta }$的方差,即$\mathrm{Var}(\hat{\theta })$. Bootstrap方法步骤如下：
(i) 考虑一个较大的整数 $A$, 对于$\alpha =$ $1,\dots ,A$,每次从 { Y 1 , Y 2 , ⋯   , Y n } \{Y_1,Y_2,\cdots,Y_n\} {Y1​,Y2​,⋯,Yn​}进行$n^{\star }$次有放回抽样,得到$Y_{\alpha 1}^{\ast },\dots ,Y_{\alpha n}^{\ast }$;
(ii)对每次有放回抽样得到的样本,计算相应的估计量${\hat{\theta }}_{\alpha }^{\ast }$;
(iii)计算 ${\hat{\theta }}_{\alpha }^{\ast }$之间的方差:
$$\begin{array}{c}v(\hat{\theta })=\frac{1}{A-1}\sum _{\alpha =1}^A{\left({\hat{\theta }}_{\alpha }^{\ast }-{\hat{\theta }}^{\ast }\right)}^2,\\ {\hat{\theta }}^{\ast }=\frac{1}{A}\sum _{\alpha =1}^A{\hat{\theta }}_{\alpha }^{\ast }.\end{array}$$
自然地, $v(\hat{\theta })$便可以作为估计量$\hat{\theta }$的方差估计,一般取$n^{\star }=n$, 大佬Efron 和 Tibshirani (1986) 指出 $A$取50到200能满足大多数情况,但是我看好像一般取1000比较多? 毕竟$A$多多益善.
Bootstrap是非常实用的Model-Free非参方法,除了估方差之外,还可以估偏差,估分布,算置信区间,做假设检验等等,这里只是简单介绍一下估方差的基本用法,后续可以写篇文章专门介绍Bootstrap.

Taylor展开线性近似

在很多情况下,估计量可能非线性的,一些常见例子比如比值、相关系数、回归系数等都是这种情况. 非线性估计量的方差通常没有显性表达式, 一种估计非线性估计量方差的方法是利用观测值的线性函数来近似, 基于线性估计量, 就可以应用特定的方差估计方法,这虽然可能会导致一定偏差, 但通常是相合的.
如何进行线性近似?可以考虑一阶Taylor展开,在此之前先介绍一些概率意义下阶的概念.
定义. 对于一列$p$维随机向量${\mathbf{Y}}_n$,实数列$r_n$,若 $\forall \epsilon >0$,存在实正数$M_{\epsilon }$,满足 P { ∣ Y j n ∣ ≥ M ε r n } ≤ ε , j = 1 , … , p , ∀ n ∈ N + , P\left\{\left|Y_{j n}\right| \geq M_{\varepsilon} r_n\right\} \leq \varepsilon, \quad j=1, \ldots, p, \quad \forall n \in \mathbb{N}^{+} , P{∣Yjn​∣≥Mε​rn​}≤ε,j=1,…,p,∀n∈N+,
则称
$${\mathbf{Y}}_n=O_p\left(r_n\right).$$
这个定义主要表明一个随机向量的阶由其每个维度决定.
考虑$g\left({\mathbf{Y}}_n\right)$在$g(\mathbf{a})$点的泰勒展开:
$$g\left({\mathbf{Y}}_n\right)=g(\mathbf{a})+{\mathbf{a}}^{\top }\nabla g(\mathbf{a})+R_n\left({\mathbf{Y}}_n,\mathbf{a}\right),$$
其中余项
$$R_n\left({\mathbf{Y}}_n,\mathbf{a}\right)=\sum _{j=1}^p\sum _{i=1}^p\frac{1}{2!}\frac{{\partial }^{2}g(\ddot{\mathbf{a}})}{\partial y_j\partial y_i}\left(Y_{jn}-a_j\right)\left(Y_{in}-a_i\right),$$
$\partial g(\mathbf{a})/\partial y_j$ 是 $g(\mathbf{y})$ 关于 $\mathbf{y}$ 的第 $j$ 个分量在 $\mathbf{y}=\ddot{\mathbf{a}}$处的偏导数，$\mathbf{y}=\ddot{\mathbf{a}}$是 $g(\mathbf{y})$ 关于 $y_j$ 和 $y_i$ 在 $\ddot{\mathbf{a}}$ 处的二阶偏导数，而 $\ddot{\mathbf{a}}$是连接${\mathbf{Y}}_n$ 和$\mathbf{a}$的线段上的某一点.
定理4.1 令${\mathbf{Y}}_n=\mathbf{a}+O_p\left(r_n\right)$, 且$r_n\to 0,n\to \infty$,则$g\left({\mathbf{Y}}_n\right)$ 可以由上式的泰勒展开形式表示,并且余项$R_n\left({\mathbf{Y}}_n,\mathbf{a}\right)=O_p\left(r_n^2\right)$.
注：该定理给出了${\mathbf{Y}}_n$的余项和泰勒展开余项$R_n\left({\mathbf{Y}}_n,\mathbf{a}\right)$的关系,
在该定理基础上,很容易得到下面的定理.
定理4.2.令${\mathbf{Y}}_n=\mathbf{a}+O_p\left(r_n\right),$ 其中 $r_n\to O,n\to \infty$,令
E { Y n } = a , E { ( Y n − a ) ⊤ ( Y n − a ) } = Σ n < ∞ . \begin{aligned} \mathrm{E}\left\{\mathbf{Y}_n\right\} & =\mathbf{a}, \\ \mathrm{E}\left\{\left(\mathbf{Y}_n-\mathbf{a}\right)^{\top}\left(\mathbf{Y}_n-\mathbf{a}\right)\right\} & ={\Sigma}_n<\infty . \end{aligned} E{Yn​}E{(Yn​−a)⊤(Yn​−a)}​=a,=Σn​<∞.​
则
$$\mathrm{E}\left\{{\left(g\left({\mathbf{Y}}_n\right)-g(\mathbf{a})\right)}^2\right\}=\nabla g(\mathbf{a}{)}^{\top }{\Sigma }_n\nabla g(\mathbf{a})+O_p\left(r_n^3\right),$$
当${\mathbf{Y}}_n$是一列独立同分布随机变量的均值时,可以得到更强的结论
定理4.2.若${\mathbf{Y}}_n$是一列独立同分布随机变量的均值,其中每个随机变量都是$p$维的,期望为$\mathbf{a}$ ,协方差矩阵为 $\Sigma$，且四阶矩存在,$\mathbf{g}(\mathbf{y})$在$\mathbf{a}$的的邻域内具有连续三阶导数的函数,则
$$\mathrm{E}\left\{{\left(g\left({\mathbf{Y}}_n\right)-g(\mathbf{a})\right)}^2\right\}=\frac{1}{n}\mathbf{d}{\Sigma }_n{\mathbf{d}}^{\prime }+O_p\left(n^{-\frac{1}{2}}\right),$$
注：这个定理很好地解决了$\mathbf{g}\left({\mathbf{Y}}_n\right)$的方差估计问题,即只需要知道 ${\mathbf{Y}}_n$的一阶矩和二阶矩信息,然后求$\mathbf{g}$的梯度就可以了. 看到这里,可以发现$Delta$方法和这个定理结果非常类似,不同点主要在于$Delta$方法在渐近正态假设下,不仅可以得到方差估计,还可以得到渐近正态的结果.
Delta方法. 设 $\sqrt{\mathrm{n}}\left({\mathbf{T}}_{\mathrm{n}}-\mathbf{\theta }\right)\stackrel{\mathrm{d}}{\to }{\mathrm{N}}_{\mathrm{p}}(0,\Sigma (\mathbf{\theta }))$. 令 $\mathbf{g}:{\mathbb{R}}^{\mathrm{p}}\mapsto {\mathbb{R}}^{\mathrm{m}}$ 在 $\mathbf{\theta }$ 可微且有非零梯度 $\nabla \mathrm{g}(\mathbf{\theta })$ ,则 n { g ( T n ) − g ( θ ) } → d N m ( 0 , ∇ g ( θ ) ⊤ Σ ( θ ) ∇ g ( θ ) ) \sqrt{\mathrm{n}}\left\{\boldsymbol{g}\left(\boldsymbol{T}_{\mathrm{n}}\right)-\boldsymbol{g}(\boldsymbol{\theta})\right\} \xrightarrow{\mathrm{d}} \mathrm{N}_{\mathrm{m}}\left(\mathbf{0}, \nabla \boldsymbol{g}(\boldsymbol{\theta})^{\top} \Sigma(\boldsymbol{\theta}) \nabla \boldsymbol{g}(\boldsymbol{\theta})\right) n ​{g(Tn​)−g(θ)}d ​Nm​(0,∇g(θ)⊤Σ(θ)∇g(θ)).
注：Delta方法的证明也是通过泰勒展开,感兴趣的同学可以尝试证明一下.
下面我们考虑一些具体例子,试试如何用泰勒展开来估方差.
例1.
设$X_1,X_2,\cdots ,X_n$是来自同一总体独立同分布的样本,$Y_1,Y_2,\cdots ,Y_n$亦是来自同一总体独立同分布的样本.对于估计量$\hat{\theta }=\bar{X}/\bar{Y}$,令$g(x,y)=x/y$,则$\nabla g(x,y{)}^{\top }=(1/y,-x/y^2)$.令$\mathbf{a}=(\bar{X},\bar{Y}{)}^{\top }$,对于$\Sigma =\mathrm{Cov}(\mathbf{a},\mathbf{a})$,可知${\sigma }_{11}=\mathrm{Var}(\bar{X})$,${\sigma }_{12}={\sigma }_{21}=\mathrm{Cov}(\bar{X},\bar{Y})$,${\sigma }_{22}=\mathrm{Var}(\bar{Y})$,这里的${\sigma }_{ij}$都可以用样本$X_1,X_2,\cdots ,X_n$和$Y_1,Y_2,\cdots ,Y_n$估出来,记为$\widehat{{\sigma }_{ij}}$,从而可以得到 $\hat{\Sigma }$.此时
$$\begin{array}{rl}v(\hat{\theta })& =\nabla g(\mathbf{a}{)}^{\top }\hat{\Sigma }\nabla g(\mathbf{a})\\ & =(1/\bar{Y},-\bar{X}/{\bar{Y}}^2)\hat{\Sigma }(1/\bar{Y},-\bar{X}/{\bar{Y}}^2{)}^{\top }\\ & =\frac{{\hat{\sigma }}_{11}}{{\bar{Y}}^2}+\frac{{\hat{\sigma }}_{22}{\bar{X}}^2}{{\bar{Y}}^4}-\frac{2{\hat{\sigma }}_{12}\bar{X}}{{\bar{Y}}^3}.\end{array}$$
其他类似估计量,如 X ˉ Y ˉ , X ˉ 1 / Y ˉ 1 − X ˉ 2 / Y ˉ 2 , exp ⁡ { X ˉ } \bar{X}\bar{Y},\bar{X}_1 / \bar{Y}_1 - \bar{X}_2 / \bar{Y}_2,\exp\{\bar{X}\} XˉYˉ,Xˉ1​/Yˉ1​−Xˉ2​/Yˉ2​,exp{Xˉ},甚至${\left(X_1{X}_2\cdots X_N\right)}^{1/N}$等,都可以利用同样方法估方差.

总结

本文主要介绍了几种常用的方差估计方法,包括Random Groups,Jacknife,Bootstrap,Delta方法等,Random Groups其实和我们常用的方差公式差不多.Jacknife和Bootstrap在样本量较大时计算量可能有点大,相对而言,Delta方法普适性较强,可以很好地处理非线性估计量的方差问题,后续将介绍一些降方差的方法和Bootstrap原理.

参考文献：Wolter K M. Introduction to variance estimation[M]. New York: Springer, 2007.