天天看點
傳回

高斯混合模型EM算法

k-means算法對于資料點和clusters之間的關系，是all-or-nothing的關系，這是一個hard decision，往往會導緻的局部最小值，這不是理想的求解。一種常見的做法，是學習這個協方差矩陣，而不是固定它們為機關矩陣。

ＧＭＭ模型及算法流程

ＧＭＭ的全稱是Gaussian Mixture Model，即高斯混合模型。

假設我們有一個訓練集 {x1,...,xm} ,在非監督學習中，這些資料都是沒有标簽的。我們希望可以對資料進行模組化，通過定義一個聯合分布 p(x(i),z(i))=p(x(i)|z(i))p(z(i)) ,這裡 z(i)∼Multinomial(ϕ) ,其中 ϕj≥0,∑kj=1ϕj=1 , 這個參數 ϕj 代表着 p(z(i)=j) 的機率。

而 xi|zi=j∼N(μj,Σj) ,ｋ代表着标簽 z(i) 可以有的取值的數目。是以我們的模型實際上是當 z(i) 從 {1,...,k} 随機取值時，從ｋ個高斯分量中選一個，再根據相應的 z(i) 生成對應的 x(i) 。

z(i) 是一個隐變量，就是說它是隐藏的，不可見的。

我們的似然函數可以這樣寫：

L(ϕ,μ,Σ)=∑i=1mlog p(x(i);ϕ,μ,Σ)=∑i=1mlog ∑z(i)=1kp(x(i)|z(i);μ,Σ)p(z(i);ϕ)

如果直接通過令偏導數為０，對其進行求解是無法得到closed form的。

這個隐變量告訴了我們，資料 x(i) 是從ｋ個高斯分量中的哪一個來的。如果我們一開始就知道是哪個，那麼最大似然問題是很簡單的。具體而言，我們就可以将似然函數展開：

L(ϕ,μ,Σ)=∑i=1mlog p(x(i)|z(i);μ,Σ)+log p(z(i);ϕ)

随後分别對 ϕ , μ , Σ 求導，可以得到：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪ϕj=1m∑mi=1I{z(i)=j}μj=∑mi=1I{z(i)=j}x(i)∑mi=1I{z(i)=j}Σj=∑mi=1I{z(i)=j}x((i)−μj)(x(i)−μj)T∑mi=1I{z(i)=j}

當我們知道類别資訊，事實上這就是一個高斯判别分析的參數估計問題。但是我們知道，是以我們可以分兩個步驟:

一個步驟猜測這個 zi 的值，也就是Ｅ-step;

一個步驟利用已知的類别資訊來估計高斯的參數 μ,Σ ，以及各個高斯分量的權重 ϕj 。

那麼這個标簽值要如何猜測呢? 我們可以在固定參數 ϕ,μ,Σ 的情況下，在給定資料 x(i) 時 z(i) 後驗機率的來判定，再結合貝葉斯定律：

p(z(i)=j|x(i);ϕ,μ,Σ)=p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)∑kl=1p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)

EM算法在k-means 聚類中也有用到，但是用的是hard decision， 而利用GMM，我們使用的其實是soft的 assignments w(i)j . 同樣它也可能會收斂到局部極小值，是以需要嘗試不同的初始參數。

用GMM來進行聚類，K個高斯component實際上就對應着K個cluster，我們根據資料來推算機率密度density estimation。

利用Mixtures of Gaussians模型來實作聚類算法的流程總結如下：

輸入：

{x(i)}mi=1 (資料)， k(聚類的數目)

初始化：

∀μi←1k

Σi←var({x(i)})

Repeat until convergence :

E-step:

p(z(i)=j|x(i);ϕ,μ,Σ)=p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)∑kl=1p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)

w(i)j:=p(z(i)=j|x(i);ϕ,μ,Σ)

M-step:

ϕj:=1m∑i=1mw(i)j

,

μj:=∑mi=1w(i)jx(i)∑mi=1w(i)j

