逆变换采样：从均匀分布到任意分布的通用随机数生成方法

做蒙特卡洛模拟的朋友肯定都遇到过这个问题我需要生成一批服从指数分布、韦布尔分布或者其他自定义分布的随机数但手头的编程语言翻来翻去只给我一个uniform(0,1)均匀分布随机数生成器。逆变换采样就是这种情况下的第一选择它不依赖任何分布的近似公式只要你能写出累积分布函数的反函数哪怕这个函数是分段定义的都能老老实实把均匀随机数变换成目标分布。这篇东西我主要写给三类人刚入门仿真、正在写强化学习环境或者排队系统的人还有工作中需要造随机数的数据分析师。内容从最基本的思想讲起到具体的代码实现和我在实际项目中踩过的坑。如果你只是想快速抄一段代码跳到第4节但如果你想搞得明白为什么这招好用、什么时候会失效建议整篇读完。1. 逆变换采样到底解决了什么问题1.1 随机数生成的最底层只有均匀分布先说个很多人忽略的事实计算机里几乎所有的随机数都是先从均匀分布来的。rand()、numpy.random.random()这些函数生成的是 $[0,1)$ 区间上的均匀分布随机数底层是线性同余生成器、梅森旋转算法或者PCG这些伪随机数算法。其他所有分布——正态、指数、泊松、Gamma——本质上都是在这个均匀随机数的基础上施加某种数学变换得来的。这就带来一个核心需求给我一个从均匀分布到目标分布的通用通道。没有这个通道你想做一点超出标准库范围的分布比如双指数分布、Gumbel分布或者某个业务场景里自定义的经验分布就只能在网上搜“有没有现成的包”。逆变换采样就是这条最通用的通道。1.2 为什么不能用均匀分布直接硬凑有一种常见的错误思路既然目标是某个形状的分布那我直接把均匀随机数乘个系数、加个偏移不行吗在连续分布上基本不行。乘系数和加偏移对应的只是线性变换它改变的只是位置的尺度和位置改不了分布的形状。比如说我想生成一个服从指数分布的随机数均匀分布的形状是一条平线指数分布的形状是一条快速衰减的曲线两者之间不是简单的平移缩放关系。这就需要一个能“重塑形状”的变换而逆变换采样做的就是这件事它把均匀分布这个“原材料”通过累积分布函数的反函数重新映射成目标分布的形状。1.3 逆变换采样适用的场景清单根据我的实际经验逆变换采样最适合以下场景CDF的反函数能写出解析表达式或者至少能通过数值方法求出来需要大批量生成随机样本性能敏感分布的CDF比PDF更容易获得比如生存分析里常见的分布需要生成服从经验分布的样本直接用原始数据的累积分布函数一句话总结它是一个在“知道CDF”和“能求出反函数”这两个条件下效率最高、实现最直接的采样方法。2. 核心原理为什么取CDF的反函数就能采样2.1 一个看似神奇但证明很短的定理逆变换采样的数学核心是一条非常简洁的定理设随机变量 $X$ 的累积分布函数为 $F(x)$$U \sim \text{Uniform}(0,1)$。如果记 $F^{-1}(u)$ 为 $F$ 的反函数那么 $X F^{-1}(U)$ 的分布和 $X$ 完全一样。证明只有三步$$ P(X \le x) P(F^{-1}(U) \le x) $$因为 $F^{-1}$ 是单调递增函数所以括号里的不等式可以等价地“作用” $F$ 两边$$ P(U \le F(x)) $$$U$ 是区间 $(0,1)$ 上的均匀分布所以落在任何区间 $(0, a)$ 内的概率就等于 $a$于是$$ F(x) $$这个结果正是 $X$ 的累积分布函数。也就是说$X$ 和 $X$ 同分布。整个推导不涉及复杂的积分逻辑上绕了一个弯但每一步都非常干净。2.2 从概率积分变换的角度找直觉想理解这个定理建议和“概率积分变换”合在一起看。概率积分变换说的是另一件事如果 $X$ 有连续的CDF $F$那么 $F(X)$ 服从 $[0,1]$ 上的均匀分布。这两个定理互为逆操作。概率积分变换是把任意分布“压”回均匀分布逆变换采样是把均匀分布“拉”回任意分布。打个比方概率积分变换就像把一张不规则形状的地图按某种规则摊平成标准尺寸逆变换采样就是把这个摊平过程倒过来按同样的规则把平面上的点映射回原来的地形。2.3 处理CDF不严格单调的情况广义逆函数连续分布里CDF一般是严格单调的直接用反函数没问题。但离散分布或者某些混合分布里CDF是阶梯状的碰到这种情况就必须用广义逆函数$$ F^{-1}(u) \inf{ x : F(x) \ge u } $$翻译成人话从所有“CDF值大于等于 $u$”的点里取最小的那个 $x$。注意这里是“大于等于”不是“大于”。实际写代码的时候这个区别会造成索引偏移如果你用二分查找在CDF数组里搜一定要搜 u的第一个位置而不是 u的。这个细节我后面还会再提一次因为真踩过坑。3. 连续分布实战解析解与代码实现3.1 最经典的例子指数分布指数分布的CDF是$$ F(x) 1 - e^{-\lambda x}, \quad x \ge 0 $$反解 $x$ 的过程很常规$$ u 1 - e^{-\lambda x} $$$$ e^{-\lambda x} 1 - u $$$$ x -\frac{\ln(1-u)}{\lambda} $$直接按这个式子写代码即可import numpy as np def inverse_exp_sample(size, lam1.0, rngNone): if rng is None: rng np.random.default_rng() u rng.random(size) return -np.log(1 - u) / lam这里有个常见的简化因为 $1-U$ 和 $U$ 服从完全相同的均匀分布所以很多人直接写成-np.log(u) / lam。数学上完全没问题而且能少一次减法运算。但我个人习惯于保留1 - u的写法原因有两个第一它保留了公式推导的原始脉络别人看代码时更容易理解是怎么来的第二在某些数值边界情况下面会讲下1-u的行为更可控。当然这属于风格问题没有绝对的对错能跑就行。3.2 韦布尔分布、柯西分布套路完全一样知道套路之后你会发现连续分布的逆变换采样本质上就是查CDF公式、求反函数、写代码这三步。韦布尔分布的CDF是 $F(x) 1 - e^{-(x/\lambda)^k}$反解得到def inverse_weibull_sample(size, lam1.0, k1.0, rngNone): if rng is None: rng np.random.default_rng() u rng.random(size) return lam * (-np.log(1 - u)) ** (1.0 / k)柯西分布的CDF是 $F(x) \frac{1}{\pi}\arctan(x) 0.5$反解得到 $x \tan\big(\pi(u - 0.5)\big)$def inverse_cauchy_sample(size, rngNone): if rng is None: rng np.random.default_rng() u rng.random(size) return np.tan(np.pi * (u - 0.5))这三个例子覆盖了三种典型形态单调衰减、带形状参数、厚尾分布。代码风格上我把rng作为参数传进去而不是直接在函数里调用全局的np.random这是因为现代 NumPy 推荐使用default_rng()创建独立的随机数生成器在并行计算和多线程场景里这样做能避免全局状态污染。这是我后来做大规模仿真时才意识到的细节决定成败。3.3 碰到没有解析反函数的CDF怎么办不是所有连续分布的CDF都有漂亮的解析反函数。最典型的就是正态分布它的CDF本身是个误差函数反函数没法写成初等表达式。这时候有几个选择查表法预先算出CDF的数值表在表上做插值用近似公式比如正态分布CDF反函数的Beasley-Springer-Moro近似换用Box-Muller采样专门为正态分布设计的变换方法数值插值的思路在实战里非常有用。有时候我拿到的是一个分布的历史样本不知道它属于哪个已知分布族这时候可以直接把样本的CDF经验估计出来然后在经验CDF上做线性插值再对这个插值函数做逆变换。这种方法叫经验分布逆变换采样效果出奇的好。def sample_from_empirical_cdf(data, size, rngNone): if rng is None: rng np.random.default_rng() sorted_data np.sort(data) n len(sorted_data) u rng.random(size) # 注意 np.searchsorted 默认搜最左边的位置等价于 F(x) u indices np.searchsorted(np.arange(1, n 1) / n, u, sideleft) indices np.clip(indices, 0, n - 1) return sorted_data[indices]这段代码里最核心的是np.searchsorted的用法。np.arange(1, n 1) / n构造的是经验CDF在每一个样本点上的累积值searchsorted找到第一个大于等于u的位置。这里有个取舍直接返回sorted_data[indices]得到的是离散的原始样本值如果希望输出更平滑可以在两个相邻样本之间做线性插值但那样会引入人为的平滑误差具体用哪种看你对保真度的要求。这种从真实数据直接采样的方法在很多离线强化学习数据集处理里我会用到简单且稳定。4. 离散分布和混合分布的逆变换实现4.1 离散情形的标准算法离散分布的逆变换采样实现起来比连续情形更直观因为CDF是阶梯函数反函数退化成一张查找表。算法分两步计算累积概率表cumsum生成均匀随机数 $u$在累积概率表里找到第一个满足cumsum[i] u的索引代码如下import numpy as np def sample_discrete(probabilities, size, rngNone): if rng is None: rng np.random.default_rng() cumsum np.cumsum(probabilities) # 归一化防止概率和不为1导致的误差 cumsum cumsum / cumsum[-1] u rng.random(size) indices np.searchsorted(cumsum, u, sideleft) return indices这里np.searchsorted默认就是找左边界正好对应广义逆函数里 u的定义。二分类、多分类的采样本质上都是这个逻辑。如果你去看很多机器学习库的multinomial或categorical采样源码内部实现思路和这段代码是一致的。4.2 混合分布的两个采样思路混合分布是多个分布按权重组合例如 70% 来自正态分布30% 来自指数分布。对混合分布采样有两条路第一条路是层化采样先按混合权重采样一个组件编号再从对应的组件分布里采样。这是最常用的做法效率高而且组件之间的独立性天然具备。第二条路是真用逆变换采样写出整体CDF求出反函数然后套用连续或离散的逆变换。这条路只在整体CDF有解析意义时才推荐。我在实际工作里很少直接用第二条路因为它把简单问题复杂化了。只有当你的后续流程明确需要用到整体CDF的数值时才值得这么干。4.3 一个常见的索引偏移坑离散采样里最容易犯的错误是把index np.searchsorted(cumsum, u)的结果直接当成样本值但忘记了searchsorted返回的是位置而非原始类别标签。如果类别标签不是0,1,2...而是业务编码必须拿着位置去类别数组里查。还有一个边界问题如果u恰好等于0searchsorted返回0这没问题但如果浮点数精度导致cumsum最后一个值不是精确的1.0那么当u很接近1时可能返回len(probabilities)越界。所以上面代码里我先做了归一化这不仅是数学上的严谨更是防御性编程。5. 实测教训数值边界、性能优化与工程化细节5.1 浮点边界u0是最大的敌人逆变换采样公式里几乎都有ln或者1/(...)这意味着当 $u0$ 时$\ln(0)$ 会得到负无穷或者直接报错。常见的规避办法有两个生成随机数时把范围限定在 $(0,1]$u 1 - rng.random(size)因为rng.random生成的是 $[0,1)$所以1 - u就在 $(0,1]$ 上用np.nextafter(0, 1)作为下界把0这个极端值抬离实战中我倾向于第一种因为1 - rng.random(size)同时保留了指数分布推导公式的形式。但这里有个我踩过的坑如果你同时用1 - u变换又用的浮点数u恰好是 1e-17 这类极小值那么-np.log(1-u)算出来的值可能极大达到几百甚至上千这在某些仿真里会产生异常大的离群点。处理方式是不要在所有代码路径里无脑用1-u。如果你明确知道u的分布最好判断一下应用场景能否容忍极端大值。如果容忍不了可以对采样结果做截断或者用别名采样等方法绕开对数变换。5.2 性能优化向量化优先逆变换采样最大的性能优势在于它天然支持向量化。不要写这样的Python循环# 慢Python 循环逐个采样 samples [] for _ in range(100000): u rng.random() samples.append(-np.log(1 - u) / lam)而应该一次性生成一大batch# 快一次性生成十万个 u rng.random(100000) samples -np.log(1 - u) / lamNumPy 的向量化运算底层调用的是 C 实现性能差距往往是几十倍。如果你的采样分布特别复杂无法直接向量化可以考虑分批向量化或者用 Numba 的njit装饰器把采样函数编译成机器码。这两种方法我都用过在项目里把千万级样本的生成时间从秒级压到了几百毫秒以内。5.3 大批量离散采样排序法替代二分查找前面说的离散采样用np.searchsorted单次查找复杂度是 $O(\log n)$生成 $m$ 个样本就是 $O(m \log n)$。如果 $m$ 和 $n$ 都很大这个复杂度就不理想了。有一个经典优化技巧先对 $m$ 个均匀随机数排序然后利用单调性线性扫描累积概率表把总复杂度降到 $O(m \log m n)$。思路是这样的排好序的 $u$ 是递增的所以在累积概率表上的查找位置也是单调递增的。你可以用一个指针在cumsum上向右移动不需要每次都从头二分。这个技巧在处理几万类别、几百万样本的场景下非常有效我在构造大规模离散状态转移矩阵的样本时用它做过优化速度提升肉眼可见。5.4 工程化CDF单调性校验如果你要实现一个通用的逆变换采样工具库我建议在初始化时校验CDF的单调性。我见过不少因为数据预处理出错导致CDF出现非单调跳变采样结果完全错乱的情况。校验代码很简单def validate_cdf(cdf_values, tol1e-9): diffs np.diff(cdf_values) if np.any(diffs -tol): raise ValueError(CDF 必须单调不减) if not np.isclose(cdf_values[0], 0.0, atoltol) or not np.isclose(cdf_values[-1], 1.0, atoltol): raise ValueError(CDF 边界值必须接近 0 和 1)这个校验函数放在库里平时不会触发但一旦遇到脏数据它能帮你省下好几个小时的排查时间。6. 逆变换采样的边界条件和替代方案6.1 什么时候逆变换采样不是好选择逆变换采样很强但不是万能钥匙。根据我的经验它有几个明显的盲区只知道概率密度函数不知道CDF。比如某些贝叶斯模型里的后验分布只有一个未归一化的密度函数CDF的数值积分成本太高这时候逆变换采样就很尴尬。高维联合分布。逆变换采样本质上是为一维分布设计的扩展到高维时需要对每个维度做条件采样需要知道完整的条件CDF链这通常很难得到。CDF反函数计算代价太高。有些分布的CDF反函数没有解析形式每算一个样本都要跑一次数值求根生成大量样本时性能会很差。分布有截断或尖峰。比如截断正态分布虽然也可以用逆变换采样先采正态的 $u$再映射到截断区间但需要额外处理CDF在截断点的归一化一不小心就出错。6.2 我有几个备选方案遇到上面的情况我常用的替代方案是。拒绝采样只需要知道目标分布的PDF不必归一化从一个提议分布里采样按接受概率决定保留还是丢弃。实现简单但接受率低时效率很差。Box-Muller变换专门生成标准正态分布通过两个均匀随机数构造一对独立的正态样本。它比逆变换采样更快也不用求反函数。但如果要生成的不是正态分布它帮不上忙。Ziggurat算法这是个查表算法被很多标准库用来生成正态分布和指数分布速度非常快但实现复杂度高不适合你自己现写直接用标准库就行。选择建议如果目标分布的CDF反函数写起来不费劲或者你能接受数值求根选逆变换采样如果只有PDF选拒绝采样如果追求极致的性能比如每秒采样千万级直接用现成库的底层算法别自己造轮子。6.3 逆变换采样在高阶话题里的身影虽然它本身是个基础算法但很多高阶采样方法都挂着逆变换的思想。比如copula方法里要把相关结构映射到目标分布的边缘分布上用的就是边缘CDF的逆变换再比如分位数回归里的样本生成本质上也是逆变换的变体。理解了这个原理你去看那些复杂方法时会有一种“原来你也在这里”的感觉。最后分享一个我实际遇到的教训。之前做一个可靠性仿真需要生成大量服从极小形状参数的韦布尔分布样本结果发现样本里偶尔出现天文数字级的极大值把平均值拉得不成样子。排查半天才发现是 $\ln(1-u)$ 里的 $u$ 太小导致对数变换产生了极端值。后来我改用对样本做上限截断再在截断区间内做采样校准问题才解决。这种极端值不是随机性的锅而是数值边界和分布厚尾叠加产生的效应。做仿真的人尤其是模型输出会直接影响决策的场合一定不能忽略这种边界情况。