c语言em算法实现,统计学习方法ｃ++实现之八 EM算法与高斯混合模型

EM算法与高斯混合模型

前言

EM算法是一种用于含有隐变量的概率模型参数的极大似然估计的迭代算法。如果给定的概率模型的变量都是可观测变量，那么给定观测数据后，就可以根据极大似然估计来求出模型的参数，比如我们假设抛硬币的正面朝上的概率为p(相当于我们假设了概率模型)，然后根据n次抛硬币的结果就可以估计出p的值，这种概率模型没有隐变量，而书中的三个硬币的问题(先抛A然后根据A的结果决定继续抛B还是C)，这种问题中A的结果就是隐变量，我们只有最后一个硬币的结果，其中的隐变量无法观测，所以这种无法直接根据观测数据估计概率模型的参数，这时就需要对隐变量进行估计，进而得到概率模型的参数，这里要注意，概率模型是已知的(已经假定好了)，包括隐变量的模型也是假设好的，只是具体的参数未知，这时候就需要用EM算法求解未知参数，这里我用EM算法估计了高斯混合模型的参数，并用高斯混合模型实现了聚类，代码地址。

EM算法

EM算法中文名称是期望极大算法，EM是expectation maximization的缩写，从名字就可以窥视算法的核心，求期望，求极大。求谁的期望？求似然函数对隐变量的期望，所以，首先必须确定隐变量是什么。其次，对谁求极大？当然是求出概率模型的参数使得上一步的期望最大。算法如下：

输入：观测变量数据Y，隐变量数据Z(这里也是知道的？其实这里我的理解是，这里不是已知的，但是却是可以根据假设的隐变量的参数得到的)，联合分布\(P(Y,Z|\theta)\), 条件分布\(P(Z|Y,\theta)\)

输出：模型参数\(\theta\)

E步:记\(\theta^{(i)}\)为第i次迭代参数\(\theta\)的估计值，i+1次迭代的E步， 计算

\(Q(\theta, \theta^{(i)})=E_Z[logP(Y,Z|\theta)|Y, \theta^{(i)}]=\sum_{Z}logP(Y,Z|\theta)P(Z|Y,\theta^{(i)})\)

M步，求使\(Q(\theta, \theta^{(i)})\)极大的\(\theta\),作为下一次迭代的\(\theta^{(i+1)}\)

重复2,3直到收敛

可以看出最重要的在于求\(Q(\theta, \theta^{(i)})\)，那么为什么每一次迭代最大化\(\sum_{Z}logP(Y,Z|\theta)P(Z|Y,\theta^{(i)})\)就能使观测数据的似然函数最大(这是我们的最终目的)？这里书上有证明，很详细，就不赘述了。先看一下我们要最大化的极大似然函数，然后这里主要引用西瓜书中的解释来从理解EM算法：

\(L(\theta)=logP(Y|\theta)=log(\sum_ZP(Y,Z|\theta))=log(\sum_ZP(Y|Z,\theta)P(Z|\theta))\)

在迭代过程中，若参数\(\theta\)已知，则可以根据训练数据推断出最优隐变量Z的值(E步)；反之，若Z的值已知，则可以方便的对参数\(\theta\)做极大似然估计(M步)。

可以看出就是一种互相计算，一起提升的过程。

这部分如果推导看不懂了，可以结合下面的二维高斯混合模型的EM算法来理解。

c++实现

这里我使用EM算法来估计高斯混合模型的参数来进行聚类，高斯混合模型还有一个很大的作用是进行前景提取，这里仅仅用二维混合高斯模型进行聚类。

代码结构

关键代码

这里实现起来其实没什么难点，难点在于推导参数的更新公式，详情参考西瓜书p206。

void GMM::EMAlgorithm(vector &alphaOld, vector>> &sigmaOld,

vector> &muOld) {

// compute gamma

for (int i = 0; i < trainDataF.size(); ++i) {

double probSum = 0;

for (int l = 0; l < alpha.size(); ++l) {

double gas = gaussian(muOld[l], sigmaOld[l], trainDataF[i]);

probSum += alphaOld[l] * gas;

}

for (int k = 0; k < alpha.size(); ++k) {

double gas = gaussian(muOld[k], sigmaOld[k], trainDataF[i]);

gamma[i][k] = alphaOld[k] * gas / probSum;

}

}

// update mu, sigma, alpha

for (int k = 0; k < alpha.size(); ++k) {

vector muNew;

vector> sigmaNew;

double alphaNew;

vector muNumerator;

double sumGamma = 0.0;

for (int i = 0; i < trainDataF.size(); ++i) {

sumGamma += gamma[i][k];

if (i==0) {

muNumerator = gamma[i][k] * trainDataF[i];

}

else {

muNumerator = muNumerator + gamma[i][k] * trainDataF[i];

}

}

muNew = muNumerator / sumGamma;

for (int i = 0; i < trainDataF.size(); ++i) {

if (i==0) {

auto temp1 = gamma[i][k]/ sumGamma * (trainDataF[i] - muNew);

auto temp2 = trainDataF[i] - muNew;

sigmaNew = vecMulVecToMat(temp1, temp2);

}

else {

auto temp1 = gamma[i][k] / sumGamma * (trainDataF[i] - muNew);

auto temp2 = trainDataF[i] - muNew;

sigmaNew = sigmaNew + vecMulVecToMat(temp1, temp2);

}

}

alphaNew = sumGamma / trainDataF.size();

mu[k] = muNew;

sigma[k] = sigmaNew;

alpha[k] = alphaNew;

}

}

