《统计学习方法》Notes(1)

前言

本篇博客基于《统计学习方法》前半部分监督学习的相关内容，主要记录一些书中的重点内容以及自己的一点总结，以及一些方法的代码实现举例。博客中所有相关代码在https://github.com/Aron00123/Statistical-Learning-Methods-Codes.git

感知机

感知机的学习算法

• 原始形式

输入：训练数据集，其中；学习率​。

输出：w，b；感知机模型。

（1）选取初值

（2）在训练集中选取数据

（3）若，

$$w\leftarrow w+\eta y_ix_i$$ $$b\leftarrow b+ \eta y_i$$

（4）转至（2），直至训练机中没有误分类点。

• 对偶形式

输入：训练数据集，其中；学习率​。

输出：；感知机模型。其中。

（1）

（2）在训练集中选取数据

（3）如果

$$\alpha_i\leftarrow \alpha_i + \eta$$ $$b\leftarrow b+\eta y_i$$

（4）转至（2）直到没有误分类数据。

K邻近算法

算法内容

• k近邻法

输入：训练数据集，其中​；

输出：实例x所属的类y。

（1）根据给定的距离度量，在训练集T中找出与x最邻近的k个点，涵盖这k个点的x的邻域记作；

（2）在中根据分类决策规则决定x的类别y

$$y=\arg \max_{c_j} \sum_{x_i\in N_k(x)} I(y_j=c_j),\quad i=1,2,...,N;j=1,2,...,K$$

其中I为指示函数，当时函数值为1，否则为0。

kd树

• 构造平衡kd树算法

输入：k维空间数据集，其中；

输出：kd树；

（1）构造根结点，根结点对应于包含T的k维空间的超矩形区域。选择为坐标轴，以T中所有实例的坐标的中位数为切分点，将根结点对应的超矩形区域切分为两个子区域。切分由通过切分点并与坐标轴垂直的超平面实现。由根结点生成深度为1的左右子结点：左子结点对应坐标小于切分点的子区域，右子结点对应于坐标大于切分点的子区域。将落在超平面上的实例点保存在根结点。

（2）对深度为j的点，选择为切分的坐标轴，，重复（1）的步骤。

（3）直到两个字区域没有实例存在时停止，从而形成kd树划分。

• 用kd树的最近邻搜索

输入：已构造的kd树，目标点x；

输出：x的最近邻。

（1）在kd树中找出包含目标点x的叶结点：从根结点出发，递归地向下访问kd树。若目标点x当前维的坐标小于切分点的坐标，则移动到左子结点，否则移动到右子结点。直到子结点为叶结点为止；

（2）以此叶结点为“当前最近点”；

（3）递归地向上回退，在每个结点进行以下操作：

（a）若该结点保存的实例点比当前最近点距离目标点更近，则以该实例点为“当前最近点”；

（b）当前最近点一定存在于该结点的一个子结点对应的区域。检查该子结点的父结点的另一子结点对应的区域是否有更近的点。具体地，检查另一子结点对应的区域是否与以目标点为球心，以目标点与“但前最近点”的距离为半径的超球体相交。若相交，则移动到另一个子结点递归地进行最近邻搜索；若不相交，则向上回退；

（4）回退至根结点时，搜索结束。最后的“当前最近点”即为x的最近邻点。

朴素贝叶斯法

​朴素贝叶斯算法是一种基于贝叶斯定理的分类算法。贝叶斯定理是一种用于计算条件概率的方法，它描述了在已知某些条件下，某事件发生的概率。在分类问题中，我们希望根据输入数据的特征将其分为不同的类别。朴素贝叶斯算法的“朴素”之处在于它假设输入数据的特征之间是相互独立的，这意味着算法会假设每个特征对于确定类别的影响是独立的。尽管这个假设在实际情况中并不总是成立，但朴素贝叶斯算法的简单性和高效性使它成为许多分类问题的一个常用选择。

• 朴素贝叶斯算法

输入：训练数据，其中有共n个特征，（第j个特征可以取个值），；实例x；

输出：实例x的分类；

（1）计算先验概率及条件概率

$$P(Y=c_k)=\frac{\overset{N}{\underset{i=1}{\sum}}I(y_i=c_k)}{N},\quad k=1,2,...,K$$ $$P(X^{(j)}=a_{jl}|Y=c_k)=\frac{\overset{N}{\underset{i=1}{\sum}}I(x_i^{(j)}=a_{jl},y_i=c_k)}{\overset{N}{\underset{i=1}{\sum}}I(y_i=c_k)},\quad j=1,2,...,n;l=1,2,...,S_j;k=1,2,...,K$$

（2）对于给定的实例，计算

$$P(Y=c_k)\prod_{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k),\quad k=1,2,...,K$$

（3）确定实例x的类

$$y=\arg \max_{c_k}P(Y=c_k)\prod_{j=1}^nP(X^{(j)}=x^{(j)}|Y=c_k)$$

决策树

决策树由节点和边组成。每个节点代表一个特征属性上的测试，每条边表示一个测试结果的输出，从根节点开始沿着路径向下，直到达到叶子节点，叶子节点包含最终的输出值。决策树学习通常包括3个步骤：特征选择、决策树的生成和决策树的修剪。

特征选择

• 信息增益算法

输入：训练数据集D和特征A；

输出：特征A对训练数据集D的信息增益；

（1）计算数据集D的经验熵

$$H(D)=-\sum_{k=1}^K\frac{|C_k|}{|D|}\log_2 \frac{|C_k|}{|D|}$$

（2）计算特征A对数据集D的经验条件熵

$$H(D|A)=\sum_{i=1}^n \frac{|D_i|}{|D|}H(D_i)=-\sum_{i=1}^n\frac{|D_i|}{|D|}\sum_{k=1}^K\frac{|D_{ik}|}{|D_i|}\log_2 \frac{|D_{ik}|}{|D_i|}$$

（3）计算信息增益

$$g(D,A)=H(D)-H(D|A)$$

通常根据信息增益准则选择最佳特征，即选择信息增益最大的特征。

• 信息增益比

以信息增益作为划分训练数据集的特征，存在偏向于选择取值较多的特征的问题。使用信息增益比可以很好的解决这个问题：

$$g_R(D,A)=\frac{g(D,A)}{H_A(D)}$$

其中，，n是特征A取值的个数。

生成决策树

• ID3算法

输入：训练数据集D，特征集A阈值；

输出：决策树T；

（1）若D中所有实例属于同一类，则T为单结点树，并将类作为该结点的类标记，返回T；

（2）若，则T为单结点树，并将D中实例数最大的类作为该结点的类标记，返回T；

（3）否则，计算A中各特征对D的信息增益，选择信息增益最大的特征 ：

（4）如果的信息增益小于阈值，则置T为单结点树，并将D中实例数最大的类作为该结点的类标记，返回T：

（5）否则，对的每一可能值，依将D分割为若干非空子集，将中实例数最大的类作为标记，构建子结点，由结点及其子结点构成树T，返回T;

（6）对第i个子结点，以为训练集，以为特征集，递归调用（1）~（5），得到子树 ,返回。

