【数理统计基础】 03 - 参数估计

　　数理统计要解决的问题是，根据样本的信息猜测随机变量的信息。随机变量的分布可能完全未知，也可能已经判定为某类分布\(f(x,\theta_1,\cdots,\theta_k)\)，但有未知参数\(\bar{\theta}=(\theta_1,\cdots,\theta_k)\)，这是数理统计中最常研究的情景。

1. 点估计

　　一类最简单的问题是，要求给出参数函数\(g(\bar{\theta})\)的一个统计量估计\(\hat{g}(X_1,\cdots,X_n)\)。因为对于某次试验，估计量\(\hat{g}\)是确定的值，它被称为\(g\)的点估计。需要估计的函数\(g\)往往是参数\(\bar{\theta}\)本身，这里只讨论常用分布的参数估计。

1.1 矩估计

　　点估计就是值的近似，而分布的矩和样本矩正是一对现成的近似关系，而且大数定理也保证了它们的极限关系。被估函数\(g\)和估计量\(\hat{g}\)分别取矩和样本矩的方法就叫矩估计法，只需要有\(k\)个不同的联立方程（1）便能得到参数的估算值\(\hat{\theta}_1,\cdots,\hat{\theta}_k\)（但不一定所有时候都要求估计参数）。而常见的分布只有一个或两个参数（正态分布），往往用样本均值\(\bar{X}\)或样本方差\(S^2\)就能得到很好的估计。

\[\hat{g}_i(X_1,\cdots,X_n)=g_i(\theta_1,\cdots,\theta_k),\;\;(i=1,\cdots,k)\tag{1}\]

　　对于常用分布的矩估计，这里就不一一列举了，计算并无本质困难。这里仅给出均匀分布\([\theta_1,\theta_2]\)的矩估计（均值和方差）结果（式（2）），请注意和其它方法的比较。

\[\hat{\theta}_1=\bar{X}-\sqrt{3}S;\;\;\hat{\theta}_2=\bar{X}+\sqrt{3}S\tag{2}\]

1.2 最大似然估计

　　矩估计法是一个非常符合直观的方法，但并不是唯一的方法，换一个思维，也能找到别的“合乎情理”的方法。站在概率论的角度，我们希望能找到参数“可能性最大”的值，但关于可能性的度量还比较含糊。教材上直接告诉我们，如果分布的密度函数是\(f(x,\bar{\theta})\)，则整个样本的密度函数为式（3），对于确定的采样\((x_1,\cdots,x_n)\)，合适的\(\hat{\theta}_i\)应当使得\(L(\bar{\theta})\)达到最大值。用这样的\(\hat{\theta}_i\)作为\(\bar{\theta}\)的点估计的方法叫做最大似然估计法，\(L(\bar{\theta})\)则称为似然函数。

\[L(x_1,\cdots,x_n,\theta_1,\cdots,\theta_k)=\prod_{i=1}^n f(x_i,\theta_1,\cdots,\theta_k)\tag{3}\]

　　这个方法看起来很合理，但和矩估计法差别太大，甚至让人担心会得出相矛盾的结论。并且如果仔细推敲，这样的似然函数其实并不合理，因为\(L(\bar{\theta})\)最大和最可能的\(\bar{\theta}\)还是有很大差别的，\(L(\bar{\theta})\)并不能轻易看成\(\bar{\theta}\)的密度函数。讲到这里，其实我们已经触及到数理统计中的一个争论点了，就是\(\bar{\theta}\)究竟应该看成确定值还是随机变量？

　　矩估计法就是把\(\bar{\theta}\)看做确定值，这也是符合直觉的。但如果非要讨论最大可能性的\(\bar{\theta}\)，就不得不把它当做随机变量看待，这就是贝叶斯思想。似然函数本质上应当是个条件概率\(P(A|B)\)，条件\(B\)就是观察值\(x_1,\cdots,x_n\)，但初始概率\(P(A)\)是什么？这就是问题的关键，\(\bar{\theta}\)应当有个初始分布，答案很简单，初始的\(\bar{\theta}\)默认是均匀分布的。这正是最大似然法适用的场合与实际意义，使用时一定要确保这个假设是成立的，详细的贝叶斯法在下一段展开讨论。

　　为了计算方便，一般是求解\(\ln L(\bar{\theta})\)的极值，即求联立方程（4）的解，得到的是函数不动点，有时还要论证是否为最值。常用分布的最大似然估计大多与矩估计法的结果一样，这是一种巧合，但也说明了最大似然估计法的有效性。不过也有与矩估计法不同的，比如正态分布的方差，得到的估计是\(m_2\)，而非修正的样本方差。再比如均匀分布\([\theta_1,\theta_2]\)，由于密度函数非零的部分是\((\theta_2-\theta_1)^{-n}\)，显然在\(\theta_2-\theta_1\)最小时取得最大值，故有式（5）的估计。