总结

前面的代码一直用vector来实现向量，但是这里用到了矩阵，矩阵的相关计算都添加的计算函数。最正规的应该是写个类，实现矩阵运算，但是这里偷懒了，以后写代码一定要考虑周到，这样添添补补的太低效了。

《统计学习方法》笔记九 EM算法及其推广

本系列笔记内容参考来源为李航 EM算法是一种迭代算法,用于含有隐变量的概率模型参数的极大似然估计或极大后验概率估计.迭代由 (1)E步:求期望 (2)M步:求极大 组成,称 ...

EM 算法求解高斯混合模型python实现

注:本文是对EM算法的一个简单总结. 1. 什么是EM算法? 引用书上的话: 概率模型有时既含有观测变量,又含有隐变量或者潜在变量.如果概率模型的变量都是观测变量,可以直接 ...

斯坦福大学机器学习，EM算法求解高斯混合模型

斯坦福大学机器学习,EM算法求解高斯混合模型.一种高斯混合模型算法的改进方法---将聚类算法与传统高斯混合模型结合起来的建模方法, 并同时提出的运用距离加权的矢量量化方法获取初始值,并采用衡量相似度的 ...

机器学习第三课(EM算法和高斯混合模型)

极大似然估计,只是一种概率论在统计学的应用,它是参数估计的方法之一.说的是已知某个随机样本满足某种概率分布,但是其中具体的参数不清楚,参数估计就是通过若干次试验,观察其结果,利用结果推出参数的大概值. ...

机器学习算法总结&lpar;六&rpar;——EM算法与高斯混合模型

极大似然估计是利用已知的样本结果,去反推最有可能(最大概率)导致这样结果的参数值,也就是在给定的观测变量下去估计参数值.然而现实中可能存在这样的问题,除了观测变量之外,还存在着未知的隐变量,因为变量未 ...

EM算法求高斯混合模型參数预计——Python实现

EM算法一般表述:       当有部分数据缺失或者无法观察到时,EM算法提供了一个高效的迭代程序用来计算这些数据的最大似然预计.在每一步迭代分为两个步骤:期望(Expectation)步骤和最大化( ...

EM算法和高斯混合模型GMM介绍

EM算法 EM算法主要用于求概率密度函数参数的最大似然估计,将问题$\arg \max _{\theta_{1}} \sum_{i=1}^{n} \ln p\left(x_{i} | \theta_{ ...

统计学习方法笔记--EM算法--三硬币例子补充

本文,意在说明第九章EM算法的三硬币例子,公式(9.5-9.6如何而来) 下面是(公式9.5-9.8)的说明, 本人水平有限,怀着分享学习的态度发表此文,欢迎大家批评,交流 ...

学习笔记——EM算法

EM算法是一种迭代算法,用于含有隐变量(hidden variable)的概率模型参数的极大似然估计,或极大后验概率估计.EM算法的每次迭代由两步组成:E步,求期望(expectation):M步,求 ...

随机推荐

eclipse jetty启动内存溢出

一.eclipse jetty启动内存溢出, 异常信息 Exception in thread "ConfigClientWorker-Default" java.lang.Out ...

java遍历map方法

java 代码: import java.util.HashMap; import java.util.Iterator; import java.util.Map; public class Map ...

PHP字符编码问题-总结

今天在网上看到一个人的对于php开发中字符编码的总结,感觉不错,摘录如下: 一,php编码转换        1.通过iconv()函数实现编码转换                语法:iconv(s ...

&period;net EF 事物 订单流水号的生成 &lpar;一&rpar;

首先需要 添加 System.Transactions 程序集 数据表: create table SalesOrder ( ID ,) primary key not null, OrderNo ) ...

formValidator

formValidator输入验证.异步验证实例 + licenseImage验证码插件实例应用   实例技术:springmvc 实现功能:完整用户登录流程.输入信息规则校验.验证码异步校验. 功能 ...

二叉树的序列化和反序列化(Java)

请实现两个函数,分别用来序列化和反序列化二叉树 序列化就是将二叉树以字符串输出,反序列化:根据自己输出的字符串,构建二叉树. 这里先序遍历输出,且为了方便反序列化,各个节点","隔 ...

Web设计快速入门

在基本顺利完成功能的基础上,就需要考虑美观的问题了,在眼球经济的当下,一个面向用户的产品,如果没有好的UI,那么它就是不合格的.这部分内容算是初出茅庐,会持续更新. "一个人的外貌决定我是否 ...

iOS边练边学--自定义等高的cell

一.storyboard自定义cell <1>创建一个继承自UITableViewCell的子类,比如ChaosDealCell <2>在storyboard中 <2.1 ...

gulp &amp&semi; webpack整合

为什么需要前端工程化? 前端工程化的意义在于让前端这个行业由野蛮时代进化为正规军时代,近年来很多相关的工具和概念诞生.好奇心日报在进行前端工程化的过程中,主要的挑战在于解决如下问题:✦ 如何管理多个项 ...

LeetCode第&lbrack;11&rsqb;题&lpar;Java&rpar;：Container With Most Water (数组容器盛水)——Medium

题目难度:Medium Given n non-negative integers a1, a2, ..., an, where each represents a point at coordina ...