Think Bayes Note 8: 骰子点数之和

Author: nex3z 2018-07-29

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()

图像为：