從EM算法到高斯混合模型（GMM）與K-means算法（理論分析+仿真分析）

目錄：

• EM算法
• GMM模型
• K-means算法

EM算法

基本問題

  首先介紹一下EM算法要解決的基本問題。EM算法通常解決的是不同分布下的參數估計問題，屬于參數估計理論的一種。在另一個角度，EM算法也可以看作是極大似然估計的加強版。

  由于極大似然估計要求采樣資料屬于同一個分布，因為隻能估計同一個分布下的參數。

  舉個例子：在學校随機地選取100個男生，可以通過極大似然估計來得到男生的平均身高。倘若随機選取100個未知性别的學生（有男有女），這時，這100個人的平均身高既不代表男生也不代表女生。

  此時，正确地做法應該是先把女生和男生分開，然後分别用極大似然估計，最後分别得到了男生和女生的平均身高。

數學描述

  若樣本 X = { X 1 , X 2 . . . . X n } X=\{X_1,X_2....X_n \} X={X1​,X2​....Xn​}，一共有 N N N個不同的樣本，其中樣本 X i X_i Xi​屬于 K ( k = 1 , 2... K ) K(k=1,2...K) K(k=1,2...K)個不同的分布。即 X i ∼ p ( x ∣ θ k ) X_i \sim p(x|\theta_k) Xi​∼p(x∣θk​)， X i X_i Xi​是從第 k k k個分布裡面抽取出來的樣本。現在的問題是需要得到 p ( x ∣ θ k ) p(x|\theta_k) p(x∣θk​)。下面就是EM算法要解決的問題了，做法很簡單，思想确實很精巧。

1.極大似然估計

  EM算法還是從極大似然出發的，前面說過，這玩意就是極大似然的加強版本。首先寫出樣本 X = { X 1 , X 2 . . . . X n } X=\{X_1,X_2....X_n \} X={X1​,X2​....Xn​}的似然函數:

L ( θ ) = ∏ i = 1 N p ( x i ∣ θ k ) (1) L(\theta)=\prod_{i=1}^Np(x_i|\theta_k)\tag{1} L(θ)=i=1∏N​p(xi​∣θk​)(1)

  觀察 ( 1 ) (1) (1)式子可以看出似然函數中含有 θ k \theta_k θk​, k k k是未知的，是以這裡引入一個叫隐變量的東西。首先隐變量有幾個特點（個人了解）：

• 隐變量是一個随機變量
• 每個樣本 X i X_i Xi​都對應一個隐變量值 Z j Z_j Zj​，但是不意味隐變量可以取 N N N個值，同一分布内（同一類）的樣本對應隐變量的值是一樣的。
• 隐變量描述的是 X i X_i Xi​歸屬于哪個分布的現象。

  主要解釋一下後面兩點，想想最開始提到的100個學生，其實每個學生我們不僅僅能測量他的身高，我們還應該知道他的性别。是以，完整地描述一個人的身高應該是 ( X i , Z j ) (X_i,Z_j) (Xi​,Zj​)，前者表示升高，後者表示性别。

  利用隐變量，對 ( 1 ) (1) (1)式進行改寫：

L ( θ ) = ∏ i = 1 N ∑ k = 1 K p ( x i ∣ θ k ) (2) L(\theta)=\prod_{i=1}^N \sum_{k=1}^Kp(x_i|\theta_k)\tag{2} L(θ)=i=1∏N​k=1∑K​p(xi​∣θk​)(2)

  這裡的解釋是既然不知道每個樣本 X i X_i Xi​所對應的分布 p ( x ∣ θ k ) p(x|\theta_k) p(x∣θk​)，那麼幹脆下面主要處理

X i ∼ ∑ k = 1 K p ( x i ∣ θ k ) X_i\sim\sum_{k=1}^Kp(x_i|\theta_k) Xi​∼k=1∑K​p(xi​∣θk​)

  這個其實是代表這每個分布都對 X i X_i Xi​有貢獻，隻是大小不一樣，那麼就直接求和取最大就好了。

  對 ( 2 ) (2) (2)式取對數，有：

