【机器学习实验】概率编程及贝叶斯方法

引言

贝叶斯方法是天生用来做推断的方法，然而它常隐藏在课本的数学分析的背后。

随着近年来贝叶斯方法在机器学习竞赛中成功应用，其重要性又引起了学习者的兴趣，但是其难点在于贝叶斯数学和概率编程之间的衔接。《Probabilistic
Programming and Bayesian Methods for Hackers》一书试图弥补以上的遗憾。

有关概率编程和贝叶斯方法的实验，我将以该在线书籍作为学习资料，希望从中能将理论知识和编程实践想结合，以加深对理论的理解，以提高解决具体问题和数学建模的能力。

我们期望使用贝叶斯推断来进行统计分析，利用计算机的计算能力通过概率编程来解决问题。该方法去除了原本数学分析的干预。PyMC是一个马尔可夫链蒙特卡罗采样(Markov chain Monte Carlo Sampling)Python工具包，它实现了Metropolis-Hastings算法，包含了画图、拟合优度(goodness-of-fit)和聚合诊断(convergence diagnostics)的方法。

贝叶斯思想

比如测试一个算法的正确性，我们通过一个个测试用例来验证算法的正确性，每次测试都让我们对算法的准确性更加有信心，但是我们无法保证算法是没有bug的，除非测试了所有可能的问题，但这常常是不实际的。

贝叶斯推断是同样的一个道理，通过成功的结果来更新信任度，在没有验证所有可能的结果之前，我们无法绝对肯定。

现在，我们来对比一下贝叶斯和频率派对于概率的解释。

从贝叶斯的世界观中，概率是一个事件发生的可信度(believability)的度量，也就是说，我们对一件事情发生是多么有信心。

而频率派的观点，他们设定概率是事件的长期频率。比如飞机发生事故的概率，频率派解释为飞机事故的长期频率。这在很多事件的概率问题上看上去确实有逻辑合理性，但当一个事件不会经常出现，无法统计长期的发生频率的话，就很难解释通了。比如总统选举的概率，因为该选举只发生一次。

贝叶斯观点认为，概率是事件发生的可信度或者信心的测度。对于一个事件的发生，如果将这种可信度设为0，那么说明我们对该事件发生没有信心；相反，如果将这种可信度设为1，那么说明我们十分肯定该事件一定会发生。0到1之间的可信度是对结果发生的权重的考虑。这样的定义与飞机发生事故的概率的解释是一致的：通过观测了飞机事故的发生频率，在没有附加信息的情况下，飞机事故的可信度应该等于该频率。

我们设定可信度的测量是针对一个个体的，而不是对于自然而言的。这是因为，不同的个体具有不同的信息，所以其对于事件发生的信任程度也不同。

将可信度认为是一种概率的哲学观点对人类是很自然的，我们持续不断地与世界万物互动接触，总是看到一部分真相，却收集了一些证据来构成这种可信度。

贝叶斯派通过观察到更多的证据和现象来更正自己的信任度，尤其是当证据与初始认为的一致时，该证据就更不能被忽视了。我们将更新的可信度设为P(A|X)，解释为在给定证据X时，事件A发生的概率。这被称作后验概率。

先验概率P(A)：硬币有50%的机会向上；后验概率P(A|X)：你看到了硬币落到地上时正面朝上，这就是信息X，于是我们就将硬币朝上的概率设为1.0，硬币朝下的概率设为0.0。
先验概率P(A)：大型复杂的程序代码可能bug；后验概率P(A|X)：如果代码通过了所有的测试用例X，尽管仍然有bug存在的可能性，但是该可能性已经减小了。

贝叶斯推断实战

如果是频率派和贝叶斯派编写一个函数，输入一个统计问题，而返回的结果是不同的。频率派的推断方程会返回一个数来表示一个估计，而贝叶斯派的方程将返回一个概率。

以之前代码调试问题为例，当询问“代码通过了所有测试，那么该代码是不是没有bug了”，频率派的方程会返回“是”；而贝叶斯派的方程将返回“是”和“否”的概率。

