A Rejection Sampler（拒接采样）

“采样”是指从一些概率分布中生成随机数，文中介绍了如何从非标准的概率分布里进行采样，以及如何计算样本在分布里对应的概率。

介绍

通常，在计算机科学和工程中，我们会遇到无法使用方程解决的问题。这些问题通常涉及复杂的系统、嘈杂的输入，或两者兼而有之。以下是一些没有精确解析解决方案的现实世界问题的示例：

1. 您已经建立了飞机的计算机模型，并希望确定飞机在不同天气条件下的承受能力。

2. 您想根据地下水扩散模型确定拟建工厂的化学径流是否会影响附近居民的供水。

3. 你有一个机器人，它从它的相机中捕捉到嘈杂的图像，并且想要恢复这些图像所描绘的物体的三维结构。

4. 您想计算如果您采取特定行动，您在国际象棋中获胜的可能性有多大。

尽管这些类型的问题无法精确解决，但我们通常可以使用称为 Monte Carlo 方法的技术来近似解决它们。在 Monte Carlo 方法中，关键思想是采取许多样本，然后您可以估计解决方案。

蒙特·卡罗方法(Monte Carlo method)也称统计模拟方法，通过重复随机采样模拟对象的概率与统计的问题，在物理、化学、经济学和信息技术领域均具有广泛应用。拒绝采样(reject sampling)就是针对复杂问题的一种随机采样方法。

蒙特卡洛的概念

蒙特卡洛方法是一种近似推断的方法，通过采样大量粒子的方法来求解期望、均值、面积、积分等问题。

蒙特卡洛原来是一个赌场的名称，用它作为名字大概是因为蒙特卡洛方法是一种随机模拟的方法，这很像赌博场里面的扔骰子的过程。最早的蒙特卡洛方法都是为了求解一些不太好求解的求和或者积分问题。

例如下图是一个经典的用蒙特卡洛求圆周率的问题，用计算机在一个正方形之中随机的生成点，计数有多少点落在1/4圆之中，这些点的数目除以总的点数目即圆的面积，根据圆面积公式即可求得圆周率。

接受-拒绝采样

对于累积分布函数未知的分布，我们可以采用接受-拒绝采样。如下图所示，p(z)是我们希望采样的分布，q(z)是我们提议的分布(proposal distribution)，令kq(z)>p(z)，我们首先在kq(z)中按照直接采样的方法采样粒子，接下来判断这个粒子落在途中什么区域，对于落在灰色区域的粒子予以拒绝，落在红线下的粒子接受，最终得到符合p(z)的N个粒子

什么是采样？

术语采样意味着从某种概率分布中生成随机值。例如，您通过滚动六面骰子获得的价值是一个样本。您在洗牌后从牌组顶部抽出的牌是一个样本。飞镖击中棋盘的位置也是一个样本。这些不同样本之间的唯一区别是它们是从不同的概率分布生成的。在骰子的情况下，分布在六个值上放置相等的权重。在卡片的情况下，分布在 52 个值上放置相等的权重。在飞镖板的情况下，分布将重量放置在一个圆形区域（尽管它可能不是均匀分布的，这取决于您作为飞镖玩家的技能）。

我们通常想要使用样本的方式有两种。第一个只是生成一个随机值供以后使用：例如，在扑克的电脑游戏中随机抽牌。使用样本的第二种方式是用于估计。例如，如果您怀疑您的朋友正在玩装满骰子的游戏，您可能需要多次掷骰子以查看某些数字是否比您预期的更频繁。或者，您可能只想描述可能性的范围，如上面的飞机示例。天气是一个相当混乱的系统，这意味着无法准确计算飞机能否在特定的天气情况下幸存下来。相反，您可以多次模拟飞机在许多不同天气条件下的行为，这将让您了解飞机在哪些条件下最有可能发生故障。

使用样本和概率编程

与计算机科学中的大多数应用程序一样，您可以在使用样本和概率进行编程时做出设计决策，这些样本和概率会影响代码的整体清洁度、连贯性和正确性。在本章中，我们将通过一个简单的示例来说明如何在计算机游戏中对随机项目进行采样。特别是，我们将专注于特定于概率的设计决策，包括用于采样和评估概率的函数、使用对数、允许再现性以及将生成样本的过程与特定应用程序分开。

