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

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

生成学习

处理分类问题，我们可以使用逻辑回归、Softmax。这两种方法都属于“判别学习”，也就是给定 $(x^{(i)},y^{(i)})$，我们学习 $P(y|x)$，并对于给定的 $x$，计算 ${\mathrm{argmax}}_y\{P(y|x)\}$。

GDA属于另一种方法——生成学习。在判别学习中，我们并不关注 $x$ 本身的分布，而在生成学习中，我们基于一些事实假设 $x$ 的分布，例如在GDA中，我们假设对于相同的 $y$，$x$ 符合高斯分布。然后基于假设，学习 $P(x|y)$ 的参数（高斯分布则学习均值 $\mu$ 和协方差 $\mathrm{\Sigma }$），以及 $P(y)$ 的参数（例如两点分布的 $\varphi$），从而得到联合概率 $P(x,y)$。在预测时，我们仍然是找到最大的 $P(y|x)$，但是是用贝叶斯公式计算：

$$P(y|x)=\frac{P(x,y)}{P(x)}$$

高斯判别分析

作出的假设

首先假设 $Y$ 服从两点分布 $Bernolli(\varphi )$，于是可以写作 $P(y)={\varphi }^y(1-\varphi {)}^{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 (-\frac{(x_i-{\mu }_i{)}^2}{2{\sigma }_i^2})$$

又 $x_i$ 相互独立，可以得到：

$$\begin{array}{rl}P(x|y=0)& =\prod _{i=1}^n\frac{1}{\sqrt{2\pi }{\sigma }_i}\exp (-\frac{(x_i-{\mu }_i{)}^2}{2{\sigma }_i^2})\\ & =\frac{1}{(2\pi {)}^{\frac{n}{2}}{\sigma }_1\dots {\sigma }_n}\exp (-\frac{1}{2}\sum _{i=1}^n\frac{(x_i-{\mu }_i{)}^2}{{\sigma }_i^2})\end{array}$$

我们定义矩阵 $\mathrm{\Sigma }=\left(\begin{array}{c}{\sigma }_1^2\\ & \ddots \\ & & {\sigma }_n^2\end{array}\right)$，可以将上式化简：

$$P(x|y=0)=\frac{1}{(2\pi {)}^{\frac{n}{2}}|\mathrm{\Sigma }{|}^{\frac{1}{2}}}\exp (-\frac{1}{2}(x-\mu ){\mathrm{\Sigma }}^{-1}(x-\mu {)}^T)$$

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

但是实际上GDA并没有假设 $\mathrm{\Sigma }$

最大似然估计

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

$$\begin{array}{rl}\mathcal{L}(\varphi ,\mu ,{\mu }^{\prime },\mathrm{\Sigma })& =\ln \prod _{i=1}^{m}P(x^{(i)},y^{(i)};\varphi ,\mu ,{\mu }^{\prime },\mathrm{\Sigma })\\ & =\ln \prod _{i=1}^{m}P(x^{(i)}\mid y^{(i)};\mu ,{\mu }^{\prime },\mathrm{\Sigma })P(y^{(i)};\varphi )\\ & =\sum _{i=1}^m\ln P(x^{(i)}\mid y^{(i)};\mu ,{\mu }^{\prime },\mathrm{\Sigma })+\sum _{i=1}^m\ln P(y^{(i)};\varphi )\end{array}$$

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

计算 $\varphi$

$$\begin{array}{rl}\frac{\partial \mathcal{L}}{\partial \varphi }& =\frac{\partial }{\partial \varphi }\sum _{i=1}^m\ln P(y^{(i)};\varphi )\\ & =\frac{\partial }{\partial \varphi }\sum _{i=1}^m{y}^{(i)}\ln \varphi +(1-y^{(i)})\ln (1-\varphi )\\ & =\sum _{i=1}^m\frac{y^{(i)}}{\varphi }-\frac{1-y^{(i)}}{1-\varphi }\end{array}$$

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

计算 $\mu ,{\mu }^{\prime }$

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

$$\begin{array}{rl}{\mathrm{\nabla }}_{\mu }\mathcal{L}& ={\mathrm{\nabla }}_{\mu }\sum _{i=1}^m(1-y^{(i)})\ln P(x^{(i)}\mid y=0;\mu ,\mathrm{\Sigma })+y^{(i)}\ln P(x^{(i)}\mid y=1;{\mu }^{\prime },\mathrm{\Sigma })\end{array}$$

其中

$$\begin{array}{rl}\frac{\partial }{\partial {\mu }_j}\ln P(x^{(i)}\mid y=0;\mu ,\mathrm{\Sigma })& =\frac{\partial }{\partial {\mu }_j}(C-\frac{1}{2}\ln |\mathrm{\Sigma }|-\frac{1}{2}(x^{(i)}-\mu ){\mathrm{\Sigma }}^{-1}(x^{(i)}-\mu {)}^T)\\ & =x_j^{(i)}-{\mu }_j\end{array}$$

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

$${\mathrm{\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 }^{\prime }=\frac{\sum y^{(i)}{x}^{(i)}}{\sum y^{(i)}}$。

计算 $\mathrm{\Sigma }$

下面用到了两个 $\mathrm{\nabla }$ 算符的性质：

• ${\mathrm{\nabla }}_{\mathrm{\Sigma }}|\mathrm{\Sigma }|=|\mathrm{\Sigma }|{\mathrm{\Sigma }}^{-1}$；
• ${\mathrm{\nabla }}_{\mathrm{\Sigma }}{\mathrm{\Sigma }}^{-1}=-{\mathrm{\Sigma }}^{-2}$。

$$\begin{array}{rl}{\mathrm{\nabla }}_{\mathrm{\Sigma }}\mathcal{L}& =\sum _{i=1}^m-\frac{1}{2|\mathrm{\Sigma }|}\cdot {\mathrm{\nabla }}_{\mathrm{\Sigma }}|\mathrm{\Sigma }|+\frac{1}{2}{\mathrm{\nabla }}_{\mathrm{\Sigma }}((x^{(i)}-\mu ){\mathrm{\Sigma }}^{-1}(x^{(i)}-\mu {)}^T)\\ & =\sum _{i=1}^m-\frac{1}{2}{\mathrm{\Sigma }}^{-1}+\frac{1}{2}(x^{(i)}-\mu {)}^T(x^{(i)}-\mu ){\mathrm{\Sigma }}^{-2}\end{array}$$

同样地令上式为 $0$，计算得

$$\mathrm{\Sigma }=\frac{1}{m}\sum _{i=1}^m(x^{(i)}-\mu {)}^T(x^{(i)}-\mu )$$

注意到，按照我们的假设，$\mathrm{\Sigma }$ 应该是一个对角矩阵（如果了解协方差矩阵的话，由 $x_i$ 相互独立，可以推出 $\mathrm{\Sigma }$ 应该是对角矩阵，对角元就是每个变量的方差），但是这里非常显然不总是对角矩阵。

（以下仅是我个人的猜测）在GDA的实现中，我们并没有关注 $x_i$ 相互独立的性质，而是直接学习了一个普遍的协方差矩阵。实际上这种定义的限制更低，更符合现实情况（现实中的变量之间存在联系比较普遍）。

代码实现

直接计算这几个参数的代码不需要解释：

Copym, 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

最后我们发现计算概率需要计算 ${\mathrm{\Sigma }}^{-1}$，然而我们算出来的可能是一个奇异矩阵。这时候可以对 $\mathrm{\Sigma }$ 进行微扰——将对角元加上一个较小的偏差。

最后我们对 $y=0,1$ 分别计算出 $P(x|y)$。

$$P(x|y;\mu ,\mathrm{\Sigma })=\frac{1}{(2\pi {)}^{\frac{n}{2}}|\mathrm{\Sigma }{|}^{\frac{1}{2}}}\cdot \exp (-\frac{(x-\mu ){\mathrm{\Sigma }}^{-1}(x-\mu {)}^T}{2})$$

但是实际上比较 $P(x|y=0),P(x|y=1)$ 不需要真正计算出概率，而是比较两者不同的地方，也就是后面的 $\exp$。

（奇异矩阵的处理参考 https://blog.csdn.net/qq_30091945/article/details/81508055 ，对其中的 Gaussian 函数有修改）

最后是实现的一个类：

Copyclass 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