904数据结构与机器学习考试复习[机器学习部分]——贝叶斯分类器
904数据结构与机器学习考试复习[机器学习部分]——贝叶斯分类器
Nothing is more practical than a good theory
(一)贝叶斯决策论
(二)极大似然估计
(三)朴素贝叶斯分类器
(四)EM算法
(零)前置知识
0 基本概念
0.1 条件概率公式
设\(A,B\)是两个事件，且\(P(B)>0\)，则在事件\(B\)已经发生的条件下，事件\(Ａ\)发生的概率，称为事件Ｂ发生下事件Ａ的条件概率(conditional probability)。
其基本求解公式为：
0.2 乘法公式
0.2.1 乘法公式的形式
由条件概率公式可得:
上式即为乘法公式
0.2.2 乘法公式的推广
对于任何正整数\(n\geq2\)，当\(P(A_1A_2...A_{n-1})>0\)时，有
0.3 全概率公式
0.3.1 全概率公式
(1)如果事件组\(B_1,B_2,...,B_n\)两两互斥，即\(B_i \cap B_j = \varnothing,i\neq j\),\(i,j=1,2,...n\),且\(P(B_i)>0\),\(i=1,2,...,n\);
(2)\(B_1\cup B_2=...=\Omega\)，则称事件组\(B_1,B_2,...,B_n\)是样本空间\(\Omega\)一个划分
设\(B_1,B_2,...,B_n\)是样本空间\(\Omega\)的一个划分，A为任一事件，则
上式即为全概率公式(formula of total probability)
0.3.2 全概率公式的意义
全概率公式的意义在于，当直接计算\(P(A)\)较为困难，而\(P(B_i)，P(A|B_i)(i=1,2,...,n)\)的计算较为简单时，可以利用全概率公式计算\(P(A)\)。其思想就是，将事件Ａ分解成几个小事件，通过求小事件的概率，然后相加从而求得事件Ａ的概率，而将事件Ａ进行分割的时候，不是直接对Ａ进行分割，而是先找到样本空间\(\Omega\)的一个个划分\(B_1,B_2,...,B_n\),这样事件Ａ就被事件\(AB_1,AB_2,...,AB_n\)分解成了n部分，即\(A=AB_1,AB_2,...,AB_n\),每一\(B_i\)发生都可能导致Ａ发生相应的概率是\(P(A|B_i)\)，由加法公式得
例题1 计算车间次品率
某车间用甲、乙、丙三台机床进行生产，各台机床次品率分别为5%，4％，2％,它们各自的产品分别占总量的25％，35％，40％，将它们的产品混在一起，求任取一个产品是次品的概率。
解 设P(A)为任取一个产品都为次品的概率，则 P(A) = 25%×5% + 4%×35% + 2%×40% = 0.0345
0.4 贝叶斯公式
与全概率公式解决的问题相反，贝叶斯公式是建立在条件概率的基础上寻找事件发生的原因（即大事件\(A\)已经发生的条件下，分割中的小事件\(B_i\)的概率），设\(B_1,B_2,...,B_n\)是样本空间\(\Omega\)的一个划分，则对任一事件\(A\)（\(P(A)>0\)),有
上式即为贝叶斯公式(Bayes formula)，\(B_i\)常被视为导致试验结果Ａ发生的“原因”，
\(P(B_i)(i=1,2,...,n)\)表示各种原因发生的可能性大小，故称为先验概率；
\(P(B_i|A)(i=1,2,...,n)\)则反映当试验产生了结果Ａ之后，再对各自原因概率的新认识，故称为后验概率。
例题2 吸毒者检测
贝叶斯定理在检测吸毒者时很有用。假设一个常规的检测结果的敏感度与可靠度均为99%，也就是说，当被检者吸毒时，每次检测呈阳性（+）的概率为99%。而被检者不吸毒时，每次检测呈阴性（-）的概率为99%。从检测结果的概率来看，检测结果是比较准确的，但是贝叶斯定理却可以揭示一个潜在的问题。假设某公司将对其全体雇员进行一次鸦片吸食情况的检测，已知0.5%的雇员吸毒。我们想知道，每位医学检测呈阳性的雇员吸毒的概率有多高。令“D”为该公司雇员吸毒事件，“N”为该公司雇员不吸毒事件，“+”为该公司雇员检测呈阳性事件。可得
- P(D)代表雇员吸毒的概率，不考虑其他情况，该值为0.005。因为公司的预先统计表明该公司的雇员中有0.5%的人吸食毒品，所以这个值就是D的先验概率。
- P(N)代表雇员不吸毒的概率，显然，该值为0.995，也就是1-P(D)。
- P(+|D)代表吸毒者阳性检出率，这是一个条件概率同时也是先验概率，由于阳性检测准确性是99%，因此该值为0.99。
- P(+|N)代表不吸毒者阳性检出率，也就是出错检测的概率，该值为0.01，因为对于不吸毒者，其检测为阴性的概率为99%，因此，其被误检测成阳性的概率为1-99%。
- P(+)代表不考虑其他因素的影响的阳性检出率。该值为0.0149或者1.49%。我们可以通过全概率公式计算得到：此概率 = 吸毒者阳性检出率(0.5% × 99% = 0.00495)+ 不吸毒者阳性检出率(99.5% × 1% = 0.00995)。P(+)=0.0149是检测呈阳性的先验概率。用数学公式描述为：
根据上述描述，我们可以计算某人检测呈阳性时确实吸毒的条件概率P(D|+)：
尽管我们的检测结果可靠性很高，但是只能得出如下结论：如果某人检测呈阳性，那么此人是吸毒的概率只有大约33%，也就是说此人不吸毒的可能性比较大。我们测试的条件（本例中指D，雇员吸毒）越难发生，发生误判的可能性越大。
但如果让此人再次复检（相当于P(D)=33.2215%，为吸毒者概率，替换了原先的0.5%），再使用贝叶斯定理计算，将会得到此人吸毒的概率为98.01%。但这还不是贝叶斯定理最强的地方，如果让此人再次复检，再重复使用贝叶斯定理计算，会得到此人吸毒的概率为99.98%（99.9794951%）已经超过了检测的可靠度
(一)贝叶斯决策论
贝叶斯决策论(Bayesian decision theoty)是概率框架下实施决策的基本方法。
对于分类任务来说，在所有相关概率都已知的理想情形下，贝叶斯决策论考虑如何基于这些概率和误判损失来选择最优的类别标记。
下面我们以多分类任务为例来解释其基本原理。
假设有N种可能的类别标记，即\(\mathcal{Y} = \left\{c_1,c_2,...,c_N\right\}\)，\(\lambda_{ij}\)是将一个真实标记为\(c_j\)的样本误分类为\(c_i\)所产生的损失。基于后验P(c_i|x)可获得将样本\(x\)分类为\(c_i\)所产生的期望损失(expected loss)，即在样本\(x\)上的“条件风险”(conditional risk)。
等式(7.1)的解释
等式左侧\(R(c_i|x)\)表示将样本\(x\)分类为\(c_i\)所产生的期望损失；这也就是说，现在已知样本\(x\)被分类为\(c_i\)，想要知道的是这一事实的期望损失有多少。
那么，根据样本\(x\)真实类别标记的不同，\(x\)被分类为\(c_i\)所产生的损失肯定也不同。当样本\(x\)真实类别标记为\(c_1\)时候误分类为\(c_i\)所产生的损失为\(\lambda_{i1}\)，当样本\(x\)真实标记为\(c_{2}\)时误分类为\(c_i\)所产生的损失为\(\lambda_{i2},...\)，当样本\(x\)真实类别标记为\(c_2\)时误分类为\(c_i\)所产生的损失为\(\lambda_{i2}，...\)，当样本\(x\)真实类别标记为\(c_j\)时误分类为\(c_i\)所产生的损失为\(\lambda_{ij}.....\)。
因此，只要知道样本\(x\)真实类别标记，就可以知道将其分类为\(c_i\)所产生的损失。现在待求的条件风险\(R(c_i|x)\)是将样本\(x\)分类为\(c_i\)所产生的期望损失，即所有损失的平均值。根据数据期望的定义，每个损失乘以其对应的概率，再求和即可。
对于样本\(x\)来说，它的标记为\(c_j\)的概率为\(P(c_j|x)\)，即已知\(x\)的情况下，类别标记为\(c_j\)的后验概率(在给定\(x\)的条件下，类别标记为\(c_j\)的(条件)概率)，因此得式(7.1)：
\(R(c_i|x) = \sum \limits_{j=1}^N \lambda_{ij}P(c_j|x)\).
其中\(N\)为可能的类别标记个数。
我们的任务是寻找一个判定准则\(h:\mathcal{X}\longmapsto \mathcal{Y}\)以最小化总体风险
等式(7.2)的解释
式(7.1)是针对单个样本\(x\)的，而式(7.2)则是针对整个数据集D所有样本的期望，即
\(R(h) = \mathbb{E}_x[R(h(x))|x] = \sum \limits_{x \in D} R(h(x)|x)P(x)\).
其中\(P(x)\)为样本\(x\)出现的概率，且满足\(\sum_{x \in D}P(x) = 1\);\(h(x)\)为判定准则\(h: \mathcal{X} \mapsto \mathcal{Y}\)对样本\(x\)预测的类别标记；则\(R(h)\)表示判定准则\(h\)的总体风险。
显然，对每个样本\(x\)，若h能最小化条件风险\(R(h(x)|x)\)，则总体风险\(R(h)\)也将被最小化。
这就产生了贝叶斯判定准则(Bayes decision rules)：为最小化总体风险，只需在每个样本上选择那个能使条件风险\(R(c|x)\)最小的类别标记，即
式(7.3)的解释
这里解释以下符号"arg min"，其中"arg"是"argument"的前三个字母，“min”是"minimum"的前三个字母。维基百科中有对"arg max"的解释:参考连接。是其反义符号。概况起来，式(7.3)表示求出使目标函数\(R(c|x)\)最小的类别标记\(c\)，并将该\(c\)返回给\(h^*(x)\)作为输出，注意此处\(x\)为已知常量，\(R(c|x)\)表达式参见式(7.1)。
此式共涉及三个概念：贝叶斯判定准则(Bayes decision rule)，贝叶斯最优分类器(Bayes optimal classifier)，贝叶斯风险(Bayes risk)。式(7.3)即为贝叶斯判定准则，所得分类器\(h^*(x)\)为贝叶斯最优分类器，此时的总体风险\(R(h^*)\)称为贝叶斯风险。
此时，\(h^*\)称为贝叶斯最优分类器(Bayes)，与之对应的总体风险\(R(h^*)\)称为贝叶斯风险(Bayes risk)。
\(1-R(h^*)\)反映了分类器所能达到的最好性能，即通过机器学习所能产生的模型精度的理论上限。
具体来说，若目标是最小化分类错误率，则误判损失\(\lambda_{ij}\)可写为
式(7.4)的解释
该式所表达的即为0/1损失，也就是说分类正确时(i=j)损失为0，否则损失为1。
此时条件风险
式(7.5)的解释
将式(7.4)代入式(7.1)，得
\(\begin{align} &R(c_i|x) = \sum \limits_{j = 1}^N \lambda_{ij}P(c_j|x) \\ &= \sum \limits_{j=1}^{i-1} P(c_j|x) + \sum \limits_{j = i+1}^N P(c_j|x) \\ &= \sum \limits_{j=1}^{i-1} P(c_j|x) + P(c_i|x) + \sum \limits_{j = i+1}^N P(c_j|x) - P(c_j|x) \\ &= \sum \limits_{j=1}^N P(c_j|x) - P(c_i|x) \\ &= 1 - P(c_i|x) \end{align}\)
其中第二个等号利用了式(7.4)，最后一个等号利用了\(\sum_{j=1}^N P(c_j|x) = 1\)。
于是，最小化分类错误率的贝叶斯最优分类器为
式(7.6)的解释
式(7.3)的贝叶斯最优分类器是最小化式(7.5)，因此等价于最大化\(P(c|x)\)。
即对每个样本\(x\)，选择能使后验概率\(P(c|x)\)最大的类别标记。
不难看出，欲使用贝叶斯判定准则来最小化决策风险，首先要获得后验概率\(P(c|x)\)。
然而，在现实任务中这通常难以直接获得。从这个角度来看，机器学习所要实现的是基于有限的训练样本集尽可能准确地估计出后验概率\(P(c|x)\)。大体来说，主要有两种策略：
- 给定\(x\)，可通过直接建模\(P(c|x)\)来预测\(c\)，这样得到的是“判别式模型”(discriminative models)；
- 先对联合概率分布\(P(x,c)\)建模，然后再由此获得\(P(c|x)\)，这样得到的是“生成式模型”(generative models)。
判别式模型与生成式模型
可以简单理解为判别式模型直接建模后验概率\(P(c|x)\)来预测类别c，而生成式模型则是对联合概率\(P(x,c)\)建模来预测类别\(c\)。
具体来说：
对于判别式模型来说，就是已知\(x\)的条件下判别其类别标记，即求后验概率\(P(c|x)\)，显然，我们已知的决策树、BP神经网络、支持向量机等，都可归入判别式模型的范畴。其中尤其以3.3节介绍的对率回归最为直接，详见式(3.23)和式(3.24)。
对于生成式模型来说，理解起来比较抽象。思考两个问题：
(1)对于数据集\(D\)来说，其中的样本是如何生成的呢?
答：当然是按照联合概率分布\(P(x,c)\)采样而得，也可以描述为根据\(P(x,c)\)生成的。(2)若已知样本\(x\)和联合概率分布\(P(x,c)\)，如何预测类别\(c\)呢？
答：在样本\(x\)已知的情况下，可以分别求出\(P(x,c_1),P(x,c_2),P(x,c_N)\)，即生成样本\((x,c_1),(x, c_2), ...,(x,c_N)\)的概率(特征\(x\)已知，可变部分只有类别\(c\))，当然是选择与\(x\)联合概率最大的类别标记，即为 \(h^*(x) = \arg \min \limits_{c \in \mathcal{Y}} P(x,c)\)
因此，之所以称为“生成式”模型，是因为所求概率\(P(x,c)\)是样本生成的概率。
对生成式模型来说，必然考虑
基于贝叶斯定理，\(P(c|x)\)可写为
其中，\(P(c)\)是类“先验”(prior)概率；\(P(x|c)\)是样本\(x\)相当于类标记\(c\)的类条件概率(class-conditional probability)，或称为“似然”(likehood)；\(P(x)\)是用于归一化的“证据”(evidence)因子。
对给定样本\(x\)，证据因子\(P(x)\)与类标记无关，因此估计\(P(c|x)\)的问题就转化为如何基于训练数据Ｄ来估计先验\(P(c)\)和似然\(P(x|c)\)。
类先验概率\(P(c)\)表达了样本空间中各类样本所占的比例，根据大数定律，当训练集包含充足的独立同分布样本时，\(P(c)\)可通过各类样本出现频率来进行估计。
对类条件概率\(P(x|c)\)来说，由于它涉及关于\(x\)所有属性的联合概率，直接根据样本出现的频率来估计将会遇到严重的困难。例如，假设样本的\(d\)个属性都是二值的，则样本空间将有\(2^d\)种可能的取值，在现实应用中，这个值往往大于训练样本数\(ｍ\)，也就是说，很多样本取值在训练集中根本没有出现，直接使用频率来估计\(P(x|c)\)显然不可行，因为“未被观测到”与“出现概率为零”通常是不同的。
式(7.7)和式(7.8)的解释
\(P(c|x) ＝ \frac{P(x,c)}{P(x)} \quad \quad (7.7)\)
\(P(c|x) = \frac{P(c)P(x|c)}{P(x)} \quad \quad (7.8)\)这两个式子本质上相同，大学学过概率论的人都应该知道有以下的概率关系：
\(P(x,c) = P(x|c)P(c) = P(c|x)P(x)\)
前面在分析生成式模型的概念时提到，若已知样本\(x\)和联合概率分布\(P(x,c)\)，预测类别\(c\)时选择与\(x\)联合概率最大的类别标记，而这实际上与选择使后验概率\(P(c|x)\)最大的类别标记等价；因为给定样本\(x\)，概率\(P(x)\)为常量(而且概率肯定为大于零)，因此，最大化\(P(x,c)\)和\(\frac{P(x,c)}{P(x)}\)等价(常量并不影响极值点的存在位置，例如\(f(x) = x^2\)与\(g(x) = 3x^2\)极值点均为\(x = 0\))，即最大化后验概率\(P(c|x)\)与最大化联合概率\(P(x,c)\)等价。由式(7.8)，估计\(P(c|x)\)的问题转化为如何基于训练数据\(D\)来估计\(P(c)\)和\(P(x|c)\)。
以西瓜书的表4.1数据集2.0为例，基于训练数据\(D\)来估计\(P(c)\)相对容易
\(P(好瓜=是) = \frac{8}{17}\)，\(P(好瓜=否) = \frac{9}{17}\)
而基于训练数据\(D\)估计\(P(x|c)\)则比较困难。
若估计当\(x = (色泽＝青绿,根蒂=蜷缩，敲声=浊响，纹理=清晰,脐部=凹陷，触感=硬滑), c=(好瓜＝是)\)的概率\(P(x|c)\)，则在前8个“好瓜=是”的样本中，仅编号为1的瓜的特征等于\(x\)，此时若将\(P(x|c)\)估计为\(\frac{1}{8}\)，这实际上是不可靠的，因为此时的样本集合太小了；
再例如若\(x = (色泽=青绿，根蒂=蜷缩，敲声=浊响，纹理=清晰,脐部＝凹陷，触感＝软粘)\)，则在前8个“好瓜＝是”的样本中根本没出现此\(x\)，但“未被观测到”与“出现概率为零”通常是不同的。
这两个式子就是贝叶斯定理，它将估计后验概率\(P(c|x)\)转化为估计先验概率\(P(c)\)和似然概率\(P(x|c)\)，而这正式朴素贝叶斯分类器、半朴素贝叶斯分类器等的理论基础。
先验概率和条件概率、后验概率和似然概率
对于条件概率\(P(A|B)\)，表示事件Ａ在另外一个事件B已经发生下发生的概率，读作为“在B的条件下A的概率”；
而\(P(A)\)称为先验概率，即在无任何已知条件下事件A的概率。
若给定B，则\(P(A|B)\)称为后验概率；此时，对于多个候选A，选择使后验概率\(P(A|B)\)最大的那一个，称为最大后验(maximum a posteriori,MAP)估计。
若给定A，则\(P(A|B)\)称为似然概率；此时，对于多个候选B，选择使似然概率\(P(A|B)\)最大的那一个，称为最大似然(maximum likelihood,ML)估计。(即7.2节)
后验概率\(P(A|B)\)可以理解为已知\(B\)之后(即后验)，A发生的(条件)概率；而似然概率\(P(A|B)\)可以理解为已知\(A\)之后，估计它更像是(似然)在哪个B的条件下得到的。
一般来说，后验概率是\(P(果|因)\)，似然概率是\(P(因|果)\)，而并不像本条解释中的A和B任意指定已知哪个。
在机器学习任务中，已知样本特征向量\(x\)来预测样本类别\(y\)，因此\(P(y|x)\)是后验概率，\(P(x|y)\)是似然概率。
另外，\(P(A|B)\)和\(P(B|A)\)的差别是很大的，正如百度百科词条[条件概率]上所述，很多人经常会犯错误，误以为\(P(A|B)\)大致等于\(P(B|A)\)，此即条件概率的谬论。
例如，A表示事件喝了口开水，B表示事件被烫着了，则\(P(A|B)\)表示该人被烫着了的条件下喝口开水的概率(已知被烫着了，求此时喝口开水的概率)，而\(P(A|B)\)表示该人喝了口开水的条件下被烫着的概率(已知喝了口开水，求此时被烫着的概率)，显然这两个概率是不同的。
(二)极大似然估计
估计类条件概率的一种常用策略是先假定其具有某种确定的概率分布形式，再基于训练样本对概率分布的参数进行估计。具体地，记关于类别\(c\)的类条件概率为\(P(x|c)\)，假设\(P(x|c)\)具有确定的形式并且被参数向量\(\theta_c\)唯一确定，则我们的任务就是利用训练集D估计参数\(\theta_c\)。为明确起见，我们将\(P(x|c)\)记为\(P(x|\theta_c)\)。
事实上，概率模型的训练过程就是参数估计(parameter estiimation)过程。
对于参数估计，统计学界的两个学派分别提供了不同的解决方案：
-
频率主义学派认为(Frequentist)认为参数虽然未知，但却是客观存在的固定值，因此，可通过优化似然函数等准则来确定参数值；
-
贝叶斯学派(Bayesian)则认为参数是未观察到的随机变量，其本身也可有分布，因此，可假定参数服从一个先验分布，然后基于观测到的数据来计算参数的后验分布。
极大似然估计(Maximum Likelihood Estimation,简称MLE)，是源自于频率主义学派，根据数据采样来估计概率分布参数的经典方法。
在朴素贝叶斯法中，学习意味着估计\(P(Y=c_k)\)和\(P(X^{(j)} = x^{(j)}|Y=c_k)\)。可以应用极大似然估计法估计相应的概率。先验概率\(P(Y=c_k)\)的极大似然估计是
设第\(j\)个特征\(x^{(j)}\)可能取值的集合为{\(a_{j1},a_{j2},...,a_{j}S_j\)}，条件概率\(P(X^{(j)}=a_{jl}|Y=c_k)\)的极大似然估计是
式中，\(x_i^{(j)}\)是第\(i\)个样本的第\(j\)个特征，\(\alpha_{jl}\)是第\(j\)个特征可能取得的第\(l\)个值；\(I\)为指示函数。
令\(D_c\)表示训练集Ｄ中第\(c\)类样本组成的集合，假设这些样本是独立同分布的，则参数\(\theta_c\)对于数据集\(D_c\)的似然是
对\(\theta_c\)进行极大似然估计，就是去寻找能最大化似然\(P(D_c|\theta_c)\)的参数值\(\hat \theta_c\)。直观上看，极大似然估计是试图在\(\theta_c\)所有可能的取值中，找到一个能使数据出现的“可能性”最大的值。
\(P(D_c|\theta_c)=\underset{x∈D_c}{\prod}P(x|\theta_c)\)中的连乘操作易造成下溢，通常使用对数似然(log-likelohood)
此时参数\(\theta_c\)的极大似然估计\(\hat \theta_c\)为
例如，在连续属性情形下，假设概率密度函数\(p(x|c)～\mathcal{N}(\mu_c,\sigma_c^2)\)，则参数\(\mu_c\)和\(\sigma_c^2\)的极大似然估计为
也就是说，通过极大似然法得到的正态分布均值就是样本均值，方差就是\((x-\hat \mu_c)(x-\hat \mu_c)^T\)的均值，这显然是一个复合直觉的结果。在离散属性情形下，也可通过类似的方式估计类条件概率。
需注意的是，这种参数化的方法虽能使类条件概率估计变得相对简单，但估计结果的准确性严重依赖于所假设的概率分布形式是否符合潜在的真实数据分布。在现实应用中，欲做出能较好地接近潜在真实分布的假设，往往需在一定程度上利用关于应用任务本身的经验知识，否则若仅凭“猜测”来假设概率分布形式，很可能产生误导性的结果。
对于极大似然法，关键在于找出似然概率的表达式。式(7.9)给出了最原始的似然表达式，注意该式中\(\theta_{c}\)是待求模型参数，极大似然要做的是寻找当当前数据集更像是由哪组模型参数\(\theta_{c}\)生成的；连乘容易造成下溢，这是由于式(7.9)的概率均小于1，例如100个0.1连乘，则乘积非常小，此时计算机可能无法表达如此小的数(这里仅为距离，具体问题请具体分析)，此即为下溢(下溢出的反面为上溢，例如100个10连乘，则乘积非常大)；解决下溢的常用方法就是取对数，将连乘操作变为连加操作，此即为式(7.10)；得到了(对数)似然的表达式之后，在\(\theta_{c}\)的可行域内最大化似然函数即得极大似然估计，此即为式(7.11)。
例:投掷硬币5次，结果依次是正面、正面、反面、正面、反面，试基于此观察结果估计硬币正面朝上的概率\(\theta\)。
解：设正面朝上的概率为\(\theta\)，各次投掷结果相互独立，则似然为
对数似然为
求导并令导数等于零
解得\(\theta = \frac{3}{5}\)，也就是说，基于观测结果对概率\(\theta\)的极大似然估计就是正面出现的比例。
极大似然估计要解决的问题
- 给定一个数据分布 \(P_{data}(x)\).
- 给定一个由参数\(\theta\)定义的数据分布\(P_{G}(x;\theta)\).
- 我们希望，求得参数\(\theta\)使得\(P_{G}(x;\theta)\)尽可能接近\(P_{data}(x)\).
可以理解为：
\(P_{G}(x;\theta)\)是某一具体的分布(比如简单的高斯分布)，而\(P_{data}(x)\)是未知的(或者及其复杂，我们可能找到一个方式去表示它)，因此我们希望通过极大似然估计的方法来确定\(\theta\)，让\(P_{G}(x;\theta)\)能大体表示\(P_{data}(x)\).
极大似然估计的解决方案
1.从\(P_{data}(x)\)采样\(m\)个样本{\(x^{1},x^{2},...,x^{m}\)}.
2.计算采样样本的似然函数\(L=\prod \limits_{i=1}^m P_G(x^i;\theta)\).
3.计算使得似然函数\(L\)最大的参数\(\theta\)：
这里再讨论一下极大似然估计这么做的必要性：
\(P_{data}(x)\)可以理解成是非常复杂的分布，不可能用某个数学表达精确表示，因此我们只能通过抽象，使用一个具体的分布模型\(P_G(x;\theta)\)近似\(P_{data}(x)\)。
所以，求\(P_G(x;\theta)\)的参数\(\theta\)的策略就变成了：
我们认为来自\(P_{G}(x;\theta)\)的样本{\(x^{1},x^{2},...,x^{m}\)}在\(P_G(x^i;\theta)\)分布中出现的概率越高，也就是\(\prod \limits_{i=1}^m P_G(x^i;\theta)\)越大，
\(P_G(x;\theta)\)和\(P_{data}(x)\)就越接近。
因此，我们期待的\(\theta\)就是使得\(\prod \limits_{i=1}^m P_G(x^i;\theta)\)最大的\(\theta\).
即：
例题3 身高统计
抽样某中学高一年级的200名男生，统计身高数据为：
step 1. 采样
我们将该校或者该地区所有高一男生的身高理解成某一分布\(P_{data}(x)\) ，以上数据便是我们从这一分布中的采样。
样本数据表示为：\(\left\{x^{(1)},x^{(2)},...,x^{(200)}\right\}\).
step 2. 建立极大似然函数
这个分布肯定是及其复杂的，会受到众多因素影响，用数学表达精准表达这个分布完全不可能！
但我们可以通过一个数学上已知并且明确的简单分布来抽象近似这个\(P_{data}(x)\).
我们期望使用数学模型解决实际问题时，第一步就是化简实际问题，去掉实际问题中很多琐碎、影响不是很大的因素，然后使用明确的数学模型去抽象近似实际问题。
这种例子实在太多了，比如我们将太阳抽象成一个完美的球体，在力学问题中将物体抽象成质点等。
根据直觉，我们认为身高分布可能非常近似高斯分布，我们不妨就把身高分布抽象成高斯分布吧！
高斯分布
因此，我们可以据此建立极大似然函数
step 3. 求解参数
求解使\(L\)尽可能大的参数\(\mu,\sigma\).
此时我们认为\(P_G\)和\(P_{data}\)是最接近的。
即：
具体求解方法：
因为这里\(L\)是连乘，并且有\(exp\)的存在，所以很自然地想到对数\(log\)来进行化简
即：
由此可得
不妨令
我们的目标是
下面选择求偏导求多元函数的最值
解得
我们发现，\(\mu^*\)实际就是样本均值，\(\sigma^*\)实际也就是样本标准差。这也符合高斯分布的特点。
接下来，将\({}\)\(\left\{x^{(1)},x^{(2)},...,x^{(200)}\right\}\)具体数值代入计算可得
由此，便得到目标的\(P_G(x;\mu^*=170.08, \sigma^* = 5.04)\).
我们将样本在不同区间出现的频率(图中散点)与\(P_G(x;\mu^*=170.08,\sigma^*=5.04)\)的概率密度函数(图中曲线)绘制出来，如下图：
发现大致散点与曲线，也即频率与频率密度函数曲线是契合的，这也进一步验证了我们的模型复合。
(三)朴素贝叶斯分类器
3.1 朴素贝叶斯分类器
根据以上讨论，我们不难发现，基于贝叶斯公式\(P(c|x) = \frac{P(c)P(x|c)}{P(x)}\)来估计后验概率\(P(c|x)\)的主要困难在于：类条件概率\(P(x|c)\)是所有属性上的联合概率，难以从有限的训练样本直接估计而得。
为了避开这个障碍，朴素贝叶斯分类器(naive Bayes classifier)采用了“属性条件独立性假设”(attribute conditional independence assumption):对已知类别，假设所有属性相互独立。换言之，假设每个属性独立地对分类结果发生影响。
基于属性条件独立性假设，\(P(c|x) = \frac{P(c)P(x|c)}{P(x)}\)可重写为
其中\(d\)为属性数目，\(x_i\)为\(x\)在第\(i\)个属性上的取值。
由于对所有类别来说\(P(x)\)相同，因此基于式\(h^*(x) =\underset{c∈\mathcal{Y}}{arg max} R(c|x)\)的贝叶斯判定准则可得，朴素贝叶斯分类器的表达式：
式(7.15)的解释
朴素贝叶斯分类器即式(7.15)，相比于式(7.14)少了\(P(x)\)，这是因为\(P(x)\)与类标记\(c\)无关，而常量并不影响极值点的存在位置，例如\(f(x) = x^2\)和\(g(x) = 3x^2\)极小值点位置相同，均为\(x = 0\)。式(7.15)是多个概率连乘的形式，可以取对数转化为连加的形式以防止下溢。
显然，朴素贝叶斯分类器的训练过程就是基于训练集D来估计类先验概率P(c)，并为每个属性估计条件概率\(P(x_i|c)\)。
令\(D_c\)表示训练集\(D\)中第\(c\)类样本组成的集合，若有充足的独立同分布样本，则可容易地估计出类先验概率：
对离散属性而言，令\(D_{c,x_i}\)表示\(D_c\)中在第\(i\)个属性上取值为\(x_i\)的样本组成的集合，则条件概率\(P(x_i|c)\)可估计为
离散变量的概率估计较为简单，直接按频率计数即可，如式(7.16)的先验概率估计和式(7.17)的条件概率估计，注意式中的\(|\cdot|\)表示求集合元素的个数。
对于连续属性可考虑概率密度函数，假定\(p(x_i|c)～\mathcal{N}(\mu_{c,i}, \sigma_{c,i}^2)\)，其中\(\mu_{c,i}\)和\(\sigma_{c,i}^2\)分别是第\(c\)类样本在第\(i\)类样本在第\(i\)个属性上取值的均值和方差，则有
而连续变量必须假设一种概率频率分布，常用的就是高斯分布；注意：此处概率密度并非概率，即取值范围不是0到1之间，对于高斯概率密度函数来说，式(7.18)的\(p(x_i|c) \in (0, \frac{1}{\sqrt{2 \pi}\sigma_c}]\)，其中当\(x_i = \mu_{c,i}\)时取得最大值\(\frac{1}{\sqrt{2\pi}}\sigma_{c,i}\)，进一步地当\(\sigma_{c,i} < \frac{1}{\sqrt{2\pi}}\sigma_{c,i}\)时，\(\frac{1}{\sqrt{2\pi}\sigma_{c,i}}>1\)，即会出现\(p(x_i|c) > 1\)的情形。
在离散变量的概率估计中，由于训练集样本的不充分性，可能会出现分子为零的情形，但“未被观测到”与“出现概率为零”是不同的，因此要进行平滑(smoothing)，常用的是拉普拉斯修正。
算法1 (朴素贝叶斯算法)
输入：训练数据\(T= \left\{(x_1,y_1),(x_2,y_2),...(x_N,y_N)\right\}\)，其中\(x_i = (x_i^{(1)},x_i^{(2)},...,x_i^{(n)})^T\),\(x_i^{(j)}\)是第\(i\)个样本的第\(j\)个特征，\(x_i^{j}∈\left\{a_{j1},a_{j2},...,a_{jS_j}\right\}\),\(a_{jl}\)是第\(j\)个特征可能取的第\(l\)个值，\(j=1,2,...,n, \quad l=1,2,...,S_j,\quad y_i∈ \left\{c_1,c_2,...,c_k \right\}\)；实例\(x\);
输出：实例\(x\)的分类。
(1) 计算先验概率及条件概率
(2)对于给定的实例\(x=(x^{(1)},x^{(2)},...,x^{(n)})^T\)，计算
(3)确定实例\(x\)的类
例题4 通过朴素贝叶斯分类器确定类标记\(y\)
试由下表的训练数据学习一个朴素贝叶斯分类器并确定\(x=(2,S)^T\)的类标记\(y\)。
表中\(X^{(1)},X^{(2)}\)为特征，取值的集合分别为\(A_1= \left\{1,2,3\right\}\),\(A_2= \left\{S,M,L\right\}\)，
\(Y\)为类标记，\(Y∈C= \left\{1,-1\right\}\)。
解 根据朴素贝叶斯算法，由上表，容易计算下列概率：
对给定的\(x=(2,S)^T\)计算：
因为\(P(Y=-1)P(X^{(1)}=2|Y=-1)P(X^{(2)}=S|Y=-1)\)最大，所以\(y=-1\)。
例题5 判别是否为好瓜
下面用西瓜数据集训练一个朴素贝叶斯分类器，对测试例“测1”进行分类
首先估计类先验概率\(P(c)\)，显然有
然后，为每个属性估计条件概率\(P(x_i|c)\)：
根据概率乘法公式，分别计算是否为好瓜的概率，可得
由于，\(0.038>6.80×10^{-5}\)，因此，朴素贝叶斯分类器将测试样本“测1”判别为“好瓜”。
需注意的是，若某个属性值在训练集中没有与某个类同时出现过，则直接基于\(P(x_i|c) = \frac{|D_{c,x_i}|}{|D_c|}\)进行判别将出现问题。例如，在使用西瓜数据集训练朴素贝叶斯分类器时，对一个“敲声=清脆”的测试例，有
由于朴素贝叶斯表达式\(h_{nb}(x) = \underset{c∈\mathcal{Y}}{arg max}P(c) \prod \limits_{i=1}^d P(x_i|c)\)中的连乘式计算出的概率值为零，因此，无论该样本的其他属性值是树木，哪怕在其他属性上明显像好瓜，分类的结果都将是“好瓜=否”，这显然不合理。
为避免其他属性携带的信息被训练集中未出现的属性值“抹去”，在估计概率值时通常要进行“平滑”(smoothing)，常用“拉普拉斯修正”(Laplacian correction)。具体来说，令N表示训练集D中可能的类别数，\(N_i\)表示第\(i\)个属性可能的取值数，则式(7.16)和式(7.17)分别修正为
例如，在上面的例子总，类先验概率可估计为
类似地，\(P_{青绿|是}\)和\(P_{青绿|否}\)以及\(P_{清脆|是}\)可估计为
显然，拉普拉斯修正避免了训练集样本不充分而导致概率估值为零的问题，并且在训练集变大时，修正过程所引入的先验(prior)的影响也会逐渐变得可忽略，使得渐趋向于实际概率值。
在现实任务中朴素贝叶斯分类器有多种使用方式。例如，若任务对预测速度要求较高，则对给定训练集，可将朴素贝叶斯分类器涉及的所有概率值估计事情先计算好存储起来，这样在进行预测时只需“查表”即可进行判别；若任务数据更替频繁，则可以采用“懒惰学习”(lazy learning)方式，先不进行任何训练，待收到预测请求时再根据当前数据集进行概率估值；若数据不断增加，则可以在现有估值基础上，仅对新增样本的属性值所涉及的概率估值进行计数修正即可实现增量学习。
3.2 贝叶斯估计
用极大似然估计可能会出现所要估计的概率值为0的情况。这时会影响到后验概率的计算结果，使分类产生偏差。解决这一问题的方法是采用贝叶斯估计。具体地，条件概率的贝叶斯估计是
其中\(\lambda\geq0\).等价于在随机变量各个取值的频数上赋予一个正数\(\lambda>0\)。
当\(\lambda=0\)时就是极大似然估计。常取\(\lambda=1\)，这时称为拉普拉斯平滑(Laplacian smoothing)。
显然，对任何\(l=1,2,...,S_j, k=1,2,...,K\)，有
同样的，先验概率的贝叶斯估计是
为了避免其他属性携带的信息被训练集中未出现的属性“抹去”，在估计概率值时通常要进行“平滑”(smoothing)，常用“拉普拉斯修正”(Laplacian correction)。具体来说，令\(N\)表示训练集\(D\)中可能的类别数，\(N_i\)表示第\(i\)个属性可能的取值数，则类先验概率\(P(c)=\frac{|D_{c,x_i}|}{|D_c|}\)和条件概率\(P(x_i|c) = \frac{|D_{c,x_i}|}{|D_c|}\)分别修正为
例如，在上面的西瓜数据例子中，类先验概率可估计为
类似地，\(P_{青绿|是}\)和\(P_{青绿|否}\)可以估计为
同时，概率\(\hat P_{清脆|是}\)可估计为
显然，拉普拉斯修正避免了因训练集样本不充分而导致概率估值为零的问题，并且在训练集变大时，修正过程所引入的先验(prior)的影响也会逐渐变得可忽略，使得估值渐趋向于实际概率值。
在现实任务中，朴素贝叶斯分类器有多种使用方式。
例如，若任务对预测速度要求较高，则对给定训练集，可将朴素贝叶斯分类器涉及的所有概率估值事先计算好存储起来，这样在进行预测时只需“查表”即可进行判别；
若任务数量更替频繁，则可采用“懒惰学习”(lazy learning)方式，先不进行任何训练，待收到预测请求时再根据当前数据集进行概率估值；
若数据不断增加，则可以在现有估值基础上，仅对新增样本的属性值所涉及的概率估值进行计数修正即可实现增量学习。
例题6 按照拉普拉斯平滑估计概率，即取\(\lambda=1\).
试由下表的训练数据学习一个朴素贝叶斯分类器并确定\(x=(2,S)^T\)的类标记\(y\)。
表中\(X^{(1)},X^{(2)}\)为特征，取值的集合分别为\(A_1= \left\{1,2,3\right\}\),\(A_2= \left\{S,M,L\right\}\)，
\(Y\)为类标记，\(Y∈C= \left\{1,-1\right\}\)。
解 \(A_1= \left\{1,2,3\right\}\),\(A_2=\left\{S,M,L\right\}\),\(C=\left\{1,-1\right\}\).根据\(P_{\lambda} \left(X^{(j)}=a_{jl}|Y=c_k\right) = \frac{\sum \limits_{i=1}^N I(x_i^{(j)} = a_{jl}, y_i =c_k) + \lambda}{\sum \limits_{i=1}^N I(y_i=c_k) + S_j\lambda}\)和\(P_{\lambda}(Y=c_k) = \frac{\sum \limits_{i=1}I(y_i=c_k)+\lambda}{N+K\lambda}\)计算下列概率：
对于给定的\(x＝(2,S)^T\)，计算
由于\(P(Y=-1)P(X^{(1)}=2|Y=-1)P(X^{(2)}=S|Y=-1)\)最大，所以\(y=-1\)。
3.3 半朴素贝叶斯分类器
为了降低贝叶斯公式(7.8)中估计后验概率\(P(c|x)\)的困难，朴素贝叶斯分类器采用了属性条件独立性假设，但在现实任务中这假设很难成立。于是，人们尝试对属性条件独立性假设进行一定程度的放松，由此产生了一类称为“半朴素贝叶斯分类器”(sermi-naive Bayes classifier)的学习方法。
半朴素贝叶斯分类器的基本想法是适当考虑一部分属性间的相互依赖信息，从而既不需进行完全联合概率计算，又不至于彻底忽略了比较强的属性依赖关系。“独依赖估计”(One-Dependent Estimator，简称ODE)是半朴素贝叶斯分类器最常用的一种策略。顾名思义，所谓“独依赖”就是假设每个属性在类别之外最多仅依赖于一个其他属性，即
式(7.21)的解释
在朴素贝叶斯分类器中，假设每个属性独立地对分类结果发生影响，因此
\(P(x|c) = \prod \limits_{i=1}^d P(x_i|c)\)
“独依赖估计”策略假设每个属性在类别之外最多仅依赖一个其他属性，因此
\(P(x|c) = \prod \limits_{i=1}^d P(x_i|c,pa_i)\)其中\(pa_i\)为属性\(x_i\)所依赖的属性。前者只依赖类别\(c\)，而后者同时依赖于\(c\)和\(pa_i\)。
根据贝叶斯公式，类似于式(7.14)到(7.15)，舍掉分母常数项\(P(x)\)，即可得式(7.21)
\(P(c|x) \propto P(c)P(x|c) = P(c) \prod \limits_{i=1}^d P(x_i|c,pa_i)\)
在上一节中详细展示了\(P(x_i|c)\)，即先挑出类别为\(c\)的样本，若是离散属性则按计数法估计\(P(x_i|c)\)，若是连续属性则求这些样本的均值和方差，按高斯分布估计\(p(x_i|c)\)。现在估计\(P(x_i|c, pa_i)\)，则先挑出类别为\(c\)、属性\(x_i\)所依赖的属性值为\(pa_i\)的样本，剩下步骤与估算\(P(x_i|c)\)。
其中\(pa_i\)为属性\(x_i\)所依赖的属性，称为\(x_i\)的父属性。此时，对每个属性\(x_i\)，若其父属性\(pa_i\)已知，则可采用类似式(7.20)的办法来估计概率值\(P(x_i|c, pa_i)\)，于是，问题的关键就转化为如何确定每个属性的父属性，不同的做法产生不同的独依赖分类器。
最直接的做法是假设所有属性都依赖于同一个属性，称为“超父”(super-parent)，然后通过交叉验证等模型选择方法来确定超父属性，由此形成了SPODE(Super-Parent ODE)方法。
例如，在图7.1(b)中，\(x_1\)是超父属性。
图7.1的解释
图中箭头指向表示依赖关系，如所有属性\(x_1,x_2,...,x_d\)均依赖于类别变量\(y\)，因此存在由\(y\)指向\(x_1,x_2,...,x_d\)的箭头。图(a)中的NB各属性之间不存在依赖关系，因此无箭头；图(b)假设所有属性均依赖于\(x_1\)(\(x_1\)是超父属性)，因此存在由\(x_1\)指向\(x_2,x_3,...,x_d\)的箭头；图(c)中属性之间依赖关系是树形结构(不存在回路，则每个结点有且仅有一个父结点)。
TAN(Tree Augmented naive Bayes)[Friedman et al., 1997]则是在最大带权生成树(maximum weighted spanning tree)算法的基础上，通过以下步骤将属性间依赖关系约简为如图7.1(c)所示的树形结构：
- (1) 计算任意两个属性之间的条件互信息(conditional mutual information)\[I(x_i, x_j|y) = \sum \limits_{x_i,x_j; c \in \mathcal{Y}} P(x_i, x_j|c)log \frac{P(x_i,x_j|c)}{P(x_i|c)P(x_j|c)}; \quad \quad (7.22) \]
式(7.22)的解释
该式写为如下形式可能更容易理解:
\(I(x_i,x_j|y) = \sum \limits_{n=1}^N P(x_i,x_j|c_n)log \frac{P(x_i, x_j|c_n)}{P(x_i|c_n)P(x_j|c_n)}\)
其中\(i,j=1,2,...,d\)且\(i \neq j\)，\(N\)为类别个数。该式共可得到\(d(d-1) / 2\)个\(I(x_i, x_j|y)\)，即每对\((x_i, x_j)\)均有一个条件互信息\(I(x_i, x_j|y)\)
-
(2) 以属性为节点构建完全图，任意两个节点之间边的权重设为\(I(x_i,x_j|y)\);
-
(3) 构建此完全图的最大带权生成树，挑选根变量，将边置为有向；
-
(4) 加入类别节点\(y\)，增加从\(y\)到每个属性的有向边。
TAN算法的解释
该算法共分为四步，第一步即式(7.22)，每对\((x_i, x_j)\)得到一个条件互信息\(I(x_i, x_j|y)\)；
第二步以\(I(x_i, x_j|y)\)为权重，构建完全图，即每个属性为节点，节点之间边的权重由第一步计算所得；
第三步条用最大带权生成树算法，去除第二步完全图中的部分边，生成一颗树，并且该树所有边的权重之和应该是所有可能生成树中最大的，即最大带权的生成树，此时的树仍是无向图，挑选根变量使之变为有向图(从根节点开始，由父节点指向孩子节点，至于如何挑选根变量，书中没说，可以将树中所有节点均作为根节点，最后将学习结果做一次集成即可)；
第四步加入节点\(y\)，即得类似于图7.1(c)所示TAN。
容易看出，条件互信息\(I(x_i,x_j|y)\)刻画了属性\(x_i\)和\(x_j\)在已知类别情况下的相关性，因此，通过最大生成树算法，TAN实际上仅保留了强相关属性之间的依赖性。
AODE(Averaged One-Dependent Estimator)是一种基于集成学习机制、更为强大的独依赖分类器。与SPODE通过模型选择确定超父属性不同，AODE尝试将每个属性作为超父来构建SPODE，然后将那些具有足够训练数据支撑的SPODE集成起来作为最终结果，即
式(7.23)的推导
贝叶斯定理式(7.8)将联合概率\(P(x,c)\)换为等价形式\(P(x|c)P(c)\)；实际上，将向量\(x\)拆开，把\(P(x,c)\)写为\(P(x_1, x_2,...,x_d,c)\)形式，此时
\(P(x,c) = P(x_1, x_2,...,x_d, c) \\ = P(x_1, x_2,...,x_d|c)P(c) \\ = P(x_1,...,x_{i-1},x_{i+1},...,x_d|c,x_i)P(c, x_i)\)
其中第二个等号即\(P(x,c) = P(x|c)P(c)\)，而第三个等号是第二个等号的推广；更通俗地来说，\(P(A,B) = P(A|B)P(B)\)，对于第二个等号，\(A = x, B = c\)，而对于第三个等号，\(A = (x_1,...,x_{i-1}, x_{i+1},...,x_d)\)，\(B = (c, x_i)\)。
类似于式(7.14)的属性条件独立性假设，则
\(P(x_1, ...,x_{i-1}, x_{i+1},...,x_d|c, x_i) = \prod \limits_{j=1 and j \neq i}^d P(x_j|c, x_i)\)。
实际上，根据式(7.25)，若不考虑平滑项，\(P(x_j|c, x_i) = 1\)，因此在连乘中不起作用，可去除上式中的\(j \neq i\)约束(当\(j =i\)时，式(7.25)的\(|D_{c,x_i,x_j}|\)与\(|D_{c, x_i}|\)相等)，即\(P(x_1,...,x_{i-1}, x_{i+1},...,x_d|c,x_i) = \prod_{j=1}^d P(x_j|c, x_i)\)
综上所述：
\(P(c|x) = \frac{P(x,c)}{P(x)} = \frac{P(c, x_i)P(x_1,...,x_{i-1}, x_{i+1},x_d)|c,x_i}{P(x)} \\ \propto P(c, x_i)P(x_1,...,x_{i-1}, x_{i+1},...,x_d|c, x_i) \\ =P(c, x_i)\prod_{j=1}^d P(x_j|c, x_i)\)上式是将属性\(x_i\)作为超父属性的，AODE尝试将每个属性作为超父来构建SPODE，然后将那些具有足够虚拟那数据支撑的SPODE集成起来作为最终结果。具体来说，对于总共\(d\)个属性来说，共有\(d\)个不同的上式，集成直接求和即可，因为对于不同的类别标记\(c\)均有\(d\)个不同的上式；至于“足够训练数据支撑的SPODE”条件，注意(7.24)和(7.25)使用到了\(|D_{c,x_i,x_j}|\)与\(|D_{c, x_i}|\)，若\(D_{x_i}\)集合中样本数量过少，则\(|D_{c,x_i}|\)与\(|D_{c,x_i,x_j}|\)将会更少，因此在式(7.23)中要求\(D_{x_i}\)集合中样本数量不少于\(m'\)。
其中\(D_{x_i}\)是在第\(i\)个属性上取值为\(x_i\)的样本的集合，\(m'\)为阈值常数。
显然，AODE需估计\(P(c,x_i)\)和\(P(x_j|c,x_i)\)。类似式(7.20)，有
其中\(N_i\)是第\(i\)个属性可能的取值数，\(D_{c,x_i}\)是类别为\(c\)且在第\(i\)个属性上取值为\(x_i\)的样本集合，\(D_{c, x_i,x_j}\)是类别为\(c\)且在第\(i\)和第\(j\)个属性上取值分别为\(x_i\)和\(x_j\)的样本集合。
式(7.24)和式(7.25)的解释
这两个式子本身就是频率计数，它们都应用了类似于式(7.19)和式(7.20)的拉普拉斯修正。其中，式(7.24)分母之所以加\(N×N\)是因为在样本集合\(D(分母)\)中类别\(c\)与属性\(x_i(分子)\)的可能组合个数共有\(N×N_i\)种，式(7.25)分母之所以加\(N_j\)是因为样本集合\(D_{c,x_i}(分母)\)中属性\(x_j(分子)\)的可能取值个数共有\(N_j\)种。
例如，对西瓜数据集3.0有，
不难看出，与朴素贝叶斯分类器类似，AODE的训练过程也是“计数”，即在训练数据集上对复合条件的样本进行计数的过程。与朴素贝叶斯分类器相似，AODE无需模型选择，既能通过预计算节省预测时间，也能采取懒惰学习方式在预测时再进行计数，并且易于实现增量学习。
既然将属性条件独立性假设放松为独依赖假设可能获得泛化性能的提升，那么，能否通过考虑属性间的高阶依赖来进一步提升泛化性能呢？也就是说，将式(7.23)中的属性\(pa_i\)替换为包含\(k\)个属性的集合\(pa_i\)，从而将ODE拓展为kDE。需注意的是，随着\(k\)的增加，准确估计概率\(P(x_i|y, pa_i)\)所需的训练样本数量将以指数级增加。因此，若训练数据非常充分，泛化性能有可能提升；但在有限样本条件下，则又陷入估计高阶联合概率的泥沼。
3.4 贝叶斯网
贝叶斯网(Bayesian network)亦称“信念网”(belief network)，它借助有向无环图(Directed Acyclic Graph，简称DAG)来刻画属性之间的依赖关系，并使用条件概率表(Conditional Probability Table，简称CPT)来描述属性的联合概率分布。
具体来说，一个贝叶斯网B由结构G和参数\(\Theta\)两部分构成，即\(B = \left \langle G,\Theta \right \rangle\)。网络结构\(G\)是一个有向无环图，其每个节点对应于一个属性，若两个属性有直接依赖关系，则它们由一条边连接起来；参数\(\Theta\)定量描述这种依赖关系，假设属性\(x_i\)在\(G\)中的父节点集为\(\pi_i\)，则\(\Theta\)包含了每个属性的条件概率表\(\theta_{x_i|\pi_i} = P_B(x_i|\pi_i)\)。
作为一个例子，下图给出了西瓜问题的一种贝叶斯网结构和属性"根蒂"的条件概率表。从图中网中网络结构可看出，“色泽”直接依赖于“好瓜”和“甜度”，而“根蒂”则直接依赖于“甜度”；进一步从条件概率表能得到“根蒂”对“甜度”量化依赖关系，如\(P(根蒂=硬挺|甜度=高) = 0.1\)等。
3.4.1 结构
贝叶斯网结构有效地表达了属性间的条件独立性。给定父节点集，贝叶斯网假设每个属性与它的非后裔属性独立，于是\(B = \left \langle G,\Theta \right \rangle\)将属性\(x_1, x_2,...,x_d\)的联合概率分布定义为
式(7.26)的解释
注意本式前面的一段话：“给定父节点集，贝叶斯网假设每个属性与它的非后裔属性独立”，该式正是基于此而来。本式给出的是贝叶斯网所有节点的联合概率分布，以图7.2为例，节点\(x_1\)无父节点，因此在式(7.26)中为\(P(x_1)\)；节点\(x_2\)亦无父节点，因此在式(7.26)中为\(P(x_2)\)；节点\(x_3\)父节点为\(x_1\)，因此在式(7.26)中为\(P(x_3|x_1)\)；节点\(x_4\)父节点为\(x_1,x_2\)，因此在式(7.26)中为\(P(x_4|x_1, x_2)\)；节点\(x_5\)父节点为\(x_2\)，因此在式(7.26)中为\(P(x_5|x_2)\)，故图(7.2)的联合概率定义为:
\(P(x_1, x_2, x_3, x_4, x_5) = P(x_1)P(x_2)P(x_3|x_1)P(x_4|x_1, x_2)P(x_5|x_2)\).
以上面图7.2为例，联合概率分布定义为
显然，\(x_3\)和\(x_4\)在给定\(x_1\)的取值时独立，\(x_4\)和\(x_5\)在给定\(x_2\)的取值时独立，分别简记为\(x_3\perp x_4|x_1\)和\(x_4 \perp x_5 |x_2\).
下图7.3显示出贝叶斯网中三个变量之间的典型依赖关系，其中前两种在式(7.26)中已有所体现。
在“同父”(common parent)结构中，给定父节点\(x_1\)的取值，则\(x_3\)与\(x_4\)条件独立。在“顺序”结构中，给定\(x\)的值，则\(y\)与\(z\)条件独立。\(V\)型结构(V-structure)亦称“冲撞”结构，给定子节点\(x_4\)的取值，\(x_1\)与\(x_2\)必不独立；奇妙的是，若\(x_4\)的取值完全未知，则\(V\)型结构下\(x_1\)和\(x_2\)却是相互独立的。
这里我们可以做一个简单的验证:
式(7.27)的解释
该式不能基于概率论去推导，而是应该基于式(7.26)去推导。式(7.27)针对图7.3中的\(V\)型结构，这也是一个贝叶斯网，因此根据式(7.26)，\(V\)型结构三个节点的联合概率为
\(P(x_1, x_2, x_4) = P(x_1)P(x_2)P(x_4|x_1, x_2)\)
这正是式(7.27)中第二 个等号的由来。第一个等号就是概率论中的边际化，即要消去联合概率中的某个变量，只需对该变量积分(或求和)即可。而\(\sum_{x_4}P(x_4|x_1, x_2) = 1\)，因此有第三个等号。
其中\(\sum_{x_4}P(x_4|x_1, x_2) = 1\)这个结论是显而易见的，因为在给定\(x_1, x_2\)下，遍历了\(x_4\)所有状态。为便于理解，我们基于表4.1西瓜数据集2.0举例，可假设\(x_1,x_2\)分别为属性色泽中的青绿和属性根蒂中的蜷缩，\(x_4\)为纹理属性的可能取值(清晰、稍糊、模糊)，显然：\(P(纹理=清晰|色泽=青绿,根蒂=蜷缩)+P(纹理=清晰|色泽=青绿,根蒂=蜷缩) \\ + P(纹理=清晰|色泽=青绿,根蒂=蜷缩) = 1\),各概率计算参见7.3节各离散属性概率的估计。
本式推导的关键在于为什么\(P(x_1, x_2, x_4) = P(x_1)P(x_2)P(x_4|x_1, x_2)\)，因为根据概率论有\(P(x_1, x_2, x_4) = P(x_1, x_2)P(x_4|x_1, x_2)\)。
这样的独立性称为"边际独立性"(marginal independence)，记为\(x_1 \perp \!\!\! \perp x_2\).
事实上，一个变量取值的确定与否，能对另两个变量间的独立性发生影响，这个现象并非\(V\)型结构所特有。例如在同父结构中，条件独立性\(x_3 \perp \!\!\! \perp x_4 |x_1\)成立，但若\(x_1\)的取值未知，则\(x_3\)和\(x_4\)就不独立，即\(x_3 \perp \!\!\! \perp x_4\)不成立；在循序结构中，\(y \perp z | x\)，但是\(y \perp \!\!\! \perp z\)不成立。
为了分析有向图中变量间的条件独立性，可使用“有向分离”(D-separation)。
我们先把有向图转变为一个无向图:
- 找出有向图中的所有\(V\)型结构，在\(V\)型结构的两个父节点之间再加上一条无向边；
- 将所有有向边改为无向边.
由此产生的无向图称为“道德图”(moral graph)，令父节点相连的过程称为“道德化”(moralization)。
基于道德图能直观、迅速地找到变量间的条件独立性。假定道德图中有变量\(x,y\)和变量集合\(z = \left\{z_i\right\}\)，若变量\(x\)和\(y\)能在图上被\(z\)分开，即从道德图中将变量集合\(z\)去除后，\(x\)和\(y\)分属两个连通分支，则称变量\(x\)和\(y\)被\(z\)有向分离，\(x \perp y |z\)成立。
例如，上图7.2所对应的道德图如图7.4所示，从图中能容易地找出所有的条件独立关系：\(x_3 \perp x_4 |x_1\)，\(x_4 \perp x_5|x_2\)，\(x_3 \perp x_2 | x_1\)，\(x_3 \perp x_5 |x_1\)，\(x_3 \perp x_5 |x_2\)等。
图7.4的解释
注意：摘除某个变量后，所有以其为端口的边全部删除。例如，摘除\(x_1\)后，与\(x_1\)相连的三条边全部删除，此时\(x_3\)与其余部分完全隔开。
由图7.2的贝叶斯网到图7.4的道德图，目标是为了分析有向图中变量间的条件独立性，使用的方法是“有向分离”。
英文D-separation和moralization在论文中经常见到，应该理解其含义。
3.4.2 学习
若网络结构已知，即属性间的依赖关系已知，则贝叶斯网的学习过程相对简单，只需通过对训练样本“计数”，估计出每个节点的条件概率表即可。但在现实应用中我们往往并不知晓网络结构，于是，贝叶斯网学习的首要任务就是根据训练数据集来找出结构最“恰当”的贝叶斯网。
‘’评分搜索“是求解这一问题的常用办法。具体来说，我们先定义一个评分函数(score function)，以此来评估贝叶斯网与训练数据的契合程度，然后基于这个评分函数来寻找结构最优的贝叶斯网。显然，评分函数引入了关于我们希望获得什么样的贝叶斯网的归纳偏好。
“评分搜索”的解释
简单来说，将所有可能的贝叶斯网络结构当作定义域，每种结构都用事先定义好的评分函数映射得到一个评分，即\(y = f(x)\)形式，\(x\)当自变量(即某种贝叶斯网络结构)，函数值\(y\)为该结构的评分，贝叶斯网络结构学习的过程相当于在定义域寻找函数最大值的过程。
接下来的式(7.28)到式(7.31)都在介绍评分函数，了解概念即可，不必细究。
常用评分函数通常基于信息论准则，此类准则将学习问题看作一个数据压缩任务，学习的目标是找到一个能以最短编码长度描述训练数据的模型，此时编码的长度包括了描述模型自身所需的字节长度和使用该模型描述数据所需的字节长度。
对贝叶斯网学习而言，模型就是一个贝叶斯网，同时，每个贝叶斯网描述了一个在训练数据上的概率分布，自有一套编码机制能使那些经常出现的样本更短的编码。于是，我们应选择那个综合编码长度(包括描述网络和编码数据)最短的贝叶斯网，这就是“最小描述长度”(Minimal Description Length，简MDL)准则。
给定训练集\(D = \left\{x_1, x_2, ...,x_m\right\}\)，贝叶斯网\(B = \left \langle G,\Theta \right \rangle\)在D上的评分函数可写为
其中，\(|B|\)是贝叶斯网的参数个数；\(f(\theta)\)表示每个参数\(\theta\)所需的字节数；而
\(LL(B|D)\)是贝叶斯网\(B\)的对数似然，显然，式(7.28)的第一项是计算编码贝叶斯网B所需的字节数，第二项是计算B所对应的概率分布\(P_B\)需多少字节来描述D。于是，学习任务就转化为一个优化任务，即寻找一个贝叶斯网B使评分函数\(s(B|D)\)最小。
若\(f(\theta) = 1\)，即每个参数用1字节描述，则得到\(AIC\)(Akaike Information Criterion)评分函数
若\(f(\theta) = \frac{1}{2}logm\)，即每个参数用\(\frac{1}{2}log m\)字节描述，则得到BIC(Bayesian Information Criterion)评分函数
显然，若\(f(\theta) = 0\)，即不计算对网络进行编码的长度，则评分函数退化为负对数似然，相应的，学习任务退化为极大似然估计。
不难发现，若贝叶斯网\(B = \left \langle G,\Theta \right \rangle\)的网络结构\(G\)固定，则评分函数\(s(B|D)\)的第一项为常数，此时，最小化\(s(B|D)\)等价于对参数\(\Theta\)的极大似然估计。由式(7.29)和(7.26)可知，参数\(\theta_{x_i|\pi_i}\)能直接在训练数据\(D\)上通过经验估计获得，即
其中\(\hat P_D(\cdot)\)是\(D\)上的经验分布。因此，为了最小化评分函数\(s(B|D)\)，只需对网络结构进行搜索，而候选结构的最优参数可直接在训练集上计算得到。
不幸的是，从所有可能的网络结构空间搜索最优贝叶斯网结构是一个NP难问题，难以快速求解。
有两种常用的策略能在有限时间内求得近似解：
-
第一种策略是贪心法，例如从某个网络结构出发，每次调整一条边(增加、删除或调整方向)，直到评分函数值不再降低为止；
-
第二种策略是通过给网络结构施加约束来削减搜索空间，例如将网络结构限定为树形结构等。
3.4.3 推断
贝叶斯网训练好之后就能用来回答"查询"(query)，即通过一些属性变量的观测值来推测其他属性变量的取值。例如在西瓜问题中，若我们观测到西瓜色泽青绿、敲声浊响、根蒂蜷缩，想知道它是否成熟、甜度如何。这样通过已知变量观测值来推测待查询变量的过程称为“推断”(inference)，已知变量观测值称为“证据”(evidence)。
最理想的是直接根据贝叶斯网定义的联合概率分布来精确计算后验概率，不幸的是，这样的”精确推断“已被证明是NP难的；换言之，当网络节点较多、连接稠密时，难以进行精确推断，此时需借助”近似推断“，通过降低精度要求，在有限时间内求得近似解。
3.4.3.1 吉步斯采样
在现实应用中，贝叶斯网的近似推断常使用吉布斯采样(Gibbs sampling)来完成，这是一种随机采样方法，我们来看看它是如何工作的。
令\(Q = \left\{Q_1, Q_2,...,Q_n \right\}\)表示待查询变量，\(E = \left\{E_1,E_2,...,E_k \right\}\)为征集变量，已知其取值为\(e = \left\{e_1, e_2, ...,e_k \right\}\)，目标是计算后验概率\(P(Q=q|E=e)\)，其中\(q = \left\{q_1,q_2,...,q_n \right\}\)是待查询变量的一组取值。以西瓜问题为例，待查询变量为\(Q = \left\{好瓜,甜度\right\}\)，证据变量为\(E = \left\{色泽,敲声,根蒂 \right\}\)且已知其取值为\(e = \left\{青绿,浊响,蜷缩\right\}\)，查询的目标值是\(q = \left\{是,高\right\}\)，即这是好瓜且甜度高的概率有多大。
如下图7.5所示，吉布斯采样算法先随机产生一个与征集\(E = e\)一致的样本\(q^0\)作为初始点，然后每步从当前样本出发产生下一个样本。具体来说，在第\(t\)次采样中，算法先假设\(q^t = q^{t-1}\)，然后对非证据变量逐个进行采样改变其取值，采样概率根据贝叶斯网\(B\)和其他变量的当前取值(即\(Z=z\))计算获得。
假定经过\(T\)次采样得到的与\(q\)一致的样本共有\(n_q\)个，则可近似估算出后验概率
实质上，吉布斯采样是在贝叶斯网所有变量的联合状态空间与证据\(E = e\)一致的子空间中进行”随机漫步“(random walk)。每一步仅依赖于前一步的状态，这是一个”马尔可夫链“(Markov chain)。在一定条件下，无论从什么初始状态开始，马尔可夫链第\(t\)步的状态分布在\(t \to \infty\)时必收敛于一个平稳分布(stationary distribution)；对于吉布斯采样来说，这个分布恰好是\(P(Q|E = e)\)。因此，在\(T\)很大时，吉布斯采样相当于根据\(P(Q|E = e)\)采样，从而保证了式(7.33)收敛于\(P(Q=q|E=e)\)。
吉布斯采样算法描述如下:
需注意的是，由于马尔可夫链通常需很长时间才能趋于平稳分布，因此吉布斯采样算法的收敛速度较慢。此外，若贝叶斯网中存在极端概率“0”或“1”，则不能保证马尔可夫链存在平稳分布，此时吉布斯采样会给出错误的估计结果。
吉布斯采样算法的解释
首先，去除图中的第11行到第13行，才是标准的吉布斯采样算法，该三行仅是利用吉布斯采样结果得到式(7.33)的后验概率而已。这里参考机器学习经典专著《Pattern Recognition and Machine Learning》(简称PRML)中的吉布斯采样算法：
\(\begin{align} &Gibbs Sampling \\ &1. Initialize \left\{z_i: i=1,...,M \right\} \\ &2. For \tau = 1,...,T: \\ & - Sample z_1^{(\tau+1)} \sim p(z_1|z_2^{(\tau)},z_3^{(\tau)},...,z_M^{(\tau)}) \\ & - Sample z_2^{(\tau+1)} \sim p(z_2|z_1^{(\tau)},z_3^{(\tau)},...,z_M^{(\tau)}) \\ & \vdots \\ & - Sample z_j^{(\tau+1)} \sim p(z_j|z_1^{(\tau+1)},z_{j-1}^{(\tau+1)},...,z_M^{(\tau)}) \\ & \vdots \\ & - Sample z_M^{(\tau+1)} \sim p(z_M|z_1^{(\tau + 1)},z_2^{(\tau+1)},...,z_M^{(\tau+1)}). \end{align}\)
可以看出，第二步的For循环即算法描述过程中第3行到第14行的For循环(去除第11行至第13行)，区别在于PRML中将上面的第4行到第10行的For循环逐行给出的。
在上述过程中，大写黑体字母表示变量，如\(E,Q,Z\)，相应的小写字母\(e,q,z\)表示变量的取值，例如西瓜数据集2.0中，属性纹理为变量，则共有三个可能取值(清晰、稍糊、模糊)。
第5行是去除变量\(Q_i\)外的其他变量(每次仅更新一个变量)，第6行则是相应的变量取值；第7行则是根据贝叶斯网\(B\)计算分布\(P_B(Q_i|Z=z)\)，也就是根据除去变量\(Q_i\)外的所有其他变量来计算条件概率分布\(P_B(Q_i|Z=z)\)(例如，假设此处\(Q_i\)表示属性纹理，则该概率表示其他变量取值计算条件概率分布，包含三个概率值，分别是纹理为清晰、稍糊、模糊时的条件概率)；第8行根据\(P_B(Q_i|Z=z)\)采样获取\(Q_i\)取值(假如根据清晰、稍糊、模糊三个概率选择其中一个，每个属性值被选到的概率对应其概率\(P_B(Q_i|Z=z)\))；第9行更新待查询变量中\(Q_i\)的取值(如果只有一个待查询变量\(Q_i\)，那么第4行到第10行的\(For\)循环实际也仅包含一次循环)。
第4行到第10行的\(For\)循环就是每次更新一个待查询变量，该\(For\)循环结束则所有待查询变量被更新一次；然后第3行到第14行的\(For\)循环结束则所有待查询变量被更新\(T\)轮。我们要估计后验概率\(P(Q=q|E=e)\)，而吉布斯采样算法对\(Q\)随机赋初值(第2行)，那么每执行一次第3行到第14行的\(For\)循环则\(Q\)就会被更新一次，迭代\(T\)次\(Q\)会被更新\(T\)次，在这\(T\)次循环当众，\(Q\)的取值可能会等于待查询的\(q\)(例如\(Q\)为是否好瓜、甜度如何，而\(q\)为坏瓜、不甜)，假设有\(n_q\)次等于待查询的\(q\)，则待解后验概率\(P(Q=q|E=e) = \frac{n_q}{T}\)。
吉布斯采样算法能有效的依据是，经过若干次第3行到第14行的\(For\)循环后，\(Q\)各种可能的取值的出现概率收敛于一个平稳分布，该分布恰是\(P(Q=q|E=e)\)。例如(坏瓜、不甜)的后验概率为0.2，则第3行到第14行的For循环执行100次当中，会有约20次得到的\(q\)取值为(坏瓜、不甜)。
吉布斯采样是马尔可夫链蒙特卡罗(Markov Chain Monte Carlo,MCMC)方法的一种实现，是MCMC方法的代表算法Metropolis-Hastings(简称MH)的一种特殊形式。
(四)EM算法
4.1 EM算法介绍
在前面的讨论中，我们一直假设训练样本所有属性变量的值都已被观察到，即训练样本是“完整”的。但在现实应用中往往会遇到“不完整”的训练样本，例如由于西瓜的根蒂已脱落，无法看出是“蜷缩”还是“硬挺”，则训练样本的“根蒂”属性变量值未知。在这种存在“未观测”变量的情形下，是否仍能对模型参数进行估计呢？
未观察变量的学名是“隐变量”(latent variable)。令\(X\)表示已观测变量集，\(Z\)表示隐变量集，\(\Theta\)表示模型参数。若欲对\(\Theta\)做极大似然估计，则应最大化对数似然
式(7.34)的解释
本式等号右侧\(ln P(X,Z|\Theta)\)就是类似式(7.10)的似然概率，不同的是，等号左侧\(LL(\Theta)|X,Z)\)与式(7.10)表达方式不同：为了表示对\(\Theta\)做最大似然估计是在已知\(X,Z\)条件下，使用了类似于条件概率的表达方式。这里用\(LL(\Theta)\)作为表示符号简记。
然而由于\(Z\)是隐变量，上式无法直接求解，此时我们可通过对\(Z\)计算期望，来最大化已观测数据的对数“边际似然”(marginal likelihood)
式(7.35)的解释
在式(7.34)的解释中提到，对\(\Theta\)做最大似然估计需要已知\(X,Z\)。但是\(Z\)是隐变量，也就是未知的，因此只需考虑已知\(X\)即可，即最大化似然\(lnP(X|\Theta)\)。
概率\(P(X|\Theta)\)可以对概率\(P(X,Z|\Theta)\)做边际化(marginalization)而得，类似于式(7.27)解释中的\(\sum_{x_4}P(x_4|x_1, x_2) = 1\)，此处有\(P(X|\Theta) = \sum_{Z}P(X,Z|\Theta)\)。
概率模型有时既含有观测变量(observable variable)，又含有隐变量或潜在变量(latent variable)。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计法估计模型参数。但是，当模型含有隐变量时，就不能简单地使用这些估计方法。
EM(Expectation-Maximization)算法是常用的估计参数隐变量的利器，它是一种迭代式的方法，其基本想法是：若参数\(\Theta\)已知，则可根据训练推断出最优隐变量\(Z\)的值(E步)；反之，若\(Z\)的值已知，则可方便地对参数\(\Theta\)做极大似然估计(E步)。
于是，以初始值\(\Theta^0\)为起点，对式\(LL(\Theta|X) = lnP(X|\Theta) = ln \sum_{Z}P(X,Z|\Theta)\)，可迭代执行以下步骤直至收敛：
- 基于\(\Theta^t\)推断隐变量\(Z\)的期望，记为\(Z^t\)；
- 基于已观测变量\(X\)和\(Z^t\)对参数\(\Theta\)做极大似然估计，记为\(\Theta^{t+1}\)。
这就是EM算法的原型。
进一步，若我们不是取\(Z\)的期望，而是基于\(\Theta^t\)计算隐变量\(Z\)的概率分布\(P(Z|X,\Theta^t)\)，则EM算法的两个步骤是：
- E步(Expectation)：求期望(以当前参数\(\Theta^t\)推断隐变量分布\(P(Z|X,\Theta^t)\)，并计算对数似然\(LL(\Theta|X,Z)\)关于\(Ｚ\)的期望)\[Q(\Theta|\Theta^t) = \mathbb{E}_{Z|X,\Theta^t}LL(\Theta|X,Z) \quad (7.36) \]
式(7.36)的解释
我们通过数学期望公式来区分离散型变量和连续型变量，此处以离散型来解释式(7.36)。
设隐变量\(Z\)可能的取值个数为Ｎ，分别记为\(Z_1,Z_2,...,Z_N\)。在已知\(X、\Theta = \Theta^t\)以及\(Z = Z_i\)时，可以计算出条件概率\(P(Z_i|X,\Theta^t)\)；而若已知\(X\)和\(Z＝Z_i\)，则可由式(7.34)得到对数似然\(LL(\Theta|X,Z_i) = lnP(X,Z_i|\Theta)\)，因此可得式(7.36):\(\mathbb{E}_{Z|X, \Theta^t}LL(\Theta|X,Z) = \sum \limits_{i=1}^N P(Z_i|X, \Theta^t)LL(\Theta|X, Z_i) \\ = \sum \limits_{i=1}^N P(Z_i|X, \Theta^t)ln P(X,Z_i|\Theta)\)
可以发现，最后的表达式中\(X,Z_i, \Theta^t\)均为已知量，尤其是\(X\)为已观测变量，\(Z_i\)为隐变量\(Z\)可能的取值，是客观已知的量，与参数本身无关；而\(\Theta^t\)是参数\(\Theta\)的当前取值，因此上式可以表达为\(Q(\Theta|\Theta^t)\)，表示已知\(\Theta^t\)条件下的表达式。
- M步(Maximization)：求极大(寻找参数最大化期望似然)，即\[\Theta^{t+1} = \underset{\Theta}{arg max}Q(\Theta|\Theta^t) \quad (7.37) \]
简要来说，EM算法使用两个步骤计算：
- 第一步是期望\((E)\)步，利用当前估计的参数值来计算对数似然的期望值；
- 第二步是最大化(M)步，寻找能使\(E\)步的似然期望最大化的参数值。然后，新得到的参数值重新被用于\(E\)步，直至收敛到局部最优解。
事实上，隐变量估计问题也可通过梯度下降等优化算法求解，但由于求和的项数将随着隐变量的数目以指数级上升，会给梯度计算带来麻烦；而\(EM\)算法则可看作一种非梯度优化方法。
一般地，用Ｙ表示观测随机变量的数据，\(Z\)表示隐随机变量的数据。\(Y\)和\(Z\)连接在一起称为完全数据(complete-data)，观测数据\(Y\)又称为不完全数据(incomplete-data)。假设给定观测数据\(Y\)，其概率分布是\(P(Y|\theta)\)，其中\(\theta\)是需要估计的模型参数，那么不完全数据\(Y\)的似然函数是\(P(Y|\theta)\)，对数似然函数\(L(\theta)=logP(Y|\theta)\);
假设Y和Ｚ的联合概率分布是\(P(Y,Z|\theta)\)，那么完全数据的对数似然函数是\(logP(Y,Z|\theta)\).
EM算法通过迭代求\(L(\theta)=logP(Y|\theta)\)的极大似然估计，每次迭代包含两步:E步(求期望)；M步(求极大化)。
算法3 (EM算法)
输入：观测变量数据\(Y\)，隐变量数据\(Z\)，联合分布\(P(Y,Z|\theta)\)，条件分布\(P(Z|Y,\theta)\);
输出：模型参数\(\theta\)。
(1) 选择参数的初值\(\theta^{(0)}\)，开始迭代；
(2) \(E\)步：记\(\theta^{(i)}\)为第\(i\)次迭代参数\(\theta\)的估计值，在第\(i+1\)次迭代的\(E\)步，计算
这里，\(P(Z|Y,\theta^{(i)})\)是在给定观测数据\(Y\)和当前的参数估计\(\theta^{(i)}\)下隐变量数据\(Z\)的条件概率分布；
(3) \(M\)步：求使\(Q(\theta,\theta^{(i)})\)极大化的\(\theta\)，确定第\(i+1\)次迭代的参数的估计值\(\theta^{(i+1)}\)。
(4)重复第(2)步和第(3)步，直到收敛。
其中，\(Q(\theta,\theta^{(i)}) = E_z[logP(Y,Z|\theta)|Y,\theta^{(i)}]
= \sum \limits_{Z}logP(Y,Z|\theta)P(Z|Y, \theta^{(i)})\)中的函数\(Q(\theta,\theta^{(i)})\)是EM算法的核心，称为Q函数(Q function)。
\(Q\)函数的定义
完全数据的对数似然函数\(logP(Y,Z|\theta)\)关于在给定观测数据\(Y\)和当前参数\(\theta^{(i)}\)下对未观测数据\(Z\)的条件概率分布\(P(Z|Y,\theta^{(i)})\)的期望称为\(Q\)函数，即
下面关于EM算法作几点说明：
步骤(1) 参数的初值可以任意选择，但需注意EM算法对初值是敏感的。
步骤(2) E步求\(Q(\theta,\theta^{(i)})\)。Ｑ函数式中Ｚ是未观测数据，Ｙ是观测数据。
注意，Q(\(\theta,\theta^{(i)}\))的第1个变元表示要极大化的参数，第2个变元表示参数的当前估计值。每次迭代实际在求Q函数及其极大。
步骤(3) M步求\(Q(\theta,\theta^{(i)})\)的极大化，得到\(\theta^{(i+1)}\)，完成一次迭代\(\theta^{(i)}\to\theta^{(i+1)}\)。每次迭代使似然函数增大或达到局部极值。
步骤(4) 给出停止迭代的条件，一般是对较小的正数\(\varepsilon_1,\varepsilon_2\)，若满足
则停止迭代。
4.2 实例理解 抛硬币模型
4.2.1 引入
假设现在有两枚硬币1和2，,随机抛掷后正面朝上概率分别为P1，P2。为了估计这两个概率，做实验，每次取一枚硬币，连掷5下，记录下结果，如下：
| 硬币 | 结果 | 统计 |
|---|---|---|
| 1 | 正正反正反 | 3正-2反 |
| 2 | 反反正正反 | 2正-3反 |
| 1 | 正反反反反 | 1正-4反 |
| 2 | 正反反正正 | 3正-2反 |
| 1 | 反正正反反 | 2正-3反 |
可以很容易地估计出\(P1\)和\(P2\)，如下：
到这里，一切似乎很美好，下面我们加大难度。
4.2.2 加入隐变量\(z\)
还是上面的问题，现在我们抹去每轮投掷时使用的硬币标记，如下：
| 硬币 | 结果 | 统计 |
|---|---|---|
| Unknown | 正正反正反 | 3正-2反 |
| Unknown | 反反正正反 | 2正-3反 |
| Unknown | 正反反反反 | 1正-4反 |
| Unknown | 正反反正正 | 3正-2反 |
| Unknown | 反正正反反 | 2正-3反 |
好了，现在我们的目标没变，还是估计\(P1\)和\(P2\)，要怎么做呢？
显然，此时我们多了一个隐变量\(z\)，可以把它认为是一个5维的向量\((z1,z2,z3,z4,z5)\)，代表每次投掷时所使用的硬币，比如\(z_1\)，就代表第一轮投掷时使用的硬币是1还是2。但是，这个变量\(z\)不知道，就无法去估计\(P1\)和\(P2\)，所以，我们必须先估计出z，然后才能进一步估计\(P1\)和\(P2\)。
但要估计\(z\)，我们又得知道\(P1\)和\(P2\)，这样我们才能用最大似然概率法则去估计z，这不是鸡生蛋和蛋生鸡的问题吗，如何破？
答案就是先随机初始化一个\(P1\)和\(P2\)，用它来估计\(z\)，然后基于\(z\)，还是按照最大似然概率法则去估计新的P1和\(P2\)，如果新的\(P1\)和\(P2\)和我们初始化的\(P1\)和\(P2\)一样，请问这说明了什么？
这说明我们初始化的\(P1\)和\(P2\)是一个相当靠谱的估计！
就是说，我们初始化的\(P1\)和\(P2\)，按照最大似然概率就可以估计出\(z\)，然后基于\(z\)，按照最大似然概率可以反过来估计出\(P1\)和\(P2\)，当与我们初始化的\(P1\)和\(P2\)一样时，说明是\(P1\)和\(P2\)很有可能就是真实的值。这里面包含了两个交互的最大似然估计。
如果新估计出来的\(P1\)和\(P2\)和我们初始化的值差别很大，怎么办呢？就是继续用新的\(P1\)和\(P2\)迭代，直至收敛。这就是下面的EM初级版。
4.2.3 EM初级版
我们不妨这样，先随便给P1和P2赋一个值，比如：
\(P1 = 0.2\)，\(P2 = 0.7\)。
然后，我们看看第一轮抛掷最可能是哪个硬币。
如果是硬币1，得出3正2反的概率为 \(0.2*0.2*0.2*0.8*0.8 = 0.00512\).
如果是硬币2，得出3正2反的概率为\(0.7*0.7*0.7*0.3*0.3=0.03087\).
然后依次求出其他4轮中的相应概率。做成表格如下：
| 轮数 | 若是硬币1 | 若是硬币2 |
|---|---|---|
| 1 | 0.00512 | 0.03087 |
| 2 | 0.02048 | 0.01323 |
| 3 | 0.08192 | 0.00567 |
| 4 | 0.00512 | 0.03087 |
| 5 | 0.02048 | 0.01323 |
按照最大似然法则：
第1轮中最有可能的是硬币2
第2轮中最有可能的是硬币1
第3轮中最有可能的是硬币1
第4轮中最有可能的是硬币2
第5轮中最有可能的是硬币1
我们就把上面的值作为z的估计值。然后按照最大似然概率法则来估计新的P1和P2。
设想我们是全知的神，知道每轮抛掷时的硬币就是如本文第001部分标示的那样，那么，P1和P2的最大似然估计就是0.4和0.5（下文中将这两个值称为P1和P2的真实值）。那么对比下我们初始化的P1和P2和新估计出的P1和P2：
| 初始化的P_1 | 估计出的P_1 | 真实的P_1 | 初始化的P_2 | 估计出的P_2 | 真实的P_2 |
|---|---|---|---|---|---|
| 0.2 | 0.33 | 0.4 | 0.7 | 0.6 | 0.5 |
看到没？我们估计的\(P1\)和\(P2\)相比于它们的初始值，更接近它们的真实值了！
可以期待，我们继续按照上面的思路，用估计出的\(P1\)和\(P2\)再来估计\(z\)，再用\(z\)来估计新的\(P1\)和\(P2\)，反复迭代下去，就可以最终得到\(P1 = 0.4，P2=0.5\)，此时无论怎样迭代，\(P1\)和\(P2\)的值都会保持0.4和0.5不变，于是乎，我们就找到了P1和P2的最大似然估计。
这里有两个问题：
1、新估计出的\(P1\)和\(P2\)一定会更接近真实的\(P1\)和\(P2\)？
答案是：没错，一定会更接近真实的P1和P2，数学可以证明，但这超出了这里的主题，请参阅其他书籍或文章。
2、迭代一定会收敛到真实的\(P1\)和\(P2\)吗？
答案是：不一定，取决于\(P1\)和\(P2\)的初始化值，上面我们之所以能收敛到\(P1\)和\(P2\)，是因为我们幸运地找到了好的初始化值。
4.2.4 EM进阶版
下面，我们思考下，上面的方法还有没有改进的余地？
我们是用最大似然概率法则估计出的z值，然后再用\(z\)值按照最大似然概率法则估计新的\(P1\)和\(P2\)。也就是说，我们使用了一个最可能的\(z\)值，而不是所有可能的\(z\)值。
如果考虑所有可能的\(z\)值，对每一个z值都估计出一个新的\(P1\)和\(P2\)，将每一个\(z\)值概率大小作为权重，将所有新的\(P1\)和\(P2\)分别加权相加，这样的\(P1\)和\(P2\)应该会更好一些。
所有的\(z\)值有多少个呢？显然，有2^5=32种，需要我们进行32次估值？
不需要，我们可以用期望来简化运算。
| 轮数 | 若是硬币1 | 若是硬币2 |
|---|---|---|
| 1 | 0.00512 | 0.03087 |
| 2 | 0.02048 | 0.01323 |
| 3 | 0.08192 | 0.00567 |
| 4 | 0.00512 | 0.03087 |
| 5 | 0.02048 | 0.01323 |
利用上面这个表，我们可以算出每轮抛掷中使用硬币1或者使用硬币2的概率。
比如第1轮，使用硬币1的概率是： \(P_1(硬币=硬币1) = \frac{0.00512}{(0.00512+0.03087)}=0.14\).
使用硬币2的概率是\(P_2(硬币=硬币2) = 1-0.14 = 0.86\).
依次可以算出其他4轮的概率，如下：
上表中的右两列表示期望值。看第一行，0.86表示，从期望的角度看，这轮抛掷使用硬币2的概率是0.86。相比于前面的方法，我们按照最大似然概率，直接将第1轮估计为用的硬币2，此时的我们更加谨慎，我们只说，有0.14的概率是硬币1，有0.86的概率是硬币2，不再是非此即彼。这样我们在估计P1或者P2时，就可以用上全部的数据，而不是部分的数据，显然这样会更好一些。
这一步，我们实际上是估计出了z的概率分布，这步被称作E步。
结合下表：
| 硬币 | 结果 | 统计 |
|---|---|---|
| Unknown | 正正反正反 | 3正-2反 |
| Unknown | 反反正正反 | 2正-3反 |
| Unknown | 正反反反反 | 1正-4反 |
| Unknown | 正反反正正 | 3正-2反 |
| Unknown | 反正正反反 | 2正-3反 |
我们按照期望最大似然概率的法则来估计新的\(P1\)和\(P2\)：
以\(P1\)估计为例，第1轮的3正2反相当于
\(0.14*3=0.42\)正
\(0.14*2=0.28\)反
依次算出其他四轮，列表如下：
| 轮数 | 正面 | 反面 |
|---|---|---|
| 1 | 0.42 | 0.28 |
| 2 | 1.22 | 1.83 |
| 3 | 0.94 | 3.76 |
| 4 | 0.42 | 0.28 |
| 5 | 1.22 | 1.83 |
| 总计 | 4.22 | 7.98 |
可以看到，改变了\(z\)值的估计方法后，新估计出的\(P1\)要更加接近0.4。原因就是我们使用了所有抛掷的数据，而不是之前只使用了部分的数据。
这步中，我们根据E步中求出的z的概率分布，依据最大似然概率法则去估计\(P1\)和\(P2\)，被称作M步。
例题7 三硬币模型
假设有3枚硬币，分别记作\(A,B,C\)。这些硬币正面出现的概率分别是\(\pi,p\)和\(q\)。进行如下掷硬币实验：先掷银币Ａ,根据其结果选出硬币Ｂ或者硬币Ｃ，正面选择硬币B，反面选硬币Ｃ；然后掷选出的硬币，掷硬币的结果出现正面记作1，出现反面记作0；独立地重复n次试验(这里，n=10)，观测结果如下：
假设只能观测到掷硬币的结果，不能观测掷硬币的过程。问如何估计三硬币正面出现的概率，即三硬币模型的参数。
解 三硬币模型可以写作
考虑求模型参数\(\theta＝(\pi,p,q)\)的极大似然估计，即
这个问题没有解析解，只有通过迭代的方法求解。
EM算法就是可以用于求解这个问题的迭代算法。
下面给出针对以上问题的EM算法，其推导过程省略。
EM算法首先选取参数的初值，记作\(\theta^{(0)}=(\pi^{(0)},p^{(0)},q^{(0)})\)，然后通过下面的步骤迭代计算参数的估计值，直至收敛为止。第\(i\)次迭代参数的估计值为\(\theta^{(i)}=(\pi^{(0)},p^{(0)},q^{(0)})\)。
EM算法的第\(i+1\)迭代如下：
E步：计算在模型参数\(\pi^{(i},p^{(i)},q^{(i)}\)下观测数据\(y_j\)来自掷硬币B的概率
M步：计算模型参数的新估计值
进行数值计算。假设模型参数的初值取为
由\(\mu_j^{(i+1)} = \frac{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}}
{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j} + (1-\pi^{(i)})(q^{(i)})^{y_j}(1-q^{(i)})^{1-y_j}}\),对\(y_j=1\)与\(y_j=0\)均有\(\mu_j^{(1)}=0.5\).
利用迭代公式，可得
由于\(\mu_j^{(i+1)}\)更新迭代公式，可得
继续迭代，得
于是得到模型参数\(\theta\)的极大似然估计：
\(\pi＝0.5\)表示硬币A是均匀的，这一结果容易理解。
如果取初值\(\pi^{(0)}=0.4,p^{(0)}=0.6,q^{(0)}=0.7\),那么得到的模型参数的极大似然估计是\(\hat \pi=0.4064,\hat p = 0.5368,\hat q = 0.6432\)。这就是说，EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。
4.2 EM算法在无监督学习中的应用
监督学习是由训练数据\(\left\{(x_1,y_1),(x_2,y_2),...,(x_N,y_N)\right\}\)学习条件概率分布\(P(Y|X)\)或决策函数\(Y=f(x)\)作为模型，用于分类、回归、标注等任务。这时训练数据中的每个样本点由输入和输出对组成。
有时训练数据只有输入没有对应的输出\(\left\{(x_1,\cdot),(x_2,\cdot),...,(x_N,\cdot)\right\}\)，从这样的数据学习模型称为无监督模型学习问题。EM算法可以用于生成模型的无监督学习。生成模型由联合概率分布\(P(X,Y)\)表示，可以认为无监督学习训练数据是联合概率分布产生的数据。\(X\)是观测数据，\(Y\)为未观测数据。
参考材料：
[1] 周志华《机器学习》
[2] 李航《统计学习方法》
[3] 邱锡鹏《神经网络与深度学习》
[4] 通俗理解GAN（一）：把GAN给你讲得明明白白
[5] 如何感性地理解EM算法?