Σj:=∑mi=1w(i)j(x(i)−μj)(x(i)−μj)T∑mi=1w(i)j

輸出：

{z(i)←argmaxj∈{1,...,k}w(i)j}mi=1

算法代碼demo及關鍵步驟注解

大費苦心寫了編輯了這麼多公式，還沒有進入正題。本文一大重點是要研讀scikit包中gaussian mixture 部分的代碼，已進一步熟悉算法和這個工具包。

官網給出的密度估計的源碼可見：

http://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_pdf.html#example-mixture-plot-gmm-pdf-py

算法就是随機産生兩個不同中心的高斯分布的資料點，再用mixture.GMM進行拟合。

關鍵部分代碼拆分如下。

1. 準備資料

産生以20,20 為中心的滿足标準正态分布的300個樣本點。

shifted_gaussian = np.random.randn(n_samples, ) + np.array([, ])

産生原點處的stretch後高斯分布的300個樣本點。通過乘以一個矩陣來完成。

C = np.array([[0., -0.7], [3.5, .7]])
stretched_gaussian = np.dot(np.random.randn(n_samples, ), C)

疊加兩部分資料：

訓練資料是一個600*2的矩陣。

2. 初始化高斯混合模型

clf = mixture.GMM(n_components=2, covariance_type=’full’)

第二個參數這樣設定，實際得到的協方差的參數 covars_是

(n_components, n_features, n_features)

這樣允許了不同的高斯分量有不同的協方差矩陣。

初始化就做了這些事情：

def __init__(self, n_components=, covariance_type='diag',
                 random_state=None, thresh=, min_covar=,
                 n_iter=, n_init=, params='wmc', init_params='wmc'):
        self.n_components = n_components
        self.covariance_type = covariance_type
        self.thresh = thresh
        self.min_covar = min_covar
        self.random_state = random_state
        self.n_iter = n_iter
        self.n_init = n_init
        self.params = params
        self.init_params = init_params

預設疊代的次數為100次

self.weights_ = np.ones(self.n_components) / self.n_components

每個分量一開始指派的權重是相同的，這裡的例子中有兩個分量，是以得到的是一個行向量[0.5,0.5].

密度估計

随後要做的事情，是生成模型。

調用方式為：

clf.fit(X_train)

裡面的關鍵步驟，就看Expection Step和Maximization Step。

Expection Step:

curr_log_likelihood, responsibilities = self.score_samples(X)
      log_likelihood.append(curr_log_likelihood.sum())

      # Check for convergence.
      if i >  and abs(log_likelihood[-] - log_likelihood[-]) < \
              self.thresh:
          self.converged_ = True
          break

第一個函數就是計算後驗機率的函數，

其中作者對

p(z(i)=j|x(i);ϕ,μ,Σ)=p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)∑kl=1p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)

表示成exp(log(.))的形式，于是把除法就轉換成了減法。

下面這一步是計算分子。相當于計算 log(p(x(i)|z(i)=j;μ,Σ))+log(p(z(i)=j;ϕ))

lpr = (log_multivariate_normal_density(X, self.means_, self.covars_,
                                               self.covariance_type)
               + np.log(self.weights_))

上面得到的矩陣為600*2，一列對應着一個高斯分量。下面這一步是計算分母，先用exp得到原值，再求和，随後再取log相當于計算 log∑kl=1p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)

萬分精簡的隻用了一個函數就完成了：

計算結果減法完成：

responsibilities = np.exp(lpr - logprob[:, np.newaxis])

函數的傳回值中，

responsibilities : array_like, shape (n_samples, n_components)

Posterior probabilities of each mixture component for each

observation

而curr_log_likelihood則是分母的值，作者通過記錄這個值所有600個樣本的和在每一次疊代中的變化來與門檻值self.thresh比較來判斷是否達到了收斂。

Maximization step

