高斯混合模型的C++实现

原理

GMM将数据的分布通过多个高斯模型进行拟合。GMM是一种聚类算法，每个component就是一个聚类中心。高斯混合模型可以得到每个数据属于每个模型的概率，是一种软聚类算法。这是来自《统计学习方法》中的定义：

过程

高斯混合模型使用EM算法估计模型参数。
1. 初始化模型的个数和每个高斯模型的参数，设定迭代结束条件（迭代次数，误差阈值）
2. 迭代：对于每一个数据，计算在每一个高斯模型中的概率
3. 根据计算得到的概率更新每个模型的参数（均值，方差）
4. 当超过迭代次数或者更新小于阈值时结束迭代。

代码实现

代码参考自网上大神，具体出处忘记了，加以修改和添加注释，如有侵权请联系~

头文件

#ifndef _GMM_H#define _GMM_H#include <vector>#include <cmath>using namespace std;class GMM{public:    void Init(const vector<double> &inputData, const int clustNum = 5, double eps = 0.01, double max_steps = 20);    void train();    int predicate(double x);//预测输入的数据属于哪一类    void print();protected:    int clusterNum;             // 限制    vector<double> means;    vector<double> means_bkp;       // 上一次的迭代数据    vector<double> sigmas;    vector<double> sigmas_bkp;      // 上一次的迭代数据    vector<double> probilities;    vector<double> probilities_bkp;    vector<vector<double>> memberships;     // 存储属于哪一个类别    vector<vector<double>> memberships_bkp;    vector<double> data;    int dataNum;        // 数据数量    double epslon;      // 相差的阈值    double max_steps;   // 迭代次数private:    double gauss(const double x, const double m, const double sigma);};#endif

实现

#include "GMM.h"#include <iostream>#include <fstream>#include <stdlib.h>#include <Windows.h>using namespace std;void GMM::Init(const vector<double> &inputData, const int clustNum, double eps, double max_steps){    /*获取输入数据*/    this->data = inputData;    this->dataNum = data.size();    /*存储最终需要的结果*/    this->clusterNum = clustNum;        // 聚类数量    this->epslon = eps;                 // 阈值    this->max_steps = max_steps;        // 最大的迭代次数    /*保留每一个类别的均值，方差参数，保留上一个的参数*/    this->means.resize(clusterNum);    this->means_bkp.resize(clusterNum);    this->sigmas.resize(clusterNum);    this->sigmas_bkp.resize(clusterNum);    /*保留每一个数据对于每一个类别下的概率，由在这个类别下的概率除以到各个类别的总概率得到*/    this->memberships.resize(clusterNum);    this->memberships_bkp.resize(clusterNum);    for (int i = 0; i < clusterNum; i++)    {        memberships[i].resize(data.size());        memberships_bkp[i].resize(data.size());    }    /*每一个类别的可能性*/    this->probilities.resize(clusterNum);    this->probilities_bkp.resize(clusterNum);    //initialize mixture probabilities 初始化每个类别的参数    for (int i = 0; i < clusterNum; i++)    {        probilities[i] = probilities_bkp[i] = 1.0 / (double)clusterNum;        //init means        means[i] = means_bkp[i] = 255.0*i / (clusterNum);        //init sigma        sigmas[i] = sigmas_bkp[i] = 50;    }}void GMM::train(){    //compute membership probabilities    int i, j, k, m;    double sum = 0, sum2;    int steps = 0;    bool go_on;    do          // 迭代    {        for (k = 0; k < clusterNum; k++)        {            //compute membership probabilities            for (j = 0; j < data.size(); j++)            {                //计算p(k|n)，计算每一个数据在每一个类别的值的加权和                sum = 0;                for (m = 0; m < clusterNum; m++)                {                    sum += probilities[m] * gauss(data[j], means[m], sigmas[m]);                }                //求分子，第j个数据在第k类中的所占的比例                memberships[k][j] = probilities[k] * gauss(data[j], means[k], sigmas[k]) / sum;            }            //求均值            //求条件概率的和，将每个数据在第k类中的概率进行累加            sum = 0;            for (i = 0; i < dataNum; i++)            {                sum += memberships[k][i];            }            //得到每个数据在属于第k类的加权值之和            sum2 = 0;            for (j = 0; j < dataNum; j++)            {                sum2 += memberships[k][j] * data[j];            }            //得到新的均值 由概率加权和除以总概率作为均值            means[k] = sum2 / sum;            //求方差  由到均值的平方的加权和作为新的方差            sum2 = 0;            for (j = 0; j < dataNum; j++)            {                sum2 += memberships[k][j] * (data[j] - means[k])*(data[j] - means[k]);            }            sigmas[k] = sqrt(sum2 / sum);            //求概率            probilities[k] = sum / dataNum;        }//end for k        //check improvement        go_on = false;        for (k = 0; k<clusterNum; k++)        {            if (means[k] - means_bkp[k]>epslon)            {                go_on = true;                break;            }        }        //back up        this->means_bkp = means;        this->sigmas_bkp = sigmas;        this->probilities_bkp = probilities;    } while (go_on&&steps++ < max_steps); //end do while}double GMM::gauss(const double x, const double m, const double sigma){    return 1.0 / (sqrt(2 * 3.1415926)*sigma)*exp(-0.5*(x - m)*(x - m) / (sigma*sigma));}// 预测int GMM::predicate(double x){    double max_p = -100;    int i;    double current_p;    int bestIdx = 0;    for (i = 0; i < clusterNum; i++)    {        current_p = gauss(x, means[i], sigmas[i]);        if (current_p > max_p)        {            max_p = current_p;            bestIdx = i;        }    }    return bestIdx;}void GMM::print(){    int i;    for (i = 0; i < clusterNum; i++)    {        cout << "Mean: " << means[i] << " Sigma: " << sigmas[i] << " Mixture Probability: " << probilities[i] << endl;    }}