为有限离散分布建立 Alias 抽样器并核对频率

0310

如果模拟中的取值只有有限个，而同一个分布需要被反复抽样，可以先把分布整理成概率向量，再用 SciPy 的 DiscreteAliasUrn 建立抽样器。建立表格需要时间和内存，但表建好以后，逐个产生样本的成本很低。这正是 Discrete Alias Urn（DAU，离散 Alias-Urn 方法）的使用场景。

本文依据 SciPy 官方 DAU 教程完整译写，并用同版 API 文档核对参数。检索时页面标为 SciPy 1.18.0；原页未给出首次发表日期，也没有单独的个人作者署名，署名保留为 SciPy 文档贡献者。示例需要 NumPy、SciPy 和 Matplotlib。下列代码只经过静态审核，本次未执行抽样或绘图。

DAU 抽样流程示意：有限质量向量经一次建表后反复生成样本，再以正确类别边界核对频率；右下列出平方权重的归一化规则
编辑绘制的流程示意图。图中公式是理论分布，不是本次运行得到的样本或性能结果。

先判断分布是否适合 DAU

DAU 接受两种输入：有限长度的概率向量（probability vector，简称 PV），或者带有有限定义域的概率质量函数（PMF）。它适用于任意形状的有限离散分布，而不是只适用于某一种命名分布。原算法基于 A. J. Walker 的方法，需要一张大小至少为 N 的表，N 是概率向量长度。教程描述每次产生随机变量时需要一个随机数和一次比较；建表时间为 O(N)。因此，初始化相对较慢，抽样很快，表格占用也随分布大小增长。

这意味着“一次建表，多次使用”比“每产生一个样本就重新构造抽样器”更符合它的设计。有限支持仍可能很大，不能把快速抽样理解为初始化没有成本。

从三个类别的概率向量开始

下面的向量对应三个取值，质量分别为 0.18、0.02 和 0.8。先创建 NumPy 随机数生成器，再把它传给 DAU：

import numpy as np
from scipy.stats.sampling import DiscreteAliasUrn

pv = [0.18, 0.02, 0.8]
urng = np.random.default_rng()
rng = DiscreteAliasUrn(pv, random_state=urng)
sample = rng.rvs()

默认的取值从 0 开始，所以这里返回的是 0、1 或 2，而不是向量中的概率本身。原文的一次交互输出是 0，但它只是可能结果；每次运行可以不同。random_state 在这里接收的是已有的 Generator 对象。

作为补充的输入检查，API 要求概率向量是一维、非负、有限的浮点数序列，不能含有 NaN 或无穷值，且至少有一个非零元素。输入可以是未归一化的非负权重，算法会据此建立分布；绘制理论 PMF 时则要自行归一化。

用 domain 平移取值，不要把它误当成类别数量

如果实际取值应为 10、11、12，可传入下面的定义域：

rng = DiscreteAliasUrn(pv, domain=(10, 13), random_state=urng)
sample = rng.rvs()

原文的示例输出是 12，同样可能变化。对于概率向量，domain[0] 只是将索引起点从 0 平移到指定整数；实际取值是从起点到“起点加向量长度减一”。此时 domain[1] 被忽略。上例只会产生 10、11、12，不会因右端写了 13 就多出第四个类别。

如果类别原本是名称或不连续整数，应在抽样后把返回的连续索引映射回自己的标签；domain 的这个参数用途不是任意标签映射。

只有质量函数时，要同时提供有限支持

不必先手工写出向量。也可以给出带 pmf 方法的对象，再以显式 domain 参数或者对象的 support 方法说明有限范围。教程采用后者：

class Distribution:
    def __init__(self, c):
        self.c = c

    def pmf(self, x):
        return x**self.c

    def support(self):
        return (0, 10)

dist = Distribution(2)
rng = DiscreteAliasUrn(dist, random_state=urng)
sample = rng.rvs()

这个例子把 0 到 10 的整数作为支持，端点都包含在内。c=2 时，未归一化的质量是 x²：0 的质量为 0，10 的质量为 100，其余整数按平方增加。原文的一次输出是 10，但不是保证返回值。pmf 在此返回的是相对质量；它不需要预先把总和调整成 1。

注意这里与传入 PV 时的区别：用 PMF 建表必须知道在哪些整数上计算质量，support() 的上界 10 有实际意义。把两种输入的 domain 解释混为一谈，会产生漏掉端点或者错误平移的问题。

