05 EM算法 - 高斯混合模型 - GMM

Posted 2023-05-06

参考技术A

04 EM算法 - EM算法收敛证明

GMM (Gaussian Mixture Model, 高斯混合模型)是指该算法由多个高斯模型线性叠加混合而成。每个高斯模型称之为component。

GMM算法 描述的是数据的本身存在的一种分布，即样本特征属性的分布，和预测值Y无关。显然GMM算法是无监督的算法，常用于聚类应用中，component的个数就可以认为是类别的数量。

回到昨天说的例子：随机选择1000名用户，测量用户的身高；若样本中存在男性和女性，身高分别服从高斯分布N(μ1,σ1)和N(μ2,σ2)的分布，试估计参数:μ1,σ1,μ2,σ2；

1、如果明确的知道样本的情况(即男性和女性数据是分开的)，那么我们使用极大似然估计来估计这个参数值。

2、如果样本是混合而成的，不能明确的区分开，那么就没法直接使用极大似然估计来进行参数的估计。

我们可以认为当前的1000条数据组成的集X，是由两个高斯分布叠加而成的(男性的分布和女性的分布)。

如果能找到一种办法把每一个高斯分布对应的参数π、 μ、σ求出来，那么对应的模型就求解出来了。

如果模型求解出来后，如何对数据进行聚类？

这个公式求出来的分别是男性和女性身高分布的概率密度，如果把π、 μ、σ都求出来，以后我们可以构建出一个 能够根据样本特征 计算出样本属于男性或女性的可能性。

实际做样本分类的时候，我们把样本X的特征x1~xn分别代入两个公式中，求出来的两个结果分别是：样本X的性别是男、是女的可能性。如果是男的可能性大于是女的可能性，我们就把样本X归入男性的分类。

假定 GMM 由k个Gaussian分布线性叠加而成，那么概率密度函数如下：

分析第1个等式：
p(x)： 概率密度函数，k个Gaussian分布线性叠加而成的概率密度函数。
∑p(k)p(x|k)： k个某种模型叠加的概率密度函数。
p(k)： 每个模型占的权重，即上面提到的π。
p(x|k)： 给定类别k后，对应的x的概率密度函数。

分析第2个等式： 目标 - 将公式写成高斯分布的样子。
π k ： 即p(k)
p(x;μ k ,∑ k )： 多元高斯(正态)分布。有了观测数据x后，在 给定了条件 下的高斯分布。这个 条件 是 1、第k个分类的均值μ k ; 2、第k个分类的方差∑ k ;

深入分析p(x;μ k ,∑ k )的参数：
如果样本有n个特征，所有的特征x1~xn一起服从一个多元的高斯分布(正态分布)，所有特征的均值应该是一个向量 (μ 1 ~μ n )；
μ k ： 第k个分类的情况下(第k个高斯分布的情况下对应的每一列的均值)；μ k = (μ k1 ~μ kn )

∑ k ： 协方差矩阵(对称阵)。现在有n个特征，协方差矩阵是一个n×n的矩阵。现在我们要算的是：

cov(x1,x1)，cov(x1,x2)，...，cov(x1,xn)

cov(x2,x1)，cov(x2,x2)，...，cov(x2,xn)
....
cov(xn,x1)，cov(x1,x2)，...，cov(xn,xn)

其中， 对角线 cov(x1,x1)、cov(x2,x2)， ... ，cov(xn,xn)中，x1和x1的协方差 = x1的方差；即cov(x1,x1) = var(x1)；所以 对角线上两个特征的协方差 = 对应的特征的方差。

协方差 （Covariance）在 概率论 和 统计学 中用于衡量两个变量的总体 误差 。而 方差 是协方差的一种特殊情况，即当两个变量是相同的情况。

协方差表示的是两个变量的总体的 误差 ，这与只表示一个变量误差的 方差 不同。 如果两个 变量 的变化趋势一致，也就是说如果其中一个大于自身的期望值，另外一个也大于自身的期望值，那么两个变量之间的协方差就是正值。 如果两个变量的变化趋势相反，即其中一个大于自身的期望值，另外一个却小于自身的期望值，那么两个变量之间的协方差就是负值。

