基于高斯混合模型案例的EM算法教学设计
Teaching Design of EM Algorithm Based on Gaussian Mixture Model Case
DOI: 10.12677/ces.2025.131053, PDF, HTML, XML,   
作者: 杜 芳, 方晓峰:火箭军工程大学基础部数学教研室,陕西 西安
关键词: EM算法高斯混合模型教学设计EM Algorithm Gaussian Mixture Model Teaching Design
摘要: EM (Expectation Maximization)算法是统计学中的核心算法,也是本校近代数理统计课程教学过程中的一个重难点。论文采用案例式、启发式、研讨式教学方法,以基于高斯混合模型(GMM)的轴承退化阶段划分问题为例,引导学生发现隐变量模型极大似然估计(MLE)存在的困难,设计问题链启发学生探寻参数估计的数值方法,并总结出EM算法的一般过程。基于matlab编程可视化EM算法下的GMM模型参数更新过程,对比MLE目标函数和EM迭代目标函数,分析EM算法的内涵思想并结合图形进行直观展示,并且挖掘其中蕴含的思政元素,在知识传授的同时实现价值塑造。
Abstract: Expectation maximization (EM) algorithm is a core algorithm in statistics and also a key and difficult point in the teaching process of modern mathematical statistics courses in our school. The paper adopts a case-based and heuristic teaching method, taking the Gaussian Mixture Model (GMM) based bearing degradation stage division problem as an example, guiding students to discover the difficulties of maximum likelihood estimation (MLE) in the latent variable model, designing a problem chain to inspire students to explore numerical methods for parameter estimation, and summarizing the general process of EM algorithm. Based on Matlab programming, the parameter update process of GMM based on EM algorithm is visualized. Comparing the MLE objective function and EM iteration objective function, the intrinsic thought of EM algorithm is analyzed and visually displayed with graphics. The ideological and political elements are also explored, so as to achieve value shaping while knowledge transmission.
文章引用:杜芳, 方晓峰. 基于高斯混合模型案例的EM算法教学设计[J]. 创新教育研究, 2025, 13(1): 394-402. https://doi.org/10.12677/ces.2025.131053

1. 引言

随着国防和军队改革的深入推进,全军兴起加强军事人才建设热潮,也对军校研究生教育提出了新的要求:针对服务实战中人才培养的需要,从学术型人才培养向应用型人才培养转变,保证国防和军队改革后的战斗力[1]。而当前军校研究生培养面临着课程内容理论性较强、科学研究工作与人才培养的结合不紧、学生自身积极主动性不够等问题。如何从科学研究和服务部队两个方面出发,建设以学生能力素质提升为主线,以创新能力培养为核心的课程,是军校研究生教育的重点任务[2]

近代数理统计是面向我校研究生开设的一门选修课程,课程聚焦于统计学中的概念、理论和方法,为学生依据实际问题进行建模或通过统计推断发现问题的内在规律提供科学方法,支撑研究生的专业方向课程学习,启迪研究生分析问题和解决问题的思维[3]。在学情方面,选修该门课程的学生已经选定研究方向,具备了一定的文献查阅能力、问题分析能力和科研探索能力,但是数理统计基础参差不齐,容易对数理统计中的抽象原理望而生畏,同时非常关心如何应用数理统计方法解决实际问题,希望从课程中学习到能应用于专业领域的知识和方法[4]

在工程实践中,存在大量对观测不完整数据或者含有隐变量的模型做参数估计的问题。直接采用极大似然估计方法通常无法获得参数估计的解析表达式,EM算法是针对此类问题的一种有效的数值求解方法[5]。EM算法是统计学中的核心算法,是高斯混合模型、隐马尔可夫模型、深度生成模型等众多机器学习方法的理论基础,也是近代数理统计课程教学中的一个重难点。

基于这样的国防需求、课程目标和学情特点,我们采用案例式、启发式、研讨式教学方法进行EM算法的教学设计,在提高课程的专业应用性的同时,注重学生多维能力的培养。本次课以基于GMM的轴承退化阶段划分问题为引例,引导学生发现在对不完整观测数据或隐变量模型做极大似然估计时面临的困境以及设计数值解法的必要性。拆解问题并组织学生进行课前思考和课上讨论,使其逐步探索总结出EM算法的迭代更新过程,进而应用算法解决引例问题,再结合matlab编程使学生对迭代过程有更加直观的认识。最后,对比MLE目标函数和EM迭代目标函数,深入分析EM算法的思想内涵,证明EM算法的有效性和收敛性。具体实施过程在2~5节展开介绍。

