看了很多介紹EM算法的文章，但是他們都沒有代碼，所以在這里寫出來。

Jensen 不等式

參考期望最大算法

Jensen不等式在優(yōu)化理論中大量用到，首先來回顧下凸函數(shù)和凹函數(shù)的定義。假設(shè)

是定義域?yàn)閷?shí)數(shù)的函數(shù)，如果對于所有的 ， 的二階導(dǎo)數(shù)大于等于0，那么 是凸函數(shù)。當(dāng) 是向量時，如果hessian矩陣 是半正定（即 是凸函數(shù)）。如果， 的二階導(dǎo)數(shù)小于0或者 就是凹函數(shù)。

Jensen不等式描述如下：

如果 是凸函數(shù)， 是隨機(jī)變量，則 是嚴(yán)格凸函數(shù)時，則

如果 是凹函數(shù)， 是隨機(jī)變量，則 ，當(dāng) 是（嚴(yán)格）凹函數(shù)當(dāng)且僅當(dāng) 是（嚴(yán)格）凸函

EM思想

極大似然函數(shù)法估計(jì)參數(shù)的一般步驟：

寫出似然函數(shù)

取對數(shù)

求導(dǎo)數(shù)，并令導(dǎo)數(shù)為0

解似然方程

給定

個訓(xùn)練樣本 ，假設(shè)樣本之間相互獨(dú)立，要擬合模型 。根據(jù)分布我們可以得到如下的似然函數(shù)：

需要對每個樣本實(shí)例的每個可能的類別

求聯(lián)合概率分布之和，即 。

如果

是已知的，那么使用極大似然估計(jì)參數(shù) 會很容易。

然而上式存在一個不確定的隱含變量(latent random variable)

，這種情況下EM算法就派上用場了。

由于不能直接最大化

，所以只能不斷地建立 的下界(E-step)，再優(yōu)化下界。一直迭代直到算法收斂到局部最優(yōu)解。

EM算法通過引入隱含變量，使用MLE(極大似然估計(jì))進(jìn)行迭代求解參數(shù)。通常引入隱含變量后會有兩個參數(shù)，EM算法首先會固定其中一個參數(shù)，然后使用MLE計(jì)算第二個參數(shù)；然后固定第二個參數(shù)，再使用MLE估計(jì)第一個參數(shù)的值。依次迭代，直到收斂到局部最優(yōu)解。

• E-Step: 通過觀察到的狀態(tài)和現(xiàn)有模型估計(jì)參數(shù)估計(jì)值 隱含狀態(tài)
• M-Step: 假設(shè)隱含狀態(tài)已知的情況下，最大化似然函數(shù)。

由于算法保證了每次迭代之后，似然函數(shù)都會增加，所以函數(shù)最終會收斂

EM推導(dǎo)

對于每個實(shí)例

，用 表示樣本實(shí)例 隱含變量 的某種分布，且 滿足條件 ，如果 是連續(xù)的，則 表示概率密度函數(shù)，將求和換成積分。

上式最后的變換用到了Jensen不等式：

由于

函數(shù)的二階導(dǎo)數(shù)為 ，為凹函數(shù)，所以使用 ；

把所以上式寫成

，那么我們可以通過不斷的優(yōu)化 的下界，來使得 不斷提高，最終達(dá)到它的最大值。

在Jensen不等式中，當(dāng)

，即為常數(shù)時，等號成立。在這里即為：

變換并對

求和得到：

因?yàn)?

，概率之和為1，所以：

因此：

可以看出，固定了參數(shù)

之后，使下界拉升的 的計(jì)算公式就是后驗(yàn)概率，一并解決了 如何選擇的問題。

EM完整的流程如下：

初始化參數(shù)分布

重復(fù)E-Step和M-Step直到收斂

E-Step: 根據(jù)參數(shù)的初始值或者上一次迭代的模型參數(shù)來計(jì)算出隱含變量的后驗(yàn)概率，其實(shí)就是隱含變量的期望值，作為隱含變量的當(dāng)前估計(jì)值:

M-Step: 最大化似然函數(shù)從而獲得新的參數(shù)值:

多維高斯分布

一元高斯分布的概率密度函數(shù)為：

因?yàn)?

是標(biāo)量，所以 等價(jià)于 ，所以上式等價(jià)于

推廣到多維得到多元高斯分布，得到K維隨機(jī)變量

的概率密度函數(shù)：

和 都是K維向量 是協(xié)方差陣的行列式，協(xié)方差陣 是 的正定矩陣，稱 服從K元正態(tài)分布，簡記為：

多元高斯分布的極大似然估計(jì)

對于

個樣本 ，其似然函數(shù)為：

分別對

和 求偏導(dǎo)，參考多元正態(tài)分布的極大似然估計(jì)得到:

# 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

生成二維分布需要乘以協(xié)方差矩陣（協(xié)方差矩陣是正定的，所以可以分解(Cholesky)成下三角矩陣),

二維高斯分布的參數(shù)分析

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)

兩個二元高斯分布混合的分布

定義分量數(shù)目 ，對每個分量設(shè)置 ，然后計(jì)算下式的對數(shù)似然函數(shù)

2. E-Step，根據(jù)當(dāng)前的

計(jì)算后驗(yàn)概率 ， 是先驗(yàn)概率， 表示點(diǎn) 屬于聚類 的后驗(yàn)概率。

3. M-Step，更新

4. 檢查是否收斂，否則轉(zhuǎn)#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)

二元高斯分布的概率密度函數(shù)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)

對數(shù)似然函數(shù)的收斂過程

參考

WIKI多元正態(tài)分布

期望最大算法

二維高斯分布的參數(shù)分析

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。