关于符号的简要说明

1. 表明 是值的概率密度函数(PDF)或概率质量函数(PMF)超过值 的随机变量。

2. PDF 是一个连续函数；而 PMF 是一个离散函数，其中 ℤ 是所有整数的集合。

3. 飞镖盘的概率分布是连续的 PDF，而骰子的概率分布是离散的 PMF。在这两种情况下，。

4. 给定一个值 ，我们可能想要评估该位置的概率密度（或质量）是多少。在数学符号中，我们将其写为 。

5. 给定 PDF 或 PMF，我们可能还想采样一个值 以与分布成正比的方式（这样我们更有可能在概率较高的地方获得样本）。在数学符号中，我们将其写为, 表明 采样成正比

采样魔法物品（例子）

假设我们正在编写角色扮演游戏 (RPG)，我们想要一种为怪物随机掉落的魔法物品生成奖励统计数据的方法。我们希望一个项目给予的最大奖励是 +5，并且较高的奖励比较低的奖励的可能性较小。如果 是奖金值的随机变量，则：

我们还可以指定应该在六个属性（敏捷、体质、力量、智力、智慧和魅力）之间分配我们的奖金。因此，具有 +5 奖励的项目可以将这些点数分布在不同的属性上（例如，+2 智慧和 +3 智力）或集中在一个属性中（例如，+5 魅力）。

我们将如何从这个分布中随机抽样？最简单的方法可能是首先对整体项目奖金进行抽样，然后对奖金在统计数据中的分配方式进行抽样。方便的是，奖金的概率分布及其分配方式都是多项式分布的实例。

多项分布

定义：当您有多个可能的结果，并且您想要表示这些结果中的每一个发生的概率。用于解释多项式分布的经典示例是ball 和 urn。这个想法是你有一个装有不同颜色球的骨灰盒（例如，30% 的红色、20% 的蓝色和 50% 的绿色）。你拿出一个球，记录它的颜色，把它放回瓮中，然后重复多次。在这种情况下，结果对应于绘制特定颜色的球，每个结果的概率对应于该颜色球的比例（例如，对于绘制蓝色球的结果，概率为。然后使用多项分布来描述抽取多个球时可能的结果组合（例如，两个绿色和一个蓝色）。

书写 MultinomialDistribution 类

通常，分布有两个用例：我们可能想要从该分布中采样，并且我们可能想要在该分布的 PMF 或 PDF 下评估一个（或多个）样本的概率。虽然执行这两个函数所需的实际计算完全不同，但它们依赖于一条共同的信息：分布的参数是什么。在多项分布的情况下，参数是事件概率， （对应于上述骨灰盒示例中不同颜色球的比例）

最简单的解决方案是简单地创建两个函数，它们都采用相同的参数，但在其他方面是独立的。但是，我通常会选择使用一个类来表示我的分布。这样做有几个好处：

1. 创建类时，您只需传入参数一次。

2. 关于分布，我们可能还想知道其他属性：均值、方差、导数等。一旦我们拥有了一些对公共对象进行操作的函数，使用类比传递函数更方便。许多不同功能的相同参数。

3. 检查参数值是否有效通常是一个好主意（例如，在多项式分布的情况下，向量 事件概率的总和应为 1)。在类的构造函数中执行一次此检查要高效得多，而不是每次调用其中一个函数时执行一次。

4. 有时计算 PMF 或 PDF 涉及计算常数值（给定参数）。使用类，我们可以在构造函数中预先计算这些常量，而不必在每次调用 PMF 或 PDF 函数时计算它们。

在实践中，这是有多少统计包可以工作，包括 SciPy 自己的发行版，它们位于scipy.stats模块中。然而，当我们使用其他 SciPy 函数时，我们没有使用它们的概率分布，这既是为了说明，也是因为 SciPy 中目前没有多项式分布。

这是该类的构造函数代码：

import numpy as np
from scipy.special import gammaln


