EM算法在高斯混合模型中的应用

EM 算法在高斯混合模型中的应用1.定义对于一个随机信号生成器，假设他的模型参数为Θ，我们能观测到的数据输出为X ，不能观测到的数据输出为Y ，且随机系统模型结构的概率密度函数为(,|)p x y Θ （1）能够观测到的一部分数据输出数据12{,,...,}N x x x ，模型的另一部分输出数据 未知,模型的参数Θ也未知。

EM 算法就是要求我们从观测数据12{,,...,}N x x x 中估计出参数Θ。

2.EM 算法的描述假设每一对随机系统的输出样本(,)n n x y 对于不同的n 相互独立，这样当(,,)p x y Θ，x 和y 都已知的情况下，概率(,,)p x y Θ也已知。

未观测的输出y 的概率分布也属于待求参数Θ。

根据独立性假设有：1(,|)(,|)Nn n n p x y p x y =Θ=Θ∏ （2）3.EM 算法的基本思路基本问题是求解下面的方程的解：arg max (,|)p x y Θ=Θ （3） 由于X 是确定量，Y 是未知的，因此即使给定了Θ，也无法求得(,|)p x y Θ的值，因此我们只能退一步求：arg max (|)p x Θ=Θ （4）其中(|)(,|)[(|),(|,)]y Y y Y p x p x y p y p x y ∈∈Θ=Θ=ΘΘ∑∑ （5） 表示考虑了未知数据y 的所有可能的取值Y 后对(|,)p x y Θ求平均值。

最后根据log 函数的单调性得到（4）的等效形式：arg max log (|)p x Θ=Θ （6） 对于(6)给出的最优化问题，考虑用下面的递推算法解决，即：先给定一个估值k Θ并计算(|)k p x Θ，然后更新k Θ得到1k +Θ并且有1log (|)log (|)k k p x p x +Θ>Θ （7）()log (|)log [(|)(|,)]|(|,)log (|,)(|,)(|)(|,)(|,)log (|,)(,)y Y k ky Y k ky Y k p x p y p x y p y p x y p y x p y x p y p x y p y x p y x B ∈∈∈Θ=ΘΘΘΘ⎡⎤=Θ⎢⎥Θ⎣⎦⎧⎫⎡⎤ΘΘ≥Θ⎨⎬⎢⎥Θ⎣⎦⎩⎭=ΘΘ∑∑∑ （8） 其中，等号在(,)k k B ΘΘ时成立，即：(,)log (|)k k k B p x ΘΘ=Θ （9）于是对log (|)p x Θ的递推算法（7）可通过(,)k B ΘΘ进行，步骤为： 1） 令k=0，先给出估值 k Θ2） 然后找出1k +Θ满足 1(,)(,)k k k k B B +ΘΘ>ΘΘ （10） 3） k 更新为k+1并返回步骤2)直到收敛令 1arg max (,)k k B +Θ=ΘΘ (11) 处理后[]{}[]{}1arg max (,)(|)(|,)arg max (|,)log (|,)arg max (|,)log (|)(|,)(|,)log (|,)arg max (|,)log (|)(|,)arg max (,)k k k ky Y k k k y Y k y Y k B p y p x y p y x p y x P y x p y p x y p y x p y x p y x p y p x y C +∈∈∈Θ=ΘΘ⎧⎫⎡⎤ΘΘ=Θ⎨⎬⎢⎥Θ⎣⎦⎩⎭=ΘΘΘ-ΘΘ=ΘΘΘ=ΘΘ∑∑∑ (12)其中[]{}(,)(|,)log (|)(|,)k k y Y C p y x p y p x y ∈ΘΘ=ΘΘΘ∑ (13) 4.EM 算法与高斯混合模型在随机系统模型中，假设m θ是通道m 的随机信号生成器的概率密度函数的参数，()p y m =是选中通道m 的概率。

记为m a 。

假设M 个随机信号生成器和通道选择随机生成器是相互独立的，从通道m 输出的数据x 的概率是：(|)m m m a p x θ (14)不考虑通信信息，输出x 的概率为：1(|)(|)Mm m m m p x a p x θ=Θ=∑ (15)其中：m θ：是第m 个通道随机信号生成器的参数。

Θ：参数集合{}1,2...,,m m m M a θ=。

观测数据为一批随机产生的输出信号，并且每个输出都是相互独立的，而每个输出来自哪个通道不可测。

于是系统模型参数估计问题就变为通过有限的输出样本12{,,...,}N x x x 估计M 个通道参数{},(1,2,...,)m m a m M θ=.应用（12）求解，其中(,)k C ΘΘ可以简化为：111112(,)log()(|,)log((|))(|,)(,)(,)M N M Nkkk m n m m m n m n m n k k C a p m x p x p m x C C θ====ΘΘ=Θ+Θ=ΘΘ+ΘΘ∑∑∑∑（16）其中：111(,)log()(|,)M Nkk m n m n C a p m x ==ΘΘ=Θ∑∑211(,)log((|))(|,)MNk k m m m n m n C p x p m x θ==ΘΘ=Θ∑∑这样我们把m a 和m p 分别放在两项里面，他们不相关，可以独立考虑。