\[\frac{\partial[\ln L(\theta_1,\cdots,\theta_k)]}{\partial\theta_i}=0,\;\;(i=1,\cdots,k)\tag{4}\]

\[\hat{\theta}_1=\min\{X_i\};\;\;\hat{\theta}_2=\max\{X_i\}\tag{5}\]

1.3 优良性准则

　　对于同一个参数，可以有不同的点估计，在具体的场景下应当如何选择，是很重要的问题。在制定准则时，有两点需要注意：一是判定准则是根据具体需求制定，好坏并不是绝对的；二是判定准则往往是针对全体样本的，某个具体样本的好坏不足以说明问题。

　　最简单的准则就是无偏性，它要求\(E(\hat{g}(\theta))=g(\theta)\)。无偏性适用于多次误差可以补偿的情况，比如买东西的重量，误差造成的双方损失可以互补。但对于精度要求高的场景，还希望\(\hat{g}(\theta)\)尽量聚拢在\(g(\theta)\)周围，也就是说它的方差还要尽量小。在所有无偏估计中，方差最小的称为最小方差无偏估计，简称MVU估计。

　　MVU估计比较难找，甚至根本不存在。有一个朴素的思想是，如果能得到\(D(\hat{g}(\theta))\)的一个下界，并且正好找到了这样的\(\hat{g}(\theta)\)，那么就是找到了MVU估计。这个看似异想天开的方法，居然还真有比较好的结论，下面来看看。结论的灵感来自于不等式\(\text{Cov}^2(\xi,\eta)\leqslant D(\xi)D(\eta)\)，其中等式成立的充要条件是\(\xi,\eta\)（中心化后）有简单的线性关系。

　　构造的思路是这样的，需要选择一个统计量\(G\)，它使得\(\text{Cov}(\hat{g},G)\)是与\(\hat{g}\)无关的常量，然后所有\(\hat{g}\)中还存在与\(G\)有简单线性关系的统计量。对\(\hat{g}\)的唯一限制条件是其期望为\(g\)，它也是唯一可利用的等式，因此\(G\)要取与样本分布密度\(p=\prod f(X_i,\theta)\)有关的函数。这个问题正面求解似乎很难，我们不妨从最简单的场景入手，就拿正态分布的均值估计\(\bar{X}\)为例，与\(p\)有关且与\(\bar{X}\)有线性关系的量是式（6）中的\(G\)。

\[G(X_1,\cdots,X_n)=\sum_{i=1}^n[\ln f(X_i,\theta)]'_{\theta}=\sum_{i=1}^n\dfrac{f'_{\theta}(X_i,\theta)}{f(X_i,\theta)}\tag{6}\]

