快速上手 — TensorFlow Probability 内置概率编程教材

文 / Mike Shwe，Google TensorFlow Probability 产品经理、Josh Dillon，Google TensorFlow Probability 软件工程师、Bryan Seybold、Matthew McAteer 和 Cam Davidson-Pilon Google 软件工程师

刚开始接触概率编程？对 TensorFlow Probability (TFP) 还不熟悉？快来看看我们为您准备了什么。《黑客的贝叶斯方法》(Bayesian Methods for Hackers) 是一本入门级实践教程，现内置于 TFP 中，而且提供有许多示例。作为所有人可用的开源资源，该 TFP 版本对此前 PyMC3 中编写的版本做了补充。

教程中提供有大量关于《黑客的贝叶斯方法》的精华知识，而且演示了如何运用概率编程解决现实问题，即便是概率编程初学者也可以轻松上手。

人人可学的概率编程

贝叶斯方法虽非概率编程之必需，但它提供了一个直观的框架，可以用于表征信念并根据新数据对这些信念进行升级。《黑客的贝叶斯方法》以 TFP 作为基础，采用注重实践的方式教授这些方法。这本书是在 Google Colab 中编写的，欢迎您运行和修改其中的 Python 示例。

TensorFlow 团队构建的 TFP 面向希望对领域知识进行编码以理解数据并进行预测的数据科学家、统计学家、ML 研究人员和从业者。TFP 是在 TensorFlow 的基础上构建的 Python 库，可以在现代硬件上轻松结合概率模型和深度学习。TFP 让您能够：

• 以交互方式探索您的数据

• 快速评估不同模型

• 自动利用现代矢量化硬件加速器

• 轻松自信地发布。TFP 经过专业构建和测试，支持 Google Cloud 环境，还有活跃的开源社区提供支持

我们在相关博文中讨论过，概率编程的应用非常广泛，包括金融和油气行业等。为什么会这样？ 因为不确定性无处不在。现实世界的各种现象经常受到我们无法控制甚或并未意识到的种种外部因素的影响，就连我们完全了解的现象也是如此。如果忽略这些因素，就可能得出错误或至少是有误导性的结论。我们构建了人人可用的 TFP，用于模拟遍布我们周围的不确定性。

解决现实问题

许多贝叶斯教程侧重于解决有分析解法的简单问题，比如抛硬币和掷骰子。而《黑客的贝叶斯方法》虽然也是从这些问题开始，但却迅速转向更加现实的问题。其中的示例范围广泛，从了解宇宙现象到探测在线用户的行为变化等都有涉猎。

在本篇文章的余下篇幅中，我们将概括介绍一个著名的现实问题，该问题在书中第 2 章 “1986 年挑战者号航天飞机失事” 中有详细表述。

1986 年 1 月 28 日，在美国航天飞机的第 25 次飞行中，挑战者号航天飞机的两个固态火箭助推器中有一个由于 O 型密封圈失效而发生了爆炸。虽然参与任务的工程师就此前飞行中出现的损坏问题与 O 型密封圈制造商进行过多次沟通，但制造商认为这种风险在可接受的范围内 [1]。

下图绘制的是对此前航天飞机任务中七次 O 型密封圈损坏事故的观测结果，这是损坏事故与环境温度的函数（在 70 华氏度时，发生过两起事故）。

可以看到，随着温度降低，O 型密封圈的损坏比例会显著增加，但是没有明确的温度阈值，即低于某个温度，O 型密封圈很可能会失效。与大多数现实世界的现象一样，这个问题存在不确定性。我们想要确定在给定温度 t 时，O 型密封圈失效的概率是多少？

我们可以用 一个逻辑函数，对 O 型密封圈在温度为 t 时 的损坏概率 p 建模，具体如下：

其中 β 决定逻辑函数的形状，𝛼 是偏移项，可以使函数从左向右移动。由于这些参数既可为正值也可为负值，并且没有特定的边界或大小偏差，我们可以将其作为高斯随机变量建模：

在 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 版本的第 2 章 ，运行完整的航天飞机示例）。

请注意，我们在第 8 行得出 p(t) 的实际值为 0 或 1，这是使用之前在第 6 和第 7 行中 𝛼 和 β 的采样值对逻辑函数进行采样得出的结果。此外，请注意 evaluate() 辅助函数让我们可以在图表和 eager 模式间无缝转换，同时将张量值转化为 NumPy。我们在第 2 章的开始部分对 eager 和图表模式，以及此辅助函数做了更详尽的描述。

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

我们想要根据此生成模型找到模型参数，这样模型就可以解释我们观测到的数据，这正是概率推理的目标！

TFP 使用未归一化联合对数概率函数来评估模型，从而执行概率推理。此 joint_log_prob 的参数是数据和模型状态。该函数会返回参数化模型生成所观测数据的联合概率对数。如需了解有关 joint_log_prob 的详细信息，请参阅 https://github.com/tensorflow/probability/blob/master/discussion/joint_log_prob.md

在挑战者号示例中，我们定义 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 类型，称作 Hamiltonian Monte Carlo：

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，因为所有示例的结果都是负值。

正如前文所言，我们真正想要知道的是： 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% 确信失效概率在 1.0 至 0.80 之间。但讽刺的是，许多人错误地解释为置信区间具有这一属性。

挑战者号失事当天的环境温度为 31 华氏度，这证明 O 型密封圈失效的后验分布将引导我们得出这样的结论：当天发生问题的置信度较高。

从这个很简单易懂的概率分析可以窥见 TFP 和贝叶斯方法的强大效用：它们可以用于提供有价值的数据分析，并对会产生重大后果的现实问题作出预测。

请继续往下读 ！

您会在此书中发现大量现实示例。您可以将一段时间内的文本消息量分析广泛用于制造和生产系统中的故障问题检测。在我们首次完成此章节 TFP 版本草稿的几周后，Google 的软件工程师就采用这些方法进行文本消息分析，用于理解生产软件的片状文本。

您也可以找到一种分析方法，用于寻找宇宙中的暗物质。还可以预测上市公司的股份回报，即股票收益。

希望您会钻研此书中的概念演示，并运用这些方法解决您所处领域中的问题。请在 github 中写出评论并拉取请求，帮助我们将这本书保持在顶部位置！

致谢

感谢 TensorFlow Probability 团队协助编写 TFP 版《黑客的贝叶斯方法》。

参考文献

[1] https://en.wikipedia.org/wiki/Space_Shuttle_Challenger_disaster

更多 AI 相关阅读：