抽取一千个样本，与理论质量比较

下面保留原文完整的绘图流程，去掉交互提示符，整理为可阅读的脚本。这里重新创建随机数生成器，因此不会延续上一段已经消耗的随机数流。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats.sampling import DiscreteAliasUrn

class Distribution:
    def __init__(self, c):
        self.c = c

    def pmf(self, x):
        return x**self.c

    def support(self):
        return (0, 10)

dist = Distribution(2)
urng = np.random.default_rng()
rng = DiscreteAliasUrn(dist, random_state=urng)
rvs = rng.rvs(1000)

fig = plt.figure()
ax = fig.add_subplot(111)
x = np.arange(1, 11)
fx = dist.pmf(x)
fx = fx / fx.sum()
ax.plot(x, fx, 'bo', label='true distribution')
ax.vlines(x, 0, fx, lw=2)
ax.hist(
    rvs,
    bins=np.r_[x, 11] - 0.5,
    density=True,
    alpha=0.5,
    color='r',
    label='samples',
)
ax.set_xlabel('x')
ax.set_ylabel('PMF(x)')
ax.set_title('Discrete Alias Urn Samples')
plt.legend()
plt.show()

这张图比较的是两种量：蓝色点和竖线表示归一化后的理论质量；红色直方图表示样本的经验频率。原文配图展示了这样的比较，但本文没有重新运行代码，因而没有把自绘配图伪装成一次抽样结果。

x = np.arange(1, 11) 给出 1 到 10，尽管支持包含 0，这里仍不丢失正质量，因为 0²=0。平方和为 385，所以理论概率是 x²/385。fx / fx.sum() 这一行不能省略，否则会拿原始权重与概率比较。

直方图边界从 0.5 到 10.5，相邻边界间隔为 1，恰好让每个整数处于一个箱子的中央。对于宽度为 1 的箱子，density=True 的高度与相对频率数值一致。如果任意改变箱宽，再把密度高度当作每个类别的概率，就会失去这种对应关系。

一千次抽样的频率一般不会精确等于理论值，小概率类别的波动尤其明显。图形核对可以发现明显错位、漏项或未归一化，不是对随机数质量或全部算法正确性的证明。原文没有固定随机种子；需要比较两次修改时，可自行把 default_rng() 改为带种子的调用，同时记录依赖版本。这是便于复核的编辑建议，不是本次实测结果。

SciPy 原教程1000个DAU样本的红色经验直方图与蓝色理论质量x平方除以385比较。
SciPy 原教程1000个DAU样本的红色经验直方图与蓝色理论质量x平方除以385比较。 来源：SciPy 文档贡献者；保留原文示例，未重新运行生成。

已有向量化 PMF 时，先算完整向量

教程说明 DAU 面向标量 PMF 回调，内部会用 np.vectorize 对支持上的点逐个计算。教程的说明把签名写为浮点参数，同版 API 则把离散自变量写为整数；实际编写时应保证函数能接受支持中的单个整数并返回标量质量。

如果 PMF 本身已经能够一次处理整段数组，先向量化求值，再把结果作为 PV 传入，通常可以避免反复经过标量回调。SciPy 内置离散分布的 pmf 就支持数组输入。原文以二项分布为例：

from scipy.stats import binom
from scipy.stats.sampling import DiscreteAliasUrn

dist = binom(10, 0.2)
domain = dist.support()
x = np.arange(domain[0], domain[1] + 1)
pv = dist.pmf(x)
rng = DiscreteAliasUrn(pv, domain=domain)

domain[1] + 1 是为了让 arange 包含支持的最后一个整数。传入 PV 时仍应保留 domain，这样索引与真实取值的位置保持对应。本例的起点正好是 0，平移没有可见效果，但这个写法也适用于起点不是 0 的有限分布。

urn_factor：用更多表空间换取抽样速度

urn_factor 控制表格大小相对于概率向量长度的倍数，默认是 1，不能小于 1。原文给出把表格扩大为两倍的用法：

urn_factor = 2
rng = DiscreteAliasUrn(pv, urn_factor=urn_factor, random_state=urng)
sample = rng.rvs()

该段延续上一节二项分布的 pv，它从 0 开始，所以省略 domain 不会改变取值位置。如果换成非零起点的向量，应继续传入相应的 domain。原文可能返回 2，实际结果可以不同。

