C++元编程——EM算法求解高斯混合模型

腾昵猫

已于 2023-05-24 14:23:15 修改

阅读量378

点赞数

文章标签： c++ 算法数据结构

于 2023-05-16 15:20:22 首次发布

本文链接：https://blog.csdn.net/Dr_Jack/article/details/130706454

版权

EM算法跟K-means差不多，步骤就是：1）首先随机固定个数类的均值和方差矩阵；2）将数据集按照目前的分类计算概率，分派到每个类上；3）对各类中的数据计算其均值和方差；4）误差达标则退出，未达标则返回第2步继续执行。下面是实现代码，模拟的是多高斯分布使用EM算法进行优化。

其中第2）3）步为核心。这里有几个公式要讲一下。

$\gamma_{ik}=\frac{\pi_{k}N(x_{i};\mu_{k}, \sigma_{k})}{\sum_{j=0}^{K}\pi_{j}N(x_{i};\mu_{j}, \sigma_{j})}$

这个值代表的是在样本xi的情况下选中第k类的概率，其中pi表示的是第k类的先验概率；所以分子表示的是在k发生且x为k类时的联合概率，分母表示的是在xi发生的边缘概率，所以，根据贝叶斯公式可得这个公式表示的是在xi的情况下选中k的概率有多大，这个就是说为的E-step，就是根据目前的参数情况判断各个样本的概率。

根据这个可以得到第k类的均值和方差(这就就是M-step，目的就是调整概率参数使得在样本集的情况下联合的概率最大化)。

$\mu_{k}=\sum_{i=0}^{N}\frac{\gamma_{ik}*x_{i}}{\sum_{j=0}^{N}\gamma_{jk}}$

