简介

本文章实现了Haojun Sun提出的一种计算高斯混合模型（GMM）重叠率的方法（论文：Measuring the component overlapping in the Gaussian mixture model）。这篇文论提出的方法可以计算任意两个混合高斯分布之间的重叠度。该方法可以用来评价GMM模型的好坏，我在我的论文中使用了这个算法，用来评价高斯混合模型聚类的可分性。

关于高斯混合模型（GMM）的相关概念可以参考另一篇博文：高斯混合模型及其EM算法的理解

使用GMM聚类或分析两个高斯混合分布的数据时，我们有时会希望两个高斯分布离得越远越好，这样表示数据才有可分性。但很多情况下两个高斯分布会有重叠。一维和二维的重叠情况如下所示（图片取自作者论文）。

我们可以计算一些指标来间接反映两个高斯分布的重叠情况。比如可以计算Mahalanobis距离，Bhattacharyya距离或Kullback-Leibler (KL)距离，可以衡量两个高斯分布的相似性。但是Mahalanobis距离预设两个分布具有相同的协方差，Bhattacharyya距离和KL距离都考虑了协方差，但却没有考虑高斯混合分布的混合系数（mixing coefficient）。而且KL距离对高维的正态分布没有解析解，计算复杂。

这篇论文提出的计算OLR的方法考虑了高斯混合分布中的所有参数，包括均值，协方差和混合系数。

OLR计算

假设有nnn个ddd维的样本X={X1,...,Xn}\boldsymbol{X} = \{X_1,..., X_n\}X={X1,...,Xn}. 其中XiX_iXi是一个ddd维向量。一个混合高斯模型的pdf可以表示为：
p(X)=∑i=1kαiGi(X,μi,Σi)(1)p(X) = \sum_{i=1}^k \alpha_iG_i(X, \mu_i, \Sigma_i) \tag{1}p(X)=i=1∑kαiGi(X,μi,Σi)(1)
其中αi\alpha_iαi是混合系数，满足αi>0\alpha_i > 0αi>0且∑i=1kαi=1\sum_{i=1}^k\alpha_i=1∑i=1kαi=1.

Gi(X)G_i(X)Gi(X)是一个ddd维高斯分布，可以表示为下面的形式：
Gi(X)=1(2π)d/2∣Σi∣1/2exp⁡(12(X−μi)TΣi−1(X−μi))(2)G_i(X) = \frac{1}{(2\pi)^{d/2} |\Sigma_i|^{1/2}} \exp \left( \frac{1}{2} (X-\mu_i)^T \Sigma_i^{-1}(X-\mu_i)\right) \tag{2}Gi(X)=(2π)d/2∣Σi∣1/21exp(21(X−μi)TΣi−1(X−μi))(2)

以二维高斯分布为例。当两个高斯分布有重叠时，会形成鞍状。如上图的d和e，二维高斯分布混合时会出现两个峰和一个鞍部；当两个分布几乎完全混合时，鞍部可能消失，但峰还在，此时明显的峰只有一个，如上图中的f。

