高斯判別分析GDA推導與代碼實現

生成學習

處理分類問題，我們可以使用邏輯回歸、Softmax。這兩種方法都屬于“判別學習”，也就是給定 \((x^{(i)}, y^{(i)})\)，我們學習 \(P(y|x)\)，并對于給定的 \(x\)，計算 \(\operatorname{argmax}_{y}\{P(y|x)\}\)。

GDA屬于另一種方法——生成學習。在判別學習中，我們并不關注 \(x\) 本身的分布，而在生成學習中，我們基于一些事實假設 \(x\) 的分布，例如在GDA中，我們假設對于相同的 \(y\)，\(x\) 符合高斯分布。然后基于假設，學習 \(P(x|y)\) 的參數（高斯分布則學習均值 \(\mu\) 和協方差 \(\Sigma\)），以及 \(P(y)\) 的參數（例如兩點分布的 \(\phi\)），從而得到聯合概率 \(P(x,y)\)。在預測時，我們仍然是找到最大的 \(P(y|x)\)，但是是用貝葉斯公式計算：

高斯判別分析

作出的假設

首先假設 \(Y\) 服從兩點分布 \(Bernolli(\phi)\)，于是可以寫作 \(P(y)=\phi^{y}(1-\phi)^{1-y}\)。

（盡管可以直接用協方差定義多元變量的高斯分布，但是這里采用另一種方法，在特殊的情況下得到等式而不需要完全理解協方差矩陣。）

然后假設給定 \(y\) 后，\(x\) 的每一個分量 \(x_i\) 都服從高斯分布且相互獨立。不妨設 \(y=0\)，即 \(x_i|y=0\sim N(\mu_i,\sigma_i)\)：

又 \(x_i\) 相互獨立，可以得到：

我們定義矩陣 \(\Sigma=\begin{pmatrix}\sigma_1^2\\&\ddots\\&&\sigma_n^2\end{pmatrix}\)，可以將上式化簡：

同樣，我們也假設 \(x|y=1\) 符合高斯分布 \(N(\mu',\Sigma)\)，需要注意的是，兩個分布都采用了同一個 \(\Sigma\)（存疑，不知道目的）。

但是實際上GDA并沒有假設 \(\Sigma\)

最大似然估計

在判別學習中，我們以 \(P(y^{(i)}\mid x^{i})\) 為似然函數，而在GDA這一類生成學習中，我們以聯合概率 \(P(x^{(i)},y^{(i)})\) 為似然函數。也即最大化如下對數似然函數：

不同于邏輯回歸，我們可以直接用導數為 \(0\) 求解參數。

計算 \(\phi\)

令上式為 \(0\)，得 \(\phi=\frac{\sum y^{(i)}}{m}\)。

計算 \(\mu,\mu'\)

一些無關緊要的常數用 \(C\) 代替，比如密度函數中的系數。

其中

則 \(\nabla_{\mu}\ln P(x^{(i)}\mid y=0;\mu,\Sigma)=x^{(i)}-\mu\)，進而

令上式為 \(0\)，可得 \(\mu=\frac{\sum(1-y^{(i)})x^{(i)}}{\sum 1-y^{(i)}}\)，同理可得 \(\mu'=\frac{\sum y^{(i)}x^{(i)}}{\sum y^{(i)}}\)。

計算 \(\Sigma\)

下面用到了兩個 \(\nabla\) 算符的性質：

• \(\nabla_{\Sigma}|\Sigma|=|\Sigma|\Sigma^{-1}\)；
• \(\nabla_{\Sigma}\Sigma^{-1}=-\Sigma^{-2}\)。

同樣地令上式為 \(0\)，計算得

注意到，按照我們的假設，\(\Sigma\) 應該是一個對角矩陣（如果了解協方差矩陣的話，由 \(x_i\) 相互獨立，可以推出 \(\Sigma\) 應該是對角矩陣，對角元就是每個變量的方差），但是這里非常顯然不總是對角矩陣。

（以下僅是我個人的猜測）在GDA的實現中，我們并沒有關注 \(x_i\) 相互獨立的性質，而是直接學習了一個普遍的協方差矩陣。實際上這種定義的限制更低，更符合現實情況（現實中的變量之間存在聯系比較普遍）。

代碼實現

直接計算這幾個參數的代碼不需要解釋：

m, n = xs.shape
# Calculate mu
self.mu = np.zeros((2, n))
classes_size = np.zeros(2)
for i in range(m):
    self.mu[ys[i]] += xs[i]
    classes_size[ys[i]] += 1
self.mu /= np.transpose([classes_size])
# Calculate Sigma
self.sigma = np.zeros((n, n))
for i in range(m):
    temp_array = xs[i] - self.mu[ys[i]]
    self.sigma += np.dot(temp_array.reshape(n, 1), temp_array.reshape(1, n))
self.sigma /= m
# Calculate phi
self.phi = np.sum(ys) / m

最后我們發現計算概率需要計算 \(\Sigma^{-1}\)，然而我們算出來的可能是一個奇異矩陣。這時候可以對 \(\Sigma\) 進行微擾——將對角元加上一個較小的偏差。

最后我們對 \(y=0,1\) 分別計算出 \(P(x|y)\)。

但是實際上比較 \(P(x|y=0),P(x|y=1)\) 不需要真正計算出概率，而是比較兩者不同的地方，也就是后面的 \(\exp\)。

（奇異矩陣的處理參考 https://blog.csdn.net/qq_30091945/article/details/81508055 ，對其中的 Gaussian 函數有修改）

最后是實現的一個類：

class GaussianDiscriminantAnalysis:
    def __init__(self):
        self.mu = None
        self.phi = None
        self.sigma = None

    def fit(self, xs, ys, **others):
        m, n = xs.shape
        # Calculate mu
        self.mu = np.zeros((2, n))
        classes_size = np.zeros(2)
        for i in range(m):
            self.mu[ys[i]] += xs[i]
            classes_size[ys[i]] += 1
        self.mu /= np.transpose([classes_size])
        # Calculate Sigma
        self.sigma = np.zeros((n, n))
        for i in range(m):
            temp_array = xs[i] - self.mu[ys[i]]
            self.sigma += np.dot(temp_array.reshape(n, 1), temp_array.reshape(1, n))
        self.sigma /= m
        # Calculate phi
        self.phi = np.sum(ys) / m

    def evaluate(self, x, mean, cov):
        dim = np.shape(cov)[0]
        # cov的行列式為零時的措施
        cov_inv = np.linalg.inv(cov + np.eye(dim) * 0.001)
        xdiff = (x - mean).reshape((1, dim))
        # 概率密度
        prob = np.exp(-0.5 * xdiff.dot(cov_inv).dot(xdiff.T))[0][0]
        return prob

    def predict(self, xs):
        predict = []
        for x in xs:
            evaluate_0 = self.evaluate(x, self.mu[0], self.sigma) * (1 - self.phi)
            evaluate_1 = self.evaluate(x, self.mu[1], self.sigma) * self.phi
            if evaluate_0 > evaluate_1:
                predict.append(0)
            else:
                predict.append(1)
        return predict