2. 基于GMM的轴承退化阶段划分问题

2.1. 问题阐述

轴承寿命预测对于开展机械设备健康管理,提高旋转机械运行安全性、可靠性具有重要意义,其中轴承退化阶段的划分能够有效提高其寿命预测的精度[6]。现实应用中通常缺乏大量故障标记数据,经常采用无监督的聚类算法进行轴承退化阶段的划分,高斯混合模型(GMM)就是一种有效的聚类算法:

p( x,θ )= k=1 K P( y k =1 ) p X|Y ( x| y k =1 ) (1)

其中随机变量Y表示轴承所属退化阶段,其取值为K维one-hot编码向量 ( y 1 , y 2 ,, y k ) y k =1 代表轴承属于第k个退化阶段,每个退化阶段的轴承服从多元高斯分布,即

p X|Y ( x| y k =1 )= ( 2π ) n 2 | Σ k | 1 2 e 1 2 ( x μ k ) T Σ k 1 ( x μ k ) ,k=1,2,,K (2)

为简化后续表达,将(1)式记为

p( x,θ )= k=1 K α k p k ( x, θ k ) (3)

其中 θ=( α 1 , θ 1 , α 2 , θ 2 ,, α K , θ K ) 为总体参数,且 k=1 K α k =1 ,于是 Y~PN( 1: α 1 , α 2 ,, α K )

轴承退化阶段划分的目的是通过观测数据 x ˜ =( x 1 , x 2 ,, x N ) ,推断出每个成分的参数 θ k 以及每个轴承 x i ( i=1,2,,N ) 所属的退化阶段 y i 。事实上,只要(3)式中的所有参数已知,即可得到 y i 的条件分布为:

P Y|X ( y|x,θ )= P( y,θ ) P X|Y ( x|y,θ ) y P( y,θ ) P X|Y ( x|y,θ ) (4)

所以解决轴承退化阶段划分的核心问题还是要通过观测 x ˜ =( x 1 , x 2 ,, x N ) 估计总体参数 θ=( α 1 , θ 1 , α 2 , θ 2 ,, α K , θ K )

2.2. GMM参数估计的困难

课前向学生发布预习问题1~3,课堂中组织学生交流汇报。

问题1. 已知观测数据 x ˜ =( x 1 , x 2 ,, x N ) 时,请尝试采用MLE方法估计总体参数。

不难发现对数似然函数具有如下形式:

i=1 N logp( x i ,θ ) = i=1 N log k=1 K α k p k ( x i , θ k ) (5)

其中包含求和项的对数,在关于参数求偏导时,导数形式复杂,很难或者无法求得参数估计的解析表达式。而当问题无法解析求解时,我们通常会尝试迭代的数值算法,EM算法就是一种针对不完整观测数据或隐变量模型参数估计的数值迭代算法。

3. EM算法迭代过程探索

数值迭代算法,通常预先给定一个初始值,通过分析目标函数在当前取值下的性质,构造出下一步的参数迭代公式。这里我们采用类似的方式,考虑下面两个问题。

3.1. 模型参数给定时隐变量(轴承退化阶段)后验分布的推断

问题2. 当已知参数为 θ t =( α 1 t , θ 1 t , α 2 t , θ 2 t ,, α K t , θ K t ) 时,分析观测为 x ˜ =( x 1 , x 2 ,, x N ) 时,退化阶段变量Y的后验分布是什么?

对学生的讨论结果进行总结,对于问题2,有:

P( y i k =1| x i , θ t )= α k t p k ( x i , θ k t ) k=1 K α k t p k ( x i , θ k t ) (6)

P( y i k =1| x i , θ t ) 记为 y kt ,则 Y i ~PN( 1: y i 1t , y i 2t ,, y i Kt ) ,即隐变量的后验分布和先验分布一样,属于多项分布。

3.2. 隐变量(轴承退化阶段)已知时(3)式的MLE

问题3. 当隐变量也可观测时,根据观测数据 ( x ˜ , y ˜ )=( x 1 ,, x N , y 1 ,, y N ) 能否对(3)式中的参数做MLE?

总结学生讨论结果,轴承特征随机变量X和退化阶段随机变量Y的联合分布如下:

p( x,y,θ )=P( y,θ ) p X|Y ( x|y,θ )= α k y k k=1 K [ p k ( x, θ k ) ] y k = k=1 K [ α k p k ( x, θ k ) ] y k (7)

已知观测数据 ( x ˜ , y ˜ )=( x 1 ,, x N , y 1 ,, y N ) 时,似然函数为:

L( θ )=p( x ˜ , y ˜ ,θ )= i=1 N p( x i , y i ,θ ) = i=1 N k=1 K [ α k p k ( x i , θ k ) ] y i k (8)

对似然函数求对数,得对数似然函数为:

logL( θ )= i=1 N k=1 K y i k log α k + y i k log p k ( x i , θ k ) (9)

对数似然函数关于参数求偏导,可得:

logL( θ ) θ k = i=1 N y i k log p k ( x i , θ k ) θ k (10)

其中 y i k 取值为0或1,相当于只利用样本中属于第k个退化阶段的轴承特征,来估计GMM中第k个成分的参数 θ k

对于成分占比参数 α k ,极大化(9)式相当于求解如下约束优化问题:

max α 1 ,, α K k=1 K i=1 N y i k log α k    s.t. k=1 K α k =1 (11)

利用拉格朗日乘子法,可得 α k = 1 N i=1 N y i k ,即样本中属于第k个退化阶段的个体所占比例。

3.3. EM算法的迭代过程

问题2和问题3是比较容易解决的,也就是说在总体参数给定的情况下,可以根据观测轴承推断其所属退化阶段的后验分布,在轴承退化阶段也可观测的情况下,能够通过MLE得到总体参数的显式估计。如果交替进行这两步,就可以得到参数的一种迭代更新过程,但是问题2中给出的是分布,而问题3中用到的是观测值,如何衔接这两个过程呢?

联系随机变量的分布和随机变量的观测值,有两种简单途径:第一种是从分布中随机抽样得到观测值,第二种是按照分布对所有可能得观测值求期望。两种途径各有优缺点,抽样法能保证可行的数值计算,但是需要设计抽样方法缓解以偏概全现象,期望法比较全面,但是当分布形式复杂时,难以获得显式表示。当采用第二种途径时,对应的迭代过程就是EM算法:

E步:当参数 θ t =( α 1 t , θ 1 t , α 2 t , θ 2 t ,, α K t , θ K t ) 时,计算完整观测下的似然函数 p( x,y,θ ) 或者对数似然函数 logp( x,y,θ ) 关于条件分布 P( y|x, θ t ) 的期望:

E Y ( p( x,Y,θ )|x, θ t ) E Y ( logp( x,Y,θ )|x, θ t ) (12)

上述期望可以看作是参数 θ 的函数,记为 Q( θ, θ t )

M步: Q( θ, θ t ) 极大化,得到参数迭代公式

θ t+1 =arg max θ Q( θ, θ t ) (13)

4. 基于EM算法的引例问题求解

4.1. GMM参数迭代公式推导

回到GMM中,(8)式中的似然函数关于条件分布 P( y ˜ | x ˜ , θ t ) 的期望为:

E Y ( L( θ )| x ˜ , θ t )= i=1 N k=1 K [ α k p k ( x i , θ k ) ] E Y ( y i k | x i , θ t ) = i=1 N k=1 K [ α k p k ( x i , θ k ) ] y i kt (14)

(9)式中的对数似然函数关于条件分布 P( y ˜ | x ˜ , θ t ) 的期望为:

E Y ( logL( θ )| x ˜ , θ t )= i=1 N k=1 K y i kt log α k + y i kt log p k ( x i , θ k ) (15)

E Y ( logL( θ )| x ˜ , θ t ) Q( θ, θ t ) ,则其关于参数的偏导为:

Q( θ, θ t ) θ k = i=1 N y i kt log p k ( x i , θ k ) θ k (16)

p k ( x i , θ k ) 替换为(2)式中的多元高斯分布可得:

Q( θ, θ t ) μ k = i=1 N y i kt Σ k 1 ( μ k x ) (17)

Q( θ, θ t ) Σ k = i=1 N y i kt [ Σ k 1 + Σ k 1 ( x μ k ) ( x μ k ) T Σ k 1 ] (18)

令偏导为0,即可获得如下参数更新公式:

μ k t+1 i=1 N y i kt x i i=1 N y i kt                 Σ k t+1 i=1 N y i kt ( x i μ k t+1 ) ( x i μ k t+1 ) T i=1 N y i kt (19)

参数 α k 更新方式与3.1中类似, α k t+1 1 N i=1 N y i kt

这些迭代公式和单高斯分布参数极大似然估计的公式非常类似,差别是GMM迭代公式中为每个观测值赋予了权重 y i kt

4.2. 基于matlab的GMM参数迭代过程展示

通过matlab编程实现上述迭代过程,图1展示了包含两个成分的一维高斯混合模型的参数和类别划分迭代过程。其中黑色曲线代表真实的GMM模型概率密度曲线,红色和蓝色曲线分别代表学习到的两个一维高斯分布成分,底部的红色“*”和蓝色“。”代表两种数据划分结果。可以看到随着迭代次数的增加,学习到的两个成分期望逐渐逼近真实数据分布中的期望,数据划分结果也逐渐接近真实数据类型。

5. EM算法的理论分析

前面几部分启发式地探寻了EM算法的迭代更新步骤,并且通过一维高斯混合模型示例直观展示了EM算法的效果。下面我们从理论角度分析证明EM算法的思想内涵和有效性。

5.1. EM算法的思想

首先,比较EM迭代优化目标 Q( θ, θ t )= E Y ( logp( x,Y,θ )|x, θ t ) 和极大似然估计的优化目标 logp( x,θ )

Figure 1. Iterative update process of one-dimensional GMM

1. 一维混合高斯模型的迭代更新过程

以连续型隐变量为例,也就是比较 logp( x,θ ) logp( x,y,θ )p( y|x, θ t )dy ,我们选择从 logp( x,θ ) 出发逐步向 logp( x,y,θ )p( y|x, θ t )dy 靠拢:

logp( x,θ )=log p( x,y,θ )dy =log p( x,y,θ ) p( y|x, θ t ) p( y|x, θ t )dy (20)

其中 logx 为凹函数,根据Jensen不等式,可得:

logp( x,θ ) log p( x,y,θ ) p( y|x, θ t ) p( y|x, θ t )dy = logp( x,y,θ )p( y|x, θ t )dy logp( y|x, θ t )p( y|x, θ t )dy =Q( θ, θ t )+H[ p( y|x, θ t ) ] (21)

其中 H[ p( y|x, θ t ) ] 和优化参数 θ 无关,(21)式说明EM迭代优化目标 Q( θ, θ t ) 是MLE优化目标 logp( x,θ ) 的下界,也就是说EM算法中每个E步,都是在当前参数 θ t 下求 logp( x,θ ) 的下界函数,再用下界函数的极大值点来更新参数,过程如图2所示。

虽然无法直接求得 logp( x,θ ) 的极大值点,但是按照这样的学习过程,在迭代过程中参数会朝着极大值点不断靠近,最终能够达到极大值点的任意小邻域内。这样的学习过程也蕴含着人生哲理:不积跬步无以至千里,不积小流无以成江海,在学习和工作中,每天积累一点,我们就会朝着目标不断靠近。

5.2. EM算法的收敛性

图2中直观来看,EM算法之所以有效,是因为随着参数的迭代更新, logp( x,θ ) 在不断提高,这一结果是可证明的。

定理:EM算法在每一次迭代后,都有 logp( x, θ t+1 )logp( x, θ t )

Figure 2. Parameter learning process of EM algorithm

2. EM算法参数学习过程

证明:根据条件概率公式,有:

p( x,θ )= p( x,y,θ ) p( y|x,θ ) (22)

两边取对数,则有:

logp( x,θ )=logp( x,y,θ )logp( y|x,θ ) (23)

两边同时关于条件分布 P( y|x, θ t ) 求期望,可得:

logp( x,θ )p( y|x, θ t )dy =logp( x,θ ) = logp( x,y,θ )p( y|x, θ t )dy logp( y|x,θ )p( y|x, θ t )dy =Q( θ, θ t ) logp( y|x,θ )p( y|x, θ t )dy (24)

θ t θ t+1 代入 logp( x,θ ) 再相减,可得:

logp( x, θ t+1 )logp( x, θ t )=[ Q( θ t+1 , θ t )Q( θ t , θ t ) ] [ logp( y|x, θ t+1 )p( y|x, θ t )dy logp( y|x, θ t )p( y|x, θ t )dy ] (25)

