高斯混合模型

混合模型是潜变量模型的一种，是最常见的形式之一。而高斯混合模型(Gaussian Mixture Models, GMM)是混合模型中最常见的一种。zzz代表该数据点是由某一个高斯分布产生的。π\piπ在这里是指该点属于哪一个高斯分布的先验概率。除次之外，我们还需要找到每一个高斯分布的参数，即均值和协方差矩阵。

p(x)=∑k=1Kπkpk(x)(1)p(x)=∑k=1KπkN(x∣μk,Σk)(2)p(x) = \sum_{k=1}^K \pi_k p_k(x) \qquad \qquad (1)\\\\ p(x) = \sum_{k=1}^K \pi_k \mathcal N(x| \mu_k, \Sigma_k) \qquad (2) p(x)=k=1∑K​πk​pk​(x)(1)p(x)=k=1∑K​πk​N(x∣μk​,Σk​)(2)

我们对混合模型的一般形式即(1)进行拓展，已知每一种分布单独的期望和协方差矩阵，求出xxx的期望和协方差矩阵。

E[x]=∑xxp(x)=∑xx∑kπkpk(x)=∑kπk∑xxpk(x)=∑kπkE[pk(x)]=∑kπkμkE[x] = \sum_x x p(x) \\\\ = \sum_x x \sum_k \pi_k p_k(x) \\\\ = \sum_k \pi_k \sum_x x p_k(x)\\\\ = \sum_k \pi_k E[p_k(x)] \\\\ = \sum_k \pi_k \mu_k E[x]=x∑​xp(x)=x∑​xk∑​πk​pk​(x)=k∑​πk​x∑​xpk​(x)=k∑​πk​E[pk​(x)]=k∑​πk​μk​

Σx=E[x2]−(E[x])2=E[xxT]−(∑kπkμk)2=∫xxTp(x)dx−(∑kπkμk)2=∫xxT∑kπkpk(x∣k)dx−(∑kπkμk)2=∑kπk∫xxTpk(x∣k)dx−(∑kπkμk)2=∑kπkE[xxT∣k]−(∑kπkμk)2=∑kπk(E[xxT∣k]−μkμkT+μkμkT)−(∑kπkμk)2=∑kπk(Σk+μkμkT)−(∑kπkμk)2\Sigma_x = E[x^2] - (E[x])^2 \\\\ = E[xx^T] - (\sum_k \pi_k \mu_k)^2 \\\\ = \int xx^T p(x) dx - (\sum_k \pi_k \mu_k)^2 \\\\ = \int xx^T \sum_k \pi_k p_k(x|k) dx - (\sum_k \pi_k \mu_k)^2 \\\\ = \sum_k \pi_k \int xx^T p_k(x|k) dx - (\sum_k \pi_k \mu_k)^2 \\\\ = \sum_k \pi_k E[xx^T|k] - (\sum_k \pi_k \mu_k)^2 \\\\ = \sum_k \pi_k (E[xx^T|k] - \mu_k\mu_k^T + \mu_k\mu_k^T) - (\sum_k \pi_k \mu_k)^2 \\\\ = \sum_k \pi_k (\Sigma_k + \mu_k\mu_k^T) - (\sum_k \pi_k \mu_k)^2 \\\\ Σx​=E[x2]−(E[x])2=E[xxT]−(k∑​πk​μk​)2=∫xxTp(x)dx−(k∑​πk​μk​)2=∫xxTk∑​πk​pk​(x∣k)dx−(k∑​πk​μk​)2=k∑​πk​∫xxTpk​(x∣k)dx−(k∑​πk​μk​)2=k∑​πk​E[xxT∣k]−(k∑​πk​μk​)2=k∑​πk​(E[xxT∣k]−μk​μkT​+μk​μkT​)−(k∑​πk​μk​)2=k∑​πk​(Σk​+μk​μkT​)−(k∑​πk​μk​)2

适用场景

多分类。

优点

鲁棒性较好

缺点

需要数据服从分布（具有严格假设）参数较多，计算相对比较复杂

K-平均演算法

我们可以通过K-Means (KMeans)得到每个类的聚类中心μk\mu_kμk​。

过程

初始化聚类中心μk\mu_kμk​将每个数据点分配到离自己最近的聚类中心zi=argmink∣∣xi−μk∣∣22z_i = arg min_k ||x_i - \mu_k||_2^2zi​=argmink​∣∣xi​−μk​∣∣22​更新聚类中心μk=1Nk∑i:zi=kxi\mu_k = \frac 1 {N_k} \sum_{i:z_i=k} x_iμk​=Nk​1​∑i:zi​=k​xi​

由于KMeans更擅长发现凸类型簇，聚类结果有可能出现偏差。

最大期望算法

Expectation Maximization (EM)。在概率模型中寻找参数最大似然估计或者最大后验估计的算法，其中概率模型依赖于无法观测的隐变量。

在经过推导之后，可以发现当类别先验相同且方差为单位矩阵的时候，EM算法和KMeans实际上是一样的。

过程

初始化先验概率均等，协方差矩阵为单位矩阵(Identity Matrix)。

E-step:

如果是软分类(soft assignment)，计算每个点到每个聚类的概率(responsibility，注意不是probability)，

