高斯判别分析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\)（存疑，不知道目的）。
最大似然估计
在判别学习中，我们以 \(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}\)。
同样地令上式为 \(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 函数有修改）
def Gaussian(self, x, mean, cov):
"""
这是自定义的高斯分布概率密度函数
:param x: 输入数据
:param mean: 均值向量
:param cov: 协方差矩阵
:return: x的概率
"""
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(covinv).dot(xdiff.T))[0][0]
return prob
最后是实现的一个类：
import numpy as np
class GaussianDiscrimitiveAnalysis:
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 Gaussian(self, x, mean, cov):
"""
这是自定义的高斯分布概率密度函数
:param x: 输入数据
:param mean: 均值向量
:param cov: 协方差矩阵
:return: x的概率
"""
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(covinv).dot(xdiff.T))[0][0]
return prob
def predict(self, xs):
sigma_i = np.linalg.inv(self.sigma)
predict = []
for x in xs:
evaluate_0 = self.Gaussian(x, self.mu[0], self.sigma)
evaluate_1 = self.Gaussian(x, self.mu[1], self.sigma)
if evaluate_0 > evaluate_1:
predict.append(0)
else:
predict.append(1)
return predict