GMM混合高斯模型

高斯混合模型（Gaussian Mixture Model）通常简称GMM，是一种业界广泛使用的聚类算法，该方法使用了高斯分布作为参数模型，并使用了期望最大（Expectation Maximization，简称EM）算法进行训练。

参考资料：
知乎：如何通俗的理解概率密度函数
知乎：如何理解条件概率
知乎：如何通俗的理解高斯混合模型（Gaussian Mixture Models）
知乎：高斯混合模型（GMM）推导及实现

最近研究经典论文：GrabCut, 里面讲到混合高斯模型GMM，特地学习一下。以下一维高斯混合模型的代码段来自：https://blog.csdn.net/u012436149/article/details/53557008

import numpy as np
def gaussian(x,mu,sigma):
    temp = -np.square(x-mu)/(2*sigma)
    return np.exp(temp)/(np.sqrt(2.0*np.pi*sigma)) # sigma = sigma^2
    
def e_step(data, pis, mus, sigmas):
    Qs = []
    # Qs的每一行表示每个data属于某一个分布的概率，概率之和不为1
    for i in range(len(data)):
        q = [pi*gaussian(data[i],mu,sigma) for pi,mu,sigma in zip(pis,mus,sigmas)]
        Qs.append(q)

    # 将Qs归一化：每个data属于每个分布的概率，概率之和为1
    Qs = np.array(Qs)
    Qs = Qs / np.sum(Qs,axis=1).reshape(-1,1) #  axis =1, 在每一行求和
    return Qs
    
def m_step(data, pis, mus, sigmas, Qs):
    data = np.array(data)
    gama_j = np.sum(Qs,axis=0)          # 每一类的概率之和
    new_pis = gama_j/len(data)        # 所有data属于每一个分布的平均概率

    mu_temp = np.sum(Qs*(data.reshape(-1,1)),axis=0) # 把每个data按概率切分到每一类，再计算每一类的总和
    new_mus = mu_temp/gama_j            # 计算均值

    X_i_mu_j = np.square(np.array([data]).reshape(-1,1)-np.array([mus]))
    new_sigmas = np.sum(Qs*X_i_mu_j,axis=0)/gama_j

    return new_pis, new_mus, new_sigmas
    
def EM(data,k=1):
    # 设置均值
    pis = [1.0/k for i in range(k)] 	# 各高斯模型的初始系数
    mus = [i for i in range(k)]       	# 混合高斯的初始均值
    sigmas = [1 for i in range(k)]    	# 混合高斯的初始标准差
    # 开始学习
    for i in range(30):
        Qs = e_step(data,pis,mus,sigmas)
        pis, mus, sigmas= m_step(data,pis,mus,sigmas,Qs)
        pis = np.around(pis,5)
        mus = np.around(mus,5)
        sigmas = np.around(sigmas,5)
        print(i+1,"---- N1:[",pis[0],mus[0],sigmas[0],"], N2:[",pis[1],mus[1],sigmas[1],"]")

if __name__ == "__main__":
    s1 = np.random.normal(10,1,10000) 	# 均值，标准差，数目
    s2 = np.random.normal(20,1,30000)	# 均值，标准差，数目
    s3 = np.concatenate((s1,s2),axis=0)	# concat
    np.random.shuffle(s3)               # 数组置乱

    EM(s3,2)

输出：

