高斯判别分析（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训练代码如下：
1 %mu=[mu0 mu1] 2 function [mu sigma phi]=GDA_train(Sample) 3 [m, n] = size(Sample); %m个样本，每个n维 4 Y = Sample(:, end); 5 X = [Sample(:,1:end-1)]; 6 7 idx = find(Y==0); 8 mu(:,1)=mean(X(idx,:)); 9 10 idx2 = find(Y==1); 11 mu(:,2)=mean(X(idx2,:)); 12 13 phi = size(idx2,1)/m; 14 15 sigma = zeros(n-1); 16 for i = 1:m 17 x = X(i, :)'; 18 muc = mu(:, Y(i) + 1); 19 sigma = sigma + (x - muc) * (x - muc)'; 20 end 21 sigma = sigma / m; 22 end