class MultinomialDistribution(object):

    def __init__(self, p, rso=np.random):
        """初识化多项式类

        参数
        ----------
        p: 存有k的数组，表示总概率
        rso: 随机数的对象，表示随机数构造器

        """

        # 检查总概率是否为1，如果不是，抛出错误
        if not np.isclose(np.sum(p), 1.0):
            raise ValueError("outcome probabilities do not sum to 1")

        # 存储输入的参数
        self.p = p
        self.rso = rso

        # log-PMF，预先计算对数概率,  (函数`np.log` 在数组上p上按元素操作)
        self.logp = np.log(self.p)

在我们进入课程的其余部分之前，让我们回顾一下与构造函数相关的两点。

描述性变量名与数学变量名

通常，鼓励程序员使用描述性变量名称：例如，使用名称independent_variableanddependent_variable，而不是x。一个标准的经验法则是永远不要使用只有一两个字符的变量名。但是，您会注意到，在我们MultinomialDistribution类的构造函数中，我们使用了 的变量名p，这违反了典型的命名约定。

虽然我同意这种命名约定应该适用于几乎所有领域，但有一个例外：数学。编写数学方程的难点在于这些方程的变量名通常只有一个字母：等，所以，如果你直接将它们转换为代码，最简单的变量名是和。显然，这些不是信息量最大的变量名称（不能传达太多信息），但是具有更多描述性的变量名称也会使代码和方程之间的切换变得更加困难。

我认为在编写直接实现方程的代码时，应使用与方程中相同的变量名。这使得很容易看到代码的哪些部分正在实现等式的哪些部分。当然，这会使代码更难孤立地理解，因此注释很好地解释各种计算的目标是特别重要的。如果方程列在学术论文中，注释应引用方程编号，以便于查找。

导入 NumPy

我们将numpy模块导入为np. 这是数值计算领域的标准做法，因为 NumPy 提供了大量有用的函数，其中许多甚至可以在单个文件中使用。

import numpy as np

从多项分布抽样

NumPy 为我们提供了一个函数：np.random.multinomial，我们可以围绕它做出一些设计决策。

播种随机数生成器

尽管我们确实想抽取一个随机样本，但有时我们希望我们的结果是可重复的：即使数字看起来是随机的，如果我们再次运行程序，我们可能希望它使用相同的“随机”数字序列.

为了允许生成这种“可重复的随机”数，我们需要告诉我们的采样函数如何生成随机数。我们可以通过使用 NumPyRandomState对象来实现这一点，它本质上是一个可以传递的随机数生成器对象。它具有大部分相同的功能np.random；不同之处在于我们可以控制随机数的来源。我们创建它如下：

>>> import numpy as np
>>> rso = np.random.RandomState(230489)

其中传递给RandomState构造函数的数字是随机数生成器的种子。只要我们用相同的种子实例化它，一个RandomState对象就会以相同的顺序产生相同的“随机”数，从而确保可复制性：

>>> rso.rand()
0.5356709186237074
>>> rso.rand()
0.6190581888276206
>>> rso.rand()
0.23143573416770336
>>> rso.seed(230489)
>>> rso.rand()
0.5356709186237074
>>> rso.rand()
0.6190581888276206

早些时候，我们看到构造函数接受一个名为 的参数rso。这个rso变量是一个RandomState已经初始化的对象。我喜欢将RandomState对象作为可选参数：偶尔不强制使用它很方便，但我确实希望有使用它的选项（如果我只使用np.random模块，我将无法去做）。

因此，如果rso未给出变量，则构造函数默认为np.random.multinomial。否则，它使用来自RandomState对象本身的多项式采样器。

什么是参数？

    def sample(self, n):
        """从结果概率为“self.p”的多项分布中抽取“n”个事件的样本.

        Parameters
        ----------
        n: 样本数

        Returns
        -------
        存有k的数组,每个结果的抽样出现次数

        NumPy 为我们提供了一个函数：np.random.multinomial，从多项分布抽样

        """
        x = self.rso.multinomial(n, self.p)
        return x

计算多项式 PMF

尽管我们不需要明确计算我们生成的魔法物品的概率，但编写一个可以计算分布的概率质量函数 (PMF) 或概率密度函数 (PDF) 的函数几乎总是一个好主意。为什么？

