opencv3/C++ 机器学习-EM算法/Expectation Maximization

EM算法/Expectation Maximization

EM算法包含两步：E，求期望（Expectation），利用概率模型参数的现有估计值，计算隐藏变量的期望；M，求极大（Maximization），利用E 步上求得的隐藏变量的期望，对参数模型进行最大似然估计。所得参数估计值用于下个E步的计算，重复至收敛。

期望最大化/EM算法以具有一定数量混合物的高斯混合分布的形式估计多变量概率密度函数的参数。

考虑从高斯混合模型画出的d维欧几里得空间中的N个特征向量{ x1，x2，...，xNx1，x2，...，xN\ x_1，x_2，...，x_N}的集合：

p(x;ak;Sk;πk)=∑mk=1πkpk(x),πk⩾0,∑mk=1πk=1p(x;ak;Sk;πk)=∑k=1mπkpk(x),πk⩾0,∑k=1mπk=1\ p(x;a_k;S_k;\pi_k)=\sum_{k=1}^m\pi_kp_k(x), \pi_k\geqslant 0, \sum_{k=1}^{m}\pi_k=1

pk(x)=φ(x;ak,Sk)=1(2π)d/2|Sk|1/2exp{−12(x−ak)TS−1k(x−ak)}pk(x)=φ(x;ak,Sk)=1(2π)d/2|Sk|1/2exp{−12(x−ak)TSk−1(x−ak)}\ p_k(x)=\varphi (x;a_k,S_k)=\frac{1}{(2\pi)^{d/2}|S_k|^{1/2}}exp\left \{ -\frac{1}{2}(x-a_k)^TS_k^{-1}(x-a_k)\right \}

其中m是高斯混合模型的数量， pkpk\ p_k是具有均值 akak\ a_k和协方差矩阵 SkSk\ S_k的正态分布密度， πkπk\ π_k是第k个高斯混合模型的权重。给定高斯混合模型个数M和样本 xi，i=1...Nxi，i=1...N\ x_i，i = 1...N，算法找到所有高斯混合模型参数的最大似然估计（MLE），即 ak，Skak，Sk\ a_k，S_k和 πkπk\ π_k：

L(x,θ)=logp(x,θ)=∑ni=1log(∑mk=1πkpk(x))→maxθ∈Θ,L(x,θ)=logp(x,θ)=∑i=1nlog(∑k=1mπkpk(x))→maxθ∈Θ, L(x,θ)=logp(x,θ)=\sum_{i=1}^nlog(\sum_{k=1}^m\pi_kp_k(x))\rightarrow max_{θ\in \Theta ,}

Θ=Θ= \Theta =
{(ak,Sk,πk):ak∈Rd,Sk=STk>0,Sk∈Rd×d,πk≥0,∑mk=1πk=1}{(ak,Sk,πk):ak∈Rd,Sk=SkT>0,Sk∈Rd×d,πk≥0,∑k=1mπk=1} \left \{ (a_k,S_k,\pi_k):a_k\in R^d, S_k=S_k^T>0,S_k\in R^{d\times d},\pi_k\geq 0, \sum_{k=1}^m\pi_k=1\right \}

EM算法是一个迭代过程。每次迭代包括两个步骤。在第一步E步即预期步骤中，使用当前可用的混合参数估计值，可以找出样本i属于混合模型k的概率 pikpik\ p_{ik}（在下面的公式中表示为 αikαik\ α_{ik}）：

aki=πkφ(x;ak,Sk)∑mj=1πjφ(x;aj,Sj)aki=πkφ(x;ak,Sk)∑j=1mπjφ(x;aj,Sj) a_{ki}=\frac{\pi_k\varphi(x;a_k,S_k)}{\sum_{j=1}^m\pi_j\varphi(x;a_j,S_j)}

在第二步M步即最大化步骤中，使用计算出的概率对高斯混合模型的参数估计值进行细化：

πk=1N∑Ni=1akiπk=1N∑i=1Naki \pi_k=\frac{1}{N}\sum_{i=1}^Na_{ki};

πk=∑Ni=1akixi∑Ni=1akiπk=∑i=1Nakixi∑i=1Naki \pi_k=\frac{\sum_{i=1}^Na_{ki}x_i}{\sum_{i=1}^Na_{ki}} ;

Sk=∑Ni=1aki(xi−ak)(xi−ak)T∑Ni=1akiSk=∑i=1Naki(xi−ak)(xi−ak)T∑i=1Naki S_k=\frac{\sum_{i=1}^Na_{ki}(x_i-a_k)(x_i-a_k)^T}{\sum_{i=1}^Na_{ki}}

或者，当提供 pikpik\ p_{ik}的初始值时，该算法可以从M步开始。当 pikpik\ p_{ik}未知时的另一种选择是使用更简单的聚类算法对输入采样进行预先聚类，从而获得初始的 pikpik\ p_{ik}（通常用k-means算法实现）。

