EM （Expectation Maximization）算法是一種迭代算法，由 Dempster 等人于1977年提出。EM 算法最大的價(jià)值在于解決含隱變量的概率模型參數(shù)的極大似然估計(jì)問題，通過 E步更新后驗(yàn)概率，通過 M步求參數(shù)極大似然估計(jì)，EM算法常用于無監(jiān)督學(xué)習(xí)，如聚類問題。
EM 算法的難點(diǎn)在于要理解目標(biāo)函數(shù)的解并非問題的最終解，這一點(diǎn)與許多經(jīng)典機(jī)器學(xué)習(xí)算法不同，問題的最終解是隱變量 z 的后驗(yàn)分布 Q 分布，而目標(biāo)函數(shù)的解是參數(shù) theta 和 P 分布，theta 往往用正態(tài)分布的參數(shù)表示。Q 分布無法通過傳統(tǒng)的極大似然估計(jì)求解，所以只能先通過目標(biāo)函數(shù)求 theta 和 P，再間接的求出 Q，更新 Q 之后再刷新目標(biāo)函數(shù)更新 theta 和 P，循環(huán)往復(fù)直至收斂。所以，theta 和 P 就像一組橋梁，連接了目標(biāo)函數(shù)與問題假設(shè)。
0、先導(dǎo)知識：Jeson 不等式
1、數(shù)據(jù)集與特征空間
2、目標(biāo)函數(shù)（極大似然估計(jì)）
此刻，如果能使 >= 變?yōu)?= ，那么目標(biāo)函數(shù)就明確了，因此繼續(xù)往下推：
此時(shí)，等號得以成立，更新目標(biāo)函數(shù)：
此時(shí)，可見目標(biāo)函數(shù)是關(guān)于 mu, sigma, P 的函數(shù)，P 是每個(gè)隱變量所占的權(quán)重，因此 P 的和 = 1。下面自然就是求極大似然估計(jì)：
3、假設(shè)空間
EM 算法的難點(diǎn)在于要理解目標(biāo)函數(shù)的解并非問題的最終解，這一點(diǎn)與許多經(jīng)典機(jī)器學(xué)習(xí)算法不同，問題的最終解是隱變量 z 的后驗(yàn)分布 Q 分布，而目標(biāo)函數(shù)的解是參數(shù) theta 和 P 分布，theta 往往用正態(tài)分布的參數(shù)表示。Q 分布無法通過傳統(tǒng)的極大似然估計(jì)求解，所以只能先通過目標(biāo)函數(shù)求 theta 和 P，再間接的求出 Q，更新 Q 之后再刷新目標(biāo)函數(shù)更新 theta 和 P，循環(huán)往復(fù)直至收斂。所以，theta 和 P 就像一組橋梁，連接了目標(biāo)函數(shù)與問題假設(shè)。
問題的假設(shè)可以理解為： Q 分布是關(guān)于 theta 和 P 分布的函數(shù)，Q* 是問題的最終解，也是假設(shè)函數(shù)的最優(yōu)值；theta* 和 P* 是目標(biāo)函數(shù)的最優(yōu)解，也是假設(shè)函數(shù)的最優(yōu)解。
4、優(yōu)化步驟
4.1、參數(shù)初始化（EM 對初始值敏感）
4.2、更新 Q 分布（E步）
4.3、更新 mu, sigma, P（M步）
4.4、循環(huán)4.2、4.3，直到誤差足夠小，算法停止
4.5、得到最終的 Q 分布，取其大者為隱變量最終的類別標(biāo)記
5、算法收斂性證明
手寫 GMM 并與 Sklearn 對比
# -*- coding: utf-8 -*-
import os
import copy
#import random as rd
import numpy as np
import pandas as pd
#from sklearn.utils import shuffle
#from sklearn.model_selection import train_test_split
#from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture
from scipy.stats import norm
from sklearn.metrics import silhouette_score
class GMM(object):
def __init__(self, z, tol, max_iter, init_params, verbose):
self.z = z
self.tol = tol
self.max_iter = max_iter
self.init_params = init_params
self.verbose = verbose
def fit(self, x):
'''模型訓(xùn)練'''
# 參數(shù)初始化
N, n = x.shape
Q = np.zeros([N, self.z])
pdf = np.ones([N, self.z])
P = np.random.uniform(low=0, high=1, size=(self.z, 1))
if self.init_params == "kmeans":
km = KMeans(n_clusters=self.z, init="k-means++")
km.fit(x)
mu = km.cluster_centers_
elif self.init_params == "random":
mu = np.random.randn(self.z, n) * 0.01
else:
raise ValueError("init_params must be 'kmeans' or 'random'")
std = np.std(x, axis=0)
std = np.tile(std, self.z).reshape(self.z, n)
mu_update = np.zeros_like(mu)
std_update = np.ones_like(std)
err_mu_res = []
err_std_res = []
for epoch in range(self.max_iter):
# E
for j in range(self.z):
# 隱變量外循環(huán)
for col in range(n):
# 列內(nèi)循環(huán)，計(jì)算聯(lián)合分布
pdf[:, j] *= norm.pdf(x[:, col], mu[j, col], std[j, col])
Q[:, j] = pdf[:, j] * P[j]
# 計(jì)算Q分布
Q = Q / np.sum(Q, axis=1).reshape(-1, 1)
# M
for j in range(self.z):
# 更新每個(gè)隱變量的參數(shù)值
mu_update[j] = x.T.dot(Q[:, j]) / np.sum(Q[:, j])
std_update[j] = np.sqrt(Q[:, j].dot((x - mu[j])**2) / np.sum(Q[:, j]))
P[j] = np.mean(Q[:, j])
# 停止條件
err_mu = np.sum(np.abs(mu - mu_update))
err_std = np.sum(np.abs(std - std_update))
if self.verbose:
print(f"epoch {epoch} err_mu {err_mu} err_std {err_std}")
err_mu_res.append(err_mu)
err_std_res.append(err_std)
if err_mu <= self.tol and err_std <= self.tol:
break
else:
mu = copy.deepcopy(mu_update)
std = copy.deepcopy(std_update)
return Q, P, err_mu_res, err_std_res
def predict(self, Q):
'''判斷類別'''
labels = np.argmax(Q, axis=1)
return labels
if __name__ == "__main__":
file_path = os.getcwd()
dataSet = pd.read_csv(file_path + "/swiss.csv")
x = dataSet[["Fertility", "Agriculture", "Catholic", "InfantMortality"]].values
# 手寫模型
model = GMM(z=3, tol=0.001, max_iter=1000, init_params="kmeans", verbose=True)
Q, P, err_mu_res, err_std_res = model.fit(x)
labels = model.predict(Q)
score = silhouette_score(x, labels, metric="euclidean")
print(f"手寫模型輪廓系數(shù) {score}")
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(err_mu_res)
plt.xlabel("epoch")
plt.ylabel("error")
plt.title("mu")
plt.show()
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(err_std_res)
plt.xlabel("epoch")
plt.ylabel("error")
plt.title("std")
plt.show()
# sklearn
gmm = GaussianMixture(n_components=4, tol=0.001, max_iter=1000, init_params="kmeans")
gmm.fit(x)
gmm.predict_proba(x)
labels = gmm.predict(x)
print(f"sklearn 模型輪廓系數(shù) {score}")
