高斯判别分析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

最新文章

  1. GridView点击行触发SelectedIndexChanged事件
  2. (一)NOR FALSH 讲解
  3. statusbar 样式
  4. jquery隐藏table表格的某一列
  5. IIS-Server is too busy _解决方法
  6. HDU 5898:odd-even number(数位DP)
  7. emoji表情符处理替换成空格
  8. STL容器与配接器
  9. The sum - SGU 122(斐波那契前N项和)
  10. aop切入点表达式
  11. String[255]在高版本Delphi里还是被解释成Byte,总体长度256,使用StrPCopy可以给Array String拷贝字符串(内含许多实验测试)
  12. java-成员方法/变量、类方法/变量等区别
  13. seajs模块标识命名和解析规则
  14. 在asp.net core2.1中添加中间件以扩展Swashbuckle.AspNetCore3.0支持简单的文档访问权限控制
  15. Spark实战
  16. Lily_music 网页音乐播放器 -可搜索(附歌词联动播放效果解说)
  17. [C++]Qt 如何处理密集型耗时的事情(频繁调用QApplication::processEvents)
  18. Python3 下实现 腾讯人工智能API 调用
  19. php手撸轻量级开发(一)
  20. mysql INSERT ... ON DUPLICATE KEY UPDATE语句

热门文章

  1. CSS 3 所有的选择器整理(2023.2)
  2. Kubernetes(k8s)控制器(四):ReplicaSet
  3. Javaweb-Tomcat(安装+配置环境)
  4. ASP.NET Core开发者指南(2022版路线图)
  5. 微信小程序分类菜单激活状态跟随列表滚动自动切换
  6. 如何优化线上WebAssembly
  7. 深入理解JavaScript对象
  8. C# 通过反射获取类字段名和值并加入到字典中(包含递归获取)
  9. sqlserver生成指定的序列
  10. ThreadLocal最终版本