一维GMM的Python代码实现

import numpy as np


class GMM(object):
    """Gaussian Mixture Model 
    """

    def __init__(self, data, K):
        """
        K: the number of gaussian models
        alpha: the weight for corresponding gaussian model
        mu: the vector of means
        sigma2: the vector of variances
        N: the number of samples
        K: the number of gaussian models
        """

        self.data = data
        self.K = K
        self.alpha = (np.ones(K) / K)
        self.mu = (np.arange(K) - K // 2) * (data.max() - data.min()) / K
        self.sigma2 = (np.ones(K))
        self.N = len(data)
        self.gamma = np.ones((self.N, K)) / K

    def phi(self):
        # phi.shape(K, N)
        phi = (1 / np.sqrt(2 * np.pi * self.sigma2.reshape(self.K, 1)) * np.exp(- (self.data - self.mu.reshape(self.K, 1)) ** 2 / (2 * self.sigma2.reshape(self.K, 1))))
        return phi

    def fit(self):
        sigma2_ = self.sigma2
        mu_ = self.mu
        while True:
            # gamma.shape(N, K)
            self.gamma = (0.1 * self.gamma + 0.9 * self.phi().T * self.alpha / (self.phi().T * self.alpha).sum(axis=1).reshape(self.N, 1))
            # mu.shape(1, K)
            self.mu = (0.1 * self.mu + 0.9 * np.matmul(self.data, self.gamma) / self.gamma.sum(axis=0))
            # sigma2.shape(1,K)
            self.sigma2 = (0.1 * self.sigma2 + 0.9 * (self.gamma * (data.reshape(self.N, 1) - self.mu) ** 2).sum(axis=0) / self.gamma.sum(axis=0))
            # alpha.shape(1, K)
            self.alpha = (0.1 * self.alpha + 0.9 * self.gamma.sum(axis=0) / self.N)
            if (np.sum((self.mu - mu_) ** 2) + np.abs(self.sigma2 - sigma2_).sum()) < 10 ** (-10):
                break

            mu_ = self.mu
            sigma2_ = self.sigma2
            print(self.gamma.argmax(axis=1))

        return self.gamma.argmax(axis=1)


data = np.concatenate((np.random.normal(-4, 1, 2000), np.random.normal(4, 1, 2000)))
gmm = GMM(data, 2)
label = gmm.fit()