一种通用随机数生成方式 -- 逆变换采样

2021-08-09 16:57 作者:nyasyamorina 0人读过 | 我要投稿

前言

这里讨论的随机数(random number)为随机实数, 对于随机整数(random integer)会叫做...随机整数, 嗯.

实际上自己并没有学过概率论什么的, 自己本质上就是一个小会高数的高中生罢了. 以下内容赌的成分偏多 (

在实际应用中, 生成任意随机数都可以分为三步: 1) 生成随机整数. 2) 从随机整数得出均匀分布随机数. 3) 从均匀分布随机数里得出其他分布随机数.

随机整数

生成随机整数的方法有很多, 有简单的线性同余发生器(Linear congruential generator, LCG), 有复杂的梅森旋转器(Mersenne twister), 甚至可以使用现实电路产生真随机数. 这里并不会详细介绍产生随机整数的方法, 有兴趣的可以去翻翻wiki(见底部).

大部分随机整数生成器会指定整数分布的区间 [MIN, MAX], 给出的随机整数将会均匀分布在区间里. 常用的随机整数生成器区间为 $%5B0%2C%202%5E%7B32%7D-1%5D$ , 或是 $%5B0%2C%202%5E%7B64%7D-1%5D$ , 刚好是uint32和uint64可以表示的全部范围. 确保随机整数分布的均匀性是非常重要的, 因为这是为了方便后续步骤变换随机数的分布方式.

均匀分布随机数

常用的均匀分布随机数的分布区间为 [0, 1), 需要注意的是半开区间, 也就是说实数1是不能取到的.

把均匀分布的随机整数转化为均匀分布的随机数是非常简单的, 把随机数减去MIN再除以MAX-MIN+1就可以得到分布在[0,1)里的随机数了. 因为一开始生成的随机整数不是连续分布的, 所以使用这种方法转化的随机数也不是连续分布的, 但是对于实际使用来说已经足够好了.

了解了如何生成均匀分布的随机数, 接下来就是生成其他分布形式的随机数了. 因为这里是说通用的生成方法, 所以有必要了解什么是随机数的分布.

概率密度函数

概率密度函数(probability density function, pdf)是一个用来描述随机数分布在给定区间的可能性的函数. 对指定区间下的pdf进行积分即为随机数分布在这个区间里的概率. 即

$Prob(a%5Cleq%20x%3Cb)%3D%5Cint_a%5Eb%5Crho(x)dx$

例子: 给出一个随机数生成方式: 2 * random(), 对应的pdf为

$%5Crho(x)%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D0.5%2Cx%5Cin%5B0%2C2)%5C%5C0%2Cotherwise%5Cend%7Bmatrix%7D%5Cright.$

则生成随机数会以50%概率分布在0到1里, 以100%概率分布在-2到2里, 以此类推.

另外一个著名的例子就是正态分布, 下面给出标准正态分布的pdf:

$%5Crho(x)%3D%5Cfrac%7B1%7D%7B%5Csqrt%7B2%5Cpi%7D%7Dexp(-%5Cfrac%7Bx%5E2%7D%7B2%7D)$

符合标准正态分布的随机数会以68.27%概率分布在-1到1之间, 95.45%分布在-2到2之间, 并且以50%分布在0到+∞之间.

累积分布函数

累积分布函数(Cumulative Distribution Function, cdf)定义为随机分布小于给定数字的概率, 根据pdf的定义, cdf为

$F(x)%3D%5Cint_%7B-%5Cinfty%7D%5Ex%5Crho(t)dt$

根据cdf, 随机数分布在特定区间下的概率可以表示为

$Prob(a%5Cleq%20x%3Cb)%3DF(b)-F(a)$

例子: 以2*random()生成的随机数的cdf为

$F(x)%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D0%2Cx%3C0%5C%5C0.5x%2Cx%5Cin%5B0%2C2)%5C%5C1%2Cx%5Cgeq2%5Cend%7Bmatrix%7D%5Cright.$

而标准正态分布的cdf为

$F(x)%3D%5Cfrac%7B1%7D%7B%5Csqrt%7B2%5Cpi%7D%7D%5Cint_%7B-%5Cinfty%7D%5Exexp(-%5Cfrac%7Bt%5E2%7D%7B2%7D)dt$

逆累积分布函数

因为pdf为随机数的概率密度函数, 所以pdf会符合关于"概率"的规则: pdf必定大于等于0 (不存在负数概率), pdf在全体实数上积分为1(概率总和为1). 因为cdf是pdf的积分, pdf的第一条特性定义了cdf必定是单调递增的, pdf的第二条特性定义了cdf有下界0和上界1.

因为cdf是单调递增的, 所以cdf存在逆函数, 并且cdf的逆也是单调递增的. cdf有下界0和上界1说明cdf的逆定义域为[0, 1]

从均匀随机数里生成任意分布的随机数

暴论: 已知某种特定分布的逆累积分布函数inv-cdf, 则把均匀分布的随机数经过inv-cdf处理后得到的随机数符合这种分布, 即

$Y%3DF%5E%7B-1%7D(X)$

其中X是均匀分布随机数, 则Y符合特定分布.

这是十分容易证明的: $Prob(F%5E%7B-1%7D(X)%5Cleq%20Y)%3DProb(X%5Cleq%20F(Y))$ , 因为X是均匀分布的, 所以有 $Prob(X%5Cleq%20F(Y))%3DF(Y)$ , 如此可知Y满足特定分布.

