计算机中随机数产生方法

随机数的产生方法

均匀分布

线性同余法

线性同余法(LCG) 是根据递归公式：N_j+1≡（A x N_j + 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即为取余运算，针对随机生成的数列N_k，显然，有0<=N_k<M,且Nk为整数，又因为递推公式N_j+1≡（A x N_j + B）(mod M)，则不难理解，唯一的N_k可决定唯一的N_k+1。数列N_k必是存在周期性的，而且周期T<=M。
为产生（0,1）之间的均匀分布，可令M=134456，A=8121,C=28411,此时满足上述条件，即所生成数列的周期T=M=134456。由于数列中每个数的出现概率相等，都为1/134456；而且本文为探究性质，故可将其视作符合精度要求。为实现输出(0,1)之间的均匀分布，需将①式修改为：N_j+1≡（A x N_j + 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，σ²=1/12，则有u=sqrt(12n)*[(Σr_i)/n -0.5]
为编写程序方便，取n=12，则有

编程实现，运行结果如下：


样本的均值为x=0.0003357 方差σ² =1.006824 标准正态分布N(0,1)的均值、方差理论值为EX=0,VAR(X)=1, 相差在可以接受的范围内，模拟成功
事实上，根据所生成的标准正态分布，可以随之生成很多相关分布，比如令Y=e^X，则根据定义容易看出Y服从对数正态分布Y ~ LN(0,1),而令Y=X²,例题中也讲述过此时Y服从卡方分布Y~χ²(1)

指数分布

依据定理：若随机变量的分布函数Fx（x）为严格单调递增的连续函数，其反函数F_x^(-1)(y)存在，则Y=Fx（X）服从（0，1）上的均匀分布U(0,1)。
此定理给出了任意随机变量和均匀分布之间R之间的关系，根据此定理，一些随机变量可以根据其分布函数的逆变换来通过均匀分布生成出。

求解出x=-λ₁ln（1-R），不妨令λ₁=5，编写代码，运行结果如下


根据λ₁=5，可计算出此指数分布参数λ=1/5=0.2，样本的均值为x=5.287232 方差σ²=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 方差σ²=4.947011 而指数分布的均值、方差理论值为EX=5,VAR(X)=5, 可以看出，样本与理论相差在可接受的范围内，故模拟成功。