令N为我们手上数据证据的数量，当我们收集了无限多个证据，即N→∞，贝叶斯的结果和频率派的结果是一致的。所以，对于很大的N，统计推断或多或少有一定客观性；而对于较小的N，该推断就不太稳定了：频率派估计就有较大的方差和较大的置信区间(confidence intervals)。这却是贝叶斯分析擅长的，贝叶斯派引入了先验概率，返回一个概率，这样就保有了对于小数据集情况下统计推断所反映出来的不确定性(uncertainty)。

有人会认为，当N很大的时候，人们会不太在意这两种技术的区别，因为此时它们提供了类似的推断，甚至更倾向于计算简单的频率派方法。

而有专家[Andrew Gelman (2005)]认为，样本数量永远不会很大。如果N很小，为了得到一个足够精确的推论，你需要获得更多的数据。但是一旦N足够大，你便开始细分数据以学习更多（比如，在一个公开的民意投票中，一旦你对于全国的情况有一个很好的估计，你便会划分不同的人群来进行更多的估计，比如按照性别，不同的年龄，不同的地域等）。N不会很充足，因为一旦你认为有充足的数据，你便转向更加细致的问题，这样你便需要更多的数据。

频率派方法是错误的吗？

频率派方法依然很有用，并且在很多领域都是先进水平(state-of-the-art)。像最小二乘线性回归、LASSO回归、期望最大算法都很有用、很快速。贝叶斯方法补充了这些技术的不足，或者用弹性建模(flexible modeling)来阐述潜在的系统。

关于大数据

反常的是，大数据的预测分析问题常
4000
常使用相对简单的算法，因此，我们辩解大数据预测的困难点不在于算法的使用，而是存储和使用大数据的难度。

更多的难以分析的问题包含中数据(medium data)以至于小数据(small data)。

贝叶斯框架

入门实验

抛硬币实验

假设你不知道抛硬币正面朝上的的概率，你定义该真实比率为p，但是并没有关于p的任何先验知识。 于是，我们尝试抛硬币并记录其观测结果。有趣的问题来了，随着我们观察到越来越多的数据，我们的推断是如何发生变化的？

%matplotlib inline
from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
figsize(11, 9)

import scipy.stats as stats

dist = stats.beta
#多次试验
n_trials = [0, 1, 2, 3, 4, 5, 8, 15, 50, 500]
#构造500个符合伯努利分布的随机抽样
data = stats.bernoulli.rvs(0.5, size=n_trials[-1])
x = np.linspace(0, 1, 100)

#使用伯努利分布的共轭先验——Beta分布
for k, N in enumerate(n_trials):
sx = plt.subplot(len(n_trials) / 2, 2, k + 1)
plt.xlabel("$p$, probability of heads") \
if k in [0, len(n_trials) - 1] else None
plt.setp(sx.get_yticklabels(), visible=False)#不显示y轴刻度
heads = data[:N].sum()#统计1的个数
y = dist.pdf(x, 1 + heads, 1 + N - heads)#构造Beta分布
plt.plot(x, y, label="observe %d tosses,\n %d heads" % (N, heads))
plt.fill_between(x, 0, y, color="#FF99CC", alpha=0.4)#填充色彩到曲线下区域
plt.vlines(0.5, 0, 4, color="k", linestyles="--", lw=1)#在0.5处画虚线

leg = plt.legend()
leg.get_frame().set_alpha(0.4)
plt.autoscale(tight=True)

plt.suptitle("Bayesian updating of posterior probabilities",
y=1.02,
fontsize=14)

plt.tight_layout()

上图，我们随着观测抛硬币的数据的增多，绘制了一系列后验概率的更新情况。
该后验概率由这个曲线来表征，其不确定程度正比于曲线的宽度。我们开始观察数据的过程中，后验概率开始平移，不断变化。最终，该概率变得紧缩且越来越接近其真实值p=0.5。

但注意到，最终图像的峰值不一定在0.5处，这是很正常的。回顾之前，我们假设我们对p的值没有一个先验知识。实际上，如果我们观测了8次的结果，其中只有1次硬币朝上，那么得到的分布情况与峰值在0.5周围的情况相比偏倚很大。随着更多数据的积累，我们会发现该概率逐渐被分配在p=0.5处。

