第 2 章 估计

2.1 参数估计

在实际问题中,对于一个总体\(X\)往往是仅知其分布的类型\(F_\theta\),而参数\(\theta=(\theta_1,\dots,\theta_m)\in \Theta \subset \mathbb{R}^m\)是未知的。对任给的实值函数

\[g:\ \mathbb{R}^m\to \mathbb{R},\]

如何根据\(X\)的样本\(x_1,\dots,x_n\)估计\(g( \theta)\)的值呢?这就是统计推断中的“参数估计”问题。

点估计:寻找一个统计量\(\hat{\theta} = T(X_1,\dots,X_n)\)作为\(\theta\)的点估计

区间估计:寻找两个统计量\(\hat{ \theta}_1 = T_1(X_1,\dots,X_n)\), \(\hat{ \theta}_2 = T_2(X_1,\dots,X_n)\),所构成的区间\([\hat{\theta}_1,\hat{\theta}_2]\)作为 \(\theta\)的区间估计

定义 2.1 假设总体\(X\)的分布是\(F_\theta\). 如果存在\(\theta_1\neq\theta_2\)使得\(F_{\theta_1}=F_{\theta_2}\),则称\(\theta\)对于\(X\)是不可识别的(unidentifiable)。

例 2.1 \(X\sim N(\mu_1+\mu_2,\sigma^2)\), 则不难发现\(\theta=(\mu_1,\mu_2,\sigma^2)\)是不可识别的。估计参数\(\theta\)是没有统计意义的。

考虑一个测量问题(如温度、血压、距离等),记测量对象真实值为\(\theta\)\(n\)次测量数据为\(x_1,x_2,\dots,x_n\)。相信大部分人都会用样本均值

\[\bar x = \frac{x_1+\dots+x_n}{n}\] 来估计\(\theta\)。但是,为什么用样本均值估计是合理的呢?这是最好的估计方式吗?这些是本章重点讨论的问题。

2.1.1 矩估计法

矩估计的想法来源于大数定理。如果总体\(X\)存在\(k\)阶矩,对任意\(\epsilon>0\),

\[ \lim_{n\to \infty} P(|\frac 1 n\sum_{i=1}^n X_i^k-E[X^k]|\ge \epsilon )=0. \]

这说明,当样本容量\(n\)较大时,样本\(k\)阶矩与总体\(k\)阶矩差别很小。矩法估计就是用样本\(k\)阶矩代替总体的\(k\)阶矩。通常用\(\hat{\theta}_M\)表示。一般步骤如下:

  • 列出估计式\(E[X^k]=g_k(\theta_1,\dots,\theta_m),\ k=1,\dots,m.\)

  • 求解关于估计量的方程组\(\theta_k = \theta_k(E[X^1],\dots,E[X^m])\)

  • \(M_k=\frac 1 n\sum_{i=1}^n X_i^k\)替代\(E[X^k]\)得到矩估计\(\hat\theta_k = \theta_k(M_1,\dots,M_m)\)

例 2.2 求总体\(X\)的期望\(\mu=E[X]\)与方差\(\sigma^2=Var[X]\)的矩估计。

解. (1)列出估计式

\[ \begin{cases} E[X] &= \mu\\ E[X^2] &= \mu^2+\sigma^2 \end{cases} \]

(2)求解关于估计量的方程组

\[ \begin{cases} \mu &= E[X]\\ \sigma^2 &= E[X^2]-(E[X])^2 \end{cases} \]

所以,\(\hat{\mu}_M = \bar X\), \(\hat{\sigma}^2_M = \frac{1}{n}\sum_{i=1}^n X_i^2-(\bar X)^2 = S_n^2.\)

:不难证明,总体的各阶中心矩的矩估计就是样本各阶中心矩。

例 2.3 设总体\(X\sim U[a,b]\), 求\(a,b\)的矩估计。

解. 易知,\(E[X]=(a+b)/2,\ Var[X]= (b-a)^2/12\).

所以,

\[ \begin{cases} a &= E[X]-\sqrt{3Var[X]}\\ b &= E[X]+\sqrt{3Var[X]} \end{cases} \]

\[ \begin{cases} \hat a_M &= \bar{X}-\sqrt{3}S_n\\ \hat b_M &= \bar{X}+\sqrt{3}S_n \end{cases} \]

例 2.4 设总体\(X\)的分布密度为

\[ f(x)=\frac{\theta}{2}e^{-\theta|x|},\ x\in\mathbb{R}, \theta>0. \]

\(\theta\)的矩估计。

解. \[ E[X]= 0,\ E[X^2]=\int_{-\infty}^{\infty}x^2\frac{\theta}{2}e^{-\theta|x|}d x=\theta\int_{0}^{\infty}x^2e^{-\theta x}d x=\frac{2}{\theta^2} \]

\[\hat{\theta}_M=\sqrt{\frac{2n}{\sum_{i=1}^n X_i^2}}.\]

除外,还可以由\(E[|X|]=1/\theta\)得到另一种矩估计。

2.1.2 最大似然估计法

最大似然估计法最早由高斯(C.F.Gauss)提出,后来被 Fisher完善。最大似然估计这一名称也是Fisher给的。这是一个目前仍得到广泛应用的方法。它是建立在最大似然原理基础上的一个统计方法。

最大似然原理:最先出现的是概率最大的

例 2.5 设有外形完全相同的两个箱子,甲箱中有99个白球和1个黑球,乙箱中有99个黑球和1个白球,今抽取一箱并从中随机抽取一球,结果取得白球,问这球是从哪个箱子中取出?

定义 2.2 假设总体\(X\)为离散随机变量,其分布函数记为\(P(X=x)=f(x;\theta)\),与参数\(\theta\)相关。设\(X_1,\dots,X_n\)为其样本,\(x_1,\dots,x_n\)为该样本的观测值。样本的似然函数(likelihood function)定义为观测到样本\(x_1,\dots,x_n\)的概率

\[L(x_1,\dots,x_n;\theta)=P(X_1=x_1,\dots,X_n=x_n)=\prod_{i=1}^{n}f(x_i;\theta). \]

固定参数\(\theta\),似然函数\(L(x_1,\dots,x_n;\theta)\)为样本的概率质量函数(Probability Mass Function, PMF)。另一方面,给定样本观测值\(x_1,\dots,x_n\),似然函数\(L(x_1,\dots,x_n;\theta)\)是一个关于\(\theta\)的函数,其中\(\theta\in\Theta\),有时简记为\(L(\theta)\)

例 2.6 设总体\(X\sim B(1,p)\), 从中抽取样本的观测值为\(1,1,0,0,1\). 不难计算似然函数为

\[L(p)=p^3(1-p)^2,\ p\in (0,1).\]

图像如下

对于该数据,不同\(p\)的值得到不同的概率。现在要估计\(p\)的值,一种合理的方式是找出使得该概率最大对应\(p\)的值作为估计值。这就是最大似然估计的核心思想。通过简单的计算,可以发现\(L(p)\)的最大值点发生在\(p=0.6\). 因此,\(0.6\)可以作为\(p\)的估计值。

更一般地,给定样本观测值\((x_1,\dots,x_n)\), 记\(L(x_1,\dots,x_n;\theta)\)的最大值点为\(\theta=T(x_1,\dots,x_n)\). 则\(\theta\)的最大似然估计量(MLE, maximum likelihood estimator)为 \[\hat{\theta}_L=T(X_1,\dots,X_n).\]

对于上述例子,如果观测数据是\(x_1,\dots,x_n\),似然函数则为

\[ L(x_1,\dots,x_n;p)=\prod_{i=1}^{n}p^{x_i}(1-p)^{1-x_i}=p^{\sum_{i=1}^nx_i}(1-p)^{n-\sum_{i=1}^nx_i}.\]

\(y=\sum_{i=1}^nx_i\)。为了便于计算,对似然函数取对数变换,得到对数似然函数为:

\[\ln L = y \ln p + (n-y)\ln (1-p).\] 对数似然函数的极大值点与似然函数的极大值点一致。故求其求导可得,对数似然方程为:

\[\frac{d \ln L}{d p} = y/p - (n-y)/(1-p)=0.\] 解得\(p= y/n=\frac{1}{n}\sum_{i=1}^nx_i\). 因为\(\frac{d^2\ln L}{d p^2}<0\), 所以\(p= y/n\)是极大值。最大似然估计量为\(\hat{p}_L = \bar X.\)

如果\(X\)是连续的总体,似然函数该如何定义?此时,若沿用上述定义,由于连续型随机变量在某点发生的概率为零,则有\(P(X_1=x_1,\dots,X_n=x_n)=0\). 然而,在该点一个领域的概率不为零。不妨考虑样本落在观测值点一个小邻域的概率。记\(O(x,\delta) = (x-\delta,x+\delta)\)\(x\)\(\delta\)邻域,则当\(\delta\)比较小时,

\[P(X\in O(x,\delta)) = F(x+\delta;\theta)-F(x-\delta;\theta) \approx f(x;\theta)\delta,\] 其中\(F\)\(f\)分别为X的分布函数和密度函数。现考虑样本落在\(x_1,\dots,x_n\)的附近的概率,

\[P(X_1\in O(x_1,\delta_1),\dots,X_n\in O(x_n,\delta_n))\approx (\prod_{i=1}^n\delta_i)\prod_{i=1}^n f(x_i;\theta),\] 其中\(\delta_i\)为比较小的常数。从中可以看出,该概率的大小与\(\prod_{i=1}^n f(x_i;\theta)\)相关。我们把它定义成连续总体下样本的似然函数,即

\[L(x_1,\dots,x_n;\theta)=\prod_{i=1}^n f(x_i;\theta).\] 此时,样本的似然函数为样本的联合密度函数。

更一般的情形,样本不一定是独立同分布,此时似然函数同样可以定义为该样本的联合密度函数。本质上,似然函数是刻画样本在给定观测值处的“可能性”。PMF/PDF用来衡量这种“可能性”。

最大似然估计的一般步骤归纳如下

第一步:写出似然函数\(L(x_1,\dots,x_n;\theta)\)

第二步:若似然函数\(L\)\(\theta\)的可微函数,则最大值必然满足似然方程

\[\frac{d L}{d \theta}=0\]

解出\(\theta\), 并验证其是否是极大值:\[\frac{d^2 L}{d \theta^2}<0.\]

注1:为方便求导,一般求对数似然函数\(\ln L(x_1,\dots,x_n;\theta)\)求极大值点

注2:若有多个参数\(\theta_1,\dots,\theta_m\),对每个变量求偏导,联立\(m\)个方程求解

注3:(对数)似然方程的解称为“驻点”(stationary point),可能为(局部)极大或者极小值点,也可能为鞍点(saddle point),为了进一步区分需要求 Hessian矩阵并分析其正定性。

三种类型的驻点

图 2.1: 三种类型的驻点

例 2.7 设总体\(X\sim N(\mu,\sigma^2)\), 从中抽取样本\(X_1,\dots,X_n\)的观测值为\(x_1,\dots,x_n\). 求参数\(\mu,\sigma^2\)的最大似然估计。

解. 似然函数为

\[ L(x_1,\dots,x_n;\mu,\sigma^2)=\prod_{i=1}^{n}f(x_i)=\prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}\sigma}e^{-(x_i-\mu)^2/(2\sigma^2)}\]

\(\theta_1=\mu,\theta_2=\sigma^2\), 对数似然函数为:

\[\ln L = (n/2)\ln (2\pi)-(n/2)\ln\theta_2-\frac{\sum_{i=1}^n(x_i-\theta_1)^2}{2\theta_2}\]

对数似然方程组为:

\[ \begin{cases} \frac{\partial \ln L}{\partial \theta_1} &=\frac{\sum_{i=1}^n(x_i-\theta_1)}{\theta_2}=0\\ \frac{\partial \ln L}{\partial \theta_2} &=-\frac{n}{2\theta_2}+\frac{\sum_{i=1}^n(x_i-\theta_1)^2}{2\theta_2^2}=0 \end{cases} \]