论文中的两个高斯分布的OLR定义如下：
OLR(G1,G2)={1if p(X)has one peakp(Xsaddle)p(Xsubmax)if p(X)has two peaks(3)OLR(G_1, G_2) = \begin{cases} 1 &\text{if $p(X)$ has one peak} \\ \frac{p(X_{saddle})}{p(X_{submax})} &\text{if $p(X)$ has two peaks} \end{cases} \tag{3} OLR(G1,G2)={1p(Xsubmax)p(Xsaddle)if p(X) has one peakif p(X) has two peaks(3)

其中XsaddleX_{saddle}Xsaddle是pdf中的鞍点（saddle point），XsubmaxX_{submax}Xsubmax是pdf中的较低的峰（lower peak point）。OLR的示意图如下图所示。OLR计算的是鞍点的pdf与较低峰的pdf的比值。这么做是因为鞍点的pdf与混合系数αi\alpha_iαi有关。注意到OLR并不是落在重叠区域内数据的比例，因此跟数据量无关，只跟数据的分布有关。定义中的p(Xsubmax)p(X_{submax})p(Xsubmax)容易求，只需将两个均值带入(1)式，取较小的值即可。但是p(Xsaddle)p(X_{saddle})p(Xsaddle)不容直接求得。下面介绍如何计算p(Xsaddle)p(X_{saddle})p(Xsaddle)。

注意到两个峰点和鞍点在整个曲面上都应该是极值点。因此XsaddleX_{saddle}Xsaddle和XsubmaxX_{submax}Xsubmax应该满足下式：

{∂p∂x1=Ax1α1G1+Bx1α2G2∂p∂x2=Ax2α1G1+Bx2α2G2(4)\begin{cases} \frac{\partial p}{\partial x_1} = A_{x_1}\alpha_1G_1 + B_{x_1}\alpha_2 G_2 \\ \frac{\partial p}{\partial x_2} = A_{x_2}\alpha_1G_1 + B_{x_2}\alpha_2 G_2 \end{cases} \tag{4} {∂x1∂p=Ax1α1G1+Bx1α2G2∂x2∂p=Ax2α1G1+Bx2α2G2(4)

其中，

(Ax1Ax2)=∇∣∣X−μ1∣∣Σ1−12=−Σ1−1(X−μ1)(Bx1Bx2)=∇∣∣X−μ2∣∣Σ1−12=−Σ1−1(X−μ2)(5)\begin{aligned} \left( \begin{matrix} A_{x_1} \\ A_{x_2} \end{matrix} \right) = \nabla ||X-\mu_1||_{\Sigma_1^{-1}}^2 = -\Sigma_1^{-1}(X - \mu_1) \\ \left( \begin{matrix} B_{x_1} \\ B_{x_2} \end{matrix} \right) = \nabla ||X-\mu_2||_{\Sigma_1^{-1}}^2 = -\Sigma_1^{-1}(X - \mu_2) \end{aligned} \tag{5} (Ax1Ax2)=∇∣∣X−μ1∣∣Σ1−12=−Σ1−1(X−μ1)(Bx1Bx2)=∇∣∣X−μ2∣∣Σ1−12=−Σ1−1(X−μ2)(5)

(4)式式一条曲线。如果XXX已知，(5)式可求出。论文接下来证明，峰点和鞍点会在同一条曲线上，曲线方程如下：
Ax1Bx2−Bx1Ax2=0(6)A_{x_1} B_{x_2} - B_{x_1} A_{x_2} = 0 \tag{6}Ax1Bx2−Bx1Ax2=0(6)

而且，鞍点会在以两个峰点（均值处的pdf）之间的曲线段上。因此只要从第一个均值开始，沿着曲线(4)一直找到另一个均值，这个过程中的极小值点就是鞍点。得到鞍点的坐标，带入(1)式，就可以求得鞍点的pdf值。(6)式中的曲线称为Ridge Curve (RC).

OLR的算法如下：

输入混合高斯分布的参数(μ1,μ2,Σ1,Σ2,α1,α2)(\mu_1, \mu_2, \Sigma_1, \Sigma_2, \alpha_1, \alpha_2)(μ1,μ2,Σ1,Σ2,α1,α2)
计算RC： Ax1Bx2−Bx1Ax2=0A_{x_1} B_{x_2} - B_{x_1} A_{x_2} = 0Ax1Bx2−Bx1Ax2=0
沿着RC，从μ1\mu_1μ1到μ2\mu_2μ2按步长δ\deltaδ找到RC中p(X)p(X)p(X)取得最大值和最小值的点
3.1 令X0=μ1X_0 = \mu_1X0=μ1，X0X_0X0的下一个点Xi+1X_{i+1}Xi+1的第一维（x坐标）Xi+11={Xi+δ(μ1−μ2)}1X_{i+1}^1 = \{X_i + \delta(\mu_1-\mu_2)\}^1Xi+11={Xi+δ(μ1−μ2)}1.
3.2 将Xi+11X_{i+1}^1Xi+11带入RC方程(6)，求得Xi+1X_{i+1}Xi+1的第二维（y坐标）Xi+12X_{i+1}^2Xi+12
3.3 根据(1)式计算p(Xi)p(X_i)p(Xi)
3.4 if p(Xi)−p(Xi−1)>0p(X_i) - p(X_{i-1}) > 0p(Xi)−p(Xi−1)>0 and p(Xi)−p(Xi+1)>0p(X_i) - p(X_{i+1}) > 0p(Xi)−p(Xi+1)>0, XiX_iXi is maximum point (peak)
3.5 if p(Xi)−p(Xi−1)<0p(X_i) - p(X_{i-1}) < 0p(Xi)−p(Xi−1)<0 and p(Xi)−p(Xi+1)<0p(X_i) - p(X_{i+1}) < 0p(Xi)−p(Xi+1)<0, XiX_iXi is minimum point
根据(3)式计算OLR

上述算法δ\deltaδ可以取δ=∣∣μ1−μ2∣∣/1000\delta = ||\mu_1 - \mu_2|| / 1000δ=∣∣μ1−μ2∣∣/1000. 作者认为当OLR小于0.6时，两个类别可分性良好（visually well separated），当OLR大于0.8时，两个类别严重重叠（strongly overlapping）。

如果是有多个类别的情况，可以计算所有任意两个类别的重叠度，最后对所有重叠度求均值作为整体的重叠度。

算法实现

求p(Xi)p(X_i)p(Xi)可以用python第三方统计包scipy.stats中的multivariate_normal计算。输入两个高斯分布的参数可以求出pdf值。

完整代码可以参考GMM Overlap Rate。论文中给出的算法有一些问题。

import math
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import inv
from scipy.stats import multivariate_normalclass BiGauss(object):"""docstring for BiGauss"""def __init__(self, mu1, mu2, Sigma1, Sigma2, pi1, pi2, steps = 100):super(BiGauss, self).__init__()self.mu1      = mu1self.mu2      = mu2self.Sigma1   = Sigma1self.Sigma2   = Sigma2self.pi1      = pi1self.pi2      = pi2self.biGauss1 = multivariate_normal(mean = self.mu1, cov = self.Sigma1, allow_singular = True)self.biGauss2 = multivariate_normal(mean = self.mu2, cov = self.Sigma2, allow_singular = True)self.steps    = stepsself.inv_Sig1 = -inv(self.Sigma1)self.inv_Sig2 = -inv(self.Sigma2)# variables to calculate RCself.A_1 = self.inv_Sig1[0][0]self.B_1 = self.inv_Sig1[0][1]self.C_1 = self.inv_Sig1[1][0]self.D_1 = self.inv_Sig1[1][1]self.A_2 = self.inv_Sig2[0][0]self.B_2 = self.inv_Sig2[0][1]self.C_2 = self.inv_Sig2[1][0]self.D_2 = self.inv_Sig2[1][1]

计算pdf

def pdf(self, x):return self.pi1 * self.biGauss1.pdf(x) + self.pi2 * self.biGauss2.pdf(x)

根据xxx求出yyy，使得(x,y)(x,y)(x,y)在RC上

def RC(self, x):E = self.A_1 * (x - self.mu1[0])F = self.C_1 * (x - self.mu1[0])G = self.A_2 * (x - self.mu2[0])H = self.C_2 * (x - self.mu2[0])I = E * self.D_2 - F * self.B_2J = H * self.B_1 - G * self.D_1K = self.B_1 * self.D_2 - self.B_2 * self.D_1M = F * G - E * HP = KQ = I + J - K * (self.mu2[1] + self.mu1[1])S = -(M + I * self.mu2[1] + J * self.mu1[1])if Q**2 - 4*P*S < 0:return Noney = max((-Q + math.sqrt(Q**2 - 4*P*S)) / (2*P), (-Q - math.sqrt(Q**2 - 4*P*S)) / (2*P))return y

求OLR

def OLR(self):e      = math.sqrt((self.mu1[0] - self.mu2[0])**2 + (self.mu1[1] - self.mu2[1])**2) / float(self.steps)x_step = e*(self.mu1[0]-self.mu2[0]) # each step for xy_step = e*(self.mu1[1]-self.mu2[1]) # each step for yp_x    = self.mu1[0] - x_stepwhile self.RC(p_x) == None:p_x = p_x - x_stepp_y   = self.RC(p_x)p     = [p_x, p_y]p_pre = self.mu1p_min = min(self.pdf(p), self.pdf(p_pre))p_max = max(self.pdf(p), self.pdf(p_pre))index = 0while index < self.steps:if self.RC(p[0] - x_step) != None:p_next = [p[0] - x_step, self.RC(p[0] - x_step)] # next point on ridge curveif self.pdf(p) > self.pdf(p_pre) and self.pdf(p) > self.pdf(p_next):p_max = self.pdf(p)if self.pdf(p) < self.pdf(p_pre) and self.pdf(p) < self.pdf(p_next):p_min = self.pdf(p)p_pre = pp     = p_nextindex += 1pdf_mu1 = self.pdf(self.mu1)pdf_mu2 = self.pdf(self.mu2)return p_min / min(pdf_mu1, pdf_mu2) if p_min < min(pdf_mu1, pdf_mu2) else 1.0

上述代码有时会计算出OLR大于1的情况，还没有分析原因，因为草稿丢了，不知道代码中的A~S变量代表什么意思……因此代码中做了限制，如果求出的OLR大于1，那么只会返回1.

论文中探讨了混合系数、均值间距离和协方差对OLR的影响。论文中给出了一个例子，如下。当α1=0.46\alpha_1 = 0.46α1=0.46时该例子可以取到最小的OLR rminr_{min}rmin. 论文没有给出rminr_{min}rmin的具体数值，但是给出了OLR随α1\alpha_1α1取值变化的曲线图。上述代码算出来的结果是rmin=0.660r_{min} = 0.660rmin=0.660，也确实在α1=0.46\alpha_1 = 0.46α1=0.46处取得。与曲线图中的位置吻合。论文中提到当α1=0.3\alpha_1 = 0.3α1=0.3时，ORL等于0.7288，上述代码给出的结果是0.7270.

示例代码中也画出了OLR随α1\alpha_1α1变化的曲线图和OLR随两个均值之间距离变化的曲线图。曲线走势与论文中的图示一致，但具体数值有些差别。

这个示例代码只能计算二维混合高斯模型，更高维的无法计算，但是理论上，这个算法是适用于任何维度的GMM的。

计算高斯混合模型的可分性和重叠度（Overlap Rate, OLR）相关推荐

高斯混合模型（GMM）及其EM算法的理解
一个例子高斯混合模型(Gaussian Mixed Model)指的是多个高斯分布函数的线性组合,理论上GMM可以拟合出任意类型的分布,通常用于解决同一集合下的数据包含多个不同的分布的情况(或者是同 ...
高斯-赛得尔迭代式 c++_高斯混合模型(Gaussian Mixture Model)与EM算法原理(一)
高斯混合模型(Gaussian Mixture Model)是机器学习中一种常用的聚类算法,本文介绍了其原理,并推导了其参数估计的过程.主要参考Christopher M. Bishop的<Pa ...
机器学习（四）高斯混合模型
高斯混合模型聚类算法原文地址:http://blog.csdn.net/hjimce/article/details/45244603 作者:hjimce 高斯混合算法是EM算法的一个典型的应用,E ...
python混合高斯分布_python 高斯混合模型
高斯混合模型高斯混合模型(Gaussian Mixed Model)是使用多个(多元)高斯分布函数实现数据分类对目的. 在介绍高斯混合模型之前,我们需要了解以下几个概念: 一元高斯分布函数协方差矩 ...
声纹识别-2.GMM-UBM（高斯混合模型-通用背景模型）
声纹识别-2.GMM-UBM(高斯混合模型-通用背景模型) 前言声纹识别-1.绪论中回顾了声纹识别的类别,性能评价指标和算法.本篇博文介绍声纹识别算法中较为传统的GMM-UBM(Gaussian M ...
【数据挖掘】高斯混合模型 ( 与 K-Means 每个步骤对比 | 初始参数设置 | 计算概率 | 计算平均值参数 | 计算方差参数 | 计算高斯分布概率参数 | 算法终止条件 )
文章目录 I . 高斯混合模型 ( 样本 -> 模型 ) II . 高斯混合模型 ( 模型 -> 样本 ) III . 高斯混合模型与 K-Means 迭代过程对比 IV . 高斯混合模 ...
K-means与高斯混合模型
K-means http://blog.pluskid.org/?p=17 Clustering 中文翻译作"聚类",简单地说就是把相似的东西分到一组,同 Classificati ...
其他算法-高斯混合模型
目录高斯模型单高斯模型高斯混合模型GMM 参数估计单高斯模型参数估计-极大似然高斯混合模型参数估计-EM算法 GMM与k-means 高斯模型单高斯模型当样本数据 x∈Rx\in\mat ...
使用高斯混合模型对不同的股票市场状况进行聚类
来源:DeepHub IMBA 本文约2300字,建议阅读5分钟本文为你介绍如何将 GMM 应用于金融市场和经济. 介绍通过过去的十年的发展,普通人越来越容易进入股票市场,每天进出市场的资金量创历 ...

计算高斯混合模型的可分性和重叠度（Overlap Rate, OLR）

文章目录

简介

OLR计算

算法实现

计算高斯混合模型的可分性和重叠度（Overlap Rate, OLR）相关推荐

最新文章

热门文章