斯坦福機器學習實現與分析之五（高斯判別分析）

高斯判別分析（GDA）簡介

---

　　首先，高斯判別分析的作用也是用於分類。對於兩類樣本，其服從伯努利分布，而對每個類中的樣本，假定都服從高斯分布，則有:

\( y\;\sim\;Bernouli(\phi) \)

\( x|y=0\;\sim\;N(\mu_0, \Sigma) \)

\( x|y=1\;\sim\;N(\mu_1, \Sigma) \)

　　這樣，根據訓練樣本，估計出先驗概率以及高斯分布的均值和協方差矩陣（注意這里兩類內部高斯分布的協方差矩陣相同），即可通過如下貝葉斯公式求出一個新樣本分別屬於兩類的概率，進而可實現對該樣本的分類。

　\(\begin{aligned} p(y|x)=\frac{p(x|y)p(y)}{p(x)} \end{aligned}\)

\( \begin{aligned} y=\underset{y}{argmax}\;{p(y|x)} =  \underset{y}{argmax}{ \;{\frac{p(x|y)p(y)}{p(x)}}}= \underset{y}{argmax}{\;p(x|y)p(y)} \end{aligned}\) 

 

GDA詳細推導

---

　　那么高斯判別分析的核心工作就是估計上述未知量\( \phi, \mu_0, \mu_1, \Sigma \)。如何來估計這些參數？又該最大似然估計上場了。其對數似然函數為：

\( \begin{aligned} l(\phi,\mu_0,\mu_1,\Sigma) &=log{\prod_{i=1}^m{p(x^{(i)},y^{(i)})}}=log{\prod_{i=1}^m{p(x^{(i)}|y^{(i)})p(y^{(i)})}} \\  &=\sum_{i=1}^m{log\;p(x^{(i)}|y^{(i)})}+\sum_{i=1}^m{log\;p(y^{(i)})} \\ &=\sum_{i=1}^m{log\;(\;p(x^{(i)}|y^{(i)}=0)^{1-y^{(i)}}*p(x^{(i)}|y^{(i)}=1)^{y^{(i)}}\;)}+\sum_{i=1}^m{log\;p(y^{(i)})} \\ &=\sum_{i=1}^m{(1-y^{(i)})log\;p(x^{(i)}|y^{(i)}=0)}+\sum_{i=1}^m{{y^{(i)}}log\;p(x^{(i)}|y^{(i)}=1)}+\sum_{i=1}^m{log\;p(y^{(i)})} \end{aligned}\)

　　注意此函數第一部分只和\(\mu_0,\Sigma\)有關，第二部分只和\(\mu_1,\Sigma\)有關，最后一部分只和\(\phi\)有關。最大化該函數，首先求\(\phi\),先對其求偏導數：

\(\begin{aligned} \frac{\partial\;l(\phi,\mu_0,\mu_1,\Sigma)}{\partial\phi}&=\frac{\sum_{i=1}^m{log\;p(y^{(i)})}}{\partial\phi} \\&= \frac{\partial\sum_{i=1}^m{log\;\phi^{y^{(i)}}(1-\phi)^{1-y^{(i)}})}}{\partial\phi} \\&=\frac{\partial\sum_{i=1}^m{y^{(i)}\;log\;\phi+(1-y^{(i)})log(1-\phi)}}{\partial\phi} \\&=\sum_{i=1}^m{(y^{(i)}\frac{1}{\phi}-(1-y^{(i)})\frac{1}{1-\phi})} \\&=\sum_{i=1}^m{(I(y^{(i)}=1)\frac{1}{\phi}-I(y^{(i)}=0)\frac{1}{1-\phi})} \end{aligned} \)

　　此處\(I\)為指示函數。令其為0，可求解出：

\(\begin{aligned} \phi=\frac{I(y^{(i)}=1)}{I(y^{(i)}=0)+I(y^{(i)}=1)}=\frac{I(y^{(i)}=1)}{m}\end{aligned}\) 

　　同樣地，對\(\mu_0\)求偏導數：

