一文看懂 “极大似然估计” 与 “最大后验估计” —— 极大似然估计篇

• 参考：
  1. 唐宇迪《人工智能数学基础》第8章
  2. Richard O. Duda 《模式分类》第三章
  3. 白板机器学习 P2 - 频率派 vs 贝叶斯派
  4. 频率学派还是贝叶斯学派？聊一聊机器学习中的MLE和MAP
• 本文符号说明：
  1. Θ = { Θ 1 , Θ 2 , . . . , Θ k } \Theta=\{\Theta_1,\Theta_2,...,\Theta_k\} Θ={Θ1​,Θ2​,...,Θk​}：待估参数随机向量，上加 $\hat{}$ 表示估计值随机向量（${\Theta }_i$ 可表示向量或单一变量）
  2. θ = { θ 1 , θ 2 , . . . , θ k } \theta=\{\theta_1,\theta_2,...,\theta_k\} θ={θ1​,θ2​,...,θk​}：待估参数向量的真实值，上加 $\hat{}$ 表示估计值（${\theta }_i$ 可表示向量值或单一变量值）
  3. $X=x$：一次随机试验对应的随机变量 $X$ 取值到 $x$，理解为某次样本采样得到 $x$
  4. $w_i$：代表某个样本标签为类别 $w_i$
  5. $\mathcal{D}$：训练数据集
• 本文历次修订后全长 2万8000余字，受到 CSDN 博文字数限制，故切分两篇发布，所以现在是两文看懂了… 本篇介绍参数估计背景和极大似然估计；下一篇介绍最大后验估计和两种方法对比
• 后文链接：一文看懂 “极大似然估计” 与 “最大后验估计” —— 最大后验估计篇

---

1. 引入

1.1 参数估计

• 参数估计：根据从总体中抽取的随机样本来估计总体分布中未知参数的过程。参数估计在基于统计学的机器学习中特别重要，因为学习的目标就是从假设空间中选出一个最优模型，而模型本质上就是一个从样本空间到标记空间的空间的映射，可以把它看作一个由 $\Theta$ 参数化的条件概率分布。某种程度上说，机器学习的过程，就是参数估计的过程，具体而言分以下两种情况

  1. 无标签情况（这时模型表示样本取值的分布，做参数估计就好比重复投一个骰子来估计其各个面出现的概率）：关注以下问题

     根据从分布 $p(x|\Theta )$ 中 i.i.d 采样的数据集 $\mathcal{D}$ 估计参数 $\Theta$ 的参数值

  2. 有标签情况（这时模型表示某样本值下标签的条件分布，一般的机器学习任务都是这个形式）：这里我们按照各个样本所属类别标签 $w_i$ 将数据集 $\mathcal{D}$ 进行划分为 { D i } \{\mathcal{D}_i\} {Di​}，并且假设每一个类别下的类条件概率 $p(x|w_i)$ 由一个参数向量 ${\Theta }_i$ 参数化，且 ${\Theta }_i$ 对和其他类别 $j\ne i$ 的样本不起任何作用（例如每个类对应一个独立的多维高斯分布），假设标记空间中一共有 $c$ 个类，则问题转化为 $c$ 个相互独立的问题

     根据从分布 $p(x|w_i,{\Theta }_i)$ 中 i.i.d 采样的数据集 ${\mathcal{D}}_i$ 估计参数 ${\Theta }_i$ 的参数值

     注意由于各个问题都是独立针对某一个类的，所以每一个问题内部可以忽略类别的存在，还原到 1 的情形从而简化符号，这时可以描述为以下 $c$ 个相互独立的问题

     根据从分布 $p(x|\Theta )$ 中 i.i.d 采样的数据集 $\mathcal{D}$ 估计参数 $\Theta$ 的参数值

  这里还有一点需要注意，我们的最终目标是不是通过 $\mathcal{D}$ 估计出 $\hat{\Theta }$，而是利用估计的 $\hat{\Theta }$ 还原样本分布 $p(X|\hat{\Theta })$ 得到模型，参数估计只是求解最优模型的方法。如果我们把参数估计看作内部过程隐去，模型可以表示为 $p(X|\mathcal{D})$

  1. 对于无标签情况，模型 $p(X|\mathcal{D})$ 给出某次掷骰子结果 $x$ 发生的概率
  2. 对于有标签情况，模型 $p(X|\mathcal{D})$（这里是 $p(x|w_i,{\mathcal{D}}_i)$ 的简写） 给出属于类别 $w_i$ 的情况下出现某个样本 $x$ 的概率，对每个类都算一下再取最大值就能实现分类任务