def _do_mstep(self, X, responsibilities, params, min_covar=):
        """ Perform the Mstep of the EM algorithm and return the class weihgts.
        """
        weights = responsibilities.sum(axis=) 
        weighted_X_sum = np.dot(responsibilities.T, X)
        inverse_weights =  / (weights[:, np.newaxis] +  * EPS)

        if 'w' in params:
            self.weights_ = (weights / (weights.sum() +  * EPS) + EPS)
        if 'm' in params:
            self.means_ = weighted_X_sum * inverse_weights
        if 'c' in params:
            covar_mstep_func = _covar_mstep_funcs[self.covariance_type]
            self.covars_ = covar_mstep_func(
                self, X, responsibilities, weighted_X_sum, inverse_weights,
                min_covar)
        return weights

前面已經說了：

ϕj 代表着 p(z(i)=j) 的機率，之類對應着weights_, 每個分量的權重。

計算的方法是

ϕj:=1m∑i=1mw(i)j

,

也等于600個樣本的分量1的後驗機率和/（600個樣本的分量1的後驗機率和+600個樣本的分量2的後驗機率和）。

weights = responsibilities.sum(axis=)
    if 'w' in params:
        self.weights_ = (weights / (weights.sum() +  * EPS) + EPS)

注意每一行兩個分量的權重求和為1,600個全部加起來為600，正好等于左右兩列歸總值的和。是以這麼計算也是對的。

更新均值：

μj:=∑mi=1w(i)jx(i)∑mi=1w(i)j

%得到*的矩陣，（,）代表分量第一維，（,）代表分量第二維
   weighted_X_sum = np.dot(responsibilities.T, X) 
   weights =  / (weights[:, np.newaxis] +  * EPS)
   if 'm' in params:
        self.means_ = weighted_X_sum * inverse_weights % 同樣是*的矩陣

更新協方差矩陣：

Σj:=∑mi=1w(i)j(x(i)−μj)(x(i)−μj)T∑mi=1w(i)j

if 'c' in params:
       covar_mstep_func = _covar_mstep_funcs[self.covariance_type]
       self.covars_ = covar_mstep_func(
           self, X, responsibilities, weighted_X_sum, inverse_weights,
           min_covar)

這裡實際調用的函數為：

def _covar_mstep_full(gmm, X, responsibilities, weighted_X_sum, norm,
                      min_covar):
    """Performing the covariance M step for full cases"""
    # Eq. 12 from K. Murphy, "Fitting a Conditional Linear Gaussian
    # Distribution"
    n_features = X.shape[]
    cv = np.empty((gmm.n_components, n_features, n_features))
    for c in range(gmm.n_components):
        post = responsibilities[:, c]
        # Underflow Errors in doing post * X.T are  not important
        np.seterr(under='ignore')
        avg_cv = np.dot(post * X.T, X) / (post.sum() +  * EPS)
        mu = gmm.means_[c][np.newaxis]
        cv[c] = (avg_cv - np.dot(mu.T, mu) + min_covar * np.eye(n_features))
    return cv

具體可見參考文獻2的公式12. 其實我看參考文獻2脫離上下文不怎麼容易懂，基本上應該說這是協方差矩陣的另一種寫法。https://en.wikipedia.org/wiki/Covariance_matrix

Σ=E[(X−EX)(X−EX)T]=E(XXT)−μμT

是以上面的公式就進一步寫成：

Σj:=∑mi=1w(i)jx(i)x(i)T∑mi=1w(i)j−μjμTj

下一次再來梳理關鍵步驟的推導吧，最後的結果是好看的。

高斯混合模型EM算法

參考文獻：

１. andrew ng cs229 handouts

http://cs229.stanford.edu/notes/cs229-notes7b.pdf

2. K. Murphy, “Fitting a Conditional Linear Gaussian Distribution”

http://www.cs.ubc.ca/~murphyk/Papers/learncg.pdf

機器學習算法
上一篇: EM算法--二維高斯混合模型(GMMs)
下一篇: EM算法過程推導

繼續閱讀