• C4.5算法

输入：训练数据集D，特征集A阈值；

输出：决策树T；

（1）若D中所有实例属于同一类，则T为单结点树，并将类作为该结点的类标记，返回T；

（2）若，则T为单结点树，并将D中实例数最大的类作为该结点的类标记，返回T；

（3）否则，计算A中各特征对D的信息增益比，选择信息增益比最大的特征 ：

（4）如果的信息增益比小于阈值，则置T为单结点树，并将D中实例数最大的类作为该结点的类标记，返回T：

（5）否则，对的每一可能值，依将D分割为若干非空子集，将中实例数最大的类作为标记，构建子结点，由结点及其子结点构成树T，返回T;

（6）对第i个子结点，以为训练集，以为特征集，递归调用（1）~（5），得到子树 ,返回。

决策树的剪枝

• 树的剪枝算法

输入：生成算法产生的整个树T；

输出：修剪后的子树；

（1）计算每个结点的经验熵；

（2）递归地从树的叶结点向上回缩。设一组叶结点回缩到其父结点之前的整体树的损失函数值大于之后的整体树的损失函数值，则进行剪枝，即将父结点变为新的叶结点；

（3）返回（2），直至不能继续为止，得到损失函数最小的子树。

设树T的业绩点个数为，t是树T的叶结点，该叶结点有个样本点，其中k类的样本点有个，则决策树学习的损失函数可以定义为

$$C_\alpha(T)=\sum_{t=1}^{|T|}N_tH_t(T)+\alpha|T|$$

式中经验熵为

$$H_t(T)=-\sum_k\frac{N_{tk}}{N_t}\log \frac{N_{tk}}{N_t}$$

可以将等式右侧第一项记为，表示模型对训练数据的预测误差，即模型与训练数据的拟合程度，表示模型复杂度，参数控制两者之间的影响。

CART算法

• 基尼指数

假设有K个类，样本点属于第k类的概率为，则概率分布的基尼指数定义为

$$Gini(p) = \sum_{k=1}^K p_k(1-p_k)=1-\sum_{k=1}^Kp_k^2$$

对于二类分类问题，若样本点属于第1个类的概率是p，则概率分布的基尼指数为

$$Gini(p)=2p(1-p)$$

对于给定的样本集合D，其基尼指数为

$$Gini(p)=1-\sum_{k=1}^K\left(\frac{|C_k|}{|D|}\right)^2$$

其中是D中属于第k累的样本子集，K是类的个数。

• CART生成算法

输入：训练数据集D，停止计算的条件（结点中的样本个数小于预定阈值或样本集的基尼指数小于预定阈值或没有更多的特征）；

输出：CART决策树；

（1）设结点的训练数据集为D，计算现有特征对该数据集的基尼指数。对每一个特征A，对其可能取的每个值，根据的结果是或否将D分割成和两部分，从而计算时的基尼指数；

（2）选取基尼指数最小的特征A及其对应的切分点作为最优特征与最优切分点，再依这两者生成两个子结点；

（3）对两个子结点递归调用（1）和（2），直至满足停止条件；

（4）生成CART决策树。

• CART剪枝算法

输入：CART算法生成的决策树；

输出：最优决策树；

（1）设；

（2）设；

（3）自下而上地对各内部结点t计算，以及

$$g(t)=\frac{C(t)-C(T_t)}{|T_t|-1}$$ $$\alpha=\min (\alpha,g(t))$$

其中表示以t为根结点的子树，是对训练数据的预测误差，是的叶结点个数。

（4）对的内部结点t进行剪枝，并对叶结点t以多数表决法决定其类，得到树T；

（5）设；

（6）如果不是由根结点及两个叶结点构成的树，则回到步骤（2）；否则令；

（7）采用交叉验证法在子树序列中选取最优子树。

支持向量机

支持向量机是一种二类分类模型可分为线性可分支持向量机、线性支持向量机以及非线性支持向量机。下面分别对这三类支持向量机以及核函数、SMO进行分析。

线性可分支持向量机

• 线性可分支持向量机的学习算法

输入：线性可分训练数据集，其中；

输出：最大间隔分离超平面和分类决策函数。

（1）构造并求解约束最优化问题：

$$\max_{w,b}\quad \frac{1}{2}||w||^2$$ $$s.t.\quad y_i(w\cdot x_i+b)-1\geq 0,\quad i=1,2,...,N$$

（2）由此得到分离超平面：

$$w^{\ast}\cdot x+b^{\ast}=0$$

分类决策函数

$$f(x)=sign(w^{\ast}\cdot x+b^{\ast})$$

• 学习的对偶算法

输入：线性可分训练数据集，其中；

输出：分离超平面和分类决策函数。

（1）构造并求解约束最优化问题：

$$\min_\alpha\quad \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^N\alpha_i\alpha_jy_iy_j(x_i\cdot x_j)-\sum_{i=1}^N\alpha_i$$ $$s.t.\quad \sum_{i=1}^N\alpha_iy_i=0$$ $$\alpha_i\geq 0,\quad i=1,2,...,N$$

求得最优解。

（2）计算

$$w^{\ast}=\sum_{i=1}^N\alpha_i^{\ast}y_ix_i$$

并选择的一个正分量，计算

$$b^{\ast}=y_j-\sum_{i=1}^N\alpha_i^{\ast}y_i(x_i\cdot x_j)$$

（3）求得分离超平面

$$w^{\ast}\cdot x + b^{\ast}=0$$

分类决策函数：

$$f(x)=sign(w^{\ast}\cdot x + b^{\ast})$$

线性支持向量机

• 线性支持向量机学习算法

输入：线性可分训练数据集，其中；

输出：分离超平面和分类决策函数。

（1）选择惩罚参数，构造并求解凸二次规划问题

$$\min_{\alpha} \quad \frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_j(x_i\cdot x_j)-\sum_{i=1}^N\alpha_i$$ $$s.t.\quad \sum_{i=1}^N\alpha_iy_i=0$$ $$0\leq \alpha_i\leq C,\quad i=1,2,...,N$$

求得最优解。

（2）计算

$$w^{\ast}=\sum_{i=1}^N\alpha_i^{\ast}y_ix_i$$

并选择的一个正分量，计算

$$b^{\ast}=y_j-\sum_{i=1}^N\alpha_i^{\ast}y_i(x_i\cdot x_j)$$

（3）求得分离超平面

$$w^{\ast}\cdot x + b^{\ast}=0$$

分类决策函数：

$$f(x)=sign(w^{\ast}\cdot x + b^{\ast})$$

非线性支持向量机与核函数

• 核函数

设是输入空间，又设是特征空间，如果存在一个从到的映射

$$\phi(x):\mathcal{X}\to\mathcal{H}$$

使得对所有的，函数满足条件

$$K(x,z)=\phi(x)\cdot \phi(z)$$

则称为核函数，为映射函数，式中为和的内积。

一些常用的核函数如下：

（1）多项式和函数

$$K(x,z)=(x\cdot z+1)^p$$

（2）高斯核函数