1. 我们可以用它进行测试：如果我们用我们的采样函数采集很多样本，那么它们应该近似于精确的 PDF 或 PMF。如果经过多次采样后，近似值很差或明显错误，那么我们就知道我们的代码中某处存在错误。
2. 您经常在以后需要它。例如，我们可能希望将随机生成的项目分类为common、uncommon和rare，这取决于它们生成的可能性。为了确定这一点，我们需要能够计算 PMF。

多项式 PMF 方程

形式上，多项分布具有以下等式：

\

当 是一个长度向量 指定每个事件发生的次数，以及是一个向量，指定每个事件发生的概率。

上面方程中的阶乘实际上可以用一个特殊的函数来表示，，称为伽马函数。当我们开始编写代码时，使用伽马函数而不是阶乘会更方便和有效，因此我们将使用：

使用对数值

科普：在计算机里面数据都是以二进制的形式存储的，如果数据超过了计算机所能存储的最大范围，就会发生溢出。传入的值是很小的负数时，就会出现下溢，输出结果就是0。

在进入实现上述等式所需的实际代码之前，我想强调一下编写概率代码时最重要的设计决策之一：使用对数值。这意味着，而不是直接使用概率，我们应该使用对数概率，对数. 这是因为概率会很快变得非常小，从而导致下溢错误。

为了激发这一点，请考虑概率必须介于 0 和 1（含）之间。NumPy 有一个有用的函数 ，finfo它会告诉我们系统浮点值的限制。例如，在 64 位机器上，我们看到可用的最小正数（由 给出tiny）是

>>> import numpy as np
>>> np.finfo(float).tiny
2.2250738585072014e-308

虽然这看起来很小，但遇到这种数量级甚至更小的概率并不罕见。此外，乘以概率是一种常见的操作，但是如果我们尝试以非常小的概率来执行此操作，我们会遇到下溢问题：

>>> tiny = np.finfo(float).tiny
>>> # if we multiply numbers that are too small, we lose all precision
>>> tiny * tiny
0.0

然而，取对数有助于缓解这个问题，因为我们可以用对数表示比通常情况下更广泛的数字。正式地，对数值范围为。实际上，它们的范围从min返回的值（finfo可以表示的最小数字）到零。该min值远小于该值的对数tiny（如果我们不在对数空间中工作，这将是我们的下限）：

>>> # this is our lower bound normally
>>> np.log(tiny)
-708.39641853226408
>>> # this is our lower bound when using logs
>>> np.finfo(float).min
-1.7976931348623157e+308

因此，通过使用对数值，我们可以极大地扩展可表示数字的范围。此外，我们可以使用加法对对数进行乘法，因为 因此，如果我们用对数进行上面的乘法，我们不必担心由于下溢而导致的精度损失：

>>> # the result of multiplying small probabilities
>>> np.log(tiny * tiny)
-inf
>>> # the result of adding small log probabilities
>>> np.log(tiny) + np.log(tiny)
-1416.7928370645282

当然，这个解决方案不是灵丹妙药。如果我们需要从对数导出数字（例如，将概率相加，而不是相乘），那么我们又回到下溢：

>>> tiny*tiny
0.0
>>> np.exp(np.log(tiny) + np.log(tiny))
0.0

尽管如此，使用对数进行所有计算可以省去很多麻烦。如果我们需要回到原始数字，我们可能会被迫失去这种精度，但我们至少保留了一些关于概率的信息——例如，足以比较它们——否则会丢失。

编写 PMF 代码

既然我们已经看到了使用对数的重要性，我们实际上可以编写我们的函数来计算 log-PMF

def log_pmf(self, x):
    """计算多项式的​​对数概率质量函数（log-PMF），其结果概率为“self.p”，用“x”。

    Parameters
    ----------
    x:  存有k的数组,每个结果的抽样出现次数

    Returns
    -------
    The evaluated log-PMF for draw `x`

    """
    # Get the total number of events
    n = np.sum(x)

    # equivalent to log(n!)
    log_n_factorial = gammaln(n + 1)
    # equivalent to log(x1! * ... * xk!)
    sum_log_xi_factorial = np.sum(gammaln(x + 1))

    # 如果self.p的概率为0, 那么self.logp对应的值为-inf
    # 如果初次发生的x事件为0时，log_pi_xi为0 then multiplying them together will give nan, but
    # we want it to just be 0.
    log_pi_xi = self.logp * x
    log_pi_xi[x == 0] = 0
    # equivalent to log(p1^x1 * ... * pk^xk)
    sum_log_pi_xi = np.sum(log_pi_xi)

    # Put it all together
    log_pmf = log_n_factorial - sum_log_xi_factorial + sum_log_pi_xi
    return log_pmf

