【Mento Carlo 3】给定分布生成随机数

3 years ago

source link: https://www.guofei.site/2017/08/21/randomgenerator1.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.

2017年08月21日

Author: Guofei

文章归类: A蒙特卡洛方法，文章编号: 10003

版权声明：本文作者是郭飞。转载随意，但需要标明原文链接，并通知本人
原文链接：https://www.guofei.site/2017/08/21/randomgenerator1.html

Edit

另一篇博客随机数发生器，讲了如何利用迭代法生成随机的U(0,1)U(0,1)均匀分布。
那么，如何利用均匀分布来生成特定的分布呢？

CDF的逆

目的：生成随机变量X对应的随机数。
已知：
R∼U(0,1)R∼U(0,1) 随机变量X的CDF是F(x)F(x),

方法：F(X)=RF(X)=R,解出X=F−1(R)X=F−1(R)就是符合要求的随机变量。

例子1：指数分布

X是指数分布，f(x)=λe−λx,x≥0f(x)=λe−λx,x≥0

解出F(x)=1−e−λxF(x)=1−e−λx

于是，X=−1λln(1−R)X=−1λln⁡(1−R)

由于1−R1−R与RR是相同的分布，所以也可以写成X=−1λln(R)X=−1λln⁡(R)

例子2：离散的情况

px=px(1−p)1−xpx=px(1−p)1−x

X=[lnRlnp+1]X=[ln⁡Rln⁡p+1]

例子3：不存在解析解的情况

对norm用逆变换法时，
∫X−∞12π−−√e−u2/2du∫−∞X12πe−u2/2du
上式无法求出解析解 ，虽然可以用数值方法解出来，但是仿真中需要大量这样的随机数。需要生成随机数的场景，一般规模很大，对运算速度要求很高。所以没法用了。
用以下方法改进：

包围取舍采样法

问题提出

X∼f(x)X∼f(x)无法求出F−1(x)F−1(x)的解析解
用数值方法运算效率又太低。

算法过程

下面的算法中有很巧妙的思想，请反复读。

step1 ：对f(x)f(x)的定义域分段，
例如，对正态分布可以只生成正实数区域，然后依概率取相反数即可。
也可以有其它分段方法，例如分n段，每段上生成随机数，根据预先求出的每段的概率，每次依概率选择对应的区间段就行了。

step2：只研究step1中的某一段，找到另一个随机变量Y，使得：

fX(x)≤fY(y)fX(x)≤fY(y)在这一区间段上恒成立。
Y上的随机数容易用其它方法生成

step3：

生成一个Y的随机数y
生成一个U(0,1)U(0,1)的随机数r
如果r>fX(y)fY(y)r>fX(y)fY(y)，返回1，否则输出y

经过以上算法后，生成的y就是服从X分布（在选定区间段上的）的随机数

注意：step1中的分段，不是随意分的，要满足以下条件：

使得step2可以实现，也就是说，可以找到满足条件的随机变量Y
预先知道落在每个区间段上的概率。例如，对于正态分布，如果按中心分2段，那么预先知道落在2个区间段上的概率各为0.5。
研究算法，发现fXfX的乘数系数并不重要。
λfX(x)≤fY(y)λfX(x)≤fY(y)，找到这样的λλ即可。
当然，λλ越大，step3中返回1的概率就越小，算法效率越高