聚类-混合高斯模型（GMM）

EM算法

参考：
http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html
http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006924.html
http://blog.csdn.net/gugugujiawei/article/details/45583051

给定的训练样本是{}，样例间独立，我们想找到每个样例隐含的类别z，使得p(x,z)最大。
p(x,z)的最大似然估计如下：

第一步是对极大似然取对数，第二步是对每个样例的每个可能类别z求联合分布概率和。但是直接求一般比较困难，因为有隐藏变量z存在，但是一般确定了z后，求解就容易了。

EM是一种解决存在隐含变量优化问题的有效方法。虽然不能直接最大化，我们可以不断地建立L的下确界（E步），然后优化下界（M步），但是一般确定了z，求解就容易了。

对于每一个样例i，让Qi表示该样例隐含变量z的某种分布，Qi满足的条件是，
根据Jenson不等式，得到如下公式：

求当等式成立时的Qi（zi）。
可以看做是对求了下界。
可求出公式：。
至此，推出了在固定后，每个样本属于各个z（分布）的概率。这一步即是E步，接下来M步，在给定的每个样本属于各个分布的情况下，调整，去极大化 的下界。

Gmm的代码示例如下：

from numpy import *

def loadDataSet(fileName):
    dataSet = []
    fr = open(fileName)
    for line in fr.readlines():
        curLine = line.strip().split('\t')
        fltLine = list(map(float, curLine))
        dataSet.append(fltLine)
    return dataSet


class GMM:

    def __init__(self, k=4, eps=0.00001):
        self.k = k
        self.eps = eps

    def fit_EM(self, X, max_iters=1000):

        # n = number of data-points, d = dimension of data points
        n, d = X.shape

        # randomly choose the starting centroids/means
        ## as 3 of the points from datasets
        mu = X[np.random.choice(n, self.k, False), :]

        # initialize the covariance matrices for each gaussians
        Sigma = [np.eye(d)] * self.k

        # initialize the probabilities/weights for each gaussians
        w = [1. / self.k] * self.k

        # responsibility matrix is initialized to all zeros
        # we have responsibility for each of n points for eack of k gaussians
        R = np.zeros((n, self.k))

        ### log_likelihoods
        log_likelihoods = []

        P = lambda mu, s: np.linalg.det(s) ** -.5 * (2 * np.pi) ** (-d / 2) \
                          * np.exp(-.5 * np.einsum('ij, ij -> i', \
                                                   X - mu, np.dot(np.linalg.inv(s), (X - mu).T).T))

        # Iterate till max_iters iterations
        while len(log_likelihoods) < max_iters:

            # E - Step

            ## Vectorized implementation of e-step equation to calculate the
            ## membership for each of k -gaussians
            for k in range(self.k):
                a = w[k] * P(mu[k], Sigma[k])
                R[:, k] = a

            ### Likelihood computation
            log_likelihood = np.sum(np.log(np.sum(R, axis=1)))

            log_likelihoods.append(log_likelihood)

            ## Normalize so that the responsibility matrix is row stochastic
            R = (R.T / np.sum(R, axis=1)).T

            ## The number of datapoints belonging to each gaussian
            N_ks = np.sum(R, axis=0)

            # M Step
            ## calculate the new mean and covariance for each gaussian by
            ## utilizing the new responsibilities
            for k in range(self.k):
                ## means
                mu[k] = 1. / N_ks[k] * np.sum(R[:, k] * X.T, axis=1).T
                x_mu = np.matrix(X - mu[k])

                ## covariances
                Sigma[k] = np.array(1 / N_ks[k] * np.dot(np.multiply(x_mu.T, R[:, k]), x_mu))

                ## and finally the probabilities
                w[k] = 1. / n * N_ks[k]
            # check for convergence
            if len(log_likelihoods) < 2: continue
            if np.abs(log_likelihood - log_likelihoods[-2]) < self.eps: break

        ## bind all results together
        from collections import namedtuple
        self.params = namedtuple('params', ['mu', 'Sigma', 'w', 'log_likelihoods', 'num_iters'])
        self.params.mu = mu
        self.params.Sigma = Sigma
        self.params.w = w
        self.params.log_likelihoods = log_likelihoods
        self.params.num_iters = len(log_likelihoods)

        return self.params

    def plot_log_likelihood(self):
        import pylab as plt
        plt.plot(self.params.log_likelihoods)
        plt.title('Log Likelihood vs iteration plot')
        plt.xlabel('Iterations')
        plt.ylabel('log likelihood')
        plt.show()

    def predict(self, x):
        p = lambda mu, s: np.linalg.det(s) ** - 0.5 * (2 * np.pi) ** \
                                                      (-len(x) / 2) * np.exp(-0.5 * np.dot(x - mu, \
                                                                                           np.dot(np.linalg.inv(s),
                                                                                                  x - mu)))
        probs = np.array([w * p(mu, s) for mu, s, w in \
                          zip(self.params.mu, self.params.Sigma, self.params.w)])
        return probs / np.sum(probs)