一种精简的GMM模型实现

sklearn提供了一种好用的GMM实现，但是它对修改不友好。

from sklearn.mixture import GaussianMixture

这里提供了一种简化的GMM实现，在保留原来风格的基础上实现了类似的功能

class GaussianMixture:
    """ Simplified Gaussian Mixture Model for direct use without sklearn dependency."""

    def __init__(self, n_components=1, covariance_type='full', tol=1e-3,
                 reg_covar=1e-6, max_iter=100, n_init=1, init_params='kmeans',
                 random_state=None):
        self.n_components = n_components
        self.covariance_type = covariance_type
        self.tol = tol
        self.reg_covar = reg_covar
        self.max_iter = max_iter
        self.n_init = n_init
        self.init_params = init_params
        self.random_state = random_state
        self.weights_ = None
        self.means_ = None
        self.covariances_ = None
        self.precisions_cholesky_ = None
        self.converged_ = False
        self.lower_bound_ = -np.inf
        self.n_iter_ = 0
        self.random_state_ = np.random.RandomState(self.random_state)

    def fit(self, X):
        """Fit the Gaussian mixture model to the data X using EM algorithm."""
        self._initialize_parameters(X)
        log_likelihood = -np.inf
        for n_iter in range(self.max_iter):
            prev_log_likelihood = log_likelihood
            log_resp = self._e_step(X)
            self._m_step(X, log_resp)
            log_likelihood = self._compute_lower_bound(log_resp)
            change = log_likelihood - prev_log_likelihood
            if abs(change) < self.tol:
                self.converged_ = True
                break
        self.n_iter_ = n_iter + 1
        self.lower_bound_ = log_likelihood
        return self

    def _initialize_parameters(self, X):
        n_samples, n_features = X.shape
        if self.init_params == 'kmeans':
            from sklearn.cluster import KMeans
            labels = KMeans(n_clusters=self.n_components, n_init=1, random_state=self.random_state).fit(X).labels_
            resp = np.zeros((n_samples, self.n_components))
            resp[np.arange(n_samples), labels] = 1
        else:
            resp = self.random_state_.rand(n_samples, self.n_components)
            resp /= resp.sum(axis=1, keepdims=True)
        self._initialize(X, resp)

    def _initialize(self, X, resp):
        nk = resp.sum(axis=0) + 10 * np.finfo(float).eps
        means = np.dot(resp.T, X) / nk[:, np.newaxis]
        covariances = self._estimate_covariances(resp, X, nk, means)
        weights = nk / nk.sum()
        self.weights_ = weights
        self.means_ = means
        self.covariances_ = covariances
        self.precisions_cholesky_ = self._compute_precision_cholesky(covariances)

    def _e_step(self, X):
        weighted_log_prob = self._estimate_weighted_log_prob(X)
        log_prob_norm = logsumexp(weighted_log_prob, axis=1)
        log_resp = weighted_log_prob - log_prob_norm[:, np.newaxis]
        return log_resp

    def _m_step(self, X, log_resp):
        resp = np.exp(log_resp)
        nk = resp.sum(axis=0) + 10 * np.finfo(float).eps
        means = np.dot(resp.T, X) / nk[:, np.newaxis]
        covariances = self._estimate_covariances(resp, X, nk, means)
        weights = nk / nk.sum()
        self.weights_ = weights
        self.means_ = means
        self.covariances_ = covariances
        self.precisions_cholesky_ = self._compute_precision_cholesky(covariances)

    def _estimate_covariances(self, resp, X, nk, means):
        n_samples, n_features = X.shape
        covariances = np.empty((self.n_components, n_features, n_features))
        for k in range(self.n_components):
            diff = X - means[k]
            covariances[k] = np.dot(resp[:, k] * diff.T, diff) / nk[k]
            covariances[k] += self.reg_covar * np.eye(n_features)
        return covariances

    def _compute_precision_cholesky(self, covariances):
        precisions_cholesky = np.empty_like(covariances)
        for k, cov in enumerate(covariances):
            try:
                prec_chol = linalg.cholesky(cov, lower=True)
                precisions_cholesky[k] = linalg.inv(prec_chol).T
            except linalg.LinAlgError:
                raise ValueError("Covariance matrix is singular.")
        return precisions_cholesky

    def _estimate_weighted_log_prob(self, X):
        return self._estimate_log_prob(X) + np.log(self.weights_)

    def _estimate_log_prob(self, X):
        n_samples = X.shape[0]
        log_prob = np.empty((n_samples, self.n_components))
        for k, (mean, prec_chol) in enumerate(zip(self.means_, self.precisions_cholesky_)):
            y = np.dot(X, prec_chol) - np.dot(mean, prec_chol)
            log_prob[:, k] = np.sum(np.square(y), axis=1)
        return -0.5 * (np.log(2 * np.pi) + log_prob)

    def _compute_lower_bound(self, log_resp):
        return log_resp.sum()

    def score_samples(self, X):
        """Calculate the weighted log probabilities for each sample."""
        log_prob = self._estimate_weighted_log_prob(X)
        return logsumexp(log_prob, axis=1)

    def predict_proba(self, X):
        """Calculate the probabilities for each sample belonging to each component."""
        weighted_log_prob = self._estimate_weighted_log_prob(X)
        log_prob_norm = logsumexp(weighted_log_prob, axis=1)
        return np.exp(weighted_log_prob - log_prob_norm[:, np.newaxis])