• 模型参数化：为了估计模型参数，我们必须了解模型抽取样本的概率形式，即 $p(X|\Theta )$ 的数学形式（否则就是一个黑箱，无从分析）。模型参数化是依靠先验信息简化模型形式，减少待估参数的一种重要方法

  1. 在最简单的情况下，模型是一个表格形式，对于输入空间中的任意样本 $x$ 都有一个对应的参数 $\theta$ 代表采样到它的概率。假设样本有 p 个维度，每个维度 q 个离散取值，则每一个类别下估计 $p(x|\Theta )$ 的参数数量为 $p^q$；当任意维度是连续取值时，参数数量变为无穷，显然这样做是不合理的
  2. 如果我们事先知道参数的个数，并且先验知识允许我们把条件概率密度进行参数化，那么问题的难度就可以显著降低。例如我们可以假设正确的 $p(x|w_i,{\Theta }_i)$ 是一个多元正态分布，这样我们只需要估计其均值 ${\mu }_i$ 和协方差矩阵 ${\sum }_i$ 即可，这样就大大减少了待估参数量。机器学习中常见的神经网络就是一种典型的参数化模型

• 本文讨论参数估计的两种经典方法：极大似然估计（MLE）和最大后验估计（MAP），他们分别体现了概率论中频率派和贝叶斯派的思想

1.2 频率派 vs 贝叶斯派

• 抽象一点说，频率学派和贝叶斯学派代表了对世界的两种认知：
  1. 频率学派认为世界是确定的，他们认为待估参数是参数空间中某个确定的取值，只是其取值未知，数据都是在这个确定的参数值下产生的，最佳估计就是使得训练样本出现概率最大的值，所以他们的方法论是从 “哪个值最有可能是真实值” 这个角度出发的，于是就有了最大似然（maximum likelihood）以及置信区间（confidence interval）这样的东西。注意一下，频率学派的人们认为参数本身是一个定值，不是随机变量，但置信区间是随着样本数量而改变的，是一个随机变量。频率派关心的就是有多大把握去圈出那个唯一的真实参数
  2. 贝叶斯学派认为世界是不确定的（或不能确定的），他们认为参数空间里的每个值都有可能是真实模型使用的值，区别只是概率不同而已。于是他们引入先验分布（prior distribution）和后验分布（posterior distribution）这样的概念来设法找出参数空间上的每个值的概率。对样本进行观测的过程，就是把先验概率密度转化成后验概率密度的过程，每得到一个样本，都使得后验概率密度函数变得更加尖锐，使其在待估参数的真实值附近形成最大的尖峰
• 具体应用上
  1. 频率学派的目标是求解使得数据产生的那个唯一且确定的模型参数 $\Theta$，因此这一派产生的方法往往转换为优化问题，其结构一般是：确定模型 —— 设计损失函数 —— 优化损失。这个结构很眼熟吧，目前大部分基于统计的机器学习方法都是属于频率派的
  2. 贝叶斯学派中，由于模型参数 $\Theta$ 是一个分布，在使用贝叶斯派的模型时（比如做贝叶斯预测时）需要对 $\theta$ 求积分，因此这一派方法最终往往转换为求积分的问题。$\Theta$ 的积分需要在整个参数空间中求，难以得到解析解，通常使用蒙特卡洛采样的方式求其数值解。由贝叶斯学派衍生的机器学习方法有概率图模型、mcmc方法等
• 目前教科书上大部分讲的都是频率派的内容，比如参数点估计、假设检验等，这里面有一些历史原因，下面引用一段 知乎用户Xiangyu Wang的回答

  如果从概率的角度看，贝叶斯学派的想法其实更为自然，这也是为什么贝叶斯学派的产生远早于频率学派。但是贝叶斯方法本身有很多问题，比如当先验选的不好或者模型不好的时候你后验分布的具体形式可能都写不出来，跟别说做统计推断了。在当年电子计算机还没发展出来的时候，对这些情况做分析几乎是不可能的，这也就大大限制了贝叶斯方法的发展。而频率学派主要使用最优化的方法，在很多时候处理起来要方便很多。所以在频率学派产生后就快速地占领了整个统计领域。直到上世纪90年代依靠电子计算机的迅速发展，以及抽样算法的进步（Metropolis-hastings, Gibbs sampling）使得对于任何模型任何先验分布都可以有效地求出后验分布，贝叶斯学派才重新回到人们的视线当中。就现在而言，贝叶斯学派日益受到重视当然是有诸多原因的，所以这并不意味这频率学派就不好或者不对。两个学派除了在参数空间的认知上有区别以外，方法论上都是互相借鉴也可以相互转化的。当代学术领域批评的最多的仅仅是频率学派里的Hypothesis testing的问题，尤其是对于p-value的误用造成了很多问题，不过这只是Hypothesis testing这种研究方法本身的问题。对应于Hypothesis testing，贝叶斯学派有自己的一套方法称为 Bayes factor。虽然Bayes factor本身比p-value要合理很多（个人见解），但是我并不觉得单靠Bayes factor的方法就可以有效解决当下p-value滥用导致的问题，因为Bayes factor同样可以导致Multiple comparisons problem。

