谨慎选择我们的先验 - 对最优拟合概率分布搜索空间的一种约束

1. 哲学上先验的概念

先验，是康德哲学中的重要概念，它并不是单纯在机器学习或者说AI中特定的专有名词，实际上非常多学科中都包含有先验（prior）和后验（posterior）的概念。

在学习贝叶斯推断统计的先验概念之前，我们这章先来了解一些简单的历史。

0x1：人类知性的组成

我们的知性有两方面的功能：

• 一是逻辑功能：逻辑功能在知识的基础上规定着思维的判断形式，这是形式逻辑的范围；
• 二是认识功能：认识功能则为我们提供新的知识，这是“先验逻辑”的领域；

自亚里士多德以来，我们的形式逻辑已经相当完备了，没有人怀疑形式逻辑是普遍必然的亦即“先天的”，因而为形式逻辑提供知识内容的先验逻辑也一定是先天的，因为它是形式逻辑的基础。

所以，形式逻辑与先验逻辑是对应的，在每一个判断形式的背后都有某种先验的要素作为它的基础，这种先验要素就是“范畴”。

0x2：康德”哥白尼式的革命“

康德之前，西方近代哲学主要分为理性主义和经验主义两派，康德哲学的出现初步调和了两派在认识论上的分歧。
康德解决的办法是，颠倒以往主观去符合客观的基本认识，认为是客观来符合主观，这就是康德自称的”哥白尼式的革命“。

1. 客观符合主观

所谓客观符合主观，并不是说对事物的认识完全由认识主体决定。客观符合主观，是说经验得到的东西要经过主观固有认识形式的”检验“，才能形成主体的认识。
这种固有认识形式不是哪个主体自己随便决定的，而是思维本身固有的结构。
从某种程度上来说，客观符合主观，这里的”主观“也是一种客观（固有、不变）。如康德本人所说，他不是主观唯心主义者。

2. 认识的来源

康德认为认识的来源包括经验，但又不全是经验，经验必须经过主观固有认识形式，即一张”知识之网“的过滤、判决后，才能形成认识。
这张知识之网是思维本身固有结构，这就保证了认识会具有普遍性。如此一来，经验主义和理性主义一定程度上被康德糅合了。

3. 先验的核心思想

先验是先天的、、普遍的、固有的、绝对的。先验先于经验，有涉及经验，是先天的先天，是知识的知识。
先验固有在思维本身之中，给予了先天、知识可能性的讨论，由此才有经验、知识。
先验的思想，表现了康德”哥白尼式的革命“的核心思想，即是有一种”先验“来讨论知识的可能性，那么从外物得到的经验就必须通过”先验“这一关，这才有了客观服从主观、对象服从主体，即”哥白尼式的革命“。

4. 先验调和了理性主义和经验主义

先验是既先于经验（先天），又与经验有关，这是康德哲学的一个重要概念。康德在一定程度上调和了经验主义和理性主义，其处理方法中最重要的一个就是”先验“的概念。
先于经验的”先验“，是作为经验的裁决者，作为知识的知识出现的，因为它考察的是”知识如何可能“。有了先验来考察知识如何可能，直观就可以对外物进行择别。
知识是离不开经验的，离不开外物的，这承认了经验主义的观点，另一方面，主观却要用先验去裁决，去作为知识结构之网、知识结构只固有框架。

笔者思考：从某种程度上讲，笔者认为先验思想也调和了频率统计推断和贝叶斯统计推断之间的鸿沟，让基于数据驱动的机器学习和基于统计进行推断的贝叶斯学习之间不在界限分明。先验是一个伟大的理论创新。

Relevant Link:

https://baike.baidu.com/item/%E5%85%88%E9%AA%8C%E5%88%86%E6%9E%90%E5%88%A4%E6%96%AD/7101736 
https://wenku.baidu.com/view/18cce387e53a580216fcfedd.html
《实践理性批判》康德

 

2. 主观和客观先验

贝叶斯先验可以分为两类：

• 客观先验：皆在让数据最大程度地影响后验；
• 主观先验：让领域专家来表达自己对先验的个人看法；

主观先验和客观先验是什么？如何精确定义某个先验是否是客观先验呢？这个小节我们来一起看几个例子，体悟一下主观和客观先验的涵义。

0x1：客观先验

1. 扁平先验

这是一种在整个未知参数范围内的均匀分布。用扁平先验意味着我们给每一个可能的值相等的权重，选择这种类型的先验我们称之为无差别原理，我们没有理由偏好某个具体数值。

扁平先验在机器学习范畴中还有另一个名词，叫最大熵原理，即当我们对未知参数的形式不确定时，应该将其初始化扁平均值形式。

在实际工程项目中，我们也一定常常有意无意地使用这个这种先验，例如你将未知参数初始化为相等的均值，或者初始化为0。

2. 杰弗里斯先验 - 建立一个不会因偶然改变变量位置而大幅变化的先验

上一个小节，我们谈到说扁平先验是一个客观先验。但他真的是绝对客观的吗？

在某种程度上，我们所说的客观是，一个不偏向后验估计的先验。扁平先验看起来是一个合理的选择，因为它对所有的参数赋予相同的概率。

但是扁平先验不是变换不变的。假设我们有一个来自于伯努利分布 θ 的随机变量 x，我们定义在 p(θ) = 1 上的先验，如下图所示：

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats


x = np.linspace(0.000, 1, 150)
y = np.linspace(1.0, 1.0, 150)
lines = plt.plot(x, y, color="#A60628", lw=3)
plt.fill_between(x, 0, y, alpha=0.2, color=lines[0].get_color())
plt.autoscale(tight=True)
plt.ylim(0, 2)

plt.show()

现在，让我们用函数来变换θ，这个变换函数的功能仅仅是在实轴上拉伸 θ。我们来看下，通过拉伸变换，先验函数会变成什么样呢？

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats


psi = np.linspace(-10, 10, 150)
y = np.exp(psi) / (1 + np.exp(psi)) ** 2
lines = plt.plot(psi, y, color="#A60628", lw=3)
plt.fill_between(psi, 0, y, alpha=0.2, color=lines[0].get_color())
plt.autoscale(tight=True)
plt.ylim(0, 1)

plt.show() 

从上图中可以看到，原本的扁平先验函数不再是平的！事实证明，扁平先验也含有信息。

杰弗里斯先验的意义是，建立一个不会因偶然改变变量位置而大幅变化的先验。 

0x2：主观先验

如果我们对先验的特定区域增大概率可能性，而对其他区域相应减小，这样便将我们的推断向具有更大可能性区域的参数偏倚。这被称为一个主观先验，或信息先验。

1. 对参数空间中特定区域有主观偏倚的先验分布

下图中，主观先验描述了一个信念，即未知参数可能位于0.5附近，而不是在极点，而客观鲜艳对此是不敏感的。

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats


colors = ["#348ABD", "#A60628", "#7A68A6", "#467821"]

x = np.linspace(0, 1)
y1, y2 = stats.beta.pdf(x, 1, 1), stats.beta.pdf(x, 10, 10)

p = plt.plot(x, y1,
             label='An objective prior \n(uninformative, \n"Principle of Indifference" )')
plt.fill_between(x, 0, y1, color=p[0].get_color(), alpha=0.3)

p = plt.plot(x, y2,
             label="A subjective prior \n(informative)")
plt.fill_between(x, 0, y2, color=p[0].get_color(), alpha=0.3)

p = plt.plot(x[25:], 2 * np.ones(25), label="another subjective prior")
plt.fill_between(x[25:], 0, 2, color=p[0].get_color(), alpha=0.3)

plt.ylim(0, 4)

plt.ylim(0, 4)
leg = plt.legend(loc="upper left")
leg.get_frame().set_alpha(0.4)
plt.title("Comparing objective vs. subjective priors for an unknown probability")

plt.show()