l o g [ L ( θ ) ] = ∑ i = 1 N l o g ∑ k = 1 K p ( x i ∣ θ k ) (3) log[L(\theta)]=\sum_{i=1}^N log\sum_{k=1}^Kp(x_i|\theta_k)\tag{3} log[L(θ)]=i=1∑N​logk=1∑K​p(xi​∣θk​)(3)

  按照極大似然的做法，接下來要對 ( 3 ) (3) (3)式求導了，由于對數内求和再求導，将會無比複雜，是以不能求導。這裡需要用到一個不等式。

2.琴聲不等式

  對于一個凹函數有以下不等式成立：

f ( E ( x ) ) ≥ E ( f ( x ) ) f(E(x))\ge E(f(x)) f(E(x))≥E(f(x))

  如果 x x x是離散的随機變量，那麼有：

f ( ∑ x p ( x ) ) ≥ ∑ f ( x ) p ( x ) f(\sum xp(x))\ge \sum f(x)p(x) f(∑xp(x))≥∑f(x)p(x)

  前面提到，每個樣本 Z i Z_i Zi​對于一個隐變量的取值 Z j Z_j Zj​,設隐變量 Z j Z_j Zj​對應的機率密度值為 Q ( Z j ) Q(Z_j) Q(Zj​)， ( 3 ) (3) (3)式可寫成如下形式：

l o g [ L ( θ ) ] = ∑ i = 1 N l o g ∑ k = 1 K Q ( Z j ) p ( x i ∣ θ k ) Q ( Z j ) (4) log[L(\theta)]=\sum_{i=1}^N log\sum_{k=1}^KQ(Z_j)\frac{p(x_i|\theta_k)}{Q(Z_j)}\tag{4} log[L(θ)]=i=1∑N​logk=1∑K​Q(Zj​)Q(Zj​)p(xi​∣θk​)​(4)

對 ( 4 ) (4) (4)應用琴聲不等式。

f − − − > l o g f--->log f−−−>log

x − − − > p ( x i ∣ θ k ) Q ( Z j ) x--->\frac{p(x_i|\theta_k)}{Q(Z_j)} x−−−>Q(Zj​)p(xi​∣θk​)​

p ( x ) − − − > Q ( Z j ) p(x)--->Q(Z_j) p(x)−−−>Q(Zj​)

  是以有：

∑ i = 1 N l o g ∑ k = 1 K Q ( Z j ) p ( x i ∣ θ k ) Q ( Z j ) ≥ ∑ i = 1 N ∑ k = 1 K l o g ( p ( x i ∣ θ k ) Q ( Z j ) ) Q ( Z j ) = J ( Z , θ ) (5) \sum_{i=1}^N log\sum_{k=1}^KQ(Z_j)\frac{p(x_i|\theta_k)}{Q(Z_j)} \ge\sum_{i=1}^N\sum_{k=1}^Klog(\frac{p(x_i|\theta_k)}{Q(Z_j)})Q(Z_j)=J(Z,\theta)\tag{5} i=1∑N​logk=1∑K​Q(Zj​)Q(Zj​)p(xi​∣θk​)​≥i=1∑N​k=1∑K​log(Q(Zj​)p(xi​∣θk​)​)Q(Zj​)=J(Z,θ)(5)

  由 ( 5 ) (5) (5)式可以看出 l o g L ( θ ) ≥ J ( Z , θ ) logL(\theta)\ge J(Z,\theta) logL(θ)≥J(Z,θ)，顯然 J ( Z , θ ) J(Z,\theta) J(Z,θ)是似然函數的下界，隻要不斷提高 J ( Z , θ ) J(Z,\theta) J(Z,θ)，就能不斷提高 l o g L ( θ ) logL(\theta) logL(θ)。這裡有個問題：

• 如何提高 J ( Z , θ ) J(Z,\theta) J(Z,θ)？？？

