算法思想：含有隱變量的極大似然估計(jì)
我們經(jīng)常會(huì)從樣本觀察數(shù)據(jù)中烹棉，找出樣本的模型參數(shù)慷蠕。 最常用的方法就是極大化模型分布的對(duì)數(shù)似然函數(shù)掸冤。
但是在一些情況下恳守，我們得到的觀察數(shù)據(jù)有未觀察到的隱含數(shù)據(jù)，此時(shí)我們未知的有隱含數(shù)據(jù)和模型參數(shù)贩虾，因而無法直接用極大化對(duì)數(shù)似然函數(shù)得到模型分布的參數(shù)催烘。怎么辦呢？這就是EM算法可以派上用場(chǎng)的地方了缎罢。那么先復(fù)習(xí)一下極大似然估計(jì)伊群。
極大似然估計(jì)（MLE）
直接舉個(gè)例子：
某位同學(xué)與一位獵人一起外出打獵，一只野兔從前方竄過策精。只聽一聲槍響舰始，野兔應(yīng)聲到下，如果要你推測(cè)咽袜，這一發(fā)命中的子彈是誰打的丸卷？你就會(huì)想，只發(fā)一槍便打中询刹，由于獵人命中的概率一般大于這位同學(xué)命中的概率谜嫉，看來這一槍是獵人射中的。
這個(gè)例子所作的推斷就體現(xiàn)了極大似然法的基本思想凹联。
MLE就是利用已知的樣本結(jié)果沐兰，反推最有可能(最大概率)導(dǎo)致這樣結(jié)果的參數(shù)值 的計(jì)算過程。直白來講蔽挠，就是給定了一定的數(shù)據(jù)住闯，假定知道數(shù)據(jù)是從某種分布中 隨機(jī)抽取出來的，但是不知道這個(gè)分布具體的參數(shù)值，即“模型已定比原，參數(shù)未 知”插佛，MLE就可以用來估計(jì)模型的參數(shù)。MLE的目標(biāo)是找出一組參數(shù)(模型中的 參數(shù))量窘，使得模型產(chǎn)出觀察數(shù)據(jù)的概率最大朗涩。
求最大似然函數(shù)估計(jì)值的一般步驟：
（1）寫出似然函數(shù)；
（2）對(duì)似然函數(shù)取對(duì)數(shù)绑改，并整理；
（3）求導(dǎo)數(shù)兄一，令導(dǎo)數(shù)為0厘线，得到似然方程；
（4）解似然方程出革，得到的參數(shù)即為所求造壮；
EM算法解決這個(gè)的思路是使用啟發(fā)式的迭代方法，既然我們無法直接求出模型分布參數(shù)骂束，那么我們可以先猜想隱含數(shù)據(jù)（EM算法的E步）耳璧，接著基于觀察數(shù)據(jù)和猜測(cè)的隱含數(shù)據(jù)一起來極大化對(duì)數(shù)似然，求解我們的模型參數(shù)（EM算法的M步)展箱。由于我們之前的隱藏?cái)?shù)據(jù)是猜測(cè)的旨枯，所以此時(shí)得到的模型參數(shù)一般還不是我們想要的結(jié)果。不過沒關(guān)系混驰，我們基于當(dāng)前得到的模型參數(shù)攀隔，繼續(xù)猜測(cè)隱含數(shù)據(jù)（EM算法的E步），然后繼續(xù)極大化對(duì)數(shù)似然栖榨，求解我們的模型參數(shù)（EM算法的M步)昆汹。以此類推，不斷的迭代下去婴栽，直到模型分布參數(shù)基本無變化满粗，算法收斂，找到合適的模型參數(shù)愚争。
從上面的描述可以看出映皆，EM算法是迭代求解最大值的算法，同時(shí)算法在每一次迭代時(shí)分為兩步轰枝，E步和M步劫扒。一輪輪迭代更新隱含數(shù)據(jù)和模型分布參數(shù)，直到收斂狸膏，即得到我們需要的模型參數(shù)沟饥。
K-Means算法
下面回顧一下K-Means算法的基本思想。
K-means算法，也稱為k-均值聚類算法贤旷，是一種非常廣泛使用的聚類算法之一广料。
假定輸入樣本為S=x1,x2,x3,...,xm,則算法步驟為：
假設(shè)輸入樣本為T=X1,X2…,Xm則算法步驟為(使用歐幾里得距離公式)
選擇初始化的k個(gè)類別中心a1,...ak
對(duì)于每個(gè)樣本X,將其標(biāo)記位距離類別中心aj最近的類別
更新每個(gè)類別的中心點(diǎn)aj為隸屬該類別的所有樣本的均值
更新后的中心點(diǎn)為：
表示第j個(gè)簇中的樣本點(diǎn)的均值為新的中心點(diǎn)。
- 重復(fù)上面兩步操作,直到達(dá)到某個(gè)中止條件
中止條件：
迭代次數(shù)幼驶、最小平方誤差MSE艾杏、簇中心點(diǎn)變化率。
一個(gè)最直觀了解EM算法思路的是K-Means算法盅藻，見之前寫的K-Means聚類算法原理购桑。在K-Means聚類時(shí)，每個(gè)聚類簇的質(zhì)心是隱含數(shù)據(jù)氏淑。我們會(huì)假設(shè)KK個(gè)初始化質(zhì)心勃蜘，即EM算法的E步；然后計(jì)算得到每個(gè)樣本最近的質(zhì)心假残，并把樣本聚類到最近的這個(gè)質(zhì)心缭贡，即EM算法的M步。重復(fù)這個(gè)E步和M步辉懒，直到質(zhì)心不再變化為止阳惹，這樣就完成了K-Means聚類。
當(dāng)然眶俩，K-Means算法是比較簡(jiǎn)單的莹汤，實(shí)際中的問題往往沒有這么簡(jiǎn)單。上面對(duì)EM算法的描述還很粗糙颠印，我們需要用數(shù)學(xué)的語言精準(zhǔn)描述体啰。
EM算法
先從一個(gè)簡(jiǎn)單的例子開始：
現(xiàn)在一個(gè)班里有50個(gè)男生，50個(gè)女生嗽仪，且男生站左荒勇，女生站右。我們假定男生的身高服從正態(tài)分布μ1闻坚，σ1沽翔，女生的身高則服從另一個(gè)正態(tài)分布μ2，σ2 窿凤。這時(shí)候我們可以用極大似然法（MLE）仅偎，分別通過這50個(gè)男生和50個(gè)女生的樣本來估計(jì)這兩個(gè)正態(tài)分布的參數(shù)。
但現(xiàn)在我們讓情況復(fù)雜一點(diǎn)雳殊，就是這50個(gè)男生和50個(gè)女生混在一起了橘沥。我們擁有100個(gè)人的身高數(shù)據(jù)，卻不知道這100個(gè)人每一個(gè)是男生還是女生夯秃。
這時(shí)候情況就有點(diǎn)尷尬座咆，因?yàn)橥ǔ碚f痢艺，我們只有知道了精確的男女身高的正態(tài)分布參數(shù)我們才能知道每一個(gè)人更有可能是男生還是女生。但從另一方面去考量介陶，我們只有知道了每個(gè)人是男生還是女生才能盡可能準(zhǔn)確地估計(jì)男女各自身高的正態(tài)分布的參數(shù)堤舒。
下面以一個(gè)知乎上的例子來計(jì)算說明EM算法，地址：
https://www.zhihu.com/question/27976634/answer/39132183
背景：公司有很多領(lǐng)導(dǎo)=[A總哺呜，劉總舌缤，C總]，同時(shí)有很多漂亮的女職員=[小甲某残，小章国撵，小乙]。（請(qǐng)勿對(duì)號(hào)入座）你迫切的懷疑這些老總跟這些女職員有問題玻墅。為了科學(xué)的驗(yàn)證你的猜想介牙，你進(jìn)行了細(xì)致的觀察。于是椭豫，
觀察數(shù)據(jù)：
1）A總，小甲旨指，小乙一起出門了赏酥；
2）劉總，小甲谆构，小章一起出門了裸扶；
3）劉總，小章搬素，小乙一起出門了呵晨；
4）C總，小乙一起出門了熬尺；
收集到了數(shù)據(jù)摸屠，你開始了神秘的EM計(jì)算：
初始化，你覺得三個(gè)老總一樣帥粱哼，一樣有錢季二，三個(gè)美女一樣漂亮，每個(gè)人都可能跟每個(gè)人有關(guān)系揭措。所以胯舷，每個(gè)老總跟每個(gè)女職員“有問題”的概率都是1/3;
這樣，（E step）
1） A總跟小甲出去過了 1/2 * 1/3 = 1/6 次绊含，跟小乙也出去了1/6次桑嘶；（所謂的fractional count）
2）劉總跟小甲，小章也都出去了1/6次
3）劉總跟小乙躬充，小章又出去了1/6次
4）C總跟小乙出去了1/3次
總計(jì)逃顶，A總跟小甲出去了1/6次讨便，跟小乙也出去了1/6次 ; 劉總跟小甲，小乙出去了1/6次口蝠，跟小章出去了1/3次器钟；C總跟小章出去了1/3次；
你開始跟新你的八卦了(M step),
A總跟小甲妙蔗，小乙有問題的概率都是1/6 / (1/6 + 1/6) = 1/2傲霸；
劉總跟小甲，小乙有問題的概率是1/6 / (1/6+1/6+1/6+1/6) = 1/4; 跟小章有問題的概率是(1/6+1/6)/(1/6 * 4) = 1/2;
C總跟小乙有問題的概率是 1眉反。
然后昙啄，你有開始根據(jù)最新的概率計(jì)算了；（E-step）
1）A總跟小甲出去了 1/2 * 1/2 = 1/4 次寸五，跟小乙也出去 1/4 次梳凛；
2）劉總跟小甲出去了1/2 * 1/4 = 1/12 次， 跟小章出去了 1/2 * 1/2 = 1/4 次梳杏；
3）劉總跟小乙出去了1/2 * 1/4 = 1/12 次韧拒， 跟小章又出去了 1/2 * 1/2 = 1/4 次；
4）C總跟小乙出去了1次十性；
重新反思你的八卦（M-step）:
A總跟小甲叛溢，小乙有問題的概率都是1/4/ (1/4 + 1/4) = 1/2劲适；
B總跟小甲楷掉，小乙是 1/12 / (1/12 + 1/4 + 1/4 + 1/12) = 1/8 ; 跟小章是 3/4 ;
C總跟小乙的概率是1。
你繼續(xù)計(jì)算霞势，反思烹植，總之，最后愕贡，你得到了真相草雕！
EM算法過程
通過上面的計(jì)算我們可以得知，EM算法實(shí)際上是一個(gè)不停迭代計(jì)算的過程固以， 根據(jù)我們事先估計(jì)的先驗(yàn)概率A促绵，得出一個(gè)結(jié)果B，再根據(jù)結(jié)果B嘴纺，再計(jì)算得到結(jié)果 A败晴，然后反復(fù)直到這個(gè)過程收斂。
可以想象飯店的后方大廚栽渴，炒了兩盤一樣的菜尖坤，現(xiàn)在，菜炒好后從鍋中倒入盤闲擦， 不可能一下子就分配均勻慢味，所以先往兩盤中倒入场梆，然后發(fā)現(xiàn)B盤菜少了，就從A中 勻出一些纯路，A少了或油，從B勻.......不停迭代。
例如驰唬，小時(shí)候顶岸，老媽給一大袋糖果給你，叫你和你姐姐等分叫编，然后你懶得去點(diǎn)糖果的個(gè)數(shù)辖佣，所以你也就不知道每個(gè)人到底該分多少個(gè)。咱們一般怎么做呢搓逾？先把一袋糖果目測(cè)的分為兩袋卷谈，然后把兩袋糖果拿在左右手，看哪個(gè)重霞篡，如果右手重世蔗，那很明顯右手這代糖果多了，然后你再在右手這袋糖果中抓一把放到左手這袋朗兵，然后再感受下哪個(gè)重污淋，然后再?gòu)闹氐哪谴ヒ恍“逊胚M(jìn)輕的那一袋，繼續(xù)下去矛市，直到你感覺兩袋糖果差不多相等了為止芙沥。呵呵诲祸，然后為了體現(xiàn)公平浊吏，你還讓你姐姐先選了。
初始化分布參數(shù) 重復(fù)下列兩個(gè)操作直到收斂：
E步驟：估計(jì)隱藏變量的概率分布期望函數(shù)救氯；
M步驟：根據(jù)期望函數(shù)重新估計(jì)分布參數(shù)找田。
參數(shù)求解過程
假設(shè)我們有一個(gè)樣本集{x(1),…,x(m)}，包含m個(gè)獨(dú)立的樣本着憨。但每個(gè)樣本i對(duì)應(yīng)的類別z(i)是未知的（相當(dāng)于聚類）墩衙，也即隱含變量。故我們需要估計(jì)概率模型p(x,z)的參數(shù)θ甲抖，但是由于里面包含隱含變量z漆改，所以很難用最大似然求解，但如果z知道了准谚，那我們就很容易求解了挫剑。
對(duì)于參數(shù)估計(jì)，我們本質(zhì)上還是想獲得一個(gè)使似然函數(shù)最大化的那個(gè)參數(shù)θ柱衔，現(xiàn)在與最大似然不同的只是似然函數(shù)式中多了一個(gè)未知的變量z樊破，見下式（1）愉棱。也就是說我們的目標(biāo)是找到適合的θ和z讓L(θ)最大。那我們也許會(huì)想哲戚，你就是多了一個(gè)未知的變量而已啊奔滑，我也可以分別對(duì)未知的θ和z分別求偏導(dǎo)，再令其等于0顺少，求解出來不也一樣嗎朋其？
但是可以看到里面有“和的對(duì)數(shù)”，求導(dǎo)后形式會(huì)非常復(fù)雜（自己可以想象下log(f1(x)+ f2(x)+ f3(x)+…)復(fù)合函數(shù)的求導(dǎo)）祈纯，所以很難求解得到未知參數(shù)z和θ令宿。那OK，我們可否對(duì)（1）式做一些改變呢腕窥？我們看（2）式粒没，（2）式只是分子分母同乘以一個(gè)相等的函數(shù)，還是有“和的對(duì)數(shù)”啊簇爆，還是求解不了癞松，那為什么要這么做呢？咱們先不管入蛆，看（3）式响蓉，發(fā)現(xiàn)（3）式變成了“對(duì)數(shù)的和”，那這樣求導(dǎo)就容易了哨毁。我們注意點(diǎn)枫甲，還發(fā)現(xiàn)等號(hào)變成了不等號(hào)，為什么能這么變呢扼褪？這就是Jensen不等式的大顯神威的地方想幻。
步驟三用到了Jenson不等式，不等式的原理如下：
設(shè)f是定義域?yàn)閷?shí)數(shù)的函數(shù)话浇，如果對(duì)于所有的實(shí)數(shù)x脏毯。如果對(duì)于所有的實(shí)數(shù)x，f(x)的二次導(dǎo)數(shù)大于等于0幔崖，那么f是凸函數(shù)食店。當(dāng)x是向量時(shí)，如果其hessian矩陣H是半正定的赏寇，那么f是凸函數(shù)吉嫩。如果只大于0，不等于0嗅定，那么稱f是嚴(yán)格凸函數(shù)自娩。
Jensen不等式表述如下：
如果f是凸函數(shù)，X是隨機(jī)變量露戒，那么：E[f(X)]>=f(E[X])
特別地椒功，如果f是嚴(yán)格凸函數(shù)捶箱，當(dāng)且僅當(dāng)X是常量時(shí)，上式取等號(hào)动漾。
Jensen不等式如下圖所示：
OK丁屎，到這里，現(xiàn)在式（3）就容易地求導(dǎo)了旱眯，但是式（2）和式（3）是不等號(hào)啊晨川，式（2）的最大值不是式（3）的最大值啊，而我們想得到式（2）的最大值删豺，那怎么辦呢共虑？
上面的式（2）和式（3）不等式可以寫成：似然函數(shù)L(θ)>=J(z,Q)，那么我們可以通過不斷的最大化這個(gè)下界J呀页，來使得L(θ)不斷提高妈拌，最終達(dá)到它的最大值。
我們固定θ蓬蝶，調(diào)整Q(z)使下界J(z,Q)上升至與L(θ)在此點(diǎn)θ處相等尘分，然后固定Q(z)，調(diào)整θ使下界J(z,Q)達(dá)到最大值（θt到θt+1）丸氛，然后再固定θ培愁，調(diào)整Q(z)……直到收斂到似然函數(shù)L(θ)的最大值處的θ*。
在Jensen不等式中說到缓窜，當(dāng)自變量X是常數(shù)的時(shí)候定续，等式成立。而在這里禾锤，即：
再推導(dǎo)下私股，由于
（因?yàn)镼是隨機(jī)變量z(i)的概率密度函數(shù)），則可以得到：分子的和等于c（分子分母都對(duì)所有z(i)求和：多個(gè)等式分子分母相加不變时肿，這個(gè)認(rèn)為每個(gè)樣例的兩個(gè)概率比值都是c）庇茫，則：
至此港粱，我們推出了在固定參數(shù)θ后螃成，使下界拉升的Q(z)的計(jì)算公式就是后驗(yàn)概率，解決了Q(z)如何選擇的問題查坪。這一步就是E步寸宏，建立L(θ)的下界。接下來的M步偿曙，就是在給定Q(z)后氮凝，調(diào)整θ，去極大化L(θ)的下界J（在固定Q(z)后望忆，下界還可以調(diào)整的更大）罩阵。那么一般的EM算法的步驟如下：
EM算法（Expectation-maximization）：
期望最大算法是一種從不完全數(shù)據(jù)或有數(shù)據(jù)丟失的數(shù)據(jù)集（存在隱含變量）中求解概率模型參數(shù)的最大似然估計(jì)方法竿秆。
EM的算法流程：
初始化分布參數(shù)θ；
重復(fù)以下步驟直到收斂：
E步驟：根據(jù)參數(shù)初始值或上一次迭代的模型參數(shù)來計(jì)算出隱性變量的后驗(yàn)概率稿壁，其實(shí)就是隱性變量的期望幽钢。作為隱藏變量的現(xiàn)估計(jì)值：
M步驟：將似然函數(shù)最大化以獲得新的參數(shù)值：
這個(gè)不斷的迭代，就可以得到使似然函數(shù)L(θ)最大化的參數(shù)θ了傅是。
示例
簡(jiǎn)單給出一個(gè)EM算法的示例匪燕，代碼如下：
#-*- conding:utf-8 -*-
# --encoding:utf-8 --
"""
實(shí)現(xiàn)GMM高斯混合聚類
根據(jù)EM算法流程實(shí)現(xiàn)這個(gè)流程
http://scipy.github.io/devdocs/generated/scipy.stats.multivariate_normal.html#scipy.stats.multivariate_normal
"""
import numpy as np
from scipy.stats import multivariate_normal
def train(x, max_iter=100):
"""
進(jìn)行GMM模型訓(xùn)練，并返回對(duì)應(yīng)的μ和σ的值(假定x數(shù)據(jù)中的簇類別數(shù)目為2)
:param x: 輸入的特征矩陣x
:param max_iter: 最大的迭代次數(shù)
:return: 返回一個(gè)五元組(pi, μ1喧笔， μ2帽驯，σ1，σ2)
"""
# 1. 獲取樣本的數(shù)量m以及特征維度n
m, n = np.shape(x)
# 2. 初始化相關(guān)變量
# 以每一列中的最小值作為mu1书闸，mu1中的元素?cái)?shù)目就是列的數(shù)目(n)個(gè)
mu1 = x.min(axis=0)
mu2 = x.max(axis=0)
sigma1 = np.identity(n)
sigma2 = np.identity(n)
pi = 0.5
# 3. 實(shí)現(xiàn)EM算法
for i in range(max_iter):
# a. 初始化多元高斯分布（初始化兩個(gè)多元高斯混合概率密度函數(shù)）
norm1 = multivariate_normal(mu1, sigma1)
norm2 = multivariate_normal(mu2, sigma2)
# E step
# 計(jì)算所有樣本數(shù)據(jù)在norm1和norm2中的概率
tau1 = pi * norm1.pdf(x)
tau2 = (1 - pi) * norm2.pdf(x)
# 概率做一個(gè)歸一化操作
w = tau1 / (tau1 + tau2)
# M step
mu1 = np.dot(w, x) / np.sum(w)
mu2 = np.dot(1 - w, x) / np.sum(1 - w)
sigma1 = np.dot(w * (x - mu1).T, (x - mu1)) / np.sum(w)
sigma2 = np.dot((1 - w) * (x - mu2).T, (x - mu2)) / np.sum(1 - w)
pi = np.sum(w) / m
# 返回最終解
return (pi, mu1, mu2, sigma1, sigma2)
if __name__ == '__main__':
np.random.seed(28)
# 產(chǎn)生一個(gè)服從多元高斯分布的數(shù)據(jù)（標(biāo)準(zhǔn)正態(tài)分布的多元高斯數(shù)據(jù)）
mean1 = (0, 0, 0) # x1\x2\x3的數(shù)據(jù)分布都是服從正態(tài)分布的尼变，同時(shí)均值均為0
cov1 = np.diag((1, 1, 1))
data1 = np.random.multivariate_normal(mean=mean1, cov=cov1, size=500)
# 產(chǎn)生一個(gè)數(shù)據(jù)分布不均衡
mean2 = (2, 2, 3)
cov2 = np.array([[1, 1, 3], [1, 2, 1], [0, 0, 1]])
data2 = np.random.multivariate_normal(mean=mean2, cov=cov2, size=200)
# 合并兩個(gè)數(shù)據(jù)
data = np.vstack((data1, data2))
pi, mu1, mu2, sigma1, sigma2 = train(data, 100)
print("第一個(gè)類別的相關(guān)參數(shù):")
print(mu1)
print(sigma1)
print("第二個(gè)類別的相關(guān)參數(shù):")
print(mu2)
print(sigma2)
print("預(yù)測(cè)樣本屬于那個(gè)類別(概率越大就是那個(gè)類別)：")
norm1 = multivariate_normal(mu1, sigma1)
norm2 = multivariate_normal(mu2, sigma2)
x = np.array([0, 1, 0])
print(pi * norm1.pdf(x)) # 屬于類別1的概率為:0.0275 => 0.989
print((1 - pi) * norm2.pdf(x))# 屬于類別1的概率為:0.0003 => 0.011
運(yùn)行結(jié)果如下：
第一個(gè)類別的相關(guān)參數(shù):
[-0.03813579 -0.0265211 0.06435217]
[[ 0.89477712 -0.01408823 -0.0630486 ]
[-0.01408823 1.04926061 0.05625028]
[-0.0630486 0.05625028 1.01136971]]
第二個(gè)類別的相關(guān)參數(shù):
[1.89859865 1.9345104 2.77821259]
[[0.84837261 1.11138647 1.12554675]
[1.11138647 2.19189624 1.37848616]
[1.12554675 1.37848616 3.43711088]]
預(yù)測(cè)樣本屬于那個(gè)類別(概率越大就是那個(gè)類別)：
0.02758258224417439
0.00039921371592994276
Process finished with exit code 0