\(\begin{aligned} \frac{\partial\;l(\phi,\mu_0,\mu_1,\Sigma)}{\partial\mu_0}&=\frac{\partial\sum_{i=1}^m{(1-y^{(i)})log\;p(x^{(i)}|y^{(i)}=0)}}{\partial\mu_0} \\&=\frac{\partial\sum_{i=1}^m{(1-y^{(i)})(log\frac{1}{\sqrt{(2\pi)^n|\Sigma|}}-\frac{1}{2}(x^{(i)}-\mu_0)^T\Sigma^{-1}(x^{(i)}-\mu_0))}}{\partial\mu_0} \\&=\sum_{i=1}^m{(1-y^{(i)})(\Sigma^{-1}(x^{(i)}-\mu_0))} \\&=\sum_{i=1}^m{I(y^{(i)}=0)\Sigma^{-1}(x^{(i)}-\mu_0)}
\end{aligned}\) 

　　令其為0，可求解得：

\(\begin{aligned} \mu_0=\frac{\sum_{i=1}^m{I(y^{(i)}=0)x^{(i)}}}{\sum_{i=1}^m{I(y^{(i)}=0)}} \end{aligned}\) 

　　根據對稱性可直接得出：

\(\begin{aligned} \mu_1=\frac{\sum_{i=1}^m{I(y^{(i)}=1)x^{(i)}}}{\sum_{i=1}^m{I(y^{(i)}=1)}} \end{aligned}\) 

　　下面對\( \Sigma \)求偏導數，由於似然函數只有前面兩部分與\( \Sigma \)有關，則將前兩部分改寫如下：

\( \begin{aligned} &\sum_{i=1}^m{(1-y^{(i)})log\;p(x^{(i)}|y^{(i)}=0)}+\sum_{i=1}^m{{y^{(i)}}log\;p(x^{(i)}|y^{(i)}=1)}\\&=\sum_{i=1}^m{(1-y^{(i)})(log\frac{1}{\sqrt{(2\pi)^n|\Sigma|}}-\frac{1}{2}(x^{(i)}-\mu_0)^T\Sigma^{-1}(x^{(i)}-\mu_0))}+\sum_{i=1}^m{{y^{(i)}}(log\frac{1}{\sqrt{(2\pi)^n|\Sigma|}}-\frac{1}{2}(x^{(i)}-\mu_1)^T\Sigma^{-1}(x^{(i)}-\mu_1))}\\&=\sum_{i=1}^m{log\frac{1}{\sqrt{(2\pi)^n|\Sigma|}}}-\frac{1}{2}\sum_{i=1}^m{(x^{(i)}-\mu_{y^{(i)}})^T\Sigma^{-1}(x^{(i)}-\mu_{y^{(i)}})}\\&=\sum_{i=1}^m{(-\frac{n}{2}log(2\pi)-\frac{1}{2}log(|\Sigma|))}-\frac{1}{2}\sum_{i=1}^m{(x^{(i)}-\mu_{y^{(i)}})^T\Sigma^{-1}(x^{(i)}-\mu_{y^{(i)}})}   \end{aligned}\)

　　進而有：

\( \begin{aligned} \frac{\partial\;l(\phi,\mu_0,\mu_1,\Sigma))}{\partial\Sigma}&=-\frac{1}{2}\sum_{i=1}^m(\frac{1}{|\Sigma|}|\Sigma|\Sigma^{-1})-\frac{1}{2}\sum_{i=1}^m(x^{(i)}-\mu_{y^{(i)}})(x^{(i)}-\mu_{y^{(i)}})^T\frac{\partial\Sigma^{-1}}{\partial\Sigma}\\&=-\frac{m}{2}\Sigma^{-1}-\frac{1}{2}\sum_{i=1}^m(x^{(i)}-\mu_{y^{(i)}})(x^{(i)}-\mu_{y^{(i)}})^T(-\Sigma^{-2}))   \end{aligned}\)

　　這里推導用到了：

\( \begin{aligned} \frac{\partial|\Sigma|}{\partial\Sigma}=|\Sigma|\Sigma^{-1}\end{aligned}\)