J ( Z , θ ) = ∑ i = 1 N ∑ k = 1 K l o g ( p ( x i ∣ θ k ) Q ( Z j ) ) Q ( Z j ) J(Z,\theta)=\sum_{i=1}^N\sum_{k=1}^Klog(\frac{p(x_i|\theta_k)}{Q(Z_j)})Q(Z_j) J(Z,θ)=i=1∑N​k=1∑K​log(Q(Zj​)p(xi​∣θk​)​)Q(Zj​)

  顯然需要把 Q ( z j ) Q(z_j) Q(zj​)先求出來（EM算法的E步），然後求導估計 θ k \theta_k θk​(EM算法的M步)。

3.E步

  因為 Z Z Z是隐變量，也是一個随機變量，是以有:

∑ z Q ( Z j ) = 1 \sum_{z}Q(Z_j)=1 z∑​Q(Zj​)=1

  當琴聲不等式取"="号時，有：

p ( x i ∣ θ k ) Q ( Z j ) = c = ∑ k = 1 K p ( x i ∣ θ k ) ∑ z Q ( Z j ) = ∑ k = 1 K p ( x i ∣ θ k ) \frac{p(x_i|\theta_k)}{Q(Z_j)}=c=\frac{\sum_{k=1}^Kp(x_i|\theta_k)}{\sum_{z}Q(Z_j)}=\sum_{k=1}^Kp(x_i|\theta_k) Q(Zj​)p(xi​∣θk​)​=c=∑z​Q(Zj​)∑k=1K​p(xi​∣θk​)​=k=1∑K​p(xi​∣θk​)

  是以 Q ( Z j ) = p ( x i ∣ θ k ) ∑ k = 1 K p ( x i ∣ θ k ) Q(Z_j)=\frac{p(x_i|\theta_k)}{\sum_{k=1}^Kp(x_i|\theta_k)} Q(Zj​)=∑k=1K​p(xi​∣θk​)p(xi​∣θk​)​

  從這裡可以看出，如果 p ( x i ∣ θ k ) = p ( x i + 1 ∣ θ k ) p(x_i|\theta_k)=p(x_{i+1}|\theta_k) p(xi​∣θk​)=p(xi+1​∣θk​),那麼 Q ( Z j ) = Q ( Z j + 1 ) Q(Z_j)=Q(Z_{j+1}) Q(Zj​)=Q(Zj+1​)。即對于同分布的樣本，其隐變量取值相同的機率是一樣的。

  通過初始化 θ k \theta_k θk​，那麼就可先估計出 Q ( Z j ) Q(Z_j) Q(Zj​)。這是E步。

4.M步

  當 Q ( Z j ) Q(Z_j) Q(Zj​)已知，帶入 J ( Z , θ ) J(Z,\theta) J(Z,θ)，令 d J ( z , θ ) d θ = 0 \frac{dJ(z,\theta)}{d{\theta}}=0 dθdJ(z,θ)​=0，估計 θ \theta θ。

  重複E步與M步，直到算法收斂，這就是EM算法啦~

高斯混合模型（GMM）

  高斯混合模型是利用多個高斯分布對任意的分布進行拟合。最大的問題也是參數估計，用的方法也是EM算法。

1.極大似然估計

  設有樣本 X = { X 1 , X 2 . . . . X n } X=\{X_1,X_2....X_n\} X={X1​,X2​....Xn​}，對樣本 X i X_i Xi​

X i ∼ ∑ k = 1 K π k N ( μ k , Σ k ) X_i\sim\sum_{k=1}^K\pi_kN(\mu_k,\Sigma_k) Xi​∼k=1∑K​πk​N(μk​,Σk​)

  即利用 K K K個高斯分布進行拟合。

  對于樣本 X i X_i Xi​引入 K K K維隐變量 Z j = [ z j 1 z j 2 ⋮ z j k ] Z_j=\begin{bmatrix} z_{j1} \\ z_{j2}\\\vdots \\z_{jk}\end{bmatrix} Zj​=⎣⎢⎢⎢⎡​zj1​zj2​⋮zjk​​⎦⎥⎥⎥⎤​，其中 ∑ k = 1 K z j k = 1 \sum_{k=1}^Kz_{jk}=1 ∑k=1K​zjk​=1，且隻有一個分量能等于 1 1 1，代表該高斯分量的機率最大。

  下面要混合高斯分布寫成包含隐變量的形式，需要用到機率公式。

  對于 Q ( z j k ) Q(z_{jk}) Q(zjk​)

