看了很多介紹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直到收斂

E-Step: 根據參數的初始值或者上一次迭代的模型參數來計算出隱含變量的后驗概率，其實就是隱含變量的期望值，作為隱含變量的當前估計值:

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)

對數似然函數的收斂過程

參考

WIKI多元正態分布

期望最大算法

二維高斯分布的參數分析

總結

以上是生活随笔為你收集整理的em算法 实例 正态分布_EM算法解GMM的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。