EM算法 - 簡書

緣起

最近室友都找實(shí)習(xí)。面試時矾削，常被問到的是壤玫，一些基礎(chǔ)的機(jī)器學(xué)習(xí)知識和最新的深度學(xué)習(xí)模型，比如邏輯回歸哼凯、正則項(xiàng)欲间、f1的計(jì)算、類別不平衡的處理断部、LSTM猎贴、Attention、Transformer蝴光、word2vec她渴、glove、ELMo蔑祟、BERT等等趁耗。因此，我打算系統(tǒng)地疆虚、扎實(shí)地將這些知識學(xué)習(xí)一遍苛败，同時將學(xué)習(xí)的筆記寫成博客右冻，以供自己以后查看。如果不幸?guī)椭藙e人著拭，還可以收獲一些虛榮:-)。對吧~

萬事開頭難牍帚。所以儡遮，先從學(xué)過的EM算法開始吧。

從貝葉斯公式到極大似然估計(jì)

貝葉斯估計(jì)

已知樣本的特征向量 $x$ 暗赶，為了求樣本的類別 $w$ 鄙币，我們可以使用經(jīng)典的貝葉斯公式，求解已知特征向量 $x$ 時不同類別的概率 $P(w_i|x)$ 蹂随，
$P(w_i|x)=\frac{P(x|w_i)P(w_i)}{P(x)}.$
然后十嘿，選擇使得 $P(w_i|x)$ 中最大的 $w_i$ 作為樣本類別的估計(jì)值，即有
$\begin{align} \hat{w}&=\arg \max_{i} P(w_i|x)\\ &=\arg \max_{i} P(x|w_i)P(w_i). \end{align}$
其中岳锁， $P(w_i)$ 被稱為先驗(yàn)概率绩衷，可通過統(tǒng)計(jì)標(biāo)記樣本中不同類別樣本的比例得到，即
$P(w_i)=\frac{\# w_i}{N}.$
$P(x|w_i)$ 被稱為類條件概率激率，當(dāng)樣本的特征數(shù)為1且是為離散值時咳燕，
$P(x|w_i)=\frac{\#(x,w_i)}{\# w_i}.$
其中， $\#(x,w_i)$ 表示特征為 $x$ 且類別為 $w_i$ 的樣本數(shù)目乒躺。樣本特征數(shù)不為1的情況招盲，這里暫且不討論。當(dāng)特征 $x$ 的取值是連續(xù)值的時候嘉冒，問題就復(fù)雜了起來曹货。我們可以使用極大似然估計(jì)來解決這個問題。

極大似然估計(jì)

在極大似然估計(jì)中讳推，我們首先假設(shè) $P(x|w_i)$ 服從一個概率分布顶籽，常用的是正態(tài)分布，即
$P(x|w_i;\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2\right].$
記參數(shù) $\theta=(\mu,\sigma )$ 娜遵，當(dāng)我們得到一組 $\theta$ 時蜕衡，類條件概率 $P(x|w_i)$ 隨之確定。

為了求解 $\theta$ 设拟，令所有類別為 $w_i$ 的樣本集合為 $D$ 慨仿，假設(shè) $D$ 中樣本 $x_1,\cdots,x_n$ 之間相互獨(dú)立，那么有
$\begin{align} P(D|\theta)&=P(x_1,x_2,\cdots,x_n|\theta)=\prod_{i=1}^n p(x_i|\theta) \end{align}$
$P(D|\theta)$ 被稱為似然函數(shù)纳胧，它描述了參數(shù) $\theta$ 的合理程度镰吆。一組合適的參數(shù) $\theta$ 應(yīng)能使得 $P(D|\theta)$ 盡量大。所以跑慕，求解 $\theta$ 最直觀的想法是找到一組 $\theta$ 使得 $P(D|\theta)$ 最大万皿。為了計(jì)算的方便摧找，我們定義對數(shù)似然函數(shù)
$l(\theta)=\ln P(D|\theta)=\sum_{i=1}^n \ln P(x_i|\theta).$
于是， $\theta$ 的估計(jì)值為
$\hat{\theta}=\arg \max_{\theta} l(\theta).$
由于極值點(diǎn)往往在導(dǎo)數(shù)等于0的地方取得牢硅，我們常常通過求導(dǎo)的方法求解上式蹬耘，即
$\frac{\partial l(\hat{\theta})}{\partial \mu} = 0, \frac{\partial l(\hat{\theta})}{\partial \sigma}=0.$
解得
$\begin{align} \hat{\mu}&=\frac{1}{n} \sum_{i=1}^n x_i,\\ \hat{\sigma}&=\frac{1}{n}\sum_{i=1}^n (x_i-\hat{\mu})^2. \end{align}$

引入EM算法

在上一節(jié)，我們看到極大似然估計(jì)通過最大化似然函數(shù)的方法减余，找到模型的最佳參數(shù)综苔。在這個過程中，我們需要對對數(shù)似然函數(shù)求偏導(dǎo)位岔，令導(dǎo)數(shù)等于0如筛，進(jìn)而得到參數(shù)的估計(jì)值。然而抒抬，在某些復(fù)雜的模型中杨刨，由于存在“隱變量”，我們難以通過對對數(shù)似然函數(shù)求偏導(dǎo)的方法得到參數(shù)的估計(jì)值擦剑。

我們考慮如下的三硬幣問題

假設(shè)有3枚硬幣妖胀，分別記作A、B抓于、C做粤，這些硬幣正面朝上的概率分別為 $\pi,p,q$ 。現(xiàn)進(jìn)行如下的拋擲實(shí)驗(yàn)：首先擲硬幣A捉撮，根據(jù)其結(jié)果選擇硬幣B和硬幣C怕品，如果是正面選擇硬幣B，如果是反面選擇硬幣C巾遭；然后擲選出的硬幣肉康，記錄擲硬幣的結(jié)果，正面記作1灼舍，反面記作0吼和；獨(dú)立地重復(fù)n次實(shí)驗(yàn)（這里，n=10）骑素，觀測結(jié)果如下：
$1,1,0,1,0,0,1,0,1,1.$
假設(shè)只能觀測到擲硬幣的結(jié)果炫乓，不能觀測擲硬幣的過程。問如何估計(jì)三硬幣模型的參數(shù) $\pi,p,q$ 献丑。

記 $\theta=(\pi,p,q)$ 末捣。首先，讓我們嘗試使用極大似然估計(jì)的方法创橄，寫出對數(shù)似然函數(shù)箩做，然后求解 $\theta$ 。不妨令每次實(shí)驗(yàn)中硬幣A的拋擲結(jié)果為 $z$ 妥畏，硬幣B和C的拋擲結(jié)果為 $y$ 邦邦。于是
$P(y|\theta)=\sum_z P(y,z|\theta)=\sum_z P(z|\theta)P(y|z,\theta).$
顯然安吁， $z$ 的取值有0、1兩種燃辖，并且 $P(z=0|\theta)=1-\pi,P(z=1|\theta)=\pi$ 鬼店。而當(dāng) $z=0$ 時，意味著選擇了硬幣C黔龟，此時 $P(y|z,\theta)=q^y(1-q)^{1-y}$ 薪韩；同理，當(dāng) $z=1$ 時捌锭，有 $P(y|z,\theta)=p^y(1-p)^{1-y}$ 。于是罗捎，上式可展開為
$P(y|\theta)=\pi p^y(1-p)^{1-y}+(1-\pi)q^y(1-q)^{1-y}.$
對于三硬幣模型多次拋擲的結(jié)果 $Y=\left[y_1,y_2,\cdots,y_n\right]$ 观谦，有
$P(Y|\theta)=\prod_{i=1}^n P(y_i|\theta).$
相應(yīng)的對數(shù)似然函數(shù)為
$l(\theta)=\ln P(Y|\theta)=\sum_{i=1}^n \ln \left[ \pi p^y_i(1-p)^{1-y_i}+(1-\pi)q^{y_i}(1-q)^{1-y_i} \right].$
事實(shí)證明對上式求偏導(dǎo)，令導(dǎo)數(shù)等于0桨菜，是無法得到解析解的豁状。也就是極大似然估計(jì)無法解決上述問題。

使用EM算法求解

一方面倒得，假設(shè)我們知道模型參數(shù) $\pi,p,q$ 泻红，我們就可以預(yù)測某次實(shí)驗(yàn)中的拋擲結(jié)果來源于硬幣B的概率，不妨記 $\mu=P(z=1|y,\theta)$ 霞掺，有
$\mu = \frac{P(y|z=1,\theta)}{P(y|z=1,\theta)+P(y|z=0,\theta)} = \frac{\pi p^y(1-p)^{1-y}}{\pi p^y(1-p)^{1-y} + (1-\pi) q^y(1-q)^{1-y}}.$
另一方面谊路，假設(shè)我們已知每次實(shí)驗(yàn)中拋擲結(jié)果來源于硬幣B的概率，我們就可以估計(jì)模型參數(shù) $\pi,p,q$ 菩彬，即
$\begin{align} \hat{\pi} &= \frac{1}{n} \sum_{i=0}^n \mu_i,\\ \hat{p} &= \frac{\sum_{i=0}^n \mu_iy_i}{\sum_{i=0}^n \mu_i},\\ \hat{q} &= \frac{\sum_{i=0}^n (1-\mu_i)y_i}{\sum_{i=0}^n (1-\mu_i）}. \end{align}$
也就是說缠劝，當(dāng)我們已知模型參數(shù)的時候，我們可以計(jì)算 $\mu$ 的值骗灶；當(dāng)我們知道 $\mu$ 的值時惨恭，我們可以估計(jì)模型參數(shù)。對于這種蛋雞問題耙旦，使用迭代法求解是最好不過了脱羡。也就是，先假設(shè)一組模型參數(shù)免都，比如 $\pi=0.5,p=0.5,q=0.5$ 锉罐，然后計(jì)算 $\mu$ 的值；接著使用 $\mu$ 的值琴昆，去估計(jì)模型參數(shù)氓鄙。重復(fù)上述過程多次之后，我們就可以得到 $\hat{\pi}=0.5,\hat{p}=0.6,\hat{q}=0.6$ 业舍。

在上述過程中抖拦，我們稱 $\mu$ 為模型的隱變量升酣，也即是無法觀測到但也不是模型參數(shù)的變量。估計(jì) $\mu$ 的過程稱為E步态罪，使用 $\mu$ 估計(jì)模型參數(shù)的過程稱為M步噩茄。下面使用python實(shí)現(xiàn)了對三硬幣模型的求解。

def E_step(pi,p,q,Y):
    Mu = []
    for y in Y:
        pyz1 = pi * p**y * (1-p)**(1-y)
        pyz0 = (1-pi) * q**y * (1-q)**(1-y)
        mu = pyz1 / (pyz1 + pyz0)
        Mu.append(mu)
    return Mu

def M_step(Mu,Y):
    n = len(Mu)
    pi = sum(Mu)/n
    p = sum([Mu[i]*Y[i] for i in range(n)]) / sum(Mu)
    q = sum([(1-Mu[i])*Y[i] for i in range(n)]) / (n-sum(Mu))
    return pi,p,q

def solve_tree_coin(Y,iter_num):
    pi = p = q = 0.5
    for i in range(iter_num):
        Mu = E_step(pi,p,q,Y)
        pi,p,q = M_step(Mu,Y)
        print(pi,p,q)


if __name__ == "__main__":
    Y = [1,1,0,1,0,0,1,0,1,1]
    solve_tree_coin(Y,5)

說句題外話复颈，三硬幣問題是沒有辦法正確估計(jì)模型的參數(shù)绩聘。如 $\pi = 0.6,p = 1,q = 0$ 和 $\pi=0.5,p=0.6,q=0.6$ 在輸出結(jié)果上是等價的。畢竟耗啦，對于這個模型觀測數(shù)據(jù)太少了凿菩，任何方法對此都無能為力。

下面我們將正式地對EM算法進(jìn)行介紹帜讲，并運(yùn)用EM算法求解混合高斯模型衅谷。

EM算法的思想

在這里不介紹EM算法的形式化定義（因?yàn)闆]有必要），只是概括一下EM算法的核心思想似将。

EM算法的全稱是Expectation Maximization Algorithm获黔，也就是期望極大算法。EM算法適用于含有隱變量的模型在验。求解過程分為E步和M步：在E步中求解隱變量的期望玷氏；在M步中使用隱變量的期望代替隱變量的值，求解模型參數(shù)腋舌。

理論上（有興趣的見李航的《統(tǒng)計(jì)學(xué)習(xí)方法）盏触，EM算法并不能保證得到全局最優(yōu)值，而且很依賴所選取初始值的好壞块饺，但EM算法簡化了很多復(fù)雜問題的求解耻陕。計(jì)算機(jī)科學(xué)中的許多方法都是基于這樣的思路。

使用EM求解GMM

首先介紹一下GMM刨沦。

GMM的全稱是Gauss Mixture Model诗宣，即混合高斯模型，也就是將多個高斯分布疊加在一起想诅。假設(shè)GMM中高斯分布的數(shù)量為M召庞，那么相應(yīng)的概率密度函數(shù)為
$p(x) = \sum_{i=1}^M a_iN(x;\mu,\Sigma).$
如下圖就是一個GMM概率密度函數(shù)的函數(shù)圖像。

概率密度函數(shù)P(x)=0.7N(-10,2)+0.3N(5,3)的函數(shù)圖像

GMM產(chǎn)生樣本的過程和三硬幣模型有些類似来破，首先按照先驗(yàn)概率 $a_i$ 選擇一個高斯分布篮灼，然后產(chǎn)生一個滿足這個高斯分布的樣本。如下圖即為GMM產(chǎn)生的2維樣本數(shù)據(jù)徘禁。

為了使用EM算法诅诱，我們首先要確定模型中的隱變量。在GMM中送朱，隱變量是樣本屬于不同高斯分布的概率 $P(y=i)$ 娘荡。于是在E步干旁，我們這樣估計(jì)隱變量
$P(y_t=i)=\frac{a_iN(x_t;\mu_i,\Sigma_i)}{\sum_{i=1}^M a_iN(x_t;\mu_i,\Sigma_i)}.$
在M步，模型參數(shù)的估計(jì)過程如下
$\begin{align} a_i & =\frac{1}{n}\sum_{t=1}^n P(y_t=i)\\ \mu_i &= \frac{\sum_{t=1}^n P(y_t=i)x_t}{\sum_{t=1}^n P(y_t=i)}\\ \Sigma_t &= \frac{\sum_{t=1}^n P(y_t=i)(x_t-\mu_i)(x_t-\mu_i)^t}{\sum_{t=1}^n P(y_t=i)} \end{align}$

介紹EM算法的大多文章都是拿GMM舉例的炮沐，但作為機(jī)器學(xué)習(xí)中的一類重要的思想争群，EM算法的應(yīng)用非常廣泛，如IBM翻譯模型大年、主題模型等等换薄。

先寫到這里吧。