Q ( z j k ) = π k , z j k = 1 Q(z_{jk})=\pi_k,z_{jk}=1 Q(zjk​)=πk​,zjk​=1

Q ( z j k ) = 1 , z j k = 0 Q(z_{jk})=1,z_{jk}=0 Q(zjk​)=1,zjk​=0

是以 Q ( z j k ) = π k z j k (6) Q(z_{jk})=\pi_k^{z_{jk}}\tag{6} Q(zjk​)=πkzjk​​(6)

對于 p ( X i , Z j ) p(X_i,Z_j) p(Xi​,Zj​)有：

p ( X i , Z j ) = ∏ k = 1 K p ( X i ∣ z j k ) Q ( z j k ) (7) p(X_i,Z_j)=\prod_{k=1}^Kp(X_i|z_{jk})Q(z_{jk})\tag{7} p(Xi​,Zj​)=k=1∏K​p(Xi​∣zjk​)Q(zjk​)(7)

( 7 ) (7) (7)式用到了全機率公式，且：

p ( X i ∣ z j k ) = N ( X i ∣ μ k , Σ k ) z j k (8) p(X_i|z_{jk})=N(X_i|\mu_k,\Sigma_k)^{z_{jk}}\tag{8} p(Xi​∣zjk​)=N(Xi​∣μk​,Σk​)zjk​(8)

把 ( 6 ) ( 7 ) (6)(7) (6)(7)代入 ( 8 ) (8) (8),帶隐變量的GMM分布為：

p ( X i , Z j ) = ∏ k = 1 K N ( X i ∣ μ k , Σ k ) z j k π k z j k (9) p(X_i,Z_j)=\prod_{k=1}^KN(X_i|\mu_k,\Sigma_k)^{z_{jk}}\pi_k^{z_{jk}}\tag{9} p(Xi​,Zj​)=k=1∏K​N(Xi​∣μk​,Σk​)zjk​πkzjk​​(9)

帶隐變量的似然函數為：

L ( X , Z , μ , Σ ) = ∏ i = 1 N ∏ k = 1 K N ( X i ∣ μ k , Σ k ) z j k π k z j k (10) L(X,Z,\mu,\Sigma)=\prod_{i=1}^N\prod_{k=1}^KN(X_i|\mu_k,\Sigma_k)^{z_{jk}}\pi_k^{z_{jk}}\tag{10} L(X,Z,μ,Σ)=i=1∏N​k=1∏K​N(Xi​∣μk​,Σk​)zjk​πkzjk​​(10)

取對數得：

l n ( L ( X , Z , μ , Σ ) ) = ∑ k = 1 K ( ∑ j = 1 N z j k ) l n ( π k ) + ∑ i = 1 N z j k l n ( N ( X i ∣ μ k , Σ k ) ) (11) ln(L(X,Z,\mu,\Sigma))=\sum_{k=1}^K(\sum_{j=1}^Nz_{jk})ln(\pi_k)+\sum_{i=1}^Nz_{jk}ln(N(X_i|\mu_k,\Sigma_k))\tag{11} ln(L(X,Z,μ,Σ))=k=1∑K​(j=1∑N​zjk​)ln(πk​)+i=1∑N​zjk​ln(N(Xi​∣μk​,Σk​))(11)

2.E步

  估計 Z j Z_j Zj​,根據EM算法的E步有：

