高斯混合模型及python代码

单高斯模型

高斯模型是一种常用的变量分布模型，一维高斯分布的概率密度函数如下：

$f(x)=\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(x-\mu {)}^2}{2{\sigma }^2})$

$\mu$和${\sigma }^2$分别是高斯分布的均值和方差。

似然函数可写为：
$P(x|\theta )={\prod }_{i=1}^N[\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(x_i-\mu {)}^2}{2{\sigma }^2})]$

$\theta =\left\{\begin{array}{c}\mu ,{\sigma }^2\end{array}\right\}$

对数似然函数：
$\log P(x|\theta )=\log {\prod }_{i=1}^N[\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(x_i-\mu {)}^2}{2{\sigma }^2})]$

$={\sum }_{i=1}^N[-\log \sqrt{2\pi }-{\sum }_{i=1}^N-\log \sigma -\frac{1}{2{\sigma }^2}(x_i-\mu {)}^2]$

$\mu$的mle[最大似然估计]为：

$$\frac{\partial {\sum }_{i=1}^N[-\log \sqrt{2\pi }-{\sum }_{i=1}^N-\log \sigma -\frac{1}{2{\sigma }^2}(x_i-\mu {)}^2]}{\partial \mu }=\frac{\partial {\sum }_{i=1}^N-\frac{1}{2{\sigma }^2}(x_i-\mu {)}^2}{\partial \mu }=0$$

$⟹-{\sum }_{i=1}^N\frac{(x_i-u)}{{\sigma }^2}=0$

$⟹{\mu }_{MLE}=\frac{1}{N}{\sum }_{i=1}^N{x}_i$

同理可求：
${\sigma }_{MLE}^2=\frac{1}{N}{\sum }_{i=1}^N(x-\mu {)}^2$

多维引入

一组数据

如下图，我们有一组数据：

X1 = np.linspace(-5, 5, num=20)

# loc: mean 均值， scale: standard deviation 标准差
gauss1=norm.pdf(X1, loc=0, scale=2)

for i in range(len(X1)):
    plt.scatter(X1[i], 0,marker ='x',color='blue')
plt.plot(X1,gauss1)
plt.show()

很显然，我们可以看出这是一个高斯分布，而且我们很容易就可以求出其分布的参数

两组数据

但是如果我们有下面的数据的话：

X2=np.linspace(8, 16, num=20)
#数据合并
X = np.hstack((X1, X2))
gauss=norm.pdf(X, loc=np.mean(X), scale=np.std(X))

for i in range(len(X1)):
    plt.scatter(X1[i], 0,marker ='x',color='blue')
for i in range(len(X2)):
    plt.scatter(X2[i], 0,marker ='x',color='red')
plt.plot(X,gauss)
plt.show()

虽然我们绘出了其高斯图，也可以求出其相关的参数，比如这里的np.mean(X)=6
但是在图上我们可以发现，在x=6的位置，数据很稀疏，这与高斯分布的原理相悖，即在均值附近的点是最密集的。

组合

所以我们如果使用两个高斯分布来拟合的数据的话就很合理了
如下：

X2=np.linspace(8, 16, num=20)
gauss2=norm.pdf(X2, loc=12, scale=2)
plt.figure(figsize=(12, 6))
for i in range(len(X1)):
    plt.scatter(X1[i], 0,marker ='x',color='blue')
plt.plot(X1,gauss1)

for i in range(len(X2)):
    plt.scatter(X2[i], 0,marker ='x',color='red')
plt.plot(X2,gauss2)

plt.show()

直观图

我们接下来用二维数据来更直观的观测

