GMM高斯混合模型EM算法代码实现 与 sklearn的比较

因为要用到GMM的EM算法，所以自己写了一个（本人很菜，这一点代码写了两天）。写完发现sklearn的包里有（心情复杂）。。那不如检测一下自己写的对不对吧。嘿嘿。做了下实现，用同一组参数初始化，看下来是没问题的。

关于GMM 的理论，有两篇博客写的很好，放下链接：
高斯混合模型（GMM）及其EM算法的理解
详解EM算法与混合高斯模型(Gaussian mixture model, GMM)

import numpy as np
from numpy import linalg
import math
from sklearn.mixture import GaussianMixture


class GMM(object):
    def __init__(self, X, K):
        self.N, self.D = X.shape
        self.K = K
        self.Mu, self.Sigma, self.Pi, self.Gamma = self.init_params()
    
    def init_params(self):
        Mu = np.random.rand(self.K, self.D)  # K*D
        Sigma = np.array([np.eye(self.D)] * self.K) * 0.1  # K*D*D
        Pi = np.array([1.0 / self.K] * self.K)  # K
        Gamma = np.random.rand(self.N, self.K)
        return Mu, Sigma, Pi, Gamma
    
    def e_step(self, X):
        prob = []
        for i in range(self.K):
            prob.append(get_Gaussian_prob(X, self.Mu[i], self.Sigma[i]))
        prob = np.array(prob).T
        Pi_Gussian = prob * self.Pi
        sum_K = np.sum(Pi_Gussian, axis=1)
        return (Pi_Gussian.T / sum_K).T  # Gamma
    
    def m_step(self, X):
        NK = np.sum(self.Gamma, axis=0)
        _1_NK = 1 / NK
        # get new Mu
        m = np.dot(self.Gamma.T, X)
        newMu = (_1_NK * m.T).T
        # get new Sigma
        newSigma = []
        for i in range(self.K):
            newSigma.append(_1_NK[i] * np.dot((self.Gamma[:, i] * (X - newMu[i]).T), (X - newMu[i])))
        newSigma = np.array(newSigma)
        # get Pi
        newPi = NK / self.N
        return newMu, newSigma, newPi
    
    def train(self, X, times):
        for i in range(times):
            self.Gamma = self.e_step(X)
            self.Mu, self.Sigma, self.Pi = self.m_step(X)
        return self.Mu, self.Sigma, self.Pi, self.Gamma


def get_Gaussian_prob(X, mu, sigma):
    """已知一个高斯分布的参数，求其概率"""
    D = len(sigma)
    if X.size == D:
        n = 1
    else:
        n, _ = X.shape
    det = linalg.det(sigma)
    
    norm_const = 1.0 / (math.pow(2 * math.pi, float(D) / 2) * math.pow(det, 1.0 / 2))
    x_mu = np.array(X - mu)
    inv = linalg.inv(sigma)
    result = np.power(math.e, -0.5 * np.dot(np.dot(x_mu, inv), x_mu.T))
    if n != 1:
        result = np.diag(result)
    return norm_const * result

把自己写的高斯混合模型EM算法 和 用sklearn包里的高斯混合模型API一起，使用同一组参数，用来确定自己写的代码没问题的代码在GitHub