Gaussian Mixture Model

• 算法特征:
  ①. 高斯分布作为基函数; ②. 多个高斯分布进行凸组合; ③. 极大似然法估计概率密度.
• 算法推导:
  GMM概率密度形式如下:
  egin{equation}
  p(x) = sum_{k=1}^{K}pi_kN(x|mu_k, Sigma_k)
  label{eq_1}
  end{equation}
  其中, $pi_k$、$mu_k$、$Sigma_k$分别表示第$k$个高斯分布的权重、均值及协方差矩阵, 且$sumlimits_{k=1}^{K}pi_k = 1, forall pi_k geq 0$.
  令样本集合为${x^{(1)}, x^{(2)}, cdots, x^{(n)}}$, 本文拟采用EM(Expectation-Maximization)算法求解上述优化变量${pi_k, mu_k, Sigma_k}_{k=1sim K}$.
  $step1$:
  随机初始化${pi_k, mu_k, Sigma_k}_{k=1sim K}$.
  $step2 sim E step$:
  计算第$i$个样本落在第$k$个高斯的概率:
  egin{equation}
  gamma_k^{(i)} = frac{pi_kN(x^{(i)}|mu_k, Sigma_k)}{sumlimits_{k=1}^{K} pi_k N(x^{(i)}|mu_k, Sigma_k)}
  label{eq_2}
  end{equation}
  $step3 sim M step$:
  计算第$k$个高斯的样本数:
  egin{equation}
  N_k = sum_{i=1}^{n}gamma_k^{(i)}
  label{eq_3}
  end{equation}
  更新第$k$个高斯的权重:
  egin{equation}
  pi_k = frac{N_k}{N}
  label{eq_4}
  end{equation}
  更新第$k$个高斯的均值:
  egin{equation}
  mu_k = frac{sumlimits_{i=1}^{n}gamma_k^{(i)}x^{(i)}}{N_k}
  label{eq_5}
  end{equation}
  更新第$k$个高斯的协方差矩阵:
  egin{equation}
  Sigma_k = frac{sumlimits_{i=1}^{n}gamma_k^{(i)}(x^{(i)} - mu_k)(x^{(i)} - mu_k)^{mathrm{T}}}{N_k}
  label{eq_6}
  end{equation}
  $step4$:
  回到$step2$, 直到优化变量${pi_k, mu_k, Sigma_k}_{k=1sim K}$均收敛.