解得\(\hat{\mu}_L=\bar X,\ \hat{\sigma}^2_L = S_n^2\). (可以验证二阶导函数非正定,即取得极大值。)
例 2.8 设总体\(X\sim U[a,b]\), 从中抽取样本\(X_1,\dots,X_n\)的观测值为\(x_1,\dots,x_n\). 求参数\(a,b\)的最大似然估计。
解. 似然函数为

\[ L(x_1,\dots,x_n;a,b)=\frac{1}{(b-a)^n}\prod_{i=1}^{n} 1\{a\le x_i\le b\}\]

注意到\(L\)关于\(a,b\)不可微。容易观察到,当\(a=\min_{i=1,\dots,n}\{x_i\},\ b=\max_{i=1,\dots,n}\{x_i\}\)\(L\)取得最大值。故

\[\hat{a}_L = X_{(1)},\ \hat{b}_L = X_{(n)}.\]

关于最大似然估计的一些说明:

  1. 最大似然估计的不变性:如果\(\hat{\theta}\)\(\theta\)的最大似然估计,则对任一函数\(g(\theta)\), 其最大似然估计为\(g(\hat{\theta})\).

  2. 当分布中有多余的参数或者数据为截尾或缺失时,似然函数的求极大值比较困难。针对这种问题,文献

Dempster, A.P.; Laird, N.M.; Rubin, D.B. (1977). Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society, Series B. 39 (1): 1–38. (cited by 54539, 2018/8/18)

提出了一种有效的Expectation–Maximization (EM)算法。

2.1.3 矩估计与最大似然估计的对比

  1. 矩估计法(也称数字特征法)
  • 直观意义比较明显,但要求总体\(k\)阶矩存在。
  • 缺点是不唯一,此时尽量使用样本低阶矩。
  • 观测值受异常值影响较大,不够稳健,实际中避免使用样本高阶矩。
  • 估计值可能不落在参数空间
  1. 极大似然估计法
  • 具有一些理论上的优点(不变性、渐近正态性)
  • 缺点是如果似然函数不可微,没有一般的求解法则。
分布名称 记号 期望 方差 矩估计 极大似然估计
0-1分布 \(B(1,p)\) \(p\) \(pq\) \(\hat p_M=\bar{X}\) \(\hat p_L=\bar{X}\)
泊松分布 \(Pois(\lambda)\) \(\lambda\) \(\lambda\) \(\hat{\lambda}_M=\bar{X}\) \(\hat{\lambda}_L=\bar{X}\)
几何分布 \(Geo(p)\) \(1/p\) \(q/p^2\) \(\hat p_M=1/\bar{X}\) \(\hat p_L=1/\bar{X}\)
均匀分布 \(\mathbb{U}[a,b]\) \((a+b)/2\) \((b-a)^2/12\) \(\hat{a}_M=\bar X-\sqrt{3}S_n\) \(\hat{a}_L=X_{(1)}\)
\(\hat{b}_M=\bar X+\sqrt{3}S_n\) \(\hat{b}_L=X_{(n)}\)
指数分布 \(Exp(\lambda)\) \(1/\lambda\) \(1/\lambda^2\) \(\hat{\lambda}_M=1/\bar{X}\) \(\hat{\lambda}_L=1/\bar{X}\)
正态分布 \(N(\mu,\sigma^2)\) \(\mu\) \(\sigma^2\) \(\hat{\mu}_M=\bar X\) \(\hat{\mu}_L=\bar X\)
\(\hat{\sigma}^2_M = S_n^2\) \(\hat{\sigma}^2_L = S_n^2\)

2.1.4 混合正态分布的参数估计

假设总体\(X\)的分布为:以概率\(\lambda\)服从\(N(\mu_1,\sigma_1^2)\), 以概率\(1-\lambda\)服从\(N(\mu_2,\sigma_2^2)\)。该混合分布的密度函数为

\[f(x;\lambda,\mu_1,\sigma_1^2,\mu_2,\sigma_2^2)=\frac{\lambda}{\sqrt{2\pi}\sigma_1}e^{-\frac{(x-\mu_1)^2}{2\sigma_1^2}}+\frac{1-\lambda}{\sqrt{2\pi}\sigma_2}e^{-\frac{(x-\mu_2)^2}{2\sigma_2^2}}.\]

样本似然函数为:

\[L(\lambda,\mu_1,\sigma_1^2,\mu_2,\sigma_2^2) = \prod_{i=1}^n \left[\frac{\lambda}{\sqrt{2\pi}\sigma_1}e^{-\frac{(x_i-\mu_1)^2}{2\sigma_1^2}}+\frac{1-\lambda}{\sqrt{2\pi}\sigma_2}e^{-\frac{(x_i-\mu_2)^2}{2\sigma_2^2}}\right].\]

不难发现,似然函数是无界的(见习题),所以最大似然估计不存在。然而,如果我们有先验的信息:\(\sigma_1=k\sigma_2\),其中\(k\)是已知的数(比如1)。似然函数则可表示为:

\[L(\lambda,\mu_1,\mu_2,\sigma_2^2) =\prod_{i=1}^n\left[ \frac{\lambda}{\sqrt{2\pi}k\sigma_2}e^{-\frac{(x_i-\mu_1)^2}{2k^2\sigma_2^2}}+\frac{1-\lambda}{\sqrt{2\pi}\sigma_2}e^{-\frac{(x_i-\mu_2)^2}{2\sigma_2^2}}\right]. \]

此时,似然函数是有界的,所以最大似然估计存在。

\(Y_i=1\)表示\(X_i\)来自\(N(\mu_1,\sigma_1^2)\)分布,\(Y_i=2\)表示\(X_i\)来自\(N(\mu_2,\sigma_2^2)\)分布。假设我们可以观测\(Y_i\)的值,基于样本\((X_i,Y_i),i=1,\dots,n\),我们可以得到\(\lambda,\mu_1,\sigma_1^2,\mu_2,\sigma_2^2\)的最大似然估计。

\(I = \{i=1,\dots,n|y_i=1\}\), 则似然函数为

\[\begin{align*} L(\lambda,\mu_1,\sigma_1^2,\mu_2,\sigma_2^2)&=\prod_{i\in I} \frac{\lambda}{\sqrt{2\pi}\sigma_1}e^{-\frac{(x_i-\mu_1)^2}{2\sigma_1^2}} \prod_{i\notin I}\frac{1-\lambda}{\sqrt{2\pi}\sigma_2}e^{-\frac{(x_i-\mu_2)^2}{2\sigma_2^2}}\\ &=\lambda^{|I|}(1-\lambda)^{n-|I|}\prod_{i\in I} \frac{1}{\sqrt{2\pi}\sigma_1}e^{-\frac{(x_i-\mu_1)^2}{2\sigma_1^2}} \prod_{i\notin I}\frac{1}{\sqrt{2\pi}\sigma_2}e^{-\frac{(x_i-\mu_2)^2}{2\sigma_2^2}}. \end{align*}\]

则只需分别求出上式中三部分的最大值点即可。对于第一部分,不难发现\(\hat\lambda = \frac 1 n\sum_{i=1}^n1\{Y_i=1\}\). 后两部分等价于求样本为\(\{X_i,i\in I\}\)\(\{X_i,i\notin I\}\)时,对应正态总体的最大似然估计,所以最大似然估计分别为

\[\hat{\mu}_1=\frac{1}{|I|}\sum_{i\in I} X_i=\frac{\sum_{i=1}^nX_i\cdot1\{Y_i=1\}}{\sum_{i=1}^n1\{Y_i=1\}},\]

\[\hat{\sigma_1^2} = \frac{1}{|I|} \sum_{i\in I}(X_i-\hat{\mu}_1)^2=\frac{\sum_{i=1}^n(X_i-\hat\mu_1)^2\cdot1\{Y_i=1\}}{\sum_{i=1}^n1\{Y_i=1\}},\]

\[\hat{\mu}_2=\frac{\sum_{i=1}^nX_i\cdot1\{Y_i=2\}}{\sum_{i=1}^n1\{Y_i=2\}},\ \hat{\sigma_2^2} =\frac{\sum_{i=1}^n(X_i-\hat\mu_2)^2\cdot1\{Y_i=2\}}{\sum_{i=1}^n1\{Y_i=2\}}.\] 上述三式要求分母,否则相应部分的估计量可以为任意常数。

上述问题不难推广到\(K\ge 3\)个不同正态分布的混合的情形。然而,大部分问题,\(Y_i\)是不可观测的(即无标签)。Kiefer (1978) 证明了样本\(X_i\)的似然方程的某个解(对应局部极大值)也是有效的估计量,同样会收敛到真值。这说明了,似然方程的极大值点也可以作为一个有效的估计量。如何找到似然方程的极大值点?我们可以利用EM算法找出似然函数局部极大值点。参考:

Andrew Ng’s lecture notes 1

Andrew Ng’s lecture notes 2

另一方面,对于混合正态分布,我们可以通过矩法得到五个参数的估计量。Cohen (1967)给出了矩法估计的一般公式,转化成求一个9次多项式方程的负根。

2.1.5 EM算法

Expectation-Maximization (EM) 算法是由Dempster et al. (1977)提出。该算法的推导用到Jensen不等式。

定理 2.1 (Jensen不等式) \(A\)\(\mathbb{R}^k\)中凸集,\(f(x)\)\(A\)上凸函数,即对任意\(\lambda\in[0,1], x,y\in A\),恒有

\[ f(\lambda x+(1-\lambda)y)\le \lambda f(x)+(1-\lambda)f(y). \] 如果\(k\)维随机向量\(X\)满足\(P(X\in A)=1\),则有

\[f(E[X])\le E[f(X)].\]

进而,如果\(f\)为严格凸函数,\(f(E[X])=E[f(X)]\)当且仅当\(P(X=E[X])=1\),即\(X\)为常数向量。 换言之,如果\(f\)为严格凸函数且\(X\)不为常数向量,则\(f(E[X])<E[f(X)]\).

如果\(f\)是凹函数,则定理中的不等式变号。注意到\(\ln x\)\((0,\infty)\)上的严格凹函数,应用Jensen不等式得到下面一个结果。

例 2.9 如果\(X\)为一个取值为正的非常数随机变量,则有\(E[\ln X]<\ln (E[X])\).

EM算法用于求解不完备数据的极大似然估计。 设完备数据为\(x=(y,z)\),其中\(y\)为观测数据(向量),\(z\)不可观测,称为潜变量。直接对观测数据求最大似然估计比较困难。 注意到观测数据的似然函数为

\[\begin{align*} L(\theta|y) &= \int p(y,z|\theta)dz=\int q(z|\eta)\frac{p(y,z|\theta)}{q(z|\eta)}dz\\ &=E_{q(z|\eta)}\left[\frac{p(y,z|\theta)}{q(z|\eta)}\right], \end{align*}\]

其中\(p(y,z|\theta)\)为完全数据的密度函数, \(q(z|\eta)\)为密度函数且满足:如果\(q(z|\eta)=0\),则\(p(y,z|\theta)=0\)。对数似然函数

\[\ell(\theta|y)= \ln \left(E_{q(z|\eta)}\left[\frac{p(y,z|\theta)}{q(z|\eta)}\right]\right)\ge E_{q(z|\eta)}\left[\ln \left(\frac{p(y,z|\theta)}{q(z|\eta)}\right)\right]=:Q(q(z|\eta),\theta).\]

上式用到Jensen不等式,当且仅当\(p(y,z|\theta)/q(z|\eta)\)为常数(即不依赖\(z\))时,上式等号成立,由于\(\int q(z|\eta)dz =1\),不难得到

\[q(z|\eta)=\frac{p(y,z|\theta)}{p(y|\theta)}=p(z|y,\theta)=:q^*(z|\theta).\] 此时,\(\ell(\theta|y) = Q(q^*(z|\theta),\theta).\) 否则\(\ell(\theta|y)>Q(q(z|\eta),\theta)\). 假设第\(t\)步的估计值为\(\theta_t\),构造如下迭代算法

\[\theta_{t+1} = \arg\max_{\theta\in\Theta} Q(q^*(z|\theta_t),\theta).\]