$\sigma_{k}=\sum_{i=0}^{N}\frac{\gamma_{ik}*(x_{i}-\mu_{k})\cdot(x_{i}-\mu_{k})'}{\sum_{j=0}^{N}\gamma_{jk}}$

$\pi_{k}=\frac{\sum_{j=0}^{N}\gamma_{jk}}{N}$

最终的停止条件就是均值和方差变化足够小就停止。

#include <vector>
#include <functional>

#include "mat.hpp"
#include "base_function.hpp"

template<int dim_size>
double poss(const mat<dim_size, 1, double>& x, const mat<dim_size, 1, double>& u, const mat<dim_size, dim_size, double>& sigma)
{
	mat<dim_size, dim_size, double> E_1 = inverse(sigma);
	mat<dim_size, 1, double> x_u = x - u;
	double _E_1_2 = (sqrt(det(sigma))*pow(2.*3.1415926535897932384626, dim_size / 2.));
	return exp(x_u.t().dot(E_1.dot(x_u))*-0.5)[0] / _E_1_2;
}

template<int dim_size>
struct em_class 
{
	std::vector<int>						vec;
	mat<dim_size, 1, double>				u;
	mat<dim_size, dim_size, double>			sigma;
	double									p;
};

template<int dim_size>
std::function<double(const mat<dim_size, 1, double>&)> gen_gama(const std::vector<mat<dim_size, 1, double> >& vec_x, const std::vector<em_class<dim_size> >& vec_cls, const int& k)
{
	return [vec_x, vec_cls, k](const mat<dim_size, 1, double>& x) 
	{
		double poss_all = 0.;
		for (int j = 0; j < vec_cls.size(); ++j) 
		{
			poss_all += (vec_cls[j].p * poss(x, vec_cls[j].u, vec_cls[j].sigma));
		}
		return vec_cls[k].p * poss(x, vec_cls[k].u, vec_cls[k].sigma) / poss_all;
	};
}

template<int dim_size>
bool update_class(const std::vector<mat<dim_size, 1, double> >& vec_x, std::vector<em_class<dim_size> >& vec_cls)
{
	double delta_u_max = 0.;
	double delta_sigma_max = 0.;
	double N = vec_x.size();
	std::vector<em_class<dim_size> > vec_cls_new(vec_cls.size());
	for (int k = 0; k < vec_cls.size(); ++k)
	{
		mat<dim_size, 1, double> uk = 0.;
		mat<dim_size, dim_size, double> sigmak = 0.;
		auto gama = gen_gama(vec_x, vec_cls, k);
		double Nk = 0.;
		for (int n = 0; n < vec_x.size(); ++n)
		{
			auto xn = vec_x[n];
			auto r_nk = gama(xn);					// xn对k的概率
			uk = uk + r_nk * xn;
			Nk += r_nk;							// 所有样本集中的数据对于类型k的概率和
		}
		uk = uk / Nk;
		for (int n = 0; n < vec_x.size(); ++n)
		{
			auto xn = vec_x[n];
			sigmak = sigmak + gama(xn) * (vec_x[n] - uk).dot((vec_x[n] - uk).t());
		}
		double pk = Nk / N;
		
		sigmak = sigmak / Nk;
		auto delta_u = (vec_cls[k].u - uk).max_abs();
		auto delta_sigma = (vec_cls[k].sigma - sigmak).max_abs();
		delta_u_max = (delta_u > delta_u_max ? delta_u : delta_u_max);
		delta_sigma_max = (delta_sigma > delta_sigma_max ? delta_sigma : delta_sigma_max);
		vec_cls_new[k].u = uk;
		vec_cls_new[k].sigma = sigmak;
		vec_cls_new[k].p = pk;
		//vec_cls[k].vec.clear();
	}
	vec_cls.swap(vec_cls_new);
	if (delta_u_max < 0.0001 && delta_sigma_max < 0.001)
		return true;
	return false;
}

template<int dim_size>
void classify(const std::vector<mat<dim_size, 1, double> >& vec_x, std::vector<struct em_class<dim_size> >& vec_cls) 
{
	for (int lp = 0; lp < 500; ++lp) 
	{
		if (update_class(vec_x, vec_cls)) 
		{
			printf("%d			exit\r\n", lp);
			break;
		}
	}
}

#include <random>
#include <weight_initilizer.hpp>

int main(int argc, char** argv) 
{
	constexpr int dim_size = 2;
	using mat_type = mat<dim_size, 1, double>;
	std::vector<mat_type> vec_x;

	/* 生成数据 */
	std::default_random_engine ge;
	std::normal_distribution <double> ud(0., 5.);
	const int set_size = 1000;
	for (int i = 0; i < set_size; ++i)
	{
		mat_type mt;
		for (int i = 0; i < mt.size(); ++i) 
		{
			mt.get(i, 0) = ud(ge);
		}
		vec_x.push_back(mt);
	}
	std::normal_distribution <double> ud1(20., 8.);
	for (int i = 0; i < set_size; ++i)
	{
		mat_type mt;
		for (int i = 0; i < dim_size; ++i)
		{
			mt.get(i, 0) = ud1(ge);
		}
		vec_x.push_back(mt);
	}
	std::normal_distribution <double> ud2(30., 4.);
	for (int i = 0; i < set_size; ++i)
	{
		mat_type mt;
		for (int i = 0; i < dim_size; ++i)
		{
			mt.get(i, 0) = ud2(ge);
		}
		vec_x.push_back(mt);
	}
	/* 初始化类别 */
	std::random_device rd;
	std::mt19937 rng(rd());
	std::shuffle(vec_x.begin(), vec_x.end(), rng);

	std::vector<struct em_class<dim_size>> vec_cls(3);
	class def_init;
	for (int lp = 0; lp < vec_cls.size(); ++lp)
	{
		weight_initilizer<def_init>::cal(vec_cls[lp].u, 0., 50.);
		weight_initilizer<def_init>::cal(vec_cls[lp].sigma, 0., 50.);
		vec_cls[lp].sigma = vec_cls[lp].sigma + vec_cls[lp].sigma.t();
		vec_cls[lp].p = (1./vec_cls.size());
	}
	/* 执行分类 */
	classify(vec_x, vec_cls);
	/* 展示分类结果 */
	for (int lp = 0; lp < vec_cls.size(); ++lp) 
	{
		printf("---------- %d -----------\r\n", lp);
		vec_cls[lp].u.print();
		vec_cls[lp].sigma.print();
	}
	return 0;
}

结果如下：

20次迭代就收敛了，结果也正确，识别出3类的均值和方差差别不大。

参考资料：

EM算法在高斯混合模型中的应用(非常简单，完全按照EM算法的节奏)_watermelon12138的博客-CSDN博客

腾昵猫

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
打赏
2
评论
C++元编程——EM算法求解高斯混合模型

EM算法跟K-means差不多，步骤就是：1）首先随机固定个数类的均值和方差矩阵；2）将数据集按照目前的分类计算概率，分派到每个类上；3）对各类中的数据计算其均值和方差；4）误差达标则退出，未达标则返回第2步继续执行。下面是实现代码，模拟的是多高斯分布使用EM算法进行优化。20次迭代就收敛了，结果也正确，识别出3类的均值和方差差别不大。
复制链接

扫一扫