在大多数情况下，这是上述多项式 PMF 等式的实现。该gammaln函数来自scipy.special，并计算对数伽马函数， 如上所述，使用 gamma 函数比阶乘函数更方便；这是因为 SciPy 给了我们一个 log-gamma 函数，而不是 log-factorial 函数。我们可以自己计算一个对数阶乘，使用类似的东西：

log_n_factorial = np.sum(np.log(np.arange(1, n + 1)))
sum_log_xi_factorial = np.sum([np.sum(np.log(np.arange(1, i + 1))) for i in x])

但如果我们使用 SciPy 内置的伽马函数，它更容易理解、更容易编码并且计算效率更高。

我们需要解决一种边缘情况：当， 那么 这很好，除了当无穷大乘以零时的以下行为：、

>>> # it's fine to multiply infinity by integers...
>>> -np.inf * 2.0
-inf
>>> # ...but things break when we try to multiply by zero
>>> -np.inf * 0.0
nan # 不是一个数字

nan意味着“不是数字”，而且处理起来几乎总是很痛苦，因为大多数计算都会nan产生另一个nan. 所以，如果我们不处理这种情况和 ，我们将得到一个nan. 这将与其他数字相加，产生另一个nan，这是没有用的。为了解决这个问题，我们专门检查以下情况，并设置结果 。

让我们暂时回到我们对使用对数的讨论。即使我们真的只需要 PMF，而不是 log-PMF，通常最好先用对数计算它，然后在需要时取幂。

def pmf(self, x):
    """Evaluates the probability mass function (PMF) of a multinomial
    with outcome probabilities `self.p` for a draw `x`.

    Parameters
    ----------
    x: numpy array of length `k`
        The number of occurrences of each outcome

    Returns
    -------
    The evaluated PMF for draw `x`

    """
    pmf = np.exp(self.log_pmf(x))
    return pmf

测试，执行。