• 代码实现:
  Part Ⅰ:
  现以如下数据集为例进行算法实施:
  

  # 生成聚类数据集
  
  import numpy
  from matplotlib import pyplot as plt
  
  
  numpy.random.seed(3)
  
  
  def generate_data(clusterNum):
      mu = [0, 0]
      sigma = [[0.03, 0], [0, 0.03]]
      data = numpy.random.multivariate_normal(mu, sigma, (1000, ))
      
      for idx in range(clusterNum  - 1):
          mu = numpy.random.uniform(-1, 1, (2, ))
          arr = numpy.random.uniform(0, 1, (2, 2))
          sigma = numpy.matmul(arr.T, arr) / 10
          tmpData = numpy.random.multivariate_normal(mu, sigma, (1000, ))
          data = numpy.vstack((data, tmpData))
      
      return data
  
  
  def plot_data(data):
      fig = plt.figure(figsize=(5, 3))
      ax1 = plt.subplot()
      
      ax1.scatter(data[:, 0], data[:, 1], s=1)
      ax1.set(xlim=[-2, 2], ylim=[-2, 2])
      # plt.show()
      fig.savefig("plot_data.png", dpi=100)
      plt.close()
  
          
  
  if __name__ == "__main__":
      X = generate_data(5)
      plot_data(X)

  

  Part Ⅱ:
  GMM实现如下:
  

  # Gaussian Mixture Model之实现
  
  import numpy
  from matplotlib import pyplot as plt
  
  
  numpy.random.seed(3)
  
  
  
  def generate_data(clusterNum):
      mu = [0, 0]
      sigma = [[0.03, 0], [0, 0.03]]
      data = numpy.random.multivariate_normal(mu, sigma, (1000, ))
      
      for idx in range(clusterNum  - 1):
          mu = numpy.random.uniform(-1, 1, (2, ))
          arr = numpy.random.uniform(0, 1, (2, 2))
          sigma = numpy.matmul(arr.T, arr) / 10
          tmpData = numpy.random.multivariate_normal(mu, sigma, (1000, ))
          data = numpy.vstack((data, tmpData))
      
      return data
  
  
  
  class GMM(object):
      
      def __init__(self, X, clusterNum):
          self.__X = X                                      # 聚类数据集, 1行代表1个样本
          self.__clusterNum = clusterNum                    # 类簇数量
          
          
      def get_clusters(self, pi, mu, sigma):
          '''
          获取类簇
          '''
          X = self.__X
          clusterNum = self.__clusterNum
          gamma = self.__calc_gamma(X, pi, mu, sigma)
          maxIdx = numpy.argmax(gamma, axis=0)
          
          clusters = list(list() for idx in range(clusterNum))
          for idx, x in enumerate(X):
              clusters[maxIdx[idx]].append(x)
          clusters = list(numpy.array(cluster) for cluster in clusters)
          return clusters
          
          
      def PDF(self, x, pi, mu, sigma):
          '''
          GMM概率密度函数
          '''
          val = sum(list(self.__calc_gaussian(x, mu[k], sigma[k]) * pi[k] for k in range(mu.shape[0])))
          return val
          
          
      def optimize(self, maxIter=1000, epsilon=1.e-6):
          '''
          maxIter: 最大迭代次数
          epsilon: 收敛判据, 优化变量趋于稳定则收敛
          '''
          X = self.__X
          clusterNum = self.__clusterNum    # 簇数量
          N = X.shape[0]                    # 样本数量
          
          pi, mu, sigma = self.__init_GMMOptVars(X, clusterNum)
          gamma = numpy.zeros((clusterNum, N))
          for idx in range(maxIter):
              print("iterCnt: {:3d},   pi = {}".format(idx, pi))
              gamma = self.__calc_gamma(X, pi, mu, sigma)
              
              Nk = numpy.sum(gamma, axis=1)
              piNew = Nk / N
              muNew = self.__calc_mu(X, gamma, Nk)
              sigmaNew = self.__calc_sigma(X, gamma, mu, Nk)
              
              deltaPi = piNew - pi
              deltaMu = muNew - mu
              deltaSigma = sigmaNew - sigma
              if self.__converged(deltaPi, deltaMu, deltaSigma, epsilon):
                  return piNew, muNew, sigmaNew, True
              pi, mu, sigma = piNew, muNew, sigmaNew
              
          return pi, mu, sigma, False
              
              
      def __calc_sigma(self, X, gamma, mu, Nk):
          sigma = list()
          for gamma_k, mu_k, Nk_k in zip(gamma, mu, Nk):
              term1 = X - mu_k
              term2 = gamma_k.reshape((-1, 1)) * term1
              term3 = numpy.matmul(term2.T, term1)
              sigma_k = term3 / Nk_k
              sigma.append(sigma_k)
          return numpy.array(sigma)
              
              
      def __calc_mu(self, X, gamma, Nk):
          mu = list()
          for gamma_k, Nk_k in zip(gamma, Nk):
              term1 = gamma_k.reshape((-1, 1)) * X
              term2 = numpy.sum(term1, axis=0)
              mu_k = term2 / Nk_k
              mu.append(mu_k)
          return numpy.array(mu)
              
              
      def __calc_gamma(self, X, pi, mu, sigma):
          gamma = numpy.zeros((pi.shape[0], X.shape[0]))
          for k in range(gamma.shape[0] - 1):
              for i in range(gamma.shape[1]):
                  gamma[k, i] = self.__calc_gamma_ik(X[i], k, pi, mu, sigma)
          term1 = numpy.sum(gamma[:-1, :], axis=0)
          gamma[-1] = 1 - term1
          return gamma
          
              
      def __converged(self, deltaPi, deltaMu, deltaSigma, epsilon):
          val1 = numpy.linalg.norm(deltaPi)
          val2 = numpy.linalg.norm(deltaMu)
          val3 = numpy.linalg.norm(deltaSigma)
          if val1 <= epsilon and val2 <= epsilon and val3 <= epsilon:
              return True
          return False
          
          
      def __calc_gamma_ik(self, x_i, k, pi, mu, sigma):
          pi_k = pi[k]
          mu_k = mu[k]
          sigma_k = sigma[k]
          
          upperVal = pi_k * self.__calc_gaussian(x_i, mu_k, sigma_k)
          lowerVal = self.PDF(x_i, pi, mu, sigma)
          gamma_ik = upperVal / lowerVal
          return gamma_ik
      
          
      def __calc_gaussian(self, x, mu, sigma):
          '''
          x: 输入x - 1维数组
          mu: 均值
          sigma: 协方差矩阵
          '''
          d = mu.shape[0]
          term0 = (x - mu).reshape((-1, 1))
          term1 = 1 / numpy.math.sqrt((2 * numpy.math.pi) ** d * numpy.linalg.det(sigma))
          term2 = numpy.matmul(numpy.matmul(term0.T, numpy.linalg.inv(sigma)), term0)[0, 0]
          term3 = numpy.math.exp(-term2 / 2)
          val = term1 * term3
          return val
          
          
      def __init_GMMOptVars(self, X, clusterNum):
          pi = self.__init_pi(clusterNum)                         # GMM权重初始化
          mu = self.__init_mu(X, clusterNum)                      # GMM均值初始化
          sigma = self.__init_sigma(X.shape[1], clusterNum)       # GMM协方差矩阵初始化 - 采用单位矩阵初始化之
          return pi, mu, sigma
          
          
      def __init_sigma(self, n, clusterNum):
          sigma = list()
          for idx in range(clusterNum):
              sigma.append(numpy.identity(n))
          return numpy.array(sigma)
          
          
      def __init_mu(self, X, clusterNum, epsilon=1.e-5):
          '''
          采用K-means方法进行初始化
          '''
          lowerBound = numpy.min(X, axis=0)
          upperBound = numpy.max(X, axis=0)
          oriMu = numpy.random.uniform(lowerBound, upperBound, (clusterNum, X.shape[1]))
          traMu = self.__updateMuByKMeans(X, oriMu)
          while numpy.linalg.norm(traMu - oriMu) > epsilon:
              oriMu = traMu
              traMu = self.__updateMuByKMeans(X, oriMu)
          return traMu    
          
              
      def __updateMuByKMeans(self, X, oriMu):
          distList = list()
          for mu in oriMu:
              distTerm = numpy.linalg.norm(X - mu, axis=1)
              distList.append(distTerm)
          distArr = numpy.vstack(distList)
          distIdx = numpy.argmin(distArr, axis=0)
          
          clusLst = list(list() for i in range(oriMu.shape[0]))
          for xIdx, dIdx in enumerate(distIdx):
              clusLst[dIdx].append(X[xIdx])
          
          traMu = list()
          for clusIdx, clus in enumerate(clusLst):
              if len(clus) == 0:
                  traMu.append(oriMu[clusIdx])
              else:
                  tmpClus = numpy.array(clus)
                  mu = numpy.sum(tmpClus, axis=0) / tmpClus.shape[0]
                  traMu.append(mu)
          return numpy.array(traMu)
              
          
          
      def __init_pi(self, clusterNum):
          pi = numpy.ones((clusterNum, )) / clusterNum
          return pi
          
          
          
  class GMMPlot(object):
      
      @staticmethod
      def profile_plot(X, gmmObj, pi, mu, sigma):
          surfX1 = numpy.linspace(-2, 2, 100)
          surfX2 = numpy.linspace(-2, 2, 100)
          surfX1, surfX2 = numpy.meshgrid(surfX1, surfX2)
          
          surfY = numpy.zeros((surfX2.shape[0], surfX1.shape[1]))
          for rowIdx in range(surfY.shape[0]):
              for colIdx in range(surfY.shape[1]):
                  surfY[rowIdx, colIdx] = gmmObj.PDF(numpy.array([surfX1[rowIdx, colIdx], surfX2[rowIdx, colIdx]]), pi, mu, sigma)
          
          clusters = gmmObj.get_clusters(pi, mu, sigma)
          
          fig = plt.figure(figsize=(10, 3))
          ax1 = plt.subplot(1, 2, 1)
          ax2 = plt.subplot(1, 2, 2)
          
          ax1.contour(surfX1, surfX2, surfY, levels=16)
          ax1.scatter(X[:, 0], X[:, 1], marker="o", color="tab:blue", s=1)
          ax1.set(xlim=[-2, 2], ylim=[-2, 2], xlabel="$x_1$", ylabel="$x_2$")
          
          for idx, cluster in enumerate(clusters):
              ax2.scatter(cluster[:, 0], cluster[:, 1], s=1, label="$cluster {}$".format(idx))
          ax2.set(xlim=[-2, 2], ylim=[-2, 2], xlabel="$x_1$", ylabel="$x_2$")
          ax2.legend()
          
          fig.tight_layout()
          # plt.show()
          fig.savefig("profile_plot.png", dpi=100)
          plt.close()
          
          
          
  if __name__ == "__main__":
      clusterNum = 5
      X = generate_data(5)
      
      gmmObj = GMM(X, clusterNum)
      pi, mu, sigma, _ = gmmObj.optimize()
      
      GMMPlot.profile_plot(X, gmmObj, pi, mu, sigma)

• 结果展示:
  左侧为聚类轮廓, 右侧为簇划分. 很显然, 高斯混合模型适用于高斯型数据分布, 此时聚类效果较为明显.
• 使用建议:
  ①. 由于EM算法只能保证局部最优, 因此选择一个合适的初值较为重要, 本文以K-means结果作为初值进行优化;
  ②. 根据基函数特征, 高斯混合模型适用于高斯型数据分布.
• 参考文档:
  https://baike.baidu.com/item/%E9%AB%98%E6%96%AF%E6%B7%B7%E5%90%88%E6%A8%A1%E5%9E%8B/8878468?fr=aladdin