$$K(x,z)=\exp (-\frac{||x-z||^2}{2\sigma^2})$$

（3）字符串核函数

设是长度大于等于n的字符串的集合，s是的元素。下面建立字符串集合到特征空间的映射。表示定义在上的实数空间，其每一维对应一个字符串，映射将字符串s对应于空间的一个向量，其在u维上的取值为

$$[\phi_n(s)]_u=\sum_{i:s(i)=u}\lambda^{l(i)}$$

其中是一个衰减参数，表示字符串i的长度，求和在s中所有与u相同的子串上进行。两个字符串s和t上的字符串核函数是基于映射的特征空间中的内积：

$$k_n(s,t)=\sum_{u\in\sum^n}[\phi_n(s)]_u[\phi_n(t)]_u=\sum_{u\in\sum^n}\sum_{(i,j):s(i)=t(j)=u}\lambda^{l(i)}\lambda^{l(j)}$$

直观上来看，两个字符串相同的子串越多，它们就越相似，字符串核函数的值就越大。字符串核函数可以由动态规划快速计算

关于核函数的通俗理解可以看https://zhuanlan.zhihu.com/p/47541349

• 非线性支持向量机

从非线性分类训练机，通过核函数与软间隔最大化，或凸二次规划，学习得到的分类决策函数

$$f(x)=sign(\sum_{i=1}^N\alpha_i^{\ast}y_iK(x,x_i)+b^{\ast})$$

称为非线性支持向量，是正定核函数。

• 非线性支持向量机学习算法

输入：线性可分训练数据集，其中；

输出：分类决策函数。

（1）选取适当的核函数和适当的参数C，构造并求解最优化问题

$$\min\quad \frac{1}{2}\sum_{i=1}^N\sum_{i=1}^N\alpha_i\alpha_jy_iy_jK(x_i,x_j)-\sum_{i=1}^N\alpha_i$$ $$s.t.\quad \sum_{i=1}^N\alpha_iy_i=0$$ $$0\leq \alpha_i\leq C,\quad i=1,2,...,N$$

求得最优解。

（2）选择的一个正分量，计算

$$b^{\ast}=y_j-\sum_{i=1}^N\alpha_i^{\ast}y_iK(x_i, x_j)$$

（3）构造决策函数：

$$f(x)=sign(\sum_{i=1}^N\alpha_i^{\ast}y_iK(x, x_i)+b^{\ast})$$

当K是正定核函数时，算法中的问题是凸二次规划问题，解释存在的。

序列最小最优化算法

SMO算法要解决一下凸二次规划的对偶问题：

$$\min_\alpha\quad \frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_jK(x_i,x_j)-\sum_{i=1}^N\alpha_i$$ $$s.t.\quad \sum_{i=1}^N\alpha_iy_i=0$$ $$0\leq\alpha_i\leq C,\quad i=1,2,...,N$$

变量是拉格朗日乘子，一个变量对应于一个样本点；变量的总数等于训练样本容量N。SMO算法是一种启发式算法：如果所有变量的解都满足此最优化问题的KKT条件，那么这个最优化问题的解就得到了。否则，选择两个变量，固定其他变量，正对这两个变量构建一个二次规划问题。这两个变量一个是违反KKT条件最严重的，另一个由约束条件自动确定，即​，如果两个变量中一个确定了，另一个也就随之确定了。如此，SMO算法将原问题不断分解为子问题并对子问题求解，进而达到求解原问题的目的。

• SMO算法内容

输入：线性可分训练数据集，其中，精度；

输出：近似解

（1）取初值，令；

（2）选取优化变量，解析求解两个变量的最优化问题，求得最优解，更新为；

（3）若在精度范围内满足停机条件

$$\sum_{i=1}^N\alpha_iy_i=0$$ $$0\leq \alpha_i\leq C,\quad i=1,2,...,N$$ $$y_i\cdot g(x_i)= \begin{cases} \geq1,\quad \{x_i|\alpha_i=0 \} \\ =1,\quad \{x_i|0 < \alpha_i < C \} \\ \leq1,\quad \{x_i|\alpha_i=C \} \end{cases}$$

其中，

$$g(x_i)=\sum_{j=1}^N\alpha_jy_jK(x_j,x_i)+b$$

则转（4）；否则令，转（2）；

（4）取​。

scikit-learn库中的SVM

• svm.SVC

用于分类问题的支持向量机。适用于非线性分类任务。

常用参数：

C: 正则化参数，默认为1.0

kernel: 核函数类型，如 ‘linear’, ‘poly’, ‘rbf’, ‘sigmoid’

degree: 多项式核函数的次数（’poly’时有效），默认为3

gamma: 核函数系数，’scale’ (默认) 或 ‘auto’

• svm.LinearSVC

适用于线性可分的分类任务，使用线性核。

常用参数：

C: 正则化参数，默认为1.0

• svm.NuSVC

与svm.SVC类似，但使用不同的正则化参数nu，控制支持向量的比例。

常用参数：

nu: 支持向量的比例（0到1之间），默认为0.5

kernel: 核函数类型

• svm.SVR

用于回归问题的支持向量机。

常用参数：

C: 正则化参数

kernel: 核函数类型

degree: 多项式核函数的次数（’poly’时有效）

gamma: 核函数系数

。。。

提升方法

AdaBoost算法

假设给定一个二类分类的训练数据集，其中，

输入：训练数据集；

输出：最终分类器G(x)；

（1）初始化训练数据的权值分布；

（2）对

（a）使用具有权值分布的训练数据集学习，得到基本分类器；

（b）计算在训练数据集上的分类误差率；

（c）计算的系数；

（d）更新训练数据记得权值分布

$$D_{m+1}=(w_{m+1,1},...,w_{m+1,i},...,w_{m+1,N})$$ $$w_{m+1,i}=\frac{w_{mi}}{Z_m}\exp{(-\alpha_my_iG_m(x_i))},\quad i=1,2,...,N$$

其中，是规范化因子，使得构成一个概率分布。

（3）构建基本分类器的线性组合，得到最终分类器。

Adaboost算法可以简述为三个步骤：
（1）首先，是初始化训练数据的权值分布D1。假设有N个训练样本数据，则每一个训练样本最开始时，都被赋予相同的权值：w1=1/N。
（2）然后，训练弱分类器。具体训练过程中是：如果某个训练样本点，被弱分类器准确地分类，那么在构造下一个训练集中，它对应的权值要减小；相反，如果某个训练样本点被错误分类，那么它的权值就应该增大。权值更新过的样本集被用于训练下一个分类器，整个训练过程如此迭代地进行下去。
（3）最后，将各个训练得到的弱分类器组合成一个强分类器。各个弱分类器的训练过程结束后，加大分类误差率小的弱分类器的权重，使其在最终的分类函数中起着较大的决定作用，而降低分类误差率大的弱分类器的权重，使其在最终的分类函数中起着较小的决定作用。
换而言之，误差率低的弱分类器在最终分类器中占的权重较大，否则较小。

训练误差分析

• AdaBoost的训练误差界

AdaBoost算法最终分类器的训练误差界为