更大的表可能略微改善抽样性能，但会增加初始化时间与内存需求。教程和 API 均建议将这个参数保持在 2 以下；示例取 2 是为了演示参数含义，不等于建议把 2 作为通用默认值。本文没有测量速度，也不提供性能提升比例。实际选择应同时考虑支持长度、总抽样次数和可用内存。

静态检查与适用限制

本文核对了导入、类方法、有限支持、归一化、直方图边界以及 PV 平移规则。示例不含外部数据下载、命令执行、明文密钥或网络连接；但很大的支持会导致建表资源消耗，负权重、非有限值和全零向量会破坏输入条件。Python 与 NumPy 的一般随机数生成器也不应被当作密码学用途的安全随机源。

本次未安装依赖、未执行上述代码、未生成随机样本、未验证跨版本可重现性。这些静态结论不能替代目标环境中的运行检查，也不是“没有漏洞”的承诺。

来源与参考

原文版权声明：© Copyright 2008, The SciPy community；原文版权归 SciPy community 及贡献者；代码相关许可见 SciPy v1.18.0 LICENSE.txt。中文译写及原创示意图由未完纪整理制作，保留原文署名、来源及许可链接。

上游代码版权与完整许可

以下保留对应上游版本的完整声明；第三方例外与附加通知按原文保留。

Copyright (c) 2001-2002 Enthought, Inc. 2003, SciPy Developers.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:

1. Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above
   copyright notice, this list of conditions and the following
   disclaimer in the documentation and/or other materials provided
   with the distribution.

3. Neither the name of the copyright holder nor the names of its
   contributors may be used to endorse or promote products derived
   from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
© 版权声明
文章版权归作者所有，未经允许请勿转载。
THE END
技术译文
# SciPy# NumPy# 随机抽样# 离散分布# Alias# 统计模拟
喜欢就支持一下吧
点赞0 分享
彩雨伞的头像-未完纪钻石会员
Self-introduction-未完纪
Self-introduction
Self-introduction
7天前 72
使用宝塔为 AList 配置反向代理与 HTTPS-未完纪
使用宝塔为 AList 配置反向代理与 HTTPS
使用宝塔为 AList 配置反向代理与 HTTPS
前天 59
新一代垃圾回收器 ZGC 的探索与实践：重读美团 JDK 11 调优案例-未完纪
新一代垃圾回收器 ZGC 的探索与实践：重读美团 JDK 11 调优案例
新一代垃圾回收器 ZGC 的探索与实践：重读美团 JDK 11 调优案例
5小时前 57
WordPress 接入 Microsoft Clarity Cookie 同意与偏好设置-未完纪
WordPress 接入 Microsoft Clarity Cookie 同意与偏好设置
WordPress 接入 Microsoft Clarity Cookie 同意与偏好设置
前天 56
注册区块变体，同时在插入器中隐藏它-未完纪
注册区块变体，同时在插入器中隐藏它
注册区块变体，同时在插入器中隐藏它
前天 56
AingDesk：在电脑上部署 DeepSeek-未完纪
AingDesk：在电脑上部署 DeepSeek
AingDesk：在电脑上部署 DeepSeek
2天前 56
相关推荐
新一代垃圾回收器 ZGC 的探索与实践：重读美团 JDK 11 调优案例-未完纪
新一代垃圾回收器 ZGC 的探索与实践：重读美团 JDK 11 调优案例
新一代垃圾回收器 ZGC 的探索与实践：重读美团 JDK 11 调优案例
5小时前 57
注册区块变体，同时在插入器中隐藏它-未完纪
注册区块变体，同时在插入器中隐藏它
注册区块变体，同时在插入器中隐藏它
前天 56
用内部审计器测试数百个网页的 Cookie 行为-未完纪
用内部审计器测试数百个网页的 Cookie 行为
用内部审计器测试数百个网页的 Cookie 行为
昨天 55
从用户序列到扩展规律：Meta 广告排序的多阶段架构-未完纪
从用户序列到扩展规律：Meta 广告排序的多阶段架构
从用户序列到扩展规律：Meta 广告排序的多阶段架构
昨天 55
用 Ansible 定制 Loki 部署-未完纪
用 Ansible 定制 Loki 部署
用 Ansible 定制 Loki 部署
昨天 55
用维度表加速字符串聚合-未完纪
用维度表加速字符串聚合
用维度表加速字符串聚合
前天 54
评论 抢沙发

请登录后发表评论

    暂无评论内容