蒙特卡罗方法
目录
1 随机数生成器 2
1.1 随机数生成器的要求 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Lehmer 线性同余法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 16807 产生器 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Tausworthe 位移计数器法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.1 R250 产生器 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Fibonacci 延迟产生器 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4.1 带载减法产生器 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4.2 带载减法 Weyl 产生器 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.5 Marsaglia 产生器 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.6 伪随机数的统计检验 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.6.1 独立性与相关系数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.6.2 均匀性检验－频率检验 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.6.3 多维频率检验 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1
1 随机数生成器
1.1 随机数生成器的要求
1. 周期长: 伪随机数一定具有周期, 希望周期尽可能长, 如大于 int 的最大值 2147483647
2. 随机性好: 生成的随机数列之间没有明显的相关性
3. 生成速度快: 生成随机数的速度要快
1.2 Lehmer 线性同余法
线性同余法的迭代公式为
𝐼
𝑛+1
= (𝑎𝐼
𝑛
+ 𝑏) mod 𝑚
为了得到 [0, 1] 之间的随机数, 可以将 𝐼
𝑛
除以 𝑚, 即
𝑥
𝑛
=
𝐼
𝑛
𝑚
其中 𝑚 的取值通常为 𝑖𝑛𝑡32 的最大值 2147483647
1.2.1 16807 产生器
16807 产生器是线性同余法的一种特例, 它已经通过了许多理论测试和大量实用的考验, 被认为是 32 位
机上最有效的产生器. 其参数为
𝑎 = 16807, 𝑏 = 0, 𝑚 = 2147483647
可惜当 𝐼
𝑛
很大时, 计算 𝑎𝐼
𝑛
会溢出. 因而需要 Schrage 方法进行取模. 考虑将 𝑚 分解
𝑚 = 𝑎 · 𝑞 + 𝑟, 0 ≤ 𝑟 < 𝑎
它们都是整数, 那么
𝑞 =
𝑚
𝑎
, 𝑟 = 𝑚 mo d 𝑎
考虑上一个迭代的结果为 𝑧, 它的上界为 𝑚, 那么考虑将下一次迭代变形
𝑎𝑧 =
𝑧
𝑞
(𝑎𝑞 + 𝑟) +
𝑟𝑧
𝑞
=
𝑧
𝑞
(𝑎𝑞 + 𝑟) + 𝑎(𝑧 mod 𝑞)𝑞 −
𝑧
𝑞
𝑟
由于 𝑎𝑞 + 𝑟 = 𝑚, 所以取模可以将其舍去, 那么
𝑎𝑧 ≡ 𝑎(𝑧 mod 𝑞) −
𝑧
𝑞
𝑟 mod 𝑚
最后根据符号判断是否应该加上 𝑚
𝐼
𝑛+1
=
𝑎(𝑧 mod 𝑞) −
𝑧
𝑞
𝑟, if 𝑎(𝑧 mod 𝑞) −
𝑧
𝑞
𝑟 ≥ 0
𝑎(𝑧 mod 𝑞) −
𝑧
𝑞
𝑟 + 𝑚, otherwise
1.3 Tausworthe 位移计数器法
Tausworthe 方法又称一般反馈位移计数器法, 它用于消除线性同余法中的关联. 假定已经有了一个随机
数列 {𝐼
1
, 𝐼
2
, · · · , 𝐼
𝑛
}, 那么就可以
𝐼
𝑛
= 𝐼
𝑛− 𝑝
⊕ 𝐼
𝑛−𝑞
其中 ⊕ 表示异或运算, 而 𝑝, 𝑞 是两个整数, 它们的选择是关键的. 一般来说, 它们的选择应该使得
𝑝
2
+ 𝑞
2
+ 1 𝑖𝑠 𝑝𝑟𝑖𝑚𝑒
这被称为梅森素数, 常见的可以列举如下
[31, 3], [98, 27], [250, 103], [1279, 216418], [9689, 84471183624444187]
一般来说,𝑝, 𝑞 越大产生的随机数质量越好. 通常使用线性同余法产生具有 𝑝 个数的初始序列, 然后再使
用 Tausworthe 方法产生随机数列
1.3.1 R250 产生器
R250 产生器是一种 Tausworthe 方法的特例, 它的参数为
𝑝 = 103, 𝑞 = 250
1.4 Fibonacci 延迟产生器
位移计数器是一种特殊的 Fibonacci 延迟产生器. 一般的 Fibonacci 延迟产生器的迭代公式为
𝐼
𝑛
= 𝐼
𝑛− 𝑝
⊗ 𝐼
𝑛−𝑞
mod 𝑚
其中 ⊗ 可以表示加减乘除异或等运算, 整数对 [𝑝, 𝑞] 表示延迟.Fibonacci 延迟产生器的周期非常长, 在
32
位机上可以达到
(
2
𝑝
− 1
)
2
31
1.4.1 带载减法产生器
带载减法产生器是一种复杂的类 Fibonacci 延迟产生器, 它的迭代公式为
𝐼
𝑛
= 𝐼
𝑛
−
22
− 𝐼
𝑛
−
43
− 𝐶
其中 𝐶 由以下公式确定
𝐶 =
0, if 𝐼
𝑛
≥ 0
1, 𝐼
𝑛
= 𝐼
𝑛
+
2
32
− 5
otherwise
1.4.2 带载减法 Weyl 产生器
带载减法 Weyl 产生器是一种复杂的类 Fibonacci 延迟产生器. 它先通过带载减法产生器产生一个序列 𝐽
𝑛
𝐽
𝑛
= 𝐽
𝑛−22
− 𝐽
𝑛−43
− 𝐶, 𝐶 =
0, if 𝐽
𝑛
≥ 0
1, 𝐽
𝑛
= 𝐽
𝑛
+
2
32
− 5
otherwise
然后再由 𝐽
𝑛
产生 𝐼
𝑛
𝐾
𝑛
=
(
𝐾
𝑛−1
− 362436069
)
mod 2
32
, 𝐼
𝑛
= (𝐽
𝑛
− 𝐾
𝑛
) mod 2
32
1.5 Marsaglia 产生器
Marsaglia 产生器使用两种随机数生成器生成随机数列. 首先使用减法操作的 Fibonacci 延迟产生器
𝑥
𝑛
=
𝑥
𝑛−97
− 𝑥
𝑛−33
, if 𝑥
𝑛−97
− 𝑥
𝑛−33
≥ 0
𝑥
𝑛−97
− 𝑥
𝑛−33
+ 1, otherwise
第二个随机数生成器是
𝑎 ◦ 𝑏 =
𝑎 − 𝑏, if ≥ 0
𝑎 − 𝑏 + 16777213/16777216, otherwise
迭代公式为
𝑦
𝑛
= 𝑦
𝑛−1
◦ (16777213/16777216)
最后组合两个随机数列
𝑧
𝑛
= 𝑥
𝑛
◦ 𝑦
𝑛
1.6 伪随机数的统计检验
1.6.1 独立性与相关系数
为了验证随机数序列中 𝑥
𝑛
与 𝑥
𝑛+𝑙
之间是否具有线性关系, 定义相关系数
𝐶
𝑙
=
⟨
𝑥
𝑛
𝑥
𝑛+𝑙
⟩
−
⟨
𝑥
𝑛
⟩
𝑥
2
𝑛
−
⟨
𝑥
𝑛
⟩
2
其中
⟨
𝑥
⟩
意为 𝑥 的平均值. 若 𝑥
𝑛
与 𝑥
𝑛+𝑙
毫不相关, 那么期望就应当有
⟨
𝑥
𝑛
𝑥
𝑛+𝑙
⟩
=
⟨
𝑥
𝑛
⟩ ⟨
𝑥
𝑛+𝑙
⟩
由于它们只是同一个序列中的不同位置, 显然都应等于
⟨
𝑥
𝑛
⟩
, 进而
𝐶
𝑙
= 0
因而期望中 𝐶
𝑙
应当越小越好
1.6.2 均匀性检验－频率检验
期望中理想的随机数生成器的序列应当是均匀的, 即它们都应服从均匀分布. 将随机数的区间划分为 𝐾
个小区间, 每个小区间的频数 𝑛
𝑘
应当都区域理论值 𝑚
𝑘
𝑚
𝑘
=
𝑁
𝐾
采用
拟合优度检验
确定随机数序列是否服从均匀分布. 取原假设
𝐻
0
: 随机数序列服从均匀分布
那么取统计量
𝑇 =
𝐾
𝑘=1
(𝑛
𝑘
− 𝑚
𝑘
)
2
𝑚
𝑘
∼ 𝜒
2
𝑘−1
它就应该服从卡方分布.𝑇 = 0 时显然支持原假设, 因而取 𝛼 上分位数, 拒绝域为
𝑇 > 𝜒
2
𝑘−1
(1 − 𝛼)
1.6.3 多维频率检验
对于多维的情形, 不妨以二维为例. 将二维平面划分为 𝐾 × 𝐾 个小区间, 那么每个小区间的频数 𝑛
𝑖 𝑗
应当
都接近于理论值 𝑚
𝑖 𝑗
𝑚
𝑖 𝑗
=
𝑁
𝐾
2
取统计量
𝑇 =
𝐾
𝑖=1
𝐾
𝑗=1
(𝑛
𝑖 𝑗
− 𝑚
𝑖 𝑗
)
2
𝑚
𝑖 𝑗
∼ 𝜒
2
(𝐾 −1)
2
同样可以构造卡方变量进行检验