需要注意的是，使用主观先验不是总是意味着我们采用领域专家的主观意见。更多时候，主观先验是对之前问题的后验，可以理解为一个逐轮迭代优化的过程。在这轮统计中，我们基于新的数据更新这个后验。

从某种程度上来说，像SGD、最大熵增优化算法等优化算法，其运行过程往往都是由很多轮”设定先验-基于后验结果更新下一轮的先验“的子过程组成的。

2. 有限空间上的扁平先验

如果我们对二项式里的先验分布进行区间限制，即只有（0.5，1】内是均匀扁平分布，则这种先验就不能称之为客观先验，即使分布在（0.5，1】内为真正的”扁平“。

扁平先验必须在整个参数范围内扁平。

3. 从领域专家处获得先验分布

指定主观经验，是指将问题所在的领域知识结合到我们的数学框架中，融入领域知识可以带来许多用处：

• 有助于优化算法收敛：例如如果我们知道未知参数是严格为正的，那么我们可以缩小关注空间，从而省下在负值查找的时间；
• 允许更准确的推断：通过加权靠近真值的先验，我们收窄我们的最终推断；
• 更好地表达了我们的不确定性；

0x3：如何决策选择主观还是客观先验

选择客观或者主观先验主要取决于需要解决的问题，但也有少数情况会优先选择某种先验。
在科学研究中，科学家面对的往往是完全陌生的新领域，基于不会有任何的领域先验，这时候选择客观先验是显而易见的，因为这消除了结论中的任何偏见。它应该让两个对研究主题有不同的信念的研究人员仍然觉得同一个客观先验是”公平“的。

我们必须记住：选择先验，无论是主观的或者客观的，扔是建模过程的一部分。
引用格尔曼的话：在模型已经拟合之后，应该检查后验分布，看看它是否有意义的。如果后验分布没有意义，这意味着样本中额外的知识（信息）尚未包括在模型中，而且违背了已经使用的先验分布的假设。这时候合适的做法是再回头去改变先验分布，使之与外部知识更加一致。

如果后验分布看起来没有意义，那么显然你对后验分布应该是怎样的有了明确的认知，这意味着当前的先验分布不包括所有的初始信息，应该更新。此时我们可以放弃当前的先验，而选择一个更能反映我们所有初始信息的先验。
格尔曼同时建议，使用一个大边界的均匀分布在大多数情况下是一个很好的客观先验选择。
然后，人们应该警惕使用具有大边界的均匀客观先验，因为它们会对极端不敏锐的数据点分配过大的先验概率（因为实际情况可能是这些极端数据点应该被分配接近0的先验概率）。
带有大方差的正态随机变量、或者在严格为正（或负）的情况下选带有宽尾的指数变量可能是更好的选择。

0x4：经验贝叶斯 - 基于数据驱动的先验选择

经验贝叶斯是结合了频率论和贝叶斯推断的技巧。
贝叶斯方法和频率论方法之前显著的差别在于，贝叶斯方法有一个带超参数先验分布，而经验方法不具有任何先验的概念。
经验贝叶斯结合了这两种方法，即使用频率论方法来选择先验分布的超参数，然后用贝叶斯方法来进行后验推断。

1. 一个简单的例子 - 估计模型超参数

假设我们希望估计方差 σ = 5 的正态分布的参数 μ，因为 μ 的范围可以是所有实数，我们可以使用正态分布作为 μ 的先验（这个正态分布和要估计的模型没有关系）。接下来，我们必须选择先验的超参数，表示为，参数可以反映我们对先验的不确定性，对于，我们有两个选项

1. 经验贝叶斯：使用经验的样本均值，这将使先验居中于经验样本均值：

2. 传统贝叶斯推断：使用主观先验知识，或者更客观的先验（0均值和大的标准差）

相比于客观贝叶斯推断，经验贝叶斯可以说是半客观的，因为当前先验模式（模型）的选择是由我们确定的（因此是主观的），而模型参数则仅由数据来确定（因此是客观的）。

需要注意的是，使用这种方法需要非常慎重。从某种程度上来说，经验贝叶斯方法是对数据重复计数。
这么理解这句话，我们使用了曾经在先验中的数据，数据同时影响了先验和后验推断，因而将影响针对观测数据的结果，进而影响MCMC的推断引擎中的结果。
这种重复计数会低估我们真正的不确定性，为了减少这种重复计数，我们只能在有很多观测样本时才使用经验贝叶斯（大数定理），否则先验将有过于强烈的影响，同时，如果可能的话，要保持适当的不确定性（通过设置大方差或其他等价方法）

笔者插入： 经验贝叶斯也违反了贝叶斯推断的哲学，教科书中的贝叶斯算法为：先验 => 观察数据 -> 后验。而经验贝叶斯的算法为：观测数据 => 先验 => 观测数据 => 后验。

理想情况下，所有的先验应在观测数据之前决定，以使数据不影响我们先验的观察。

Relevant Link:

https://www.zhihu.com/question/67846877/answer/257216646
《思考：快与慢》丹尼尔·卡尼曼 - 锚定相关主题

 

3. 先验的表现形式

搞机器学习的同学对先验这个词应该不会陌生，它大量地出现在问题分析和模型设计过程中，例如：

• 在模型设计中引入正则化技术：从贝叶斯的角度来看，正则化等价于对模型参数引入先验分布。
• 在特征工程中结合领域经验设计特殊的特征向量：例如设计【身高、体重、握力、百米跑步速度】这几位特征维度用于检测学校男生还是女生，显然样本集在这种特征空间上是明显可分的。在特征工程的阶段就得到一个明显可分的样本集是非常值得高兴的，因为这意味着我们的项目已经成功了90%，接下来的工作即使用随机森林可能也能得到非常好的分类结果。
• 在贝叶斯推断中使用特定先验概率分布函数
• 模型参数初始化时人工设定特定的初始化值

笔者认为：先验的本质是约束。不管何种形式的先验，其本质都是对参数空间搜索函数的牵引和回拉作用，使得最终的搜索结果在一定程度上受初始先验的“影响”。

 

4. 贝叶斯分析中一些常用的先验

这个章节我们将介绍在贝叶斯分析和方法中常用的一些先验分布。 在开始讨论之前，笔者希望先抛出一个问题和大家一起思考。为什么gamma分布、泊松分布可以作为先验分布？有哪些标准决定了一个函数可以作为先验分布？

笔者认为这个问题可以从几个方面展开思考：

1. 思考gamma函数这类函数被发明的最早的原因，是因为其可以作为某个物理现象或当时的实际问题的数学模型；
2. 理论上，任意一个非负的可积函数，你都可以通过变换变成某个分布的密度函数；
3. 数学领域有句名言，数学是解开宇宙秘密的一把钥匙，很多数学公式，你第一眼看到就会从心底产生一个感觉，哇！好美啊！好工整好对称啊！这类数学模型往往可以抽象概括一大类物理现象，因此很自然地被作为先验模型使用。

0x1：二项分布 - N次二项实验（只有两种互斥结果的实验）中的离散概率密度函数分布

二项分布就是重复n次独立的伯努利试验，在每次试验中只有两种可能的结果，而且两种结果发生与否互相对立，并且相互独立，与其它各次试验结果无关，事件发生与否的概率在每一次独立试验中都保持不变，则这一系列试验总称为n重伯努利实验，当试验次数为1时，二项分布服从0-1分布。

1. 数学公式 - N次实验中发生K次事件的”累计“概率

二项分布（Binomial Distribution），即重复n次的伯努利试验（Bernoulli Experiment）。

如果在一次伯努利实验中事件发生的概率是p，则不发生的概率1-p。

N次独立重复试验中发生k次的概率是：

记作：ξ~B(n，p)。期望：Eξ = np；方差：Dξ = np(1-p)；

上图中可以很清楚看到，n、p 这2个参数是如何影响二项分布函数的形状的。

