（本文假设读者懂贝叶斯统计分析的基本框架，主要是为了记录自己的学习，因此可能很多地方比较简略，并不适合当作入门读物。代码部分可以跳过。）

1986年，美国挑战者号航天飞机升空失败，导致数名宇航员死亡。事后调查发现，事故原因是 O 型环的设计问题，导致其对一些外在条件非常敏感，其中一个影响条件是外部温度。在此前的 24 次类似发射中，有 23 次的发射数据保留了下来，其中 7 次发生了O型环出故障。在挑战者号发射前夜，虽然工作人员讨论了 O 型环问题，但是认为这过往的 7 次故障没有体现出明显的规律，因此没有留意。

下图是过往 23 次发射记录中的O型环故障情况图，横轴是外部温度（华氏度），纵轴表示是否发生故障。
从图中来看，一个大致的趋势是温度越低，越可能发生故障（65度左右是一个可能的分水岭）。我们是否可以从中得到一个关于温度和故障几率的函数 p(t)，给定一个温度 t，返回在该温度下O型环会故障的概率（或概率分布）？

让我们试着用贝叶斯统计来算出这个函数。

确定模型

很多函数可以做到给定一个值，返回一个 [0, 1] 之间的值，其中最常见的就是 logistic 函数：

其中 β 和 α 是需要通过贝叶斯统计计算出的参数。如果没有 α，则 t=0 会是 p(t) 的转折点，为了让分水岭发生在 t=65 度左右，所以需要一个 α 参数来平移函数图像

计算模型参数

α 和 β 可以取任何值，但是肯定有一个相对比较合理的区间，只是我们不知道这个区间，因此可以使用正态分布来模拟这两个参数。

import pymc3 as pm
import theano as ttwith pm.Model() as model:beta = pm.Normal('beta', mu=0, tau=0.001, testval=0)alpha = pm.Normal('alpha', mu=0, tau=0.001, testval=0)p = pm.Deterministic("p", 1.0/(1. + tt.exp(beta*temperature + alpha)))

参数 tau 就是表示正态分布中的方差，准确来说 tau = 1/variance²，因此 alpha 和 beta 都符合平均数 0，方差为 101010\sqrt{10}1010 的正态分布。另外 alpha 和 beta 的初始值（testval 参数决定）都是 0，这仅仅是为了避免如果让这两个参数的初始值随机决定，可能会非常大，这两个值非常大就会导致 p 非常近似 0 或 1，而 pymc3 的伯努利函数在处理概率是 0 或 1 时的情况下表现很不好，因此设置 testval =0 纯粹只是一个为了照顾 pymc3 计算能力的一个设定

有了 p（即给定温度后的故障概率），则任意一个O型环的故障概率就是一个伯努利事件 Ber(p(t))

这个伯努利分布如果展开的话，公式比较复杂，因此我们使用 MCMC 方法进行抽样，通过抽样结果来还原这个分布（其中 observered=D 中的 D 是已有的那 23 条数据）

with model:observed = pm.Bernoulli("bernoulli_obs", p, observed=D)start = pm.find_MAP()step = pm.Metropolis()trace = pm.sample(120000, step=step, start=start)burned_trace = trace[100000::2]

posterior distribution of α and β

使用 MCMC 方法后，现在我们已经有了 α，β 和 p 的 posterior distribution。

下图中的 α 和 β 的分布都是从 burned_trace 中数据制作的图。可以看到，β 始终大于 0，且小于0.5。alpha则小于 -5，大于 -35

β 始终大于 0 正说明温度和故障率有关系，否则 β 就会 =0，让温度这个变量从公式中消失。当然由于 prior distribution 信息不足，数据量又较小，因此 α 的 β 的分布范围都很宽

现在我们有了 α 的 samples，也有了 β 的samples（都来自 burned_trace），则每一对 α 和 β 都对应一个 logistic 函数，每一个函数给定一个温度 t 都会产生一个故障概率，因此给定一个温度 t 会有很多个故障概率，将这些概率取平均，就是这个温度 t 下预期的故障概率。

下图紫色范围就是给定一个温度，故障率有 95% 可能落在紫色区间内。可以看到，大约在 60 度左右时，紫色范围最宽，大约在 0.25-0.75 之间，过了 70 度，紫色范围收窄。这意味着我们应该在 60-65 度这个区间内对O型环做更多测试