注意到,\(\ell(\theta_{t+1}|y)\ge Q(q^*(z|\theta_t),\theta_{t+1})\ge Q(q^*(z|\theta_t),\theta_{t})=\ell(\theta_t|y)\). 如果\(p(y,z|\theta_{t+1})/p(y,z|\theta_t)\)\(z\)相关(即不为常数向量),则有\(\ell(\theta_{t+1}|y)>\ell(\theta_t|y)\). 所以,该迭代算法使得似然函数单调递增。如果\(\theta_t\)收敛到\(\theta^*\),那么在满足一定条件下,\(\theta^*\)为似然函数的驻点。

注意到

\[Q(q^*(z|\theta_t),\theta)=E_{p(z|y,\theta_t)}\left[\ln \left(p(y,z|\theta)\right)\right]-E_{p(z|y,\theta_t)}\left[\ln \left(p(z|y,\theta_t)\right)\right].\] 上式等式最后一项与\(\theta\)无关,所以优化问题等价于

\[\theta_{t+1} = \arg\max_{\theta\in\Theta} E_{p(z|y,\theta_t)}\left[\ln \left(p(y,z|\theta)\right)\right].\]

EM算法包含两步:

  1. 第一步求期望\(E_{p(z|y,\theta_t)}\left[\ln \left(p(y,z|\theta)\right)\right]\), 该期望称为预期的对数似然函数。

  2. 第二步则求预期的对数似然函数最大值。

以下定理保证该迭代算法收敛到似然函数的驻点。

定理 2.2 如果\(E_{p(z|y,\theta_t)}\left[\ln \left(p(y,z|\theta)\right)\right]\)关于\(\theta\)\(\theta_t\)连续,则\(\theta_{t}\)收敛到似然函数\(L(\theta|y)\)的某一驻点(局部极大值点或者鞍点)。
证明. 证明见Lehmann and Casella的《点估计理论》P460。

例 2.10 假设完备数据\((y,z)\)服从指数型分布族

\[p(y,z|\theta) = c(\theta)\exp\left(\sum_{i=1}^k T_i(y,z)c_i(\theta)\right)h(y,z).\]

则对数似然函数为: \[ \ln p(y,z|\theta) = \ln c(\theta)+\sum_{i=1}^k T_i(y,z)c_i(\theta)+\ln h(y,z). \]

\[\begin{align*} \theta_{t+1} &= \arg\max_{\theta\in\Theta} E_{p(z|y,\theta_t)}\left[\ln c(\theta)+\sum_{i=1}^k T_i(y,z)c_i(\theta)+\ln h(y,z)\right]\\ &=\arg\max_{\theta\in\Theta} \left\lbrace\ln c(\theta)+\sum_{i=1}^k c_i(\theta)E_{p(z|y,\theta_t)}[T_i(y,z)]\right\rbrace. \end{align*}\]

要求解该优化问题,只须求\(E_{p(z|y,\theta_t)}[T_i(y,z)]=E_{\theta_t}[T_i(y,z)|y]\)的表达式即可。

例 2.11 再次考虑混合正态分布。设完全数据\(x=((y_i,z_i),i=1,\dots,n)=(y,z)\),其中\(z_i=1\)或者\(2\), 当\(z_i=1\)是表示\(Y_i\sim N(\mu_1,\sigma_1^2)\),当\(z_i=2\)是表示\(Y_i\sim N(\mu_2,\sigma_2^2)\). 这里分类变量\(z_i\)是不可以观测。注意到直接对观测数据\(y\)构建似然函数,然后求极大值是非常困难的。下面我们通过EM算法求解。记\(\theta=(\mu_1,\sigma_1^2,\mu_2,\sigma_2^2,\lambda)\). 完全数据的对数似然函数为

\[ \ln p(y,z|\theta)=\sum_{i=1}^n \ln p(y_i,z_i|\theta) = \sum_{i=1}^n \ln p(y_i|z_i,\theta)+\ln p(z_i|\theta). \]

首先计算

\[\begin{align*} &E_{p(z|y,\theta_t)}[\ln p(y,z|\theta)]=E_{p(z|y,\theta_t)}\left[\sum_{i=1}^n \ln p(y_i|z_i,\theta)+\ln p(z_i|\theta)\right]\\ &=\sum_{i=1}^n \sum_{j=1}^2 P(z_i=j|y_i,\theta_t)[-\frac 1 2\ln (2\pi)-\ln\sigma_j-\frac{(x_i-\mu_j)^2}{2\sigma_j^2}+\ln p(z_i=j|\theta)]\\ &=-\frac 12\sum_{i=1}^n \sum_{j=1}^2 w_{ij}[\ln \sigma_j^2+(y_i-\mu_j)^2/\sigma_j^2]+\sum_{i=1}^n [w_{i1}\ln\lambda+w_{i2}\ln(1-\lambda)]\\ &=s_1(\mu_1,\mu_2,\sigma_1^2,\sigma_2^2)+s_2(\lambda). \end{align*}\]

其中,

\[\begin{align*} w_{ij}&= P(z_i=j|y_i,\theta_t)=\frac{p(y_i|z_i=j,\theta_t)p(z_i=j|\theta_t)}{\sum_{j=1}^2p(y_i|z_i=j,\theta_t)p(z_i=j|\theta_t)}\\ &=\frac{\phi(y_i;\mu_{j,t},\sigma_{j,t}^2)(\lambda_{t}1\{j=1\}+(1-\lambda_{t})1\{j=2\})}{\lambda_{t}\phi(y_i;\mu_{1,t},\sigma_{1,t}^2)+(1-\lambda_{t})\phi(y_i;\mu_{2,t},\sigma_{2,t}^2)}, \end{align*}\]

\[s_1(\mu_1,\mu_2,\sigma_1^2,\sigma_2^2)=-\frac 12\sum_{i=1}^n \sum_{j=1}^2 w_{ij}[\ln \sigma_j^2+(y_i-\mu_j)^2/\sigma_j^2],\]

\[s_2(\lambda)=\sum_{i=1}^n [w_{i1}\ln\lambda+w_{i2}\ln(1-\lambda)],\]

\(\phi(x;\mu,\sigma^2)\)\(N(\mu,\sigma^2)\)的密度函数。 注意到\(w_{ij}\)\(\theta\)无关,可视为常数。对\(E_{p(z|y,\theta_t)}[\ln p(y,z|\theta)]\)求最大值,等价于分别对\(s_1,s_2\)两部分求最大值。先考虑对第一部分\(s_1(\mu_1,\mu_2,\sigma_1^2,\sigma_2^2)\)求最大值. 对\(\mu_j\)求偏导得到

\[\frac{\partial s_1(\mu_1,\mu_2,\sigma_1^2,\sigma_2^2)}{\partial \mu_j}=\sum_{i=1}^n w_{ij}[(y_i-\mu_j)/\sigma_j^2]=0.\] 于是有

\[\mu_{j,t+1}=\frac{\sum_{i=1}^nw_{ij}y_{i}}{\sum_{i=1}^nw_{ij}}, j=1,2.\]

对求\(\sigma^2_j\)偏导得到

\[\frac{\partial s_1(\mu_1,\mu_2,\sigma_1^2,\sigma_2^2)}{\partial \sigma^2_j}=-\frac 12\sum_{i=1}^n w_{ij}\left[\frac 1{\sigma_j^2}-\frac{(y_i-\mu_j)^2}{\sigma_j^4}\right]=0.\]

于是有

\[\sigma_{j,t+1}^2=\frac{\sum_{i=1}^nw_{ij}(y_i-\mu_{j,t+1})^2}{\sum_{i=1}^nw_{ij}}.\]

考虑第二部分\(s_2(\lambda)\)求最大值. 对\(\lambda\)求导得,