笔者插入：可以把人生比作二项分布，我们无时无刻不在进行选择和决策，而选择就是一种”向左走、向右走“的伯努利实验，要获得更大的期望（收获），就要尽量扩大n，也就是一分耕耘一分收获，另一方面，要尽可能提高在单个决策中的成功率，这样，总体的方差才会小，我们才有更大的可能达到我们的人生目标。

2. 公式的原理推导和理解

我们一起用一个抛硬币游戏来理解二项分布概率公式的内涵。

1）等概率伯努利实验

抛硬币的结果有两种，结果是正面（H）；是结果是反面（T）。

在没有特定定制硬币的情况下，我们说：硬币正面向上（H）的概率是 ½ ；反面向上（T）的概率是 ½。

游戏开始后掷骰子：

很显然，四 的概率是 1/6 （六面里有一面是四）。不是四的概率是 5/6 （六面里有五面不是四）。

现在问一个问题：抛一个公平的硬币三次 …… 结果是两个正面的概率是多少？

抛开数学公式的推导，我们先从实验统计的角度来回答这个问题，抛一个公平的硬币三次（H 是正面， T 是反面）可以有 8 个结果：

"两个正面" 可以是任何次序： "HHT"、"THH" 和 "HTH" 都有两个正面（和一个反面）。所以3个结果有 "两个正面"。

每个结果的可能性是一样的，一共有 8 个可能，所以每个结果的概率是 1/8

所以结果是"两个正面"这个事件的概率是：3  ×  1/8  =  3/8 

我们继续计算一下所有可能事件的发生概率：

• P(三个正面) = P(HHH) = 1/8
• P(两个正面) = P(HHT) + P(HTH) + P(THH) = 1/8 + 1/8 + 1/8 = 3/8
• P(一个正面) = P(HTT) + P(THT) + P(TTH) = 1/8 + 1/8 + 1/8 = 3/8
• P(没有正面) = P(TTT) = 1/8

 将概率分布画在直方图上：

注意到图是对称的，这个不是偶然，我们后面会谈到，这和在单个伯努利实验中的概率p有关。

我们进一步来思考上面这个推导过程的数学化表示，我们注意到N次实验中发生K次事件本质上是一个排列组合问题，具体哪K次并不关注，所以我们可以用数列阶乘的方式来计算总共的排列组合次数：

通常这个是叫 "n取k"。

我们重新计算一次抛 3次，结果有 2个 正面的排列组合：

所以结果是"两个正面"这个事件的概率是：3  ×  1/8  =  3/8 = 0.375

和我们人工统计得到的结果一样，当然在数字小的情况下可以使用统计方式，当数字很大的时候就需要通过公式进行计算了。

2）非等概率伯努利实验

在上面讲的情况里，成功和失败的可能性是相同的。可是，当硬币有偏误（不公平的硬币）时，一面的可能性便会大于另一面。例如正面是70%，反面是30%。

还是进行3次实验，我们通过树图的方式来看一下所有可能事件的概率结果：

图中突出显示了有 "两个正面" 的结果。

将概率分布画上直方图上：

可以看到直方图出现了偏置，这个偏置是由于概率p的不均等造成的。

将3个结果相加，得到"两个正面"这个事件的概率是：3  ×  0.147  =  0.441。

接下来用数学公式来重新计算这个过程。

0.7 是正面的概率，称它为 p；
其他的结果（反面）的概率是：1-p；

一共有 n 个结果；

2 是目标的结果个数，称它为 k； 
得到我们 "两个正面" 的概率是：pk；  

其他的结果的个数是：n-k； 
"其他的结果" 的概率是：(1-p)(n-k)；

所以所有的结果一同发生的概率是：

每个结果的概率

结果的总数是：

得到"两个正面"这个事件的概率是：3  ×  0.147  =  0.441。  

3. 应用条件

• 各观察单位只能具有相互对立的一种结果，如阳性或阴性，生存或死亡等，属于两分类资料；
• 已知发生某一结果（阳性）的概率为π，其对立结果的概率为1-π，实际工作中要求π是从大量观察中获得比较稳定的数值；
• n次试验在相同条件下进行，且各个观察单位的观察结果相互独立，即每个观察单位的观察结果不会影响到其他观察单位的结果。如要求疾病无传染性、无家族性等； 

我们在使用二项分布对问题场景进行建模前需要特别思考上述条件，问题域是否满足二项分布的先决条件。

同时还有一点需要读者朋友们注意，二项分布和接下要要讨论的泊松分布等分布，都是对顺序不敏感的，即对N次实验中发生K次的顺序是不敏感的，如果你的问题场景中，还需要考虑到这K次具体发生的顺序，则无法直接套用原始的二项分布公式。

Relevant Link:

https://www.shuxuele.com/data/binomial-distribution.html

0x2：泊松分布 - n趋向无穷时的离散二项分布

Poisson分布，是一种离散概率分布，由法国数学家西莫恩·德尼·泊松（Siméon-Denis Poisson）在1838年时发表。

1. 数学公式

注意：泊松分布的公式和二项分布的公式是相同的，区别在于泊松分布在n趋向于无穷时，公式可以近似等于另一个计算等式。

泊松分布的参数λ是单位时间(或单位面积)内随机事件的平均发生次数。

泊松分布的期望和方差均为λ。特征函数为：

 

不同λ取值下，函数分布曲线图如下：

上图中可以很清楚看到，λ（样本均值） 这个参数是如何影响泊松分布函数的形状的。

2. 公式的原理推导和理解

公司楼下有家馒头店，每天早上六点到十点营业，生意挺好，就是发愁一个事情，应该准备多少个馒头才能既不浪费又能充分供应？

老板决定采用数据驱动的方式解决这个问题，先采集数据！老师统计了一周每日卖出的馒头：

均值为：

按道理讲均值是不错的选择（基于最小二乘的损失函数），但是如果每天准备5个馒头的话，从统计表来看，至少有两天不够卖， 的时间不够卖：

那怎么办呢？很显然，我们不能按照MAX值来准备，那么确实保险程度最高，但是浪费粮食的概率也很高。很显然，这里需要基于概率推断的思想，选择一个有较高概率（95%）馒头刚好够供给的数字，问题是这个数字怎么算出来呢？

老板的解决策略是：用概率统计的方式解决问题，首先先将问题转化为二项实验问题，在二项分布的框架内进行思考和计算：

老板尝试把营业时间抽象为一根线段，把这段时间用  来表示：

然后把周一的三个馒头按照销售时间放在线段上：

把  均分为四个时间段：

此时，在每一个时间段上，要不卖出了（一个）馒头，要不没有卖出：

可以看到在每个时间段，就有点像抛硬币，要不是正面（卖出），要不是反面（没有卖出）。此时，我们已经将某一天卖出多少馒头的问题，转换为了二项分布问题，我们来逐条回顾下二项分布的达成条件：

• 每次试验是独立的：T内的每个时间段是否卖出馒头都是条件独立的；
• 每个试验只有两个可能结果：卖出 or 不卖出；
• 每个试验里的 "成功" 概率是不变的：我们假定顾客购买这个行为是固定概率的；

这样， 内卖出3个馒头的概率，就和抛了4次硬币（4个时间段），其中3次正面（卖出3个）的概率一样了。

这样的概率通过二项分布来计算就是：

老板继续思考，但是，当他把周二的七个馒头放在线段上，分成四段就不够了：

从图中看，每个时间段，有卖出3个的，有卖出2个的，有卖出1个的，就不再是单纯的“卖出、没卖出”了。二项分布的结果唯一性条件不满足了，不能套用二项分布了。

解决这个问题也很简单，把  分为20个时间段，那么每个时间段就又变为了抛硬币：

这样， 内卖出7个馒头的概率就是（相当于抛了20次硬币，出现7次正面）：

老板继续思考，为了能一次性保证在一个时间段内只会发生“卖出、没卖出”，干脆把时间切成  份，而且分得越细越好，用极限来表示：

接下来的是问题，在单个时间格内的概率p是多少？也即每次实验中顾客是否购买馒头的概率是多少？

在上面的假设下，问题已经被转为了二项分布。二项分布的期望为：