$$\frac{1}{N}\sum_{i=1}^NI(G(x_i)\neq y_i)\leq \frac{1}{N}\sum_i\exp (-y_if(x_i))=\prod_mZ_m$$

Proof: When , we can get ,

then

$$\exp(-y_if_i(x))\geq 1$$

So the left part can be done. The right part needs

$$w_{mi}\exp (-\alpha_my_iG_m(x_i))=Z_mw_{m+1,i}$$

The entire proof then proceeds as follows

$$\frac{1}{N}\sum_i\exp (-y_if(x_i)) = \frac{1}{N} \sum_i \exp (-\sum_{m=1}^{M}\alpha_my_iG_m(x_i))\quad\quad\quad\quad$$ $$\quad\quad\quad\quad\quad=\sum_iw_{1i}\prod_{m=1}^{M}\exp (-\alpha_my_iG_m(x_i))$$ $$\quad\quad\quad\quad\quad\quad=Z_i\sum_iw_{2i}\prod_{m=1}^M\exp (-\alpha_my_iG_m(x_i))$$ $$\quad\quad\quad\quad\quad\quad\quad=Z_1Z_2\sum_iw_{3i}\prod_{m=3}^{M}\exp (-\alpha_my_iG_m(x_i))$$ $$=...\quad\quad\quad\quad\quad\quad\qquad$$ $$\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad=Z_1Z_2\cdot \cdot \cdot Z_{M-1}\sum_iw_{Mi}\exp (-\alpha_My_iG_M(x_i))$$ $$=\prod_{m=1}^MZ_m\quad\quad\quad\quad\quad\quad$$

• 二类分类问题AdaBoost的训练误差界

$$\prod_{m=1}^{M}Z_m=\prod_{m=1}^{M}[2\sqrt{e_m(1-e_m)}]=\prod_{m=1}^M\sqrt{(1-4\gamma_m^2)}\leq\exp (-2\sum_{m=1}^{M}\gamma_m^2)$$

这里。

Proof: According to the definition of , we can get

$$Z_m=\sum_{i=1}^Nw_{mi}\exp (-\alpha_my_iG_m(x_i))$$ $$\qquad\qquad\quad=\sum_{y_i=G_m(x_i)}w_{mi}e^{-\alpha_m}+\sum_{y_i\neq G_m(x_i)}w_{mi}e^{\alpha_m}$$ $$=(1-e_m)e^{-\alpha_m}+e_me^{\alpha_m}\quad$$ $$\qquad\quad=2\sqrt{e_m(1-e_m)}=\sqrt{(1-4\gamma_m^2)}$$

And when it comes to the enquality

$$\prod_{m=1}^M\sqrt{(1-4\gamma_m^2)}\leq\exp (-2\sum_{m=1}^M\gamma_m^2)$$

it can be prooved by the Taylor Expand of and at :

$$\sqrt{(1-4\gamma_m^2)}\leq \exp (-2\gamma_m^2)$$

EM算法

算法的引入

EM算法是一种迭代算法，用于在统计模型中估计概率模型的参数，尤其适用于模型中包含隐变量的情况。关于EM算法的解释我找到了这样一篇回答
https://www.zhihu.com/question/27976634/answer/3449882660 可以作为参考。

输入：观测变量数据Y，隐变量数据Z，联合分布，条件分布。

输出：模型参数。

（1）选择参数的初值，开始迭代；

（2）E步：记为第i次迭代参数的估计值，在第i+1次迭代的E步，计算$Q(\theta, \theta^{(i)})=E_z[\log P(Y,Z|\theta)|Y,\theta^{(i)}]=\underset{z}\sum\log P(Y,Z|\theta)P(Z|Y,\theta^{(i)})$，这里是在给定观测数据Y和当前的参数估计下隐变量数据Z的条件概率分布；

（3）M步：求使极大化的，确定第i+1次迭代的参数的估计值；

（4）重复第（2）步和第（3）步，直到收敛。

算法中的Q函数尤为重要，是EM算法的核心。Q函数是完全数据的对数似然函数关于在给定观测数据Y和当前参数下对未观测数据Z的条件概率分布的期望。

下面是一个简单的EM算法例子，用来估计投硬币正面朝上的概率：

'''
EM/EM_simple.py
'''
import numpy as np

def expectation(data, bias):
    heads_prob = bias
    tails_prob = 1 - bias
    return data * heads_prob / (data * heads_prob + (1 - data) * tails_prob)

def maximization(data, responsibilities):
    return np.sum(responsibilities) / len(data)

def EM(data, initial_bias, num_iterations):
    bias = initial_bias

    for i in range(num_iterations):
        responsibilities = expectation(data, bias)
        bias = maximization(data, responsibilities)

    return bias

# Generate synthetic data
np.random.seed(0)
data = np.random.randint(0, 2, size=100)
# Initial guess for bias
initial_bias = 0.5
# Run EM algorithm
estimated_bias = EM(data, initial_bias, num_iterations=100)

print("Estimated bias:", estimated_bias)

在高斯混合模型学习中的应用

高斯混合模型是指具有如下形式的概率分布模型：。其中，是系数（总和为1），是高斯密度分布（正态分布），。假设观测数据由上述高斯混合模型生成，接下来使用EM算法估计高斯混合模型的参数​。

• 首先可以明确对于每个观测数据，K个0-1隐变量表示此数据是否来自第k个分模型。则完全数据可以表示为​。从而可以写出完全数据的似然函数：

$$P(y,\gamma|\theta)=\overset{N}{\underset{j=1}{\prod}}P(y_j,\gamma_{j1},\gamma_{j2},...,\gamma_{jK}|\theta)=\overset{K}{\underset{k=1}{\prod}}\overset{N}{\underset{j=1}{\prod}}[\alpha_k\phi(y_j|\theta_k)]^{\gamma_{jk}}\\ =\overset{K}{\underset{k=1}{\prod}}\alpha_k^{n_k}\overset{N}{\underset{j=1}{\prod}}[\phi(y_j|\theta_k)]^{\gamma_{jk}}=\overset{K}{\underset{k=1}{\prod}}\alpha_k^{n_k}\overset{N}{\underset{j=1}{\prod}}[\frac{1}{\sqrt{2\pi}\sigma_k}\exp{(-\frac{(y_i-\mu_k)^2)}{2\sigma_k^2}}]^{\gamma_{jk}}$$

式中​。从而完全数据的对数似然函数为：

$$\log P(y,\gamma|\theta)=\overset{K}{\underset{k=1}{\sum}}\{n_k\log \alpha_k+\overset{N}{\underset{j=1}{\sum}}\gamma_{jk}[\log(\frac{1}{\sqrt{2\pi}})-\log\sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2] \}$$

• 由上，可以确定Q函数