　　可以计算得到，\(\text{Cov}(\hat{g},G)=g'(\theta)\)，然后还能得到\(E(G)=0\)，接下来由\(X_i\)的独立性可以把\(D(G)\)记作\(nI(\theta)\)。最终得到式（7）的克拉美-劳不等式，其中\(I(\theta)\)被称为费歇尔信息量（如果\(X_i\)是离散的也有类似的表达式），也许在信息论中会有更好的阐述，这里就不探究了。证明中还需满足一些一致性的要求，这里也省略了，请自行参考教材。

\[D[\hat{g}(X_1,\cdots,X_n)]\geqslant\dfrac{[g'(\theta)]^2}{nI(\theta)};\;\;I(\theta)=\int_{-\infty}^{\infty}\dfrac{[f'_{\theta}(x,\theta)]^2}{f(x,\theta)}\,\text{d}x\tag{7}\]

　　式（6）对任何统计量都成立，我们更应当关注等号成立的条件，即\(\hat{g}\)和\(G\)有线性关系。对于均值估计\(\bar{X}\)，要想它是MVU估计，只需\([\ln f(x,\theta)]'_\theta\)是\(x\)的线性函数（当然还要验证一致性，这里略去），容易验证常见分布一般都满足这个条件。另外还可以证明，方差\(S^2\)也是\(\sigma^2\)的MVU估计。

2. 区间估计

　　参数估计的目的是对参数更多的了解，点估计的结果虽然直接易用，但却丢失了太多参数的信息，使用上也没有灵活性。为了包含参数的更多信息，我们希望找到两个统计量\(\hat{g}_1,\hat{g}_2\)，以区间形式估计参数，并达到一定的概率要求。一般是对给定足够小的\(\alpha>0\)，要找到尽量小的区间\([\hat{g}_1,\hat{g}_2]\)，使它能以\(1-\alpha\)的概率包含\(g(\theta)\)（式（8））。

\[P[\hat{g}_1(X_1,\cdots,X_n)\leqslant g(\theta)\leqslant\hat{g}_2(X_1,\cdots,X_n)]=1-\alpha\tag{8}\]

　　这样的估计方法叫区间估计，其中\([\hat{g}_1,\hat{g}_2]\)叫置信区间，而\(1-\alpha\)是区间的置信系数。有两点需要强调：一个是这里仍然是把参数当做确定值，把样本当做随机变量，所以置信系数的意义是“区间能包含参数”的概率，而非“参数落在区间里”的概率；另一个是区间长度越小越好，但不做强求，因为区间本身就是随机变量，对它最小值的讨论比较困难。

　　为了构造统计量\(\hat{g}_1,\hat{g}_2\)，观察式（8），其中只包含待估参数和样本值，以及它们之间的概率不等式。一种比较方便的构造方法是这样的，找一个变量\(G(g(\theta),X_1,\cdots,X_n)\)，它服从一个比较简单的分布\(F\)。为了生成置信区间，一般把变量值限定在\(E(G)\)的两侧，每测的概率分别取\((1-\alpha)/2\)。如果用\(f(\alpha)\)表示\(F\)的\(\alpha\)分位点，则建立不等式（9），整理后便能得到式（8）置信区间。

\[E(G)-f(1-\frac{\alpha}{2})\leqslant G(g(\theta),X_1,\cdots,X_n)\leqslant E(G)+f(\frac{\alpha}{2})\tag{9}\]

　　这里的\(G\)就是前面提到过的枢轴变量，因此该方法也叫枢轴变量法。很多分布的枢轴变量比较难构造或者计算量大，甚至有时对分布完全未知，这时如果样本足够大，可以利用中心极限定理，以标准正态分布作为枢轴变量。这里我们只讨论正态分布，它的常用枢轴变量还有上篇介绍的三大变量，请先回顾相关性质。以下讨论仅给出枢轴变量，具体置信区间请自行计算，并没有本质困难。

　　先讨论单样本的正态分布\(X\sim N(\mu,\sigma^2)\)。估计均值\(\mu\)时，\(\sigma\)可能已知也可能未知，上篇的公式（8）和（15）便是对这两种情况的枢轴变量。再估计方差\(\sigma^2\)，同样分为\(\mu\)已知和未知两种情况，\(\mu\)已知的情况比较简单，未知时上篇的公式（17）便是我们要的枢轴变量。

　　再讨论两样本的正态分布，一般是根据两个随机变量的观察值，比较它们的参数。设两个随机变量为\(X\sim N(\mu_1,\sigma_1^2)\)和\(X\sim N(\mu_2,\sigma_2^2)\)，样本分别是\(X_1,\cdots,X_m\)和\(Y_1,\cdots,Y_n\)。一种是要考察\(\mu_1-\mu_2\)的大小，一般的做法当然是用\(\bar{X}-\bar{Y}\)去估计它。当\(\sigma_i\)都已知时，\(\bar{X}-\bar{Y}\)的方差为\(\sigma^2=\dfrac{\sigma_1^2}{m}+\dfrac{\sigma_2^2}{n}\)，枢轴变量比较显然。

　　当\(\sigma_i\)都未知时，暂时没办法把\(\sigma\)消除（即使用\(S_1^2+S_2^2\)也不行），这里只讨论\(\sigma_1=\sigma_2\)的场景。为了能使用\(t\)分布，直接使用式（10）中的\(S^2\)来近似方差，容易得到枢轴变量（11）。当\(\sigma_1\ne \sigma_2\)时，暂时没有完美的解决方法，该问题称为贝伦斯-费歇尔问题。

\[S^2=\dfrac{\sum\limits_{i=1}^m(X_i-\bar{X})^2+\sum\limits_{i=1}^n(Y_i-\bar{Y})^2}{m+n-2}\tag{10}\]

\[\dfrac{(\bar{X}-\bar{Y})-(\mu_1-\mu_2)}{\sqrt{1/m+1/n}\cdot S}\sim t_{m+n-2}\tag{11}\]

　　最后来比较\(X,Y\)的方差，直接作差比较难处理，而且意义也不明显。一般是估计\(\sigma_1^2/\sigma_2^2\)的大小，它可以直接使用上一篇的式（19）作为枢轴变量。

3. 贝叶斯估计

　　现在来正式讨论贝叶斯估计，它的模型直接从事件的条件概率扩展而来，只不过由事件概率扩展为分布密度（同样适用于离散分布）。贝叶斯法的最大特点就是把参数\(\theta\)看做一个随机变量，如何理解这一点非常关键。现实中参数\(\theta\)一定是确定的，只不过我们不知道它的信息。但根据过去的认识或合理的假设，对\(\theta\)的所有可能值会有个评估，这样的评估就使得\(\theta\)有了随机变量的性质。需要着重强调的是，随机变量不是变量，它只是对不同可能性的一种描述。

　　脑洞再开大一点，则可以认为我们以一定概率处在不同的平行时空中，而\(\theta\)在每一个时空中都有一个确定的值。在得到观察\(X_1,\cdots,X_n\)后，我们需要重新评估处在不同时空的概率。这是典型的条件概率问题，但要注意，这时讨论的样本空间是\(\theta,x_1,\cdots,x_n\)。假设随机变量\(X\)的密度函数为\(f(x,\theta)\)，参数\(\theta\)的先验分布的密度函数为\(h(\theta)\)，容易得到\(\theta\)的后验概率的密度函数（式（12））。

\[h(\theta\,|\,x_1,\cdots,x_n)=\dfrac{p(\theta)}{\int p(\theta)\,\text{d}\theta},\;\;p(\theta)=h(\theta)\prod_{i=1}^nf(x_i,\theta)\tag{12}\]

　　前面说过，最大似然估计本质上也是贝叶斯思想，只是先验分布采用的是均匀分布。这里有个很现实的问题，如何在无限区间（比如整个实数域、所有正数等）上定义均匀分布？这个对我们还是太困难，也许测度论中会有完美解释？这就不得而知了。但还是可以换个思路，\(h(\theta)\)在条件分布中本质上起到的是“权重”的作用，也就是说它的根本意义在于表明概率之间的“比重”。比如对于均匀分布，只需取\(h(\theta)=1\)就足以说明“均匀”的性质，不必要求\(h(\theta)\)是一个严格的密度函数（但由式（12）易知后验概率一定是密度函数）。

　　但无论如何，在式（12）的使用过程中，必须先给出一个先验分布\(h(\theta)\)，这个分布的选择非常影响估计结果。很多时候先验分布难以确定，只能凭主观经验，这给贝叶斯方法带来了很多诟病。但由于贝叶斯思想的有效性和方便性，它在数理统计中仍然大行其道，甚至形成了所谓的贝叶斯学派，以区别于坚持频率方法的学者。一种和解的方法是承认两个模型本质的不同，并且互相补充、互相学习。但以个人粗浅的了解，我觉得贝叶斯思想是对传统模型的扩充，它是用先验概率把传统模型补充完整而已。这个补充就如同虚数于实数系统一样，是打破直觉却非常必要的抽象，是现代数学所具有的特征。

　　贝叶斯模型是完整的，且逻辑自洽的，方法本身不应该被诟病。既然问题出在先验概率的选择，那么在使用时挑选最合适的即可。这就是另外一个问题了，需要更多的理论分析和支持，不能把这部分工作的欠缺怪罪于贝叶斯模型本身。在大部分场合，\(h(\theta)\)一般遵循“同等无知”原则，这个原则的缺点也是显然的：如果对\(\theta\)是同等无知的，则对它的函数则基本不满足。这时选择对谁同等无知就很关键，在正态分布中一般取\(h(\sigma)=\sigma^{-1}\)，在指数分布中一般取\(h(\lambda)=\lambda^{-1}\)。

　　后验概率可以看做包含了参数的所有信息，它可以被用作点估计和区间估计。在点估计时，最合理的应当是取期望值，而非最大似然法中的最大值。伯努利分布（先验概率取均匀分布）在贝叶斯方法中的期望值是\(\hat{p}=\dfrac{N+1}{n+2}\)，这在\(n\)很小时明显更合理（最大似然法得到\(\dfrac{N}{n}\)）。

　　后验概率上的区间估计实现起来非常方便，它有统一的计算过程，而不依赖于具体分布。针对指定的置信系数，寻找最小的置信区间，是存粹的分析学计算问题。后验分布有对称轴时，最小置信区间一般也是对称的。其它情况下，可以先固定一个边界以确定另一个边界，然后变动第一个边界寻找最小区间。实在复杂的计算，也可以直接交给计算机完成。最后需要提醒一下，这里的置信区间和第二段中的置信区间有着本质的区别，一个是\(\theta\)自身（随机变量）的取值区间，一个是可以包含\(\theta\)（确定值）的区间，请仔细体会。