那么：

那么 μ 怎么计算得到呢？

一个简单的回答是，可以用历史样本均值（期望）计算得到：

可以用它来近似：

笔者插入：

上面这样计算的假设前提是，我们可以基于历史数据来得出一个概率值，但是这里其实隐含了一个问题，只有在历史数据满足大数定理的时候（N很大，历史统计数据足够多）这个公式结果才是无偏的，否则这个计算结果可能是一个有偏结果，这个问题不影响我们理解概念，但是读者朋友在实际项目中要注意。大数定理！大数定理！大数定理！

有了 了之后，我们继续推导上面的二项分布公式，看看如何从二项公式演进到泊松分布公式：

我们来算一下这个极限：

其中：

所以：

上面就是泊松分布的概率密度函数，也就是说，在  时间内卖出  个馒头的概率为：

一般来说，我们会换一个符号，让  ，所以：

这就是教科书中的泊松分布的概率密度函数，公式中λ是一个形状参数，也等于均值和方差，可以这么理解，先验样本中的均值和方差，决定对应泊松分布的形状。

带入之前计算得到的 λ = 5，得到泊松分布公式：

画出概率密度函数的曲线：

接下来就是概率推断的问题了，需要反复强调的是，在概率推断领域，我们无法得到一个100%的实值结果，而只能得到在某个结果下的概率，我们基于这个概率（例如大于95%）进行一个推断和决策。

所以，基于上述泊松公式，老板接下来需要思考的问题是：我需要知道在93%的概率下，我准备的馒头足够供给，这就足够满足我店铺的正常盈利了，另外7%算作小概率事件，可以忽略不计。

可以看到，如果每天准备8个馒头的话，那么足够卖的概率就是把前8个的概率加起来：

最终的决策就是：每天准备8个馒头！问题得以解决。

3. 物理意义

泊松分布适合于描述单位时间（或空间）内随机事件发生的次数的概率。如某一服务设施在一定时间内到达的人数，电话交换机接到呼叫的次数，汽车站台的候客人数，机器出现的故障数，自然灾害发生的次数，一块产品上的缺陷数，显微镜下单位分区内的细菌分布数等等。
在实际事例中，当一个随机事件，例如某电话交换台收到的呼叫、来到某公共汽车站的乘客、某放射性物质发射出的粒子、显微镜下某区域中的白血球等等，以固定的平均瞬时速率λ（或称密度）随机且独立地出现时，那么这个事件在单位时间（面积或体积）内出现的次数或个数就近似地服从泊松分布P(λ)。因此,泊松分布在管理科学、运筹学以及自然科学的某些问题中都占有重要的地位。（在早期学界认为人类行为是服从泊松分布，2005年在nature上发表的文章揭示了人类行为具有高度非均匀性。）

Relevant Link:

https://blog.csdn.net/ccnt_2012/article/details/81114920
https://baike.baidu.com/item/%E6%B3%8A%E6%9D%BE%E5%88%86%E5%B8%83/1442110?fr=aladdin
https://blog.csdn.net/ccnt_2012/article/details/81114920 

0x3：正态分布 - 连续概率密度分布

正态分布（Normal distribution），也称“常态分布”，又名高斯分布（Gaussian distribution），最早由A.棣莫弗在求二项分布的渐近公式中得到。C.F.高斯在研究测量误差时从另一个角度导出了它。

正态曲线呈钟型，两头低，中间高，左右对称因其曲线呈钟形，因此人们又经常称之为钟形曲线。

1. 数学公式

若随机变量 X 服从一个位置参数为 μ、尺度参数为 σ 的概率分布，且其概率密度函数为：

则这个随机变量就称为正态随机变量，正态随机变量服从的分布就称为正态分布，记作

数学期望为μ、方差为σ^2。

特别的，当时，正态分布就成为标准正态分布：

2. 图形特征

上图中可以很清楚看到，μ， σ 这两个参数是如何影响正态分布函数的形状的。

• 集中性：正态曲线的高峰位于正中央，即均数所在的位置。
• 对称性：正态曲线以均数为中心，左右对称，曲线两端永远不与横轴相交。
• 均匀变动性：正态曲线由均数所在处开始，分别向左右两侧逐渐均匀下降。

3. 物理意义

正态分布有极其广泛的实际背景，生产与科学实验中很多随机变量的概率分布都可以近似地用正态分布来描述。一般来说，如果一个量是由许多微小的独立随机因素影响的结果，那么就可以认为这个量具有正态分布。

在笔者所在的网络安全领域，我们可以假设一台机器在历史上一个月的网络QPS满足正态分布，因为在每个时刻，网络QPS是由很多因素累计决定的（网络抖动、应用负载变化、网络入侵等），通过正态分布对历史的网络QPS进行建模，从而对未来的网络外联QPS进行异常检测。 

Relevant Link: 

https://www.zhihu.com/question/26854682
https://baike.baidu.com/item/%E6%AD%A3%E6%80%81%E5%88%86%E5%B8%83/829892?fr=aladdin

0x4：Gamma（伽马）分布 - 连续概率函数

伽玛分布（Gamma Distribution）是统计学的一种连续概率函数，记为

它是一个正实数的随机变量，它是概率统计中一种非常重要的分布。

1. 数学公式

密度函数为：

期望和方差分别为：

Gamma分布中的参数α称为形状参数（shape parameter），β称为逆尺度参数（scale parameter）。

2. 变化趋势

特殊的，

当 α=1 时，伽马分布就是参数为 β 的指数分布，X~Exp（β）;

当 α=n/2，β=2时，伽马分布就是自由度为n的卡方分布，X^2(n);

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats

figsize(12.5, 5)
gamma = stats.gamma

parameters = [(1, 0.5), (9, 2), (3, 0.5), (7, 0.5)]
x = np.linspace(0.001, 20, 150)
for alpha, beta in parameters:
    y = gamma.pdf(x, alpha, scale=1. / beta)
    lines = plt.plot(x, y, label="(%.1f,%.1f)" % (alpha, beta), lw=3)
    plt.fill_between(x, 0, y, alpha=0.2, color=lines[0].get_color())
    plt.autoscale(tight=True)

plt.legend(title=r"$\alpha, \beta$ - parameters")

plt.show()

Relevant Link: 

https://www.cnblogs.com/coshaho/p/9653460.html
https://www.cnblogs.com/alps/p/5601300.html
https://www.zhihu.com/question/34866983/answer/60541847

0x5：威沙特分布

除了标量的随机变量之外，随机变量还可以是矩阵形式的！具体地说，威沙特分布是所有半正定矩阵的分布。

合适的协方差矩阵是正定的，因此该威沙特分布是一个协方差矩阵的适当的先验，我们在下图中，绘制一个来自 4 x 4 和 15 x 15 威沙特分布的某些实现

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats

figsize(12.5, 5)

n = 4
for i in range(10):
    ax = plt.subplot(2, 5, i + 1)
    if i >= 5:
        n = 15
    plt.imshow(stats.wishart.rvs(n + 1, np.eye(n)), interpolation="none",
               cmap="hot")
    ax.axis("off")

plt.suptitle("Random matrices from a Wishart Distribution")

plt.show()

有一点需要注意的是这些矩阵的对称性，它反映了协方差的对称性。

0x6：Beta分布 

1. 数学公式 

在概率论中，贝塔分布，也称B分布，是指一组定义在（0，1）区间的连续概率分布，有两个参数α，β。

Β分布的概率密度函数是：

2. 变化趋势

α，β参数决定了Beta分布的形状：

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats

figsize(12.5, 5)

params = [(2, 5), (1, 1), (0.5, 0.5), (5, 5), (20, 4), (5, 1)]

x = np.linspace(0.01, .99, 100)
beta = stats.beta
for a, b in params:
    y = beta.pdf(x, a, b)
    lines = plt.plot(x, y, label = "(%.1f,%.1f)"%(a,b), lw = 3)
    plt.fill_between(x, 0, y, alpha = 0.2, color = lines[0].get_color())
    plt.autoscale(tight=True)