$$Q(\theta,\theta^{(i)})=E[\log P(\mathbf{y} \mid \gamma, \theta^{(i)})] =E\left \{ \overset{K}{\underset{k=1}{\sum}} n_k \log \alpha_k + \overset{N}{\underset{j=1}{\sum}} \gamma_{jk} \left[ \log \left( \frac{1}{\sqrt{2\pi}} \right) - \log \sigma_k - \frac{1}{2\sigma_k^2} (y_j - \mu_k)^2 \right]\right \}\\ =\overset{K}{\underset{k=1}{\sum}} \left \{ \overset{N}{\underset{j=1}{\sum}} (E\gamma_{jk}) \log \alpha_k + \overset{N}{\underset{j=1}{\sum}} (E\gamma_{jk}) \left[ \log \left( \frac{1}{\sqrt{2\pi}} \right) - \log \sigma_k - \frac{1}{2\sigma_k^2} (y_j - \mu_k)^2 \right] \right \}$$

接着需要计算，记为。

$$\hat{\gamma}_{jk}=E(\gamma_{jk}|y,\theta)=P(\gamma_{jk}=1|y,\theta)=\frac{P(\gamma_{jk}=1,y_j|\theta)}{\overset{K}{\underset{k=1}{\sum}}P(\gamma_{jk}=1,y_j|\theta)}=\frac{P(y_j \mid \gamma_{jk} = 1, \theta) P(\gamma_{jk} = 1 \mid \theta)}{\overset{K}{\underset{k=1}{\sum}} P(y_j \mid \gamma_{jk} = 1, \theta) P(\gamma_{jk} = 1 \mid \theta)}\\ =\frac{\alpha_k\phi(y_j|\theta_k)}{\overset{K}{\underset{k=1}{\sum}}\alpha_k\phi(y_j|\theta_k)},j=1,2,...,N;k=1,2,...,K$$

是当前模型参数下第j个观测数据来自第k个分模型的概率，称为分模型k对观测数据的相应度。将和​代回上式，有

$$Q(\theta, \theta^{(i)}) = \overset{K}{\underset{k=1}{\sum}} \left\{ n_k \log \alpha_k + \overset{N}{\underset{j=1}{\sum}} E[\hat{\gamma}_{jk}] \left[ \log \left( \frac{1}{\sqrt{2\pi}} \right) - \log \sigma_k - \frac{1}{2\sigma_k^2} (y_j - \mu_k)^2 \right]\right\}$$

• 接着确定EM算法的M步，用，和表示的各参数。求前两者只需对上式求偏导并令其为0即可，得到；后者在​条件下求偏导并令其为0得到：

$$\mu_k = \frac{\sum_{j=1}^N \hat{\gamma}_{jk} y_j}{\sum_{j=1}^N \hat{\gamma}_{jk}},\quad \hat{\sigma}_k^2 = \frac{\sum_{j=1}^N \hat{\gamma}_{jk} (y_j - \mu_k)^2}{\sum_{j=1}^N \hat{\gamma}_{jk}},\quad\hat{\alpha}_k = \frac{\sum_{j=1}^N \hat{\gamma}_{jk}}{N}, \quad k = 1, 2, \ldots, K$$

重复以上计算直到对数似然函数值无明显变化为止。

下面展示一个简单的高斯混合模型应用EM算法的代码实现：

'''
EM/EM.py
'''
import numpy as np
from scipy.stats import norm


def expectation(data, means, covariances, weights):
    num_components = len(means)
    num_points = len(data)
    responsibilities = np.zeros((num_points, num_components))

    for i in range(num_points):
        for j in range(num_components):
            responsibilities[i, j] = weights[j] * norm.pdf(data[i], loc=means[j], scale=np.sqrt(covariances[j]))

    responsibilities /= np.sum(responsibilities, axis=1, keepdims=True)
    return responsibilities

def maximization(data, responsibilities):
    num_components = responsibilities.shape[1]
    num_points = len(data)

    total_weight = np.sum(responsibilities, axis=0)
    weights = total_weight / num_points
    means = np.sum(data[:, np.newaxis] * responsibilities, axis=0) / total_weight

    covariances = np.zeros(num_components)
    for j in range(num_components):
        diff = data - means[j]
        covariances[j] = np.sum(responsibilities[:, j] * diff ** 2) / total_weight[j]

    return means, covariances, weights

def EM(data, initial_means, initial_covariances, initial_weights, num_iterations):
    means = initial_means
    covariances = initial_covariances
    weights = initial_weights

    for i in range(num_iterations):
        responsibilities = expectation(data, means, covariances, weights)
        means, covariances, weights = maximization(data, responsibilities)

    return means, covariances, weights

# Generate synthetic data
np.random.seed(0)
data = np.concatenate([np.random.normal(-10, 1, 1000), np.random.normal(10, 1, 1000)])

# Initial guesses for parameters
initial_means = [-8, 8]
initial_covariances = [1, 1]
initial_weights = [0.4, 0.6]

# Run EM algorithm
estimated_means, estimated_covariances, estimated_weights = EM(data, initial_means, initial_covariances,
                                                               initial_weights, num_iterations=100)

print("Estimated means:", estimated_means)
print("Estimated covariances:", estimated_covariances)
print("Estimated weights:", estimated_weights)

EM算法推广

EM算法还可以解释为F函数的极大-极大算法，基于此可以推广至GEM算法（广义期望极大算法）等。

F函数的极大-极大算法

假设隐变量数据Z的概率分布为，定义分布与参数的函数，此函数即为F函数。式中是分布的熵。关于F函数有几个重要的引理：

• 对于固定的，存在唯一的分布极大化，这时，并且随连续变化。
• 若，则。

由上，可以得到关于EM算法用F函数的极大-极大算法的解释：

设为观测数据的对数似然函数，为EM算法得到的参数估计序列，如果函数在和有局部最大值，那么也在有局部最大值。类似地，如果函数在和达到全局最大值，那么也在达到全局最大值。

EM算法的一次迭代可由F函数的极大-极大算法实现：

设为第i次迭参数的估计，为第i次迭代函数的估计。在第i+1次迭代的两步为：

（1）对固定的，求使极大化；

（2）对固定的，求极大化。

证明略。

GEM算法

• GEM算法1:

输入：观测数据，F函数

输出：模型参数。

（1）选择参数的初值，开始迭代；

（2）第i+1次迭代，第1步：求使极大化；

（3）第2步：求使极大化；

（4）重复第（2）步和第（3）步，直到收敛。

• GEM算法2:

输入：观测数据，F函数

输出：模型参数。

（1）选择参数的初值，开始迭代；

（2）第i+1次迭代，第1步：计算；

（3）第2步：求使；

（4）重复第（2）步和第（3）步，直到收敛。

• GEM算法3:

输入：观测数据，F函数

输出：模型参数。

（1）选择参数的初值，开始迭代；

（2）第i+1次迭代，第1步：计算；

（3）第2步：进行d次条件极大化：

在保持不变的条件下求使达到极大的；然后在的条件下求使达到极大的；如此继续，一共执行d次，得到使得；

（4）重复第（2）步和第（3）步，直到收敛。

隐马尔可夫模型

基本概念