• 总而言之，极大似然估计和最大后验估计，都是总体的（某个指标的）分布类型已知前提下的一种参数估计方法。在对事物建模时，用 Θ = { θ 1 , θ 2 , . . . } \Theta = \{\theta_1,\theta_2,...\} Θ={θ1​,θ2​,...} 表示模型的参数，这两种方法的目标就是对分布的未知参数 $\Theta$ 进行估计，以得到确定的分布。区别在于极大似然估计仅根据观测到的结果（观测值/样例）进行估计，而最大后验估计进一步考虑了参数 $\theta$ 的先验分布 $p(\theta )$。总的来说
  1. 频率派比较简单粗暴，上来就是大量重复实验、中心极限定理什么的，在已知模型且允许大量重复实验时可以得到比较精确的概率评估，但在难以大量实验或缺乏数据时可能不适用或产生严重偏差

     极大似然估计是代表频率学派思想的参数估计方法！

  2. 贝叶斯派更符合人对概率的理解，通过引入先验概率，可以在少量数据的情况下给出比较准确的估计，当然这对先验概率的准确性有要求。随着新数据的获取，贝叶斯派会不断修正先验概率，这间接地调整了 $\theta$ 在已知事件下的的条件概率（后验概率），最终在大样本情况下取得和频率派一样的结果

     最大后验估计是代表贝叶斯学派思想的参数估计方法！

2. 先验概率、后验概率、似然函数

• 假设某个随机变量 $X$ 的取值 $x$ 受到参数 $\Theta$ 取值 $\theta$ 的影响，例如抛硬币的结果受到硬币均匀性参数的影响，我们希望通过试验找出硬币均匀性参数。这里用 “有多大可能性投出正面” 表示均匀性，即 $\theta =p(X=正面)$
  1. 先验概率$p(\Theta )$：所谓 “先验”，顾名思义就是指试验之前，不依靠观测数据的未知参数概率分布。在没有做任何试验时，人们可以利用经验或者统计结果给出先验概率，比如我们可能会猜测硬币是均匀理想的，所以先验概率是 $p(\Theta =0.5)=1$；也可能我们认为硬币有一定概率不均匀，先验概率是 $p(\Theta =0.5)=0.9,p(\Theta =0.6)=0.1$。我们通常说的概率默认指先验概率
  2. 似然函数$L(\Theta |X)$：代表试验或事件 $X=x$ 已经发生之后，在 $\Theta$ 的各个取值下得到这个结果的概率。这是一个关于参数 $\Theta$ 的函数，其在数值上等于关于 $X$ 的条件概率分布 $p(X|\Theta )$。假设我们抛硬币10次，有7次正面朝上，则这个试验结果的似然函数为 $L(\Theta |X=x)=p(X|\Theta )={\Theta }^7(1-\Theta {)}^3$，图像如下
     
     因此我们有 $\underset{\theta }{\mathrm{argmax}}L(\Theta =\theta |X=x)=0.7$，这个其实就是极大似然估计
  3. 后验概率$p(\Theta |X)$：所谓 “后验”，顾名思义就是指试验或事件 $X=x$ 已经发生之后，在得到这个结果的前提下的未知参数的条件概率分布。我们可以使用贝叶斯公式在上述三者间建立联系
     $$p(\Theta |X)=\frac{p(\Theta ,X)}{p(X)}=\frac{p(\Theta )p(X|\Theta )}{p(X)}$$ 所谓最大后验估计，就是要找到 $\underset{\theta }{\mathrm{argmax}}p(\Theta =\theta |X=x)$

3. 极大似然估计（MLE）

• 考虑这样一种情境：一个菜鸟和一个高手运动员一起打靶，现在观察到靶纸上一个十环的弹孔，最可能是谁打出的呢？显然我们会认为是高手打出的，因为他更有可能打出十环。在观察到某一现象探求其发生原因时，选择可能性最大的（概率最大的）原因，这就是最大似然估计的思想。
• 极大似然估计的示意图如下
  

3.1 似然函数