在挑战者号发生事故的当天，温度是 31 度，根据上面的 logistic 函数计算，我们可以画出在 31 度下发生故障的可能性。只需要每一对 α 和 β 产生的 logistic 函数都带入 31 度，就能得到故障概率的频数

下图中可以看到，在数千对 α 和 β 对应的数千个 logistic 函数中，有差不多 3000 个函数都显示故障率达到了 100%，概率非常高

因此，如果单纯从这个分析来看，O 型环的故障并非毫无迹象可言。

但是我个人相信，航天飞机的发射并不会只考虑这一个因素，而是很多因素的综合考虑。以上分析不但有事后诸葛亮之嫌，而且即便事先得出了同样的结论，也未必就意味着绝对不应该发射。此案例主要供学习使用而已。

以上所有图片和代码来自《Bayesian Methods for Hackers》一书。仅供参考

贝叶斯小数据分析—— 23 条数据决定宇航员生死（使用 PyMC3）相关推荐

机器学习实战：朴素贝叶斯算法在新闻文本数据上的分类表现
https://www.toutiao.com/a6647102437532369421/ 2019-01-17 08:01:00 大家好,今天跟大家学习一下通过sklearn的朴素贝叶斯模型实战.前 ...
机器学习（8）朴素贝叶斯算法（20条新闻分类）
目录一.基础理论二.实战:20条新闻分类 1.读取数据 2.训练集划分 3.特征工程(文本特征提取) 4.朴素贝叶斯算法训练 5.模型评估方法一:预测值与真实值比对方法二:计算准确率总代码 ...
【视频】线性回归中的贝叶斯推断与R语言预测工人工资数据|数据分享
最近我们被客户要求撰写关于线性回归的研究报告,包括一些图形和统计输出. 在这个视频中,我们转向简单线性回归中的贝叶斯推断. 我们将使用一个参照先验分布,它提供了频率主义解决方案和贝叶斯答案之间的联系. ...
清华团队通过监督贝叶斯嵌入，对单细胞染色质可及性数据进行细胞类型注释...
本文约3200字,建议阅读9分钟本文介绍了清华团队在单细胞技术的最新进展. 单细胞技术的最新进展使得能够在细胞水平上表征表观基因组异质性.鉴于细胞数量呈指数增长,迫切需要用于自动细胞类型注释的计算方 ...
自编程实现朴素贝叶斯算法，Navie Bayes程序（python），并对鸢尾花数据进行分类。
自编程实现朴素贝叶斯算法,Navie Bayes程序(python),并对鸢尾花数据进行分类. 目录自编程实现朴素贝叶斯算法,Navie Bayes程序(python),并对鸢尾花数据进行分类. 朴 ...
贝叶斯分位数回归、lasso和自适应lasso贝叶斯分位数回归分析免疫球蛋白、前列腺癌数据...
原文链接:http://tecdat.cn/?p=22702 贝叶斯回归分位数在最近的文献中受到广泛关注,本文实现了贝叶斯系数估计和回归分位数(RQ)中的变量选择,带有lasso和自适应lasso惩罚 ...
[zt]数学之美番外篇：平凡而又神奇的贝叶斯方法
数学之美番外篇:平凡而又神奇的贝叶斯方法 Tags: 数学, 机器学习与人工智能, 计算机科学 save it69 saved tags: 贝叶斯 math bayesian algorithm 数学 ...
04机器学习实战之朴素贝叶斯
朴素贝叶斯概述贝叶斯分类是一类分类算法的总称,这类算法均以贝叶斯定理为基础,故统称为贝叶斯分类.本章首先介绍贝叶斯分类算法的基础--贝叶斯定理.最后,我们通过实例来讨论贝叶斯分类的中最简单的一种: ...
NB贝叶斯平凡而又神奇的贝叶斯方法
转自:http://mindhacks.cn/2008/09/21/the-magical-bayesian-method/ 概率论只不过是把常识用数学公式表达了出来. --拉普拉斯目录 0. 前言 ...
非常全面的贝叶斯网络介绍非常多的例子说明
0. 前言这是一篇关于贝叶斯方法的科普文,我会尽量少用公式,多用平白的语言叙述,多举实际例子.更严格的公式和计算我会在相应的地方注明参考资料.贝叶斯方法被证明是非常 general 且强大的推理框架 ...

贝叶斯小数据分析—— 23 条数据决定宇航员生死（使用 PyMC3）

确定模型

计算模型参数

posterior distribution of α and β

贝叶斯小数据分析—— 23 条数据决定宇航员生死（使用 PyMC3）相关推荐

最新文章

热门文章