第四章 随机数产生原理

第四章随机数产生原理一、引言二、伪随机数产生原理三、[0，1]均匀分布随机数的算法四、其他分布随机数的产生五、正态分布随机数的产生六、MATLAB统计库中的随机数发生器七、随机数的检验八、案例3九、习题一、引言以随机数产生为基础的Monte-Carlo方法已成为现代科研的重要手段之一。其意义早以超出了概率论与数理统计的范畴。广泛应用于计算方法、随机数规划、管理科学、物理化学中的高分子结构的研究等领域。我们来看一些例子。1、数值计算的研究数值计算的研究可以说是一切Monte-Carlo应用的基础，在计算数学领域我们遇到了很多的复杂计算，一个典型的例子是计算积分。对于一维、二维的问题早已获得解决。但是当遇到高维积分问题时，我们传统的数值方法都由于计算量太大而陷于了困境。但是高维积分问题又偏偏是物理、高分子化学、运筹学和最优化问题迫切而必须解决的问题。我们来看一个例子。},,,{21nxxxx这里dxxf)(这是一个众所周知的积分公式，我们当然也可以把它一般的看为是一个高维积分，如果从传统的数值计算方法来看待，则高维问题是随着维数的增加计算量成指数增加的，计算很快就失去控制。但是如果我们换一个角度来看待这个问题，从概率论的角度，它实际是：dxxfxfE)()]([即是f（x）的均值，对于均值我们有一个很好的估计，即NxfxfENii1)()]([ˆ【例4.1.1】用Monte-Carlo对积分102dxx解：将积分区域和值域看成是一个边长为一的正方形。利用均匀分布随机数将点撒在正方形中，计算小于函数的个数并除全部点数。这就是积分的近似值。%利用Monte-Carlo方法计算定积分x=rand(1,1000);x_2=x.^2;JF=mean(x_2)%作Monte-Carlo积分示意图fori=1:1000xx=rand(1,100);yy=rand(1,100);endx1=linspace(0,1,50);y1=x1.^2;plot(xx,yy,'.',x1,y1,'linewidth',2)axisequalh=legend('x-y','x^2');JF=0.3346面积计算结果为：s=0.3482【例4.1.2】利用Monte-Carlo方法计算定积分。10,1022)sin(yxdxdyyx解：抽两组随机数，求每组元素的平方代入给定的函数，然后求平均值即得积分的近似值。%Monte-Carlo方法积分二重积分并与数值方法的结果进行比较Q=dblquad('sin(x.^2+y.^2)',0,1,0,1)%数值求积分命令x=rand(2,100000);%产生两组随机数Sin_xy=sin(x(1,:).^2+x(2,:).^2);%代入函数JF_Sin_xy=mean(Sin_xy)%用Monte-Carlo方法求积分值计算结果为：Q=0.5613JF_Sin_xy=0.5612当抽样数很大时结果很接近。我们可以从Monte-Carlo方法中看出，如果维数增加实际计算难度并不增加，因此是处理高维问题的有效方法。这里x是积分定义域中的均匀分布随机数，这是革命性的一个视角。从这个视角，我们把繁难的积分计算变成了简单的算术平均，而大量的抽样计算又是计算机的拿手好戏，更重要的是当维数增加时并不增加计算难度，从而用Monte-Carlo方法研究高维积分问题已是当今计算数学界的热门课题。2、管理科学的系统仿真研究管理科学中的系统仿真研究，如服务系统、库存系统等。其共性就是研究的对象是随机数，如顾客到达时间一般是一个服从指数分布的随机数，而服务时间也可以看成是服从某种分布的随机数，当一个系统是多队多服务体系时，问题就变的相当复杂了。我们很难用数学的解析式来表达。这时Monte-Carlo方法也是有利的武器。对于这个领域的已有各种比较成熟的专用软件如GPSS、SIMULATION等可以使用。3.物理化学中的分子领域50年代科学家已经在高分子领域使用Monte-Carlo方法了。这一领域所研究的问题更加复杂，计算量非常之大。高分子材料是由几乎是无穷的高分子链组成，而每一个链的长度又是10的好几次方。而链的构象又是千差万别，而且是随机游动的。如何在中微观上几乎是无规律的现象中去判断其宏观的性质？用数学的解析式来说明这样的现象是苍白无力的。Monte-Carlo方法是一个很好的工具，它使得科学家用Monte-Carlo方法去探索高分子运动的规律。一个典型的例子是：对于高分子多链体的研究这是一个很复杂的问题，直到标度理论和重正化群理论方法的引入，才使得单链构象统计问题得到了较好的解决。例：用计算机模拟高分子链链的末端距末端距：空间一链的末端与始端的距离为末端距，由于我们将始端放在坐标原点，所以末端距的计算公式为：末端距=（X2+Y2+Z2)1/2这里X，Y，Z为链的末端点的坐标。显然末端距随链的不同而不同，即为随机变量。-0.500.51-1.5-1-0.50012-1.5-1-0.50三根链的起点(0,0,0)蓝链的末端距绿链的末端距红链的末端距二、伪随机数产生原理前面Monte-Carlo方法的例子是以高质量的随机数为基础的。通过完全的随机抽样或调查可以产生随机序列。当我们用Monte-Carlo方法研究一个实际问题时，我们需要快速地获得大量的随机数。用计算机产生这样的随机数是非常方便的，用数学方法在计算机上产生的随机数称为伪随机数。nxxx,,,21基本定理：如果随机变量X的分布函数F（x）连续，则R=F（x）是[0，1]上的均匀分布的随机变量。证：因为分布函数F（x）是在（0，1）上取值单调递增的连续函数，所以当X在（－∞，x）内取值时，随机变量R则在（0，F(x)）上取值，且对应于（0，1）上的一个R值，至少有一个x满足，见图4.2.1xPxFr)(以表示随机变量R的分布函数，则有)(1rFrfPrRPxF)()(1(4.2.1)1110)(001rrrrFXPr当当当=证毕图4.2.1(4.2.2)基本定理给出了任一随机变量和均匀分布R之间的关系。而有些随机变量可以通过分布函数的逆变换来获得，因此如果我们可以产生高质量的均匀分布，我们就可以通过变换获得高质量的其他分布。见公式(4.2.3))(1rF(4.2.3)例4.2.1求指数分布的随机数。令)1(0xxtedteRxeR1)1ln()1ln(RxxR从而我们用服从[0，1]上的随机数R，通过上面的公式就可以得到指数分布的随机数了。例4.2.1产生1000个均匀分布随机数，利用变换产生λ=6的指数分布并进行拟合优度检验。clc,clearx=linspace(0,20,100);R=rand(1,1000);%产生1000个（0，1）均匀随机数ex=-6*log(1-R);%变换为指数分布随机数ex=ex';m=mean(ex)v=var(ex)subplot(1,2,1),cdfplot(ex)subplot(1,2,2),hist(ex)kstest(ex,[exexpcdf(ex,6)])%拟合优度检验结果为：H=0,接受原假设，变换后的确为λ=6的指数分布三、（0，1）均匀分布随机数的产生1、算法要求（1）产生的数值序列要具有均匀总体简单子样的一些概率统计特性，通常包括分布的均匀性，抽样的随机性，试验的独立性和前后的一致性。（2）产生的随机数要有足够长的周期。（3）产生随机数速度快，占用内存小。为了达到快速的要求，一般采用递推公式),,|,,,(12110knnnknxxxaaafx（4.3.1)目前最常用的方法是上述方法的一个特例：混合同余法MMyxbayynnnn/)(mod1(4.3.2)其中a，b，M以及初值y都是正整数，容易看出x满足0≤x≤1。其中modM运算定义为：任一整数y可唯一表示为公式znMyzMy)(mod则乘同余法当b=0时，有MMyxayynnnn/)(mod1(4.3.4)加同余法以下形式的同余法称为加同余法MMyxyyyynnknnnn/)(mod121(3.4.5)例4.3.1历史上比较有名的称为“菲波那西”数列为加同余法的特例。MMyxyyynnnnn/)(mod21(4.3.6)当M=8时，取初值得“菲波那西”数列。0，1，1，2，3，5，8，13，21，34，55，89，144，253……对上述数列取模得：0，1，1，2，3，5，0，5，5，7，1，1……(4.3.7)再除以模M我们可得到(0,1)之间的序列。我们知道对于一个来自均匀分布的随机序列，应该满足独立性、均匀性等统计特性，但伪随机数往往有一些缺陷，例如(4.3.7)序列到一定长度后，又开始重复下面的序列，这称为周期性是一种明显的规律，与随机性矛盾。通常我们只能选用一个周期内的序列作为我们的伪随机数。因此研究一种算法，使得其产生的随机序列的周期尽可能长，我们可以通过调节（4.3.1）的参数来实现。因此如何来获得一个周期比较长的序列，就成了我们研究的一个内容：有关伪随机数序列的周期有如下的一些结论：定理4.3.1混合同余法产生序列达最大周期M的充要条件：（1）b与M互素（2）对于M的任意素因子p，有（3）如果4是M的因子，则显然乘同余法产生的序列达不到周期M（不满足（1））。当取（k为任意整数）时，因为M只有一个素因子2，且4是M的因子，则由条件（2）、（3）有，从而混合同余发生器达到最大周期的算法为：)(mod1pa)4(mod1akm2)4(mod1akiikiiyxdycy2/)2)(mod12()14(1（3.4.8)其中c，d为任意整数。混合同余发生器是否达到最大周期M与初始值无关。),2,1)((modiMbaii对于乘同余发生器，由同余运算的定义，知其由如下性质(1)如果则有：)(mod)(mod21212211MbbaaMbaba),(modMcbca（2）如果),(modMcMba则其中（c，M）是c，M的最大公约数。利用这些性质可得到以下定理。定理4.3.2对乘同余发生器，若，则使成立得最小正整数V就是此发生器得周期。1),(0Mx)(mod1MaV在数论中称V为a关于M的阶数，对于乘同余发生器，若种子与M互素，则其周期就是关于M的阶数。这样一来，寻找达到最大周期的同余发生器的问题就转化为数论方面寻求M达到最大阶数a的问题了。Knuth对这一问题的研究作了总结。从算法上我们知道公式是递推的，因此一般的随机发生器程序都要预先赋初值，这种初值为种（Seed），有些统计软件如SPSS要求用户给出Seed.以均匀分布（0，1）随机变量R变换成的随机变量。以r,ε,u,分别表示(0,1)均匀分布,指数分布，N(0,1)标准正态分布。其他常见的分布如卡方分布、F分布等的抽样方法见表4.3.1。四、其他分布随机数的产生1、直接抽样法由基本定理我们知道，对于有些随机变量可以建立与R的函数关系，因此只需对R进行抽样，利用函数的映射关系我们就可以方便地得到该随机变量的抽样了。如前面的指数分布随机数。2、变换抽样产生随机变量的变换抽样方法，是讨论均匀分布的不同函数分布，为随机变量抽样提供一些简单可行的算法。在概率论中，从不同的角度出发，对随机变量函数进行了讨论，以下列出一些结果。设随机变量X具有密度函数，是对随机变量X的变换，且的逆函数存在，记为有一阶连续导数，则随机变量的密度函数)(xf)(XgY)(xg，)()(1yhxg)(XgY)()()('*yhyhfyf（4.4.1)例4.4.1R，1-R均为（0，1）上的均匀分布随机变量，设随机向量（X，Y）具有二维联合密度。对于随机变量X，Y进行函数变换)