\( \begin{aligned} \frac{\partial\Sigma^{-1}}{\partial\Sigma}=-\Sigma^{-2}\end{aligned}\) 

　　令其為0，從而求得：

\( \begin{aligned} \Sigma=\frac{1}{m}\sum_{i=1}^m(x^{(i)}-\mu_{y^{(i)}})(x^{(i)}-\mu_{y^{(i)}})^T\end{aligned}\) 

　　上面的推導似乎很復雜，但其結果卻是非常簡潔。通過上述公式，所有的參數都已經估計出來，需要判斷一個新樣本x時，可分別使用貝葉斯求出p(y=0|x)和p(y=1|x)，取概率更大的那個類。

　　實際計算時，我們只需要比大小，那么貝葉斯公式中分母項可以不計算，由於2個高斯函數協方差矩陣相同，則高斯分布前面那相同部分也可以忽略。實際上，GDA算法也是一個線性分類器，根據上面推導可以知道，GDA的分界線(面)的方程為：

        \( \begin{aligned} (1-\phi)exp((x-\mu_0)^T\Sigma^{-1}(x-\mu_0)={\phi}exp((x-\mu_1)^T\Sigma^{-1}(x-\mu_1)\end{aligned}\)

　　取對數展開后化解，可得：

 \( \begin{aligned} 2x^T\Sigma^{-1}(\mu_1-\mu_0)=\mu_1^T\Sigma^{-1}\mu_1-\mu_0^T\Sigma^{-1}\mu_0+log\;\phi-log(1-\phi)\end{aligned}\)

　　若\( \begin{aligned} A=2\Sigma^{-1}(\mu_1-\mu_0)=(a_1,a_2,...,a_n)\quad b=\mu_1^T\Sigma^{-1}\mu_1-\mu_0^T\Sigma^{-1}\mu_0+log\;\phi-log(1-\phi)\end{aligned}\)，則

 \( \begin{aligned} a_1x_1+a_2x_2+...+a_nx_n=b\end{aligned}\)

　　這就是GDA算法的線性分界面。

 

GDA實現

---

　　這里也采用前面講邏輯回歸生成的數據來進行實驗，直接load進來進行處理，詳見邏輯回歸。GDA訓練代碼如下： 

%mu=[mu0 mu1]
function [mu sigma phi]=GDA_train(Sample)
    [m, n] = size(Sample); %m個樣本，每個n維
    Y = Sample(:, end);
    X = [Sample(:,1:end-1)];
    
    idx = find(Y==0);
    mu(:,1)=mean(X(idx,:));
    
    idx2 = find(Y==1);    
    mu(:,2)=mean(X(idx2,:));
    
    phi = size(idx2,1)/m;
    
    sigma = zeros(n-1);
    for i = 1:m
        x = X(i, :)';
        muc = mu(:, Y(i) + 1);
        sigma = sigma + (x - muc) * (x - muc)';
    end
    sigma = sigma / m;    
end

 　　測試代碼： 

function GDA_test
    clear
    close all
    clc
    load('log_data.mat');
    [mu sigma phi]=GDA_train(Sample);
    
    %顯示結果，以下代碼不通用，樣本維數增加時顯示不可用
    figure,
    idx = find(Sample(:,3)==1);
    plot(Sample(idx,1), Sample(idx,2), 'g*');hold on
    idx = find(Sample(:,3)==0);
    plot(Sample(idx,1), Sample(idx,2), 'ro');hold on
    
    [t1 t2]=meshgrid(min(Sample(:,1)):.1:max(Sample(:,1)), min(Sample(:,2)):.1:max(Sample(:,2)));
    G1=Gaussian(mu(:,1),sigma,t1,t2,1-phi);
    G2=Gaussian(mu(:,2),sigma,t1,t2,phi);
    contour(t1,t2,G1);hold on
    contour(t1,t2,G2);hold on
    
    A=2*inv(sigma)*(mu(:,2)-mu(:,1));
    b=mu(:,2)'*inv(sigma)*mu(:,2) - mu(:,1)'*inv(sigma)*mu(:,1) + log(phi) - log(1-phi);
    x1=min(Sample(:,1)):.1:max(Sample(:,1));
    x2=(b-A(1)*x1)/A(2);
    plot(x1,x2,'m')
end

function G = Gaussian(mu,sigma,t1,t2,phi)
    for i=1:size(t1,1)
        for j=1:size(t1,2)
            x=[t1(i,j);t2(i,j)];
            z=(x-mu)'*inv(sigma)*(x-mu);
            G(i,j)=phi*exp(-z)/0.001;
        end
    end    
end

　　訓練結果如下，訓練樣本中，正負樣本均為100個，故\(\phi=0.5\)：

 

　　改變正負樣本數量，即相當於改變先驗概率，則實驗結果如下(相應的\(\phi\)的值顯示在圖像標題)：

       

 

算法分析

---

　　1.與邏輯回歸的關系

　　　　根據上面的結果以及貝葉斯公式，可有

 \( \begin{aligned} p(y=1|x)&=\frac{p(x|y=1)p(y=1)}{p(x)} \\&=\frac{N(\mu_1,\Sigma)\phi}{N(\mu_0,\Sigma)(1-\phi)+N(\mu_1,\Sigma)\phi}\\&=1/{(1+\frac{N(\mu_0,\Sigma))}{N(\mu_1,\Sigma)}\frac{1-\phi}{\phi})} \end{aligned}\)

　　　　而

\( \begin{aligned} \frac{N(\mu_0,\Sigma)}{N(\mu_1,\Sigma)}&= exp\{(x-\mu_0)^T\Sigma^{-1}(x-\mu_0)-(x-\mu_1)^T\Sigma^{-1}(x-\mu_1)\}\\&=exp\{ 2(\mu_1-\mu_0)^T\Sigma^{-1}x+(\mu_0^T\Sigma\mu_0-\mu_1^T\Sigma\mu_1)\}\end{aligned}\)

　　　　那么,令

\( \begin{aligned}
2\Sigma^{-1}(\mu_1-\mu_0) =(\theta_1,\theta_2,...,\theta_n)^T\\
\theta_0=\mu_0^T\Sigma\mu_0-\mu_1^T\Sigma\mu_1+log\frac{1-\phi}{\phi}\\
\end{aligned}\)

　　　　則

\( \begin{aligned}
p(y=1|x)=\frac{1}{1+exp(\theta_0+\theta_1x_1+\theta_2x_2+...+\theta_nx_n)}
\end{aligned}\)

　　　　這不就是邏輯回歸的形式么？

　　　　在推導邏輯回歸的時候，我們並沒有假設類內樣本是服從高斯分布的，因而GDA只是邏輯回歸的一個特例，其建立在更強的假設條。故兩者效果比較：

　　　　a.邏輯回歸是基於弱假設推導的，則其效果更穩定，適用范圍更廣

　　　　b.數據服從高斯分布時，GDA效果更好

　　　　c.當訓練樣本數很大時，根據中心極限定理，數據將無限逼近於高斯分布，則此時GDA的表現效果會非常好

 

　　2.為何要假設兩類內部高斯分布的協方差矩陣相同？

　　　　從直觀上講，假設兩個類的高斯分布協方差矩陣不同，會更加合理（在混合高斯模型中就是如此假設的），而且可推導出類似上面簡潔的結果。

　　　　假定兩個類有相同協方差矩陣，分析具有以下幾點影響：

　　　　A．當樣本不充分時，使用不同協方差矩陣會導致算法穩定性不夠；過少的樣本甚至導致協方差矩陣不可逆，那么GDA算法就沒法進行

　　　　B．使用不同協方差矩陣，最終GDA的分界面不是線性的，同樣也推導不出GDA的邏輯回歸形式

 

　　3.使用GDA時對訓練樣本有何要求？

　　　　首先，正負樣本數的比例需要符合其先驗概率。若是預先明確知道兩類的先驗概率，那么可使用此概率來代替GDA計算的先驗概率；若是完全不知道，則可以公平地認為先驗概率為　　50%。

　　　　其次，樣本數必須不小於樣本特征維數，否則會導致協方差矩陣不可逆，按照前面分析應該是多多益善。