本周机器学习小项目的主要目标在于引入概率模型。之前内容使用函数来描述模型，本周则侧重于从概率角度描述机器学习。概率与统计几乎是整个机器学习最重要的基础，甚至可能是深度学习之后新一代机器学习模型的基础。
本周主要内容包括：频率学派与贝叶斯学派、最大似然估计（MLE）、最大后验估计（MAP）、交叉熵作为损失函数与其反向传播过程实现、一个高斯环境下的贝叶斯分类器实现。最后就神经网络与贝叶斯分类器的等价关系作一说明与对比。
本周推荐阅读时间：25min
一、阅读建议
- 推荐阅读时间：25min
- 推荐阅读文章：李航. 统计学习方法[M]. 清华大学出版社, 2012.
- 本章涉及主要内容：
- 最大似然估计
- 最大后验估计
- 新的损失函数：交叉熵
- 高斯环境下的贝叶斯分类器
二、软件环境
- Python3
- Numpy
- Matplotlib
三、前置基础
- 概率与统计
四、数据描述
本章数据分为两个部分，第一部分依然来源于 Yann LeCun 文章，第二部分使用 Numpy 生成的训练样本。
4.1 数据1
- 数据来源：Yann LeCun 文章
- 数据下载：http://yann.lecun.com/exdb/mnist/
- 数据描述：MNIST（Mixed National Institute of Standards and Technology database）是一个计算机视觉数据集，它包含 70000 张手写数字的灰度图片，其中每一张图片包含 28 X 28 个像素点。保存形式为长度为 784 的向量。
4.2 数据2
- 数据来源：Numpy 产生数据
- 生成目的：用于测试高斯环境下的贝叶斯分类器算法
- 数据描述：数据有两列属性，数据总共两类，用于分类算法验证
- 生成代码：参考 7.1 节
5. 理论部分
5.1 最大似然估计与最大后验估计
举个例子来：讲抛硬币实验。记录 n 次随机试验 $\Omega=\{\omega_1,\cdots,\omega_n\} $ ，出现正面的次数为 m，假设正面出现的概率为 $\theta=p(正面)$，那么产生随机试验结果 $\Omega$ 的概率为：
$$p(\Omega|\theta)$$(1.1)
1.1式中 $\theta$ 为自变量，频率学派认为这个概率 $\theta$ 最优解应该使得 1.1 取最大值。在假设了正面出现的概率后，同时假设样本产生是有顺序，而且是独立同分布的，那么可以计算出现实验结果 $\Omega$ 概率的具体形式：
$$p(\Omega|\theta)=\theta^m (1-\theta)^{n-m}$$(1.2)
求解上述以 $\theta$ 为自变量的函数的最大值:
1.3式中为了方便计算将概率修改为了对数形式：
求其最小值，对于抛硬币问题来讲：
$$f(\theta)=m \log(\theta)+(n-m)\log(1-\theta)$$(1.5)
1.5取得最小值时 $\theta=\frac{m}{n}$，所以频率学派认为取得正面的概率为 $m/n$，这个结果可以通过频率统计得到，而求解硬币取得正面概率的过程称为最大似然估计（MLE）。
贝叶斯学派则认为概率 $\theta$ 本身符合一个先验分布。这个先验分布可以纠正采样的偏差：
$$p(\theta|m)=\frac{p(m|\theta)p(\theta)}{p(m)}=\frac{p(m|\theta)p(\theta)}{\int{p(m|\phi)d\phi}}$$(1.6)
1.6式中
- $p(\theta|m)$ 为后验
- $p(m|\theta)$ 为似然
- $p(\theta)$ 为先验
- 后验概率 = 先验概率 + 数据信息
先验概率是什么意思，假设为了估计此次投掷硬币为正面的概率，我先用的自己的硬币做了 n=4 次实验，正面出现了 m=2 次，假设此时正面概率为 x，那么出现这种情况的概率为：
$$\begin{matrix}p(x|n, m)=\frac{x^m (1-x)^{n-m}}{constant}\=\frac{\Gamma(n+2)}{\Gamma(m+1)\Gamma(n+1)}x^m (1-x)^{n-m}\=\frac{1}{B(m+1, n-m+1)}x^m(1-x)^{n-m}=30 x^{2}(1-x)^2\end{matrix}$$(1.7)
1.7 中 constant 为归一化常数，称为 Beta 函数：
$$B(m,n)=\int_0^1 x^{m-1}(1-x)^{n-1}dx$$(1.8)
到此为止，我们都在描述 x 的分布，也就是硬币为正面概率的分布，可以看到使用我自己的硬币估计时，取得 0.5 的概率是最大的：
这个概率的分布可以作为我们的先验分布，它是概率的概率。以这个先验分布，我们去估计其他硬币的实验，假设使用另一枚硬币，抛 4 次，出现正面为 p=1 次，出现反面为 q=3 次，那么在贝叶斯理论下出现正面的概率为：
$$\begin{matrix}p(\theta|p, q)=\frac{p(p, q|\theta)p(\theta)}{\int{p(p, q|\phi)d\phi}}\=\frac{\theta^{p+2}(1-\theta)^{q+2}}{\int{\phi^{p+2}(1-\phi)^{q+2}d\phi}}=\frac{\theta^{3}(1-\theta)^{5}}{\int{\phi^{3}(1-\phi)^{5}d\phi}}\end{matrix}$$(1.9)
将 1.9式绘制成图：
此时硬币概率最大值为 0.375，这里的数值很重要，我们以最大似然估计估计的硬币概率是 0.25，此时是没有先验分布的，如果引入了先验分布，取得正面概率应该为 0.375。而我们抛硬币的次数仅为 4 次，因此很大可能出现偏差，而先验概率的引入则可以纠正这种偏差。从另一个角度来看相当于在最大似然估计的基础上加入了正则化项，这称为最大后验估计（MAP）
名词：独立同分布 = Independent and identical distribution=i.i.d.
5.2 信息论与交叉熵
首先来看熵：
$$H(p)=\mathbb{E}_p(\log(1/p(x)))=-\int p(x)\log(p(x))dx$$(1.10)
这里其实没有必要对于熵进行更具体的比喻，因为数学定义本身就很直观。再来看交叉熵：
交叉熵用来衡量两个分布 p, q 的相似度，分布相似度越高 H(p, q) 越小。
5.3 概率模型下的神经网络
分类问题神经网络预测过程可以从函数描述转换为概率描述：
$$y=f(x)\rightarrow p(y|x,w)$$(1.12)
1.12式子描述了：给定样本 x 后神经网络输出属于某一类的概率 y，网络的可训练参数为 w。下面的问题就是如何将神经网络的输出转换为概率。假设神经网络输出为 h，因为概率均是大于 0 的，说以需要将 h 转换为大于 0 的数值，同时所有情况的概率之和为 1，整个的处理过程称之为 softmax：
而神经网络的训练过程可以描述为：
1.14中，d 为 one-hot 形式的向量，其代表样本属于某一类概率的标签。比如对于手写数字来说，总共有十个数字，那么 one-hot 向量长度为 10，如果手写数字为 6，那么 6 的位置上数字为 1，代表概率为 100%，其他位置为 0。可以看到最大似然估计与交叉熵作为损失函数具有相同的意义。
5.4 交叉熵损函数的反向传播过程
假设标签为 d，神经网络输出为 h，那么 loss 函数定义为：
反向传播过程：
5.5 最大似然估计下的贝叶斯分类器
假设数据分布是未知的，我们可以给出数据的近似分布形式，近似分布形式有可训练参数。这些可训练参数可以通过最大似然估计得到。再来复述一下这个过程：假定数据分布形式 $p(x;w)$-> 通过最大似然估计求解 $w=\hat{w}$-> 计算 x 的概率 $p(x;\hat{w})$。
首先我们需要假设数据的分布形式：$$p(x)=\frac{1}{\sqrt{2\pi}\sigma}\exp(-\frac{(x-\mu)^2}{2\sigma^2})$$(1.17)
这里在假设分布的时候用到了两个可训练参数，$\sigma,\mu$。使用最大似然估计对可训练参数进行求解：
使得 1.18 导数为 0：
1.18 假定的分布为高斯分布，1.19 为根据最大似然估计所求的可训练参数。可以看到两个可训练参数分别就是统计所得到的均值和方差。贝叶斯分类器实际上就是：
1.20 中 $p(y_i)$ 为样本属于第 $y_i$ 类概率，x 为输入数据。模型为给定样本估计其属于某一类的概率。$p(x_j|y_i)$ 可以假设为高斯分布，那么其求解形式为最大似然估计 1.19 式。而属于某一类的概率 $p(y_i)$ 为先验，可以通过训练数据统计得到。
5.6 贝叶斯分类器与神经网络
实际上，高斯环境下的贝叶斯分类器等同于单层神经网络的，只是需要加入数据的二次方项。具体内容请参考《神经网络与机器学习》，或者文章朴素贝叶斯与单层神经网络关系
6. 代码部分
6.1 交叉熵损失函数
def _loss_softmax_cross_entropy(self, Y, *args, **kw): return np.square(Y) def _d_loss_softmax_cross_entropy(self, Y, *args, **kw): B = np.shape(Y)[0] X = self.outputs[-2] eX = np.exp(X) norm = eX / np.sum(eX, axis=1, keepdims=True) ce = - Y + norm return ce def loss_softmax_cross_entropy(self): self.value.append([]) self.d_value.append([]) self.layer.append((self._loss_softmax_cross_entropy, None, self._d_loss_softmax_cross_entropy, None)) self.layer_name.append("loss-softmax-cross-entropy") 6.2 高斯环境下的贝叶斯分类器
import numpy as npclass GaussianNaiveBayes(): def N(self, x, mu, std): """ 标准正态分布 """ par = 1/(np.sqrt(2*np.pi)*std) return par*np.exp(-(x-mu)**2/2/std**2) def logN(self, x, class_type): """ 标准正态分布对数 """ if class_type==0: return np.log(self.N(x, self.mu1, self.std1)) else: return np.log(self.N(x, self.mu2, self.std2)) def fit(self, X, y): """ 训练过程为最大似然估计参数 """ X1 = X[y==0] X2 = X[y==1] self.mu1 = np.mean(X1, axis=0) self.mu2 = np.mean(X2, axis=0) self.std1 = np.std(X1, axis=0) self.std2 = np.std(X2, axis=0) self.pre1 = np.log(len(X1)/len(X)) self.pre2 = np.log(len(X2)/len(X)) def predict_proba(self, xx): """ 预测过程 """ prb = [] for x in xx: prb1_log = np.sum(self.logN(x, 0)) + self.pre1 prb2_log = np.sum(self.logN(x, 1)) + self.pre2 prb1 = np.e ** prb1_log prb2 = np.e ** prb2_log prb1 = prb1 / (prb1 + prb2) prb2 = prb2 / (prb1 + prb2) prb.append([prb1, prb2]) return np.array(prb)7. 程序运行
7.1 数据产生
# 引入绘图库import matplotlib.pyplot as pltimport matplotlib# 设置风格和中文plt.rcParams['font.sans-serif']=['SimHei'] plt.rcParams['axes.unicode_minus']=Falsematplotlib.style.use("ggplot")# 训练数据产生theta1 = np.random.random([1000]) * np.pi * 2theta2 = np.random.random([1000]) * np.pi * 2r1 = np.random.normal(loc=1, scale=0.2, size=[1000])r2 = np.random.normal(loc=0, scale=0.3, size=[1000])X = np.zeros([2000, 2])y = np.zeros([2000])X[:1000, 0] = np.cos(theta1) * r1X[:1000, 1] = np.sin(theta1) * r1X[1000:, 0] = np.cos(theta2) * r2X[1000:, 1] = np.sin(theta2) * r2y[1000:] = 1# 将标签转换为one-hot编码y_one_hot = np.zeros([2000, 2])y_one_hot[:1000, 0] = 1y_one_hot[1000:, 1] = 1# 训练输入数据添加二次方项X2 = np.concatenate([X, X[:, 0:1] * X[:, 0:1], X[:, 1:2] * X[:, 1:2], X[:, 0:1] * X[:, 1:2]], axis=1)# 绘图plt.scatter(X[:, 0], X[:, 1], c=y)plt.title("训练数据散点图")plt.show()产生训练图片
7.2 神经网络
定义单层神经网络，但是在数据预处理过程之中加入数据的二次方项：
w = np.random.uniform(-0.1, 0.1, size=[5, 2])b = np.zeros([2])NeuralNet = NN()NeuralNet.matmul(w)NeuralNet.bias_add(b)NeuralNet.loss_softmax_cross_entropy()# 数据加入二次方项X2 = np.concatenate([X, X[:, 0:1] * X[:, 0:1], X[:, 1:2] * X[:, 1:2], X[:, 0:1] * X[:, 1:2]], axis=1)for itr in range(200): idx = np.random.randint(0, 2000, size=[32]) NeuralNet.fit(X2[idx], y_one_hot[idx])7.3 朴素贝叶斯
朴素贝叶斯算法直接输入原始数据
NaiveBayes = GaussianNaiveBayes()NaiveBayes.fit(X, y)7.4 绘制分类边界并对比
x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.1), np.arange(y_min, y_max, 0.1))x_test = np.c_[xx.ravel(), yy.ravel()]x2_test = np.concatenate([x_test, x_test[:, 0:1] * x_test[:, 0:1], x_test[:, 1:2] * x_test[:, 1:2], x_test[:, 0:1] * x_test[:, 1:2]], axis=1)y_neural_pred = np.argmax(NeuralNet.predict(x2_test), axis=1)y_bayes_pred = np.argmax(NaiveBayes.predict_proba(x_test), axis=1)plt.subplot(131)plt.scatter(X[:, 0], X[:, 1], c=y)plt.title("训练数据散点图")plt.subplot(132)plt.scatter(X[:, 0], X[:, 1], c=y)plt.contourf(xx, yy, y_neural_pred.reshape(xx.shape), alpha=.4, cmap=plt.get_cmap("spring"))plt.title("单层神经网络预测结果")plt.subplot(133)plt.scatter(X[:, 0], X[:, 1], c=y)plt.contourf(xx, yy, y_bayes_pred.reshape(xx.shape), alpha=.4, cmap=plt.get_cmap("spring"))plt.title("高斯环境下的朴素贝叶斯分类器结果")plt.show()最终获取结果：
接下来需要修改什么
- 时序依赖问题
- 训练过程迭代缓慢
- 其他梯度迭代方法比如 Adam 引入
本文首发于GitChat，未经授权不得转载，转载需与GitChat联系。
阅读全文: http://gitbook.cn/gitchat/activity/5afce08b109cd76e3c1dc794
您还可以下载 CSDN 旗下精品原创内容社区 GitChat App ，阅读更多 GitChat 专享技术内容哦。
扫一扫
1047
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