E ( Z j ) = Q ( Z j ) Z j = π k N ( X i ) ∑ k = 1 k π k N ( X i ) E(Z_j)=Q(Z_j)Z_j=\frac{\pi_kN(X_i)}{\sum_{k=1}^k\pi_kN(X_i)} E(Zj​)=Q(Zj​)Zj​=∑k=1k​πk​N(Xi​)πk​N(Xi​)​

  其中， E ( Z j ) E(Z_j) E(Zj​)是一個 K K K維的向量。這個式子可以直接根據EM算法導出，也可以根據機率公式進行計算。

2.M步

  有了E步，即知道了 Z j Z_j Zj​，那麼可以直接對 ( 11 ) (11) (11)式子求導了，具體的求導運算不寫了，直接給結果吧。

μ k = ∑ i = 1 N z j k X i ∑ i = 1 N z j k \mu_k=\frac{\sum_{i=1}^Nz_{jk}X_i}{\sum_{i=1}^Nz_{jk}} μk​=∑i=1N​zjk​∑i=1N​zjk​Xi​​

Σ k = ∑ i = 1 N z j k ( X i − μ k ) T ( X i − μ k ) ∑ i = 1 N z j k \Sigma_k=\frac{\sum_{i=1}^Nz_{jk}(X_i-\mu_k)^T(X_i-\mu_k)}{\sum_{i=1}^Nz_{jk}} Σk​=∑i=1N​zjk​∑i=1N​zjk​(Xi​−μk​)T(Xi​−μk​)​

π k = ∑ i = 1 N z j k N \pi_k=\frac{\sum_{i=1}^Nz_{jk}}{N} πk​=N∑i=1N​zjk​​

3.matlab仿真

clear;%%以下用GMM模型進行聚類

%首先産生資料,%%産生了80個資料點，均服從二維高斯分布
X1=normrnd(3,1,[100 2]);
X2=normrnd(6,1.5,[40 2]);
X=[X1',X2'];  

%初始化相應的參數；
N=length(X);
beta=0.001;%正則化因子
% Pai=0.5*ones(2,1);%%這是混合GMM的分量
Pai(1)=0.2;Pai(2)=0.8;

mu_1=rand(1)*ones(2,1);%産生均值矩陣
mu_2=rand(1)*ones(2,1); 

Var_1=rand(1)*eye(2);%産生協方差矩陣
Var_2=rand(1)*eye(2); 

for iter=1:20
    
    
%E步
    for i=1:N
        
    N_x_1=mvnpdf(X(:,i),mu_1,Var_1);%對應機率密度的值
    N_x_2=mvnpdf(X(:,i),mu_2,Var_2);
    
    EZ(1,i)=Pai(1)*N_x_1/(Pai(1)*N_x_1+Pai(2)*N_x_2); %估計隐變量   
    EZ(2,i)=Pai(2)*N_x_2/(Pai(1)*N_x_1+Pai(2)*N_x_2); 
      
    end
   
%M步

%1)更新均值向量    
N_k1=sum(EZ(1,:));
N_k2=sum(EZ(2,:)); %隐變量對于分量的和
m1=0;m2=0;
    for i=1:N  
        m1=m1+EZ(1,i)*X(:,i); 
        
        m2=m2+EZ(2,i)*X(:,i); 
    end   
mu_1=m1/N_k1;
mu_2=m2/N_k2; 

%2)更新協方差矩陣

V1=0;V2=0; %代表臨時的方差矩陣

    for i=1:N         
    V1=V1+EZ(1,i)*(X(:,i)-mu_1)*(X(:,i)-mu_1)';
    V2=V2+EZ(2,i)*(X(:,i)-mu_2)*(X(:,i)-mu_2)'; %%更新方差矩陣               
    end
    
    Var_1=V1/N_k1+beta*eye(2);
    Var_2=V2/N_k2+beta*eye(2); 
    
%3)更新pai值
Pai(1)=N_k1/N;
Pai(2)=N_k2/N;

    
    
end
    


%%看看效果
for i=1:N
    if(EZ(1,i)>=EZ(2,i))
        rx1(:,i)=X(:,i);
    else
        rx2(:,i)=X(:,i);
    end