plt.ylim(0)
plt.legend(loc = 'upper left', title="(a,b)-parameters")

plt.show()

上图中有一点请读者朋友注意，当指定参数为（1，1）时，Beta分布等同于扁平分布，它是一种均匀分布。

因此，Beta分布是均匀分布的更一般形式。

3. Beta分布特性

1. 作为概率的概率分布，Beta(a, b)在(0, 1)上对 θ 积分必定为1；
2. Beta(a, b)同时能作为先验分布和后验分布，能够模拟各种概率分布情况；
3. Beta分布可以模拟出以(0, 1)上任意点为峰值的曲线，这表明Beta分布可以模拟极大似然法求出的任意最大值点概率值；

Relevant Link: 

https://baike.baidu.com/item/%E8%B4%9D%E5%A1%94%E5%88%86%E5%B8%83/8994021?fr=aladdin 
http://www.cnblogs.com/coshaho/p/9658135.html

0x7：共轭先验

1. 什么是共轭先验

在贝叶斯概率理论中，如果后验概率P(θ|x)和先验概率p(θ)满足同样的分布律，那么，先验分布和后验分布一起被叫做共轭分布（彼此是对方的共轭），同时，先验分布叫做似然函数的共轭先验分布。

假定 X 来自或被认为是来自一个著名的分布，称之为 fa，其中 a 可能是 f 的未知参数（f 可以是一个正态分布或者二项分布）。对于特定的 fa，如果存在先验分布 pβ，使得：

其中，是一组不同的参数，但 p 是和先验相同的分布。

满足该关系的先验 p 称为共轭先验。

2. 共轭先验的性质 

• 代数上的方便性，可以在建模时直接给出后验分布的封闭形式，否则的话只能数值计算（例如SGD启发式搜索）。
• 可以形成一个先验链，即现在的后验分布可以作为下一次计算的先验分布，如果形式相同，就可以形成一个链条。
• 一切指数族分布成员都存在一个共轭先验。
  • 在伯努利分布中，共轭先验是 Beta 分布；　　
  • 在高斯分布中，均值的共轭先验是高斯分布，精度的共轭先验是Wishart分布；

3. 共轭先验的局限

• 共轭先验是不客观的，因此，它只有当需要主观先验时才有用。不能保证共轭先验能够顾忌开发者的主观意识；
• 对于简单的一维问题，通常存在共轭先验。但是对于更大更复杂的问题，涉及更复杂的结构，基本没希望找到共轭先验；

4. 一个例子：Beta先验+二项分布数据 = Beta后验分布 

Beta分布和二项分布之间有一个有趣的关系。假设我们感兴趣的是一些未知的比例或者概率 p。我们设定它符合一个Beta(α，β) 先验分布，然后我们通过观察收集一个由二项式过程 X~B(N，p) 昌盛的数据，其中 p 仍然是未知的，后验分布仍然是Beta分布，即：p | X ~ Beta（a+X，β+N-X）

一个Beta先验分布连同二项式生成的观测数据形成一个Beta后验分布。这是一个非常有用的性质，无论是从计算的角度还是启发性的角度。

具体地说，如果我们设 p（一个均匀分布）的先验为 Beta（1，1），观察 X ~ B（N，p）的数据，我们得到的后验是 Beta（1+X，1+N-X）。

例如，如果我们在 N=25 次试验里观察到 X=10 次成功，那么我们关于 p 的后验是 Beta（1+10，1+25-10）= Beta（11，16）的分布。

Relevant Link:  

https://www.jianshu.com/p/8de4fd24e53c 

0x8：各个分布之间的推演和转化关系：

1. 二项分布和泊松分布的关系

在 n 重伯努利实验中，事件A在每次实验中发生的概率为 Pn，在小数据情况下，Pn与实验次数有关。如果= λ > 0，则泊松公式可由二项分布公式推导得到新的表达等式：

由定理可知，当二项分布 b(n, p) 的参数 n 很大，p 很小，而 λ = np 大小适中时（实际中，n=>100, p<=0.1, np<=10时），二项分布可用参数为 λ=np 的泊松分布来近似（也即满足大数定理时）。

这就是二项分布的泊松逼近，当然 n 应该尽可能地大，否则近似效果往往不佳。

二项分布的泊松近似常常被应用于研究稀有事件（即每次实验中事件出现的概率p很小），当伯努利实验的次数n很大时，事件发生的频数分布近似于泊松分布的频数分布。
实际表明，在一般情况下，当 p<0.1 时，这种近似是很好的，甚至n不必很大都可以，当然p越小，n相对就会增大，这是期望公式的限定。下图中，当 p=0.01，n=2，我们对比了二项分布和泊松分布的频数分布

通过概率分布图可以看出这个近似，当二项分布的  很小的时候，两者比较接近：

也可以从另一个方面理解，当 n 趋向无穷大的时候，才能满足将随机事件切位为单独的小的时间格，在每个时间块中，独立伯努利条件才能成立，这个时候二项分布进化为泊松分布。

二项分布和泊松分布的关系，就是从小数据到大数据时，大数定律逐渐占据主导地位的转变。在N趋近于无穷时，统计可以近似等同于概率。

2. 二项分布和正态分布的关系

如果n很大，函数分布最终成正态分布，二项分布的极限分布为正态分布。故当n很大时，二项分布的概率可用正态分布的概率作为近似值。

3. 泊松分布和正态分布的关系 

二项分布既可以用泊松分布近似，也可以用正态分布近似。显然，泊松分布和正态分布在一定条件下也具有近似关系。

对任意的 a < b：

，其中，

4. 用泊松分布还是正态分布来近似二项分布？

二项分布的泊松近似和正态近似各自适用的条件是不同的。 

• 当 p 很小时，即使 n 不是很大，用泊松分布近似二项分布，已经相当吻合。但是在这种情况下，如果用正态分布来近似就会造成较大的误差；
• 当 n 充分大，p 既不接近于0也不接近于1时（最好满足 0.1<=p<0.9），正态分布可以较好地近似二项分布；

Relevant Link:

https://wenku.baidu.com/view/5d9ac306e87101f69e3195e4.html
https://www.zybang.com/question/dcdea243cfc493eac4dde9eb4c05b373.html
https://baike.baidu.com/item/%E4%BC%BD%E9%A9%AC%E5%88%86%E5%B8%83/7245468?fr=aladdin
https://baike.baidu.com/item/Gamma%E5%88%86%E5%B8%83/1033808
http://www.cnblogs.com/JustForCS/p/5264315.html 
https://blog.csdn.net/u012279165/article/details/73693157
https://baike.baidu.com/item/%E4%BA%8C%E9%A1%B9%E5%88%86%E5%B8%83/1442377?fr=aladdin
http://www.360doc.com/content/17/1231/22/9200790_718001949.shtml 

 

5. 从针对伯努利二项分布密度估计到Beta分布密度估计的演进，看先验和后验的关系

0x1：从二元变量伯努利分布说起

首先，我们考虑一个二元随机变量x ∈ {0, 1}。例如，x可能描述了扔硬币的结果，

• x = 1表 示“正面”
• x = 0表示反面

我们可以假设由一个损坏的硬币，这枚硬币正面朝上的概率未必等于反面朝上的概率。x = 1的概率被记作参数μ，因此

其中0 ≤ μ ≤ 1。我们可以看到，p(x = 0 | μ) = 1 − μ。x的概率分布因此可以写成：

这被叫做伯努利分布(Bernoulli distribution)。很容易证明，这个分布是归一化的，并且均值和方差为：

现在我们假设我们有一个x的观测值的数据集D = {x₁, . . . , x_N}。假设每次观测都是独立地 从p(x | μ)中抽取的，因此我们可以构造关于μ的似然函数如下：

在频率学家的观点看来，我们可以通过最大化似然函数来估计μ的值，或者等价地，最大化对数似然函数。在伯努利分布的情形下，对数似然函数为：