因为 θ t+1 =arg max θ Q( θ, θ t ) ,所以等式右边第一项 Q( θ t+1 , θ t )Q( θ t , θ t )0 ,将等式右边第第二项记为K,对其做等价变换,有:

K= log p( y|x, θ t+1 ) p( y|x, θ t ) p( y|x, θ t )dy (26)

根据Jensen不等式,有:

Klog p( y|x, θ t+1 ) p( y|x, θ t ) p( y|x, θ t )dy =log p( y|x, θ t+1 )dy =log1=0 (27)

综上所述, logp( x, θ t+1 )logp( x, θ t )0logp( x, θ t+1 )logp( x, θ t )

再根据单调有界原理,若 logp( x,θ ) 有上界,则EM算法收敛到 logp( x,θ ) 的稳定点。从图2也可以看出,EM算法不能保证收敛到全局最优点,而且收敛位置受到初值影响。较为可行的办法是选取几个不同的初值进行迭代,然后在多个估计间加以选择,减轻初值选取对结果的影响[7]

6. 教学效果分析

学生通过查阅资料和课前讨论,针对教师课前给出的3个问题做了充分准备,并且在授课过程中主动发言,课堂参与的主动性和积极性被充分调动。此外学生的科研兴趣也得到了激发,在课后进一步追踪EM算法的研究进展及其在前沿机器学习方法中的应用。批改课后大作业发现,作业的结构完整度、统计方法和专业应用的契合度都有明显提升,学生基本都能通过编程实现算法。此外,通过微信群向选修课程的12名学生发出调查问卷,请他们对本节课的教学效果进行评价,统计结果分析见表1

Table 1. Feedback on class effectiveness of EM algorithm

1. EM算法课堂效果反馈

调查项目

备选项及对应占比

对教学模式是否满意?

A. 非常满意(占比91.7%)

B. 满意(占比8.3%)

C. 不满意(占比0%)

对教学内容是否满意?

A. 非常满意(占比100%)

B. 满意(占比0%)

C. 不满意(占比0%)

对能力提升是否有用?

A. 非常有用(占比83.3%)

B. 比较有用(占比16.7%)

C. 不太有用(占比0%)

对专业研究是否有帮助?

A. 非常有帮助(占比83.3%)

B. 比较有帮助(占比16.7%)

C. 帮助不大(占比0%)

7. 结语

本节课以与学生专业相关的轴承退化阶段划分问题为引例,结合课前问题,引导学生发现隐变量模型极大似然估计中的困难,从而探索迭代数值算法,总结EM算法的一般步骤,并应用其解决引例问题,培养了学生分析问题、解决问题的能力,锻炼了学生的探索思维,提高了学生的创新意识。最后通过理论推导,分析EM算法的内涵和收敛性,培养学生严谨的科研态度,并且挖掘算法中蕴含的思政元素,在知识传授、能力培养的同时实现价值塑造。学生课堂上主动讨论,课后积极探索并认真完成大作业,对于课堂效果满意度较高,为课程的进一步建设提供了方向。

致 谢

感谢相关文献对本文的启发以及审稿专家提出的宝贵意见。

参考文献

[1] 赵慎, 王红云, 郭希维, 张自宾. 军校研究生教育应用型培养模式探讨[J]. 中国教育技术装备, 2018(22): 92-93+98.
[2] 安涛, 张腾, 何宇廷. 聚焦备战打仗型军校研究生培养模式的探索与实践[J]. 大学教育, 2024(14): 111-114.
[3] 肖枝洪, 黄守成. 人工智能时代下研究生应用数理统计优质课程建设[J]. 大学数学, 2024, 40(2): 41-46.
[4] 宫春梅. 基于创新能力培养的研究生《数理统计》课程教学研究[J]. 创新创业理论研究与实践, 2020(16): 25-26.
[5] Dempster, A.P., Laird, N.M. and Rubin, D.B. (1977) Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society Series B: Statistical Methodology, 39, 1-22.
https://doi.org/10.1111/j.2517-6161.1977.tb01600.x
[6] 陈东楠, 胡昌华, 郑建飞, 郑红倩, 裴洪. 基于一种新型健康指标的轴承退化阶段自适应划分方法[J]. 哈尔滨工程大学学报, 2025(9): 1-13.
[7] 茆诗松, 王静龙, 濮晓龙. 高等数理统计[M]. 第3版. 北京: 高等教育出版社, 2022.