EM 算法及其在高斯混合模型中的应用 EM 算法在高斯混合模型中的应用 1.定义 对于一个随机信号生成器，假设他的模型参数为  ，我们能观测到的数据输 出为 X，不能观测到的数据输出为 Y，且随机系统模型结构的概率密度函数为 ) ( , p x y  | （1） { , x x 能够观测到的一部分数据输出数据 1 2 ,..., { , x x 未知,模型的参数  也未知。EM 算法就是要求我们从观测数据 1 2 计出参数  。 2.EM 算法的描述 }N x ，模型的另一部分输出数据 }N x 中估 ,..., 假设每一对随机系统的输出样本( , ) n p x y  ，x和y都已知的情况下，概率 ( , ( , 分布也属于待求参数  。 根据独立性假设有： ) x y 对于不同的n相互独立，这样当 n , p x y  也已知。未观测的输出y的概率 ) , ( , p x y |   ) N  n 1  ( , p x y n n |  ) （2） 3.EM算法的基本思路 基本问题是求解下面的方程的解： arg max ( , p x y   | )  （3） | p x y  的 ) 由于X是确定量，Y是未知的，因此即使给定了  ，也无法求得 ( , 值，因此我们只能退一步求：   arg max ( p x | )  其中 | ) ( p x   [ ( p y 表示考虑了未知数据y的所有可能的取值Y后对 ( | p x y  求平均值。 最后根据log函数的单调性得到（4）的等效形式： ( p x y ) ),  , ( , p x y   )]  y Y  y Y  ) , | | |     arg max log ( p x | )  （4） （5） （6） 对于(6)给出的最优化问题，考虑用下面的递推算法解决，即：先给定一个估值 k 并计算 ( p x  ，然后更新 k 得到 1k  并且有 )k | log ( p x | 1 k  )  log ( p x | k  ) （7） log ( p x |   ) log  [ y Y   log  y Y     ( , p y x |  k ) | ( p y  p y ( p y ( p y x | ,  k )log    k )  y Y      B    , k ( k , )]  )  , | ) ( , p x y   | ( | p x y  k | ( , ) p y x  | , ( ) p x y  k ) ( | , p y x |     )     （8）    其中，等号在 ( B   时成立，即： log ( p x )    B ( ) , k k | k  ) （9） 

EM 算法及其在高斯混合模型中的应用 | ) p x  的递推算法（7）可通过 ( 于是对 log ( 1） 令k=0，先给出估值 k ( ( B 2） 然后找出 1k  满足 3） k更新为k+1并返回步骤2)直到收敛 ) k   令 处理后 arg max ( B        B ) , , , 1 1 k k k k B   进行，步骤为： )k , k ) （10） (11)   k 1  ) y Y  ( , | p y x k   argmax argmax ( , B       , | ( P y x   argmax ( , | p y x y Y  k argmax ( , ) C   argmax y Y      k  )log k  )log k  )log ( p y     ( p y  ( p y )   ) ( | | ,  p x y    k ) ( | , p y x     , | ( | ) ( | , ) p y x p x y      ) ( | , p x y )   | k  )log ( | , p y x k   ) (12) 其中 C (    k ) ,  y Y  4.EM算法与高斯混合模型  ( , | p y x k  )log  ( p y | )   ( p x y | ,  )  (13) 在随机系统模型中，假设 m 是通道 m 的随机信号生成器的概率密度函数的 参数， ( p y m 是选中通道 m 的概率。记为 ma 。 ) 假设 M 个随机信号生成器和通道选择随机生成器是相互独立的，从通道 m 输出的数据 x 的概率是： ( 不考虑通信信息，输出 x 的概率为：    ) ( p x M | a p x  m m m | ) a p x  m m m ( | m 1  ) (14) (15) 其中： m ：是第 m 个通道随机信号生成器的参数。  ：参数集合 ,m a  1,2..., m m 。 M  观测数据为一批随机产生的输出信号，并且每个输出都是相互独立的，而每 个输出来自哪个通道不可测。于是系统模型参数估计问题就变为通过有限的输出 { , x x 样本 1 2 应用（12）求解，其中 ( x 估计 M 个通道参数 }N  m )k C   可以简化为： ma m 1,2,..., ,..., M   . ( ) , , C ( M N k , )    log(  1 1 m n   ( ( C      ) , k 2 k ) C 1 ,  其中： a ) ( p m x | , n m k   ) M N  m 1  n 1  log( ( p x m m |  m )) ( p m x | , k  ) n （16） 

EM 算法及其在高斯混合模型中的应用 C 1 (    k ) , C 2 (    k ) , M N  n 1 1 m   M N  m 1  n 1  log( a ) ( p m x | , k  ) n m log( ( p x m m |  m )) ( | p m x , k  ) n 这样我们把 ma 和 mp 分别放在两项里面，他们不相关，可以独立考虑。 在 ( C   中应用约束条件： )k , M  a m 1 m 1  用拉格朗日乘子优化 ma 得到： 1 k   m a 1 N N  n 1  ( | p m x n , k  ) 上式的含义是，选中 m 号通道的概率估计 1k 道的条件概率（根据上一次估值  估算）的平均。其中的 ( 得出。 ma  是每个观测数据 nx 来自于 m 通 p m x  通过下式 )k , | n ( p m x | , n k   ) | ( a p x n m m  ( x a p m m M ' ' k )  m k |  m n m 1  ) ' C   中的 m 的优化取决于分布函数的类型，对于 ( p x  为高斯分 2( m 布时， )k ) m m , |    m m m ,   其中 m 是分布的均值， m 是方差。再经过推导，有： a 1 k  m  ( | p m x n , k  ) N  n 1  m 通道参数 ③ m  得更新可以看作是对 nx 的加权，加权系数 ( ,m p m x  | 可以看成是根据上一次的参数估计 k 算出来得 nx 率属于 m 通道的概率。  , m )k 最后，上面的EM算法可能收敛到局部极大点，因此需要选择多个参数  的初 p x  最大的解可由下式算 p x  最大的解， ( ) ) | | 始值进行迭代计算，并选择使得 ( 出： 1 N N   1  N n 1  n N  n 1        1 k   m  1 k   m = | x p m x n n ( , k  ) ， ( p m x | n , k  ) ( p m x | n , k  ) | x n  k  m 1 2 |  N  n 1  ( p m x | , k  ) n       ① ② 1/2 

EM 算法及其在高斯混合模型中的应用 ( p x |   )  ( p y | )  ( p x y | , )   M M      N m N 1  n 1    ( a p x m n n |  m m )    y Y  M   1  M m 1 N 21 m      n 1  m 1  ( a p x m n |  m )    5.EM算法在matlab中的实现 利用上面推导出的公式①②③，我们以二个一维的高斯分布（ 1N , 2N ）来 验证EM算法的性能，首先用二个一维的高斯分布来建立一个高斯混合模型 HN 。 假设： 2 , N 1 HN ) ( N    , 1 N     2   其中 1 与 2 为混合系数，且有 1 2 , , 1 N N 1 1 2 2 各一维高斯分布的均值和方差 , 值与真实值做比较，就可以得到该算法的性能。 , 1 2 1 2 1  N   2 ( , 2 2 )  ，我们要用EM算法估计混合系数和 ( ,       。并将利用EM算法估计出的 ) 2 2 1 2 2 1 2 1 , , , , (       的真实值为（0.4,0.6,1,2,0.25,0.36） 首先我们取 这样我们得到一个混合高斯分布，他的密度函数为 HN ,然后产生1000个服从 HN 的分布的观测样本点。接下来要做的就是对这1000个样本点用EM算法进行处 理，来估计出一组 ) , 2 1 2 2 2 2 1 2 1 , , , , (       的值。 , 在使用EM算法时，要首先给 这里假设初值为 1 =0.3， 2 =0.7， 1  0.8， 2  1.8， 2 Matlab具体程序如下： 1 (       设定一组初值 ) , ) , , , , 2 2 2 2 1 2 1 2 1 2 1  0.2， 2 2  0.25 Y=zeros(1,10000); for i=1:10000 if rand(1)>0.3 Y(i)=normrnd(2,sqrt(0.36),1,1); Y(i)=normrnd(1,sqrt(0.25),1,1); else end end A=[0.3 0.7]; M=[0.8 1.8]; S=[0.2 0.25]; for n=1:1000 for j=1:2 a3=0; a4=0; a5=0; for k=1:10000 a1=0; for t=1:2 %高斯混合模型 %设置参数 的初值 %设置均值 的初值 %设置方差 2 的初值 a1=A(t)*1/sqrt(2*pi*S(t))*exp(-(Y(k)-M(t))^2/(2*S(t)))+a1; 

EM 算法及其在高斯混合模型中的应用 end f= A(j) * 1/sqrt(2*pi*S(j))*exp(-(Y(k)-M(j))^2/(2*S(j))); a2=f/a1; a3=a2+a3; %a3 对应公式 N  n 1  ( p m x | n , k  ) a4=a2*Y(k)+a4; %a4 对应公式 N  n 1  ( p m x | , k  ) x n n a5=a2*(Y(k)-M(j))^2+a5; %a5 对应公式 N  n 1  ( p m x n | , k  )( k x   n m  1 2 ) end A(j)=a3/10000; M(j)=a4/a3; S(j)=a5/a3; end end %循环更新系数值 %循环更新均值值 %循环更新方差值 运行程序，查看变量 A,M,S 的值，与真实值比较一下，就可以得到用 EM 算法估计的高 0.6937]，M=[ 1.0093 2.0013]，S=[ 0.2558 0.36] ，我们可以从中看出用 EM 算法 斯混合模型的性能了。得到参数为 A=[ 0.3063 0.3632]，而真实值为 A=[0.4 0.6]，M=[1 2]，S=[0.25 估计的高斯混合模型的相关参数与真实值是比较接近的。