在频率学派的这种观点中，对数似然函数只通过和式 ∑_nx_n 依赖于x_n的N次观察。这个和式是这个分布下数据的充分统计量(sufficient statistic)。

如果我们令 lnp(D | μ)关于μ的导数等于零，我们就得到了最大似然的估计值：

这也被称为样本均值(sample mean)。如果我们把数据集里x = 1(正面朝上)的观测的数量记作m，那么我们可以把上式写成下面的形式：

因此在最大似然的框架中，正面朝上的概率是数据集里正面向上的观测所占的比例，完全由数据决策密度估计结果。

现在假设我们扔一个硬币3次，碰巧3次都是正面朝上。那么N = m = 3，且μML = 1。这种情况下，最大似然的结果会预测所有未来的观测值都是正面向上。常识告诉我们这个是不合理的。事实上，这是最大似然中过拟合现象的一个极端例子。

0x2：Beta分布，在伯努利分布基础之上增加了先验约束

根据上一章的讨论，我们已经看到伯努利分布的参数μ的最大似然解，就是数据集里x = 1的观测所占的比例。

正如我们已经提到过的那样，这对于小规模的数据集会给出严重的过拟合结果。

为了用贝叶斯的观点重新看待这个问题，我们需要引入一个关于μ的先验概率分布p(μ)。

这里，我们考虑一种形式简单的先验分布。这种形式简单的先验分布有很多有用的性质。

为了找到这个先验分布，我们注意到似然函数是某个因子与μ^x(1 − μ)^1−x的乘积的形式。如果我们选择一个正比于μ和(1 − μ)的幂指数的先验概率分布， 那么后验概率分布（正比于先验和似然函数的乘积）就会有着与先验分布相同的函数形式。这个性质被叫做共轭性（conjugacy）。因此，我们把这种先验分布选择为Beta分布，定义为：

其中，Γ(x)是Gamma函数，上面公式保证了Beta分布式归一化，即：

Beta分布的均值和方差为：

参数a和b经常被称为超参数(hyperparameter)，因为它们控制了参数μ的概率分布。下图给出了不同的超参数值对应的Beta分布的图像：

μ的后验概率分布现在可以这样得到：把Beta先验与二项似然函数相乘，然 后归一化。只保留依赖于μ的因子，我们看到后验概率分布的形式为：

其中l = N − m，即对应于硬币“反面朝上”的样本数量。我们看到上式关于μ的函数形式与先验分布相同，这反映出先验关于似然函数的共轭性质。实际上，它仅仅是另一个Beta分布。 

从公式上来分析，我们看到，如果一个数据集里有m次观测为x = 1，有l次观测为x = 0，那么从先验概率到后验概率，

• a的值增大了m
• b的值增大了l

这让我们可以简单地把先验概率中的超参数a和b分别 看成x = 1和x = 0的有效观测数（effective number of observation）。注意，a和b不一定是整数。 

接下来，如果我们继续观测到更多的数据，那么后验概率分布可以扮演先验概率的角色。

为了说明这一点，我们可以假想每次值取一个观测值，然后在每次观测之后更新当前的后验分布。 更新方法是让当前的后验分布与新观测值的似然函数相乘，然后归一化，获得新的修正后的后验分布。

在每个阶段，后验概率是一个Beta分布，对于x = 1和x = 0的观测总数(先验的和实际的)由参数a和b给出。观测到一个x = 1仅仅对应于把a的值增加1，而观测到x = 0会使b增加1。下图说明了这个过程中的一个步骤。 

贝叶斯顺序推断中的一个步骤的例子。先验概率为Beta分布，参数为a = 2, b = 2，似然函数中，N = m = 1，对应于x = 1的一次观测，从而后验概率分布为Beta分布，参数 为a = 3, b = 2。 

我们看到，如果我们使用了贝叶斯观点，那么学习过程中的顺序(sequential)方法可以自然而然地得出。它与先验和似然函数的选择无关，只取决于数据独立同分布的假设。

顺序方法每次使用一个观测值，或者每次使用一小批观测值，然后在使用下一个观测值之前丢掉它们。

例如，顺序方法可以被用于实时学习的场景中。在实时学习的场景中，输入为一个稳定持续的数据流，模型必须在观测到所有数据之前就进行预测。由于顺序学习的方法不需要把所有的数据都存储到内存里，因此顺序方法对于大的数据集也很有用。最大似然方法也可以转化成顺序的框架。 

如果我们的目标是尽可能好地预测下一次试验的输出，那么我们必须估计给定观测数据集D的情况下，x的预测分布。根据概率的加和规则和乘积规则，这个预测分布的形式为：

这个结果可以简单地表述为对应于x = 1的观测结果（包括实际的观测值和假想的先验观测值）所占的比例。

注意，在数据集无限大的极限情况下，m, l → ∞，此时上式的结果变成了最大似然的结果。贝叶斯的结果和最大似然的结果在数据集的规模趋于无穷的情况下会统一到一起。这是一个很普遍的情况。对于有限规模的数据集，μ的后验均值总是位于先验均值和上式给出的μ的最大似然估计之间。

从下图中，我们可以看到，当观测的数量增加时，后验分布的图像变得更尖了。这通过Beta分布方差的结果也能够看出来，如果a → ∞或 者b → ∞，那么方差就趋于零。

实际上，这个性质是贝叶斯学习的一个共有的属性：

随着我们观测到越来越多的数据，后验概率表示的不确定性将会持续下降

 

6. 数据量N和先验分布对最终后验的综合决策

0x1：通过乘积公式体现数据量和先验分布对后验估计的动态制衡作用

在机器学习中有一个被广泛接受的观点，我们拥有的数据越多、数据量越多，先验就越不重要，这是符合直觉的。毕竟，我们的先验也是基于以前的信息，足够多的新信息完全可以替代我们以前信息的价值，因为规律永远蕴含在数据中，只要拥有数据，就可以不断从中提取出信息。

同时，足够多的数据对先验的修正也是有帮助的，如果我们的先验明确是错误的，那么数据的自我修正性质将呈现给我们一个不那么错的后验估计结果。

我们可以从数学上阐述上面的观点。给定数据集 X，对参数 θ 的后验分布可以写作：

写成对数形式：

对数似然函数会随着样本量而变化，因为它是数据的一个函数；

但是先验的密度函数不会随着数据而变化；

因此，当样本量增加时，的绝对值会变大，但保持不变。

因此，随着样本量增加，整体函数更多地受到的影响，所选择的先验的影响会变小。

因此，只要非零概率的区域是相同的，那么推断的收敛和先验无关。

0x2：选择退化的先验

只要先验在某个区域有非零的概率，通过N数据量的训练后，后验就可以在这个区域有任何可能的概率。

但是！当某个区域先验概率初始值为0时，无论输入多少的数据，后验都无法在这个区域得到任何概率了。从数学公式上很容易理解，这是由于乘法的性质决定的。

我们用一个小实验来说明，假设我们的数据是伯努利分布，我们希望估计p（成功的概率），我们现在选择一个”不合适“的先验 Uniform（0.5，1），这里说不合适是应该我们事先知道了数据的实际分布，真实项目场景中当然不可能有这种好事，这里仅仅是为了说明先验分布选错了会带来什么影响。

我们已知了数据的真实分布，但是我们选的先验在真实值0.35处的概率为0，我们来看下mcmc推断的结果会如何：

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats
import pymc as pm

figsize(12.5, 15)

p_actual = 0.35
x = np.random.binomial(1, p_actual, size=100)
print x[:10]

p = pm.Uniform('p', 0.5, 1)
obs = pm.Bernoulli('obs', p, value=x, observed=True)

mcmc = pm.MCMC([p, obs])
mcmc.sample(10000, 2000)

p_trace = mcmc.trace('p')[:]
plt.xlabel('Value')
plt.ylabel('Density')
plt.hist(p_trace, bins=30, histtype='stepfilled', normed=True)