隐马尔可夫模型是关于时序的概率模型，描述一个隐藏的马尔可夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测从而产生观测随机序列的过程。模型包含三个元素：状态转移概率矩阵、观测概率矩阵和初始状态概率向量。具象地说，记Q是所有可能的状态的集合，V是所有可能的观测的集合：，其中N是可能的状态数，M是可能的观测数。则

• 状态转移概率矩阵，其中是在时刻t处于状态的条件下在时刻t+1转移到状态的概率；
• 观测概率矩阵，其中是在时刻t处于状态的条件下生成观测的概率。
• 初始状态概率向量，其中是时刻t=1处于状态的概率。

状态转移概率矩阵A与初始状态概率向量​确定了隐藏的马尔可夫链，生成不可观测的状态序列。观测概率矩阵B确定了如何从状态生成观测，与状态序列综合确定了如何产生观测序列。

隐马尔可夫模型主要有3个基本问题：

• 概率计算问题。给定模型和观测序列，计算在模型下观测序列O出现的概率。
• 学习问题。已知观测序列，估计模型参数，使得在该模型下观测序列概率最大。即用极大似然估计的方法估计参数。
• 预测问题，即解码（decoding）问题。已知模型和观测序列，求对给定观测序列条件概率最大的状态序列。即给定观测序列，求最有可能的对应的状态序列。

概率计算算法

直接计算法

此方法通过概率公式直接计算出，但是由于计算量过大，所以此方法理论可行但是实际不可行。具体地，通过列举所有可能的长度为T的状态序列，求各个状态序列I与观测序列的联合概率，然后对所有可能的状态序列求和，得到：

状态序列的概率是：；

对固定的状态序列，观测序列的概率是：，从而O和I同时出现的联合概率为。然后对所有可能的状态序列I求和即可。最后的求和计算量很大，是的，从而计算上不可行。

forward-backward algorithm（前向-后向算法）

• 前向概率

给定隐马尔可夫模型，定义到时刻t部分观测序列为且状态为的概率为前向概率，记作

​。

• 前向算法

输入：隐马尔可夫模型，观测序列O；

输出：观测序列概率；

（1）初值；

（2）对，；

（3）​。

前向算法的计算复杂度为，实际上是可以接受的。

• 后向概率

给定隐马尔可夫模型，定义在时刻t状态为的条件下，从t+1到T的部分观测序列为的概率为后向概率，记作。

• 后向算法

输入：隐马尔可夫模型，观测序列O；

输出：观测序列概率。

（1）初值；

（2）对，；

（3）。

利用前向概率和后向概率可以得到：。

一些概率与期望值

• 给定模型与观测O，在时刻t处于状态的概率

$$\gamma_t(i)=P(i_t=q_i|O,\lambda)=\frac{P(i_t=q_i,O|\lambda)}{P(O|\lambda)}=\frac{\alpha_t(i)\beta_t(i)}{\overset{N}{\underset{j=1}{\sum}}\alpha_t(j)\beta_t(j)}$$

• 给定模型与观测O，在时刻t处于状态且在时刻t+1处于状态的概率

