随机数的产生方法
均匀分布
线性同余法
线性同余法(LCG) 是根据递归公式:Nj+1≡(A x Nj + B)(mod M)①来确定一随机数列,公式中A,B,M是产生器设定的常数,模运算mod的定义如下:
任一正整数y可唯一表示为y=n*M+z,M、z、n均为整数,0<=z<M,则y(mod M)=z;
LCG的最大周期为M,但大部分情况都会少于M。要令LCG达到最大周期,应符合以下条件:
- B,M互质
- M的所有质因数都可以整除A-1;
- 若M是4的倍数,A-1也是;
- A,B,N0都比M小;
- A,B是正整数。
事实上mod即为取余运算,针对随机生成的数列Nk,显然,有0<=Nk<M,且Nk为整数,又因为递推公式Nj+1≡(A x Nj + B)(mod M),则不难理解,唯一的Nk可决定唯一的Nk+1。数列Nk必是存在周期性的,而且周期T<=M。
为产生(0,1)之间的均匀分布,可令M=134456,A=8121,C=28411,此时满足上述条件,即所生成数列的周期T=M=134456。由于数列中每个数的出现概率相等,都为1/134456;而且本文为探究性质,故可将其视作符合精度要求。为实现输出(0,1)之间的均匀分布,需将①式修改为:Nj+1≡(A x Nj + B)(mod M)/M,运行程序,可得如下结果:
上图所示的“类函数分布图”中,横坐标代表数据的值x,纵坐标代表在总数为500的情况下,小于或等于值x的数据的个数,横坐标x的值从左向右递增。可以看出,“类函数分布图”曲线形状与对应分布函数图像相似,也可再一定程度上反应出样本所服从的分布的一些特征。又由样本的均值为x=0.517687, 方差σ2=0.082915,而(0,1)之间均匀分布的均值、方差理论值为EX=0.5, VAR(X)=1/12=0.83333,相差在可以接受的范围内,模拟成功
在实现(0,1)之间均匀分布的随机数的生成之后,接下来试图通过此均匀分布生成常用的分布:正态分布、泊松分布,以及指数分布
正态分布
中心极限定理
此定理不论随机变量x原先服从什么分布都成立。根据此定理,因代入μ=0.5,σ2=1/12,则有u=sqrt(12n)*[(Σri)/n -0.5]
为编写程序方便,取n=12,则有
编程实现,运行结果如下:
样本的均值为x=0.0003357 方差σ2 =1.006824 标准正态分布N(0,1)的均值、方差理论值为EX=0,VAR(X)=1, 相差在可以接受的范围内,模拟成功
事实上,根据所生成的标准正态分布,可以随之生成很多相关分布,比如令Y=eX,则根据定义容易看出Y服从对数正态分布Y ~ LN(0,1),而令Y=X2,例题中也讲述过此时Y服从卡方分布Y~χ2(1)
指数分布
依据定理:若随机变量的分布函数Fx(x)为严格单调递增的连续函数,其反函数F_x^(-1)(y)存在,则Y=Fx(X)服从(0,1)上的均匀分布U(0,1)。
此定理给出了任意随机变量和均匀分布之间R之间的关系,根据此定理,一些随机变量可以根据其分布函数的逆变换来通过均匀分布生成出。
求解出x=-λ1ln(1-R),不妨令λ1=5,编写代码,运行结果如下
根据λ1=5,可计算出此指数分布参数λ=1/5=0.2,样本的均值为x=5.287232 方差σ2=26.87719 而指数分布的均值、方差理论值为EX=5,VAR(X)=25, 可以看出,样本与理论相差还较大,简要分析,可能是均匀分布效果不稳定,以致出现此结果。
泊松分布
利用一连续型的随机分布生成离散型泊松分布,算法如下:
algorithm poisson random number (Knuth):
init:
Let L ← e^{−λ}, k ← 0 and p ← 1.
do:
k ← k + 1.
Generate uniform random number u in [0,1] and let p ← p × u.
while p > L.
return k − 1.
算法可以描述为:
第一步:给定一个参数, 生成一系列的随机数,这些随机数服从(0, 1) 之间均匀分布。
第二步:求这些随机数的乘积,当乘积小于或者等于e-λ时,程序停止。记下此时参与求乘积的随机数的个数。
第三步:程序终止时参与乘积的随机数的个数减去一即服从参数为λ的泊松分布。
此段算法的详细证明因本人也是看相关资料才搞明白的,故不再加以细述,只文末附资料链接。取λ=5,编写代码,运行程序后有如下结果:
样本的均值为x=5.044 方差σ2=4.947011 而指数分布的均值、方差理论值为EX=5,VAR(X)=5, 可以看出,样本与理论相差在可接受的范围内,故模拟成功。