plt.show()

从上图中可以看到，后验分布大量堆积在先验的下界。在数据的作用下，后验分布在”极力“靠近真值，但是因为先验分布的下界之外是0概率，后验概率无法改变。

如果在实际项目中看到了类似的情况，很有可能说明你的先验假设不太正确。

0x3：一个例子说明数据对先验的修正作用

下面通过一个例子来说明本小节观点。考察两个二项分布参数 θ 的后验的收敛，一个是扁平先验，一个是朝着 0 偏移的先验。当样本量增加时，它们的后验收敛，因此其推断也收敛。

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats
import pymc as pm

figsize(12.5, 15)

p = 0.6
beta1_params = np.array([1.,1.])
beta2_params = np.array([2,10])
beta = stats.beta

x = np.linspace(0.00, 1, 125)
data = stats.bernoulli.rvs(p, size=500)

plt.figure()
for i,N in enumerate([0,4,8, 32,64, 128, 500]):
    s = data[:N].sum()
    plt.subplot(8,1,i+1)
    params1 = beta1_params + np.array([s, N-s])
    params2 = beta2_params + np.array([s, N-s])
    y1,y2 = beta.pdf(x, *params1), beta.pdf( x, *params2)
    plt.plot(x,y1, label = r"flat prior", lw =3)
    plt.plot(x, y2, label = "biased prior", lw= 3)
    plt.fill_between(x, 0, y1, color ="#348ABD", alpha = 0.15)
    plt.fill_between(x, 0, y2, color ="#A60628", alpha = 0.15)
    plt.legend(title = "N=%d" % N)
    plt.vlines(p, 0.0, 7.5, linestyles = "--", linewidth=1)
    #plt.ylim( 0, 10)#


plt.show()

Relevant Link: 

《贝叶斯方法》 

 

7. 贝叶斯多臂游戏机问题

0x1：游戏规则

假设你面对10台游戏机机（多臂游戏机）。每台游戏机会以某种概率发奖金，每台游戏机的奖金相同，只是概率不同。有些游戏机非常大方，有些则很少。当然，游戏参与者事先不知道这些概率。

我们每次仅能选择一个游戏机，我们的任务是制定一个策略，赢取最多的奖金。

当然，如果我们知道哪台游戏机拥有最大的概率，然后总是挑这台，必定会产生最多的奖金。因此，我们的任务可以表述为”尽快找出最好的游戏机“。

该任务因为游戏机的随机性而变得复杂，在偶然情况下，次优的游戏机也可以返回许多奖金，这可能使得我们相信，这就是最优的那台。同样，在偶然的情况下，最好的游戏机也可能返回很低的奖金。

我们是应该继续尝试那台在本轮失败的机器，还是放弃挑选另一台？

一个更为棘手的问题是，如果我们发现了一台返回奖金相当不错的游戏机，我们是继续依靠它维持我们相当不错的成绩，还是尝试其他机器以期找到一个更好的游戏机？

这就是著名的探索与利用困境。

0x2：该问题的现实意义

探索与利用困境并不是数学家虚构的数字游戏，它在我们的日常生产生活中处处可见。

1. 互联网展示广告：公司有一系列可以展示给潜在客户的广告，但该公司并不清楚要遵循哪些广告策略，以最大限度地提高销售。
2. 生态学：动物只有有限的能量用于耗费，而且某些行为带来的回报是不确定的。动物如何最大化其适应度？
3. 金融：在随时间变化的回报量中，哪些股票期权能给出最高的回报？
4. 临床试验：一位研究人员希望在众多的方案中找出最好的治疗方案，同时最大限度地减少损失。
5. 心理学：赏罚如何影响我们的行为？人类如何学习？

0x3：游戏策略

1. 选择先验分布

该算法开始于一个完全无知的状态，它什么都不知道，并开始通过测试系统来获取数据。在获取数据和结果上，它可以学习什么是最好的和最差的行为。 

贝叶斯解决方案首先假定每个游戏机发奖金的先验概率。因为我们假定对这些概率完全无知，所以自然的，我们采用0到1的扁平分布（Beta分布）。

2. 算法流程

我们将10台游戏机抽象为x轴上【0，9】10个坐标数字，对应的，每个游戏机本轮的抽奖结果作为y值，这样，所有游戏机的抽奖作为就【x，y】坐标化了。

1. 首轮游戏：对所有游戏机（这里N=10台）设定一个扁平先验，也即初始化阶段是零知识的，对10台游戏机随机进行一次抽取即可；
2. 获取本轮样本数据：选择本轮抽取中，样本值最高的游戏机 b，即选择 B = argmax Xb，根据那个样本值最高的游戏机 b 的样本结果，作为本轮试验的样本数据，【x，y】，x代表第几胎游戏机，y代表对应的值。
3. 更新后验：基于本轮的样本数据【x，y】更新先验分布，这可以理解为一个后验修正过程；
4. 重复2-3过程；

这个算法包含的思想是：我们不应该直接放弃目前结果不理想的游戏机，而是随着我们建立的信念认为还有更好的游戏机，应该以一定的下降概率去选择它们。随着玩的次数逐渐增多，不好的游戏机的概率会下降，好的游戏机的概率会上升。

3. MCMC推断所有游戏机的发奖率后验分布

我们在代码中人工设定的真值隐含概率为：[0.85, 0.60, 0.75]。

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats
import pymc as pm

rand = np.random.rand


class Bandits(object):
    """
    This class represents N bandits machines.

    parameters:
        p_array: a (n,) Numpy array of probabilities >0, <1.

    methods:
        pull( i ): return the results, 0 or 1, of pulling
                   the ith bandit.
    """

    def __init__(self, p_array):
        self.p = p_array
        self.optimal = np.argmax(p_array)

    def pull(self, i):
        # i is which arm to pull
        return np.random.rand() < self.p[i]

    def __len__(self):
        return len(self.p)


class BayesianStrategy(object):
    """
    Implements a online, learning strategy to solve
    the Multi-Armed Bandit problem.

    parameters:
        bandits: a Bandit class with .pull method

    methods:
        sample_bandits(n): sample and train on n pulls.

    attributes:
        N: the cumulative number of samples
        choices: the historical choices as a (N,) array
        bb_score: the historical score as a (N,) array
    """

    def __init__(self, bandits):
        self.bandits = bandits
        n_bandits = len(self.bandits)
        self.wins = np.zeros(n_bandits)
        self.trials = np.zeros(n_bandits)
        self.N = 0
        self.choices = []
        self.bb_score = []

    def sample_bandits(self, n=1):
        bb_score = np.zeros(n)
        choices = np.zeros(n)

        for k in range(n):
            # sample from the bandits's priors, and select the largest sample
            choice = np.argmax(np.random.beta(1 + self.wins, 1 + self.trials - self.wins))

            # sample the chosen bandit
            result = self.bandits.pull(choice)

            # update priors and score
            self.wins[choice] += result
            self.trials[choice] += 1
            bb_score[k] = result
            self.N += 1
            choices[k] = choice

        self.bb_score = np.r_[self.bb_score, bb_score]
        self.choices = np.r_[self.choices, choices]
        return


figsize(11.0, 10)

beta = stats.beta
x = np.linspace(0.001,.999,200)

def plot_priors(bayesian_strategy, prob, lw = 3, alpha = 0.2, plt_vlines = True):
    ## plotting function
    wins = bayesian_strategy.wins
    trials = bayesian_strategy.trials
    for i in range(prob.shape[0]):
        y = beta(1+wins[i], 1 + trials[i] - wins[i])
        p = plt.plot(x, y.pdf(x), lw = lw)
        c = p[0].get_markeredgecolor()
        plt.fill_between(x,y.pdf(x),0, color = c, alpha = alpha,
                         label="underlying probability: %.2f" % prob[i])
        if plt_vlines:
            plt.vlines(prob[i], 0, y.pdf(prob[i]) ,
                       colors = c, linestyles = "--", lw = 2)
        plt.autoscale(tight = "True")
        plt.title("Posteriors After %d pull" % bayesian_strategy.N +\
                    "s"*(bayesian_strategy.N > 1))
        plt.autoscale(tight=True)
    return

