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

程序员文章站 2022-07-14 14:38:28

...

EM算法/Expectation Maximization

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

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

考虑从高斯混合模型画出的d维欧几里得空间中的N个特征向量{ $x_{1} ， x_{2} ， . . . ， x_{N}$ }的集合：

$p (x; a_{k}; S_{k}; π_{k}) = \sum_{k = 1}^{m} π_{k} p_{k} (x), π_{k} ⩾ 0, \sum_{k = 1}^{m} π_{k} = 1$

$p_{k} (x) = φ (x; a_{k}, S_{k}) = \frac{1}{(2 π)^{d / 2} | S_{k} |^{1 / 2}} e x p {- \frac{1}{2} (x - a_{k})^{T} S_{k}^{- 1} (x - a_{k})}$

其中m是高斯混合模型的数量， $p_{k}$ 是具有均值 $a_{k}$ 和协方差矩阵 $S_{k}$ 的正态分布密度， $π_{k}$ 是第k个高斯混合模型的权重。给定高斯混合模型个数M和样本 $x_{i} ， i = 1... N$ ，算法找到所有高斯混合模型参数的最大似然估计（MLE），即 $a_{k} ， S_{k}$ 和 $π_{k}$ ：

$L (x, θ) = l o g p (x, θ) = \sum_{i = 1}^{n} l o g (\sum_{k = 1}^{m} π_{k} p_{k} (x)) \to m a x_{θ \in Θ,}$

$Θ =$
${(a_{k}, S_{k}, π_{k}) : a_{k} \in R^{d}, S_{k} = S_{k}^{T} > 0, S_{k} \in R^{d \times d}, π_{k} \geq 0, \sum_{k = 1}^{m} π_{k} = 1}$

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

$a_{k i} = \frac{π_{k} φ (x; a_{k}, S_{k})}{\sum_{j = 1}^{m} π_{j} φ (x; a_{j}, S_{j})}$

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

$π_{k} = \frac{1}{N} \sum_{i = 1}^{N} a_{k i}$ ;

$π_{k} = \frac{\sum_{i = 1}^{N} a_{k i} x_{i}}{\sum_{i = 1}^{N} a_{k i}}$ ;

$S_{k} = \frac{\sum_{i = 1}^{N} a_{k i} (x_{i} - a_{k}) (x_{i} - a_{k})^{T}}{\sum_{i = 1}^{N} a_{k i}}$

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

EM算法的一个主要问题是需要估计大量参数。大多数参数存在于协方差矩阵中，这些矩阵大小为d×d，其中d是特征空间维度。但在许多实际问题中，协方差矩阵接近于对角线或者甚至接近 $μ_{k} * I$ ，其中 $I$ 是单位矩阵， $μ_{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;
}

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

点坐标分类

从文件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;
}

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

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

上一篇： cutadapt 截取fastq导致fq文件报错

下一篇： opencv3/C++ 机器学习-逻辑回归/Logistic Regression

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

EM算法/Expectation Maximization

OpenCV EM类

相关函数

应用示例

图像分割

点坐标分类

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

opencv3/C++ 机器学习-提升算法/Boosting