例子: 三角形分布

三角形分布因为它的pdf形如三角形而得名. 三角形分布需要给出下限a, 上限b, 和众数c.

三角形分布的pdf定义为

$%5Crho(x)%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D%5Cfrac%7B2(x-a)%7D%7B(b-a)(c-a)%7D%2Cx%5Cin%5Ba%2Cc)%5C%5C%5Cfrac%7B2(b-x)%7D%7B(b-a)(b-c)%7D%2Cx%5Cin%5Bc%2Cb)%5C%5C0%2Cotherwise%5Cend%7Bmatrix%7D%5Cright.$

对pdf积分求逆得到逆累积分布函数:

$F%5E%7B-1%7D(x)%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Da%2B%5Csqrt%7B(b-a)(c-a)x%7D%2Cx%3C%5Cfrac%7Bc-a%7D%7Bb-a%7D%5C%5Cb-%5Csqrt%7B(b-a)(b-c)(1-x)%7D%2Cotherwise%5Cend%7Bmatrix%7D%5Cright.$

如果X是均匀分布的随机数, 那么 $F%5E%7B-1%7D(X)$ 是三角形分布, 可以简单地写一个jio本进行验证

例子2: 正态分布

正态分布的pdf为

$%5Crho(x%3B%5Csigma%2C%5Cmu)%3D%5Cfrac%7B1%7D%7B%5Csigma%5Csqrt%7B2%5Cpi%7D%7D%5Cexp(-%5Cfrac%7B(x-%5Cmu)%5E2%7D%7B2%5Csigma%5E2%7D)$

那么cdf为

$F(x%3B%5Csigma%2C%5Cmu)%3D%5Cfrac%7B1%7D%7B2%7D%5Cleft(1%2Berf(%5Cfrac%7Bx-%5Cmu%7D%7B%5Csigma%5Csqrt2%7D)%5Cright)$

其中 $erf(x)%3D%5Cfrac%7B2%7D%7B%5Csqrt%5Cpi%7D%5Cint_0%5Ex%5Cexp(-t%5E2)dt$ 称为误差函数(error function), 则cdf的逆为

$F%5E%7B-1%7D(x%3B%5Csigma%2C%5Cmu)%3D%5Csigma%5Csqrt2erf%5E%7B-1%7D(2x-1)%2B%5Cmu$

需要注意到的是, erf并不是初等函数, erf的逆也不能直接求出. 但幸运的是, 目前erf的逆有许多近似的方法. 使用一种快速但低准确度的逆erf进行验证随机数的正确性:

关于逆变换采样的讨论

1) 逆变换采样具有通用性, 所有随机分布都有相应的分布函数pdf, 由pdf推导出cdf的逆即可生成任意分布的随机数. 但是需要手动计算积分和求逆, 很有可能不存在初等函数可以表示的解. 2) 大部分随机过程都是由随机整数生成器为基础, 而随机整数生成器存在周期性, 比如说LCG只有 $2%5E%7B32%7D$ 长度的周期, 而梅森旋转器的周期长达 $2%5E%7B19937%7D-1$ . 使用其他产生特定分布的方法可能会使得周期产生变化, 比如使用LCG与C++TR1里的方法产生正态分布会使得周期缩短至 $2%5E%7B16%7D$ , 这是因为在C++TR1里的生成正态分布的方法需要取两次随机数才可以生成一个正态分布随机数. 而对于逆变换采样, 不管什么分布都只需要取一次随机数. 3) 逆变换采样具有对比拒绝采样(rejection sampling)更快的运算速度. 对于拒绝采样以及其衍生算法, 需要采用大量的判断语句, 而分支对运行速度是不好的 (特别是GPU运算). 因为逆变换采样需要用到大量初等函数, 所以面对其他的优化算法就显得特别累赘了.

示例代码在下面贴出, 运行需要Python3 + numpy + matplotlib.

END.

吐槽

为什么会有一篇这样的东西出来呢. 其实是最近自己在用C艹写光追, 打算自己造轮子的时候遇上了"如何生成正态分布随机数"这个难题(baidu-less). 然后灵机一动自己发现了这种计算方法, 但是直到专栏写到一半才看到了原来这种方法叫逆变换采样. 嗨, 还以为发现了什么新大陆.

关于上面的示例代码, 原本是打算用C艹写的, 顺便做一波性能分析, 但是突然发现C艹画图实在太麻烦的(相对于matplotlib), 于是还是滚回了心爱的python下. 原本代码是打算发到gayhub上的, 但是感觉为了一个100多行的文件去开一个仓库也太小题大做了, 而且发百度云有点low, 也不知道有没有其他分享代码的方法, 就直接在文章结尾贴出了.

更多资讯:

梅森旋转器: en.wikipedia.org/wiki/Mersenne_Twister
线性同余发生器: en.wikipedia.org/wiki/Linear_congruential_generator
概率密度函数: en.wikipedia.org/wiki/Probability_density_function
累积分布函数: en.wikipedia.org/wiki/Cumulative_distribution_function
逆变换采样: en.wikipedia.org/wiki/Inverse_transform_sampling
如何产生正态分布的随机数？zhihu.com/question/29971598
逆误差函数(1) Anthony J. Strecok - On the Calculation of the Inverse of the Error Function
逆误差函数(2) Mathematics of Computation - Rational Chebyshev Approximations for the Inverse of the Error Function
逆误差函数(3) Mike Giles - Approximating the erfinv function

标签：