本文介绍了混合高斯聚类算法。首先介绍了混合高斯的类表示是一个高斯模型，相似性度量定义为服从类参数为高斯分布，其是一种典型的基于模型的密度聚类算法。然后介绍了混合高斯模型假设类间服从伯努利分布，类内服从高斯分布，结合最大似然函数给出了混合高斯模型的目标函数。最后介绍了混合高斯模型的EM求解流程。
作者 | 文杰
编辑 | yuquanle
模型聚类
高斯混合
高斯混合的类表示是一个高斯模型，相似性度量定义为服从类高斯分布的概率（Kmeans的相似度量是距离度量），所以高斯混合聚类也可以看作是有参的密度聚类。
高斯混合假设类之间服从伯努利分布，样本在某一类下服从高斯分布，也就是说每个样本独立服从多元高斯分布。为了使得所有样本的概率最大化，即最大化对数似然函数：
也就是说假设类之间服从一个伯努利分布：
样本在类下的条件概率服从高斯分布：
那么样本和类标签的联合分布为：
以上，当已知时，即标签信息已知的话，类似于高斯判别分析（当然，高斯判别分析中多个高斯分布之间具有相同的协方差），对应的且只属于一类（类标已知），那么上式有：
最大似然估计有参数：
可以看出为每一类样本所占的比例，为该类下样本的均值，为该类下样本的协方差。
考虑到高斯混合模型中的类划分是概率划分，表示第个样本属于第类的概率。所以，高斯混合模型的所有参数都需要乘上类的划分概率。
高斯混合模型流程
1）初始化参数隐类别数，模型参数
2）采用EM算法，先假设参数，期望最大化，然后更新样本的划分概率更新参数
a）E-step：对所有样本，根据参数划分每个样本类概率：
b）M-step：根据划分后的类概率更新参数
这里可以看出，当是已知的，即类标签已知，则直接进行参数估计等价于高斯判别分析，当是硬划分，同Kmenas又是一致的。
代码实战
代码语言：javascript复制Matrix GMM::E_step(const Matrix &x)
{
Matrix w;
w.initMatrix(K, x._row, 0);
Matrix wNorm;
wNorm.initMatrix(1, x._row, 0);
Matrix xOne;
double p = 0;
for(size_t i = 0; i < x._row; i )
{
for(int k = 0; k < w._row; k )
{
xOne = x.getOneRow(i);
p = Phi[k] * 1.0 / (pow(PI,double(K/2.0)) * COV2(Sigma[k])) * exp(-0.5 *(((xOne-Mu[k]) * Sigma[k].niMatrix()) * (xOne-Mu[k]).transposeMatrix())._data[0][0]);
w._data[k][i] = p;
wNorm._data[0][i] = p;
}
for(int k = 0; k < w._row; k )
{
w._data[k][i] /= wNorm._data[0][i];
}
}
return w;
}
int GMM::M_step(const Matrix &w, const Matrix &x)
{
Phi.clear();
Mu.clear();
Sigma.clear();
Matrix mu;
mu.initMatrix(1, x._col, 0);
Matrix sigma;
sigma.initMatrix(x._col, x._col, 0);
for(size_t k = 0; k < w._row; k )
{
Phi.push_back(0);
for(size_t j = 1; j < x._col; j )
mu._data[0][j] = 0;
for(size_t i = 0; i < sigma._row; i )
for(size_t j = 1; j < sigma._col; j )
sigma._data[i][j] = 0;
for(size_t i = 0; i < x._row; i )
{
Phi[k] = w._data[k][i];//计算phi
for(size_t j = 0; j < x._col; j )
mu._data[0][j] = w._data[k][i] * x._data[i][j];//计算mu
}
Phi[k] /= x._row;//更新phi
for(size_t j = 0; j < x._col; j )
mu._data[0][j] /= (Phi[k] * x._row);
Mu.push_back(mu);//更新mu
for(size_t i = 0; i < x._row; i )
{
for(size_t j = 0; j < x._col; j )
{
for(size_t j2 = 0; j2 < x._col; j2 )
{
sigma._data[j][j2] = w._data[k][i] * (x._data[i][j]-mu._data[0][j]) * (x._data[i][j2]-mu._data[0][j2]);
}
}
}
for(size_t i = 0; i < sigma._row; i )
for(size_t j = 0; j < sigma._col; j )
sigma._data[i][j] /= (Phi[k] * x._row);
Sigma.push_back(sigma);//更新Sigma
}
return 0;
}The End
