Think Bayes Note 8: 骰子点数之和
1. 点数之和
在《龙与地下城》角色扮演游戏中，玩家通过投掷 3 个 6 面骰子来决定自己角色的属性。求投掷 3 个 6 面骰子点数之和的分布。
这个问题可以从通过模拟或枚举的方法来计算。
1.1. 模拟
可以通过多次模拟投掷 3 个骰子的过程，并汇总结果，得到结果的分布。定义代表单个骰子的 Die 类如下：
class Die(Pmf):
def __init__(self, sides):
super().__init__()
for x in range(1, sides + 1):
self.set(x, 1)
self.normalize()
其中，__init__() 的参数 sides 为骰子的面数，初始化时，设置骰子每个面出现的概率相等。模拟代码为：
d6 = Die(6) dice = [d6] * 3 pmf_sum_sim = sample_sum(dice, 1000, 'sim')
d6 是一个 6 面骰子， dice 是包含 3 个 d6 的列表；sample_sum() 用于从 dice 中进行采样，代码为：
def random_sum(dists):
return sum(dist.random() for dist in dists)
def sample_sum(dists, n, name=''):
values = [random_sum(dists) for _ in range(n)]
pmf = Pmf(name=name)
for value in values:
pmf.incr(value)
pmf.normalize()
return pmf
random_sum() 的参数 dists 是一个分布的列表，它依调用每个分布的 random() 方法并对结果求和，这里相当于计算投掷 3 个 6 面骰子得到的点数之和。sample_sum() 则会对 dists 进行 n 次试验，得到一个结果列表，生成一个 Pmf。
1.2. 枚举
除了像上面一样通过模拟投掷 3 个 6 面骰子的过程来得到骰子点数和的分布，还可以枚举所有可能的结果并计算其概率。代码如下：
pmf_sum_exact = d6 + d6 + d6 pmf_sum_exact.name = 'exact'
Pmf 重载了加号操作符如下：
class Pmf(DfWrapper):
# ...
def __add__(self, other):
pmf = Pmf()
for v1, p1 in self.iter_items():
for v2, p2 in other.iter_items():
pmf.incr(v1 + v2, p1 * p2)
return pmf
# ...
__add__() 会计算自身和 other 两个 Pmf 相加，
使用加号计算 a 和 b 两个 Pmf 相加时，会得到一个新的 Pmf，该 Pmf 中包含 a 和 b 中的值相加的所有可能结果，以及对应的概率。
绘制模拟和枚举的结果：
pmf_sum_sim.plot_with([pmf_sum_exact])
图像为：
2.最大化
假设我们使用这 3 个 6 面骰子连续投掷了 6 轮，记下每一轮中三个骰子点数之和，其最大值的分布是怎样的呢？这个问题依旧可以从模拟和枚举的角度来计算。
2.1. 模拟
模拟的代码如下：
rolls = [pmf_sum_exact] * 6 pmf_max_sim = sample_max(rolls, 1000, 'sim')
这里 rolls 为 6 个 pmf_sum_exact 的列表， pmf_sum_exact 是前面得到的投掷 3 个 6 面骰子的和的分布。sample_max() 从 rolls 中采样最大值，代码为：
def random_max(dists):
return max(dist.random() for dist in dists)
def sample_max(dists, n, name=''):
values = [random_max(dists) for _ in range(n)]
# print(values)
pmf = Pmf(name=name)
for value in values:
pmf.incr(value)
pmf.normalize()
return pmf
这里 sample_max() 和 random_max() 与前面的 sample_sum() 和 random_sum() 类似，只是将用到 sum 的地方换成了 max。
2.2. 枚举
定义 pmf_max() 如下：
def pmf_max(pmf_1, pmf_2, name=''):
name = name if name else pmf_1.name
pmf = Pmf(name=name)
for v1, p1 in pmf_1.iter_items():
for v2, p2 in pmf_2.iter_items():
pmf.incr(max(v1, v2), p1 * p2)
return pmf
pmf_max() 用于求两个 Pmf 中最大值得概率，它会遍历 pmf_1 和 pmf_2 中的所有的条目，取较大的值，并更新其概率。
枚举的过程如下：
pmf_max_exact = pmf_sum_exact.copy(name='exact')
for i in range(5):
pmf_max_exact = pmf_max(pmf_max_exact, pmf_sum_exact)
这里先复制了一份 pmf_sum_exact 做为 pmf_max_exact，然后连续 5 次调用 pmf_max(pmf_max_exact, pmf_sum_exact)，即得到投掷 6 轮后最大值的分布。
绘制模拟和枚举的结果：
pmf_max_sim.plot_with([pmf_max_exact])
图像为：
2.3. 使用 CDF 计算最大值的分布
前面模拟和枚举的方法都涉及大量迭代，计算速度较慢。可以使用 CDF 来加速计算。由 CDF 的定义，有：
\begin{equation}
CDF(x) = P(X \leq x)
\end{equation}
设 $X$ 的 CDF 为 $CDF_X$，$Y$ 的 CDF 为 $CDF_Y$，现从 $X$ 中选出一个值 $x$、从 $Y$ 中选出一个值 $y$，计算最大值 $z = max(x, y)$，对于一个任意值 $z_0$，若要 $z \leq z_0$，则要求 $x \leq z_0$ 且 $y \leq a$。如果 $X$ 和 $Y$ 是独立的，则有：
\begin{equation}
CDF_Z(z) = CDF_X(z)CDF_Y(z)
\end{equation}
如果从 $CDF_X$ 中选出 $k$ 个值 $x_1, x_2,…,x_k$，则 $z = max(x_1, x_2,…,x_k)$ 的 CDF 为：
\begin{equation}
CDF_Z(z) = CDF_X(z)^k
\end{equation}
据此可以在 Cdf 类中实现 max() 方法如下：
class Cdf(DfWrapper):
# ...
def max(self, k):
cdf = self.copy()
cdf.d.prob = self.d.prob ** k
return cdf
# ...
该方法会计算重复进行 $k$ 次试验所得结果的最大值的 CDF。
回到使用这 3 个 6 面骰子投掷 6 轮求最大值分布的问题，代码如下：
cdf_sum_exact = pmf_sum_exact.make_cdf() cdf_best_attr = cdf_sum_exact.max(6) pmf_best_attr = cdf_best_attr.make_pmf() pmf_best_attr.plot()
图像为：
3. 混合分布
在前面的问题中，我们只使用了 6 面骰子，假设我们有一盒骰子，其中有：
- 2 个 4 面骰子
- 3 个 6 面骰子
- 2 个 8 面骰子
- 1 个 12 面骰子
- 1 个 20 面骰子
现在从盒子中随机选择一个投掷，得到的点数的分布是怎样的呢？
可以使用如下代码构造一个 Pmf，包含各种骰子的 Pmf 以及选中它的概率：
pmf_dice = Pmf() pmf_dice.set(Die(4), 2) pmf_dice.set(Die(6), 3) pmf_dice.set(Die(8), 2) pmf_dice.set(Die(12), 1) pmf_dice.set(Die(20), 1) pmf_dice.normalize()
定义 make_mixture() 用于计算混合概率：
def make_mixture(meta_pmf, name='mix'):
mix = Pmf(name=name)
for pmf, p1 in meta_pmf.iter_items():
for x, p2 in pmf.iter_items():
mix.incr(x, p1 * p2)
return mix
从盒子中随机选择一个投掷，得到的点数的分布为：
mix = make_mixture(pmf_dice) mix.plot_bar()
图像为：