true_Mu = [[0.5, 0.5], [5.5, 2.5], [1, 7]]
true_Var = [[1, 3], [2, 2], [6, 2]]
num1, mu1, var1 = 400, true_Mu[0], true_Var[0]
X1 = np.random.multivariate_normal(mu1, np.diag(var1), num1)
# 第二簇的数据
num2, mu2, var2 = 600, true_Mu[1], true_Var[1]
X2 = np.random.multivariate_normal(mu2, np.diag(var2), num2)
# 第三簇的数据
num3, mu3, var3 = 1000, true_Mu[2], true_Var[2]
X3 = np.random.multivariate_normal(mu3, np.diag(var3), num3)
# 合并在一起
X = np.vstack((X1, X2, X3))
# 显示数据
plt.figure(figsize=(10, 8))
plt.axis([-10, 15, -5, 15])
plt.scatter(X1[:, 0], X1[:, 1], s=5)
plt.scatter(X2[:, 0], X2[:, 1], s=5)
plt.scatter(X3[:, 0], X3[:, 1], s=5)
plt.show()

我们构建了三簇数据，这很明显是三个高斯分布

# 显示数据

fig, ax1 = plt.subplots(figsize=(10, 8))

plt.scatter(X1[:, 0], X1[:, 1], s=5, color="blue")
plt.scatter(X2[:, 0], X2[:, 1], s=5, color="green")
plt.scatter(X3[:, 0], X3[:, 1], s=5, color="red")

#副坐标轴
ax2 = ax1.twinx()
sns.distplot(X1[:, 0],hist=False, color="blue")
sns.distplot(X2[:, 0],hist=False, color="green")
sns.distplot(X3[:, 0],hist=False, color="red")

plt.show()

现在我们想把这多个高斯模型放在一个模型里面该如何解决呢？

高斯混合模型

高斯混合模型可以看作是由 K 个单高斯模型组合而成的模型，这 K 个子模型是混合模型的隐变量（Hidden variable）。一般来说，一个混合模型可以使用任何概率分布，这里使用高斯混合模型是因为高斯分布具备很好的数学性质以及良好的计算性能。
设有随机变量X，则混合高斯模型可以用下式表示：
$P(x)={\sum }_{k=1}^K{\alpha }_{k}N({\mu }_k,{\Sigma }_k)$

高斯混合模型（Gaussian Mixture Model，GMM）是一种软聚类模型。 GMM也可以看作是K-means的推广，因为GMM不仅是考虑到了数据分布的均值，也考虑到了协方差。和K-means一样，我们需要提前确定簇的个数。

1. 参数设定

$X=(x_1,x_2,\dots \dots ,x_N)$
$Z=(z_1,z_2,\dots \dots ,z_N)$
$\theta =(p,\mu ,\sum )$

其中参数$\theta$包含隐变量$Z$的概率分布，各个高斯的均值与协方差矩阵：
$p=(p_1,p_2,\dots \dots ,p_K)$
$\mu =({\mu }_1,{\mu }_2,\dots \dots ,{\mu }_K)$
$\sum =({\sum }_1,{\sum }_2,\dots \dots ,{\sum }_K)$
$K$是值高斯模型的个数

2. EM算法-E步

$Q$函数：
$Q(\theta ,{\theta }^{(t)})={\mathbb{E}}_{Z|X,{\theta }^{(t)}}[logP(X,Z|\theta )]$

$={\sum }_{Z}logP(X,Z|\theta )P(Z|X,{\theta }^{(t)})$

对上式进行分别简化：
先看右半部分

$P(Z|X,{\theta }^{(t)})={\prod }_{i=1}^{N}P(z_i|x_i,{\theta }^{(t)})$

$={\prod }_{i=1}^N[\frac{P(x_i|z_i)P(z_i)}{{\sum }_{z_i=1}^{K}P(x_i|z_i)P(z_i)}]$

几何意义

对于任意一个点$x_i$会受到上图中两个高斯分布的影响，每个高斯分布对于点的权重不一样，权重的大小即为高斯概率密度函数的值，为了方便需要对其进行归一化，即$P(z_i=left|x_i,\theta )=\frac{a}{a+b}$，$P(z_i=right|x_i,\theta )=\frac{b}{a+b}$，为了计算整体的权重需要对所有点的权重进行加总并平均。

左边[为了方便计算右边简写为$P(z_1,z_2,\dots \dots ,z_N)$]

$={\sum }_{Z}logP(X,Z|\theta )P(Z|X,{\theta }^{(t)})$