end


subplot(2,1,1);
scatter(X1(:,1),X1(:,2),'r');
hold on;
scatter(X2(:,1),X2(:,2),'g');
xlabel('實際分類');
hold off;


subplot(2,1,2);
scatter(rx1(1,:),rx1(2,:),'r');
hold on;
scatter(rx2(1,:),rx2(2,:),'g');
xlabel('GMM模型聚類')
hold off;
           

結果：

K-means算法

  K-means隻是GMM模型的退化，弄明白了EM和GMM算法，K-means很容易了解。

1.理論分析

  目标函數如下：

J = ∑ i = 1 N ∑ k = 1 K ∣ ∣ ( X i − μ k ) ∣ ∣ J=\sum_{i=1}^N \sum_{k=1}^K||(X_i-\mu_k)|| J=i=1∑N​k=1∑K​∣∣(Xi​−μk​)∣∣

  E步，經過初始後，先大概分類：

d k = ∣ ∣ X i − μ k ∣ ∣ d_k=||X_i-\mu_k|| dk​=∣∣Xi​−μk​∣∣

  對于 X i X_i Xi​， d k d_k dk​最大，則劃分為第 k k k類，這個過程也可以用GMM模型中的隐變量表示。

  M步，估計 μ k \mu_k μk​

μ k = ∑ k = 1 N k X k N k \mu_k=\frac{\sum_{k=1}^{N_k}{X_k}}{N_k} μk​=Nk​∑k=1Nk​​Xk​​

  表達了對是以的 k k k類資料求取均值。

2.matlab仿真

clear;%%以下用K-means模型進行聚類

%首先産生資料,%%産生了80個資料點，均服從二維高斯分布
X1=normrnd(3,1,[100 2]);
X2=normrnd(6,1.5,[100 2]);
X=[X1',X2'];  

%初始化相應的參數；
N=length(X);
iter=20;
x1_=[0,2]';
x2_=[3,3]'; %均值（質心的初始值）

for j=1:iter
    %E步
    for i=1:N
    d_1=norm(X(:,i)-x1_);
    d_2=norm(X(:,i)-x2_);
    if(d_1<d_2)
        z(:,i)=[1;0];
    else
        z(:,i)=[0;1];
    end
    
    end


%M步
%分别計算更新後的x1_,x2_
tmp_1=0;tmp_2=0;
N_1=0;N_2=0;
    for(i=1:N)
        if( z(:,i)==[1;0])
            tmp_1=tmp_1+X(:,i);
            N_1=N_1+1;
        else
            tmp_2=tmp_2+X(:,i);
            N_2=N_2+1;
        end

    end
x1_=tmp_1/N_1;
x2_=tmp_2/N_2;

%%看看誤差
    e(j)=0;
    for(i=1:N)  
        e(j)=e(j)+norm(X(:,i)-x1_)+norm(X(:,i)-x2_);
    end
    
end

%畫圖

%%看看效果
for i=1:N
    if(z(:,i)==[1,0]')
        rx1(:,i)=X(:,i);
    else
        rx2(:,i)=X(:,i);
    end
end


subplot(3,1,1);
scatter(X1(:,1),X1(:,2),'r');
hold on;
scatter(X2(:,1),X2(:,2),'g');
xlabel('實際分類');
hold off;


subplot(3,1,2);
scatter(rx1(1,:),rx1(2,:),'r');
hold on;
scatter(rx2(1,:),rx2(2,:),'g');
xlabel('K-means算法聚類');
hold off;
subplot(3,1,3)
m=1:1:iter;
scatter(m',e);
xlabel('誤差變化');
           

最後要感謝幾位大神的部落格，受益匪淺！！

參考文獻：

[1]https://blog.csdn.net/jinping_shi/article/details/59613054

[2]https://blog.csdn.net/lin_limin/article/details/81048411?utm_medium=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.add_param_isCf&depth_1-utm_source=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.add_param_isCf

[3]https://blog.csdn.net/zouxy09/article/details/8537620