em算法实例正态分布

看了很多介绍EM算法的文章，但是他们都没有代码，所以在这里写出来。

Jensen 不等式

参考期望最大算法

Jensen不等式在优化理论中大量用到，首先来回顾下凸函数和凹函数的定义。假设

是定义域为实数的函数，如果对于所有的

，

的二阶导数大于等于0，那么

是凸函数。当

是向量时，如果hessian矩阵

是半正定（即

是凸函数）。如果，

的二阶导数小于0或者

就是凹函数。

Jensen不等式描述如下：

如果
是凸函数，
是随机变量，则
是严格凸函数时，则
如果
是凹函数，
是随机变量，则
，当
是（严格）凹函数当且仅当
是（严格）凸函

EM思想

极大似然函数法估计参数的一般步骤：

写出似然函数
取对数
求导数，并令导数为0
解似然方程

给定

个训练样本

，假设样本之间相互独立，要拟合模型

。根据分布我们可以得到如下的似然函数：

需要对每个样本实例的每个可能的类别

求联合概率分布之和，即

。

如果

是已知的，那么使用极大似然估计参数

会很容易。

然而上式存在一个不确定的隐含变量(latent random variable)

，这种情况下EM算法就派上用场了。

由于不能直接最大化

，所以只能不断地建立

的下界(E-step)，再优化下界。一直迭代直到算法收敛到局部最优解。

EM算法通过引入隐含变量，使用MLE(极大似然估计)进行迭代求解参数。通常引入隐含变量后会有两个参数，EM算法首先会固定其中一个参数，然后使用MLE计算第二个参数；然后固定第二个参数，再使用MLE估计第一个参数的值。依次迭代，直到收敛到局部最优解。

E-Step: 通过观察到的状态和现有模型估计参数估计值隐含状态
M-Step: 假设隐含状态已知的情况下，最大化似然函数。

由于算法保证了每次迭代之后，似然函数都会增加，所以函数最终会收敛

EM推导

对于每个实例

，用

表示样本实例隐含变量

的某种分布，且

满足条件

，如果

是连续的，则

表示概率密度函数，将求和换成积分。

上式最后的变换用到了Jensen不等式：

由于

函数的二阶导数为

，为凹函数，所以使用

；

把所以上式写成

，那么我们可以通过不断的优化

的下界，来使得

不断提高，最终达到它的最大值。

在Jensen不等式中，当

，即为常数时，等号成立。在这里即为：

变换并对

求和得到：

因为

，概率之和为1，所以：

因此：

可以看出，固定了参数

之后，使下界拉升的

的计算公式就是后验概率，一并解决了

如何选择的问题。

EM完整的流程如下：

初始化参数分布
重复E-Step和M-Step直到收敛
1. E-Step: 根据参数的初始值或者上一次迭代的模型参数来计算出隐含变量的后验概率，其实就是隐含变量的期望值，作为隐含变量的当前估计值:
2. M-Step: 最大化似然函数从而获得新的参数值:

多维高斯分布

一元高斯分布的概率密度函数为：

因为

是标量，所以

等价于

，所以上式等价于

推广到多维得到多元高斯分布，得到K维随机变量

的概率密度函数：

和

都是K维向量

是协方差阵的行列式，协方差阵

是

的正定矩阵，称

服从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)