二项分布与Beta分布

二项分布

在概率论和统计学中，二项分布是n个独立的[是/非]试验中成功的次数的离散概率分布，其中每次试验的成功概率为p。

比如：
抛一次硬币出现正面的概率是0.5，抛10次硬币，出现k次正面的概率。
掷一次骰子出现六点的概率是1/6，掷 6次骰子，出现k次六点的概率。

每次抛硬币或者掷骰子都和上次的结果无关，所以每次实验都是独立的。二项分布是一个离散分布，k的取值范围为从0到n，只有n+1种可能的结果。

scipy.stats.binom为二项分布，下面用它计算抛十次硬币，出现k次正面的概率分布。

n = 10
k = np.arange(n+1)
pcoin = stats.binom.pmf(k, n, 0.5)#k -> quantiles, n and p -> shape parameters
plt.stem(k, pcoin)#绘制火柴棍图
plt.margins(0.1,0.05)#设置自动缩放的边缘

下面是投掷6次骰子，出现6点的概率分布。

n = 6
k = np.arange(n+1)
pdice = stats.binom.pmf(k, n, 1.0/6)
plt.stem(k, pdice)
plt.margins(0.1)

Beta分布

贝塔分布是一个作为伯努利分布和二项式分布的共轭先验分布的密度函数，在机器学习和数理统计学中有重要应用。贝塔分布中的参数可以理解为伪计数，伯努利分布的似然函数可以表示为，表示一次事件发生的概率，它为贝塔有相同的形式，因此可以用贝塔分布作为其先验分布。

对于硬币或者骰子这样的简单实验，我们事先能很准确地掌握系统成功的概率。然而通常情况下，系统成功的概率是未知的。为了测试系统的成功概率p，我们做n次试验，统计成功的次数k，于是很直观地就可以计算出p = k/n。然而由于系统成功的概率是未知的，这个公式计算出的p只是系统成功概率的最佳估计。也就是说实际上p也可能为其它的值，只是为其它的值的概率较小。

例如有某种特殊的硬币，我们事先完全无法确定它出现正面的概率。然后抛10次硬币，出现5次正面，于是我们认为硬币出现正面的概率最可能是0.5。但是即使硬币出现正面的概率为0.4，也会出现抛10次出现5次正面的情况。因此我们并不能完全确定硬币出现正面的概率就是0.5，所以p也是一个随机变量，它符合Beta分布。

Beta分布是一个连续分布，由于它描述概率$p$的分布，因此其取值范围为0到1。 Beta分布有alpha和beta两个参数，其中alpha为成功次数加1，beta为失败次数加1。

连续分布用概率密度函数描述，下面绘制实验10次，成功4次和5次时，系统成功概率p的分布情况。可以看出k=5时，曲线的峰值在p=0.5处，而k=4时，曲线的峰值在p=0.4处。

当n=30，k=12时，曲线的峰值还在p=0.4处，但是因为随着实验次数的增多，p取其他值的可能就会变小，对p的估计就更有信心，因此山峰就更陡峭了。

n = 10

k = 5

p = np.linspace(0, 1, 100)

pbeta = stats.beta.pdf(p, k+1, n-k+1)

plt.plot(p, pbeta, label="k=5", lw=2)

k = 4

pbeta = stats.beta.pdf(p, k+1, n-k+1)

plt.plot(p, pbeta, label="k=4", lw=2)

n = 30

k = 12

pbeta = stats.beta.pdf(p, k+1, n-k+1)

plt.plot(p, pbeta, label="k=3", lw=2)

plt.xlabel("$p$")

plt.legend(loc="best")

小结

这一小节，我们先对贝叶斯概念有个初步的了解，在下一小节中，我将进一步地介绍不同的数学分布，并引入pymc库的使用。

转载请注明作者Jason Ding及其出处
Github博客主页(http://jasonding1354.github.io/)
CSDN博客(http://blog.csdn.net/jasonding1354)
简书主页(http://www.jianshu.com/users/2bd9b48f6ea8/latest_articles)
百度搜索jasonding1354进入我的博客主页