看了很多介绍EM算法的文章，但是他们都没有代码，所以在这里写出来。
Jensen 不等式
参考期望最大算法
Jensen不等式在优化理论中大量用到，首先来回顾下凸函数和凹函数的定义。假设
Jensen不等式描述如下：
- 如果
是凸函数，
是随机变量，则
是严格凸函数时，则
- 如果
是凹函数，
是随机变量，则
，当
是（严格）凹函数当且仅当
是（严格）凸函
EM思想
极大似然函数法估计参数的一般步骤：
- 写出似然函数
- 取对数
- 求导数，并令导数为0
- 解似然方程
给定
需要对每个样本实例的每个可能的类别
如果
然而上式存在一个不确定的隐含变量(latent random variable)
由于不能直接最大化
EM算法通过引入隐含变量，使用MLE(极大似然估计)进行迭代求解参数。通常引入隐含变量后会有两个参数，EM算法首先会固定其中一个参数，然后使用MLE计算第二个参数；然后固定第二个参数，再使用MLE估计第一个参数的值。依次迭代，直到收敛到局部最优解。
- E-Step: 通过观察到的状态和现有模型估计参数估计值 隐含状态
- M-Step: 假设隐含状态已知的情况下，最大化似然函数。
由于算法保证了每次迭代之后，似然函数都会增加，所以函数最终会收敛
EM推导
对于每个实例
上式最后的变换用到了Jensen不等式：
由于
把所以上式写成
在Jensen不等式中，当
变换并对
因为
因此：
可以看出，固定了参数
EM完整的流程如下：
- 初始化参数分布
- 重复E-Step和M-Step直到收敛
- E-Step: 根据参数的初始值或者上一次迭代的模型参数来计算出隐含变量的后验概率，其实就是隐含变量的期望值，作为隐含变量的当前估计值:
- M-Step: 最大化似然函数从而获得新的参数值:
- E-Step: 根据参数的初始值或者上一次迭代的模型参数来计算出隐含变量的后验概率，其实就是隐含变量的期望值，作为隐含变量的当前估计值:
多维高斯分布
一元高斯分布的概率密度函数为：
因为
推广到多维得到多元高斯分布，得到K维随机变量
多元高斯分布的极大似然估计
对于
分别对
# use multivariate_normal to generate 2d gaussian distribution
mean = [3, 4]
cov = [[1.5, 0], [0, 3.3]]
x = np.random.multivariate_normal(mean, cov, 500)
plt.scatter(x[:, 0], x[:, 1])mu_hat = np.mean(x, axis=0)
print(mu_hat)
sigma_hat = ((x-mu_hat).T @ (x-mu_hat)) / 500
print(sigma_hat)
#[2.89991371 4.08421214]
#[[ 1.43340175 -0.01134683]#[-0.01134683 3.28850441]]
高斯混合模型
生成一维的高斯分布
sigma * np.random.randn(...) + mu
生成二维分布需要乘以协方差矩阵（协方差矩阵是正定的，所以可以分解(Cholesky)成下三角矩阵),
二维高斯分布的参数分析
from scipy.stats import multivariate_normaldef gen_gaussian(conv, mean, num=1000):points = np.random.randn(num, 2)points = points @ conv + meanreturn pointsfig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111)conv, mean = np.array([[1, 0], [0, 5]]), np.array([2, 4])
points1 = gen_gaussian(conv, mean)
plt.scatter(points1[:, 0], points1[:, 1])conv, mean = np.array([[2, 0], [0, 3]]), np.array([10, 15])
points2 = gen_gaussian(conv, mean)
plt.scatter(points2[:, 0], points2[:, 1])points = np.append(points1, points2, axis=0)K = 2
X = points
mu = np.array([[2, 4], [10, 15]])
cov = np.array([[[1, 0], [0, 5]], [[2, 0], [0, 3]]])x, y = np.meshgrid(np.sort(X[:, 0]), np.sort(X[:, 1]))
XY = np.array([x.flatten(), y.flatten()]).T
reg_cov = 1e-6 * np.identity(2)
for m, c in zip(mu, cov):c = c + reg_covmng = multivariate_normal(mean=m, cov=c)ax.contour(np.sort(X[:, 0]), np.sort(X[:, 1]), mng.pdf(XY).reshape(len(X), len(X)), colors='black', alpha=0.3)
- 定义分量数目
，对每个分量设置
，然后计算下式的对数似然函数
2. E-Step，根据当前的
3. M-Step，更新
4. 检查是否收敛，否则转#2
# https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.multivariate_normal.htmlimport math# implement my own gaussian pdf, get same result as multivariate_normal.pdf
def gaussian(X, K, mu, cov):p = 1.0 / np.sqrt(np.power(2*math.pi, K) * np.linalg.det(cov))i = (X-mu).T @ np.linalg.inv(cov) @ (X-mu)p *= np.power(np.e, -0.5*i)return pX = np.random.rand(2,)
mu = np.random.rand(2, )
cov = np.array([[1, 0], [0, 3]])mng = multivariate_normal(mean=mu, cov=cov)
gaussian(X, 2, mu, cov), mng.pdf(X)输出 (0.08438545576275427, 0.08438545576275427)
# This import registers the 3D projection, but is otherwise unused.
from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as npfig = plt.figure()
ax = fig.gca(projection='3d')# Make data.
X = np.arange(-5, 5, 0.25)
Y = np.arange(-5, 5, 0.25)
X, Y = np.meshgrid(X, Y)
Z = np.array([gaussian(x, 2, mu, cov) for x in zip(X.flatten(), Y.flatten())]).reshape(X.shape)print(X.shape, Y.shape, Z.shape)
# Plot the surface.
surf = ax.plot_surface(X, Y, Z, cmap=plt.cm.coolwarm,linewidth=0, antialiased=False)
from tqdm import tqdmX = points
N, K, D = len(points), 2, len(X[0])
mu = np.random.randint(min(X[:,0]),max(X[:,0]),size=(K, D))
d = np.max(X)
rcov = 1e-6*np.identity(D)
cov = np.zeros((K, D, D))
for dim in range(K):np.fill_diagonal(cov[dim], d)pi0 = np.random.rand()
pi = np.array([pi0, 1-pi0])rnk = np.zeros((N, K))
muh, covh, Rh = [], [], []log_likelihoods = []for i in tqdm(range(100)):muh.append(mu)covh.append(cov)# E-Steprnk = np.zeros((N, K)) for m, co, p, k in zip(mu, cov, pi, range(K)):co = co + rcovmng = multivariate_normal(mean=m, cov=co)d = np.sum([pi_k * multivariate_normal(mean=mu_k, cov=cov_k).pdf(X) for pi_k, mu_k, cov_k in zip(pi, mu, cov+rcov)], axis=0)rnk[:, k] = p * mng.pdf(X) / dRh.append(rnk)# for n in range(N):
# d = sum([pi[k]*gaussian(X[n], K, mu[k], cov[k]) for k in range(K)])
# for k in range(K):
# rnk[n, k] = pi[k] * gaussian(X[n], K, mu[k], cov[k]) / d# M-Stepmu, cov, pi = np.zeros((K, D)), np.zeros((K, D, D)), np.zeros((K, 1)) for k in range(K):nk = np.sum(rnk[:, k], axis=0)# new meanmuk = (1/nk) * np.sum(X*rnk[:, k].reshape(N, 1), axis=0)mu[k] = muk# new conv matrixcovk = (rnk[:, k].reshape(N, 1) * (X-muk)).T @ (X-muk) + rcovcov[k] = covk / nk# new pipi[k] = nk / np.sum(rnk)log_likelihoods.append(np.log(np.sum([p*multivariate_normal(mu[i], cov[k]).pdf(X) for p, i, k in zip(pi, range(len(X[0])), range(K))])))plt.plot(log_likelihoods, label='log_likelihoods')
plt.legend()fig = plt.figure()
ax = fig.add_subplot(111)
plt.scatter(points1[:, 0], points1[:, 1])
plt.scatter(points2[:, 0], points2[:, 1])for m, c in zip(mu, cov):mng = multivariate_normal(mean=m,cov=c)ax.contour(np.sort(X[:, 0]), np.sort(X[:, 1]), mng.pdf(XY).reshape(len(X), len(X)),colors='black',alpha=0.3)
参考
WIKI多元正态分布
期望最大算法
二维高斯分布的参数分析
