1. 写在之前的话
0x1：贝叶斯推断的思想
我们从一个例子开始我们本文的讨论。小明是一个编程老手，但是依然坚信bug仍有可能在代码中存在。于是，在实现了一段特别难的算法之后，他开始决定先来一个简单的测试用例，这个用例通过了。接着，他用了一个稍微复杂的测试用例，再次通过了。接下来更难的测试用例也通过了，这时，小明开始觉得这段代码出现bug的可能性大大大大降低了....
上面这段白话文中，已经包含了最质朴的贝叶斯思想了！简单来说，贝叶斯推断是通过新得到的证据不断地更新我们的信念。
贝叶斯推断很少会做出绝对的判断，但可以做出非常可信的判断。在上面的例子中，小明永远无法100%肯定代码是无缺陷的，除非他穷举测试了每一种可能出现的情形，这在实践中几乎是不可能的。但是，我们可以对代码进行大量的测试，如果每一次测试都通过了，我们“更”有把握觉得这段代码是没问题的。
贝叶斯推断的工作方式就在这里：我们会随着新的证据不断更新之前的信念，但很少做出绝对的判断，除非所有其他的可能都被一一排除。
0x2：贝叶斯思维和频率派思维的区别
1. 频率派关于概率的解释 - 统计
频率派是更古典的统计学派，它们认为概率是事件在长时间内发生的概率。例如，在频率派的哲学语境里，飞机事故的概率指的是长期来看，飞机事故的概率值。
对许多事件来说，这样解释是有逻辑的。但是对某些没有长期频率的事件来说，这样的解释是难以理解的。例如，在总统选举时，某个候选人的获选的概率是难以用统计的方式来回答的，因为这个候选人的选举可能只发生一次。频率论为了解决这个问题，提出了“替代现实”的说法，从而用在所有这些替代的“现实”中发生的概率定义了这个概率。
2. 贝叶斯派关于概率的解释 - 信息
贝叶斯派对概率有更直观的解释。它把概率解释成是对事件发生的信心。简单地说，概率是观点的概述。
如果某人把概率 0 赋给某个事件的发生，表明他完全确定此事不会发生；相反，如果某人赋的概率值是 1，则表明他十分肯定此事一定会发生。
0 和 1 之间的概率值可以表明此事可能发生的权重。
这个概率定义可以和飞机事故概率一致，如果除去所有外部信息，一个人对飞机事故发生的信心应该等同于他了解到的飞机事故的频率。同样，贝叶斯概率定义也能适用于总统选举，并且使得概率（信心）更加有意义：即你对候选人 A 获胜的信息有多少？
对我们而言，将对一个事件发生的信息等同于概率是十分自然的，这正是我们长期以来同世界打交道的方式。我们只能了解到部分的真相，但可以通过不断收集证据来完善我们对事物的观念。
3. 贝叶斯推断的先验概率和后验概率
1）先验概率
对于统计学派来说，在进行统计之前，我们对概率是一无所知的。对于统计学派和统计方法来说，所有的概率都必须在完成统计之后（包括最大似然估计）才能得到一个统计概率。
但是对于贝叶斯推断来说，在进行统计之前，我们对事物已经拥有了一个先验的判断，即先验概率，记为 P(A)。
例如：P(A)：硬币有50%的概率是正面；P(A)：一段很长的很复杂的代码可能含有bug
2）后验概率
我们用 P(A|X) 表示更新后的信念，其含义为在得到证据 X 后，A 事件的概率。考虑下面的例子：
P(A|X)：我们观察到硬币落地后是正面，那么硬币是正面的后验概率会比50%大；
P(A|X)：代码通过了所有的 X 个测试，现在代码可能仍然有bug，但是现在这个概率变得非常小了。
3）贝叶斯推断柔和地融合了我们对事物的先验领域知识以及实际的观察样本的统计结果
在上述的例子中，即便获得了新的证据，我们也并没有完全地放弃初始的信念，但是我们重新调整了信念使之符合目前的证据（训练样本），也就是说，证据让我们对某些结果更有信息。
通过引入先验的不确定性，我们事实上允许了我们的初始信念可能是错误的。在观察数据、证据或其他信息之后，我们不断更新我们的信念使得它错得不那么离谱
0x3：在大数据情况下统计思想和贝叶斯推断思想是相等的？
频率推断和贝叶斯推断都是一种编程函数，输入的是各种统计问题。频率推断返回的是一个统计值（通常是统计量，平均值是一个典型的例子），而贝叶斯推断则会返回一个概率值。
这个小节我们来讨论下，当样本数量 N 不断增大时，频率推断和贝叶斯推断会发生什么。
1. 频率推断：N越大，频率统计越接近真实概率
频率统计学派假设世界上存在一个真理，即概率。统计做的事只是通过大量的实验来揭示这个真理本身而已。当 N 趋于无穷时，频率统计得到的统计值会无限接近于概率本身，即规律本身。
2. 贝叶斯推断：N越大，先验的离谱程度会不断降低，即越靠近真实的规律本身
当我们增加 N 时，即增加更多的证据时，初始的信念会不断地被“洗刷”。这是符合期望的。例如如果你的先验是非常荒谬的信念类似“太阳今天会爆炸”，那么每一天的结果都会不断“洗刷”这个荒谬的先验，使得后验概率不断下降。
3. N趋于无穷大时，频率推断和贝叶斯推断趋近相等
让 N 表示我们拥有的证据的数量。如果 N 趋于无穷大，那么贝叶斯的结果通常和频率派的结果一致。
因此，对于足够大的 N，统计推断多多少少还是比较客观的。
另一方面，对于较小的 N，统计推断相对而言不太稳定，频率统计的结果有更大的方差和置信区间。所以才会有相关性检验、卡方检验、皮森相关系数这些评价算法。
贝叶斯推断在小数据这方面效果要更好，通过引入先验概率和返回概率结果（而不是某个固定的值），贝叶斯推断保留了不确定性，这种不确定性正式小数据集的不稳定性的反映。
4. 频率统计仍然是十分有用的！
虽然频率统计的有更大的方差和置信区间，在小数据集下表现也不一定稳定。但是频率统计仍然非常有用，在很多领域可能都是最好的办法。例如最小方差线性回归、LASSO回归、EM算法，这些都是非常有效和快速的方法。
贝叶斯方法能够在这些方法失效的场景发挥作用，或者是作为更有弹性的模型而补充。
0x4：贝叶斯推断和极大似然估计、SGD随机梯度下降是什么关系？
贝叶斯推断和极大似然估计、SGD随机梯度下降一样，都是一种求解算法模型超参数（模型分布）的方法。
贝叶斯推断更侧重于先验和证据的动态融合，得到极大后验概率估计；
极大似然估计更侧重于对训练样本的频率统计；
而SGD在各种场景下都能表现出较好的性能，通过反向梯度的思想，不断进行局部最优的求解，从而获取近似的全局最优解。
2. 贝叶斯框架
我们知道，对于我们感兴趣的估计，可以通过贝叶斯思想被解释为概率。在频率统计中，这个估计时通过统计得到的，而在贝叶斯估计中，估计时通过后验概率得到的。
我们已经从抽象概念上了解了何谓后验概率，后验概率即为：新证据对初始信念洗刷后的结果，即修正后的信念。这个概念怎么落实到数学计算上呢？
0x1：贝叶斯定理
我们需要一个方法能够将几个因素连接起来，即先验概率、我们的证据、后验概率。这个方法就是贝叶斯定理，得名于它的创立者托马斯.贝叶斯。
上面的公式并不等同于贝叶斯推断，它是一个存在于贝叶斯推断之外的数学真理。或者应该说贝叶斯推断是在贝叶斯定理的基石上建立起来的。
在贝叶斯推断里，贝叶斯定理公式仅仅被用来连接先验概率 P(A) 和更新后的后验概率 P(A|B)
0x2：一个贝叶斯推断的例子 - 图书管理员还是农民
斯蒂文是一个害羞的人，他乐于助人，但是他对其他人不太关注。他非常乐见事情处于合理的顺序，并对他的工作非常细心。从特征工程的角度来看，似乎我们更应该判定史蒂文为已个图书管理员。
但是这里忽略了一个非常重要的先验知识，即男性图书管理员的人数只有男性农民的 1/20，这个先验知识将极大改变史蒂文是否是图书管理员这个事件的后验概率。
1. 形式化定义这个问题
把问题简化，假设世上只有两种职业，图书管理员和农民，并且农民的数量确实是图书管理员的20倍。
设事件 A 为史蒂文是一个图书管理员。如果我们没有史蒂文的任何信息，那么 P(A) = 1/21。这是我们的先验。
我们假设从史蒂文的邻居们那里我们获得了关于他的一些信息，我们称它们为 X，我们想知道的就是 P(A|X)。由贝叶斯定理得：P(A|X) = P(X|A) * P(A) / P(X)。
P(X|A)被定义为在史蒂文真的是一个图书管理员的情况下，史蒂文的邻居们给出的某种描述的概率，即如果史蒂文真的是一个图书管理员，他的邻居们将他描述为一个图书管理员的概率。这个值很可能接近于1，假设它为 0.95。
P(X) 可以解释为：任何人史蒂文的描述和史蒂文邻居的描述一致的概率，我们将其做逻辑上的改造：
P(X) = P(X and A) + P(X and ~A) = P(X|A)P(A) + P(X|~A)P(~A)
P(X|A)P(A) 这两个值我们都已经知道了。另外，~A 表示史蒂文不是一个图书管理员的概率，那么他一定是一个农民。所以 P(~A) = 1 - P(A) = 20 / 21。
P(X|~A) 即在史蒂文是一个农名的情况下，史蒂文的邻居给出的某个描述的概率，我们假设其为 0.5。
综上，P(X) = 0.95 * 1/24 + 0.5 * 20/21 = 0.515773809524
带入贝叶斯定理公式，P(A|X) = 0.95 * 1/20 / 0.515773809524 = 0.0877091748413
可以看到，因为农民数量远大于图书管理员，所以后验概率并不算高，尽管我们看到 P(X|A)的概率要大于 P(X|~A)，即我们得到的特征描述更倾向于表明史蒂文是一个图书管理员。
这个例子非常好的阐述了，先验概率以及证据时如何动态影响最终的后验概率结果的。
2. 通过可视化的代码呈现该问题
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt plt.rcParams['savefig.dpi'] = 300 plt.rcParams['figure.dpi'] = 300 colours = ["#348ABD", "#A60628"] prior = [1/20., 20/21.] # 图书管理员和农民的先验概率分别是 1/20 和 20/21 posterior = [0.087, 1-0.087] # 后验概率 plt.bar([0, .7], prior, alpha=0.70, width=0.25, color=colours[0], label="prior distribution", lw="3", edgecolor=colours[0]) plt.bar([0 + 0.25, .7 + 0.25], posterior, alpha=0.7, width=0.25, color=colours[1], label="posterior distribution", lw="3", edgecolor=colours[1]) plt.ylim(0, 1) plt.xticks([0.20, 0.95], ['Librarian', 'Farmer']) plt.title("Prior and Posterior probability of Steve's occupation") plt.ylabel("Probability") plt.legend(loc="upper left") plt.show()
可以看到，在得到 X 的观测值之后，史蒂文为图书管理员的后台概率增加了。同时，农民的后验概率降低了。
但同时也看到，史蒂文为图书管理员的后验概率增加的并不是很多，史蒂文是农民的可能性依旧很大。
这表明，先验概率对贝叶斯推断影响是很大的，例如我们假设一台服务器平时不会被黑客攻陷（IOC），但是如果观测到一系列的异常事件，被攻陷的后验概率就会不断增大，但这取决于观测到的现象的数量和强度。
3. 一个从短信数据推断行为模式的例子 - 使用计算机执行贝叶斯推断
0x1：场景描述 - 我们的数据是什么样的？
下面的直方图显示了一段时间内的短信数量的统计。
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt count_data = np.loadtxt("data/txtdata.csv") n_count_data = len(count_data) plt.bar(np.arange(n_count_data), count_data, color="#348ABD") plt.xlabel("Time (days)") plt.ylabel("count of text-msgs received") plt.title("Did the user's texting habits change over time?") plt.xlim(0, n_count_data) plt.show()
数据集下载链接。
0x2：问题定义 - 我们要从数据中挖掘出什么东西？
我们的问题是：这段时间内用户的行为出现了什么变化吗？或者是说，这段时间内是否出现了异常的行为模式，又具体是几号发生的？
这特别像安全攻防领域经常面对的一个问题，即在一段日志中进行数据挖掘，从中找出是否存在异常行为模式（即入侵IOC）。
通常来说，一个直观的想法可能会是这样，我们通过ngram或者定长的时间窗口进行日志聚合，然后通过人工经验聚合中抽取一系列的特征，并希望借助有监督或者无监督的方法直接得到一个 0/1 判断，即这个窗口期间内是否存在异常点。这是一种典型的端到端判断场景。
但是其实我觉得应该更灵活地理解算法，不一定每一个场景都适合使用判别式模型进行二分或者多分判断，有时候，我们其实可以尝试使用生成式模型，基于训练数据，得到一个或多个概率分布生成式。而多个生成式模型之间的分界点，很可能就是代表了一种模式跃迁，即行为异常。
0x3：使用什么生成式模型进行数据拟合？
1. 二项分布？正态分布？泊松分布？
二项分布，正态分布，泊松分布这3个都是对自然界的事物在一定时间内的发生次数进行泛化模拟的分布模型。实际上，它们3个之间相差并不是很大。
如果忽略分布是离散还是连续的前提（二项分布和泊松分布一样都是离散型概率分布，正态分布是连续型概率分布），二项分布与泊松分布以及正态分布至少在形状上是十分接近的，也即两边低中部高。
当 n 足够大，p 足够小（事件间彼此独立，事件发生的概率不算太大，事件发生的概率是稳定的），二项分布逼近泊松分布，λ=np；
一个被广泛接受的经验法则是如果 n≥20，p≤0.05，可以用泊松分布来估计二项分布值。
至于正态分布是一个连续分布 当实验次数 n 再变大，几乎可以看成连续时二项分布和泊松分布都可以用正态分布来代替。
2. 单位时间内随机事件发生的次数趋近于一个相对固定值 - 选泊松分布
在日常生活中，大量事件我们可以预见其在时间窗口内的总数，但是无法准确知道具体某个时刻发生的次数。例如
1. 某医院平均每小时出生3个婴儿； 2. 某公司平均每10分钟接到1个电话； 3. 某超市平均每天销售4包 xx 牌奶粉； 4. 某网站平均每分钟有2次访问；
这些事件的特点就是，我们可以预估这些事件的总数，但是没法具体知道发生时间。例如医院具体下一小时的前10分钟出现几个婴儿，有可能前10分钟出现1个，也可能一下子出生3个。
泊松分布就是描述某段时间内，事件具体的发生概率。
P 表示概率，N 表示某种函数关系，t 表示时间，n 表示数量，λ 表示事件的频率。
从公式上看，我们得到几个泊松分布的基本判断：
1. 单位时间内发生 n 次数的概率随着时间 t 变化；
2. 概率受发生次数 n 的影响；
3. λ 相当于这件事的先验属性，在整个时间周期内都影响了最终的概率，λ 被称为此分布的一个超参数，它决定了这个分布的形式。
λ 可以位任意正数，随着 λ 的增加，得到大值的概率会增大。相反地，当 λ 减小时，得到小值的概率会增大。λ 可以被称为 Poisson分布的强度。
接下来两个小时，一个婴儿都不出现的概率是（已知 λ = 3）：
基本不可能发生。
接下来一个小时，至少出现两个婴儿的概率是：
即 80%。
泊松分布的图像大概是这样的：
可以看到，在频率附近，事件发生的概率最高，然后向两边对称下降。每小时出生 3 个婴儿，这是最可能的结果，出生得越多或越少，就越不可能。
Poisson分布的一个重要性质是：它的期望值等于它的参数。即：E[Z|λ] = λ。
下图展示了不同 λ 取值下的概率质量分布，可以看到，随着 λ 的增加，取到大值的概率会增加。
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt import scipy.stats as stats a = np.arange(16) poi = stats.poisson lambda_ = [1.5, 4.25] colours = ["#348ABD", "#A60628"] plt.bar(a, poi.pmf(a, lambda_[0]), color=colours[0], label="$\lambda = %.1f$" % lambda_[0], alpha=0.60, edgecolor=colours[0], lw="3") plt.bar(a, poi.pmf(a, lambda_[1]), color=colours[1], label="$\lambda = %.1f$" % lambda_[1], alpha=0.60, edgecolor=colours[1], lw="3") plt.xticks(a + 0.4, a) plt.legend() plt.ylabel("probability of $k$") plt.xlabel("$k$") plt.title("Probability mass function of a Poisson random variable; differing \ $\lambda$ values") plt.show()
3. 用单个泊松分布就足够拟合数据集了吗？
好！我们现在已经确定要使用泊松分布来对数据集进行拟合。接下来的问题是，我们是不是只要用单个泊松分布就足够对这份数据集进行完美拟合了呢？
我们再来看看我们的问题场景。
我们这里通过观察数据，我们认为在某个时间段，发生了一个行为模式的突变，即在后期，收到短信的频率好像提高了。
实际上，这已经进入到异常检测的范畴内了，我们可以基于泊松分布混合模型对这个数据集进行拟合。
4. 怎么形式化描述泊松分布参数 λ 的阶跃情况？
我们不能确定参数 λ 的真实取值，但是可以猜测在某个时段 λ 增加了（λ取大值时，更容易得到较大的结果值），也即一天收到的短信条数比较多的概率增大了。
我们怎么形式化地描述这个现象呢？假设在某一个时间段内(用τ表示)(τ可能也无法准确定义到具体的某一天)，参数?突然变大。
所以?可以取两个值，在时段τ之前有一个，在τ时段之后有另一个。这种突然发生变化的情况称之为转换点。
我们以贝叶斯推断的思想来看待这个问题：
1. 首先，存在一个 λ 阶跃点只是一个先验假设。事实上，可能并不存在这样一个 λ 阶跃点。 2. 第二，具体的阶跃点 τ 是什么时间呢？可能我们无法准确知道时间点，只能得到一个时间区间的概率分布。
1）确定泊松分布参数 λ 的先验分布
如果实际中 λ 没有变化，则 
，那么 λ 的先验概率是相等的。
在贝叶斯推断下，我们需要对不同可能取值的 λ 分配相应的先验概率。那具体什么样的先验概率对
来说才是最有可能的呢？
严格来说，我们对这件事几乎没有任何先验知识，因为我们可能对这个数据集的属主用户行为模式一无所知。
注意到 λ 可以取任意正数，所以我们可以用指数分布来表示 λ 的先验分布，因为指数分布对任意正数都存在一个连续密度函数，这对模拟 λ 来说是一个很好的选择。但是，指数分布本身也有它自己的参数，所以这个参数也需要包括在我们的模型公式里面。我们称它为 α。不同的 α 决定了该指数分布的形状。
α 称之为超参数。该参数能够影响其他参数。按照最大熵算法的原理，我们不希望对这个参数赋予太多的主观色彩。
注意到我们计算样本中计算平均值的公式：
，它正好就是我们要估计的超参数 α 的逆，所以，超参数 α 可以直接根据观测样本的计算得到。
2）确定阶跃点 τ 的先验分布
我们对于参数τ，由于受到噪声数据的影响，很难为它选择合适的先验。所以还是依照最大熵原理，我们假定每一天的先验估计都是一致的。即：
在不知道任何先验知识的情况下，等概论被证明是一种最好的先验假设。
0x4：使用PyMC进行贝叶斯推断
Cronin对概率编程有一段激动人心的描述：
换一种方式考虑这件事件，跟传统的编程仅仅向前运行不同的是，概率编程既可以向前也可以向后运行。它通过向前运行来计算其中包含的所有关于世界的假设结果（例如对模型空间的描述），但它也从数据中向后运行，以约束可能的解释。
在实践中，许多概率编程系统将这些向前和向后的操作巧妙地交织在一起，以给出有效的最佳解释。
1. 声明模型中的随机变量
按照前面我们的讨论，α 取样本中计算平均值的逆，而泊松分布的的参数 λ 由 α决定。对阶跃点先验估计的 τ 参数，使用等概率估计（最大熵原理）。
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt import pymc as pm count_data = np.loadtxt("data/txtdata.csv") n_count_data = len(count_data) plt.bar(np.arange(n_count_data), count_data, color="#348ABD") plt.xlabel("Time (days)") plt.ylabel("count of text-msgs received") plt.title("Did the user's texting habits change over time?") plt.xlim(0, n_count_data) alpha = 1.0 / count_data.mean() # Recall count_data is the # variable that holds our txt counts lambda_1 = pm.Exponential("lambda_1", alpha) lambda_2 = pm.Exponential("lambda_2", alpha) tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data) print("Random output:", tau.random(), tau.random(), tau.random())
代码输出的结果为：('tau Random output:', array(36), array(57), array(38))
表明我们的阶跃点 τ 参数是等概率随机生成的。
参数 λ 受参数 τ 影响，所以在阶跃点两边的不同的泊松分布参数 λ 也是随机的，它根据 τ 随机生成。
我们希望PyMC收集所有变量信息，以及观测数据，并从中产生一个model。PyMC可以借助一个叫马尔科夫链蒙特卡洛（MCMC）的算法。得到参数 λ1，λ2，τ 的后验估计
2. 利用MCMC根据观测数据得到对模型参数的极大后验估计结果
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt import pymc as pm count_data = np.loadtxt("data/txtdata.csv") n_count_data = len(count_data) plt.bar(np.arange(n_count_data), count_data, color="#348ABD") plt.xlabel("Time (days)") plt.ylabel("count of text-msgs received") plt.title("Did the user's texting habits change over time?") plt.xlim(0, n_count_data) alpha = 1.0 / count_data.mean() # Recall count_data is the # variable that holds our txt counts lambda_1 = pm.Exponential("lambda_1", alpha) lambda_2 = pm.Exponential("lambda_2", alpha) tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data) @pm.deterministic def lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2): out = np.zeros(n_count_data) out[:tau] = lambda_1 # lambda before tau is lambda1 out[tau:] = lambda_2 # lambda after (and including) tau is lambda2 return out print("tau Random output:", tau.random(), tau.random(), tau.random()) #print("lambda distribution:", lambda_(tau, lambda_1, lambda_2)) observation = pm.Poisson("obs", lambda_, value=count_data, observed=True) model = pm.Model([observation, lambda_1, lambda_2, tau]) mcmc = pm.MCMC(model) mcmc.sample(40000, 10000, 1) lambda_1_samples = mcmc.trace('lambda_1')[:] lambda_2_samples = mcmc.trace('lambda_2')[:] tau_samples = mcmc.trace('tau')[:] # histogram of the samples: ax = plt.subplot(311) ax.set_autoscaley_on(False) plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85, label="posterior of $\lambda_1$", color="#A60628", normed=True) plt.legend(loc="upper left") plt.title(r"""Posterior distributions of the variables $\lambda_1,\;\lambda_2,\;\tau$""") plt.xlim([15, 30]) plt.xlabel("$\lambda_1$ value") ax = plt.subplot(312) ax.set_autoscaley_on(False) plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.85, label="posterior of $\lambda_2$", color="#7A68A6", normed=True) plt.legend(loc="upper left") plt.xlim([15, 30]) plt.xlabel("$\lambda_2$ value") plt.subplot(313) w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples) plt.hist(tau_samples, bins=n_count_data, alpha=1, label=r"posterior of $\tau$", color="#467821", weights=w, rwidth=2.) plt.xticks(np.arange(n_count_data)) plt.legend(loc="upper left") plt.ylim([0, .75]) plt.xlim([35, len(count_data) - 20]) plt.xlabel(r"$\tau$ (in days)") plt.ylabel("probability") plt.show()
0x5：得到的模型参数后验分布结果说明了什么？
首先，我们可以看到，贝叶斯推断得到的是一个概率分布，不是一个准确的结果。我们使用分布描述位置的 λ 和 τ。
分布越广，我们的后验概率越小
1. 有很大概率存在异常行为模式的突变
我们也可以看到参数的合理值：λ1 大概是18，λ2 大概是23。这两个 λ 的后验分布明显是不同的，这表明该用户接受短信的行为确实发生了变化。
2. 后验分布是我们假设的指数分布吗？
我们还注意到，λ 的后验分布并不是我们在模型建立时假设的指数分布，事实上，λ 的后验分布并不是我们能从原始模型中可以辨别的任何分布。
但是这恰恰就是概率编程的好处所在，我们借助计算机绕开了复杂的数学建模过程。
3. 阶跃点时间点
我们的分析结果返回了一个 τ 的一个概率分布。在很多天中，τ 不太好确定，总体来说，我们可以认为在3~4天内可以认为是潜在的阶跃点。
我们可以看到，在45天左右，有50%的把握可以说用户的行为是有所改变的。这个点的概率分布突然出现了一个峰值，我们可以推测产生这种情况的原因是：短信费用的降低，天气提醒短信的订阅，或者是一段感情的开始。
0x6：根据模型参数得到后验样本
在开始讨论后验样本之前。首先需要明白的是，模型的后验样本并不是真实存在的事实，它更多地像是我们基于已经得到的概率分布得到的一个期望值。因为很显然，我们的模型不能很好地表达出真实物理中的随机分布以及噪音的影响。我们能得到的只能是一个对整体情况的概括。
我们用后验样本来回答一下下面这个问题：
在第 t （0 <= t <= 70）天中，期望收到多少条短信？
显然，这个问题要求我们给出一个概括统计，即求每天短信数量的期望。
Poisson分布的期望值等于它的参数 λ。因此，问题相当于：在时间 t 中，参数 λ 的期望值是多少？
# -*- coding: utf-8 -*-
import numpy as np
from matplotlib import pyplot as plt
import pymc as pm
count_data = np.loadtxt("data/txtdata.csv")
n_count_data = len(count_data)
plt.bar(np.arange(n_count_data), count_data, color="#348ABD")
plt.xlabel("Time (days)")
plt.ylabel("count of text-msgs received")
plt.title("Did the user's texting habits change over time?")
plt.xlim(0, n_count_data)
alpha = 1.0 / count_data.mean() # Recall count_data is the
# variable that holds our txt counts
lambda_1 = pm.Exponential("lambda_1", alpha)
lambda_2 = pm.Exponential("lambda_2", alpha)
tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data)
@pm.deterministic
def lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2):
out = np.zeros(n_count_data)
out[:tau] = lambda_1 # lambda before tau is lambda1
out[tau:] = lambda_2 # lambda after (and including) tau is lambda2
return out
print("tau Random output:", tau.random(), tau.random(), tau.random())
#print("lambda distribution:", lambda_(tau, lambda_1, lambda_2))
observation = pm.Poisson("obs", lambda_, value=count_data, observed=True)
model = pm.Model([observation, lambda_1, lambda_2, tau])
mcmc = pm.MCMC(model)
mcmc.sample(40000, 10000, 1)
lambda_1_samples = mcmc.trace('lambda_1')[:]
print "lambda_1_samples"
print lambda_1_samples
lambda_2_samples = mcmc.trace('lambda_2')[:]
print "lambda_2_samples"
print lambda_2_samples
tau_samples = mcmc.trace('tau')[:]
print "tau_samples"
print tau_samples
# tau_samples, lambda_1_samples, lambda_2_samples contain
# N samples from the corresponding posterior distribution
N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_count_data)
for day in range(0, n_count_data):
# ix is a bool index of all tau samples corresponding to
# the switchpoint occurring prior to value of 'day'
ix = day < tau_samples # 如果天数在转折点tau之前，取值lambda_1,否则取lambda_2,然后在取平均值
# Each posterior sample corresponds to a value for tau.
# for each day, that value of tau indicates whether we're "before"
# (in the lambda1 "regime") or
# "after" (in the lambda2 "regime") the switchpoint.
# by taking the posterior sample of lambda1/2 accordingly, we can average
# over all samples to get an expected value for lambda on that day.
# As explained, the "message count" random variable is Poisson distributed,
# and therefore lambda (the poisson parameter) is the expected value of
# "message count".
print "lambda_1_samples[ix]"
print lambda_1_samples[ix]
print "lambda_2_samples[~ix]"
print lambda_2_samples[~ix]
expected_texts_per_day[day] = (lambda_1_samples[ix].sum() # 这个期望值我们假设是服从泊松分布的，分布的期望值＝参数lambda。
+ lambda_2_samples[~ix].sum()) / N
plt.plot(range(n_count_data), expected_texts_per_day, lw=4, color="#E24A33",
label="expected number of text-messages received")
plt.xlim(0, n_count_data)
plt.xlabel("Day")
plt.ylabel("Expected # text-messages")
plt.title("Expected number of text-messages received")
plt.ylim(0, 60)
plt.bar(np.arange(len(count_data)), count_data, color="#348ABD", alpha=0.65,
label="observed texts per day")
plt.legend(loc="upper left");
plt.show()
在代码中可以看到，给定某天 t，我们对所有可能的 λ 求均值。
1. 基于后验样本，从统计学上确定两个 λ 是否真的不一样
我们基于参数 λ 的后验分布得到一个结论，λ 存在阶跃点。但是如果这不是真相呢？有一部分分布其实是重合的呢？这个可能性有多大呢？怎么量化地评估这个可能性呢？
一个方法就是计算出 P（λ1 < λ2 | data），即在获得观察数据的前提下，λ1 的真实值比 λ2 小的概率。
print (lambda_1_samples < lambda_2_samples).mean() 0.999933333333
很显然，有99.93%的概率说这两个值是不相等的。
可以看到，这本质上是计算在 λ1 和 λ2 后验概率分布的重合程度。
0x7：扩充至两个阶跃点 - 能否借助贝叶斯推理自动挖掘潜在的规律？
我们在章节的一开始就作出一个假设，即认为我们的数据集中蕴含了一个现象，即用户的行为模式在某个时间点发生了突变。
但是是否存在可能说这个转折点不只一个呢？
假设现在我们对一个转折点表示很怀疑，我们现在认为用户的行为发生了两次改变．扩充之后，用户的行为分为三个阶段，三个泊松分布对应三个λ1,λ2,λ3，两个转折点τ1,τ2。
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt import pymc as pm count_data = np.loadtxt("data/txtdata.csv") n_count_data = len(count_data) plt.bar(np.arange(n_count_data), count_data, color="#348ABD") plt.xlabel("Time (days)") plt.ylabel("count of text-msgs received") plt.title("Did the user's texting habits change over time?") plt.xlim(0, n_count_data) alpha = 1.0 / count_data.mean() # Recall count_data is the # variable that holds our txt counts lambda_1 = pm.Exponential("lambda_1", alpha) lambda_2 = pm.Exponential("lambda_2", alpha) lambda_3 = pm.Exponential("lambda_3", alpha) tau_1 = pm.DiscreteUniform("tau_1", lower=0, upper=n_count_data-1) tau_2 = pm.DiscreteUniform("tau_2", lower=tau_1, upper=n_count_data) @pm.deterministic def lambda_(tau_1=tau_1, tau_2=tau_2, lambda_1=lambda_1, lambda_2=lambda_2, lambda_3=lambda_3): out = np.zeros(n_count_data) out[:tau_1] = lambda_1 # lambda before tau is lambda1 out[tau_1:tau_2] = lambda_2 # lambda after (and including) tau is lambda2 out[tau_2:] = lambda_3 # lambda after (and including) tau is lambda2 return out observation = pm.Poisson("obs", lambda_, value=count_data, observed=True) model = pm.Model([observation, lambda_1, lambda_2, lambda_3, tau_1, tau_2]) mcmc = pm.MCMC(model) mcmc.sample(40000, 10000) lambda_1_samples = mcmc.trace('lambda_1')[:] lambda_2_samples = mcmc.trace('lambda_2')[:] lambda_3_samples = mcmc.trace('lambda_3')[:] tau_1_samples = mcmc.trace('tau_1')[:] tau_2_samples = mcmc.trace('tau_2')[:] # histogram of the samples: # lambda_1 ax = plt.subplot(311) ax.set_autoscaley_on(False) plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, label="$\lambda_1$", normed=True) plt.legend(loc="upper left") plt.grid(True) plt.title(r"""Posterior distributions of the variables $\lambda_1,\;\lambda_2,\;\tau$""") # x轴坐标范围 plt.xlim([15, 30]) plt.xlabel("$\lambda_1$ value") # lambda_2 ax = plt.subplot(312) ax.set_autoscaley_on(False) plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, label=" $\lambda_2$",color='#3009A6',normed=True) plt.legend(loc="upper left") plt.grid(True) plt.xlim([30, 90]) plt.xlabel("$\lambda_2$ value") # lambda_3 ax = plt.subplot(313) ax.set_autoscaley_on(False) plt.hist(lambda_3_samples, histtype='stepfilled', bins=30, label=" $\lambda_2$",color='#6A63A6',normed=True) plt.legend(loc="upper left") plt.grid(True) plt.xlim([15, 30]) plt.xlabel("$\lambda_3$ value") plt.show()
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt import pymc as pm count_data = np.loadtxt("data/txtdata.csv") n_count_data = len(count_data) plt.bar(np.arange(n_count_data), count_data, color="#348ABD") plt.xlabel("Time (days)") plt.ylabel("count of text-msgs received") plt.title("Did the user's texting habits change over time?") plt.xlim(0, n_count_data) alpha = 1.0 / count_data.mean() # Recall count_data is the # variable that holds our txt counts lambda_1 = pm.Exponential("lambda_1", alpha) lambda_2 = pm.Exponential("lambda_2", alpha) lambda_3 = pm.Exponential("lambda_3", alpha) tau_1 = pm.DiscreteUniform("tau_1", lower=0, upper=n_count_data-1) tau_2 = pm.DiscreteUniform("tau_2", lower=tau_1, upper=n_count_data) @pm.deterministic def lambda_(tau_1=tau_1, tau_2=tau_2, lambda_1=lambda_1, lambda_2=lambda_2, lambda_3=lambda_3): out = np.zeros(n_count_data) out[:tau_1] = lambda_1 # lambda before tau is lambda1 out[tau_1:tau_2] = lambda_2 # lambda after (and including) tau is lambda2 out[tau_2:] = lambda_3 # lambda after (and including) tau is lambda2 return out observation = pm.Poisson("obs", lambda_, value=count_data, observed=True) model = pm.Model([observation, lambda_1, lambda_2, lambda_3, tau_1, tau_2]) mcmc = pm.MCMC(model) mcmc.sample(40000, 10000) lambda_1_samples = mcmc.trace('lambda_1')[:] lambda_2_samples = mcmc.trace('lambda_2')[:] lambda_3_samples = mcmc.trace('lambda_3')[:] tau_1_samples = mcmc.trace('tau_1')[:] tau_2_samples = mcmc.trace('tau_2')[:] # histogram of the samples: # lambda_1 # tau_1 w = 1.0 / tau_1_samples.shape[0] * np.ones_like(tau_1_samples) plt.hist(tau_1_samples, bins=count_data, alpha=1, label=r"$\tau_1$",color="blue",weights=w ) plt.grid(True) plt.legend(loc="upper left") plt.xlabel(r"$\tau_1$ (in days)") plt.ylabel("probability") # tau_2 w = 1.0 / tau_2_samples.shape[0] * np.ones_like(tau_1_samples) plt.hist(tau_2_samples, bins=count_data, alpha=1, label=r"$\tau_2$",weights=w,color="red",) plt.xticks(np.arange(count_data)) plt.grid(True) plt.legend(loc="upper left") plt.ylim([0, 1.0]) plt.xlim([35, len(count_data) - 20]) plt.xlabel(r"$\tau_2$ (in days)") plt.ylabel("probability") plt.show()
贝叶斯推理展示了5个未知数的后验推理结果。我们可以看到模型的转折点大致在第45天和第47天的时候取得。这个是真实的情况吗？我们的模型有可能过拟合了吗？
确实，我们都可能对数据中有多少个转折点抱有怀疑的态度。这意味着应该有多少个转折点可以通过设置一个先验分布，并让模型自己做决定。
这是一种非常重要的思想：我们对所有未知的变量都可以用一个先验分布来假设它，通过观测样本的训练，得到最佳的后验分布。
Relevant Link:
https://blog.csdn.net/u014281392/article/details/79330286
4. 挑战者号事故案例
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] # plot it, as a function of temperature (the first column) print("Temp (F), O-Ring failure?") print(challenger_data) plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k", alpha=0.5) plt.yticks([0, 1]) plt.ylabel("Damage Incident?") plt.xlabel("Outside temperature (Fahrenheit)") plt.title("Defects of the Space Shuttle O-Rings vs temperature") plt.show()
实验数据如下：
Date,Temperature,Damage Incident 04/12/1981,66,0 11/12/1981,70,1 3/22/82,69,0 6/27/82,80,NA 01/11/1982,68,0 04/04/1983,67,0 6/18/83,72,0 8/30/83,73,0 11/28/83,70,0 02/03/1984,57,1 04/06/1984,63,1 8/30/84,70,1 10/05/1984,78,0 11/08/1984,67,0 1/24/85,53,1 04/12/1985,67,0 4/29/85,75,0 6/17/85,70,0 7/29/85,81,0 8/27/85,76,0 10/03/1985,79,0 10/30/85,75,1 11/26/85,76,0 01/12/1986,58,1 1/28/86,31,Challenger Accident
0x1：我们如何对数据分布进行建模？
1. 通过非线性回归模拟样本数据的函数表示
首先，通过简单的观察，这个样本数据的自变量和因变量明显不符合一个线性函数的关系。而且，我们甚至有理由怀疑它们之间存在一个阶跃函数的关系，即在65°时发生阶跃。但阶跃函数太过武断，不够柔和，而且也没有证据。
我们使用逻辑斯蒂线性回归，尝试对观测样本进行拟合。
# -*- coding: utf-8 -*- # Code source: Gael Varoquaux # License: BSD 3 clause import numpy as np import matplotlib.pyplot as plt from sklearn import linear_model # data set X = np.array([ [66.], [70.], [69.], [80.], [68.], [67.], [72.], [73.], [70.], [57.], [63.], [70.], [78.], [67.], [53.], [67.], [75.], [70.], [81.], [76.], [79.], [75.], [76.], [58.] ]) y = np.array([ 0., 1., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1., 0., 0., 1., 0., 0., 0., 0., 0., 0., 1., 0., 1. ]) print "X" print X.ravel() print "y" print y # and plot the data plt.figure(1, figsize=(8, 6)) plt.clf() plt.scatter(X.ravel(), y, color='black', zorder=20) X_test = np.linspace(50, 85, 300) # run the classifier clf = linear_model.LogisticRegression(C=1e5) clf.fit(X, y) def model(x): y_value = 1 / (1 + np.exp(-x)) return y_value # X_test * clf.coef_ + clf.intercept = 相当于: wx + b # model()负责转换为逻辑斯蒂回归函数的 y 值 print "clf.coef: ", clf.coef_ print "clf.intercept_: ", clf.intercept_ y_value = model(X_test * clf.coef_ + clf.intercept_).ravel() plt.plot(X_test, y_value, color='red', linewidth=3) # 画出模型拟合后的曲线 plt.ylabel('y') plt.xlabel('X') plt.xticks(range(50, 85)) plt.yticks([0, 0.5, 1]) plt.ylim(-.25, 1.25) plt.xlim(50, 85) plt.legend(('Logistic Regression Model'), loc="lower right", fontsize='small') plt.show()
我们得到了sklearn对样本数据的逻辑斯蒂回归拟合结果，我们可以看到，一个很明显的结论是：随着温度下降，出现故障的风险会不断提高，在50°时甚至接近于1，即必然发生。
同时我们也打印出了拟合的模型参数：
clf.coef: [[-0.22842526]] clf.intercept_: [ 14.77684766]
这两个参数决定这个逻辑函数的分布形态，我们后面还会看到，我们用贝叶斯推断的方式，也能获取到这2个参数的估计结果。
注意，这里coef是包含负号的，所以我们可以认为sklearn推导出来的结果是：1 / 1 + e ^ 0.2284 + 14.77
2. 通过数学建模和贝叶斯推理的方式进行数据拟合
1）温度和事故发生的概率如何建模？
? 是这个模型的参数，是我们需要从数据中进行推理得到的。下面画出? = 1,3, −5时的?(?)的曲线
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt def logistic(x, beta): return 1.0 / (1.0 + np.exp(beta * x)) x = np.linspace(-4, 4, 100) plt.plot(x, logistic(x, 1), label=r"$\beta = 1$") plt.plot(x, logistic(x, 3), label=r"$\beta = 3$") plt.plot(x, logistic(x, -5), label=r"$\beta = -5$") plt.legend() plt.show()
画这个图的意义是什么呢？其实没啥意义，这只是随便估计了一下参数，但起码我们看到，参数 ? 起码必须是正数，而且 ? 应该是一个倾向于小的值。当然，这只是我们的猜测，具体的结果要进行贝叶斯编程后推导得到。
此外，在上图中，在 ? = 0 附件概率会发生改变，但是在实际的数据中在 65 到 70 这个范围内概率发生改变。所以还需要在 logistic 函数中添加一个偏置项，函数变为：
不同的 ? 对应的 logistic 函数如下
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) x = np.linspace(-4, 4, 100) plt.plot(x, logistic(x, 1), label=r"$\beta = 1$", ls="--", lw=1) plt.plot(x, logistic(x, 3), label=r"$\beta = 3$", ls="--", lw=1) plt.plot(x, logistic(x, -5), label=r"$\beta = -5$", ls="--", lw=1) plt.plot(x, logistic(x, 1, 1), label=r"$\beta = 1, \alpha = 1$",color="#348ABD") plt.plot(x, logistic(x, 3, -2), label=r"$\beta = 3, \alpha = -2$",color="#A60628") plt.plot(x, logistic(x, -5, 7), label=r"$\beta = -5, \alpha = 7$",color="#7A68A6") plt.legend(loc="lower left"); plt.show()
? 参数控制着曲线在横轴上左右移动，这也就是把 ? 称为偏置的原因。
2）逻辑函数的参数如何建模？
读者小伙伴们注意，在贝叶斯概率编程领域，我们要习惯用概率分布的思维方式来思考问题。对于逻辑函数的参数 ? 和参数 ? 来说，并没有限定要为正数，或者在某个范围内要相对的大。基本上来说，我们不知道这两个参数的先验，因此我们用正态分布来模拟它们是一个相对比较合适的选择。
，正态分布有两个参数均值?，方差?。也可以用?2代替1⁄?，它们二者互为倒数。?越小，正态分布的曲线越平坦(不确定性越大)，反之亦然(不确定性越小)。?的取值范围为正数。
对模型参数建模之后，就可以代入上一小节温度和事故发生概率的模型中。
import pymc as pm temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha))
3）我们的观测数据如何和模型参数融合建模？
前两个小节，逻辑函数和正态分布都是针对单个样本的输入和输出进行建模，我们的实验观察是多个样本，那么我们怎么将我们的观测数据和我们的概率值（通过逻辑函数模拟的）联系起来呢？
如果按照我们在机器学习编程中的做法，这个时候，我们已经对输入和输出进行了建模（逻辑函数），那接下来的做法就是选择一个损失函数（目标函数）（这里很可能就是0-1损失），通过计算观察样本的整体损失并通过极值求导的方式来得到模型的最优拟合参数。
但是记住，我们现在是在贝叶斯概率编程领域中，贝叶斯推断理论认为，输入变量和输出变量之间是不存在绝对的因果推导关系的，我们所谓的观测结果只是在在一个概率分布下的一次（一批）具体实验具现化结果，我们对任何的因果关系都要持怀疑态度，即采集概率分布的方式来进行建模。
回到我们这个具体的问题场景，在一定的温度下，我们可以得到事件发生的概率，但是事件是否真的发生，有两种可能，发生 or 不发生。因此我们可以用 伯努利分布来对“事件发生概率 -> 是否真的发生”进行概率建模。
，其中 p(t) 是我们的逻辑函数，ti 是观测到的温度值。
4）通过pymc对观测样本进行模型参数的后验估计
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] # plot it, as a function of temperature (the first column) print("Temp (F), O-Ring failure?") print(challenger_data) plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k", alpha=0.5) plt.yticks([0, 1]) plt.ylabel("Damage Incident?") plt.xlabel("Outside temperature (Fahrenheit)") plt.title("Defects of the Space Shuttle O-Rings vs temperature") temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] # histogram of the samples: plt.subplot(211) plt.title(r"Posterior distributions of the variables $\alpha, \beta$") plt.hist(beta_samples, histtype='stepfilled', bins=35, alpha=0.85, label=r"posterior of $\beta$", color="#7A68A6", normed=True) plt.legend() plt.subplot(212) plt.hist(alpha_samples, histtype='stepfilled', bins=35, alpha=0.85, label=r"posterior of $\alpha$", color="#A60628", normed=True) plt.legend() plt.show()
0x2：我们从模型参数的后验概率估计中可以得到什么信息？
1. 温度对事故的发生肯定是有影响的
所有的 ? 样本值都大于 0，说明温度这一自变量是因变量，一定有所影响。
如果后验概率集中在 0 附件，假设? = 0，说明温度对失败的概率没有影响。
2. ? 一定是小于 0 的
所有的 ? 后验值都是负的且远离 0，可以确信 ? 一定是小于 0 的。
3. 我们对模型参数的真实取值依然无法 100% 确定
鉴于数据的分布情况，我们不确定真实的参数?, ?到底是什么。当然，考虑到样本集比较小，以及发生缺陷和没有缺陷事故的大量重叠，得到这样的结果也是正常的。
0x3：基于模型参数的后验估计反推事件发生的概率的后验估计
像之前的短信行为后验估计一样，我们现在得到了逻辑函数的模型参数的后验估计，我们现在可以反过来做，把估计得到的结果当做输入，得到模型的输出结果，即事件发生的后验概率。
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] # plot it, as a function of temperature (the first column) print("Temp (F), O-Ring failure?") print(challenger_data) plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k", alpha=0.5) plt.yticks([0, 1]) plt.ylabel("Damage Incident?") plt.xlabel("Outside temperature (Fahrenheit)") plt.title("Defects of the Space Shuttle O-Rings vs temperature") temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] # 查看对既定的温度取值的期望概率，即我们队所有后验样本取均值，得到一个最大后验 p(t) t = np.linspace(temperature.min() - 5, temperature.max() + 5, 50)[:, None] p_t = logistic(t.T, beta_samples, alpha_samples) print "p_t" print p_t mean_prob_t = p_t.mean(axis=0) plt.plot(t, mean_prob_t, lw=3, label="average posterior \nprobability of defect") # 求均值即求对特定温度的下事件发生的期望概率 plt.plot(t, p_t[0, :], ls="--", label="realization from posterior") # 在不同的模型参数后验估计下，事件发生概率的后验估计 plt.plot(t, p_t[-2, :], ls="--", label="realization from posterior") # 在不同的模型参数后验估计下，事件发生概率的后验估计 plt.scatter(temperature, D, color="k", s=50, alpha=0.5) plt.title("Posterior expected value of probability of defect; \ plus realizations") plt.legend(loc="lower left") plt.ylim(-0.1, 1.1) plt.xlim(t.min(), t.max()) plt.ylabel("probability") plt.xlabel("temperature") plt.show()
1. 95%置信区间（CI）
一个有趣的问题是：在哪个温度时，我们对缺陷发生的概率是最不确定的？
我们通过编码来看看，期望会的曲线和每个点对应的95%的置信区间（CI）
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm from scipy.stats.mstats import mquantiles np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] # plot it, as a function of temperature (the first column) print("Temp (F), O-Ring failure?") print(challenger_data) plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k", alpha=0.5) plt.yticks([0, 1]) plt.ylabel("Damage Incident?") plt.xlabel("Outside temperature (Fahrenheit)") plt.title("Defects of the Space Shuttle O-Rings vs temperature") temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] # 查看对既定的温度取值的期望概率，即我们队所有后验样本取均值，得到一个最大后验 p(t) t = np.linspace(temperature.min() - 5, temperature.max() + 5, 50)[:, None] p_t = logistic(t.T, beta_samples, alpha_samples) print "p_t" print p_t mean_prob_t = p_t.mean(axis=0) # vectorized bottom and top 2.5% quantiles for "confidence interval" qs = mquantiles(p_t, [0.025, 0.975], axis=0) plt.fill_between(t[:, 0], *qs, alpha=0.7, color="#7A68A6") plt.plot(t[:, 0], qs[0], label="95% CI", color="#7A68A6", alpha=0.7) plt.plot(t, mean_prob_t, lw=1, ls="--", color="k", label="average posterior \nprobability of defect") plt.xlim(t.min(), t.max()) plt.ylim(-0.02, 1.02) plt.legend(loc="lower left") plt.scatter(temperature, D, color="k", s=50, alpha=0.5) plt.xlabel("temp, $t$") plt.ylabel("probability estimate") plt.title("Posterior probability estimates given temp. $t$") plt.show()
95%置信区间，或者称95% CI，在图中用紫色区块显示。
我们知道，pymc得到模型参数 ?, ? 并不是一个固定的值，而是一个概率分布，而不同的 ?, ? 得到不同的逻辑函数，因此，对每一个温度 t，会有不同的事故发生概率后验值。
上图紫色区域的含义是：对每一个温度值，它都包含了95%的分布，例如在 65 度时，有 95%的确信度相信事故发生的概率范围为 0.25 到 0.75 之间。
更一般情况，当温度到达 60 度时，CI’s 快速变窄；当温度超过 70 度时 CI’s 也有类似的情况。这将给我提供接下来如何分析的思路。
1. 首先，在低温区间，尤其是50°以下，95%置信区间是很窄的，说明在低温时候发生事故的概率是非常大的，且这个事实非常确定； 2. 我们可能在温度区间 60-65 上做更多 的 O 形密封圈的测试，这样就能得到更好的估计概率。
2. 挑战者号事故当天发生了什么
挑战者号发生灾难的那天，外部的温度是 31 华氏度。在此时的温度下，发生灾难的后验分布如何?
还是一样注意，贝叶斯推断得到的后验分布是一个概率分布（这点在刚开始学习概率编程时候特别不适应），因此分布图如下。
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm from scipy.stats.mstats import mquantiles np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] # plot it, as a function of temperature (the first column) print("Temp (F), O-Ring failure?") print(challenger_data) plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k", alpha=0.5) plt.yticks([0, 1]) plt.ylabel("Damage Incident?") plt.xlabel("Outside temperature (Fahrenheit)") plt.title("Defects of the Space Shuttle O-Rings vs temperature") temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] prob_31 = logistic(31, beta_samples, alpha_samples) plt.xlim(0.995, 1) plt.hist(prob_31, bins=1000, normed=True, histtype='stepfilled') plt.title("Posterior distribution of probability of defect, given $t = 31$") plt.xlabel("probability of defect occurring in O-ring"); plt.show()
看起来挑战者号必定会发生由于 O 形密封环导致的灾难。
0x4：我们的模型适用吗？
接下来的问题是，如何验证我们的模型的吻合度是否不好？
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm from scipy.stats.mstats import mquantiles np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] # plot it, as a function of temperature (the first column) print("Temp (F), O-Ring failure?") print(challenger_data) plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k", alpha=0.5) plt.yticks([0, 1]) plt.ylabel("Damage Incident?") plt.xlabel("Outside temperature (Fahrenheit)") plt.title("Defects of the Space Shuttle O-Rings vs temperature") temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] simulated = pm.Bernoulli("bernoulli_sim", p) N = 10000 mcmc = pm.MCMC([simulated, alpha, beta, observed]) mcmc.sample(N) simulations = mcmc.trace("bernoulli_sim")[:] print(simulations.shape) plt.title("Simulated dataset using posterior parameters") for i in range(4): ax = plt.subplot(4, 1, i + 1) plt.scatter(temperature, simulations[1000 * i, :], color="k", s=50, alpha=0.6) plt.show()
我们希望从统计上评估上面的仿真数据和我们的原始观测数据是否真的很相似，即我们希望评估出模型模拟得究竟有多好。需要找一种量化的方式来量化衡量这里所谓的“好”。
1. 分离图
我们使用作图的方式完成这件事，下面的图形测试是一种用于逻辑回归拟合程度的数据可视化方法。这种图被称为分离图。分离图可以让用户用一种图形化的方法对比不同的模型并从中选出最适合的。
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm from scipy.stats.mstats import mquantiles np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] simulated = pm.Bernoulli("bernoulli_sim", p) N = 10000 mcmc = pm.MCMC([simulated, alpha, beta, observed]) mcmc.sample(N) simulations = mcmc.trace("bernoulli_sim")[:] print(simulations.shape) # 对每一个模型，给定温度，计算后验模拟产生值1的次数比例，并对所有返回的模拟值取均值，这样我们得到在每一个数据点上发生缺陷的后验可能 posterior_probability = simulations.mean(axis=0) print("posterior prob of defect | realized defect ") for i in range(len(D)): print("%.2f | %d" % (posterior_probability[i], D[i])) # 接下来我们根据后验概率对每一列进行排序 ix = np.argsort(posterior_probability) print("probb | defect ") for i in range(len(D)): print("%.2f | %d" % (posterior_probability[ix[i]], D[ix[i]]))
我们可以在一个图中更好地可视化展示这些数据。
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm from scipy.stats.mstats import mquantiles from separation_plot import separation_plot np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] simulated = pm.Bernoulli("bernoulli_sim", p) N = 10000 mcmc = pm.MCMC([simulated, alpha, beta, observed]) mcmc.sample(N) simulations = mcmc.trace("bernoulli_sim")[:] print(simulations.shape) # 对每一个模型，给定温度，计算后验模拟产生值1的次数比例，并对所有返回的模拟值取均值，这样我们得到在每一个数据点上发生缺陷的后验可能 posterior_probability = simulations.mean(axis=0) # 可视化展示 separation_plot(posterior_probability, D) plt.show()
蛇形线表示排序后的后验概率，蓝色柱子表示真实发生的缺陷（观测结果），空的地方表示没有发生缺陷。
随着概率增加，可以看到事故发生的概率越来越大，在图中的右边，说明后验概率越来越大（因为蛇形线接近 1），表示更多的事故发生。理想情况所有的蓝色柱子都要靠近右边，这里面会存在一些偏差，偏差存在的地方预测也就缺失了。
垂直的线表示在这个模型下，所能观测到的事故次数。这根垂直的线可以让我们比较模型预测出的事故次数和真实事故次数各是多少。
将这个模型的分离图和其它模型的分离图进行对比，会发现更多的有用信息。下面将上面的模型与其它三个进行对比：
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm from scipy.stats.mstats import mquantiles from separation_plot import separation_plot np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] simulated = pm.Bernoulli("bernoulli_sim", p) N = 10000 mcmc = pm.MCMC([simulated, alpha, beta, observed]) mcmc.sample(N) simulations = mcmc.trace("bernoulli_sim")[:] print(simulations.shape) # 对每一个模型，给定温度，计算后验模拟产生值1的次数比例，并对所有返回的模拟值取均值，这样我们得到在每一个数据点上发生缺陷的后验可能 posterior_probability = simulations.mean(axis=0) # Our temperature-dependent model separation_plot(posterior_probability, D) plt.title("Temperature-dependent model") plt.show() # Perfect model # i.e. the probability of defect is equal to if a defect occurred or not. p = D separation_plot(p, D) plt.title("Perfect model") plt.show() # random predictions p = np.random.rand(23) separation_plot(p, D) plt.title("Random model") plt.show() # constant model constant_prob = 7. / 23 * np.ones(23) separation_plot(constant_prob, D) plt.title("Constant-prediction model") plt.show()
1. 如果模型很完美，当事故确实发生时，它预测的后验概率是 1 2. 一个随机模型预测的结果完全是随机的，将和温度没有任何关系 3. 对于一个常数模型，即P(? = 1|?) = ?, ∀?。?的最优值是观测到的事故发生的概率，即? = 7⁄23
在随机模型中，可以看到，随着概率的增加，事故并没有向右侧聚集。常数模型也是这个情况。
完美模型中没有显示概率曲线。因为它一直位于柱子的顶端。此处的完美模型只是用于演示，并不能从中推断出什么科学依据。
在实际项目中，我们要评估的是我们的后验概率估计的分离图是如何介于随机模型和完美模拟之间的，自然，越靠近完美模型意味着我们的模型拟合优度越好。
扫一扫
6120
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
