import numpy as npclass 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()