EM算法的一个主要问题是需要估计大量参数。大多数参数存在于协方差矩阵中，这些矩阵大小为d×d，其中d是特征空间维度。但在许多实际问题中，协方差矩阵接近于对角线或者甚至接近μk∗Iμk∗I μ_k* I，其中III是单位矩阵，μk" role="presentation" style="position: relative;">μkμkμ_k是混合相关的“比例”参数。因此，一个健壮的计算方案是对协方差矩阵加较强的约束，然后用估计的参数作为较少约束优化问题的输入（通常对角协方差矩阵已经足够了）。

OpenCV EM类

应用示例

图像分割

使用EM算法对图像进行分割。

#include <opencv2/opencv.hpp>
#include <iostream>
using namespace std;
using namespace cv;
using namespace cv::ml;int main()
{Vec3b colors[] ={Vec3b(0, 0, 255), Vec3b(0, 255, 0), Vec3b(255, 100, 100), Vec3b(255, 0, 255)};Mat data, labels, src, dst;src = imread("E:/image/image/red.jpg", 1);resize(src, src, Size(src.cols/1.5,src.rows/1.5));if(src.empty()){printf("can not load image \n");return -1;}src.copyTo(dst);for (int i = 0; i < src.rows; i++)for (int j = 0; j < src.cols; j++){Vec3b point = src.at<Vec3b>(i, j);Mat tmp = (Mat_<float>(1, 3) << point[0], point[1], point[2]);data.push_back(tmp);}Ptr<EM> model = EM::create();model->setClustersNumber(4); //类个数model->setCovarianceMatrixType(EM::COV_MAT_SPHERICAL);model->setTermCriteria(TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, 300, 0.1));model->trainEM(data, noArray(), labels, noArray());int n = 0;//显示结果，不同的类别用不同的颜色for (int i = 0; i < dst.rows; i++)for (int j = 0; j < dst.cols; j++){int index = labels.at<int>(n);dst.at<Vec3b>(i, j) = colors[index];n++;}imshow("src", src);imshow("dst", dst);waitKey(0);return 0;
}

点坐标分类

从文件points.txt中读取点坐标以及对应的分类，然后使用EM算法对点所在区域进行划分。

#include <opencv2/opencv.hpp>
#include "opencv2/ml.hpp"
#include <iostream>
#include <fstream> using namespace std;
using namespace cv;
using namespace cv::ml;//EM算法
int main()
{Mat src, dst;vector<Point>  trainedPoints;vector<int>    trainedPointsMarkers;//读取文件中的点坐标FILE *fp;int flge = 0;int fpoint,flabel;fp = fopen("E:\\points.txt", "r+");if (fp == NULL){printf("Cannot open the file!\n");exit(0);}Point point;while (!feof(fp)){   fscanf(fp, "%d", &fpoint);if (feof(fp)) break;//依次为横坐标、纵坐标、分类if ((flge%3==0? point.x = fpoint: flge%3==1? point.y = fpoint:flge%3==2? flabel = fpoint : -1)<0) return -1;if (flge%3==2){trainedPoints.push_back(point);trainedPointsMarkers.push_back(flabel);}flge++;}vector<Vec3b>  colors(4);colors[0] = Vec3b(0, 255, 0);colors[1] = Vec3b(0, 0, 255);colors[2] = Vec3b(0, 255, 255);colors[3] = Vec3b(255, 0, 0);src.create( 480, 640, CV_8UC3 );src = Scalar::all(0);// 绘制点for( size_t i = 0; i < trainedPoints.size(); i++ ){Scalar c = colors[trainedPointsMarkers[i]];circle( src, trainedPoints[i], 3, c, -1 );}src.copyTo(dst);imshow( "points", src );Mat samples;Mat(trainedPoints).reshape(1, (int)trainedPoints.size()).convertTo(samples, CV_32F);int nmodels = (int)colors.size();vector<Ptr<EM> > em_models(nmodels);Mat modelSamples;for( int i = 0; i < nmodels; i++ ){modelSamples.release();for( int j = 0; j < samples.rows; j++ ){if( trainedPointsMarkers[j] == i )modelSamples.push_back(samples.row(j));}// 训练模型if( !modelSamples.empty() ){const int componentCount = 5;Ptr<EM> em = EM::create();//高斯混合模型中混合成分的数量em->setClustersNumber(componentCount);//协方差矩阵的类型。em->setCovarianceMatrixType(EM::COV_MAT_DIAGONAL);//训练模型em->trainEM(modelSamples, noArray(), noArray(), noArray());em_models[i] = em;}}Mat testSample(1, 2, CV_32FC1 );Mat logLikelihoods(1, nmodels, CV_64FC1, Scalar(-DBL_MAX));for( int y = 0; y < src.rows; y += 3 ){for( int x = 0; x < src.cols; x += 3 ){testSample.at<float>(0) = (float)x;testSample.at<float>(1) = (float)y;for( int i = 0; i < nmodels; i++ ){if( !em_models[i].empty() )logLikelihoods.at<double>(i) = em_models[i]->predict2(testSample, noArray())[0];}Point maxLoc;minMaxLoc(logLikelihoods, 0, 0, 0, &maxLoc);dst.at<Vec3b>(y, x) = colors[maxLoc.x];}}imshow( "EM", dst );waitKey();return 0;
}

文件points.txt中的内容为：
（依次为横坐标、纵坐标、分类）