\[s_2'(\lambda)=\sum_{i=1}^n\left(\frac{w_{i1}}{\lambda}-\frac{w_{i2}}{1-\lambda}\right)=0.\]

于是有

\[\lambda_{t+1}=\frac{\sum_{i=1}^nw_{i1}}{\sum_{i=1}^n\sum_{j=1}^2w_{ij}}=\frac{1}{n}\sum_{i=1}^nw_{i1}.\]

EM算法求解混合正态分布,红色为真实值,n=10000

图 2.2: EM算法求解混合正态分布,红色为真实值,n=10000

2.2 估计的优良性标准

前面介绍了点估计中两种经典的方法,对于同一个问题,两种方法得到的估计量可能不一样。自然要问,哪种最好? 本节讨论估计量的优良性质。

2.2.1 无偏性

定义 2.3 设总体\(X\sim F(x;\theta),\theta\in \Theta\), \(T(X_1,\dots,X_n)\)\(g(\theta)\)的估计量。

  • 无偏估计量:

\[E[T(X_1,\dots,X_n)]=g(\theta), \forall \theta\in \Theta\]

  • 渐近无偏估计量:
\[\lim_{n\to \infty}E[T(X_1,\dots,X_n)]=g(\theta), \forall \theta\in \Theta\]

无偏性意味着:虽然估计量\(T\)由于随机可能偏离真值\(g(\theta)\), 但取其平均值(期望)却等于\(g(\theta)\). 即没有系统偏差。

  • 样本均值是总体的均值的无偏估计,即\(E[\bar X]=E[X]\)

  • 样本方差是总体方差的渐近无偏估计,即\(\lim_{n\to \infty}E[S_n^{2}]=Var[X]\)

  • 修正样本方差是总体方差的无偏估计,即\(E[S_n^{*2}]=Var[X]\)

如果\(g(\theta)\)存在无偏估计量,则称\(g(\theta)\)是可估的。但注意不是所有的参数估计都存在一个无偏估计量。

例 2.12 假设总体\(X\sim B(k,\theta)\), 其中\(k\ge 1\)已知,\(\theta\in (0,1)\)未知。 证明\(1/\theta\)不存在无偏估计量。

证明. 假设\(T(X_1,\dots,X_n)\)\(1/\theta\)无偏估计量。不妨考虑\(n=1\)情形,

\[E[T] = \sum_{i=0}^k C_k^i \theta^i(1-\theta)^{k-i} T(i)=:h(\theta).\] 显然当\(\theta\to 0\), \(h(\theta)\to T(0)\), 但\(\frac{1}{\theta}\to \infty\). 所以,\(E[T]\neq 1/\theta\).

此外,还应当注意一点是:无偏估计量不一定比有偏的估计量好(如下图所示)。下节给出一个评判标准,并通过一个具体的例子说明有时候有偏估计量更好。

无偏VS有偏

图 2.3: 无偏VS有偏

2.2.2 均方误差

定义 2.4 \(T(X_1,\dots,X_n)\)\(g(\theta)\)的估计量,其均方误差 (mean squared error, MSE)为

\[M_{\theta}(T):=E_\theta[(T(X_1,\dots,X_n)-g(\theta))^2].\]

均方根误差 (root mean squared error, RMSE)为

\[R_{\theta}(T):=\sqrt{E_\theta[(T(X_1,\dots,X_n)-g(\theta))^2]}.\]

注意到:

\[M_{\theta}(T)=(E[T]-g(\theta))^2+Var(T)=\text{偏差}^2+\text{方差}.\]

:如果\(T\)\(g(\theta)\)的无偏估计,则\(M_{\theta}(T)=Var(T)\)

比较两个估计量的优劣

定义 2.5 \(T_1(X_1,\dots,X_n)\)\(T_2(X_1,\dots,X_n)\)都为\(g(\theta)\)的估计量,

  • 如果\(M_{\theta}(T_1)\le M_{\theta}(T_2),\forall \theta\in \Theta\), 则称\(T_1\)不次于\(T_2\)
  • 在此基础上,如果存在一个\(\theta_0\in\Theta\)使得\(M_{\theta_0}(T_1)< M_{\theta_0}(T_2)\), 则称\(T_1\)\(T_2\)有效。
例 2.13 设总体\(X\)的期望\(\mu\)方差为\(\sigma^2\), \(X_1,\dots,X_n\)为其样本(\(n>1\)),证明下列估计量\(\hat{\mu} = \sum_{i=1} C_iX_i\)\(\mu\)的无偏估计的充要条件是\(\sum_{i=1}^nC_i = 1.\) 在满足该条件前提下,\(C_i\)取何值时,\(\hat{\mu}\)的最有效。

解. \(E[\hat{\mu}]=\mu\Leftrightarrow \sum_{i=1}^nC_i = 1\)

\[Var[\hat{\mu}]=\sigma^2\sum_{i=1}^nC_i^2\ge \sigma^2\frac{(C_1+\dots+C_n)^2}{n}=\frac{\sigma^2}{n}.\]

而且唯一的最小值在\(C_i=1/n,i=1,\dots,n\)处取得。

例 2.14 \(X_1,\dots,X_n\)\(N(\mu,\sigma^2)\)分布的样本,参数\(\mu,\sigma^2\)未知。样本方差\(S_n^2\)与修正样本方差\(S_n^{*2}\)作为\(\sigma^2\)的两种估计量,哪个更有效?

解. 由于\(S_n^{*2}\)是无偏的,所以均方误差

\[M(S_n^{*2}) = Var[S_n^{*2}]=\frac{2\sigma^4}{n-1}.\]

\(S_n^{2}\), 其均方误差为

\[\begin{align*} M(S_n^{*2}) &= Var[S_n^{2}]+(E[S_n^2]-\sigma^2)^2\\ &=\frac{2(n-1)\sigma^4}{n^2}+(\frac{(n-1)\sigma^2}{n}-\sigma^2)^2 \\&=\frac{(2n-1)\sigma^4}{n^2}. \end{align*}\]

\[\frac{M(S_n^{*2})}{M(S_n^{2})}=\frac{2n^2}{(n-1)(2n-1)}>1\]

所以,\(S_n^{2}\)\(S_n^{*2}\)有效。

启发:无偏估计量不一定是最有效的。

思考:对于上题,考虑估计量\(T_k=k\sum_{i=1}^n(X_i-\bar X)^2\),其中\(k\)为给定常数。特别地,当\(k=1/n\)时,\(T_k=S_n^2\);当\(k=1/(n-1)\)时,\(T_k=S_n^{*2}\)。样本方差\(S_n^{2}\)是不是比其它的\(T_k\)更有效?如果不是,那么最优的\(k\)是多少?

2.2.3 一致最小方差无偏估计

定义 2.6 如果\(T_0(X_1,\dots,X_n)\)\(g(\theta)\)的无偏估计,如果对于\(g(\theta)\)的任意无偏估计量\(T(X_1,\dots,X_n)\)都有

\[Var[T_0]\le Var[T],\ \forall\theta\in\Theta,\]

则称\(T_0\)\(g(\theta)\)一致最小方差无偏估计量 (uniformly minimum-variance unbiased estimator, UMVUE)。 如果\(Var[T_0]\le Var[T]\)\(\theta=\theta_0\)处成立,则称\(T_0\)\(g(\theta)\)局部最小方差无偏估计量 (locally minimum-variance unbiased estimator, LMVUE)。

例 2.15 \(X\)为取值\(-1,0,1,\dots\)的离散型随机变量,分布为

\[P(X=-1)=p,\ P(X=k)=q^2p^k,\ k=0,1,\dots,\] 其中\(0<p<1,\ q=1-p\)。现考虑用\(X\)来估计\(p\)\(q^2\). 对于估计\(p\),一个简单的无偏估计为

\[T_1 = 1\{X=-1\}.\] 对于估计\(q^2\),一个简单的无偏估计为 \[T_2 = 1\{X=0\}.\] 为求解最小方差估计量,现介绍以下引理。
引理 2.1 \(T_0\)\(g(\theta)\)的任一无偏估计量,则所有的无偏估计量都可以表示为\(T=T_0-U\),其中\(U\)为“零”的无偏估计,即\(E[U]=0\).

对于上例,考虑“零”的无偏估计\(U=U(X)\)。注意到,

\[\begin{align*} E[U]&=\sum_{i=-1}^\infty U(i)P(X=i)\\ &=pU(-1)+q^2\sum_{i=0}^\infty U(i)p^k=0,\ \forall p\in(0,1). \end{align*}\]

\(p\to 0\), 可得\(U(0)=0\). 于是,

\[\sum_{k=1}^\infty U(k)p^{k-1} + \frac{U(-1)}{(1-p)^2}=0.\] 注意到\(1/(1-p)^2=\sum_{k=1}^\infty kp^{k-1}\). 则有

\[\sum_{k=1}^\infty [U(k)+kU(-1)]p^{k-1} =0.\] 由于上式对于任意\(p\in (0,1)\),则所有系数应该为零,即\(U(k)=-kU(-1),k=1,\dots,\infty\). 这意味着,\(U\)为“零”的无偏估计量当且仅当\(U(k)=ak\)对于所有的\(k=-1,0,1,\dots\)和某个\(a\)。要使方差最小等价于最小化二阶矩:

\[E[(T_i-U)^2]=\sum_{k=-1}^\infty P(X=k)[T_i(k)-ak]^2,\ i=1,2.\]

两种情况下,最小值点分别为

\[a_1^* = -\frac{p}{p+q^2\sum_{k=1}^\infty k^2p^k},\ a_2^*=0.\] 由于\(a_2^*\)不依赖\(p\),所以\(T_2^*=T_2-a_2^*X=T_2\)是UMUVE. 但\(a_1^*\)依赖\(p\),所以\(T_1^*=T_1-a_1^*X\)是LMUVE,对于估计\(p\),不存在一致最小方差无偏估计量。

通过这个例子,我们发现对于同一总体,有些参数估计问题存在UMVUE,有些则不存在。那么自然要问,什么情况下所有可估的参数一定存在UMVUE?Blackwell, Rao, Lehmann, Scheffe等统计学家获得了一系列寻求UMVUE的理论和方法。为解决这个问题,下面先引入完全统计量的概念。

定义 2.7 \(T(X_1,\dots,X_n)\)为统计量。如果对任何(Borel可测)函数\(u(\cdot)\), 只要\(E[u(T)]=0\)(对一切\(\theta\))就可以推出\(P(u(T)=0)=1\), 则称统计量\(T\)为参数\(\theta\)完全的统计量

如何理解统计量的完全性?假设非常数统计量\(T\)的分布与参数\(\theta\)无关,则对于所有的函数\(u(x)\), \(E[u(T)]\)的值为常数\(c\),与\(\theta\)无关。这表明\(E[u(T)-c]=0\)对任意的\(\theta\)成立,但由于\(u\)的任意性,\(u(T)-c\)可以不为零(比如取\(u(T)=T\)),故\(T\)不可能为完全统计量。如果统计量\(T\)的分布与参数\(\theta\)无关,我们称此类统计量为辅助统计量。这类统计量不包含参数信息,故认为是“辅助的”。完全性排除了非常数辅助统计量。

例 2.16 \(X_1,\dots,X_n\)是来自\(N(\mu,1)\)的样本。显然\(T_1=X_1-X_2\), \(T_2=X_1-\bar X\)的分布与参数\(\mu\)无关,故它们不是完全统计量。
例 2.17 \(X_1,\dots,X_n\)是来自伯努利分布\(B(1,p)\)的样本。证明\(T=\sum_{i=1}^n X_i\)是参数\(p\)的完全统计量。

证明. 注意到\(T\sim B(n,p)\)。假设对一切\(p\in(0,1)\)都有\(E[u(T)]=0\),则

\[E[u(T)]=\sum_{i=0}^{n}u(i)C_n^ip^i(1-p)^{n-i}= 0.\]\(y=p/1-p\),则对任意\(y\in \mathbb{R}\)都有\(\sum_{i=0}^n C_n^iu(i)y^i=0\)。等式左边为\(y\)的多项式,所以该多项式的所有系数均为0,即\(u(i)=0,\ i=0,\dots,n\). 这意味着\(u(T)=0\). 这就说明\(T\)为完全统计量。

定理 2.3 考虑指数型分布族\(\mathcal{F}=\{f_\theta(x);\theta\in\Theta\}\)中的分布\(f_\theta(x)\)(分布列或者密度函数)都可以表示成如下形式:

\[f_\theta(x)=c(\theta)\exp\{\sum_{j=1}^kc_j(\theta)T_j(x)\}h(x).\] 如果\(\Theta\)有内点,则该分布族的充分统计量

\[\left(\sum_{i=1}^nT_1(x_i),\dots,\sum_{i=1}^nT_k(x_i)\right)\] 是完全的。
证明. 证明略。
定理 2.4 (Black-Lehmann-Scheffe定理) 考虑参数分布族\(\mathcal{F}=\{f_\theta(x);\theta\in\Theta\}\),设\(T(X_1,\dots,X_n)\)为其完全的充分统计量。则所有的可估参数\(g(\theta)\)均存在最小方差无偏估计且可表示为\(T\)的一个函数\(\psi(T)\)(该表示在概率意义是唯一的)。

证明. \(T_1\)\(g(\theta)\)的任意无偏估计量,记\(\psi_1(t) = E[T_1|T=t]\). 由于\(T\)是充分统计量,所以\(\psi_1(t)\)与参数\(\theta\)无关,即\(\psi_1(T)=E[T_1|T]\)为一统计量。由全期望公式知,\(E[\psi_1(T)]=E[E[T_1|T]]=E[T_1]=g(\theta)\). 所以\(\psi_1(T)\)\(g(\theta)\)的无偏估计量。 由全方差公式知,

\[Var[T_1] = E[Var[T_1|T]]+Var[E[T_1|T]]\ge Var[\psi_1(T)].\]

同理,设\(T_2\)\(g(\theta)\)的另一无偏估计量,\(\psi_2(T)=E[T_2|T]\)同样是无偏的,且\(Var[\psi_2(T)]\le Var[T_2].\)\(u(x) = \psi_2(x)-\psi_1(x)\). 则\(E[u(T)]=E[\psi_2(T)]-E[\psi_1(T)]=0\). 由于\(T\)是完全统计量,所以\(P(\psi_2(T)=\psi_1(T))=1\), 这表明\(\psi_2(T)\)\(\psi_1(T)\)在概率意义上是相等的,为UMVUE.

该定理表明如果存在完全充分统计量\(T\),UMVU估计量必然可以表示为该统计量的函数\(\psi(T)\)。根据这个性质可以求解UMVU估计量,即求解方程

\[E[\psi(T)]=g(\theta), \forall \theta\in\Theta.\]

例 2.18 \(X\sim B(1,p)\),其中\(p\in(0,1)\). 求\(g(p)=p(1-p)\)的UMVUE.

解. 由定理2.3知,\(T=\sum_{i=1}^n X_i\)为完全充分统计量。所以,对任意的\(p\in(0,1)\)恒有

\[E[\psi(T)]=\sum_{i=0}^n \psi(i)C_n^i p^i(1-p)^{n-i} = p(1-p).\]

\(\rho = p/(1-p)\), 则\(p = \rho/(1+\rho)\). 上式等价于

\[\sum_{i=0}^n \psi(i)C_n^i \rho^i = \rho(1+\rho)^{n-2}=\sum_{i=1}^{n-1}C_{n-2}^{i-1}\rho^i.\]

比较系数可得

\[\psi(x)=\frac{(n-x)x}{n(n-1)}.\]

所以,\(p(1-p)\)的UMUV估计量为\(n(1-\bar X)\bar X/(n-1)\).

注意到\(X_i\)取值为\(0,1\), 所以\(X_i=X_i^2\). 于是,

\[S_n^2=\frac 1n\sum_{i=1}^n X_i^2-(\bar X)^2=\frac 1n\sum_{i=1}^n X_i-(\bar X)^2=\bar X(1-\bar X).\]

因此,\(S_n^{*2}=nS_n^2/(n-1)=n(1-\bar X)\bar X/(n-1)\)为上述推导的UMUVE.

另外还可以通过求条件期望的方式得到UMVU估计量。假设\(T_1\)\(g(\theta)\)的任意无偏估计量,\(T\)为完全充分统计量。则\(E[T_1|T]\)为UMVU估计量。这种方法的难点在于计算条件期望。

对于正态分布总体\(X\sim N(\mu,\sigma^2)\),我们知道\((\bar X,S_n^2)\)为完全充分统计量。则\(\bar X\)\(\mu\)的UMVUE, \(S_n^{*2}\)\(\sigma^2\)的UMVUE. 下面考虑其它参数的UMVUE.

例 2.19 设总体\(X\sim N(\mu,\sigma^2)\),其CDF记为\(F(x)\)

  1. 给定\(\alpha\in (0,1)\),求分位数\(F^{-1}(\alpha)\)的UMVUE.

  2. 假设\(\sigma=1\), 给定\(u\in \mathbb{R}\),求\(F(u)\)\(F'(u)\)的UMVUE.

解. (1)注意到\(F^{-1}(\alpha)=\mu+u_\alpha \sigma\). 又\(\bar X\)\(\mu\)的无偏估计。现在先求\(\sigma\)的无偏估计。由抽样分布定理知,

\[Y:=\frac{nS_n^2}{\sigma^2}\sim \chi^2(n-1).\]

对任意的\(r>1-n\), 有

\[\begin{align*} E[(S_n/\sigma)^r]&=E[n^{-r/2}Y^{r/2}]=n^{-r/2}\int_0^\infty x^{r/2}\frac{x^{\frac {n-1}2-1}e^{-\frac x2}}{2^{\frac {n-1}2}\Gamma((n-1)/2)}dx\\ &=\frac{2^{\frac {r+n-1}2}\Gamma((r+n-1)/2)}{n^{r/2}2^{\frac {n-1}2}\Gamma((n-1)/2)}\int_0^\infty \frac{x^{(r+n-1)/2-1}e^{-x/2}}{2^{\frac {r+n-1}2}\Gamma((r+n-1)/2)}dx\\ &=\frac{2^{\frac {r+n-1}2}\Gamma((r+n-1)/2)}{n^{r/2}2^{\frac {n-1}2}\Gamma((n-1)/2)}\\ &=\frac{2^{r/2}\Gamma((r+n-1)/2)}{n^{r/2}\Gamma((n-1)/2)}=:K_{n,r}. \end{align*}\]

所以,

\[E[S_n^r]=\left(\frac 2n\right)^{r/2}\frac{\Gamma((r+n-1)/2)}{\Gamma((n-1)/2)}\sigma^r=K_{n,r}\sigma^r.\]

特别地,取\(r=1\),则\(E[S_n]=K_{n,1}\sigma\). 所以\(S_n/K_{n,1}\)\(\sigma\)的无偏估计。由此可得,\(F^{-1}(\alpha)\)的一个无偏估计量为

\[\bar X+u_\alpha S_n/K_{n,1}.\]

该统计量为完全充分统计量\((\bar X,S_n^2)\)的函数,所以是UMVUE.

(2)当\(\sigma\)已知时,\(\bar X\)\(\mu\)的完全充分统计量。令\(T_1=1\{X_1\le u\}\),则有\(E(T_1)=P(X_1\le u)=F(u)\)。所以,\(T_1\)\(F(u)\)的无偏估计。考虑

\[\begin{align*} E[T_1|\bar X=\bar x]&=P(X_1\le u|\bar X=\bar x)=P(X_1-\bar X\le u-\bar x|\bar X=\bar x)\\ &=P(X_1-\bar X\le u-\bar x)=\Phi\left[\sqrt{\frac{n}{n-1}}(u-\bar x)\right]. \end{align*}\]

其中用到\(\bar X\)\(X_1-\bar X\)独立,且\(X_1-\bar X\sim N(0,(n-1)/n)\). 因此,\(\Phi\left[\sqrt{\frac{n}{n-1}}(u-\bar X)\right]\)\(F(u)\)的UMVUE. 不难发现,

\[\frac{\partial }{\partial u}\Phi\left[\sqrt{\frac{n}{n-1}}(u-\bar X)\right]=\sqrt{\frac{n}{n-1}}\phi\left[\sqrt{\frac{n}{n-1}}(u-\bar X)\right]\]

\(F'(u)\)的无偏估计量,从而为UMVUE. 这里\(\Phi\)\(\phi\)分别为标准正态分布的CDF和PDF.

上题第二问中,如果\(\sigma\)未知,经过柯尔莫哥洛夫(1950)研究,UMVUE同样存在,但推导过程比较繁琐,这里省略。结果见陈家鼎等教材P26-27.

UMVUE考虑无偏估计中最好的一种。注意到在正态总体下,总体方差的UMVUE为修正样本方差,修正样本方差的均方误差却大于样本方差的均方误差。这表明如果考虑所有类型的估计量,UMVUE不一定是最好的。那为什么不在所有的估计量中研究所谓最好的估计量呢?是否可以相应地定义“一致最小均方误差估计量”?然而,这样的估计量是不存在的。不妨\(g(\theta)\)不恒为一个常数。假设存在它的估计量\(T_0\),对于所有的估计量\(T\)满足\(M_\theta (T_0)\le M_\theta(T),\forall\theta\in\Theta\),即\(T_0\)为一致最小均方误差估计量。那么,对任意\(\theta_0\in\Theta\),取\(T=g(\theta_0)\),则\(M_{\theta_0} (T_0)\le M_{\theta_0}(T)=0\), 从而有\(P(T_0=g(\theta_0))=1\). 由于\(\theta_0\)的任意性以及\(g(\theta)\)不恒为一个常数,这样的\(T_0\)是不存在的。这也就是为什么我们在无偏估计量中考虑最优性的原因。

定理 2.5 (Cramer-Rao不等式) 设总体\(X\)的密度为\(f(x;\theta)\), 参数\(\theta\in (a,b)\). \(X_1,\dots,X_n\)\(X\)的样本,\(\psi(X_1,\dots,X_n)\)\(g(\theta)\)的一个无偏估计,且满足下列正则性条件:

  • \(X\)的支撑与\(\theta\)无关;
  • \(g'(\theta)\)\(\frac{df(x;\theta)}{d\theta}\)都存在且对一切\(\theta\)

\[\begin{align*} \int_{-\infty}^\infty \frac{df(x;\theta)}{d\theta} d x &= 0,\\ \int_{-\infty}^\infty\frac d{d\theta} L(\vec x;\theta) d \vec x&=0,\\ \frac d{d\theta}\int_{-\infty}^\infty \psi(\vec x) L(\vec x;\theta) d \vec x&=\int_{-\infty}^\infty \psi(\vec x) \frac d{d\theta}L(\vec x;\theta) d \vec x, \end{align*}\]

  • \(I(\theta):=E[(\frac {d\ln f(X;\theta)}{d\theta})^2]>0\),
则有 \[Var_\theta[\psi(X_1,\dots,X_n)]\ge \frac{[g'(\theta)]^2}{nI(\theta)}.\]

证明. \[\begin{align*} g'(\theta) &= dE[\psi(X_1,\dots,X_n)]/d\theta \\ &=\frac{d}{d\theta}\int \psi(\vec x)L(\vec x;\theta)d \vec x\\ &=\int \psi(\vec x)\frac{d}{d\theta}L(\vec x;\theta)d \vec x\\ &=\int [\psi(\vec x)-g(\theta)]\frac{d}{d\theta}L(\vec x;\theta)d \vec x\\ &=\int [\psi(\vec x)-g(\theta)]\frac{d \ln L(\vec x;\theta)}{d\theta} L(\vec x;\theta)d \vec x\\ &=E\left[[\psi(\vec X)-g(\theta)]\frac{d \ln L(\vec X;\theta)}{d\theta}\right]. \end{align*}\]

Cauchy–Schwarz不等式得,

\[\begin{align*} [g'(\theta)]^2 &\le E[[\psi(\vec X)-g(\theta)]^2] E[(\frac{d}{d\theta}\ln L(\vec X;\theta))^2]\\ &=Var[\psi(\vec X)]E[(\sum_{i=1}^n \frac{d\ln f(X_i;\theta)}{d\theta})^2]\\ &=Var[\psi(\vec X)] n I(\theta). \end{align*}\]

其中用到

\[\begin{align*} E\left[\frac{d\ln f(X_i;\theta)}{d\theta}\right]&=\int \frac{\frac{df(x;\theta)}{d\theta}}{f(x;\theta)}f(x;\theta)dx\\&=\int \frac{df(x;\theta)}{d\theta}dx=0. \end{align*}\]

C-R不等式给出无偏估计量方差的下界,如果某个无偏估计量达到这个下界且定理2.5中的条件对所有的无偏估计成立,则可以说明是该无偏估计量是一致最小方差无偏的。然而,有些情况下,C-R不等式的下界不一定达到,见陈家鼎等编著的教材例2.9, p30. \(I(\theta)\)叫做Fisher信息量。离散情形有类似的结论。

有时候定理2.5中的条件并不满足,但C-R不等式的下界还是可以用来刻画模型参数的“可估能力”。该下界越小越容易得到精度更高的估计。

例 2.20 \(X\sim N(\mu,\sigma^2)\), 其中\(\mu\)未知,\(\sigma\)已知。 Fisher信息量为

\[I(\mu) = E\left[(\frac {d\ln f(X;\mu)}{d\mu})^2\right]=\frac 1{\sigma^4}E[(X-\mu)^2]=\frac 1{\sigma^2}.\]

\[Var[\bar X] = \frac{\sigma^2}{n}=\frac{1}{nI(\mu)}\]

样本均值\(\bar X\)的方差达到了C-R不等式的下界。

试证明:若\(\mu\)已知,则\(\sigma^2\)的估计量\(\frac 1n\sum_{i=1}^n(X_i-\mu)^2\)的方差达到了C-R不等式的下界。

2.2.4 统计量的大样本性质

  1. 统计量的相合性(consistency)

(弱)相合估计:称\(T_n(X_1,\dots,X_n)\)\(g(\theta)\)的相合估计,如果对任何 \(\epsilon>0\), 有\[\lim_{n\to\infty}P(|T_n-g(\theta)|\ge \epsilon)=0.\]
也称\(T_n\)依概率收敛到\(g(\theta)\),记为\(T_n\stackrel p\to g(\theta)\).

强相合估计:称\(T_n(X_1,\dots,X_n)\)\(g(\theta)\)的强相合估计,如果

\[P(\lim_{n\to\infty}T_n=g(\theta))=1.\] 也称\(T_n=\)以概率1收到\(g(\theta)\),记为敛\(T_n\stackrel {w.p.1}\to g(\theta)\).

说明

  • 由强大数定理知,矩估计一般是强估计的
  • 最大似然估计在十分广泛的条件下也是有相合性(见下一节)

引理 2.2 如果\(T_n\stackrel p\to a\),且函数\(f(x)\)\(x=a\)处连续,则\(f(T_n)\stackrel p\to f(a)\).

证明. 因为\(f(x)\)\(x=a\)处连续,所有对任意\(\epsilon>0\)存在\(\delta>0\)使得任意\(x\)满足\(|x-a|<\delta\),均有\(|f(x)-f(a)|<\epsilon\). 注意到\(\{|f(T_n)-f(a)|\ge \epsilon\}\subset \{|T_n-a|\ge \delta\}\). 所以

\[0\le \lim_{n\to \infty} P(|f(T_n)-f(a)|\ge \epsilon)\le \lim_{n\to \infty} P(|T_n-a|\ge \delta)=0.\]

这表明\(f(T_n)\stackrel p\to f(a)\).

该引理容易推广到多元连续的情形。

引理 2.3 如果\(T_n^{(i)}\stackrel p\to a_i\), \(i=1,\dots,k\),且\(f(x_1,\dots,x_k)\)\((a_1,\dots,a_k)\)点连续,则\(f(T_n^{(1)},\dots,T_n^{(k)})\stackrel p\to f(a_1,\dots,a_k)\).

例 2.21 设总体\(X\)的期望\(\mu\)方差为\(\sigma^2\), \(X_1,\dots,X_n\)为其样本,证明
  • 样本均值\(\bar X\)\(\mu\)的相合估计量;
  • 样本\(k\)阶原点矩\(M_k\)是总体\(k\)阶原点矩\(E[X^k]\)的相合估计量;
  • 样本方差\(S_n^2\)和修正样本方差\(S_n^{2*}\)都是\(\sigma^2\)的相合估计量。

证明. 由辛钦大数定律知,\(\bar X\stackrel p\to \mu\), \(M_k\stackrel p\to E[X^k]\). 由引理2.3

\[S_n^2 = \frac{1}{n}\sum_{i=1}^nX_i^2-\bar X^2\stackrel p\to E[X^2]-E[X]^2=\sigma^2.\]

同理,\(S_n^{2*}=\frac{n-1}{n}S_n^2\stackrel p\to \sigma^2\).

下面给出判断相合估计的一个常用的充分条件。

定理 2.6 \(T(X_1,\dots,X_n)\)\(g(\theta)\)的估计量。如果

\[\lim_{n\to \infty}E[T(X_1,\dots,X_n)] = g(\theta),\ \lim_{n\to \infty}Var[T(X_1,\dots,X_n)] =0,\]

\(T(X_1,\dots,X_n)\)\(g(\theta)\)的相合估计量。

证明. \(T_n=T(X_1,\dots,X_n)\). 注意到

\[\{|T_n-g(\theta)|\ge \epsilon\}\subset \{|T_n-E[T_n]|\ge \epsilon/2\}\cup \{|E[T_n]-g(\theta)|\ge \epsilon/2\}.\]

对任意\(\epsilon>0\), 存在\(N\), 当\(n\ge N\)时,\(|E[T_n]-g(\theta)|< \epsilon/2\). 此时

\[\{|T_n-g(\theta)|\ge \epsilon \}\subset \{|T_n-E[T_n]|\ge \epsilon/2\}.\]

所以,

\[P(|T_n-g(\theta)|\ge \epsilon)\le P(|T_n-E[T_n]|\ge \epsilon/2)\le \frac{4 Var[T_n]}{\epsilon^2}\to 0.\]

此外还可以用Markov不等式(如果\(X\)为非负随机变量且\(a>0\),则\(P(X\ge a)\le E[X]/a\))证明。

所以,

\[\begin{align*} P(|T_n-g(\theta)|\ge \epsilon)&=P((T_n-g(\theta))^2\ge \epsilon^2)\le \frac{E[(T_n-g(\theta))^2]}{\epsilon^2}\\ &=\frac{Var[T_n]+(E[T_n]-g(\theta))^2}{\epsilon^2}\to 0. \end{align*}\]

  1. 统计量的渐近正态性(asymptotic normality)

定义 2.8 \(T(X_1,\dots,X_n)\)\(\theta\)的估计量。如果存在一个趋于零的正数列\(\sigma_n(\theta)\), 使得\((T-\theta)/\sigma_n(\theta)\)的分布收敛到标准正态分布,则称\(T(X_1,\dots,X_n)\)\(\theta\)渐近正态估计,或称\(T\)具备渐近正态性,记为

\[T\stackrel{\cdot}{\sim} N(\theta, \sigma_n(\theta)^2).\]

渐近正态性在构建参数的渐近置信区间中扮演着非常重要的角色。下面定理给出最大似然估计的渐近正态性。

定理 2.7 \(X\)的密度为\(f(x;\theta)\), 其参数空间\(\Theta\)是非退化区间,且满足下列正则性条件:

  • 对一切\(\theta\in\Theta\), \(\frac{\partial \ln f}{\partial\theta}, \frac{\partial^2 \ln f}{\partial\theta^2}, \frac{\partial^3 \ln f}{\partial\theta^3}\) 都存在

  • 对一切\(\theta\in\Theta\), 有\(|\frac{\partial \ln f}{\partial\theta}|<F_1(x),\ |\frac{\partial^2 \ln f}{\partial\theta^2}|<F_2(x),\ |\frac{\partial^3 \ln f}{\partial\theta^3}|<H(x),\) 其中\(F_1(x),F_2(x)\)在实数轴上可积,且\(\int_{-\infty}^\infty H(x)f(x;\theta)dx<M\), \(M\)\(\theta\)无关。

  • 对一切\(\theta\in\Theta\), 有\(0<I(\theta)=E[(\frac{\partial \ln f}{\partial\theta})^2]<\infty\).

则在参数真值\(\theta\)\(\Theta\)内点的情况下,其似然方程有一个解\(\hat{\theta}_L\)存在,且

\[\hat{\theta}_L\stackrel{p}{\to}\theta,\ \hat{\theta}_L\stackrel{\cdot}{\sim} N(\theta,[nI(\theta)]^{-1}).\]

值得注意的是,最大似然估计的渐近方差为C-R不等式的下界,从这个角度可以说明最大似然估计具有良好性质。证明参考:陈希孺. 概率论与数理统计. 中国科技大学出版社, 1992

2.3 区间估计

2.3.1 区间估计的定义

定义 2.9 设总体\(X\sim F(x;\theta),\ \theta\in \Theta\). 如果统计量\(T_1(X_1,\dots,X_n)\), \(T_2(X_1,\dots,X_n)\)使得对给定的\(\alpha\in(0,1)\)

\[P(T_1\le g(\theta)\le T_2)=1-\alpha,\ \forall \theta\in\Theta,\]

则称随机区间\([T_1,T_2]\)为参数\(g(\theta)\)置信度(置信概率)\(1-\alpha\)置信区间(Confidence Interval)\(T_1,T_2\)分别称为置信下界置信上界

说明:

  • 在一些情况下,定义中的“等式”无解,此时考虑的置信区间\([T_1,T_2]\)应满足

\[P(T_1\le g(\theta)\le T_2)\ge 1-\alpha,\ \forall \theta\in\Theta.\]

  • 这里允许\(T_1=-\infty\)或者\(T_2=\infty\),这两种情况为单侧置信区间,否则称为双侧置信区间。
置信区间示意图

图 2.4: 置信区间示意图

2.3.2 枢轴量法

目标:找到\(g(\theta)\)的区间估计,置信度为\(1-\alpha\).

Step 1: 选择恰当的枢轴量(Pivot quantity)\(G(X_1,\dots,X_n;g(\theta))\),其满足以下性质

  • \(G\)不含有其他未知参数
  • \(G\)的分布确定,即不含未知参数\(\theta\)
  • 一般地,\(G\)是关于参数\(g(\theta)\)的单调函数

Step 2: 求\(a,b\)使得\(P(a\le G\le b)=1-\alpha\)

Step 3: 转化不等式\(a\le G\le b\)为如下形式: \[ T_1 \le g(\theta) \le T_2. \]

例 2.22 若总体为指数分布\(Exp(\lambda)\),求未知参数\(\lambda\)的置信区间。

Step 1: 选择枢轴量

\[G(X_1,\dots,X_n;\lambda) = 2\lambda n\bar X\sim Ga(n,1/2)=\chi^2(2n)\]

Step 2: 求\(a,b\)使得\(P(a\le G\le b)=1-\alpha\),即

\[P(a\le 2\lambda n\bar X\le b)=1-\alpha\]

Step 3: \(\lambda\)的置信区间为\([a/(2n\bar X),b/(2n\bar X)]\).

如何选择\(a,b\)?

  • 平分法:\(a=\chi^2_{\alpha/2}(2n), b=\chi^2_{1-\alpha/2}(2n)\)

  • 最优方案?参考书p35

平分法示意图

图 2.5: 平分法示意图

以下通过R模拟来实现这个过程。

set.seed(0)
# generate data from exponential distribution
lambda <- 2
n <- 1000
X <- rexp(n,lambda)

# find out the confidence interval (CI)
alpha <- 0.05
a <- qchisq(p=alpha/2,df=2*n)
b <- qchisq(p=1-alpha/2,df=2*n)
CI <- c(a/2/sum(X),b/2/sum(X))
cat((1-alpha)*100,"% CI is [",CI[1],", ", CI[2],"]",sep="")
## 95% CI is [1.824, 2.065]
# generate R batches of the data
R <- 100
CIs <- matrix(0,R,3)

for(i in 1:R){
  X <- rexp(n,lambda)
  CIs[i,] = c(a/2/sum(X),1/mean(X),b/2/sum(X))
}
## plot the CIs
plot(0, xlim=c(0, R), ylim=c(min(CIs)-0.02,
      max(CIs)+0.02), type="n",xlab="Sample ID",ylab="")
count <- 0
for (i in 1:nrow(CIs)) {
  if (CIs[i, 1]>lambda | CIs[i, 3]<lambda){
    color = "red"
    count = count +1
    if(CIs[i, 1]>lambda)
      text(i,CIs[i, 3]+0.02,count)
    else
      text(i,CIs[i, 1]-0.02,count)
  }else{
    color = "blue"
  }
  lines(x=rep(i, 2), y=c(CIs[i, 1], CIs[i, 3]))
  points(x=i, y=CIs[i,2], pch=16, col=color)
}
abline(h=lambda,lty = 3,col="red")

2.3.3 单个正态总体的区间估计

设总体\(X\sim N(\mu,\sigma^2)\), 如何找出未知参数\(\mu\)\(\sigma^2\)的置信区间?

  • 已知\(\sigma^2\), 找出\(\mu\)的置信区间

  • 未知\(\sigma^2\), 找出\(\mu\)的置信区间

  • 已知\(\mu\), 找出\(\sigma^2\)的置信区间

  • 未知\(\mu\), 找出\(\sigma^2\)的置信区间

  1. 已知方差,求期望的置信区间

由抽样定理知,\(\bar{X}\sim N(\mu,\sigma^2/n)\). 因此 \(U = \frac{\bar{X}-\mu}{\sigma/\sqrt{n}}\sim N(0,1)\)

\[P\left(a\le \frac{\bar{X}-\mu}{\sigma/\sqrt{n}}\le b\right) = 1-\alpha\]

\(\mu\)的置信度为\(1-\alpha\)的置信区间为

\[\left[\bar{X}-b\frac{\sigma}{\sqrt{n}},\ \bar{X}-a\frac{\sigma}{\sqrt{n}}\right]\]

最优的选择\(b=-a=u_{1-\alpha/2}\), 此时置信区间为:

\[\left[\bar{X}-u_{1-\alpha/2}\frac{\sigma}{\sqrt{n}},\ \bar{X}+u_{1-\alpha/2}\frac{\sigma}{\sqrt{n}}\right]=\bar{X}\pm u_{1-\alpha/2}\frac{\sigma}{\sqrt{n}}\]

  1. 方差未知,求期望的置信区间

由抽样定理知,

\[T = \frac{\bar{X}-\mu}{S_n/\sqrt{n-1}}= \frac{\bar{X}-\mu}{S_n^*/\sqrt{n}}\sim t(n-1)\]

\[P\left(a\le \frac{\bar{X}-\mu}{S_n^*/\sqrt{n}}\le b\right) = 1-\alpha\]

\(\mu\)的置信度为\(1-\alpha\)的置信区间为:

\[\left[\bar{X}-b\frac{S_n^*}{\sqrt{n}},\ \bar{X}-a\frac{S_n^*}{\sqrt{n}}\right]=\left[\bar{X}-b\frac{S_n}{\sqrt{n-1}},\ \bar{X}-a\frac{S_n}{\sqrt{n-1}}\right]\]

最优的选择\(b=-a=t_{1-\alpha/2}(n-1)\), 此时置信区间为:

\[\left[\bar{X}-t_{1-\alpha/2}(n-1)\frac{S_n^*}{\sqrt{n}},\ \bar{X}+t_{1-\alpha/2}(n-1)\frac{S_n^*}{\sqrt{n}}\right].\] 也可以表示为\(\bar{X}\pm t_{1-\alpha/2}(n-1)S_n^*/\sqrt{n}\).

例 2.23 假设OPPO手机充电五分钟通话时间\(X\sim N(\mu,\sigma^2)\). 随机抽取6部手机测试通话时间(单位:小时)为

\[1.6,\ 2.1,\ 1.9,\ 1.8,\ 2.2,\ 2.1,\]
  • 已知\(\sigma^2=0.06\), 求\(\sigma^2\)的置信度为\(95\%\)的置信区间。

  • \(\sigma^2\)未知, 求\(\sigma^2\)的置信度为\(95\%\)的置信区间。

解. 查表知,\(u_{1-\alpha/2}=u_{0.975}=1.96,\ t_{1-\alpha/2}=t_{0.975}=2.5706\). 且\(\bar X = 1.95,\ S_n=0.206\).

  1. \(\left[1.95-1.96\frac{\sqrt{0.06}}{\sqrt{6}},\ 1.95+1.96\frac{\sqrt{0.06}}{\sqrt{6}}\right]=[1.754,\ 2.146]\).

  2. \(\left[1.95-2.5706\frac{0.206}{\sqrt{6-1}},\ 1.95+2.5706\frac{0.206}{\sqrt{6-1}}\right]=[1.713,\ 2.187]\).

一些思考

  • 分析这两种的结果会发现,由同一组样本观察值,按同样的置信概率,对\(\mu\)计算出的置信区间因为\(\sigma\)的是否已知会不一样。这因为:当\(\sigma\)为已知时,我们掌握的信息多一些,在其他条件相同的情况下,对\(\mu\)的估计精度要高一些,即表现为\(\mu\)的置信区间长度要小些。反之,当\(\sigma\)为未知时,对\(\mu\)的估计精度要低一些,即表现为\(\mu\)的置信区间长度在大一些。这是因为当\(n\)比较小时,\(t_{1-\alpha/2}(n-1)>u_{1-\alpha/2}\).

  • 还可以发现,当样本量\(n\)不断增大时,两种情况下的置信区间会慢慢接近。 也就意味着大样本信息可以弥补\(\sigma\)的缺失带来的偏差(大数定律)。

  1. 已知期望,求方差的置信区间

选择枢轴量

\[T =\sum_{i=1}^n\frac{(X_i-\mu)^2}{\sigma^2}\sim \chi^2(n)\]

\[P\left(\chi^2_{\alpha/2}(n)\le\sum_{i=1}^n\frac{(X_i-\mu)^2}{\sigma^2}\le \chi^2_{1-\alpha/2}(n)\right) = 1-\alpha\]

\(\sigma^2\)的置信度为\(1-\alpha\)的置信区间为:

\[\left[\frac{\sum_{i=1}^n(X_i-\mu)^2}{\chi^2_{1-\alpha/2}(n)},\ \frac{\sum_{i=1}^n(X_i-\mu)^2}{\chi^2_{\alpha/2}(n)}\right]\]

  1. 期望未知,求方差的置信区间

选择枢轴量

\[T =\frac{nS_n^2}{\sigma^2}=\sum_{i=1}^n\frac{(X_i-\bar X)^2}{\sigma^2}\sim \chi^2(n-1)\] \[P\left(\chi^2_{\alpha/2}(n-1)\le\sum_{i=1}^n\frac{(X_i-\bar X)^2}{\sigma^2}\le \chi^2_{1-\alpha/2}(n-1)\right) = 1-\alpha\]

\(\sigma^2\)的置信度为\(1-\alpha\)的置信区间为:

\[\left[\frac{\sum_{i=1}^n(X_i-\bar X)^2}{\chi^2_{1-\alpha/2}(n-1)},\ \frac{\sum_{i=1}^n(X_i-\bar X)^2}{\chi^2_{\alpha/2}(n-1)}\right]=\left[\frac{nS_n^2}{\chi^2_{1-\alpha/2}(n-1)},\ \frac{nS_n^2}{\chi^2_{\alpha/2}(n-1)}\right]\]

2.3.4 两个独立正态总体的区间估计

设两个独立总体\(X\sim N(\mu_1,\sigma_1^2)\), \(Y\sim N(\mu_2,\sigma^2)\), 如何找出未知参数\(\mu\)\(\sigma^2\)的置信区间?其中\(X\)的样本为\(X_1,\dots,X_m\), 样本方差为\(S_{1m}^2\); \(Y\)的样本为\(Y_1,\dots,Y_n\), 样本方差为\(S_{2n}^2\)

  • 已知\(\sigma_1^2,\sigma_2^2\), 找出\(\mu_1-\mu_2\)的置信区间

  • 以知\(\sigma_1^2=\sigma_2^2=\sigma^2\), 找出\(\mu_1-\mu_2\)的置信区间

  • 已知\(\mu_1,\mu_2\), 找出\(\sigma_1^2/\sigma_2^2\)的置信区间

  • 未知\(\mu_1,\mu_2\), 找出\(\sigma_1^2/\sigma_2^2\)的置信区间

应用场景

  • 比较男生、女生两个群体的身高/体重/成绩平均水平的差异
  1. 已知方差,求均值差的置信区间

选择枢轴量:

\[U=\frac{(\bar X-\bar Y)-(\mu_1-\mu_2)}{\sqrt{\sigma_1^2/m+\sigma_2^2/n}}\sim N(0,1).\]

\(\mu_1-\mu_2\)的置信度为\(1-\alpha\)的置信区间为:

\[\left[(\bar{X}-\bar{Y})-u_{1-\alpha/2}\sqrt{\frac{\sigma_1^2}m+\frac{\sigma_2^2}n},\ (\bar{X}-\bar{Y})+u_{1-\alpha/2}\sqrt{\frac{\sigma_1^2}m+\frac{\sigma_2^2}n}\right]\]

  1. 已知方差相同,求均值差的置信区间

选择枢轴量:

\[T=\frac{(\bar X-\bar Y)-(\mu_1-\mu_2)}{S_w\sqrt{1/m+1/n}}\sim t(m+n-2).\]

其中\(S_w =\sqrt{(mS_{1m}^2+nS_{2n}^2)/(m+n-2)}.\)

\(t_{1-\alpha/2}(m+n-2)=t_{1-\alpha/2}\)\(\mu_1-\mu_2\)的置信度为\(1-\alpha\)的置信区间为:

\[\left[(\bar{X}-\bar{Y})-t_{1-\alpha/2}S_w\sqrt{\frac 1m+\frac 1n},\ (\bar{X}-\bar{Y})+t_{1-\alpha/2}S_w\sqrt{\frac 1m+\frac 1n}\right]\]

例 2.24 假设OPPO手机充电五分钟通话时间\(X\sim N(\mu_1,\sigma_1^2)\), VIVO手机充电五分钟通话时间\(Y\sim N(\mu_2,\sigma_2^2)\). 随机抽取6部手机测试通话时间(单位:小时)为

\[\text{OPPO}:\ 1.6,\ 2.1,\ 1.9,\ 1.8,\ 2.2,\ 2.1\]

\[\text{VIVO}:\ 1.8,\ 2.2,\ 1.5,\ 1.4,\ 2.0,\ 1.7\]

\(\mu_1-\mu_2\)的置信度为\(95\%\)的置信区间:

  • 已知\(\sigma_1^2 = 0.06,\ \sigma_2^2 = 0.08\).
  • 已知\(\sigma_1^2 =\sigma_2^2\).
解. \(m=n=6\), \(\bar{X}=1.95,\ \bar{Y}=1.77\), \(S_{1m}^2=0.042, S_{2n}^2=0.064, S_w = 0.252.\) 查表知,\(u_{0.975}=1.96\), \(t_{0.975}(10)=2.23\).
  • 第一种情况为\([-0.12,\ 0.48]\)
  • 第二种情况为\([-0.14,\ 0.50]\)
  1. 已知均值,求方差比的置信区间

\[T_1 =\sum_{i=1}^m\frac{(X_i-\mu_1)^2}{\sigma_1^2}\sim \chi^2(m),\ T_2 =\sum_{i=1}^n\frac{(Y_i-\mu_2)^2}{\sigma_2^2}\sim \chi^2(n)\]

\[\frac{T_1/m}{T_2/n}=\frac{\frac 1 m\sum_{i=1}^m(X_i-\mu_1)^2}{\frac 1 n\sum_{i=1}^n(Y_i-\mu_2)^2}\frac{\sigma_2^2}{\sigma_1^2}\sim F(m,n)\]

\(\sigma_1^2/\sigma_2^2\)的置信度为\(1-\alpha\)的置信区间为:

\[\left[\frac{1}{F_{1-\alpha/2}(m,n)}\frac{\frac 1 m\sum_{i=1}^m(X_i-\mu_1)^2}{\frac 1 n\sum_{i=1}^n(Y_i-\mu_2)^2},\ \frac{1}{F_{\alpha/2}(m,n)}\frac{\frac 1 m\sum_{i=1}^m(X_i-\mu_1)^2}{\frac 1 n\sum_{i=1}^n(Y_i-\mu_2)^2} \right]\]

  1. 均值未知,求方差比的置信区间

\[T_1=\frac{(m-1)S_{1m}^{*2}}{\sigma_1^2}=\sum_{i=1}^m\frac{(X_i-\bar X)^2}{\sigma_1^2}\sim \chi^2(m-1)\]

\[T_2=\frac{(n-1)S_{2n}^{*2}}{\sigma_2^2}=\sum_{i=1}^n\frac{(Y_i-\bar Y)^2}{\sigma_2^2}\sim \chi^2(n-1)\]

\[\frac{T_1/(m-1)}{T_2/(n-1)}=\frac{S_{1m}^{*2}}{S_{2n}^{*2}}\frac{\sigma_2^2}{\sigma_1^2}\sim F(m-1,n-1)\] \(\sigma_1^2/\sigma_2^2\)的置信度为\(1-\alpha\)的置信区间为:

\[\left[\frac{1}{F_{1-\alpha/2}(m-1,n-1)}\frac{S_{1m}^{*2}}{S_{2n}^{*2}},\ \frac{1}{F_{\alpha/2}(m-1,n-1)}\frac{S_{1m}^{*2}}{S_{2n}^{*2}} \right]\]

一些说明

  • 枢轴量法的难点在于寻找枢轴量,没有统一的方法。正态总体下的应用应当熟练掌握。

  • 另外一种求置信区间方法叫统计量方法,不作要求,感兴趣陈家鼎等编著的教材pp42-46.

  • “最优的置信区间”是否存在?目前尚缺乏对置信区间的优良性讨论。

2.3.5 非正态总体参数的区间估计

\(\mu=E[X],\sigma^2=Var[X]\)分别为总体\(X\)的期望和方差。 由中心极限定理,

\[\frac{\bar X-\mu}{\sigma/\sqrt{n}}\stackrel{\cdot}\sim N(0,1).\]\(\sigma\)已知时,总体期望\(\mu\)的置信度为\(1-\alpha\)的区间估计可以近似为

\[\left[\bar X-u_{1-\alpha/2}\frac{\sigma}{\sqrt{n}}, \bar X+u_{1-\alpha/2}\frac{\sigma}{\sqrt{n}}\right].\] 如果\(\sigma\)未知,可以用样本标准差\(S_n\)(或者修正样本差\(S_n^*\))替代\(\sigma\)

\[\left[\bar X-u_{1-\alpha/2}\frac{S_n}{\sqrt{n}}, \bar X+u_{1-\alpha/2}\frac{S_n}{\sqrt{n}}\right].\]

2.4 分布的估计

本节考虑分布函数和密度函数的估计,目标是通过样本的观测值构造一种函数来近似这两种函数。本节所介绍的方法不需要知道总体的具体的分布类型,属于非参统计方法。

2.4.1 分布函数的估计

定义 2.10 设总体\(X\)的样本\((X_1,\dots,X_n)\)的一次观测值\((x_1,\dots,x_n)\), 并将它们由小到大排列\(x_{(1)}\le x_{(2)}\le \dots\le x_{(n)}\), 经验分布函数(或称样本分布函数)定义为

\[ F_n(x) =\frac{1}{n}\sum_{i=1}^n 1\{x_i\le x\} = \begin{cases} 0,&\ x<x_{(1)}\\ 1/n,&\ x_{(1)}\le x<x_{(2)}\\ 2/n,&\ x_{(2)}\le x<x_{(3)}\\ &\vdots\\ k/n,&\ x_{(k)}\le x<x_{(k+1)}\\ &\vdots\\ 1,&\ x>x_{(n)}\\ \end{cases}. \]

经验分布函数示意图

经验分布函数的性质 固定的\(x\)\(n\)\(F_n(x)\)表示事件\(\{X\le x\}\)的频率,由强大数定律知,

\[F_n(x)\to P(X\le x)=F(x),\]

\[P\left(\lim_{n\to\infty}F_n(x)=F(x)\right)=1.\]

格里汶科定理给出更强的结果(几乎处处一致收敛):

\[P\left(\lim_{n\to\infty}\sup_{x\in \mathbb{R}}|F_n(x)-F(x)|=0\right)=1.\]

:由此可见,当\(n\)相当大时,经验分布函数\(F_n(x)\)是母体分布函数\(F(x)\)的一个良好近似。数理统计学中一切都以样本为依据,其理由就在于此。

2.4.2 直方图法

只考虑一维连续型总体\(X\sim f(x)\)。设\(X_1,\dots,X_n\)为样本,\(R_n(a,b)\)表示落在区间\((a,b]\)中的个数。由中值定理得,存在\(x_0\in(a,b]\)使得

\[f(x_0)=\frac 1{b-a}\int_a^b f(x)dx\approx \frac {R_n(a,b)}{n(b-a)}\]

\(-\infty<t_0<t_1<\dots<t_m<\infty\)\(t_{i+1}-t_i=h>0\). 直方图法的密度估计为:

\[ f_n(x)= \begin{cases} \frac{R_n(t_i,t_{i+1})}{nh},\ x\in(t_i,t_{i+1}],i=0,\dots,m-1\\ 0, x\le t_0,x>t_m \end{cases} \]

实际上选取\(t_0\)为比\(X_{(1)}\)略小的数,选取\(t_m\)为比\(X_{(n)}\)略大的数。经验法则\(m\approx 1+3.322\log_{10} n.\)

案例:身高数据

直方图法的相合性

定理 2.8 \(f(\cdot)\)在点\(x\)连续且\(\lim_n h_n=0,\lim_n nh_n=\infty\), 则对任何\(\epsilon>0\)

\[\lim_{n\to\infty} P(|f_n(x)-f(x)|\ge \epsilon)=0.\]

定理 2.9 \(f(\cdot)\)\(\mathbb{R}\)上一致连续,\(\int_{-\infty}^\infty |x|^\delta d x<\infty\)(对某个\(\delta>0\)), 且\(\lim_n h_n=0,h_n\ge (\ln n)^2/n\), 则

\[P(\lim_{n\to\infty} \sup_x|f_n(x)-f(x)|=0)=1.\]

证明陈家鼎等编著的教材pp54-55.

2.4.3 核估计法

中心差分

\[f(x)\approx \frac{F(x+h)-F(x-h)}{2h}\approx \frac{F_n(x+h)-F_n(x-h)}{2h}\]

\[\hat{f}_n(x) = \frac{1}{2hn}\sum_{i=1}^n 1\{x-h<X_i\le x+h\}=\frac{1}{2hn}\sum_{i=1}^n K_0\left(\frac{x-X_i}{h}\right)\]

其中

\[K_0(x)= \frac 12 1\{-1\le x<1\}\]

核函数\(K(x)\)\(\mathbb{R}\)上的非负函数且满足\(\int_{-\infty}^\infty K(x)=1\).

核估计\(\hat{f}_n(x) = \frac{1}{2hn}\sum_{i=1}^n K\left(\frac{x-X_i}{h}\right)\)

常用的核函数

  1. 均匀核函数:

\[K_0(x)= \frac 12 1\{-1\le x\le1\}\]

\[K_1(x)= 1\{-1/2\le x\le1/2\}\]

  1. 正态核函数:

\[K_2(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}\]

核估计的相合性

定理 2.10 设核函数\(K(x)\)满足条件

\[\int_{-\infty}^\infty (K(x))^2 dx<\infty,\ \lim_{|x|\to \infty} |x|K(x)=0,\]

又密度函数\(f\)在点\(x\)处连续,且\(h_n\to 0\), \(nh_n\to\infty\), 则对一切\(\epsilon>0\), 有

\[\lim_{n\to\infty} P(|\hat{f}_n(x)-f(x)|\ge \epsilon) = 0.\]

证明见pp56-58.

例 2.25 R软件包dslabs收集了男生和女生的身高数据(单位英寸),由此估计男女身高总体的密度函数。

一些说明

  1. 收敛速度的比较:在满足一些正则性的条件(如,\(h_n\to 0\), \(nh_n\to\infty\))下,可以证明
  • 直方图法的均方误差为\(O(n^{-2/3})\)
  • 核估计的均方误差为\(O(n^{-4/5})\)
  1. 核估计的带宽(bandwidth) \(h_n\)如何选取?
  • 如果选择正态核函数,经验法则\(h_n\approx 1.06S_nn^{-1/5}\)
  1. 延伸阅读

2.5 本章习题

习题 2.1 \(X\)的分布密度函数为

\[f(x)=\frac{1}{2\sigma} e^{-|x|/\sigma}\ (\sigma>0),\]

\(X_1,\dots,X_n\)\(X\)的样本,求\(\sigma\)的最大似然估计。

习题 2.2 \(X_1,\dots,X_n\)是来自\([\theta,\theta+1]\)上均匀分布的样本,其中\(\theta\in\mathbb{R}\), 证明\(\theta\)的最大似然估计不止一个,并求出所有的最大似然估计。

习题 2.3 设随机变量\(X\)以均等机会按\(N(0,1)\)分布取值和按\(N(\mu,\sigma^2)\)分布取值,其中\(\mu\in \mathbb{R},\sigma^2>0\). 这时\(X\)的分布密度函数为这两个分布的密度的平均,即

\[f(x;\mu,\sigma^2) = \frac 12\frac{1}{\sqrt{2\pi}}e^{-x^2/2}+\frac 12\frac{1}{\sqrt{2\pi}\sigma}e^{-(x-\mu)^2/(2\sigma^2)},\]

\(X_1,\dots,X_n\)为此混合分布的简单随机样本,证明\(\mu,\sigma^2\)不存在最大似然估计。能否通过矩法估计\(\mu,\sigma^2\)

习题 2.4 (附加题I,选做)考虑上题的模型。设\(Y\)为一随机变量,\(Y=1\)表示\(X\)来自\(N(0,1)\)分布,\(Y=0\)表示\(X\)来自\(N(\mu,\sigma^2)\)分布,即\(Y\sim b(1,0.5)\). 假设我们可以观测\(Y_i\)的值,基于样本\((X_i,Y_i),i=1,\dots,n\),是否可以求出\(\mu,\sigma^2\)的最大似然估计?事实上,\(Y_i\)的值不可观测(通常称为潜变量),此时你有没有更好的办法估计\(\mu,\sigma^2\)

习题 2.5 (附加题II,选做)若考虑更一般的混合分布:

\[f(x;\lambda,\mu_1,\sigma_1^2,\mu_2,\sigma_2^2)=\frac{\lambda}{\sqrt{2\pi}\sigma_1}e^{-(x-\mu_1)^2/(2\sigma_1^2)}+\frac{1-\lambda}{\sqrt{2\pi}\sigma_2}e^{-(x-\mu_2)^2/(2\sigma_2^2)}\]

其中\(\lambda\in[0,1],\mu_1,\mu_2\in \mathbb{R},\sigma_1^2,\sigma_2^2>0\), 你能求出未知参数\(\lambda,\mu_1,\sigma_1^2,\mu_2,\sigma_2^2\)的矩估计吗?

习题 2.6 \(X_1,\dots,X_n\)是来自分布密度为

\[f(x;\theta)=\frac{\Gamma(\theta+1)}{\Gamma(\theta)\Gamma(1)}x^{\theta-1}1\{0\le x\le 1\}\]

的总体的样本,其中\(\theta>0\), 试用矩法估计\(\theta\).

习题 2.7 \(X_1,\dots,X_n\)是来自分布密度为

\[f(x;c,\theta)=\frac{1}{2\theta}1\{c-\theta\le x\le c+\theta\}\]

的总体的样本,其中\(\theta>0,c\in\mathbb{R}\), 试用矩法估计\(c,\theta\).

习题 2.8 \(X_1,\dots,X_n\)为来自参数为\(\lambda\)的Poisson分布的样本. 在下列选项中选出用于估计参数\(\lambda\)的无偏估计量。( )

A. \(\bar X\)

B. \(S_n^{*2}=\frac{1}{n-1}\sum_{i=1}^{n}(X_i-\bar X)^2\)

C. \(\frac 1 {n-1}\sum_{i=1}^{n-1}X_i\)

D. \(S_n^2=\frac{1}{n}\sum_{i=1}^n(X_i-\bar X)^2\)

E. \(\frac{1}2 \bar X + \frac 12 S_n^{*2}\)

习题 2.9 \(X_1,\dots,X_n\)为来自参数为\(\lambda\)的Poisson分布的样本, 已知\(\bar X\)是未知参数\(\lambda\)的完全统计量。在下列选项中选出用于估计参数\(\lambda\)的最有效的估计量。 ( )

A. \(\bar X\)

B. \(S_n^{*2}=\frac{1}{n-1}\sum_{i=1}^n(X_i-\bar X)^2\)

C. \(\frac 1 {n-1}\sum_{i=1}^{n-1}X_i\)

D. \(\frac{1}2 \bar X + \frac 12 S_n^{*2}\)

习题 2.10 \(X,\dots,X_n\)为来自参数为\(\lambda\)的Poisson分布的样本,求\(\lambda^2\)的无偏估计。已知\(\bar X\)是参数\(\lambda\)的完全统计量,能否找到\(\lambda^2\)的最小方差无偏估计量?

习题 2.11 \(X_1,\dots,X_n\)\(N(\mu,\sigma^2)\)分布的样本,参数\(\mu,\sigma^2\)未知。证明样本方差\(S_n^2\)与修正样本方差\(S_n^{*2}\)均为\(\sigma^2\)的弱相合估计量。

习题 2.12 \(X_1,\dots,X_n\)为总体\(N(\mu,\sigma^2)\), 其中\(\mu\)已知,\(\sigma^2\)未知。证明\(\sigma^2\)的估计量

\[T(X_1,\dots,X_n)=\frac 1n\sum_{i=1}^n(X_i-\mu)^2\] 的方差达到C-R不等式的下界。

习题 2.13 Let \(X_1,\dots,X_n\) be a simple random sample taken from the density

\[f(x;\theta)=\frac{2x}{\theta^2},\quad 0\le x\le \theta.\]

  1. Find an expression for \(\hat\theta_L\), the maximum likelihood estimator (MLE) for \(\theta\).

  2. Find an expression for \(\hat\theta_M\), the method of moments estimator for \(\theta\).

  3. For the two estimators \(\hat\theta_L\) and \(\hat\theta_M\), which one is more efficient in terms of mean squared error (MSE)?

习题 2.14 \(X_1,\dots,X_n\)\(U(0,\theta)\)的样本,求\(\theta\)的置信水平为\(1-\alpha\)的置信区间。设得到了\(5\)个样本值\(0.08,0.28,0.53,0.91,0.89\), 求\(\theta\)的置信水平为\(0.95\)的置信区间。

习题 2.15 陈家鼎等编著教材P62第27题
习题 2.16 陈家鼎等编著教材P62第28题

习题 2.17 分析R软件的dslabs包中的身高数据heights, 利用R软件完成以下问题。

  1. 假设整个总体服从正态分布,求期望和方差的95%置信区间。

  2. 为了判断“正态总体”的假设的合理性,画图比较核估计密度与正态分布密度的差异?

  3. 假设男生总体与女生总体均服从正态分布(方差相同)且独立,求这两个总体平均水平的差的95%置信区间。可否认为男生总体的平均身高大于女生总体的平均身高?你的理由是什么?

  4. 为了考察第3问中“男女总体的方差相同”的假设是否合理,不妨求这两个总体的方差比的95%置信区间。并观察该置信区间是否包含1?

相关的R语言操作见 【R画图技巧】 【R密度估计】

出现无法安装R的包“dslabs”情况:为了方便大家使用数据我把里面的数据集导出来。这样就不需要安装那个包,直接导入数据data.RData即可,具体的使用方法如下