TensorFlow Probability 概率编程入门级实操教程

雷锋网 AI 科技评论按：TensorFlow Probability（TFP）是一个基于 TensorFlow 的 Python 库，能够更容易地结合概率模型和深度学习。数据科学家、统计学以及机器学习研究者或者从业者，都可以用它编码领域知识 (Domain Knowledge)，从而理解数据并写出自己的应用。针对那些对 TFP 还不那么熟悉的入门者，日前，谷歌 TensorFlow Probability 的产品经理 Mike Shwe 及软件工程师 Josh Dillon、谷歌的软件工程师 Bryan Seybold 及 Matthew McAteer 、 Cam Davidson-Pilon 共同在 TensorFlow 官网上发布介绍 TensorFlow Probability 的入门级实操性教程——《Bayesian Methods for Hackers》的文章，雷锋网 AI 科技评论编译如下。

之前没有学过概率编程？对 TensorFlow Probability （TFP）还不熟悉？下面我们为你准备了入门级实操性教程——《Bayesian Methods for Hackers》（教程查看地址： https://camdavidsonpilon.github.io/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/#tensorflow ），这门教程的实例现在也在 TFP 中开放了。作为对所有人开放的开源资源，TFP 版本的概率编程对之前用 PyMC3 写的那版进行了补充。

《Bayesian Methods for Hackers》具备许多优势：它不仅能让概率初学者较容易上手，而且还展示了如何将概率编程应用于现实问题。

每个人都可以学的概率编程

虽然概率编程不要求贝叶斯方法（Bayesian approach），但是该方法提供了一个相对直观的框架，来表示信念（representing beliefs），并基于新的数据来更新这些信念。《Bayesian Methods for Hackers》使用 TFP 为基础，以实操的方式来教授这些技术。由于这本书由 Google Colab 所写，你可以运行并修改其中的 Python 示例。

TensorFlow 团队开发的 TFP 专门面向数据科学家、统计学家以及以及想要编码领域知识（Domain Knowledge）来了解数据并进行数据预测的机器学习研究者和从业者。TFP 是基于 TensorFlow 的 Python 开发库，能够更容易地结合概率模型和先进硬件上的深度学习。TFP 可以让你：

互动地探究数据
快速地评估不同的模型
自动地利用先进的、矢量化的硬件加速器
更轻易、更有把握地启用。TFP 经专业地开发和测试，可使用现成的 Google-Cloud，并且还有一个强大的开源社区的支持。

正如我们在相关博文中曾讨论过的，概率编程有非常多样化的应用，包括金融、石油和天然气等各行各业。为什么？——因为不确定性无处不在。现实世界的现象——即便是那些我们完全了解的现象，都受到我们无法控制甚至无法意识到的外部因素的影响。如果忽略这些因素，这些模型所得出的结论可能往好里说是误导性的或者是完全错误的。我们开发出了对所有场景都可用的 TFP，就是为了对我们周围所有的不确定性进行建模。

解决现实世界的问题

许多贝叶斯教程都是聚焦于解决那些已有分析结果的简单问题：比如掷硬币和掷骰子的问题。不过《Bayesian Methods for Hackers》一书是从这些简单的问题开始，之后迅速转向更加现实的问题，例如从理解宇宙到检测某个在线用户的行为变化等。

在这篇文章接下来的部分，我们将一起探讨一个非常有名的现实世界的问题，这个问题在 Bayesian Hackers 一书中的第二章节（查看地址： https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Chapter2_MorePyMC/Ch2_MorePyMC_TFP.ipynb ）有更详细的描述：1986 年航天飞机挑战者的灾难。

1986 年 1 月 28 日，美国挑战者号航天飞机的第 25 次飞行中，由于一处 O 形圈故障，挑战者号的两个固体火箭助推器其中的一个发生了爆炸。虽然控制中心的工程师们与 O 形圈制造商就先前飞行中的损坏进行了多次沟通，但制造商坚持认为风险是可以接受的。

下图将对先前航天飞机任务中的七次 O 形圈损坏事件的观测，描述成了一个环境温度的函数。（在 70 度时，有两次损坏事件。）

F7zeMju.png!web

你会注意到，随着温度的降低，O 形圈损坏的比例会明显增多，不过没有明显的温度阈值——低于该阈值时 O 形圈就一定会失效。与现实世界中的大多数现象一样，这个问题涉及不确定性。我们希望在给定温度 t 下，来确定 O 形圈失效的概率。

我们可以特别使用逻辑函数模拟温度 t 下 O 形环损坏的概率 p ：

VbYZBnZ.png!web

其中β确定概率函数的形状，α是偏移项，控制函数从左向右移动。由于这两个参数都既是正的或负的，也没有特定的边界或大小的偏差，我们可以将它们建模为高斯分布随机变量：

YjYruma.png!web

在 TFP 中，我们可以用 tfp.distributions.Normal 直观地表示 α 和 β，其代码片段如下：

temperature_ = challenger_data_[:, 0]    
temperature = tf.convert_to_tensor(temperature_, dtype=tf.float32)    
D_ = challenger_data_[:, 1]                # defect or not?    
D = tf.convert_to_tensor(D_, dtype=tf.float32)    
 
