本文为第一系列的原创文章之一, 目的在于总结“Simulation” -Sheldon M.Ross一书中出现的模拟随机变量的方法.
我们平时所接触的各种随机变量和这些随机变量的分布, 往往会有一个确定的表达式. 但是很多复合的随机变量, 其分布往往很难被确定. 即使是确定的表达式, 当我们需要生成这个随机变量时, 也很难直接从其表达式得到.
所以我们需要对这些随机变量的生成过程进行模拟, 根据模拟得到的变量数据, 可以进行进一步的模拟, 或者得到其均值与方差等我们需要研究的数据.
一般而言, 所有的随机变量都可以通过最简单的分布: 均匀分布得到. 接下来, 我们将从均匀分布出发, 得到各种随机变量的模拟方法.
随机数与均匀分布
要想得到均匀分布, 首先需要得到随机数. 伪随机数的生成有很多种, 比较常用的有取模法. 根据迭代公式, 第 n 个数 \( x_n \)的产生为:
\[x_n = ax_{n-1}\quad modulo\quad m\]
其中 \(modulo\) 表示取模运算. 得到的伪随机数 \(x_n\) 的取值在 \(0\) 到 \(m-1\) 之间, 将 \(x_n\) 除以 \(m\) , 当 \(m\) 足够大, 可以认为 \({x_n}/{m}\)服从 \(Uniform(0,1)\) .
书中讲到: 对于数为 32 位有符号数的机器, 我们一般取 \( m=2^{31}-1, a=7^5=16807 \) , 这是不符合逻辑的, 32 位有符号数的最大正数本身即为 \( m=2^{31}-1 \) , 对一个机器里能表示的最大的数取模相当于其本身.
经过读者考证发现, 这种 \(m\) 与 \(a\) 的取法来源于Lehmer Random Number Generator, 也被称为Park–Miller Random Number Generator (由Stephen K. Park 与 Keith W. Miller 提出), 属于比较传统的线性同余随机数参数取法. 在 Apple CarbonLib 这一 C 语言 API 中被使用. 在C99标准中, 执行的不是 32 位数, 而是 64 位数, 源码如下:
uint32_t lcg_parkmiller(uint32_t a) {
return ((uint64_t)a * 48271UL) % 2147483647UL;
}作者可能并非是计算机专业, 从而会出现这种错误而没有注意. (并不是误写, 作者还写道, 对于 36 位机器, 取值应当是 \(m=2^{35-1}\) , 且不说 36 位机器多么不常见, 这里依旧是不对的.)
总的来说, 这种随机数的方法被称为 Linear Congruential Generators, 即 LCG 方法, 有时候 我们用 \((ax_{n-1} + c)\) 取代 \(ax_{n-1}\) , 这被称为 Mixed Congruential Generators.
无论是哪种方法, LCG 不可避免的有以下几个短板:
○ 周期短, 无论怎么取值, 周期无法超过 \(m\) , 尤其是当 \(m\) 为 \(2\) 的整数次幂时, 低位的周期会更短.
○ 伪随机不是真正的随机, 如果LCG被用作在 \(n\) 维空间内随机取点, 这些点都会落到最多 \((n!m)^{\frac{1}{n}}\) 个超平面上.
所以 LCG作为一种简单的随机数算法, 当对于模拟量不大的计算, 可以考虑. 要寻求比较大的周期的随机数产生方法时, 可以考虑使用 Mersene Twister , 周期可以大到 \(2^{19937} – 1\) (这是一个梅森素数), 而且时间上并没有很大的增加, 唯一的缺点是需要占用 2.5 KB的计算空间, 这个算法的 mini 版也需要占用 127 位.
在得到随机数之后, 书中给了一些利用随机数进行简单模拟的例子, 例如计算积分, 计算 \(\pi\) , 产生随机的排列, 由于数学原理过于简单, 这里不再陈述.
接下来我们将看到如何产生泊松分布与二项分布的随机序列.