$$\xi_t(i, j) = \frac{P(i = q_i, i+1 = q_j \mid O, \lambda)}{P(O \mid \lambda)} = \frac{P(i = q_i, i+1 = q_j \mid O, \lambda)}{\overset{N}{\underset{i=1}{\sum}} \overset{N}{\underset{j=1}{\sum}} P(i = q_i, i+1 = q_j \mid O, \lambda)} =\frac{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{\overset{N}{\underset{i=1}{\sum}} \overset{N}{\underset{j=1}{\sum}} P(i = q_i, i+1 = q_j \mid O, \lambda)}$$

• • 在观测O下状态i出现的期望值：；
  • 在观测O下状态i转移的期望值：；
  • 在观测O下由状态i转移到状态j的期望值：。

学习算法

隐马尔可夫模型的学习根据训练数据是包括观测序列和对应的状态序列还是只有观测序列可分为监督学习和无监督学习。

监督学习方法

假设已给训练数据包含S个长度相同的观测序列和对应的状态序列，那么可以利用极大似然估计法来估计隐马尔可夫模型的参数。

• 转移概率​的估计：；

• 观测概率的估计：；

• 初始状态概率的估计为S个样本中初始状态为的频率。

Baum-Welch算法（无监督学习）

输入：观测数据；

输出：隐马尔可夫模型参数。

（1）初始化。对n=0，选取，得到模型；

（2）递推，对，，。右端各值按观测和模型​计算；

（3）终止。得到模型参数。

预测算法

近似算法

大致地说，近似算法的思想是在每个时刻t选择在该时刻最有可能出现的状态，从而得到一个状态序列，将它作为预测的结果。给定隐马尔可夫模型和观测序列，在时刻t处于状态的概率，在每一时刻t最有可能的状态是，从而得到状态序列。近似算法的优点是计算简单、直观，缺点是不能保证预测的状态序列整体是最有可能的序列。

维特比算法

维特比算法是一种动态规划算法。根据动态规划的原理，最优路径如果在时刻t通过节点，那么这一路径从结点到终点的部分路径必然为最优的。从而我们只需从时刻t=1开始，递推地计算在时刻t状态为i的各条部分路径的最大概率，直至得到时刻t=T状态为i的各条路径的最大概率，此时最大概率即为最优路径的概率，最优路径的终点也可以同时得到。接着可以由后向前逐步得到最优路径的所有结点，从而求出具体路径。

导入两个变量和，定义在时刻t状态为i的所有单个路径中概率最大值为

。由定义可得变量的递推公式：。定义在时刻t状态为i的所有单个路径中概率最大的路径的第t-1个结点为。

输入：模型和观测；

输出：最优路径。

（1）初始化；

（2）递推，对t=2,…,T，；

（3）终止，。

（4）最优路径回溯，对t=T-1,T-2,…,1，。

条件随机场

​条件随机场（CRF）是给定一组输入随机变量条件下另一组输出随机变量的条件概率分布模型，特点是假设输出随机变量构成马尔可夫随机场。下文只涉及线性链条件随机场。关于更加形象地理解CRF，我找到了这样一篇回答：
https://www.zhihu.com/question/35866596/answer/139485548 可以参考一下。还有一篇比较啰嗦的
https://www.zhihu.com/question/35866596/answer/236886066 也可以参考。

概率无向图模型

概率图模型是由图表示的概率分布。设有联合概率分布P(Y)，由无向图G=(V,E)表示，结点表示随机变量，边表示随机变量之间的概率依赖关系。如果联合概率分布P(Y)满足成对、局部或全局马尔可夫性，就称此联合概率分布为概率无向图模型，或马尔可夫随机场。下面对成对、局部和全局马尔可夫性作出解释：

• 成对马尔可夫性

设u和v是无向图G中任意两个没有边连接的结点，结点u和v分别对应随机变量和。其他所有结点为O，对应的随机变量组是。成对马尔可夫性指给定随机变量组的条件下随机变量和是条件独立的，即。

• 局部马尔可夫性

设是无向图G中任意一个结点，W是与v有边连接的所有结点，O是v和W以外的所有其他结点。则给在给定随机变量组的条件下随机变量与随机变量组是独立的，即，特别的在时，。

• 全局马尔可夫性

设结点集合A，B是在无向图G中被结点集合C分开的任意结点集合，“分开”指删除结点集合C后A与B不连通。则给定随机变量组条件下随机变量组和是条件独立的，即。

概率无向图模型的因子分解

将概率无向图模型的联合概率分布表示为其最大团上的随机变量的函数的乘积形式的操作，成为概率无向图模型的因子分解，即。其中C是图G上的最大团，Z是规范化因子。规范化因子保证构成一个概率分布。称为势函数，通常定义为。

条件随机场的定义与形式

设X与Y是随机变量，是在给定X的条件下Y的条件概率分布。若随机变量Y构成一个由无向图表示的马尔可夫随机场，即对成立，则称条件概率分布为条件随机场。表示图G中结点v的所有邻居，表示v以外的所有结点。下面主要考虑线性链的情况：设，均为线性链表示的随机变量序列，若在给定随机变量序列X的条件下，随机变量序列Y构成条件随机场，即满足马尔可夫性（在i=1，i=n时只考虑单边），则称为线性链条件随机场。特别地，在标注问题中，X表示输入观测序列，Y表示对应的输出标记序列或状态序列。

条件随机场的参数化形式

设​为线性链条件随机场，则在随机变量X取值为x的条件下，随机变量Y取值为y的条件概率

$$P(y|x)=\frac{1}{Z(X)}\exp{(\underset{i,k}{\sum}\lambda_kt_k(y_{i-1},y_i,x,i)+\underset{i,l}{\sum}\mu_ls_l(y_i,x,i))}$$

其中。式中，和是特征函数，和是对应的权值。是规范化因子，求和是在所有可能的输出序列上进行的。是定义在边上的特征函数，称为转移特征，依赖于当前和前一个位置；是定义在结点上的特征函数，称为状态特征，依赖于当前位置。和都依赖于位置，是局部特征函数，且通常取值为0或1（当满足条件是取值为1，否则取值为0）。条件随机场完全由特征函数，和对应的权值，​确定。

条件随机场的简化形式

设有个转移特征，个状态特征，，记

$$f_k(y_{i-1},y_i,x,i)= \begin{cases} t_k(y_{i-1},y_i,x,i) \\ s_l(y_i,x,i) \end{cases}$$

然后，对转移特征与状态特征在各个位置i求和，记作。用表示特赠你的权值，即

$$w_k= \begin{cases} \lambda_k,\quad k=1,2,...,K_1 \\ \mu_l,\quad k=K_1+l;l=1,2,...,K_2 \end{cases}$$

从而

$$P(y|x)=\frac{1}{Z(x)}\exp{\overset{k}{\underset{k=1}{\sum}}w_kf_k(y,x)},\quad Z(x)=\underset{y}{\sum}\exp{\overset{k}{\underset{k=1}{\sum}}w_kf_k(y,x)}$$

若令，则

$$P_w(y|x)=\frac{\exp{(w\cdot F(y,x))}}{Z_w(x)},\quad Z_w(x)=\underset{y}{\sum}\exp{(w\cdot F(y,x))}$$

条件随机场的矩阵形式

假设是一个线性链条件随机场，表示对给定观测序列x，相应的标记序列y的条件概率。对每个标记序列引进特殊的起点和终点状态标记和，这时标注序列的概率可以通过矩阵形式表示并计算。对观测序列x的每一个位置，定义一个m阶矩阵（m是标记的取值的个数）

$$M_i(x)=[M_i(y_{i-1},y_i|x] \\ M_i(y_{i-1}y_i|x)=\exp{(W_i(y_{i-1},y|x))} \\ W_i(y_{i-1},y_i|x)=\overset{K}{\underset{i=1}{\sum}}w_kf_k(y_{i-1},y_i,x,i)$$

这样，给定观测序列x，标记序列y的非规范化概率可以通过n+1个矩阵的乘积表示。从而

$$P_w(y|x)=\frac{1}{Z_w(x)}\overset{n+1}{\underset{i=1}{\prod}}M_i(y_{i-1},y_i|x)$$

其中，为规范化因子，是n+1个矩阵的乘积的（start,stop）元素：

$$Z_w(x)=(M_1(x)M_2(x)\cdot\cdot\cdot M_{n+1}(x))_{start,stop}$$

概率计算问题

前向-后向算法

• 前向向量

对每个指标，定义前向向量：

$$\alpha_o(y|x)= \begin{cases} 1,\quad y=start \\ 0,\quad 否则 \end{cases}$$

递推公式为：

$$\alpha_i^T(y_i|x)=\alpha_{i-1}^T(y_{i-1}|x)M_i(y_{i-1},y_i|x),\quad i=1,2,...,n+1$$

又可表示为。易知是m维列向量。

• 后向向量

类似地，定义后向向量：

$$\beta_{n+1}(y_{n+1}|x)= \begin{cases} 1,\quad y_{n+1}=stop \\ 0,\quad 否则 \end{cases}$$

递推公式为：

$$\beta_i(y_i|x)=M_i(y_i,y_{i+1}|x)\beta_{i-1}(y_{i+1}|x)$$

又可表示为。

由前向-后向向量定义不难得到：，式中是元素均为1的m维列向量。

概率计算以及期望值计算

由上述定义，已得到如下两个条件概率的表达式：

$$P(Y_i=y_i|x)=\frac{\alpha_i^T(y_i|x)\beta_i(y_i|x)}{Z(x)}$$ $$P(Y_{i-1}=y_{i-1},Y_i=y_i|x)=\frac{\alpha_{i-1}^T(y_{i-1}|x)M_i(y_{i-1},y_i|x)\beta_i(y_i|x)}{Z(x)}$$

其中，​。

同样利用前向-后向向量，可以计算特征函数关于条件分布和联合分布的数学期望。

• 关于条件分布​的数学期望是：

$$E_{P(Y|X)}[f_k] = \sum_y P(y|x)f_k(y, x) = \sum_{i=1}^{n+1} \sum_{y_{i-1},y_i} f_k(y_{i-1}, y_i, x, i) \frac{\alpha_{i-1}^T(y_{i-1} | x) M_i(y_{i-1}, y_i | x) \beta_i(y_i | x)}{Z(x)}$$

• 关于联合分布的数学期望是（假设经验分布为）：

$$E_{P(X, Y)}[f_k] = \sum_{x,y} P(x, y) \sum_{i=1}^{n+1} f_k(y_{i-1}, y_i, x, i) = \sum_{x} \tilde{P}(x) \sum_{y} P(y \mid x) \sum_{i=1}^{n+1} f_k(y_{i-1}, y_i, x, i) = \sum_{x} \tilde{P}(x) \sum_{i=1}^{n+1} \sum_{y_{i-1}, y_i} f_k(y_{i-1}, y_i, x, i) \frac{\alpha_{i-1}^T(y_{i-1} \mid x) M_i(y_{i-1}, y_i \mid x) \beta_i(y_i \mid x)}{Z(x)}$$

其中，。

学习算法

• 改进的迭代尺度法

输入：特征函数；经验分布；

输出：参数估计值；模型；

（1）对所有，取初值​；

（2）对每一，

（a）当时，令是如下方程的解：

$$\sum_{x,y} \tilde{P}(x) P(y \mid x) \sum_{i=1}^{n+1} t_k(y_{i-1}, y_i, x, i) \exp(\delta_k T(x, y)) = E_{\tilde{p}}[t_k]$$

当时，令是如下方程的解：

$$\sum_{x,y} \tilde{P}(x) P(y \mid x) \sum_{i=1}^{n} s_l(y_{i-1}, y_i, x, i) \exp(\delta_{K_1+l} T(x, y)) = E_{\tilde{p}}[s_l]$$

其中，

$$T(x, y) = \sum_k f_k(y, x) = \sum_{k=1}^{K} \sum_{i=1}^{n+1} f_k(y_{i-1}, y_i, x, i)$$

（b）更新值：​；

（3）如果不是所有​都收敛，重复（2）。

上式中，表示数据(x,y)中的特征总数，(x,y)不同特征总数也可能不同，为了解决这个问题，定义松弛特征

$$s(x, y) = S - \sum_{i=1}^{n+1} \sum_{k=1}^{K} f_k(y_{i-1}, y_i, x, i)$$

S是一个常数，选择足够大的S使得对训练数据集的所有(x,y)，成立，此时特征总数可取S。从而可以求得对于转移特征，的更新方程是：

$$\delta_k=\frac{1}{S}\log \frac{E_{\tilde{P}}[t_k]}{E_P[t_k]}$$

其中，

$$E_p(t_k) = \sum_x \tilde{P}(x) \sum_{i=1}^{n+1} \sum_{y_{i-1}, y_i} t_k(y_{i-1}, y_i, x, i) \frac{\alpha_{i-1}^T(y_{i-1} \mid x) M_i(y_{i-1}, y_i \mid x) \beta_i(y_i \mid x)}{Z(x)}$$

对于状态特征，的更新方程是：

$$\delta_{K_1+l}=\frac{1}{S}\log \frac{E_{\tilde{P}}[s_l]}{E_P[s_l]}$$

其中，

$$E_p(s_i) = \sum_x \tilde{P}(x) \sum_{i=1}^n \sum_{y_i} s_i(y_i, x, i) \frac{\alpha_i^T(y_i \mid x) \beta_i(y_i \mid x)}{Z(x)}$$

以上算法称作算法S。此算法中需要使S取足够大，这样一来每步迭代的增量向量会变大，算法收敛会变慢。下面的算法T可以解决这个问题。算法T对每个观测序列x计算其特征总数的最大值，利用前向-后向递推公式，容易计算。此时关于转移特征的参数更新方程为：

$$E_{p_L}[t_k] = \sum_{x,y} \tilde{P}(x) P(y \mid x) \sum_{i=1}^{n+1} t_k(y_{i-1}, y_i, x, i) \exp(\delta_k T(x)) \\ = \sum_x \tilde{P}(x) \sum_y P(y \mid x) \sum_{i=1}^{n+1} t_k(y_{i-1}, y_i, x, i) \exp(\delta_k T(x)) \\ = \sum_x \tilde{P}(x) \alpha_{k,t} \exp(\delta_k \cdot t) \\ = \sum_{f=0}^{T_{\text{max}}} \alpha_{k,t} \beta_k^t$$

式中是特征的期待值，，是多项式方程的唯一实根，从而可求得。关于状态特征的参数更新方程同理：

$$E_{p}[s_1] = \sum_{x,y} \tilde{P}(x) P(y \mid x) \sum_{i=1}^{n} s_l(y_i, x, i) \exp(\delta_{K_1+l} T(x)) \\ = \sum_x \tilde{P}(x) \sum_y P(y \mid x) \sum_{i=1}^{n} s_l(y_i, x, i) \exp(\delta_{K_1+l} T(x)) \\ = \sum_x \tilde{P}(x) b_{l,t} \exp(\delta_k \cdot t) \\ = \sum_{t=0}^{T_{\text{max}}} b_{l,t} \gamma_t$$

• 拟牛顿法

输入：特征函数；经验分布；

输出：最优参数值；最优模型；

（1）选定初始点，取为正定对称矩阵，置k=0；

（2）计算。若，则停止计算；否则转（3）；

（3）由求出；

（4）一维搜索：求使得

$$f(w^{(k)}+\lambda_kp_k)=\min_{\lambda\geq 0}f(w^{(k)}+\lambda p_k)$$

（5）置；

（6）计算，若则停止计算，否则按下式求出：

$$B_{k+1}=B_k+\frac{y_ky_k^T}{y_k^T\delta_k}-\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k}$$

其中，；

（7）k=k+1，转（3）。

预测算法

• 条件随机场预测的维特比算法

输入：模型特征向量和权值向量w，观测序列；

输出：最优路径；

（1）初始化

$$\delta_1(j)=w\cdot F_1(y_0=start,y_1=j,x),\quad j=1,2,...,m$$

（2）递推，对

$$\delta_i(l)=\max_{1\leq j\leq m}\{\delta_{i-1}(j)+w\cdot F_i(y_{i-1}=j,y_i=l,x) \},\quad l=1,2,...,m$$ $$\Psi_i(l)=\arg \max_{1\leq j\leq m}\{\delta_{i-1}(j)+w\cdot F_i(y_{i-1}=j,y_i=l,x) \},\quad l=1,2,...,m$$

（3）终止

$$\max_y (w\cdot F(y,x))=\max_{1\leq j\leq m}\delta_n(j)$$ $$y_n^{\ast}=\arg \max_{1\leq j\leq m}\delta_n(j)$$

（4）返回路径

$$y_i^{\ast}=\Psi_{i+1}(y_{i+1}),\quad i=n-1,n-2,...,1$$

一个简单的线性CRF例子

import sklearn_crfsuite
from sklearn import metrics

# 定义训练数据和标签
X_train = [
    [{'token': 'I', 'pos': 'PRON'}, {'token': 'am', 'pos': 'VERB'}, {'token': 'fine', 'pos': 'ADJ'}],
    [{'token': 'How', 'pos': 'ADV'}, {'token': 'are', 'pos': 'VERB'}, {'token': 'you', 'pos': 'PRON'}]
]

y_train = [['PRON', 'VERB', 'ADJ'], ['ADV', 'VERB', 'PRON']]

# 创建 CRF 模型实例
crf = sklearn_crfsuite.CRF(algorithm='lbfgs', max_iterations=100, all_possible_transitions=True)

# 训练 CRF 模型
crf.fit(X_train, y_train)

# 测试数据
X_test = [[{'token': 'You', 'pos': 'PRON'}, {'token': 'are', 'pos': 'VERB'}, {'token': 'great', 'pos': 'ADJ'}]]
y_test = [['PRON', 'VERB', 'ADJ']]

# 预测标签
y_pred = crf.predict(X_test)

# 打印预测结果
print("Predicted labels:", y_pred[0])

# 评估模型性能
print(metrics.classification_report(y_test[0], y_pred[0]))