1 ---- N1:[ 3e-05 8.9876 81.83219 ], N2:[ 0.99997 17.50295 292.15809 ]
2 ---- N1:[ 4e-05 15.61764 69.1004 ], N2:[ 0.99996 17.50276 19.81197 ]
3 ---- N1:[ 4e-05 14.56491 29.01452 ], N2:[ 0.99996 17.50279 19.8117 ]
4 ---- N1:[ 4e-05 14.00006 26.48641 ], N2:[ 0.99996 17.50284 19.8115 ]
5 ---- N1:[ 4e-05 13.43604 25.08787 ], N2:[ 0.99996 17.50286 19.81137 ]
6 ---- N1:[ 4e-05 12.86251 23.14873 ], N2:[ 0.99996 17.50289 19.81122 ]
7 ---- N1:[ 4e-05 12.24169 20.37554 ], N2:[ 0.99996 17.50292 19.81105 ]
8 ---- N1:[ 5e-05 11.51537 16.26617 ], N2:[ 0.99995 17.50296 19.81084 ]
9 ---- N1:[ 6e-05 10.63522 10.04298 ], N2:[ 0.99994 17.50311 19.81006 ]
10 ---- N1:[ 9e-05 9.82837 2.84751 ], N2:[ 0.99991 17.50337 19.80862 ]
11 ---- N1:[ 0.00023 9.66224 0.78311 ], N2:[ 0.99977 17.50451 19.80246 ]
12 ---- N1:[ 0.00087 9.63772 0.4433 ], N2:[ 0.99913 17.50954 19.77528 ]
13 ---- N1:[ 0.00361 9.62874 0.3113 ], N2:[ 0.99639 17.5312 19.65804 ]
14 ---- N1:[ 0.01476 9.62727 0.24855 ], N2:[ 0.98524 17.62069 19.17018 ]
15 ---- N1:[ 0.04924 9.63664 0.23563 ], N2:[ 0.95076 17.91003 17.53966 ]
16 ---- N1:[ 0.1066 9.66379 0.28948 ], N2:[ 0.8934 18.43801 14.21387 ]
17 ---- N1:[ 0.16483 9.71534 0.43784 ], N2:[ 0.83517 19.03959 9.66799 ]
18 ---- N1:[ 0.21921 9.81934 0.70377 ], N2:[ 0.78079 19.65987 4.33717 ]
19 ---- N1:[ 0.2482 9.97096 0.96095 ], N2:[ 0.7518 19.98927 1.2408 ]
20 ---- N1:[ 0.25 9.98935 0.98003 ], N2:[ 0.75 20.00714 1.00113 ]
21 ---- N1:[ 0.25 9.98935 0.97969 ], N2:[ 0.75 20.00714 1.00081 ]
22 ---- N1:[ 0.25 9.98935 0.97969 ], N2:[ 0.75 20.00714 1.00081 ]
23 ---- N1:[ 0.25 9.98935 0.97969 ], N2:[ 0.75 20.00714 1.00081 ]
24 ---- N1:[ 0.25 9.98935 0.97969 ], N2:[ 0.75 20.00714 1.00081 ]
25 ---- N1:[ 0.25 9.98935 0.97969 ], N2:[ 0.75 20.00714 1.00081 ]
26 ---- N1:[ 0.25 9.98935 0.97969 ], N2:[ 0.75 20.00714 1.00081 ]
27 ---- N1:[ 0.25 9.98935 0.97969 ], N2:[ 0.75 20.00714 1.00081 ]
28 ---- N1:[ 0.25 9.98935 0.97969 ], N2:[ 0.75 20.00714 1.00081 ]
29 ---- N1:[ 0.25 9.98935 0.97969 ], N2:[ 0.75 20.00714 1.00081 ]
30 ---- N1:[ 0.25 9.98935 0.97969 ], N2:[ 0.75 20.00714 1.00081 ]

我发现很多博客中绘制的高斯混合模型的示意图不对，因为混合高斯模型要求各单个高斯模型的权重之和为1，但很多图看上去单个高斯分布的权重都为1，如下面的这张图就是错误的，非常误导人。

于是这里我重新绘制了一张正确的高斯混合模型的示意图（各个高斯分布的权重和分别为：1/2,1/3,1/6,为了看的更清楚，这里将函数的值扩大了60倍），需要注意的是下面的图跟上面的代码没有关系。

使用 https://www.desmos.com/calculator?lang=zh-CN 进行函数绘制。