ri,k=πkN(x∣μk,Σk)∑k′πk′N(x∣μk′,Σk′)r_{i,k} = \frac{\pi_k \mathcal N(x| \mu_k, \Sigma_k)} {\sum_{k'} \pi_k' \mathcal N(x| \mu_k', \Sigma_k')} ri,k​=∑k′​πk′​N(x∣μk′​,Σk′​)πk​N(x∣μk​,Σk​)​

如果是硬分类(hard assignment)，直接分配到最有可能的中心，

zi=argmaxkπkN(x∣μk,Σk)∑k′πk′N(x∣μk′,Σk′)z_i = argmax_k \frac{\pi_k \mathcal N(x| \mu_k, \Sigma_k)} {\sum_{k'} \pi_k' \mathcal N(x| \mu_k', \Sigma_k')} zi​=argmaxk​∑k′​πk′​N(x∣μk′​,Σk′​)πk​N(x∣μk​,Σk​)​

M-step:

重新估计先验、聚类中心和协方差矩阵，

μk=1∑iri,k∑iri,kxiΣk=1∑iri,k∑iri,k(xi−μk)(xi−μk)Tπk=1N∑iri,k\mu_k = \frac 1 {\sum_i r_{i,k}} \sum_i r_{i,k}x_i \\\\ \Sigma_k = \frac 1 {\sum_i r_{i,k}} \sum_i r_{i,k} (x_i-\mu_k)(x_i-\mu_k)^T\\\\ \pi_k = \frac 1 N \sum_i r_{i,k} μk​=∑i​ri,k​1​i∑​ri,k​xi​Σk​=∑i​ri,k​1​i∑​ri,k​(xi​−μk​)(xi​−μk​)Tπk​=N1​i∑​ri,k​

目标函数

最大化对数似然。

logL(θ)=∑i=1nlogpθ(xi)=∑i=1nlog∑zipθ(xi∣zi)log L(\theta) = \sum_{i=1}^n log p_\theta(x_i) = \sum_{i=1}^n log \sum_{z_i}p_\theta(x_i|z_i) logL(θ)=i=1∑n​logpθ​(xi​)=i=1∑n​logzi​∑​pθ​(xi​∣zi​)

该函数有可能不是凸函数。

E-step: 对于每个数据点，在已知θt−1\theta_{t-1}θt−1​计算ziz_izi​的期望值

M-step: 对于ziz_izi​求出θt\theta_tθt​的最大似然

代码实现

使用python最基础的代码。

import numpy as np# specifying data pointsx0 = np.array([[-3], [1]])x1 = np.array([[-3], [-1]])x2 = np.array([[3], [-1]])x3 = np.array([[3], [1]])data = [x0,x1,x2,x3]n = len(data)d = len(x0) ## dimension# specifying initial conditionK = 2mu0_start = np.array([[-1], [0]])mu1_start = np.array([[1], [0]])mus_start = [mu0_start,mu1_start]Ident = np.eye(d)Sigmas_start = [Ident, Ident]priors_start = [1/K]*Kdef calc_responsibility(x_list, prior_list, mu_list, Sigma_list, i, k):# used to prevent singular matrixrand_mat = 1e-6*np.random.rand(K, K)x = x_list[i]mu = mu_list[k]Sigma = Sigma_list[k]prior = prior_list[k]r = prior*np.exp(-(x-mu).T@np.linalg.inv(Sigma+rand_mat)@(x-mu)/2)total = 0for itr in range(len(mu_list)):mu_prime = mu_list[itr]Sigma_prime = Sigma_list[itr] + rand_matprior_prime = prior_list[itr]total += prior_prime*np.exp(-(x-mu_prime).T@np.linalg.inv(Sigma_prime)@(x-mu_prime)/2)return (r/total)[0,0]def calc_mu(x_list, r_list):total = np.zeros((d,1))for i in range(len(x_list)):total += r_list[i]*x_list[i]total /= np.sum(r_list)return totaldef calc_cov_mat(x_list, r_list, mu_):size = len(r_list)res = np.zeros((K, K))for i in range(size):tmp = np.outer(x_list[i]-mu_,x_list[i]-mu_)res += tmp*r_list[i]return res / np.sum(r_list)def calc_prior(r_list):return sum(r_list)/ndef run_EM(data, priors_start, mus_start, Sigmas_start, iters=30, printing=False):priors_prev = priors_startmus_prev = mus_startSigmas_prev = Sigmas_startresponsibilities = np.zeros((K, n))priors = [0]*Kmus = [0]*KSigmas = [0]*Kfor itr in range(iters):# calculate responsbility for each data point for each classfor k in range(K):for i in range(n):responsibilities[k][i] = calc_responsibility(data, priors_prev, mus_prev, Sigmas_prev, i, k)# calculate class priorpriors = [calc_prior(responsibilities[k], n) for k in range(K)]# calculate class centermus = [calc_mu(data, responsibilities[k]) for k in range(K)]# calculate class covariance matrixSigmas = [calc_cov_mat(data, responsibilities[k], mus[k]) for k in range(K)]# update parameters priors_prev = priorsmus_prev = musSigmas_prev = Sigmasif printing:print("-----iteration {}-----".format(itr))print("Priors: {}".format(priors))print("Mean: ")for mu_ in mus:print(np.matrix(mu_))print("Covariance matrix: ")for Sigma_ in Sigmas:print(np.matrix(Sigma_))return priors, mus, Sigmas

Reference

Probability and Information Theory in Machine Learning (ECE 601), Matthew Malloy, Fall

如果觉得《机器学习基础专题：高斯混合模型和最大期望EM算法以及代码实现》对你有帮助，请点赞、收藏，并留下你的观点哦！