beta = tfd.Normal(name="beta", loc=0.3, scale=1000.).sample()    
alpha = tfd.Normal(name="alpha", loc=-15., scale=1000.).sample()    
p_deterministic = tfd.Deterministic(name="p", loc=1.0/(1. + tf.exp(beta * temperature_ + alpha))).sample()    
 
[    
prior_alpha_,    
prior_beta_,    
p_deterministic_,    
D_,    
] = evaluate([    
alpha,    
beta,    
p_deterministic,    
D,    
])

(如果要运行这个代码片段，请前往本书第二章的 Google Colab 版本，来运行整个航天飞机示例)。

要注意的是，我们在第 8 行得到 p(t) 的实际值 0 或 1，其中我们使用此前在第 6 行和第 7 行中采样的α和β值从概率函数（logistic function）中采样。此外，注意 evaluate() 辅助函数可以让我们实现图表和 eager 模式之间的无缝转换，与此同时将张量值转换为 numpy。关于 eager 和图表模型以及这一辅助函数的更详细的讲解，请查看第二章节的开头部分。

为了将温度 t 下的失效概率 p(t) 与我们观测的数据联系起来，我们可以使用带参数 p(t) 的伯努利随机变量。要注意的是，一般而言 Ber(p) 是随机变量，它的值为 1 的概率为 p，其它情况下都为 0。因此，生成模型的最后一部分是温度值为 ? 的情况下观测到的有缺陷事件的数量 D?，我们可以对其建模为：

ENrABfV.png!web

给定这一生成模型的情况下，我们希望找到模型参数从而让模型能够解释所观察到的数据——这正是是概率推理的目标。

TFP 通过使用非标准化的联合对数概率函数评估模型来执行概率推断。此 joint_log_prob 的论点是数据和模型状态。该函数返回参数化模型生成观测数据的联合概率的对数。如果要了解更多关于 joint_log_prob 的信息，可以查看这个短文介绍。

针对上述挑战者的示例，这里我们定义 joint_log_prob 的方式如下：

def challenger_joint_log_prob(D, temperature_, alpha, beta):    
"""    
   Joint log probability optimization function.    
 
   Args:    
     D: The Data from the challenger disaster representing presence or     
        absence of defect    
     temperature_: The Data from the challenger disaster, specifically the temperature on     
        the days of the observation of the presence or absence of a defect    
     alpha: one of the inputs of the HMC    
     beta: one of the inputs of the HMC    
   Returns:     
     Joint log probability optimization function.    
   """    
rv_alpha = tfd.Normal(loc=0., scale=1000.)    
rv_beta = tfd.Normal(loc=0., scale=1000.)    
logistic_p = 1.0/(1. + tf.exp(beta * tf.to_float(temperature_) + alpha))    
rv_observed = tfd.Bernoulli(probs=logistic_p)    
 
return (    
rv_alpha.log_prob(alpha)    
+ rv_beta.log_prob(beta)    
+ tf.reduce_sum(rv_observed.log_prob(D))    
)

注意 15-18 行怎样简单地对生成模型、每行的随机变量进行编码。同时也要注意 rv_alpha 和 rv_beta 表示此前对 ? 和 β 分布的随机变量。相比之下， rv_observed 表示给定一个参数为 ? 和 β 的概率分布情况下温度和 O 形圈输出的观察可能性的条件分布。

接下来，我们使用 joint_log_prob 函数，并将其发送到 tfp.mcmc 模块。马尔可夫链蒙特卡洛（MCMC）算法对未知的输入值进行有依据的猜测，并计算 joint_log_prob 函数中参数集的可能性。通过多次重复这一过程，MCMC 构建了可能的参数的分布，而构建这一分布是概率推理的目标。

因此，我们将在 challenge_joint_log_prob 函数上设置一种特定类型的称为「哈密顿蒙特卡洛」的 MCMC：

number_of_steps = 60000    
burnin = 50000    
 
# Set the chain's start state.    
initial_chain_state = [    
0. * tf.ones([], dtype=tf.float32, name="init_alpha"),    
0. * tf.ones([], dtype=tf.float32, name="init_beta")    
]    
 
# Since HMC operates over unconstrained space, we need to transform the    
# samples so they live in real-space.    
unconstraining_bijectors = [    
tfp.bijectors.Identity(),    
tfp.bijectors.Identity()    
]    
 
# Define a closure over our joint_log_prob.    
unnormalized_posterior_log_prob = lambda *args: challenger_joint_log_prob(D, temperature_, *args)    
 
# Initialize the step_size. (It will be automatically adapted.)    
with tf.variable_scope(tf.get_variable_scope(), reuse=tf.AUTO_REUSE):    
step_size = tf.get_variable(    
name='step_size',    
initializer=tf.constant(0.5, dtype=tf.float32),    
trainable=False,    
use_resource=True    
)    
 