$={\sum }_{z_1=1}^K{\sum }_{z_2=1}^K...{\sum }_{z_N=1}^K({\sum }_{i=1}^N\underset{f_i(z_i)}{\underbrace{[log{\alpha }_{z_i}+logN(x_i|{\mu }_{z_i},{\sigma }_{z_i}^2)]}})P(z_1,z_2,\dots \dots ,z_N)$

$={\sum }_{z_1=1}^K{\sum }_{z_2=1}^K...{\sum }_{z_N=1}^K(f_1(z_1)+f_2(z_2)+\dots \dots +f_N(z_N))P(z_1,z_2,\dots \dots ,z_N)$

$={\sum }_{z_1=1}^K{f}_1(z_1){\sum }_{z_2=1}^K...{\sum }_{z_N=1}^{K}P(z_1,z_2,\dots \dots ,z_N)+\dots \dots$

对于红色部分我们使用边缘概率分布特性可以得到

$={\sum }_{z_1=1}^K{f}_1(z_1)P(z_1)+......$

$={\sum }_{i=1}^N{\sum }_{z_i=1}^K{f}_i(z_i)P(z_i)$

$={\sum }_{i=1}^N{\sum }_{z_i=1}^K[log{\alpha }_{z_i}+logN(x_i|{\mu }_{z_i},{\sigma }_{z_i}^2)]P(z_i|x_i,{\theta }^{(t)})$

2. EM算法-M步

根据确定的Q函数求下一时刻的参数，首先是Z的概率分布：
${\alpha }^{(t+1)}=arg\underset{\alpha }{\underbrace{max}}{\sum }_{i=1}^N{\sum }_{z_i=1}^K[log{\alpha }_{z_i}+logN(x_i|{\mu }_{z_i},{\sigma }_{z_i}^2)]P(z_i|x_i,{\theta }^{(t)})$

这里不能直接求导，因为概率分布有一个约束：${\sum }_{i=1}^K{\alpha }_i=1$
所以需要运用拉格朗日乘子法进行求解

$L(\alpha ,\lambda )={\sum }_{i=1}^N{\sum }_{z_i=1}^{K}log{\alpha }_{z_i}P(z_i|x_i,{\theta }^{(t)})+\lambda ({\sum }_{i=1}^K{\alpha }_i-1)$

对$\alpha$求导并令导数等于0

$\frac{\partial L(\alpha ,\lambda )}{\partial \alpha }={\sum }_{i=1}^N{\sum }_{z_i=1}^K\frac{P(z_i|x_i,{\theta }^{(t)})}{{\alpha }_i}+\lambda =0$

$⟹{\sum }_{i=1}^N{\sum }_{z_i=1}^{K}P(z_i|x_i,{\theta }^{(t)})=-\lambda {\sum }_{i=1}^K{\alpha }_i$

带入条件：${\sum }_{i=1}^K{\alpha }_i=1$，$P(z_i|x_i,{\theta }^{(t)})$也是一个概率分布，所以有${\sum }_{z_i=1}^{K}P(z_i|x_i,{\theta }^{(t)})=1$

所以可得：$\lambda =-N$

带入可得：
${\alpha }^{(t+1)}=\frac{1}{N}{\sum }_{i=1}^{N}P(z_i|x_i,{\theta }^{(t)})$

$P(z_i|x_i,{\theta }^{(t)})$是T时刻的值，所以是已知的。

同理可以求其他参数：

${\mu }^{(t+1)}=\frac{{\sum }_{i=1}^{N}P(z_i|x_i,{\theta }^{(t)})x_i}{{\sum }_{i=1}^{N}P(z_i|x_i,{\theta }^{(t)})}$

${\sum }^{(t+1)}=\frac{{\sum }_{i=1}^{N}P(z_i|x_i,{\theta }^{(t)})(x_i-\mu {)}^2}{{\sum }_{i=1}^{N}P(z_i|x_i,{\theta }^{(t)})}$

\color{颜色}{文字}
\underbrace{a+b+c+d}_{Sample}

python实现

import numpy as np
from scipy import stats
from matplotlib.patches import Ellipse

