GMM高斯混合模型Python代码实现

前言

上一篇GMM高斯混合模型原理推导(二)我们已经对高斯混合模型进行了推导，接下来就对其进行代码实现。

算法流程

①随机初始化模型参数$p^t,{\mu }^t,{\Sigma }^t$。

②计算出$P(z_i=Ck|x_i,\theta )$

③依据公式计算出$p^{t+1},{\mu }^{t+1},{\Sigma }^{t+1}$

④计算$p^{t+1},{\mu }^{t+1},{\Sigma }^{t+1}$和$p^t,{\mu }^t,{\Sigma }^t$的差值，如果差值小于$ϵ$(自己设定的值)。如果小于则说明变化太小，证明收敛，结束算法。否则循环②，③步骤

代码复现

值得注意的是，EM算法本身对随机初始化的参数是极其敏感的，不同的初始值很有可能会造成不同的结果，呃，貌似基本上优化算法都有这个问题，都很容易陷入局部最优。而局部最优点的个数往往是随着维度而指数增长的。因此，对于GMM模型的初始化参数。请遵循以下原则：

不要把所有同类型的值都初始化为一样。比如对于隐变量，$p_1=1$了，$p_2$就不要再初始化为1。就算你隐变量都初始化为一样了，均值也尽量不要一样。第一类均值和第二类均值尽量不一样。那假如隐变量和均值都一样了，请不要再把协方差也设置为一样，否则天王老子来了都没用。具体原因请看推出来的公式，自己代入公式试试就知道了。

训练过程及其结果

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
plt.ion()
class GMM():
    def __init__(self):
        pass
    def train(self,x,K):
        '''
        :param x: 观测变量
        :param K: 表示聚成几类
        :return:
        '''
        x_num,x_dim=x.shape #观测变量维度
        self.p=np.ones(shape=(K,1)) #隐变量对应概率，初始化全为1
        self.μ=np.random.random(size=(K,x_dim,1)) #均值，随机初始化
        self.Σ=np.ones(shape=(K,x_dim,x_dim)) #协方差，初始化全为1
        self.ε=1e-6 #精度，0.001
        while True: #循环迭代
            #multivariate_normal.pdf是从多维正态分布中计算均值为mean，协方差为cov的关于x的值。
            #由于本文设置了两个变量，所以要循环所有的均值和标准差，每次都计算出对应的pdf值。
            #最终计算出来的，实际上就是P(z|x)的分母。
            pdf_arr=np.array([stats.multivariate_normal.pdf(x,mean=i.reshape(-1),cov=j,allow_singular=True) for i,j in zip(self.μ,self.Σ)])
            p_k_all_x =(self.p.T@pdf_arr).T

            # 用于计算是否参数有更新
            loss=0
            #迭代所有的簇的参数，更新对应簇的参数
            for k in np.arange(K):
                #计算P(z|x)的分子
                p_kx= self.p[k]*stats.multivariate_normal.pdf(x,mean=self.μ[k].reshape(-1),cov=self.Σ[k],allow_singular=True)
                p_kx=p_kx.reshape(-1,1)

                #计算出P(z|x)，但不连加，仍保留初始状态
                result=(p_kx/p_k_all_x)

                #计算P(z|x)
                result_sum=result.reshape(-1).sum()
                #依据公式计算处新的p_k
                new_p=result_sum/x_num
                #依据公式计算出新的μ
                new_μ=(x.T@result)/result_sum

                #计算出新的Σ
                mole=0
                mean_x=x-self.μ[k].T #都减去均值
                for i in np.arange(x.shape[0]):
                    x_μ=(mean_x[i].reshape(-1,1))@(mean_x[i].reshape(-1,1).T) #计算(x-μ)(x-μ)^T
                    mole=mole+result[i]*x_μ #乘以对应的P(z_i=Ck|x_i)
                #计算新的Σ
                new_Σ=mole/result_sum
                loss+=np.abs((self.Σ[k]-new_Σ).sum())
                loss+=np.abs((self.μ[k]-new_μ).sum())
                loss+=np.abs((self.p[k]-new_p).sum())
                #更新参数
                self.p[k]=new_p
                self.μ[k]=new_μ
                self.Σ[k]=new_Σ
            #计算聚类结果
            pdf=np.array([stats.multivariate_normal.pdf(x,mean=i.reshape(-1),cov=j,allow_singular=True) for i,j in zip(self.μ,self.Σ)])
            y=np.argmax(pdf,axis=0)
            #绘图
            plot_figure(x,y)

            if loss<self.ε: #检查是否更新的值符合要求
                break
        return self.p,self.μ,self.Σ
#绘图函数
def plot_figure(x,y):
    color_map={0:"r",1:"b"}
    color=[ color_map[i] for i in y]
    plt.cla()
    plt.scatter(x[:,0],x[:,1],c=color)
    plt.pause(0.1)
if __name__ == '__main__':
    x1=stats.norm.rvs(-3,1,size=(200,2)) #从均值为-3，标准差为1的正态分布中抽取200个维度为2的数据
    x2=stats.norm.rvs(3,1,size=(150,2)) #从均值为3，标准差为1的正态分布中抽取200个维度为2的数据
    x=np.concatenate((x1,x2),axis=0) #合并数据
    gmm=GMM() #初始化模型
    p,mean,cov=gmm.train(x,2) #训练模型
    print(p)
    print(mean)
    print(cov)