# Defining the HMC    
hmc=tfp.mcmc.TransformedTransitionKernel(    
inner_kernel=tfp.mcmc.HamiltonianMonteCarlo(    
target_log_prob_fn=unnormalized_posterior_log_prob,    
num_leapfrog_steps=2,    
step_size=step_size,    
step_size_update_fn=tfp.mcmc.make_simple_step_size_update_policy(),    
state_gradients_are_stopped=True),    
bijector=unconstraining_bijectors)    
 
# Sampling from the chain.    
[    
posterior_alpha,    
posterior_beta    
], kernel_results = tfp.mcmc.sample_chain(    
num_results = number_of_steps,    
num_burnin_steps = burnin,    
current_state=initial_chain_state,    
kernel=hmc)    
 
# Initialize any created variables for preconditions    
init_g = tf.global_variables_initializer()

最终，我们将真正通过 evaluate() 辅助函数进行推理：

evaluate(init_g)    
[    
posterior_alpha_,    
posterior_beta_,    
kernel_results_    
] = evaluate([    
posterior_alpha,    
posterior_beta,    
kernel_results    
])    
 
print("acceptance rate: {}".format(    
kernel_results_.inner_results.is_accepted.mean()))    
print("final step size: {}".format(    
kernel_results_.inner_results.extra.step_size_assign[-100:].mean()))    
 
alpha_samples_ = posterior_alpha_[burnin::8]    
beta_samples_ = posterior_beta_[burnin::8]

绘制 α 和 β 的分布图时，我们注意到这两个参数的分布相当宽泛，正如我们所预期的那样，数据点相当少，并且失效和非失效观测的温度出现重叠。然而，即使分布宽泛，我们也可以相当确定温度确实对 O 形圈损坏的概率有影响，因为β的所有样本都大于 0。同样，我们也可以确信 α 明显小于 0，因为所有的样本都能很好地转化为负数。

jA7riiz.png!web

正如我们上面所提到的，我们真正想知道的是：在给定温度下，O 形环损坏的预期概率是多少？为了计算这一概率，我们可以对来自后验的所有样本求平均值，从而得到概率 p(t?) 的可能值。

alpha_samples_1d_ = alpha_samples_[:, None]  # best to make them 1d    
beta_samples_1d_ = beta_samples_[:, None]    
 
beta_mean = tf.reduce_mean(beta_samples_1d_.T[0])    
alpha_mean = tf.reduce_mean(alpha_samples_1d_.T[0])    
[ beta_mean_, alpha_mean_ ] = evaluate([ beta_mean, alpha_mean ])    
 
print("beta mean:", beta_mean_)    
print("alpha mean:", alpha_mean_)    
def logistic(x, beta, alpha=0):    
"""    
   Logistic function with alpha and beta.    
 
   Args:    
     x: independent variable    
     beta: beta term     
     alpha: alpha term    
   Returns:     
     Logistic function    
   """    
return 1.0 / (1.0 + tf.exp((beta * x) + alpha))    
 
t_ = np.linspace(temperature_.min() - 5, temperature_.max() + 5, 2500)[:, None]    
p_t = logistic(t_.T, beta_samples_1d_, alpha_samples_1d_)    
mean_prob_t = logistic(t_.T, beta_mean_, alpha_mean_)    
[    
p_t_, mean_prob_t_    
] = evaluate([    
p_t, mean_prob_t    
])

我们可以在整个温度范围内计算 95％的可信区间。要注意的是，这个区间是可靠的，而不是在统计分析的频率论方法中常用到的置信区间。95％可信区间告诉我们，能够以 95％的概率确定真实值将位于此区间内。例如，在正如下图中的紫色区域所显示的，在 50 度时，我们可以 95％确定 O 形圈损坏的概率介于 1.0 和 0.80 之间。讽刺地是，许多人都错误地将置信区间解释为具有这一特征。

amYBFzI.png!web

挑战者号灾难发生的当天，其温度为 31 华氏度，这一事实证明了，O 形圈失效的后验分布将使我们高度确信挑战者号会出现损坏的问题。

VB3Q7zZ.png!web

这种非常简单的概率分析展示了 TFP 和贝叶斯方法的强大：它们可以提供有价值的分析，同时可以对可能带来严重后果的现实世界的问题进行预测。

在《Bayesian Methods for Hackers》一书中，你可以看到大量现实世界示例。对短信量随时间变化的分析，可以在制造业和生产系统中的各种故障检测问题得到广泛应用。在我们首次起草本章节的 TFP 版本的几周时间内，谷歌的软件工程师就应用了短信分析的方法来理解生产软件的文本片状（text flakiness）。

你可以找到某个分析来寻找宇宙中的暗物质。此外，书中还有预测上市公司的股票收益率的方法。

我们希望大家能深入了解本书中的概念性实操演示，并将这些技术应用到各自所在领域的问题中。欢迎大家在 Github 中发表评论和提出要求，来帮助我们对这本书不断改进，从而使其实现最佳状态！

本文参考： https://en.wikipedia.org/wiki/Space_Shuttle_Challenger_disaster

via： https://medium.com/tensorflow/an-introduction-to-probabilistic-programming-now-available-in-tensorflow-probability-6dcc003ca29e 雷锋网 (公众号：雷锋网) AI 科技评论编译。

雷锋网原创文章，未经授权禁止转载。详情见转载须知。