蒙特卡罗采样方法
目录
1 直接抽样法 2
1.1 离散型变量分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 连续型变量分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 多维离散型变量分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.4 多维连续型变量分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 变换抽样法 3
2.1 Box-Muller 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 球面均匀分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3 舍选抽样法 5
3.1 简单舍选法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.2 半圆分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.3 乘分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
4 重要抽样 8
1
1 直接抽样法
1.1 离散型变量分布
设变量 𝑥 是离散型的，取值为 𝑥
1
, 𝑥
2
, . . . , 相应值出现的几率为 𝑝
1
, 𝑝
2
, . . . . 那么取一个 [0, 1] 上均匀分布
的随机数 𝜉, 如果
𝑛−1
𝑖=1
𝑝
𝑖
< 𝜉 ≤
𝑛
𝑖=1
𝑝
𝑖
,
抽样结果为 𝑥
𝑛
. 这是因为概率相等
1.2 连续型变量分布
设变量 𝑥 的概率密度函数为 𝑝(𝑥), 那么即有
𝑥
−∞
𝑝(𝑥
′
)d𝑥
′
< 𝜉 ≤
𝑥+Δ𝑥
−∞
𝑝(𝑥
′
)d𝑥
′
也就是说
𝜉 =
𝑥
−∞
𝑝(𝑥)d𝑥
可以反解出 𝑥 的值
1.3 多维离散型变量分布
考虑 (𝑥, 𝑦) 的联合分布 𝑝(𝑥, 𝑦)
𝑝
1·
𝑝
11
𝑝
12
··· 𝑝
1𝑚
𝑝
2·
𝑝
21
𝑝
22
··· 𝑝
2𝑚
𝑝
3·
.
.
.
.
.
.
.
.
.
.
.
.
𝑝
4·
𝑝
𝑛1
𝑝
𝑛2
··· 𝑝
𝑛𝑚
𝑝
··
𝑝
·1
𝑝
·2
··· 𝑝
·𝑚
可以先抽出一行, 使得
𝑘−1
𝑖=1
𝑝
𝑖·
< 𝜉
1
≤
𝑘
𝑖=1
𝑝
𝑖·
即抽出了第 𝑘 行, 再抽出一列, 使得
𝑙−1
𝑖=1
𝑝
𝑘,𝑖
𝑝
𝑘·
< 𝜉
2
≤
𝑙
𝑖=1
𝑝
𝑘,𝑖
𝑝
𝑘·
那么抽出的结果为 (𝑥
𝑘
, 𝑦
𝑙
)
1.4 多维连续型变量分布
考虑 (𝑥, 𝑦) 的联合分布 𝑝(𝑥, 𝑦), 那么可以先抽出 𝑥. 可以由联合密度求出 𝑥 的边缘密度
𝑝(𝑥) =
∞
−∞
𝑝(𝑥, 𝑦)d𝑦
再由条件密度求出 𝑦
𝑝(𝑦|𝑥) =
𝑝(𝑥, 𝑦)
𝑝(𝑥)
也即
𝜉
1
=
𝑥
−∞
𝑝(𝑥
′
)d𝑥
′
=
𝑥
−∞
+∞
−∞
𝑝(𝑥
′
, 𝑦
′
)d𝑦
′
d𝑥
′
𝜉
2
=
𝑦
−∞
𝑝(𝑦|𝑥)d𝑦 =
𝑦
−∞
𝑝(𝑥, 𝑦
′
)d𝑦
′
∞
−∞
𝑝(𝑥, 𝑦
′
)d𝑦
′
若期望分布 𝑥, 𝑦 是独立的, 那么 𝑝(𝑥, 𝑦) = 𝑝(𝑥)𝑝(𝑦), 上式就可以化为
𝜉
1
=
𝑥
−∞
+∞
−∞
𝑝(𝑥
′
, 𝑦
′
)d𝑦
′
d𝑥
′
, 𝜉
2
=
+∞
−∞
𝑦
−∞
𝑝(𝑥
′
, 𝑦
′
)d𝑦
′
d𝑥
′
2 变换抽样法
变换抽样法即通过变换一个简单的分布函数 𝑔(𝑥) 通过变换得到一个复杂的分布函数 𝑝(𝑥). 最简单的情况
取 𝑔(𝑥) 为 [0, 1] 上的均匀分布
𝑔(𝑥) =
1, 0 ≤ 𝑥 ≤ 1
0, otherwise
希望通过蒙特卡罗方法得到一组服从 𝑔(𝑥) 的随机数
𝑥
1
, 𝑥
2
, ··· , 𝑥
𝑛
那么通过变换
𝑦
𝑖
= 𝑓 (𝑥
𝑖
)
使得 𝑦
𝑖
服从 𝑝(𝑦). 希望找到 𝑥, 𝑦 之间的关系 𝑓 (𝑥). 显然应该使得概率相等
𝑝(𝑦)𝑑𝑦 = 𝑔(𝑥)𝑑𝑥
通过求解该微分方程即可得到 𝑓 (𝑥). 代入均匀分布, 得到
𝑝(𝑦)𝑑𝑦 = 𝑔(𝑥)𝑑𝑥 = 𝑑𝑥
进行积分
𝑦
−∞
𝑝(𝑦)𝑑𝑦 =
𝑥
0
𝑑𝑥 = 𝑥
即得变换
𝑥 =
𝑦
−∞
𝑝(𝑦)𝑑𝑦
这说明直接抽样法是变换抽样法的特例
对于二维的情况, 假设 𝑔(𝑥, 𝑦 ) 为二维均匀分布, 期望中的分布为 𝑝(𝑥, 𝑦), 那么它们应当满足关系
𝑝(𝑥, 𝑦)𝑑𝑥𝑑𝑦 = 𝑔(𝑢, 𝑣)𝑑𝑢𝑑𝑣 = 𝑑𝑢𝑑𝑣
即
𝑝(𝑥, 𝑦) =
𝜕(𝑢, 𝑣)
𝜕
(
𝑥, 𝑦
)
2.1 Box-Muller 法
考虑标准正态分布的分布函数
𝑝(𝑥) =
1
√
2𝜋
𝑒
−
𝑥
2
2
直接积分是无法得到解析形式的
𝑦
−∞
𝑝(𝑦)𝑑𝑦 = 𝑥
考虑变换
𝑥 = −
√
−2 ln 𝑢 cos(2𝜋𝑣), 𝑦 = −
√
−2 ln 𝑢 sin(2𝜋𝑣)
可以得到其反变换
𝑢 = 𝑒
−
𝑥
2
+ 𝑦
2
2
, 𝑣 =
1
2𝜋
arctan
𝑦
𝑥
求其雅可比行列式
𝜕(𝑢, 𝑣)
𝜕(𝑥, 𝑦)
=
1
2𝜋
𝑒
−
𝑥
2
+ 𝑦
2
2
那么 𝑥, 𝑦 的联合分布函数 𝑝(𝑥, 𝑦) 与 𝑢, 𝑣 的联合分布函数 𝑔(𝑢, 𝑣) 满足
𝑝(𝑥, 𝑦) =
𝜕(𝑢, 𝑣)
𝜕(𝑥, 𝑦)
𝑔(𝑢, 𝑣)
这意味着如果 𝑢, 𝑣 服从均匀分布, 那么 𝑥, 𝑦 就服从正态分布
𝑝(𝑥, 𝑦) =
1
√
2𝜋
𝑒
−
𝑥
2
+ 𝑦
2
2
它可以分离为两个独立的正态分布
𝑝(𝑥, 𝑦) = 𝑝(𝑥)𝑝(𝑦), 𝑝(𝑥) =
1
√
2𝜋
𝑒
−
𝑥
2
2
, 𝑝(𝑦) =
1
√
2𝜋
𝑒
−
𝑦
2
2
2.2 球面均匀分布
二维球 (即圆面) 上的均匀分布固然可以通过 𝜙 在 [0, 2𝜋],𝑟 在 [0, 1] 上的均匀分布取
𝑥 = 𝑟 cos 𝜙, 𝑦 = 𝑟 sin 𝜙
得到, 但是三角函数需要大量的计算, 可以取一对均匀分布的随机数 (𝑢, 𝑣)
𝑝(𝑢) =
1, 0 ≤ 𝑢 ≤ 1
0, otherwise
𝑝(𝑣) =
1, 0 ≤ 𝑣 ≤ 1
0, otherwise
舍去其中 𝑢
2
+ 𝑣
2
> 1 的点, 就得到了圆面上的均匀分布
𝑝(𝑥, 𝑦) =
1
𝜋
, 𝑥
2
+ 𝑦
2
≤ 1
0, otherwise
如果要求在二维球面 (即圆) 上的均匀分布, 则可以取
(𝑥
′
, 𝑦
′
) =
(𝑥, 𝑦)
𝑥
2
+ 𝑦
2
可以进一步得到一般地高维球和高维球面上的均匀分布. 设要取 𝑛 维球上的均匀分布, 那么可以取 𝑛 个均
匀分布的随机数
𝑢
1
, 𝑢
2
, ··· , 𝑢
𝑛
, 𝑝 (𝑢
𝑖
) =
1, 0 ≤ 𝑢
𝑖
≤ 1
0, otherwise
舍去其中 𝑢
2
1
+ 𝑢
2
2
+ ··· + 𝑢
2
𝑛
> 1 的点, 就得到了 𝑛 维球上的均匀分布
𝑝(𝑥
1
, 𝑥
2
, ··· , 𝑥
𝑛
) =
Γ(
𝑛
2
)
𝜋
𝑛
2
, 𝑥
2
1
+ 𝑥
2
2
+ ··· + 𝑥
2
𝑛
≤ 1
0, otherwise
其中 Γ 为伽玛函数, 用于表示 𝑛 维球的体积. 同理, 只需要除以半径即可得到 𝑛 维球面上的均匀分布
(𝑥
′
1
, 𝑥
′
2
, ··· , 𝑥
′
𝑛
) =
(𝑥
1
, 𝑥
2
, ··· , 𝑥
𝑛
)
𝑥
2
1
+ 𝑥
2
2
+ ··· + 𝑥
2
𝑛
, 𝑝 (𝑥
′
1
, 𝑥
′
2
, ··· , 𝑥
′
𝑛
) =
Γ
(
𝑛
2
)
𝜋
𝑛
2
−1
, 𝑥
2
1
+ 𝑥
2
2
+ ··· + 𝑥
2
𝑛
= 1
0, otherwise
3 舍选抽样法
舍选法抽样即利用一个二维分布 𝑔(𝑥, 𝑦) 通过 𝑦 的抽样舍去部分 𝑥, 使得 𝑥 的分布为 𝑝(𝑥). 即希望寻找一
个
ℎ
(
𝑥
)
,
遵循以下的抽样方法
1. 从 𝑔(𝑥, 𝑦) 中抽取一个 (𝑥, 𝑦)
2. 如果 𝑦 < ℎ(𝑥), 则保留 𝑥, 否则舍去
3. 重复以上步骤
希望抽出来的 𝑥 服从分布 𝑝(𝑥). 𝑥 的密度为条件密度
𝑝(𝑥) = 𝑝
𝜉
𝑥
|𝜉
𝑦
< ℎ(𝜉
𝑥
)
意为, 𝑥 的概率密度为 𝜉
𝑦
小于 ℎ(𝜉
𝑥
) 的条件下 𝜉
𝑥
的条件概率密度, 这是 𝜉
𝑥
的条件分布. 为了求出它, 考
虑分布函数
𝐹
𝜉
𝑥
|𝜉
𝑦
<ℎ (𝜉
𝑥
)
(𝑥) =
𝑝
𝜉
𝑥
< 𝑥, 𝜉
𝑦
< ℎ(𝜉
𝑥
)
𝑝(𝜉
𝑦
< ℎ(𝜉
𝑥
))
分子分母可以分别计算
𝑝
𝜉
𝑥
< 𝑥, 𝜉
𝑦
< ℎ(𝜉
𝑥
)
=
𝑥
−∞
d𝜉
𝑥
ℎ(𝜉
𝑥
)
−∞
𝑔(𝜉
𝑥
, 𝜉
𝑦
)d𝜉
𝑦
𝑝(𝜉
𝑦
< ℎ(𝜉
𝑥
)) =
∞
−∞
d𝜉
𝑥
ℎ(𝜉
𝑥
)
−∞
𝑔(𝜉
𝑥
, 𝜉
𝑦
)d𝜉
𝑦
进而得到条件分布函数为
𝐹
𝜉
𝑥
|𝜉
𝑦
<ℎ (𝜉
𝑥
)
(𝑥) =
𝑥
−∞
d𝜉
𝑥
ℎ(𝜉
𝑥
)
−∞
𝑔(𝜉
𝑥
, 𝜉
𝑦
)d𝜉
𝑦
∞
−∞
d𝜉
𝑥
ℎ(𝜉
𝑥
)
−∞
𝑔(𝜉
𝑥
, 𝜉
𝑦
)d𝜉
𝑦
对其求导即可得到条件密度
𝑝(𝑥) =
d
d𝑥
𝐹
𝜉
𝑥
|𝜉
𝑦
<ℎ (𝜉
𝑥
)
(𝑥) =
ℎ(𝑥)
−∞
𝑔(𝑥, 𝜉
𝑦
)d𝜉
𝑦
∞
−∞
d𝜉
𝑥
ℎ(𝜉
𝑥
)
−∞
𝑔(𝜉
𝑥
, 𝜉
𝑦
)d𝜉
𝑦
3.1 简单舍选法
考察简单的情况, 期望分布 𝑝(𝑥) 定义在有限区间 [𝑎, 𝑏] 上, 则可取 𝑔(𝑥, 𝑦) 为 [𝑎, 𝑏] × [0, 1] 上的均匀分
布, 那么
𝑝(𝑥) =
ℎ(𝑥)
0
d𝑦
𝑏
𝑎
d𝜉
𝑥
ℎ(𝜉
𝑥
)
0
d𝜉
𝑦
=
ℎ(𝑥)
𝑏
𝑎
ℎ(𝜉
𝑥
)d𝜉
𝑥
那不妨就取 ℎ(𝑥) = 𝑝(𝑥) 即可满足条件. 考虑它的几何意义如下图
它实际上是取了一个 [𝑎, 𝑏] × [0, 1] 的矩形, 并在矩形中均匀取点. 若是落在 𝑝(𝑥) 的下方, 则保留 𝑥, 否则
舍去. 由于取点是均匀的, 那么保留的点也是均匀的, 因而可以使得采样的分布为 𝑝(𝑥)
不过当分布是一个” 尖峰” 的函数时, 该方法抽样的效率很低. 这时可以考虑利用变换抽样法抽取一个简
单的分布函数 𝑘 (𝑥), 然后进一步通过舍选法抽取 𝑝(𝑥)
抽样效率低的问题在于矩形区域内的点很多都被舍去了, 因而可以考虑使用一个更接近的函数, 减少无效
的点, 提高抽样效率. 一个比较好的选择是分段阶梯函数
3.2 半圆分布
希望抽样得到半圆上的分布
𝑝(𝑥) =
1
𝜋
√
1 − 𝑥
2
取两个均匀分布的随机变量
𝜉 ∼ 𝑈(0, 1), 𝜂 ∼ 𝑈(−1, 1)
那么它们的联合分布为
𝑔(𝜉, 𝜂) =
1
2
, 0 ≤ 𝜉 ≤ 1, −1 ≤ 𝜂 ≤ 1
作变换
𝑥 =
𝜉
2
− 𝜂
2
𝜉
2
+ 𝜂
2
, 𝑦 = 𝜉
2
+ 𝜂
2
则可求得 (𝑥, 𝑦) 的联合分布
𝑝(𝑥, 𝑦) = 𝑔(𝜉, 𝜂)
𝜕(𝜉, 𝜂)
𝜕(𝑥, 𝑦)
=
1
4𝜋
√
1 − 𝑥
2
那么自然
𝑝(𝑥) =
4
𝜋
1
0
𝑝(𝑥, 𝑦)d𝑦 =
1
𝜋
√
1 − 𝑥
2
正是所求的半圆分布. 注意到本节舍选法中有
𝑝(𝑥) =
ℎ(𝑥)
−∞
𝑔(𝑥, 𝜉
𝑦
)d𝜉
𝑦
∞
−∞
d𝜉
𝑥
ℎ(𝜉
𝑥
)
−∞
𝑔(𝜉
𝑥
, 𝜉
𝑦
)d𝜉
𝑦
考虑到 𝑝(𝑥, 𝑦) 的定义域, 那么希望
𝑝(𝑥) =
ℎ(𝑥)
0
𝑝(𝑥, 𝑦)d𝑦
1
−1
d𝜉
𝑥
ℎ(𝜉
𝑥
)
0
𝑝(𝜉
𝑥
, 𝜉
𝑦
)d𝜉
𝑦
注意到取 ℎ(𝑥) = 1 即可使得等式成立, 因而得到舍选法抽样
1. 从 [0, 1] × [−1, 1] 上均匀分布中抽取一个 (𝜉, 𝜂)
2. 由 𝑥 =
𝜉
2
− 𝜂
2
𝜉
2
+ 𝜂
2
, 𝑦 = 𝜉
2
+ 𝜂
2
得到 (𝑥, 𝑦)
3. 如果 𝑦 > 1, 则舍去 (𝑥, 𝑦), 否则保留 𝑥
由此得到了 cos 𝜑, 𝜑 ∈ [0, 𝜋] 的均匀分布抽样. 同样还可以取
sin 𝜑 =
2𝜂
𝜉
2
+ 𝜂
2
得到 sin 𝜑, 𝜑 ∈ [0, 𝜋] 的均匀分布抽样
3.3 乘分布
乘分布即密度函数可以拆成乘积的形式
𝑝(𝑥) = ℎ(𝑥)𝑞(𝑥)
其中 𝑞(𝑥) 是合法的密度函数, 即满足归一化
∞
−∞
𝑞(𝑥)d𝑥 = 1
那么可以先产生一个服从 𝑞(𝑥) 的随机数 𝑥, 然后产生一个 [0, ℎ
𝑚𝑎 𝑥
] 上的均匀分布的随机数 𝑦, 如果 𝑦 <
ℎ(𝑥), 则保留 𝑥, 否则舍去. 这样得到的 𝑥 服从 𝑝 (𝑥)
4 重要抽样
重要抽样希望通过对贡献大的位置 (即分布多的地方) 多抽样, 从而减小计算误差. 设期望抽样的密度函
数为 𝑓 (𝑥), 有一个与它相近的较简单的密度函数 𝑔( 𝑥) 满足
𝑓 (𝑥)
𝑔(𝑥)
≈
1
,
∞
−∞
𝑔
(
𝑥
)
d
𝑥
=
1
那么 𝑓 (𝑥) 的积分就变为了
𝑏
𝑎
𝑓 (𝑥)d𝑥 =
𝑏
𝑎
𝑓 (𝑥)
𝑔(𝑥)
𝑔(𝑥)d𝑥
此时的 𝑔(𝑥) 就可以视为权值.抽样的权值可以表现在抽样点的分布密度上. 若是以 𝑔( 𝑥) 为抽样的概率密
度函数, 则积分的求和可以化为
𝑏
𝑎
𝑓 (𝑥)d𝑥 ≈
𝑓 (𝑥)
𝑔(𝑥)
≈
1
𝑁
𝑁
𝑖=1
𝑓 (𝑥
𝑖
)
𝑔(𝑥
𝑖
)