>>> dist = MultinomialDistribution(np.array([0.25, 0.25, 0.25, 0.25]))
>>> dist.log_pmf(np.array([1000, 0, 0, 0])
-1386.2943611198905
>>> dist.log_pmf(np.array([999, 0, 0, 0])
-1384.9080667587707

采样魔法物品，重温

现在我们已经编写了多项函数，我们可以将它们用于生成我们的魔法物品。为此，我们将创建一个名为MagicItemDistribution 的类，位于文件中rpg.py：

lass MagicItemDistribution(object):

    # these are the names (and order) of the stats that all magical
    # items will have 五种属性
    stats_names = ("dexterity", "constitution", "strength",
                   "intelligence", "wisdom", "charisma")

    def __init__(self, bonus_probs, stats_probs, rso=np.random):
        """Initialize a magic item distribution parameterized by `bonus_probs`
        and `stats_probs`.

        Parameters
        ----------
        bonus_probs: numpy array of length m  奖金概率
            The probabilities of the overall bonuses. Each index in
            the array corresponds to the bonus of that amount (e.g.,
            index 0 is +0, index 1 is +1, etc.)

        stats_probs: numpy array of length 6  统计概率
            The probabilities of how the overall bonus is distributed
            among the different stats. `stats_probs[i]` corresponds to
            the probability of giving a bonus point to the ith stat;
            i.e., the value at `MagicItemDistribution.stats_names[i]`.

        rso: numpy RandomState object (default: np.random)  随机数生成器
            The random number generator

        """
        # Create the multinomial distributions we'll be using
        self.bonus_dist = MultinomialDistribution(bonus_probs, rso=rso)
        self.stats_dist = MultinomialDistribution(stats_probs, rso=rso)

我们MagicItemDistribution类的构造函数采用奖励概率、统计概率和随机数生成器的参数。尽管我们在上面指定了我们想要的奖励概率，但将参数编码为传入的参数通常是一个好主意。这为在不同分布下采样项目留下了可能性。（例如，也许奖励概率会随着玩家等级的增加而改变。）我们将统计数据的名称编码为一个类变量stats_names，尽管这很容易成为构造函数的另一个参数。

如前所述，对魔法物品进行抽样有两个步骤：首先对总体奖金进行抽样，然后对各项统计数据中的奖金分布进行抽样。因此，我们将这些步骤编码为两种方法：_sample_bonus和_sample_stats：

def _sample_bonus(self):
    """抽样总奖金的价值.

    Returns
    -------
    integer
        The overall bonus 总奖金

    """
    # The bonus is essentially just a sample from a multinomial
    # distribution with n=1; i.e., only one event occurs.
    sample = self.bonus_dist.sample(1)

    # `sample` is an array of zeros and a single one at the
    # location corresponding to the bonus. We want to convert this
    # one into the actual value of the bonus.
    bonus = np.argmax(sample)
    return bonus

def _sample_stats(self):
    """抽样总奖金及其在不同统计数据中的分布情况.

    Returns
    -------
    numpy array of length 6
        The number of bonus points for each stat

    """
    # First we need to sample the overall bonus
    bonus = self._sample_bonus()

    # Then, we use a different multinomial distribution to sample
    # how that bonus is distributed. The bonus corresponds to the
    # number of events.
    stats = self.stats_dist.sample(bonus)
    return stats

我们本可以将它们变成一个单一的方法——尤其是因为_sample_stats它是唯一依赖的函数_sample_bonus——但我选择将它们分开，既是因为它使采样例程更容易理解，而且因为将其分解成更小的部分使得代码更容易测试。

您还会注意到这些方法以下划线为前缀，表明它们实际上并不打算在类之外使用。相反，我们提供了以下功能sample：

def sample(self):
    """Sample a random magical item.

    Returns
    -------
    dictionary
        The keys are the names of the stats, and the values are
        the bonus conferred to the corresponding stat.

    """
    stats = self._sample_stats()
    item_stats = dict(zip(self.stats_names, stats))
    return item_stats

该sample函数本质上与 相同_sample_stats，除了它返回一个以统计名称作为键的字典。这为采样项目提供了一个干净且易于理解的界面——很明显哪些统计数据有多少奖励点——但它也保持开放的选项，_sample_stats如果一个人需要采集很多样本并且需要效率。

我们使用类似的设计来评估项目的概率。同样，我们公开了高级方法pmf，log_pmf并采用以下形式生成的字典sample：

def log_pmf(self, item):
    """Compute the log probability of the given magical item.

    Parameters
    ----------
    item: dictionary
        The keys are the names of the stats, and the values are
        the bonuses conferred to the corresponding stat.

    Returns
    -------
    float
        The value corresponding to log(p(item))

    """
    # First pull out the bonus points for each stat, in the
    # correct order, then pass that to _stats_log_pmf.
    stats = np.array([item[stat] for stat in self.stats_names])
    log_pmf = self._stats_log_pmf(stats)
    return log_pmf

def pmf(self, item):
    """Compute the probability the given magical item.

    Parameters
    ----------
    item: dictionary
        The keys are the names of the stats, and the values are
        the bonus conferred to the corresponding stat.

    Returns
    -------
    float
        The value corresponding to p(item)

    """
    return np.exp(self.log_pmf(item))

这些方法依赖于_stats_log_pmf，它计算统计数据的概率（但它采用数组而不是字典）：

def _stats_log_pmf(self, stats):
    """Evaluate the log-PMF for the given distribution of bonus points
    across the different stats.

    Parameters
    ----------
    stats: numpy array of length 6
        The distribution of bonus points across the stats

    Returns
    -------
    float
        The value corresponding to log(p(stats))

    """
    # There are never any leftover bonus points, so the sum of the
    # stats gives us the total bonus.
    total_bonus = np.sum(stats)

    # First calculate the probability of the total bonus
    logp_bonus = self._bonus_log_pmf(total_bonus)

    # Then calculate the probability of the stats
    logp_stats = self.stats_dist.log_pmf(stats)

    # Then multiply them together (using addition, because we are
    # working with logs)
    log_pmf = logp_bonus + logp_stats
    return log_pmf

_stats_log_pmf反过来，方法依赖于_bonus_log_pmf，它计算整体奖金的概率：

def _bonus_log_pmf(self, bonus):
    """Evaluate the log-PMF for the given bonus.

    Parameters
    ----------
    bonus: integer
        The total bonus.

    Returns
    -------
    float
        The value corresponding to log(p(bonus))

    """
    # Make sure the value that is passed in is within the
    # appropriate bounds
    if bonus < 0 or bonus >= len(self.bonus_dist.p):
        return -np.inf

    # Convert the scalar bonus value into a vector of event
    # occurrences
    x = np.zeros(len(self.bonus_dist.p))
    x[bonus] = 1

    return self.bonus_dist.log_pmf(x)

我们现在可以创建我们的分布如下：

>>> import numpy as np
>>> from rpg import MagicItemDistribution
>>> bonus_probs = np.array([0.0, 0.55, 0.25, 0.12, 0.06, 0.02])
>>> stats_probs = np.ones(6) / 6.0
>>> rso = np.random.RandomState(234892)
>>> item_dist = MagicItemDistribution(bonus_probs, stats_probs, rso=rso)

创建后，我们可以使用它来生成一些不同的项目：

>>> item_dist.sample()
{'dexterity': 0, 'strength': 0, 'constitution': 0, 
 'intelligence': 0, 'wisdom': 0, 'charisma': 1}
>>> item_dist.sample()
{'dexterity': 0, 'strength': 0, 'constitution': 1, 
 'intelligence': 0, 'wisdom': 2, 'charisma': 0}
>>> item_dist.sample()
{'dexterity': 1, 'strength': 0, 'constitution': 1, 
 'intelligence': 0, 'wisdom': 0, 'charisma': 0}

而且，如果我们愿意，我们可以评估抽样项目的概率：

>>> item = item_dist.sample()
>>> item
{'dexterity': 0, 'strength': 0, 'constitution': 0, 
 'intelligence': 0, 'wisdom': 2, 'charisma': 0}
>>> item_dist.log_pmf(item)
-4.9698132995760007
>>> item_dist.pmf(item)
0.0069444444444444441

估计攻击伤害

我们已经看到了采样的一种应用：生成怪物掉落的随机物品。我之前提到过，当您想从整个分布中估计某些内容时，也可以使用抽样，当然在某些情况下我们可以使用我们的方法MagicItemDistribution来做到这一点。例如，假设我们的 RPG 中的伤害是通过滚动一定数量的 D12（十二面骰子）来实现的。默认情况下，玩家可以掷一个骰子，然后根据他们的力量奖励添加骰子。因此，例如，如果他们有 +2 的力量加值，他们可以掷三个骰子。造成的伤害是骰子的总和。

我们可能想知道玩家在找到一定数量的武器后可能造成多少伤害；例如，作为设置怪物难度的一个因素。假设收集两个物品后，我们希望玩家能够在大约 50% 的战斗中在三击内击败怪物。怪物应该有多少生命值？

回答这个问题的一种方法是通过抽样。我们可以使用以下方案：

1. 随机选择一件魔法物品。

2. 根据物品的奖励，计算攻击时将掷出的骰子数。

3. 根据将要掷出的骰子数量，为超过 3 次命中造成的伤害生成一个样本。

4. 多次重复步骤 1-3。这将导致对损伤分布的近似。

实施损害分配

class DamageDistribution(object):

    def __init__(self, num_items, item_dist,
                 num_dice_sides=12, num_hits=1, rso=np.random):
        """Initialize a distribution over attack damage. This object can
        sample possible values for the attack damage dealt over
        `num_hits` hits when the player has `num_items` items, and
        where attack damage is computed by rolling dice with
        `num_dice_sides` sides.

        Parameters
        ----------
        num_items: int
            The number of items the player has.
        item_dist: MagicItemDistribution object
            The distribution over magic items.
        num_dice_sides: int (default: 12)
            The number of sides on each die.
        num_hits: int (default: 1)
            The number of hits across which we want to calculate damage.
        rso: numpy RandomState object (default: np.random)
            The random number generator

        """

        # This is an array of integers corresponding to the sides of a
        # single die.
        self.dice_sides = np.arange(1, num_dice_sides + 1)

        # Create a multinomial distribution corresponding to one of
        # these dice.  Each side has equal probabilities.
        self.dice_dist = MultinomialDistribution(
            np.ones(num_dice_sides) / float(num_dice_sides), rso=rso)

        self.num_hits = num_hits
        self.num_items = num_items
        self.item_dist = item_dist

    def sample(self):
        """Sample the attack damage.

        Returns
        -------
        int
            The sampled damage

        """
        # First, we need to randomly generate items (the number of
        # which was passed into the constructor).
        items = [self.item_dist.sample() for i in xrange(self.num_items)]

        # Based on the item stats (in particular, strength), compute
        # the number of dice we get to roll.
        num_dice = 1 + np.sum([item['strength'] for item in items])

        # Roll the dice and compute the resulting damage.
        dice_rolls = self.dice_dist.sample(self.num_hits * num_dice)
        damage = np.sum(self.dice_sides * dice_rolls)
        return damage

构造函数将骰子的面数、我们想要计算的伤害次数、玩家拥有的物品数量、魔法物品的分布（类型为MagicItemDistribution）和随机状态对象作为参数。默认情况下，我们设置num_dice_sides为 12，因为虽然它在技术上是一个参数，但不太可能改变。类似地，我们将num_hits默认值设置为 1，因为更可能的用例是我们只想对单次命中的伤害进行一次采样。

然后我们在sample. （注意与的结构相似性MagicItemDistribution。）首先，我们生成一组玩家可能拥有的魔法物品。然后，我们查看这些项目的强度统计数据，并从中计算要掷骰子的数量。最后，我们掷骰子（再次依赖于我们可信赖的多项式函数）并计算由此造成的伤害。

评估概率发生了什么？

您可能已经注意到，我们没有包括log_pmf或pmf在我们的功能DamageDistribution。这是因为我们实际上不知道 PMF 应该是什么！这将是等式

这个等式说的是我们需要计算每个可能的伤害量的概率，给定每一个可能的集合。我们实际上可以通过蛮力来计算它，但它不会很漂亮。这实际上是我们想要使用采样来逼近我们无法精确计算（或者很难精确计算）的问题的解决方案的完美示例。因此，我们将在下一节中展示如何使用许多样本来近似分布，而不是为 PMF 提供一种方法。

近似分布

现在我们有了机制来回答我们之前的问题：如果玩家有两个物品，并且我们希望玩家能够在 50% 的概率内在三z招内击败怪物，那么怪物应该有多少生命值？

首先，我们创建我们的分布对象，使用item_dist与rso我们之前创建的相同的对象：

>>> from rpg import DamageDistribution
>>> damage_dist = DamageDistribution(2, item_dist, num_hits=3, rso=rso)

现在我们可以绘制一堆样本，并计算第 50 个百分位数（大于 50% 样本的损坏值）：

>>> samples = np.array([damage_dist.sample() for i in xrange(100000)])
>>> samples.min()
3
>>> samples.max()
154
>>> np.percentile(samples, 50)
27.0

玩家可能造成的伤害范围相当广泛，但它有一条长尾：第 50 个百分位数为 27 点，这意味着在一半的样本中，玩家造成的伤害不超过 27 点。因此，如果我们想使用这个标准来设置怪物难度，我们会给他们 27 点生命值。

概括

在本章中，我们已经看到了如何编写代码来从非标准概率分布中生成样本，以及如何计算这些样本的概率。在研究这个例子时，我们已经涵盖了几个适用于一般情况的设计决策：

1. 使用类表示概率分布，并包括用于采样和评估 PMF（或 PDF）的函数。

2. 使用对数计算 PMF（或 PDF）。

3. 从随机数生成器对象生成样本以实现可重复的随机性。

4. 编写输入/输出清晰易懂的函数（例如，使用字典作为的输出MagicItemDistribution.sample），同时仍然暴露这些函数的不太清楚但更高效且纯数字的版本（例如，MagicItemDistribution._sample_stats）。