在本章中，我们将讨论随机性和概率。我们将通过从一组数据中选择元素，简要探讨概率的基本原理。然后，我们将学习如何使用 Python 和 NumPy 生成（伪）随机数，以及如何根据特定的概率分布生成样本。在本章的结尾，我们将介绍一些涉及随机过程和贝叶斯技术的高级主题，并使用马尔可夫链蒙特卡罗方法在一个简单模型上估计参数。
概率是特定事件发生可能性的量化。我们总是凭直觉使用概率，尽管有时形式理论可能与直觉相悖。概率论旨在描述随机变量的行为，其值未知，但该随机变量的值取某些（范围）值的概率已知。这些概率通常是几种概率分布中的一种。可以说，最著名的概率分布是正态分布，例如，正态分布可以描述某个特征在大量人群中的传播。
我们将在第 6 章处理数据和统计中再次看到概率，在这里我们讨论统计。在这里，我们将使用概率论来量化误差，并建立一个分析数据的系统理论。
在本章中，我们将介绍以下配方：
- 随机选择项目
- 生成随机数据
- 更改随机数生成器
- 生成正态分布随机数
- 处理随机过程
- 用贝叶斯技术分析转换率
- 用蒙特卡罗模拟估计参数
在本章中，我们需要标准的 scientific Python 包 NumPy、Matplotlib 和 SciPy。我们还需要 PyMC3 包作为最终配方。您可以使用最喜爱的软件包管理器安装此软件包，例如pip：
python3.8 -m pip install pymc3 此命令将安装 PyMC3 的最新版本，在编写本文时为 3.9.2。该软件包为概率规划提供了便利，其中包括执行由随机生成的数据驱动的许多计算，以了解问题解决方案的可能分布。
本章的代码可以在 GitHub 存储库的的Chapter 04文件夹中找到 https://github.com/PacktPublishing/Applying-Math-with-Python/tree/master/Chapter%2004 。
查看以下视频以查看代码的运行：https://bit.ly/2OP3FAo 。
概率和随机性的核心是从某种集合中选择一个项目。正如我们所知，从集合中选择项目的概率量化了该项目被选择的可能性。随机性描述的是根据概率从集合中选择项目，没有任何额外的偏差。随机选择的反面可以描述为确定性选择。一般来说，使用计算机复制一个纯粹的随机过程是非常困难的，因为计算机及其处理本身就是确定性的。然而，我们可以生成伪随机数序列，当正确构造时，可以证明随机性的合理近似。
在本食谱中，我们将从一个集合中选择项目，并学习与概率和随机性相关的一些关键术语，这些术语将贯穿本章。
Python 标准库包含一个用于生成（伪）随机数的模块，称为random，但在本配方中，以及在本章中，我们将使用 NumPyrandom模块。NumPyrandom模块中的例程可用于生成随机数数组，并且比标准库中的例程稍微灵活一些。像往常一样，我们以化名np进口 NumPy。
在继续之前，我们需要修正一些术语。样本空间是集合（没有重复元素的集合），事件是样本空间的子集。事件A发生的概率表示为P（A，是一个介于 0 和 1 之间的数字。概率为 0 表示事件永远不会发生，而概率为 1 表示事件肯定会发生。整个样本空间的概率必须为 1。
当样本空间是离散的，那么概率就是与每个元素相关的 0 到 1 之间的数字，其中所有这些数字的总和是 1。这意味着从集合中选择单个项（由单个元素组成的事件）的概率。我们将考虑从一个离散集合中选择项目的方法，并处理在 Ty2 T2 中生成正态分布随机数的连续的 Ty1 T1 情形。
执行以下步骤从容器中随机选择项目：
- 第一步是设置随机数生成器。目前，我们将为 NumPy 使用默认的随机数生成器，这在大多数情况下是推荐的。我们可以通过从 NumPy
random模块调用default_rng例程来实现这一点，该模块将返回一个随机数生成器的实例。我们通常在没有种子的情况下调用此函数，但对于此配方，我们将添加种子12345，以便我们的结果可以重复：
rng = np.random.default_rng(12345)
# changing seed for repeatability- 接下来，我们需要创建我们将从中选择的数据和概率。如果已经存储了数据，或者希望选择概率相等的元素，则可以跳过此步骤：
data = np.arange(15)
probabilities = np.array(
[0.3, 0.2, 0.1, 0.05, 0.05, 0.05, 0.05, 0.025,
0.025, 0.025, 0.025, 0.025, 0.025, 0.025, 0.025]
)作为一个快速的健全性测试，我们可以使用断言来检查这些概率是否确实总和为 1：
assert round(sum(probabilities), 10) == 1.0, \
"Probabilities must sum to 1"- 现在，我们可以使用随机数生成器
rng上的choice方法，根据刚刚创建的概率从data中选择样本。对于此选择，我们希望打开替换，因此多次调用该方法可以从整个data中选择：
selected = rng.choice(data, p=probabilities, replace=True)
# 0- 要从
data中选择多个项目，我们还可以提供size参数，该参数指定要选择的数组的形状。这与许多其他 NumPy 数组创建例程的shape关键字参数的作用相同。size的参数可以是整数或整数元组：
selected_array = rng.choice(data, p=probabilities, replace=True, size=(5, 5))
#array([[ 1, 6, 4, 1, 1],
# [ 2, 0, 4, 12, 0],
# [12, 4, 0, 1, 10],
# [ 4, 1, 5, 0, 0],
# [ 0, 1, 1, 0, 7]])default_rng例程创建一个新的伪随机数生成器（PRNG）实例（带或不带种子），该实例可用于生成随机数，或者如我们在配方中看到的，从预定义数据中随机选择项。NumPy 还有一个基于隐式状态的接口，用于使用直接来自random模块的例程生成随机数。但是，通常建议显式创建生成器，使用default_rng或自己创建Generator实例。以这种方式更明确是更有益的，并且应该导致更可重复的结果（在某种意义上）。
种子是传递给随机数生成器以生成值的值。生成器仅基于种子以完全确定的方式生成数字序列。这意味着具有相同种子的相同 PRNG 的两个实例将生成相同的随机数序列。如果没有提供种子，生成器通常会生成一个取决于用户系统的种子。
NumPy 的Generator类是低级伪随机位生成器的包装器，它是实际生成随机数的地方。在最新版本的 NumPy 中，默认的 PRNG 算法是 128 位置换同余生成器。相比之下，Python 的内置random模块使用 Mersenne Twister PRNG。有关 PRNG 算法不同选项的更多信息，请参见更改随机数生成器配方。
Generator实例上的choice方法根据基础BitGenerator生成的随机数进行选择。可选的p关键字参数指定与所提供数据中的每个项目相关联的概率。如果未提供此参数，则假定均匀概率，其中每个项目具有相同的被选择概率。replace关键字参数指定选择是否应进行替换。我们启用了“替换”，以便可以多次选择相同的元素。choice方法使用生成器给出的随机数进行选择，这意味着使用choice方法时，使用相同种子的两个相同类型的 PRNG 将选择相同的项目。
choice方法也可以通过传递replace=False作为参数来创建给定大小的随机样本。这保证了从数据中选择不同的项，这有利于生成随机样本。例如，这可以用于从整个用户组中选择要测试界面新版本的用户；大多数样本统计技术依赖于随机选择的样本。
许多任务都涉及生成大量的随机数，这些随机数最基本的形式是介于 0 范围内的整数或浮点数（双精度）≤ x<1。理想情况下，这些数字应该统一选择，因此如果我们绘制大量此类数字，它们应该大致均匀地分布在 0 范围内≤ x<1。
在这个配方中，我们将看到如何使用 NumPy 生成大量随机整数和浮点数，并使用直方图显示这些数字的分布。
在开始之前，我们需要从 NumPyrandom模块导入default_rng例程，并创建一个默认随机数生成器的实例，以便在配方中使用：
from numpy.random import default_rng
rng = default_rng(12345) # changing seed for reproducibility我们在随机选择项目配方中讨论了这个过程。
我们还以别名plt导入 Matplotlibpyplot模块。
执行以下步骤以生成均匀随机数据并绘制直方图以了解其分布：
- 为了生成 0 到 1 之间的随机浮点数，包括 0 而不是 1，我们在
rng对象上使用random方法：
random_floats = rng.random(size=(5, 5))
# array([[0.22733602, 0.31675834, 0.79736546, 0.67625467, 0.39110955],
# [0.33281393, 0.59830875, 0.18673419, 0.67275604, 0.94180287],
# [0.24824571, 0.94888115, 0.66723745, 0.09589794, 0.44183967],
# [0.88647992, 0.6974535 , 0.32647286, 0.73392816, 0.22013496],
# [0.08159457, 0.1598956 , 0.34010018, 0.46519315, 0.26642103]])- 为了生成随机整数，我们在
rng对象上使用integers方法。这将返回指定范围内的整数：
random_ints = rng.integers(1, 20, endpoint=True, size=10)
# array([12, 17, 10, 4, 1, 3, 2, 2, 3, 12])- 为了检查随机浮点数的分布，我们首先需要生成一个大的随机数数组，就像我们在步骤 1中所做的那样。虽然严格来说这不是必需的，但更大的样本将能够更清楚地显示分布。我们生成这些数字如下：
dist = rng.random(size=1000)- 为了显示我们生成的数字的分布，我们绘制了数据的直方图：
fig, ax = plt.subplots()
ax.hist(dist)
ax.set_title("Histogram of random numbers")
ax.set_xlabel("Value")
ax.set_ylabel("Density")结果图如图 4.1所示。如我们所见，数据大致均匀分布在整个范围内：
Figure 4.1: Histogram of randomly generated random numbers between 0 and 1
Generator接口提供了三种生成基本随机数的简单方法，不包括我们在随机选择项目配方中讨论的choice方法。除了生成随机浮点数的random方法和生成随机整数的integers方法之外，还有生成原始随机字节的bytes方法。这些方法中的每一个都在基础BitGenerator实例上调用一个相关的方法。这些方法中的每一种都允许更改生成的数字的数据类型，例如，从双精度浮点数更改为单精度浮点数。
Generator类上的integers方法通过添加endpoint可选参数，结合了旧RandomState接口上的randint和random_integers方法的功能。（在旧界面中，randint方法不包括上端点，random_integers方法包括上端点。）Generator上的所有随机数据生成方法都允许自定义生成数据的数据类型，这在旧界面中是不可能的。（此接口是在 NumPy 1.17 中引入的。）
在图 4.1中，我们可以看到我们生成的数据直方图在 0 范围内大致一致≤ x<1。也就是说，所有的钢筋都大致水平。（由于数据的随机性，它们不是完全水平的。）这是我们对均匀分布的随机数的期望，例如通过random方法生成的随机数。我们将在生成正态分布随机数配方中更详细地解释随机数的分布。
NumPy 中的random模块提供了默认 PRNG 的几种替代方案，默认 PRNG 使用 128 位置换同余生成器。虽然这是一个很好的通用随机数生成器，但它可能不足以满足您的特定需求。例如，此算法与 Python 内部随机数生成器中使用的算法非常不同。我们将遵循 NumPy 文档中关于运行可重复但适当随机的模拟的最佳实践指南。
在本食谱中，我们将向您展示如何更改为另一种伪随机数生成器，以及如何在程序中有效地使用种子。
像往常一样，我们以别名np进口 NumPy。由于我们将使用random包中的多个项目，我们也将使用以下代码从 NumPy 导入该模块：
from numpy import random您需要选择 NumPy 提供的可选随机数生成器之一（或定义您自己的随机数生成器；请参阅本配方中的*有更多…*部分）。对于这个配方，我们将使用 MT19937 随机数生成器，它使用基于 Mersenne Twister 的算法，就像 Python 内部随机数生成器中使用的算法一样。
以下步骤说明如何以可复制的方式生成种子和不同的随机数生成器：
- 我们将生成一个
SeedSequence对象，该对象可以从给定的熵源重复生成新的种子。我们可以将自己的熵作为一个整数提供，就像我们为default_rng提供种子一样，或者我们可以让 Python 从操作系统收集熵。我们将在这里使用后一种情况来演示它的用法。为此，我们不提供任何其他参数来创建SeedSequence对象：
seed_seq = random.SeedSequence()- 现在，我们有了一种方法来为会话的其余部分生成随机数生成器的种子，接下来我们将记录熵，以便稍后在必要时重现此会话。下面是熵应该是什么样子的一个例子；您的结果必然会有所不同：
print(seed_seq.entropy)
# 9219863422733683567749127389169034574- 现在，我们可以创建底层的
BitGenerator实例，为包装Generator对象提供随机数：
bit_gen = random.MT19937(seed_seq)- 接下来，我们围绕这个
BitGenerator实例创建包装Generator对象，以创建一个可用的随机数生成器：
rng = random.Generator(bit_gen)如随机选择项目配方中所述，Generator类是底层BitGenerator的包装器，实现给定的伪随机数算法。NumPy 通过BitGenerator类PCG64（默认）的不同子类提供了几种伪随机数算法的实现；MT19937（如本配方所示）；Philox；和SFC64。这些位生成器是在 Cython 中实现的。
PCG64生成器应提供具有良好统计质量的高性能随机数生成。（32 位系统上可能不是这种情况。）MT19937生成器比更现代的 PRNG 慢，不能生成具有良好统计特性的随机数。但是，这是 Python 标准库random模块使用的随机数生成器算法。Philox生成器相对较慢，但生成的随机数质量非常高，SFC64生成器速度快，质量好，但缺少其他生成器中可用的一些功能。
本配方中创建的SeedSequence对象是一种以独立且可复制的方式为随机数生成器创建种子的方法。特别是，如果您需要为多个并行进程创建独立的随机数生成器，但仍然需要能够在以后重建每个会话以调试或检查结果，那么这将非常有用。存储在该对象上的熵是从操作系统收集的 128 位整数，用作随机种子的来源。
SeedSequence对象允许我们为彼此独立的每个进程/线程创建一个单独的随机数生成器，从而消除可能导致结果不可预测的任何数据竞争问题。它还生成彼此非常不同的种子值，这有助于避免某些 PRNG 的问题（例如 MT19937，它可以生成具有两个相似的 32 位整数种子值的非常相似的流）。显然，当我们依赖于这些值的独立性时，使用两个独立的随机数生成器生成相同或非常相似的值是有问题的。
BitGenerator类用作原始随机整数生成器的公共接口。前面提到的类是在 NumPy 中通过BitGenerator接口实现的类。您还可以创建自己的BitGenerator子类，尽管这需要在 Cython 中实现。
Refer to the NumPy documentation at https://numpy.org/devdocs/reference/random/extending.html#new-bit-generators for more information.
在生成随机数据配方中，我们按照 0 和 1 之间的均匀分布生成随机浮点数，但不包括 1。然而，在大多数需要随机数据的情况下，我们需要遵循几种不同的分布中的一种。粗略地说，分布函数是描述随机变量具有低于x的值的概率的函数f（x。实际上，分布描述了随机数据在一定范围内的传播。特别是，如果我们创建一个遵循特定分布的数据直方图，那么它应该大致类似于分布函数的图形。这是最好的例子。
最常见的分布之一是正态分布，它经常出现在统计学中，并构成了许多统计方法的基础，我们将在第 6 章、处理数据和统计学中看到。在此配方中，我们将演示如何按照正态分布生成数据，并绘制此数据的直方图以查看分布的形状。
在生成随机数据配方中，我们从 NumPyrandom模块导入default_rng例程，并使用种子生成器创建一个Generator实例进行演示：
from numpy.random import default_rng
rng = default_rng(12345)像往常一样，我们将 Matplotlibpyplot模块作为plt导入，将 NumPy 作为np导入。
在以下步骤中，我们生成服从正态分布的随机数据：
- 我们在
Generator实例上使用normal方法，根据normal分布生成随机数据。正态分布有两个参数、位置和刻度。还有一个可选的size参数，用于指定生成数据的形状。（参见生成随机数据配方了解size参数的更多信息。）我们生成一个包含 10000 个值的数组，以获得一个大小合理的样本：
mu = 5.0 # mean value
sigma = 3.0 # standard deviation
rands = rng.normal(loc=mu, scale=sigma, size=10000)- 接下来，我们将绘制该数据的直方图。我们增加了直方图中
bins的数量。这并不是绝对必要的，因为默认数字（10）已经足够了，但它确实显示了稍微好一点的分布：
fig, ax = plt.subplots()
ax.hist(rands, bins=20)
ax.set_title("Histogram of normally distributed data")
ax.set_xlabel("Value")
ax.set_ylabel("Density")- 接下来，我们创建一个函数，该函数将生成一系列值的预期密度。这是通过将正态分布的概率密度函数乘以样本数（10000）得出的：
def normal_dist_curve(x):
return 10000*np.exp(-0.5*((x-
mu)/sigma)**2)/(sigma*np.sqrt(2*np.pi))- 最后，我们绘制了数据直方图上的预期分布图：
x_range = np.linspace(-5, 15)
y = normal_dist_curve(x_range)
ax.plot(x_range, y, "k--")结果如图 4.2所示。我们可以在这里看到，采样数据的分布与正态分布曲线的预期分布密切相关：
Figure 4.2: Histogram of data drawn from a normal distribution centered at 5 with a scale of 3, with the expected density overlaid
正态分布具有由以下公式定义的概率密度函数：
这与正态分布函数F（x相关，根据以下公式：
该概率密度函数在与位置参数一致的平均值处达到峰值，并且“钟形”的宽度由比例参数确定。我们可以在图 4.2中看到normal方法在Generator对象上生成的数据直方图与预期分布非常吻合。
Generator类使用 256 步 Ziggurat 方法生成正态分布的随机数据，这比同样在 NumPy 中可用的 Box-Muller 或反向 CDF 实现要快。
正态分布是连续概率分布的一个例子，因为它是为实数定义的，而分布函数是由积分（而不是和）定义的。正态分布（和其他连续概率分布）的一个有趣特征是，选择任意给定实数的概率为 0。这是合理的，因为只有测量在该分布中选择的值位于给定范围内的概率才有意义。（选择特定值的概率不应为零，这是没有意义的。）
正态分布在统计学中很重要，主要是由于*中心极限定理。*粗略地说，该定理表明，具有共同均值和方差的独立同分布（IID）随机变量的和最终类似于具有共同均值和方差的正态分布。不管这些随机变量的实际分布如何，这都是成立的。这允许我们在许多情况下使用基于正态分布的统计测试，即使变量的实际分布不一定是正态的。（然而，在引用中心极限定理时，我们确实需要非常谨慎。）
除了正态分布之外，还有许多其他的连续概率分布。我们已经在 0 到 1 的范围内遇到了均匀分布。更一般地说，在a范围内均匀分布≤ x≤ b的概率密度函数由下式给出：
连续概率密度函数的其他常见示例包括指数分布、β分布和γ分布*。这些分布中的每一个在Generator类上都有相应的方法，该方法从该分布生成随机数据。这些文件通常根据发行版的名称命名，所有文件都使用小写字母。因此，对于上述分布，相应的方法是exponential、beta和gamma。这些分布都有一个或多个参数*，如正态分布的位置和比例，这些参数决定了分布的最终形状。您可能需要查阅 NumPy 文档（https://numpy.org/doc/1.18/reference/random/generator.html#numpy.random.Generator ）或其他来源，查看每个分布需要哪些参数。NumPy 文档还列出了可以生成随机数据的概率分布。
随机过程无处不在。粗略地说，随机过程是一个由相关随机变量组成的系统，通常与时间t 相关≥ 0表示连续随机过程，或通过自然数n=1，2，表示离散随机过程。许多（离散）随机过程满足马尔可夫性质，这使得它们成为马尔可夫链。马尔可夫性是过程是无记忆的陈述，因为只有当前值对下一个值的概率很重要。
在本食谱中，我们将研究一个随机过程的简单示例，该过程模拟了一段时间内到达一个站点的公交车数量。这个过程称为泊松过程。泊松过程N（t有一个参数λ，通常称为强度或速率，以及N（t在给定时间t取N值的概率由以下公式得出：
该方程描述了n辆公交车在t时间到达的概率。从数学上讲，该方程意味着N（t具有参数λt的泊松分布。然而，有一种简单的方法可以通过取服从指数分布的到达时间之和来构造泊松过程。例如，设Xi为i-1-st到达与i*-th 到达之间的时间，与参数λ呈指数分布。现在，我们采用以下等式：*
*
这里，数字N（t）是最大的N，使得t\N<=t。这是我们将在本配方中完成的结构。我们还将通过计算到达时间的平均值来估计过程的强度。
在开始之前，我们从 NumPy 的random模块导入default_rng例程，并创建一个带有种子的新随机数生成器，用于演示：
from numpy.random import default_rng
rng = default_rng(12345)除了随机数生成器之外，我们还将 NumPy 作为np导入，将 Matplotlibpyplot模块作为plt导入。我们还需要提供 SciPy 软件包。
以下步骤显示了如何使用泊松过程对公交车到达进行建模：
- 我们的第一个任务是通过从指数分布中采样数据来创建样本到达时间。NumPy
Generator类上的exponential方法需要一个scale参数，即1/λ，其中λ是速率。我们选择 4 的速率，并创建 50 个样本到达时间：
rate = 4.0
inter_arrival_times = rng.exponential(scale=1./rate, size=50)- 接下来，我们使用 NumPy
add通用函数的accumulate方法计算实际到达时间。我们还创建一个包含整数 0 到 49 的数组，表示每个点的到达数：
arrivals = np.add.accumulate(inter_arrival_times)
count = np.arange(50)- 接下来，我们使用
step绘图方法绘制随时间变化的到达量：
fig1, ax1 = plt.subplots()
ax1.step(arrivals, count, where="post")
ax1.set_xlabel("Time")
ax1.set_ylabel("Number of arrivals")
ax1.set_title("Arrivals over time")结果如图 4.3所示，其中每条水平线的长度代表到达时间：
Figure 4.3: Arrivals over time, where inter-arrival times are exponentially distributed, which makes the number of arrivals at a time a Poisson process
- 接下来，我们定义一个函数，该函数将评估每次计数的概率分布，我们将在这里取为
1。这使用了我们在本配方介绍中给出的泊松分布公式：
from scipy.special import factorial
N = np.arange(15)
def probability(events, time=1, param=rate):
return ((param*time)**events/factorial(events))*np.exp(-
param*time)- 现在，我们绘制单位时间计数的概率分布，因为我们在上一步中选择了
time=1。稍后，我们将在此绘图中添加：
fig2, ax2 = plt.subplots()
ax2.plot(N, probability(N), "k", label="True distribution")
ax2.set_xlabel("Number of arrivals in 1 time unit")
ax2.set_ylabel("Probability")
ax2.set_title("Probability distribution")- 现在，我们继续从样本数据估计速率。我们通过计算到达间隔时间的平均值来实现这一点，对于指数分布，它是标度1/λ的估计器：
estimated_scale = np.mean(inter_arrival_times)
estimated_rate = 1.0/estimated_scale- 最后，我们用每单位时间的计数的估计率绘制概率分布图。我们在*步骤 5:*中生成的真实概率分布的基础上绘制此图
ax2.plot(N, probability(N, param=estimated_rate), "k--", label="Estimated distribution")
ax2.legend()图 4.4给出了结果图，其中我们可以看到，除了一个小差异外，估计的分布非常接近真实分布：
Figure 4.4: Poisson distribution of the number of arrivals per time unit, the true distribution, and the distribution estimated from the sampled data
泊松过程是一种计数过程，如果事件（在时间上）是随机间隔的，且具有固定参数的指数分布，则对一段时间内发生的事件（公交车到达）的数量进行计数。我们按照引言中描述的构造，通过从指数分布中采样到达时间来构造泊松过程。然而，事实证明，当泊松过程以概率的形式给出定义时，这一事实（到达时间呈指数分布）是所有泊松过程的特性。
在这个配方中，我们用给定的rate参数从指数分布中抽取 50 个点。我们必须做一个小的转换，因为从指数分布取样的 NumPyGenerator方法使用了一个相关的scale参数，即1超过rate。一旦我们有了这些点，我们就创建一个数组，其中包含这些指数分布数字的累积和。这就创造了我们的到达时间。实际泊松过程如图 4.3所示，是到达时间与当时发生的相应事件数的组合。
指数分布的平均值（期望值）与尺度参数一致，因此从指数分布中抽取样本的平均值是估计尺度（速率）参数的一种方法。这个估计并不完美，因为我们的样本相对较小。这就是为什么图 4.4中的两个图之间存在微小差异的原因。
有许多类型的随机过程描述各种各样的真实场景。在此配方中，我们使用泊松过程对到达时间进行建模。泊松过程是一个连续的随机过程，这意味着它由一个连续变量t参数化≥ 0，而不是离散变量，n=1,2，…。泊松过程实际上是马尔可夫链，根据马尔可夫链的适当广义定义，也是更新过程的一个例子。续订流程是描述一段时间内发生的事件数量的流程。这里描述的泊松过程是更新过程的一个示例。
许多马尔可夫链除了其定义的马尔可夫性质外，还满足一些性质。例如，如果以下等式适用于所有n、i和j值，则马尔可夫链为齐次：
简单地说，这意味着在一个步骤中从一个状态移动到另一个状态的概率不会随着步骤数的增加而改变。这对于检查马尔可夫链的长期行为非常有用。
构造齐次马尔可夫链的简单例子非常容易。假设我们有两种状态，A和B。在任何给定的步骤中，我们可能处于状态A或状态B。我们根据概率在状态之间移动。例如，假设从A状态过渡到A状态的概率为 0.4，从A状态过渡到B状态的概率为 0.6。同样，假设从B到B的转换概率为 0.2，从B到A的转换概率为 0.8。请注意，在这两种情况下，切换概率加上保持相同和的概率都为 1。在这种情况下，我们可以通过以下等式，以矩阵形式表示从每个状态过渡的概率：
*
该矩阵称为转移矩阵。这里的想法是，通过将包含处于状态a和B（分别为位置 0 和 1）的概率的向量相乘，给出步骤后处于特定状态的概率。例如，如果我们从状态A开始，那么概率向量将包含索引 0 处的 1 和索引 1 处的 0。然后，1 步后处于状态A的概率为 0.4，处于状态B的概率为 0.6。考虑到我们前面概述的可能性，这就是我们所期望的。但是，我们也可以使用矩阵公式编写此计算：
为了在两步之后获得处于任一状态的概率，我们将右侧向量再次乘以转移矩阵T，以获得以下结果：
我们可以继续这个过程无限来获得一系列状态向量，它们构成了我们的马尔可夫链。这种构造可以应用于建模许多简单的现实问题，如果需要，可以使用更多的状态。
贝叶斯概率允许我们通过考虑数据系统地更新我们对情况的理解（在概率意义上）。在更专业的语言中，我们使用数据更新了前分布（我们目前的理解），以获得后分布。例如，在检查浏览网站后继续购买产品的用户比例时，这一点尤其有用。我们从先验信念分布开始。为此，我们将使用贝塔分布，该分布模拟给定成功次数（完成购买）与失败次数（未购买）的成功概率。对于这个配方，我们将假设我们之前的信念是，我们期望 100 个视图中有 25 个成功（75 个失败）。这意味着我们先前的信念遵循贝塔（25，75）分布。假设我们希望计算实际成功率至少为 33%的概率。
我们的方法大致分为三个步骤。我们首先需要了解我们之前对转换率的看法，我们已经确定转换率遵循贝塔（25，75）分布。通过将先验分布的概率密度函数从 0.33 积分到 1，我们计算出转换率至少为 33%的概率。下一步是应用贝叶斯推理，用新信息更新我们的先验信念。然后，我们可以对后验信念进行同样的积分，以检查在给定新信息的情况下，转换率至少为 33%的概率。
在本食谱中，我们将看到如何使用贝叶斯技术，根据我们假设网站的新信息更新先前的信念。
像往常一样，我们需要分别以np和plt的形式导入 NumPy 和 Matplotlib 包。我们还需要 SciPy 软件包，作为sp进口。
以下步骤说明如何使用贝叶斯推理估计和更新转换率估计：
- 第一步是建立先验分布。为此，我们使用来自 SciPy
stats模块的beta分布对象，它有各种方法来处理 beta 分布。我们从别名beta_dist下的stats模块导入beta分布，然后为概率密度函数创建一个方便函数：
from scipy.stats import beta as beta_dist
beta_pdf = beta_dist.pdf- 接下来，我们需要计算在先验信念分布下，成功率至少为 33%的概率。为此，我们使用 SciPy
integrate模块中的quad例程，该例程执行函数的数值积分。我们使用它将在步骤 1中引入的贝塔分布的概率密度函数与我们先前的参数进行积分。我们根据我们之前的分发将概率打印到控制台：
prior_alpha = 25
prior_beta = 75
args = (prior_alpha, prior_beta)
prior_over_33, err = sp.integrate.quad(beta_pdf, 0.33, 1, args=args)
print("Prior probability", prior_over_33)
# 0.037830787030165056- 现在，假设我们已经收到了一些关于新时期成功与失败的信息。例如，在此期间，我们观察到 122 次成功和 257 次失败。我们创建新变量以反映这些值：
observed_successes = 122
observed_failures = 257- 为了获得具有β分布的后验分布的参数值，我们只需将观察到的成功和失败分别添加到
prior_alpha和prior_beta参数中：
posterior_alpha = prior_alpha + observed_successes
posterior_beta = prior_beta + observed_failures- 现在，我们重复我们的数值积分，使用后验分布计算成功率现在高于 33%的概率（使用我们之前计算的新参数）。同样，我们在终端中打印此概率：
args = (posterior_alpha, posterior_beta)
posterior_over_33, err2 = sp.integrate.quad(beta_pdf, 0.33, 1,
args=args)
print("Posterior probability", posterior_over_33)
# 0.13686193416281017- 我们可以在这里看到，根据更新后的后验分布，新的概率是 13%，而不是之前的 3%。这是一个显著的差异，尽管我们仍然不相信在给定这些值的情况下转换率高于 33%。现在，我们绘制先验分布和后验分布图，以可视化概率的增加。首先，我们创建一个值数组，并根据这些值评估概率密度函数：
p = np.linspace(0, 1, 500)
prior_dist = beta_pdf(p, prior_alpha, prior_beta)
posterior_dist = beta_pdf(p, posterior_alpha, posterior_beta)- 最后，我们将步骤 6中计算的两个概率密度函数绘制到一个新的图上：
fig, ax = plt.subplots()
ax.plot(p, prior_dist, "k--", label="Prior")
ax.plot(p, posterior_dist, "k", label="Posterior")
ax.legend()
ax.set_xlabel("Success rate")
ax.set_ylabel("Density")
ax.set_title("Prior and posterior distributions for success rate")结果图如图 4.5所示，我们可以看到后验分布更窄，且中心位于前验分布的右侧：
Figure 4.5: Prior and posterior distributions of a success rate following a beta distribution
贝叶斯技术通过采用先验信念（概率分布）并使用贝叶斯定理将先验信念与给定该先验信念的数据的可能性结合起来，形成后验信念。这实际上与我们在现实生活中可能理解事物的方式相似。例如，当你在某一天醒来时，你可能会相信（从天气预报或其他方面来看），外面有 40%的可能性下雨。打开百叶窗后，你会看到外面多云，这可能表明下雨的可能性更大，因此我们根据这一新数据更新了我们的信念，即 70%的可能性下雨。
为了理解这是如何工作的，我们需要理解条件概率。条件概率是指在另一事件已经发生的情况下，一个事件发生的概率。在符号中，假设事件B已经发生，则事件A的概率写为：
贝叶斯定理是一个强大的工具，可以（象征性地）写成如下：
概率P（A代表我们先前的信念。事件B代表我们收集的数据，因此PBA是我们根据先前的信念产生数据的可能性。概率P**B代表我们的数据产生的概率，PAB代表我们在给定数据的情况下的后验信念。在实践中，概率P（B）可能难以计算或以其他方式估计，因此用比例形式的贝叶斯定理替换上述强等式是很常见的：
在配方中，我们假设我们之前的产品是 beta 版的。贝塔分布的概率密度函数如下式所示：
这里，Γ（α是伽马函数。似然为二元分布，其概率密度函数由以下等式给出：
这里，k是观察次数，j是成功的观察次数之一。在配方中，我们观察到m=122成功和n=257 失败，得出k=m+n=379和j=m=122。为了计算后验分布，我们可以利用贝塔分布是二项式分布的共轭先验这一事实，看到贝叶斯定理比例形式的右侧是贝塔分布，参数为α+m**和β+n这是我们在配方中使用的。贝塔分布是二项式随机变量的共轭先验，这一事实使其在贝叶斯统计中非常有用。**
***我们在本配方中演示的方法是使用贝叶斯方法的一个相当基本的示例，但它对于以系统方式更新给定新数据的先前信念仍然有用。
贝叶斯方法可以用于各种各样的任务，使其成为一个强大的工具。在这个配方中，我们使用了贝叶斯方法，根据我们对网站表现的先前信念和从用户收集的额外数据，对网站的成功率进行建模。这是一个相当复杂的例子，因为我们之前的信念是建立在 beta 分布上的。这里是另一个使用贝叶斯定理的例子，它只使用简单的概率（0 到 1 之间的数字）来检验两个相互竞争的假设。
假设你每天回家时把钥匙放在同一个地方，但有一天早上醒来发现钥匙不在这个地方。在搜索了一小段时间后，你找不到它们，因此得出结论，它们一定已经从存在中消失了。我们把这个假设称为H1。现在，H1肯定解释了数据，D你找不到你的钥匙，因此可能PDH1=1。（如果你的钥匙不存在了，那么你就不可能找到它们。）另一种假设是，你只是在前一天晚上回家时把它们放在别的地方。我们把这个假设称为H2。现在这个假设也解释了数据，所以PDH2=1，但实际上H2远比H1可信。假设你的钥匙完全消失的概率是百万分之一——这是一个巨大的高估，但我们需要保持数字的合理性——而你估计前天晚上你把钥匙放在别处的概率是百分之一。通过计算后验概率，我们得出以下结论：
这突出了一个现实，即你只是把钥匙放错地方的可能性是它们消失的可能性的 10000 倍。果不其然，你很快就会发现你的钥匙已经在口袋里了，因为那天早上早些时候你已经捡到了。
蒙特卡罗方法广泛地描述了使用随机抽样来解决问题的技术。当潜在问题涉及某种不确定性时，这些技术尤其强大。一般方法包括执行大量模拟，每个模拟根据给定的概率分布对不同的输入进行采样，然后将结果聚合，以提供比任何单个样本解更好的真实解近似值。
马尔可夫链蒙特卡罗（MCMC）是一种特定的蒙特卡罗模拟，在这种模拟中，我们构建了一个马尔可夫链，该马尔可夫链对我们所寻求的真实分布具有连续更好的近似性。其工作原理是根据每个阶段仔细选择的接受概率接受或拒绝随机抽样的提议状态，目的是构建一个马尔可夫链，其唯一平稳分布正是我们希望找到的未知分布。
在此配方中，我们将使用 PyMC3 包和 MCMC 方法来估计简单模型的参数。该软件包将处理运行模拟的大部分技术细节，因此我们不需要进一步讨论不同 MCMC 算法实际如何工作的细节。
通常，我们将 NumPy 包和 Matplotlibpyplot模块分别作为np和plt导入。我们还导入并创建一个默认随机数生成器，其中包含一个种子，用于演示，如下所示：
from numpy.random import default_rng
rng = default_rng(12345)我们还需要 SciPy 包中的一个模块用于此配方以及 PyMC3 包，这是一个概率编程包。
执行以下步骤，使用马尔可夫链蒙特卡罗模拟使用样本数据估计简单模型的参数：
- 我们的第一项任务是创建一个函数，该函数表示我们希望识别的底层结构。在这种情况下，我们将估计二次多项式（2 次多项式）的系数。此函数采用两个参数，即范围内的点（固定）和我们希望估计的可变参数：
def underlying(x, params):
return params[0]*x**2 + params[1]*x + params[2]- 接下来，我们设置
true参数和size参数，确定我们生成的样本中有多少点：
size = 100
true_params = [2, -7, 6]- 我们生成用于估计参数的样本。这将包括由我们在步骤 1中定义的
underlying函数生成的基础数据，以及一些遵循正态分布的随机噪声。我们首先生成一系列的x值，这些值将在整个配方中保持不变，然后在我们的随机数生成器上使用underlying函数和normal方法生成样本数据：
x_vals = np.linspace(-5, 5, size)
raw_model = underlying(x_vals, true_params)
noise = rng.normal(loc=0.0, scale=10.0, size=size)
sample = raw_model + noise- 在开始分析之前，最好绘制样本数据，并覆盖底层数据。我们使用
scatter绘图方法仅绘制数据点（无连接线），然后使用虚线绘制底层二次结构：
fig1, ax1 = plt.subplots()
ax1.scatter(x_vals, sample, label="Sampled data")
ax1.plot(x_vals, raw_model, "k--", label="Underlying model")
ax1.set_title("Sampled data")
ax1.set_xlabel("x")
ax1.set_ylabel("y")结果是图 4.6，在图中我们可以看到，即使有噪声，基础模型的形状仍然可见，尽管该模型的确切参数不再明显：
Figure 4.6: Sampled data with the underlying model overlaid
- 我们已经准备好开始分析，因此我们现在导入别名为
pm的 PyMC3 包，如下所示：
import pymc3 as pm- PyMC3 编程的基本对象是
Model类，该类通常使用上下文管理器接口创建。我们还为参数创建了先验分布。在这种情况下，我们将假设我们的先验参数正态分布，平均值为 1，标准偏差为 1。我们需要 3 个参数，所以我们提供了shape参数。Normal类创建随机变量，将用于蒙特卡罗模拟：
with pm.Model() as model:
params = pm.Normal("params", mu=1, sigma=1, shape=3)- 我们为基础数据创建一个模型，可以通过将我们在步骤 6中创建的随机变量
param传递到我们在步骤 1中定义的underlying函数中来实现。我们还创建了一个变量来处理我们的观察结果。为此，我们使用Normal类，因为我们知道噪声正态分布在底层数据y周围。我们将标准偏差设置为2，并将观察到的sample数据传递到observed关键字参数中（这也在Model上下文中）：
y = underlying(x_vals, params)
y_obs = pm.Normal("y_obs", mu=y, sigma=2, observed=sample)- 为了运行模拟，我们只需要在
Model上下文中调用sample例程。我们传递cores参数以加快计算速度，但将所有其他参数保留为默认值：
trace = pm.sample(cores=4)这些模拟应该需要很短的时间来执行。
- 接下来，我们使用 PyMC3 中的
plot_posterior例程绘制后验分布。此例程从执行模拟的采样步骤中获取trace结果。我们预先使用plt.subplots例程创建自己的图形和轴，但这并不是绝对必要的。我们在一个图形上使用了三个子图，我们将Axes的axs2元组传递给ax关键字参数下的绘图路由：
fig2, axs2 = plt.subplots(1, 3, tight_layout=True)
pm.plot_posterior(trace, ax=axs2)结果图如图 4.7所示，您可以看到这些分布中的每一个都近似为正态分布，其平均值与真实参数值相似：
Figure 4.7: Posterior distributions of estimated parameters
- 现在，使用跟踪中的
params项上的mean方法从跟踪中检索每个估计参数的平均值，这只是一个 NumPy 数组。我们传递axis=0参数是因为我们需要参数估计矩阵中每一行的平均值。我们在终端中打印这些估计参数：
estimated_params = trace["params"].mean(axis=0)
print("Estimated parameters", estimated_params)
# Estimated parameters [ 2.03213559 -7.0957161 5.27045299]- 最后，通过将x值和估算参数传递给步骤 1中定义的
underlying函数，我们使用估算参数生成估算基础数据。然后，我们将这些估计的基础数据与真实的基础数据一起绘制在同一轴上：
estimated = underlying(x_vals, estimated_params)
fig3, ax3 = plt.subplots()
ax3.plot(x_vals, raw_model, "k", label="True model")
ax3.plot(x_vals, estimated, "k--", label="Estimated model")
ax3.set_title("Plot of true and estimated models")
ax3.set_xlabel("x")
ax3.set_ylabel("y")
ax3.legend()结果图如图 4.8所示，在该范围内，这两个模型之间只有很小的差异：
Figure 4.8: True model and estimated model plotted on the same axes. There is a small discrepancy between the estimated parameters and the true parameters
在Model上下文管理器中可以找到此配方中有趣的代码部分。此对象跟踪随机变量、协调模拟并跟踪状态。上下文管理器为我们提供了一种方便的方法，将概率变量从周围的代码中分离出来。
我们首先提出一个先验分布，用于表示参数的随机变量的分布，其中有三个。我们提出了正态分布，因为我们知道参数不能偏离值 1 太远。（例如，我们可以通过查看步骤 4中生成的图来判断这一点。）使用正态分布将使接近当前值的值具有更高的概率。接下来，我们添加与观测数据相关的细节，这些细节用于计算用于接受或拒绝状态的接受概率。最后，我们使用sample例程启动采样器。这将构造马尔可夫链并生成所有步骤数据。
sample例程根据要模拟的变量类型设置采样器。由于正态分布是一个连续变量，sample例程选择了无掉头取样器（螺母。这是一个合理的连续变量通用采样器。螺母的一个常见替代品是 Metropolis 取样器，它的可靠性较低，但在某些情况下比螺母更快。PyMC3 文件建议尽可能使用螺母。
采样完成后，我们绘制轨迹的后验分布（马尔可夫链给出的状态），以查看我们生成的近似值的最终形状。我们可以在这里看到，所有三个随机变量（参数）正态分布在大约正确的值附近。
在引擎盖下，PyMC3 使用 Theano 来加速计算。这使得 PyMC3 可以在图形处理单元（GPU）上执行计算，而不是在中央处理单元（CPU）上执行计算，从而大大提高了计算速度。Theano 还支持动态生成 C 代码，以进一步提高计算速度。
蒙特卡罗方法非常灵活，我们这里给出的例子是一个可以使用它的特殊情况。蒙特卡罗方法应用的一个更典型的基本示例是估计积分值，通常是蒙特卡罗积分。蒙特卡罗积分的一个非常有趣的例子是估计π的值≈ 3.1415. 让我们简单地看看它是如何工作的。
首先，我们取单位圆盘，它的半径是 1，因此有一个面积π。我们可以把这个圆盘围成一个正方形，顶点在点（1,1），（-1,1），（1，-1）和（-1，-1）上。此正方形的面积为 4，因为边长度为 2。现在我们可以在这个正方形上均匀地生成随机点。当我们这样做时，这些随机点中任何一个位于给定区域内的概率与该区域的面积成正比。因此，可以通过将区域内随机生成的点的比例乘以正方形的总面积来估计区域的面积。特别是，我们可以通过简单地将磁盘内随机生成的点的数量乘以 4，再除以生成的点的总数来估计磁盘的面积。
我们可以很容易地用 Python 编写一个函数来执行此计算，该函数可能如下所示：
import numpy as np
from numpy.random import default_rng
def estimate_pi(n_points=10000):
rng = default_rng()
points = rng.uniform(-1, 1, size=(2, n_points))
inside = np.less(points[0, :]**2 + points[1, :]**2, 1)
return 4.0*inside.sum() / n_points仅运行此函数一次将给出π的合理近似值：
estimate_pi() # 3.14224我们可以通过使用更多的点来提高估计的准确性，但我们也可以运行多次并平均结果。让我们运行这个模拟 100 次，并对结果进行平均（我们将使用 concurrent futures 对其进行并行化，以便我们可以在需要时运行更多的样本）：
from concurrent.futures import ProcessPoolExecutor, as_completed
from statistics import mean
with ProcessPoolExecutor() as pool:
fts = [pool.submit(estimate_pi) for _ in range(100)]
results = list(ft.result() for ft in as_completed(fts))
print(mean(results))运行此代码一次，将π的估计值打印为 3.1415752，这是对真实值的更好估计。
PyMC3 软件包有许多特性，这些特性通过许多示例（进行了说明 https://docs.pymc.io/ ）。还有另一个基于 TensorFlow 的概率编程库（https://www.tensorflow.org/probability ）。
下面这本书是概率和随机过程的一个很好的综合参考：
- Grimmett，G.和 Stirzaker，D.（2009）。概率与随机过程。第三版牛津：牛津大学出版社*。*
贝叶斯定理和贝叶斯统计的简单介绍如下：
- Kurt，W.（2019）。贝叶斯统计有趣的方式。CA，旧金山：没有淀粉压榨公司。*******