理解了公式后，再来看看公式在图像上是如何体现的：

如果样本X只有一个特征x1，在二维的坐标系上的表示出来。特征x1是由n个单变量样本的高斯分布叠加而成的。向量x1 k = ∑ k (x1 (1) ,x1 (2) ,~,x1 (n) )，如k=(男、女)，累加男性分类下的特征高斯分布和女性分类下的高斯分布；

图中 红色曲线 表示原有数据的分布情况，我认为这个原有数据是由多个比较的高斯分布叠加而成的， 蓝色曲线 表示单个单个高斯分布的分布情况。向量x1 = (x1 (1) ,x1 (2) ,~,x1 (n) )；

PS： 蓝1+蓝2=红 体现的就是公式 p(x) = ∑πp(x;μ,∑k)；

在得知数据的特征 x=(x1~xn) 后，如果我们想把数据合理得聚类到一个分类中，我们该如何去计算呢？

既然我已经得到了k个高斯分布对应的概率密度函数(现在设k=3，共3个分类)，将当前特征的x=(x1~xn)代入我们的概率密度函数： p(x) = ∑πp(x;μ,∑k)；

我们分别计算p(蓝1)、p(蓝2)、p(蓝3)，蓝色三条线各对应k分类中的一个，哪个数大，我认为当前的样本该分到哪一类。

GMM算法的两个前提：
1、数据服从高斯分布；
2、我们人为定义了分类个数k。

问：我们人为假定了高斯分布的分类个数k，就类似于我们聚簇时分的聚簇中心个数一样。参数π、μ、σ该如何求出来?

答：和K-Means算法一样，我们可以用 EM算法 来求解这个问题。 GMM也满足EM算法的聚类思想，首先人为得定义了聚类的个数k，从数据特征X中发掘潜在关系的一种模型。而且我还默认数据是服从多个高斯分布的。

GMM算法中的隐含条件是：第k个模型占的权重 - 、 第k个高斯分布的情况下对应的每一列的均值 - 、协方差矩阵 cov(xi,xj) - ；因为本质上我们是知道数据原有的分类状况的，只是无法观测到隐含在数据中的这些特性，使用EM的思想可以迭代得求解出这些隐含变量。

对联合概率密度函数求对数似然函数：

对联合概率密度函数求对数后，原本 连乘 的最大似然估计变成了 连加 的函数状态。

EM算法求解 - E步：

套用公式后，我们可以假定隐含变量z的分布：Q(z (i) = j)；
我们认为分布wj (i) = 第i个观测值对应的隐含分类第z (i) 类； = 以(看不见的参数π、μ、∑)为参数的情况下，输入第i观测值的特征x后得到的分类z (i) 类；

EM算法求解 - M步：
M步第1行就是上一章通过化简找到 下界 的那个函数：

1、对均值求偏导:

2、对方差求偏导:

3、对概率使用拉格朗日乘子法求解:

06 EM算法 - 案例一 - EM分类初识及GMM算法实现

6. EM算法-高斯混合模型GMM+Lasso详细代码实现

1. 前言

我们之前有介绍过4. EM算法-高斯混合模型GMM详细代码实现，在那片博文里面把GMM说涉及到的过程，可能会遇到的问题，基本讲了。今天我们升级下，主要一起解析下EM算法中GMM（搞事混合模型）带惩罚项的详细代码实现。

2. 原理

由于我们的极大似然公式加上了惩罚项，所以整个推算的过程在几个地方需要修改下。

在带penality的GMM中，我们假设协方差是一个对角矩阵，这样的话，我们计算高斯密度函数的时候，只需要把样本各个维度与对应的(mu_k)和(sigma_k)计算一维高斯分布，再相加即可。不需要通过多维高斯进行计算，也不需要协方差矩阵是半正定的要求。

我们给上面的(1)式加入一个惩罚项，
[ lambdasum_{k=1}^Ksum_{j=1}^Pfrac{|mu_k-ar{x}_j|}{s_j} ]
其中的(P)是样本的维度。(ar{x}_j)表示每个维度的平均值，(s_j)表示每个维度的标准差。这个penality是一个L1范式，对(mu_k)进行约束。