1. 似然函数和最大似然估计值

   • 现有总体 $X$ 的一个样本 $X_1,X_2,...,X_n$，观测值为 $x_1,x_2,...,x_n$，利用这一观测对未知参数 $\Theta \in I$ 进行估计，使得选定的 $\hat{\Theta }$ 最有利于 $x_1,x_2,...,x_n$ 的出现
   • 似然函数 $L(\Theta |x)$：$X_1,X_2,...,X_n$ 取值为 $x_1,x_2,...,x_n$ 的概率
     $$L(\Theta |X)=L({\theta }_1,{\theta }_2,...,{\theta }_k|x_1,x_2,...,x_n)=\{\begin{array}{rlr}& \prod _{i=1}^{n}p(x_i|{\theta }_1,{\theta }_2,...,{\theta }_k)& X是离散型随机变量\\ & \prod _{i=1}^{n}f(x_i|{\theta }_1,{\theta }_2,...,{\theta }_k)& X是连续型随机变量\end{array}$$
   • 使得似然函数 $L(\hat{\Theta }|X)$ 最大的 $\hat{\theta }=\hat{\theta }(x_1,x_2,...,x_n)\in I$，使称为参数 $\theta$ 的 最大似然估计值，而相应的统计量 $\hat{\Theta }(X_1,X_2,...,X_n)$ 称为参数 $\theta$ 的 最大似然统计量。我们的目标就是找出最大似然估计值 $\hat{\theta }$，即
     $$\begin{array}{rl}\hat{\theta }& =\underset{\theta }{\mathrm{argmax}}L(\Theta =\theta |X=x)\\ & =\underset{\theta }{\mathrm{argmax}}\prod _{i=1}^{n}p(x_i|\theta )\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{i=1}^{n}logp(x_i|\theta )\\ & =\underset{\theta }{\mathrm{argmin}}-\sum _{i=1}^{n}logp(x_i|\theta )\end{array}$$

2. 注意：

   1. 似然函数 $L(\Theta |X)$ 代表了参数 $\Theta$ 和样本 $X$ 的关联程度，其值最大时，代表此 $\Theta =\theta$ 使得样本观测值 $X=x$ 出现的概率最大，是 $\Theta$ 的最佳估计值
   2. 我们的目标是最大化 $L(\Theta |X)$，为了简化计算，可以向上面那样两边取对数，把连乘转换为连加。这样处理后得到的函数称为 对数似然函数
   3. 上面推导的最后一行所优化的函数被称为 负对数似然 Negative Log Likelihood/NLL，被很多ML方法使用，比如逻辑斯蒂回归中的参数估计就是执行的这个优化，另外这个和交叉熵损失关系密切
   4. 似然函数 $L$ 和概率分布函数 $p$、概率密度函数 $f$ 是完全不同的两个数学对象：前者是关于 $\Theta$ 的函数，后者是关于 $X$ 的函数。所以这里的等号理解为函数值形式的相等，而不是两个函数本身是同一函数

3.2 最大似然估计的步骤

1. 写出样本的似然函数 $L(\Theta |X)$
2. 如果 $p(x_i|{\theta }_1,{\theta }_2,...,{\theta }_k)$ 或 $f(x_i|{\theta }_1,{\theta }_2,...,{\theta }_k)$ 关于 ${\theta }_i(i=1,2,...,k)$ 可微。则令
   $$\frac{\partial L(\Theta |X)}{\partial {\theta }_i}=0或\frac{\partial lnL(\Theta |X)}{\partial {\theta }_i}=0$$ 由于 $L(\theta )$ 是连乘形式，又 $lnx$ 是 $x$ 的单调增函数，因此 $L(\theta )$ 和 $lnL(\theta )$ 在同一 $\theta$ 取极值，所以通常使用解对数似然方程组 $\frac{\partial lnL(\theta )}{\partial {\theta }_i}=0$ 的方法，求得 ${\theta }_i$ 的最大似然估计量 $\hat{\Theta }=\hat{\Theta }(X_1,X_2,...,X_n)(i=1,2,...,k)$
3. 若 $p(x_i|{\theta }_1,{\theta }_2,...,{\theta }_k)$ 或 $f(x_i|{\theta }_1,{\theta }_2,...,{\theta }_k)$ 不可微，或似然方程组无解，则应由定义用其他方法求得 $\hat{\Theta }$，例如当 $L(\Theta |X)$ 为 $\Theta$ 的单调函数时，$\hat{\Theta }$ 为 $\Theta$ 取值的上限或下限

3.3 示例

3.3.1 抛硬币

• 现有一个正反面不是很均匀的硬币，抛10次结果为：反，反，正，反，反，正，反，正，反，正，共6反4正。设硬币正面朝上的概率为 $\Theta =\theta$，用最大似然方法估计 $\theta$
• 已知每次抛硬币都是0-1分布，其分布函数为 $p(x=k|\theta )={\theta }^k(1-\theta {)}^{1-k}(k=0,1)$，于是似然函数和对数似然函数为
  $$\begin{gathered} L(\theta |x)=\prod _{i=1}^{n}p(x|\theta )=\prod _{i=1}^n{\theta }^{x_i}(1-\theta {)}^{1-x_i} \\ lnL(\theta |x)=\sum _{i=1}^n[x_{i}ln\theta +(1-x_i)ln(1-\theta )] \end{gathered}$$
• 令导数为0解 $\theta$
  $$\begin{array}{rl}\frac{\partial }{\partial \theta }lnL(\theta |x)& =\sum _{i=1}^n\frac{\partial }{\partial \theta }[x_{i}ln\theta +(1-x_i)ln(1-\theta )]\\ & =\sum _{i=1}^n{x}_i\frac{\partial }{\partial \theta }ln\theta +\sum _{i=1}^n(1-x_i)\frac{\partial }{\partial \theta }ln(1-\theta )\\ & =\frac{1}{\theta }\sum _{i=1}^n{x}_i-\frac{1}{1-\theta }\sum _{i=1}^n(1-x_i)\\ & \stackrel{令}{=}0\end{array}$$
  带入观测值 $x_i(i=1,2,...,10)$，解得 $\theta =0.4$

3.3.2 均匀分布

• 设样本服从正态分布 $U(a,b)$，现有样本观测值 $x_1,x_2,...,x_n$，用最大似然估计方法求参数 $a,b$
  $$f(x)=\{\begin{array}{rlr}& \frac{1}{b-a}& ,a\le x\le b\\ & 0& ,其他\end{array}$$
• 似然函数为
  $$L(a,b|x)=\{\begin{array}{rlr}& \frac{1}{(b-a{)}^n}& ,a\le x_i\le b\\ & 0& ,其他\end{array}$$
• 此时似然函数作为 $a,b$ 的二元函数是不连续的，这时不能用偏导数来求解，必须从极大似然估计的定义出发求 $L(a,b)$ 最大值
• 显然，为使 $L(a,b)$ 尽量大，b-a应尽量小，在考虑到a和b的取值限制，极大似然估计为
  { a ^ = m i n { x 1 , x 2 , . . . , x n } b ^ = m a x { x 1 , x 2 , . . . , x n } \left\{ \begin{aligned} \hat{a} = min\{x_1,x_2,...,x_n\} \\ \hat{b} = max\{x_1,x_2,...,x_n\} \end{aligned} \right. {a^=min{x1​,x2​,...,xn​}b^=max{x1​,x2​,...,xn​}​

3.3.3 一维高斯分布

• 设样本服从正态分布 $N(\mu ,{\sigma }^2)$，现有样本观测值 $x_1,x_2,...,x_n$，用最大似然估计方法求参数 $\mu ,\sigma$
• 似然函数和对数似然函数为
  $$\begin{array}{rl}L(\mu ,{\sigma }^2|x)& =\prod _{i=1}^n\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(x_i-\mu {)}^2}{2{\sigma }^2})\\ & =(2\pi {\sigma }^2{)}^{-\frac{n}{2}}exp(-\frac{1}{2{\sigma }^2}\sum _{i=1}^N(x_i-\mu {)}^2)\\ lnL(\mu ,{\sigma }^2|x)& =-\frac{n}{2}ln(2\pi )-\frac{n}{2}ln({\sigma }^2)-\frac{1}{2{\sigma }^2}\sum _{i=1}^N(x_i-\mu {)}^2\end{array}$$
• 求偏导数
  $$\{\begin{array}{rl}\frac{\partial lnL(\mu ,{\sigma }^2|x)}{\partial \mu }& =\frac{1}{{\sigma }^2}\sum _{i=1}^n(x_i-\mu )\\ \frac{\partial lnL(\mu ,{\sigma }^2|x)}{\partial {\sigma }^2}& =-\frac{n}{2{\sigma }^2}+\frac{1}{2{\sigma }^4}\sum _{i=1}^n(x_i-\mu {)}^2\end{array}$$
• 令偏导数为0，解得
  $$\{\begin{array}{rl}\hat{\mu }& =\bar{x}=\frac{1}{n}\sum _{i=1}^n{x}_i\\ \widehat{{\sigma }^2}& =\frac{1}{n}\sum _{i=1}^n(x_i-\hat{\mu }{)}^2=\frac{1}{n}\sum _{i=1}^n(x_i-\bar{x}{)}^2\end{array}$$
• 分析

  1. 易证这里 $\hat{\mu }$ 是一个无偏估计量

  2. $\widehat{{\sigma }^2}$ 是一个有偏估计量，证明如下：首先化简 $\widehat{{\sigma }^2}$
     $$\begin{array}{rl}\widehat{{\sigma }^2}& =\frac{1}{n}\sum _{i=1}^n(x_i-\bar{x}{)}^2\\ & =\frac{1}{n}\sum _{i=1}^n(x_i^2-2x_i\bar{x}+{\bar{x}}^2)\\ & =\frac{1}{n}\sum _{i=1}^n{x}_i^2-\frac{2\bar{x}}{n}\sum _{i=1}^n{x}_i+\frac{1}{n}n{\bar{x}}^2\\ & =\frac{1}{n}\sum _{i=1}^n{x}_i^2-2{\bar{x}}^2+{\bar{x}}^2\\ & =\frac{1}{n}\sum _{i=1}^n{x}_i^2-{\bar{x}}^2\end{array}$$ 然后考察是否有偏
     $$\begin{array}{rl}\mathbb{E}\widehat{{\sigma }^2}& =\mathbb{E}[\frac{1}{n}\sum _{i=1}^n(x_i-\bar{x}{)}^2]\\ & =\mathbb{E}[\frac{1}{n}\sum _{i=1}^n{x}_i^2-{\bar{x}}^2]\\ & =\frac{1}{n}\sum _{i=1}^n\mathbb{E}[x_i^2]-\mathbb{E}[{\bar{x}}^2]\\ & =\frac{1}{n}n[(\mathbb{E}{x}_i{)}^2+D{x}_i]-[(\mathbb{E}\bar{x}{)}^2+D\bar{x}]\\ & ={\mu }^2+{\sigma }^2-{\mu }^2-D[\frac{1}{n}\sum _{i=1}^n{x}_i]\\ & ={\sigma }^2-\frac{n}{n^2}D{x}_i\\ & ={\sigma }^2-\frac{1}{n}{\sigma }^2\\ & =\frac{n-1}{n}{\sigma }^2\end{array}$$ 可见一维正态分布的方差的绝对无偏估计量为 $\frac{1}{n-1}{\sum }_{i=1}^n(x_i-\bar{x}{)}^2$；最大似然估计得到的 $\widehat{{\sigma }^2}$ 是渐进无偏估计量（当样本量 $n$ 趋于无穷时趋于无偏估计），只要数据集足够大，这样的结果是可以接受的。另外，用类似的方法还能证明 $\frac{1}{n}{\sum }_{i=1}^n(x_i-\mu {)}^2$ 是无偏的。

  3. 从概率论的角度看，极大似然估计得到的估计量是关于样本的n元函数，属于一种统计量，统计量是否一致非常重要，而是否无偏完全不重要。常见的很多（也许可以说是大部分）统计量都是有偏的。

3.3.4 多维高斯分布

• 设随机向量 $X\text{X}X$ 服从期望为 $\mu \text{μ}\mu$，协方差矩阵为 $\Sigma \text{Σ}\Sigma$ 的多维高斯分布，即 $X\text{X}X\sim N(\mu \text{μ}\mu ,\Sigma \text{Σ}\Sigma )$，则 $n$ 维 r.v. ${X\text{X}X}_{n\times 1}$ 的概率密度函数为
  $$f(x\text{x}x)=f(x_1,...,x_n)=\frac{1}{(2\pi {)}^{n/2}|\Sigma {|}^{1/2}}{e}^{-\frac{1}{2}(\mathbf{x}-\mu {)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mu )}$$ 其中 $\Sigma \text{Σ}\Sigma$ 是一个实对称矩阵，指数部分的 $(\mathbf{x}-\mu {)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mu )$ 是一个二次型，利用杨本观测值 ${x\text{x}x}_1,{x\text{x}x}_2,...,{x\text{x}x}_m$ 做极大似然估计
• 首先补充一些证明会用到的定理
  1. $\Sigma \text{Σ}\Sigma$ 是实对称矩阵，${\Sigma \text{Σ}\Sigma }^{-1}$ 也是实对称矩阵，有 ${\Sigma \text{Σ}\Sigma }^{-1}=({\Sigma \text{Σ}\Sigma }^{-1}{)}^{\top }$
  2. 二次型和矩阵迹的公式 $x^{\top }Ax=\text{tr}(Ax{x}^{\top })$
  3. $\frac{\partial \left(x^{\top }Ax\right)}{\partial }=\left(A+A^{\top }\right)x$
  4. $\frac{\partial \mathrm{tr}\left(X^{-1}A\right)}{\partial X}=-X^{-1}{A}^T{X}^{-1}$
  5. $\frac{\partial \ln (|X|)}{\partial X}={\left(X^{-1}\right)}^{\top }={\left(X^{\top }\right)}^{-1}$
• 首先计算对数似然函数
  $$\begin{array}{rl}lnL(\mu \text{μ}\mu ,\Sigma \text{Σ}\Sigma |x\text{x}x)& =ln\prod _{k=1}^{m}P({x\text{x}x}_k|\mu \text{μ}\mu ,\Sigma \text{Σ}\Sigma )\\ & =\sum _{k=1}^{m}lnP({x\text{x}x}_k|\mu \text{μ}\mu ,\Sigma \text{Σ}\Sigma )\\ & =-\frac{m}{2}\ln [(2\pi {)}^m|\Sigma \text{Σ}\Sigma |]-\frac{1}{2}\sum _{k=1}^m({x\text{x}x}_k-\mu \text{μ}\mu {)}^{\top }{\Sigma \text{Σ}\Sigma }^{-1}({x\text{x}x}_k-\mu \text{μ}\mu )\\ & =-\frac{m}{2}\ln [(2\pi {)}^m|\Sigma \text{Σ}\Sigma |]-\frac{1}{2}\text{tr}[\sum _{k=1}^m{\Sigma \text{Σ}\Sigma }^{-1}({x\text{x}x}_k-\mu \text{μ}\mu )({x\text{x}x}_k-\mu \text{μ}\mu {)}^{\top }]\\ & =-\frac{m}{2}\ln [(2\pi {)}^m|\Sigma \text{Σ}\Sigma |]-\frac{1}{2}\text{tr}[{\Sigma \text{Σ}\Sigma }^{-1}\sum _{k=1}^m({x\text{x}x}_k-\mu \text{μ}\mu )({x\text{x}x}_k-\mu \text{μ}\mu {)}^{\top }]\\ & =-\frac{m}{2}\ln [(2\pi {)}^m|\Sigma \text{Σ}\Sigma |]-\frac{1}{2}\text{tr}[{\Sigma \text{Σ}\Sigma }^{-1}S]\end{array}$$ 其中 $S={\sum }_{k=1}^m({x\text{x}x}_k-\mu \text{μ}\mu )({x\text{x}x}_k-\mu \text{μ}\mu {)}^{\top }$ 是一个对称矩阵，有 $S^{\top }=S$
• 对 $\mu \text{μ}\mu$ 求偏导数
  $$\begin{array}{rl}\frac{\partial lnL(\mu \text{μ}\mu ,\Sigma \text{Σ}\Sigma |x\text{x}x)}{\partial \mu \text{μ}\mu }& =\frac{\partial -\frac{1}{2}{\sum }_{k=1}^m({x\text{x}x}_k-\mu \text{μ}\mu {)}^{\top }{\Sigma \text{Σ}\Sigma }^{-1}({x\text{x}x}_k-\mu \text{μ}\mu )}{\partial \mu \text{μ}\mu }\\ & =-1\ast (-\frac{1}{2})\sum _{k=1}^{m}2{\Sigma \text{Σ}\Sigma }^{-1}({x\text{x}x}_t-\mu \text{μ}\mu )\\ & =\sum _{k=1}^m{\Sigma \text{Σ}\Sigma }^{-1}({x\text{x}x}_t-\mu \text{μ}\mu )\\ & ={\Sigma \text{Σ}\Sigma }^{-1}\sum _{k=1}^m({x\text{x}x}_t-\mu \text{μ}\mu )\end{array}$$ 对 $\Sigma \text{Σ}\Sigma$ 求偏导数
  $$\begin{array}{rl}\frac{\partial lnL(\mu \text{μ}\mu ,\Sigma \text{Σ}\Sigma |x\text{x}x)}{\partial \Sigma \text{Σ}\Sigma }& =\frac{\partial -\frac{m}{2}\ln [(2\pi {)}^m|\Sigma \text{Σ}\Sigma |]-\frac{1}{2}\text{tr}[{\Sigma \text{Σ}\Sigma }^{-1}S]}{\partial \Sigma \text{Σ}\Sigma }\\ & =-\frac{m}{2}{\left({\Sigma \text{Σ}\Sigma }^{\top }\right)}^{-1}+\frac{1}{2}{\Sigma \text{Σ}\Sigma }^{-1}{S}^{\top }{\Sigma \text{Σ}\Sigma }^{-1}\\ & =\frac{1}{2}[-m{\Sigma \text{Σ}\Sigma }^{-1}+{\Sigma \text{Σ}\Sigma }^{-1}S{\Sigma \text{Σ}\Sigma }^{-1}]\end{array}$$
• 令两个偏导数为 0，得到极大似然估计值
  $$\{\begin{array}{rl}\widehat{\mu \text{μ}\mu }& =\overline{x\text{x}x}=\frac{1}{m}\sum _{k=1}^m{x\text{x}x}_k\\ \widehat{\Sigma \text{Σ}\Sigma }& =\frac{1}{m}S=\frac{1}{m}\sum _{k=1}^m({x\text{x}x}_k-\widehat{\mu \text{μ}\mu })({x\text{x}x}_k-\widehat{\mu \text{μ}\mu }{)}^{\top }\end{array}$$ 和一维情况类似，这里 $\widehat{\mu \text{μ}\mu }$ 是绝对无偏估计值，$\widehat{\Sigma \text{Σ}\Sigma }$ 是渐进无偏估计值，$\Sigma \text{Σ}\Sigma$ 的绝对无偏估计值为
  $$\frac{1}{m-1}\sum _{k=1}^m({x\text{x}x}_k-\widehat{\mu \text{μ}\mu })({x\text{x}x}_k-\widehat{\mu \text{μ}\mu }{)}^{\top }$$
• 最后，注意协方差矩阵的估计形式为 $\widehat{\Sigma \text{Σ}\Sigma }=\frac{1}{m}{\sum }_{k=1}^m({x\text{x}x}_k-\widehat{\mu \text{μ}\mu })({x\text{x}x}_k-\widehat{\mu \text{μ}\mu }{)}^{\top }$，它求平均的每一项都是 $n\times 1$ 的向量和自身的转置相乘，得到的结果中每一个位置的元素都是独立计算的，即任意特征间的方差/协方差都是独立计算的，如果我们假设 $n$ 个特征全部独立，相当于增加了一个非对角位置全部为 $0$ 的约束，只需把协方差矩阵非对角线位置元素全部置为 $0$ 满足约束即可，这时退化为 $n$ 个一维高斯分布情况

3.4 程序演示

• 下面演示二维高斯分布的极大似然估计，本题改编自 《模式分类（第二版）》p.127 Prob.1，先搞点数据

  import numpy as np
  import matplotlib.pylab as plt
  
  w = np.array([[0.83,1.6],
                 [1.1,1.6],
                 [-0.44,-0.41],
                 [0.047,-0.45],
                 [0.28,0.35],
                 [-0.39,-0.48],
                 [-0.34,-0.079],
                 [-0.3,-0.22],
                 [1.1,1.2],
                 [0.18,-0.11]])
  
  fig = plt.figure(figsize = (5,5))
  a0 = fig.add_subplot(1,1,1,label='a0') 
  a0.scatter(w[:,0],w[:,1],c='red',alpha=0.5)
  
  mu = np.mean(w,axis=0)
  cov = np.dot((w-mu).T,(w-mu))/N
  
  print('期望的极大似然估计为\n{}'.format(np.around(mu,4)))
  print('方差的极大似然估计为\n{}'.format(np.around(cov,4)),end='\n\n')
  
  '''
  期望的极大似然估计为
  [0.2067 0.3001]
  方差的极大似然估计为
  [[0.3346 0.4305]
   [0.4305 0.645 ]]
  '''
  

• 分别假设数据服从二维联合高斯分布和二维独立的高斯分布（把协方差矩阵除了对角线以外清零），结果如下

  %matplotlib notebook
  import numpy as np
  import matplotlib.pylab as plt
  from scipy import stats
  
  plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
  plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
  
  fig = plt.figure(figsize = (10,10))
  a0 = fig.add_subplot(2,2,1,label='a0',projection='3d') 
  a1 = fig.add_subplot(2,2,2,label='a1',projection='3d') 
  a2 = fig.add_subplot(2,2,3,label='a2') 
  a3 = fig.add_subplot(2,2,4,label='a3') 
  
  w = np.array([[0.83,1.6],
                 [1.1,1.6],
                 [-0.44,-0.41],
                 [0.047,-0.45],
                 [0.28,0.35],
                 [-0.39,-0.48],
                 [-0.34,-0.079],
                 [-0.3,-0.22],
                 [1.1,1.2],
                 [0.18,-0.11]])
  N = w.shape[0]
  
  mu = np.mean(w,axis=0)
  cov = np.dot((w-mu).T,(w-mu))/N
  
  print('期望的极大似然估计为\n{}'.format(np.around(mu,4)))
  print('方差的极大似然估计为\n{}'.format(np.around(cov,4)),end='\n\n')
  
  x, y = np.mgrid[-2:2.5:.1, -2:2.5:.1]
  pos = np.empty(x.shape + (2,))
  pos[:, :, 0] = x; pos[:, :, 1] = y
  
  rv = stats.multivariate_normal(mu, cov)   # 生成多元正态分布
  a0.plot_surface(x, y, rv.pdf(pos),alpha=0.5,cmap=plt.cm.cool)
  a0.scatter(w[:,0],w[:,1],np.zeros(w.shape[0]),c='green')
  a0.set_title('联合估计')
  a2.contourf(x, y, rv.pdf(pos))         # 等高线
  a2.grid(alpha=0.5)                     # 坐标网格
  a2.scatter(w[:,0],w[:,1],c='red',alpha=0.5)
  
  mask = np.eye(2,2,dtype = 'bool')
  cov[~mask] = 0
  
  rv = stats.multivariate_normal(mu, cov)   # 生成多元正态分布
  a1.plot_surface(x, y, rv.pdf(pos),alpha=0.5,cmap=plt.cm.cool)
  a1.scatter(w[:,0],w[:,1],np.zeros(w.shape[0]),c='green')
  a1.set_title('独立估计')
  a3.contourf(x, y, rv.pdf(pos))         # 等高线
  a3.grid(alpha=0.5)                     # 坐标网格
  a3.scatter(w[:,0],w[:,1],c='red',alpha=0.5)
  

  
  可见增加各维度独立假设后，相当于把原先的超椭圆摆正了，从 x/y 轴方向看，超椭圆的宽度和长度是一致的。关于多维高斯分布的详细介绍，参考：随机过程（2.2）—— 多维高斯分布