hidden_prob = np.array([0.85, 0.60, 0.75])
bandits = Bandits(hidden_prob)
bayesian_strat = BayesianStrategy(bandits)

draw_samples = [1, 1, 3, 10, 10, 25, 50, 100, 200, 600]

for j,i in enumerate(draw_samples):
    plt.subplot(5, 2, j+1)
    bayesian_strat.sample_bandits(i)
    plot_priors(bayesian_strat, hidden_prob)
    #plt.legend()
    plt.autoscale(tight = True)
plt.tight_layout()

plt.show()

请注意，我们并不是真正关心对隐含概率的精确估计，这点和机器学习中的回归预测是不同的。

我们更感兴趣的是选择最好的游戏机，或者更准确地说，更有信心地选择最好的游戏机。

出于这样的原因，红色游戏机的分布很宽，这代表了我们对隐含概率所知甚少，即从样本数据中提取到的信息有限，或者说样本数据给我们的先验带来的熵减很小。我们有充足的理由相信，红色游戏机不是最好的，所以选择忽略它。

另一方面，经过1000轮之后，大多数蓝色游戏机遥遥领先，因此我们几乎总是选择这台游戏机。这是一件好事，因为它经常能带来较好的回报。

4. 定义损失函数

上一小节我们得到所有游戏机的后验概率分布，也大致知道了该如何选游戏机，但这是不够的。我们的目标不是玩数字游戏，我们的目标是确确实实地给出一个可以落地执行的游戏机选择策略，类似这样的，【蓝色，蓝色，绿色，蓝色.....蓝色】这种序列。

要回答这个问题，就需要在概率分布和实际问题之间搭起一个桥梁，即损失函数，通过损失函数的数值化评估来得出最佳的后验策略。

我们需要一个指标来计算我们做的如何。理论上说，绝对最好的方法是始终挑那个获胜概率最大的游戏机。

记这台最好的游戏机的赢的概率为Wopt，我们可以定义一个理论的总遗憾，表示如果从一开始就选择最好的游戏机，和我们每轮实际选择的游戏机，这两种选择之间在收益上的差距。

在此公式中，Wb(i) 是所选游戏机在第 i 轮出奖的概率。

很显然，总遗憾为0意味着该策略获得最好的成绩，但这几乎是不太可能的，因为一开始我们的算法往往会做出错误或者不那么好的选择，只是随着轮数的增加，算法做出正确选择的概率逐渐增大。

理想情况下，总遗憾应该扁平化，因为它逐渐学习到最好的游戏机，即找到最好的后验分布对应的游戏机，这意味着我们常常能收敛到 Wb(i) = Wopt。

5. 选择游戏策略

我们已经可以得到不同游戏机的后验分布，这可以作为每轮选择的参考，但具体怎么选，依赖于我们选择的游戏策略。同时我们也有了损失评估函数，可以实时地看到每一轮选择后的损失。

在下面的代码中，我们对比了在不同的游戏策略下，总遗憾的函数曲线：

1. 随机：顾名思义，类似于丢色子，这显然不明智，如果用随机策略，就没必要费那么大劲去统计样本以及计算游戏机的后验分布了；
2. 贝叶斯的最大置信边界：选择底层概率的95%置信区间的最大上界的游戏机；
3. 贝叶斯-UCB算法：选择有最大得到的游戏机，其中得分是一个动态的后验分布的分位数；
4. 后验均值：选择具有最大后验均值的游戏机；
5. 最大比例：选择目前观测到的赢的比例最大的游戏机；

代码：

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats
import pymc as pm

rand = np.random.rand

figsize(12.5, 5)
from other_strats import *

#define a harder problem
hidden_prob = np.array([0.15, 0.2, 0.1, 0.05])
bandits = Bandits(hidden_prob)

#define regret
def regret(probabilities, choices):
    w_opt = probabilities.max()
    return (w_opt - probabilities[choices.astype(int)]).cumsum()

#create new strategies
strategies= [upper_credible_choice,
            bayesian_bandit_choice,
            ucb_bayes ,
            max_mean,
            random_choice]
algos = []
for strat in strategies:
    algos.append(GeneralBanditStrat(bandits, strat))

# train 10000 times
for strat in algos:
    strat.sample_bandits(10000)

# test and plot
for i, strat in enumerate(algos):
    _regret = regret(hidden_prob, strat.choices)
    plt.plot(_regret, label=strategies[i].__name__, lw=3)

plt.title("Total Regret of Bayesian Bandits Strategy vs. Random guessing")
plt.xlabel("Number of pulls")
plt.ylabel("Regret after $n$ pulls");
plt.legend(loc="upper left")

plt.show() 

从上图中可以看到，除了随机和后验均值策略之外，其他策略的总遗憾是逐渐收敛的，这表示了我们正在实现较优的选择。

6. 评估总遗憾期望

上个小节中，我们已经得到了3种游戏策略，都在总遗憾函数上呈现出了收敛的趋势。但是为了更科学，以消除任何可能的运气成分，我们应该看一下总遗憾期望。它定义为所有可能场景的总遗憾的期望值：

可以证明，任何次优策略的总遗憾期望都有对数形式的下界（从原始函数对数形式收敛来理解）。形式为：

因此，任何符合对数增加遗憾的策略，都可以称之为解决了多臂游戏机问题。

使用大数定理，我们可以通过进行很多次同样的实验来近似贝叶斯游戏机的总遗憾期望。

为了对不同策略间的差异性有一个更好的比较，我们在对数尺度中绘制了函数图：

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats
import pymc as pm

rand = np.random.rand

figsize(12.5, 5)
from other_strats import *

#define a harder problem
hidden_prob = np.array([0.15, 0.2, 0.1, 0.05])
bandits = Bandits(hidden_prob)

#define regret
def regret(probabilities, choices):
    w_opt = probabilities.max()
    return (w_opt - probabilities[choices.astype(int)]).cumsum()

#create new strategies
strategies= [upper_credible_choice,
            bayesian_bandit_choice,
            ucb_bayes ,
            max_mean,
            random_choice]

trials = 500
expected_total_regret = np.zeros((10000, 3))

for i_strat, strat in enumerate(strategies[:-2]):
    for i in range(trials):
        general_strat = GeneralBanditStrat(bandits, strat)
        general_strat.sample_bandits(10000)
        _regret = regret(hidden_prob, general_strat.choices)
        expected_total_regret[:, i_strat] += _regret
    plt.plot(expected_total_regret[:, i_strat] / trials, lw=3, label=strat.__name__)


[pl1, pl2, pl3] = plt.plot(expected_total_regret[:, [0,1,2]], lw = 3)
plt.xscale("log")
plt.legend([pl1, pl2, pl3],
           ["Upper Credible Bound", "Bayesian Bandit", "UCB-Bayes"],
            loc="upper left")
plt.ylabel("Exepected Total Regret \n after $\log{n}$ pulls");
plt.title( "log-scale of above" );
plt.ylabel("Exepected Total Regret \n after $\log{n}$ pulls");

plt.show()

0x4：算法扩展

1. 添加学习速率

我们能够通过加入一个学习速率项（就像深度学习中那样），促进该算法更快地更新去学习变化的环境

1. 如果rate<1，则该算法将更快地忘记先前的获胜，并且会有一个走向无知的下行压力；
2. 如果rate>1，则意味着算法将以风险较高的方式运行，而且更经常地把赌注压在早期赢的游戏机上，对不断变化的环境更有韧性；

2. 层次算法

我们可以在较小的游戏机算法之上再建立一个贝叶斯游戏机算法。即再建立一个贝叶斯游戏机模型，用于选择选择哪个子模型。原理上类似决策树和随机森林的概念。