class GMM(object):
    def __init__(self, k: int, d: int):
        '''
        k: K值即类别数
        d: 样本属性的数量即特征数量，为了可视化本例取2
        '''
        self.K = k
        # 初始化参数
        self.p = np.random.rand(k)
        self.p = self.p / self.p.sum()      # 保证所有p_k的和为1
        self.means = np.random.rand(k, d)#初始化均值
        self.covs = np.empty((k, d, d))#初始化方差
        for i in range(k):                  # 随机生成协方差矩阵，必须是半正定矩阵
            self.covs[i] = np.eye(d) * np.random.rand(1) * k
            
    # 画出聚类图像
    def plot_clusters(self,data,  Mu=None, Var=None):
        colors = ['b', 'g', 'r']
        n_clusters = len(self.means)
        plt.figure(figsize=(10, 8))
        plt.axis([-10, 15, -5, 15])
        plt.scatter(data[:, 0], data[:, 1], s=5)
        ax = plt.gca()
        for i in range(n_clusters):
            #迭代模型区域
            #Ellipse，参数：坐标，椭圆胖瘦
            ellipse = Ellipse(list(self.means[i]), 3 * self.covs[i][0][0], 3 * self.covs[i][1][1],alpha= 0.5 ,lw= 2, edgecolor=colors[i],ls='--')
            
            ax.add_patch(ellipse)
        if (Mu is not None) & (Var is not None):
            #原始数据
            for i in range(n_clusters):
                #plot_args = {'fc': 'None', 'lw': 2, 'edgecolor': colors[i], 'alpha': 0.5}
                ellipse = Ellipse(Mu[i], 3 * Var[i][0], 3 * Var[i][1],  fc='None',lw= 2, edgecolor=colors[i])
                ax.add_patch(ellipse)         
        plt.show()


    def fit(self, data: np.ndarray):
        '''
        data: 数据矩阵，每一行是一个样本，shape = (N, d)
        '''
        for _ in range(100):
            density = np.empty((len(data), self.K))
            for i in range(self.K):
                # 生成K个概率密度函数并计算对于所有样本的概率密度
                norm = stats.multivariate_normal(self.means[i], self.covs[i])
                density[:,i] = norm.pdf(data)
            # 计算所有样本属于每一类别的后验
            posterior = density * self.p#概率*权重
            posterior = posterior / posterior.sum(axis=1, keepdims=True)#归一化:N*3
            #print(posterior)
            # 计算下一时刻的参数值
            p_hat = posterior.sum(axis=0)
            #print(p_hat)
            mean_hat = np.tensordot(posterior, data, axes=[0, 0])#3*N*N*2=3*2
            #print(mean_hat)
            # 计算协方差
            cov_hat = np.empty(self.covs.shape)
            for i in range(self.K):
                tmp = data - self.means[i]
                cov_hat[i] = np.dot(tmp.T*posterior[:,i], tmp) / p_hat[i]
            # 更新参数
            self.covs = cov_hat
            self.means = mean_hat / p_hat.reshape(-1,1)
            self.p = p_hat / len(data)
            if _ % 10 ==0:
                GMM.plot_clusters(X,Mu,Var)
        print(self.p)
        print(self.means)
        print(self.covs)

阴影部分是迭代结果，初始化的数据跟原始数据【实线数据】有较大差距

随着迭代：

可以看到迭代结果跟原始数据渐渐重合。

预测结果
权重：
[0.33384088 0.44272027 0.22343885]

均值：
[[5.53894562 2.59218286]
[1.14375546 7.02480626]
[0.41072795 0.54925972]]

方差：
[[[ 2.01873055 0.02578939]
[ 0.02578939 2.01838067]]

[[ 6.10498076 0.20116635]
[ 0.20116635 2.00151559]]

[[ 0.93262483 -0.09693276]
[-0.09693276 3.0662879 ]]]

输入真实值：
Mu = [[0.5, 0.5], [5.5, 2.5], [1, 7]]
Var = [[1, 3], [2, 2], [6, 2]]

完整代码
参考