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

EM算法在高斯混合模型中的应用1.定义对于一个随机信号生成器，假设他的模型参数为，我们能观测到的数据输出为X，不能观测到的数据输出为Y，且随机系统模型结构的概率密度函数为(,|)pxy（1）能够观测到的一部分数据输出数据12{,,...,}Nxxx，模型的另一部分输出数据未知,模型的参数也未知。EM算法就是要求我们从观测数据12{,,...,}Nxxx中估计出参数。2.EM算法的描述假设每一对随机系统的输出样本(,)nnxy对于不同的n相互独立，这样当(,,)pxy，x和y都已知的情况下，概率(,,)pxy也已知。未观测的输出y的概率分布也属于待求参数。根据独立性假设有：1(,|)(,|)Nnnnpxypxy（2）3.EM算法的基本思路基本问题是求解下面的方程的解：argmax(,|)pxy（3）由于X是确定量，Y是未知的，因此即使给定了，也无法求得(,|)pxy的值，因此我们只能退一步求：argmax(|)px（4）其中(|)(,|)[(|),(|,)]yYyYpxpxypypxy（5）表示考虑了未知数据y的所有可能的取值Y后对(|,)pxy求平均值。最后根据log函数的单调性得到（4）的等效形式：argmaxlog(|)px（6）对于(6)给出的最优化问题，考虑用下面的递推算法解决，即：先给定一个估值k并计算(|)kpx，然后更新k得到1k并且有1log(|)log(|)kkpxpx（7）log(|)log[(|)(|,)]|(|,)log(|,)(|,)(|)(|,)(|,)log(|,)(,)yYkkyYkkyYkpxpypxypypxypyxpyxpypxypyxpyxB（8）其中，等号在(,)kkB时成立，即：(,)log(|)kkkBpx（9）于是对log(|)px的递推算法（7）可通过(,)kB进行，步骤为：1）令k=0，先给出估值k2）然后找出1k满足1(,)(,)kkkkBB（10）3）k更新为k+1并返回步骤2)直到收敛令1argmax(,)kkB(11)处理后1argmax(,)(|)(|,)argmax(|,)log(|,)argmax(|,)log(|)(|,)(|,)log(|,)argmax(|,)log(|)(|,)argmax(,)kkkkyYkkkyYkyYkBpypxypyxpyxPyxpypxypyxpyxpyxpypxyC(12)其中(,)(|,)log(|)(|,)kkyYCpyxpypxy(13)4.EM算法与高斯混合模型在随机系统模型中，假设m是通道m的随机信号生成器的概率密度函数的参数，()pym是选中通道m的概率。记为ma。假设M个随机信号生成器和通道选择随机生成器是相互独立的，从通道m输出的数据x的概率是：(|)mmmapx(14)不考虑通信信息，输出x的概率为：1(|)(|)Mmmmmpxapx(15)其中：m：是第m个通道随机信号生成器的参数。：参数集合1,2...,,mmmMa。观测数据为一批随机产生的输出信号，并且每个输出都是相互独立的，而每个输出来自哪个通道不可测。于是系统模型参数估计问题就变为通过有限的输出样本12{,,...,}Nxxx估计M个通道参数,(1,2,...,)mmamM.应用（12）求解，其中(,)kC可以简化为：111112(,)log()(|,)log((|))(|,)(,)(,)MNMNkkkmnmmmnmnmnkkCapmxpxpmxCC（16）其中：111(,)log()(|,)MNkkmnmnCapmx211(,)log((|))(|,)MNkkmmmnmnCpxpmx这样我们把ma和mp分别放在两项里面，他们不相关，可以独立考虑。在(,)kC中应用约束条件：11Mmma用拉格朗日乘子优化ma得到：111(|,)NkkmnnapmxN上式的含义是，选中m号通道的概率估计1kma是每个观测数据nx来自于m通道的条件概率（根据上一次估值估算）的平均。其中的(|,)knpmx通过下式得出。'''1(|)(|,)(|)kkmmnmnMknmmmmapxpmxapx2(,)kC中的m的优化取决于分布函数的类型，对于(|)mmmpx为高斯分布时，,mmm其中m是分布的均值，m是方差。再经过推导，有：111(|,)NkkmnnapmxN①111(|,)(|,)NknnknmNknnxpmxpmx，②1/212111(|,)||=(|,)NkknnmknmNknnpmxxpmx③m通道参数,mm得更新可以看作是对nx的加权，加权系数(|,)kmpmx可以看成是根据上一次的参数估计k算出来得nx率属于m通道的概率。最后，上面的EM算法可能收敛到局部极大点，因此需要选择多个参数的初始值进行迭代计算，并选择使得(|)px最大的解，(|)px最大的解可由下式算出：12111111(|)(|)(|,)(|)(|)nmNyYNMMMmnmmmmnNMmnmmnpxpypxyapxapx5.EM算法在matlab中的实现利用上面推导出的公式①②③，我们以二个一维的高斯分布（1N,2N）来验证EM算法的性能，首先用二个一维的高斯分布来建立一个高斯混合模型HN。假设：2111(,)NN,2222(,)NN1122HNNN其中1与2为混合系数，且有121，我们要用EM算法估计混合系数和各一维高斯分布的均值和方差22121212(,,,,,)。并将利用EM算法估计出的值与真实值做比较，就可以得到该算法的性能。首先我们取22121212(,,,,,)的真实值为（0.4,0.6,1,2,0.25,0.36）这样我们得到一个混合高斯分布，他的密度函数为HN,然后产生1000个服从HN的分布的观测样本点。接下来要做的就是对这1000个样本点用EM算法进行处理，来估计出一组22121212(,,,,,)的值。在使用EM算法时，要首先给22121212(,,,,,)设定一组初值这里假设初值为1=0.3，2=0.7，10.8，21.8，210.2，220.25Matlab具体程序如下：Y=zeros(1,10000);fori=1:10000ifrand(1)0.3Y(i)=normrnd(2,sqrt(0.36),1,1);elseY(i)=normrnd(1,sqrt(0.25),1,1);endend%高斯混合模型A=[0.30.7];%设置参数的初值M=[0.81.8];%设置均值的初值S=[0.20.25];%设置方差2的初值forn=1:1000forj=1:2a3=0;a4=0;a5=0;fork=1:10000a1=0;fort=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(|,)Nknnpmxa4=a2*Y(k)+a4;%a4对应公式1(|,)Nknnnpmxxa5=a2*(Y(k)-M(j))^2+a5;%a5对应公式121(|,)()NkknnmnpmxxendA(j)=a3/10000;%循环更新系数值M(j)=a4/a3;%循环更新均值值S(j)=a5/a3;%循环更新方差值endend运行程序，查看变量A,M,S的值，与真实值比较一下，就可以得到用EM算法估计的高斯混合模型的性能了。得到参数为A=[0.30630.6937]，M=[1.00932.0013]，S=[0.25580.3632]，而真实值为A=[0.40.6]，M=[12]，S=[0.250.36]，我们可以从中看出用EM算法估计的高斯混合模型的相关参数与真实值是比较接近的。