在(,)k C ΘΘ中应用约束条件：11Mmm a==∑用拉格朗日乘子优化m a 得到：111(|,)Nk k mn n a p m x N +==Θ∑上式的含义是，选中m 号通道的概率估计1k m a +是每个观测数据n x 来自于m 通道的条件概率（根据上一次估值Θ估算）的平均。

其中的(|,)k n p m x Θ通过下式得出。

'''1(|)(|,)(|)kkm m n m n Mk n m m mm a p x p m x a p x θθ=Θ=∑ 2(,)k C ΘΘ中的m θ的优化取决于分布函数的类型，对于(|)m m m p x θ为高斯分布时，{},m m m θμσ=其中m μ是分布的均值，m σ是方差。

再经过推导，有：111(|,)Nk kmnn ap m x N+==Θ∑ ①111(|,)(|,)N k nn k n m Nk nn xp m x p m x μ+==Θ=Θ∑∑ ， ②1/212111(|,)||=(|,)N k k n n m k n m Nk n n p m x x p m x μσ++==⎛⎫Θ- ⎪ ⎪ ⎪Θ ⎪⎝⎭∑∑ ③m 通道参数{},m m μσ得更新可以看作是对n x 的加权，加权系数(|,)k m p m x Θ可以看成是根据上一次的参数估计k Θ算出来得n x 率属于m 通道的概率。

最后，上面的EM 算法可能收敛到局部极大点，因此需要选择多个参数Θ的初始值进行迭代计算，并选择使得(|)p x Θ最大的解，(|)p x Θ最大的解可由下式算出：[]12111111(|)(|)(|,)(|)(|)nm N y YNM MM m n m m m m n NM m n m m n p x p y p x y ap x a p x θθ∈======Θ=ΘΘ⎡⎤=⎣⎦⎡⎤=⎢⎥⎣⎦∑∑∑∑∏∑∏5.EM 算法在matlab 中的实现利用上面推导出的公式①②③，我们以二个一维的高斯分布（1N ,2N ）来验证EM 算法的性能，首先用二个一维的高斯分布来建立一个高斯混合模型H N 。

假设：2111(,)N N μσ ,2222(,)N N μσ1122H N N N αα=+其中1α与2α为混合系数，且有121αα+=，我们要用EM 算法估计混合系数和各一维高斯分布的均值和方差22121212(,,,,,)ααμμσσ。

并将利用EM 算法估计出的值与真实值做比较，就可以得到该算法的性能。

首先我们取22121212(,,,,,)ααμμσσ的真实值为（0.4,0.6,1,2,0.25,0.36） 这样我们得到一个混合高斯分布，他的密度函数为H N ,然后产生1000个服从H N 的分布的观测样本点。

接下来要做的就是对这1000个样本点用EM 算法进行处理，来估计出一组22121212(,,,,,)ααμμσσ的值。

在使用EM 算法时，要首先给22121212(,,,,,)ααμμσσ设定一组初值这里假设初值为1α=0.3，2α=0.7，1μ=0.8，2μ=1.8，21σ=0.2，22σ=0.25 Matlab 具体程序如下：Y=zeros(1,10000); for i=1:10000if rand(1)>0.3Y(i)=normrnd(2,sqrt(0.36),1,1); elseY(i)=normrnd(1,sqrt(0.25),1,1); endend %高斯混合模型A=[0.3 0.7]; %设置参数 α的初值 M=[0.8 1.8]; %设置均值 μ的初值 S=[0.2 0.25]; %设置方差 2σ的初值 for n=1:1000 for j=1:2 a3=0; a4=0; a5=0;for k=1:10000 a1=0;for t=1:2a1=A(t)*1/sqrt(2*pi*S(t))*exp(-(Y(k)-M(t))^2/(2*S(t)))+a1;endf= A(j) * 1/sqrt(2*pi*S(j))*exp(-(Y(k)-M(j))^2/(2*S(j))); a2=f/a1;a3=a2+a3; %a3对应公式1(|,)Nknn p m x =Θ∑a4=a2*Y(k)+a4; %a4对应公式1(|,)Nk nn n p m x x =Θ∑a5=a2*(Y(k)-M(j))^2+a5; %a5对应公式121(|,)()Nk k n n m n p m x x μ+=Θ-∑endA(j)=a3/10000; %循环更新系数值M(j)=a4/a3; %循环更新均值值 S(j)=a5/a3; %循环更新方差值 end end运行程序，查看变量A,M,S 的值，与真实值比较一下，就可以得到用EM 算法估计的高斯混合模型的性能了。

得到参数为A=[ 0.3063 0.6937]，M=[ 1.0093 2.0013]，S=[ 0.2558 0.3632]，而真实值为A=[0.4 0.6]，M=[1 2]，S=[0.25 0.36] ，我们可以从中看出用EM 算法估计的高斯混合模型的相关参数与真实值是比较接近的。