高斯判别分析GDA推导与代码实现

高斯判别分析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)\)，但是是用贝叶斯公式计算：

\[P(y|x) = \frac{P(x,y)}{P(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)\)：

\[P(x_i|y=0)=\frac{1}{\sqrt{2\pi}\sigma_i}\exp\left(-\frac{(x_i-\mu_i)^2}{2\sigma_i^2}\right) \]

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

\[\begin{aligned} P(x|y=0)&=\prod_{i=1}^n\frac{1}{\sqrt{2\pi}\sigma_i}\exp\left(-\frac{(x_i-\mu_i)^2}{2\sigma_i^2}\right)\\ &=\frac{1}{(2\pi)^{\frac{n}{2}}\sigma_1\dots\sigma_n}\exp\left(-\frac{1}{2}\sum_{i=1}^n\frac{(x_i-\mu_i)^2}{\sigma_i^2}\right)\\ \end{aligned} \]

我们定义矩阵 \(\Sigma=\begin{pmatrix}\sigma_1^2\\&\ddots\\&&\sigma_n^2\end{pmatrix}\)，可以将上式化简：

\[P(x|y=0)=\frac{1}{(2\pi)^{\frac{n}{2}}|\Sigma|^{\frac{1}{2}}}\exp\left(-\frac{1}{2}(x-\mu)\Sigma^{-1}(x-\mu)^{T}\right) \]

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

但是实际上GDA并没有假设 \(\Sigma\)

最大似然估计

在判别学习中，我们以 \(P(y^{(i)}\mid x^{i})\) 为似然函数，而在GDA这一类生成学习中，我们以联合概率 \(P(x^{(i)},y^{(i)})\) 为似然函数。也即最大化如下对数似然函数：

\[\begin{aligned} \mathcal{L}(\phi,\mu,\mu',\Sigma)&=\ln\prod_{i=1}^{m}P(x^{(i)},y^{(i)};\phi,\mu,\mu',\Sigma)\\ &=\ln\prod_{i=1}^{m}P(x^{(i)}\mid y^{(i)};\mu,\mu',\Sigma)P(y^{(i)};\phi)\\ &=\sum_{i=1}^{m}\ln P(x^{(i)}\mid y^{(i)};\mu,\mu',\Sigma)+\sum_{i=1}^m\ln P(y^{(i)};\phi) \end{aligned} \]

不同于逻辑回归，我们可以直接用导数为 \(0\) 求解参数。

计算 \(\phi\)

\[\begin{aligned} \frac{\partial \mathcal{L}}{\partial \phi}&=\frac{\partial}{\partial\phi}\sum_{i=1}^m\ln P(y^{(i)};\phi)\\ &=\frac{\partial}{\partial\phi}\sum_{i=1}^my^{(i)}\ln\phi+(1-y^{(i)})\ln(1-\phi)\\ &=\sum_{i=1}^m\frac{y^{(i)}}{\phi}-\frac{1-y^{(i)}}{1-\phi} \end{aligned} \]

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

计算 \(\mu,\mu'\)

一些无关紧要的常数用 \(C\) 代替，比如密度函数中的系数。

\[\begin{aligned} \nabla_{\mu}\mathcal{L}&=\nabla_{\mu}\sum_{i=1}^m(1-y^{(i)})\ln P(x^{(i)}\mid y=0;\mu,\Sigma)+y^{(i)}\ln P(x^{(i)}\mid y=1;\mu',\Sigma) \end{aligned} \]

其中

\[\begin{aligned} \frac{\partial}{\partial \mu_j}\ln P(x^{(i)}\mid y=0;\mu,\Sigma)&=\frac{\partial}{\partial \mu_j}\left(C-\frac{1}{2}\ln|\Sigma|-\frac{1}{2}(x^{(i)}-\mu)\Sigma^{-1}(x^{(i)}-\mu)^T\right)\\ &=x^{(i)}_j-\mu_j \end{aligned} \]

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

\[\nabla_{\mu}\mathcal{L}=\sum_{i=1}^m(1-y^{(i)})(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}\)。

\[\begin{aligned} \nabla_{\Sigma}\mathcal{L}&=\sum_{i=1}^m-\frac{1}{2|\Sigma|}\cdot\nabla_{\Sigma}|\Sigma|+\frac{1}{2}\nabla_{\Sigma}\left((x^{(i)}-\mu)\Sigma^{-1}(x^{(i)}-\mu)^{T}\right)\\ &=\sum_{i=1}^m-\frac{1}{2}\Sigma^{-1}+\frac{1}{2}(x^{(i)}-\mu)^T(x^{(i)}-\mu)\Sigma^{-2} \end{aligned} \]

同样地令上式为 \(0\)，计算得

\[\Sigma=\frac{1}{m}\sum_{i=1}^m(x^{(i)}-\mu)^T(x^{(i)}-\mu) \]

注意到，按照我们的假设，\(\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;\mu,\Sigma)=\frac{1}{(2\pi)^{\frac{n}{2}}|\Sigma|^{\frac{1}{2}}}\cdot\exp\left(-\frac{(x-\mu)\Sigma^{-1}(x-\mu)^T}{2}\right) \]

但是实际上比较 \(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