加入penality后(1)变为
[ L( heta, heta^{(j)})=sum_{k=1}^Kn_k[logpi_k-frac{1}{2}(log(oldsymbol{Sigma_k})+frac{{(x_i-oldsymbol{mu}_k})^2}{oldsymbol{Sigma}_k})] - lambdasum_{k=1}^Ksum_{j=1}^Pfrac{|mu_k-ar{x}_j|}{s_j} ]

这里需要注意的一点是，因为penality有一个绝对值，所以在对(mu_k)求导的时候，需要分情况。于是(2)变成了
[ mu_k=frac{1}{n_k}sum_{i=1}^Ngamma_{ik}x_i ]
[ mu_k= left {egin{array}{cc} frac{1}{n_k}(sum_{i=1}^Ngamma_{ik}x_i - frac{lambdasigma^2}{s_j}), & mu_k >= ar{x}_j\\frac{1}{n_k}(sum_{i=1}^Ngamma_{ik}x_i + frac{lambdasigma^2}{s_j}), & mu_k < ar{x}_j end{array} ight. ]

3. 算法实现

• 和不带惩罚项的GMM不同的是，我们GMM+LASSO的计算高斯密度函数有所变化。

#计算高斯密度概率函数，样本的高斯概率密度函数，其实就是每个一维mu,sigma的高斯的和
def log_prob(self, X, mu, sigma):
    N, D = X.shape
    logRes = np.zeros(N)
    for i in range(N):
        a = norm.logpdf(X[i,:], loc=mu, scale=sigma)
        logRes[i] = np.sum(a)
    return logRes

• 在m-step中计算(mu_{k+1})的公式需要变化，先通过比较(mu_{kj})和(means_{kj})的大小，来确定绝对值shift的符号。

def m_step(self, step):
    gammaNorm = np.array(np.sum(self.gamma, axis=0)).reshape(self.K, 1)
    self.alpha = gammaNorm / np.sum(gammaNorm)
    for k in range(self.K):
        Nk = gammaNorm[k]
        if Nk == 0:
            continue
        for j in range(self.D):
            if step >= self.beginPenaltyTime:
                # 算出penality的偏移量shift，通过当前维度的mu和样本均值比较，确定shift的符号，相当于把lasso的绝对值拆开了
                shift = np.square(self.sigma[k, j]) * self.penalty / (self.std[j] * Nk)
                if self.mu[k, j] >= self.means[j]:
                    shift = shift
                else:
                    shift = -shift
            else:
                shift = 0
            self.mu[k, j] = np.dot(self.gamma[:, k].T, self.X[:, j]) / Nk - shift
            self.sigma[k, j] = np.sqrt(np.sum(np.multiply(self.gamma[:, k], np.square(self.X[:, j] - self.mu[k, j]))) / Nk)

• 最后需要修改loglikelihood的计算公式

def GMM_EM(self):
    self.init_paras()
    for i in range(self.times):
        #m step
        self.m_step(i)
        # e step
        logGammaNorm, self.gamma= self.e_step(self.X)
        #loglikelihood
        loglike = self.logLikelihood(logGammaNorm)
        #penalty
        pen = 0
        if i >= self.beginPenaltyTime:
            for j in range(self.D):
                pen += self.penalty * np.sum(abs(self.mu[:,j] - self.means[j])) / self.std[j]

        # print("step = %s, alpha = %s, loglike = %s"%(i, [round(p[0], 5) for p in self.alpha.tolist()], round(loglike - pen, 5)))
        # if abs(self.loglike - loglike) < self.tol:
        #     break
        # else:

        self.loglike = loglike - pen

4. GMM算法实现结果

用我实现的GMM+LASSO算法，对多个penality进行计算，选出loglikelihood最大的k和penality，与sklearn的结果比较。

fileName = amix1-est.dat, k = 2, penalty = 0 alpha = [0.52838, 0.47162], loglike = -693.34677
fileName = amix1-est.dat, k = 2, penalty = 0 alpha = [0.52838, 0.47162], loglike = -693.34677
fileName = amix1-est.dat, k = 2, penalty = 1 alpha = [0.52789, 0.47211], loglike = -695.26835
fileName = amix1-est.dat, k = 2, penalty = 1 alpha = [0.52789, 0.47211], loglike = -695.26835
fileName = amix1-est.dat, k = 2, penalty = 2 alpha = [0.52736, 0.47264], loglike = -697.17009
fileName = amix1-est.dat, k = 2, penalty = 2 alpha = [0.52736, 0.47264], loglike = -697.17009
myself GMM alpha = [0.52838, 0.47162], loglikelihood = -693.34677, bestP = 0
sklearn GMM alpha = [0.53372, 0.46628], loglikelihood = -176.73112
succ = 299/300
succ = 0.9966666666666667
[0 1 0 0 1 1 0 1 1 1 0 0 1 0 0 1 0 0 0 1]
[0 1 0 0 1 0 0 1 1 1 0 0 1 0 0 1 0 0 0 1]
fileName = amix1-tst.dat, loglike = -2389.1852339407087
fileName = amix1-val.dat, loglike = -358.1157431278091
fileName = amix2-est.dat, k = 2, penalty = 0 alpha = [0.56, 0.44], loglike = 53804.54265
fileName = amix2-est.dat, k = 2, penalty = 0 alpha = [0.82, 0.18], loglike = 24902.5522
fileName = amix2-est.dat, k = 2, penalty = 1 alpha = [0.82, 0.18], loglike = 23902.65183
fileName = amix2-est.dat, k = 2, penalty = 1 alpha = [0.56, 0.44], loglike = 52929.96459
fileName = amix2-est.dat, k = 2, penalty = 2 alpha = [0.82, 0.18], loglike = 22907.40397
fileName = amix2-est.dat, k = 2, penalty = 2 alpha = [0.82, 0.18], loglike = 22907.40397
myself GMM alpha = [0.56, 0.44], loglikelihood = 53804.54265, bestP = 0
sklearn GMM alpha = [0.56217, 0.43783], loglikelihood = 11738677.90164
succ = 200/200
succ = 1.0
[0 1 0 0 1 0 1 1 0 0 0 1 1 0 1 0 1 1 0 1]
[0 1 0 0 1 0 1 1 0 0 0 1 1 0 1 0 1 1 0 1]
fileName = amix2-tst.dat, loglike = 51502.878096147084
fileName = amix2-val.dat, loglike = 6071.217012747491
fileName = golub-est.dat, k = 2, penalty = 0 alpha = [0.575, 0.425], loglike = -24790.19895
fileName = golub-est.dat, k = 2, penalty = 0 alpha = [0.525, 0.475], loglike = -24440.82743
fileName = golub-est.dat, k = 2, penalty = 1 alpha = [0.55, 0.45], loglike = -25582.27485
fileName = golub-est.dat, k = 2, penalty = 1 alpha = [0.6, 0.4], loglike = -26137.97508
fileName = golub-est.dat, k = 2, penalty = 2 alpha = [0.55, 0.45], loglike = -26686.02411
fileName = golub-est.dat, k = 2, penalty = 2 alpha = [0.55, 0.45], loglike = -26941.68964
myself GMM alpha = [0.525, 0.475], loglikelihood = -24440.82743, bestP = 0
sklearn GMM alpha = [0.5119, 0.4881], loglikelihood = 13627728.10766
succ = 29/40
succ = 0.725
[0 1 0 1 0 1 0 1 0 1 0 0 0 0 0 1 1 1 0 1]
[0 1 0 1 1 1 0 0 1 1 0 1 0 1 0 0 1 1 0 0]
fileName = golub-tst.dat, loglike = -12949.606698037718
fileName = golub-val.dat, loglike = -11131.35137056415

5. 总结

通过一番改造，实现了GMM+LASSO的代码，如果读者有什么好的改进方法，或者我有什么错误的地方，希望多多指教。