计算机模拟第16课

5.4随机数与伪随机数,一、随机数与伪随机数真随机数：不可预计性，不可重复伪随机数：利用算法或公式产生的随机数。1.真随机数的发生源晶体管噪声发生器、放射粒子计数器等2.伪随机数：例如0，1区间均匀分布的随机数,二、伪随机数产生方法1、平方取中法从某个初始的2k位整数开始，求出这个数的平方，去头截尾取其中间2k位作为一个新的随机数，重复以上过程，则得到一列随机数。例1已知种子x0=3187，利用平方取中法产生四位数的随机数序列。(3187)2=10156969x1=1569(1569)2=02461761x2=4617.,例2、已知x0=44，试产生两位数的随机数序列。(44)2=1936x1=93(93)2=8649x2=64(64)2=4096x3=09(09)2=0081x4=08(08)2=0064x5=06(06)2=0036x6=03(03)2=0009x7=00(00)2=0000 x8=00例3、若利用平方取中法得到xi的一个中间值为4500，继续类推，则有以下结果：(4500)2=20250000 xi+1=2500(2500)2=06250000 xi+2=2500缺点：重复周期短；可能产生例外；较长的随机数序列可能无法通过统计性检验。,平方取中法的函数r=pfqz(k,x0,n)其中，k：随机数种子位数的一半x0：随机数种子n：产生的随机数个数r:产生的随机数序列,2.线性同余法若两整数A，B之差是m的整数倍，则称A和B按m同余。A-B=m*k记为B=Amodmxn+1=(axn+c)modm线性同余式(0xm-1)例:设a=5,c=3,m=16,取x0=7，用线性同余法产生随机数序列。x0=7x1=(57+3)mod16=6依次有x2=1x3=8x4=11x5=10 x6=5x7=12x8=15x9=14x10=9x11=0 x12=3x13=2x14=13x15=4x16=7令Ri=xi/m即可得到0，1区间分布的随机数。,介绍一种利用素数模乘同余法产生随机数的程序。xn+1=(axn)modpaxn-xn+1=kpxn+1=axn-kp=axn+k1p-k1p-kp素数p=ab+c(0bp,0ca),素数模乘同余法的函数r=primod(x0,n)其中，x0：随机数种子n：产生的随机数个数r:产生的随机数序列,5.5任意分布的伪随机变量的抽样,一、直接抽样法1.连续随机变量的直接抽样设连续随机变量的分布密度函数为f(x),则分布函数F(x)=，若F(x)存在反函数F-1(x)，令为（0，1）区间的一个均匀分布的随机数，再令F()，则F-1()证明：随机变量x的概率为:,例1:设某随机变量的分布函数由下式给出,且产生的均匀分布的随机数为0.1021,0.2162,0.7621,现将它们转化为上述分布下的随机变量。,设为(0,1)区间的均匀随机数,令F(),则,例2：产生（a,b）区间均匀分布的随机数，已知,例3：产生指数分布的随机数，已知其分布函数：,2.离散随机变量的直接抽样设离散随机变量X的可能取值为x1,x2,xk密度函数为Pk=P(X=xk),k=1,2,3其分布函数为(1)取为(0,1)区间的均匀随机数(2)求非负整数k，使得满足F(xk-1)F(xk)(3)令=xk,即为所求随机数,例4：产生取值可能为0,1,2,3,4的离散分布随机函数，其概率函数为Pk=(k+1)/15,k=0,1,2,3,4,例5：产生几何分布的随机数，其概率函数为：其分布函数为：,当x=0,qx=1;x=,qx=0,functionjihe(k,q)r=rand(k)ix=alog(r)/alog(q)jihe=ixend二、函数变换法利用随机变量函数的概率分布特性，通过函数变换关系，可以从均匀分布的随机数产生非均匀分布的随机数，或者从一种分布的随机数产生另外一种分布的随机数。1、随机变量,设连续型随机变量x的密度函数为P(x),其形式为P(x)=f(x),则其分布函数不存在反函数形式作变量替换，令x=g(y)，且存在其分布函数,若G(y)存在反函数则令=G(y)y=G-1()x=g(y)=g(G-1()2、二维随机向量设(1,2)为一个二维随机向量，其联合密度函数为f(1,2),对(1,2)作函数变换其逆变换存在,若1,2存在连续的一阶导数，则函数变换的雅可比行列式为：例6：用函数变换法产生标准正态分布的随机数,解：设r1,r2(0,1)区间的两个相互独立的均匀随机数，作函数变换其逆变换,functiontran(k1,k2,y1,y2)r1=rand(k1)r2=rand(k2)p=sqrt(-2.0*alog(r1)q=2*3.1415926*r2y1=p*cos(q)y2=p*sin(q)tran=y1end,例2：用函数变换法产生柯西分布的随机数其中为实数,特别,当设x1,x2是两个相互独立的标准正态分布的随机数,作如下变换,其逆变换为,所以关于(y1,y2)的联合密度函数为,关于y2的边缘密度函数：,subroutinecch(a,u,k1,k2,y1,y2,z)p=sqrt(-2.0*alog(rand(k1)q=2*3.1415926*rand(k2)y1=p*cos(q)y2=p*sin(q)z=y1/y2z=z*a+uend,3.随机数的检验（经验检验）1)参数检验：均值，方差2)均匀性检验3)独立性检验：Pjj阶自相关系数(j=0不相关)subroutinetest(n,k,a,b,ia,sl)dimensiona(n),b(k),ia(k),sl(9)sl(1)=0.5sl(2)=0.5do10I=3,9sl(i)=0.010continuean=float(n),ak=float(k)do20i=1,kai=float(i)ia(i)=0b(i)=ai/ak20continuedo30j=1,nyf1=a(j)if(yf1.gt.sl(1)sl(1)=yf1if(yf1.lt.sl(2)sl(2)=yf1if(yf1-0.5).gt.0.0)sl(6)=sl(6)+1.0sl(3)=sl(3)+yf1/an,sl(4)=sl(4)+yf1*2/ando40i=1,kif(yf1.gt.b(i)goto40ia(i)=ia(i)+1goto3040continue30continuesl(7)=sl(6)/ansl(5)=sl(4)-sl(3)*sl(3)do50i=1,kb(i)=ia(i)-an/akx2=b(i)*b(i),sl(8)=sl(8)+x2*ak/an50continuedo70j=1,20aj=float(n-j)pj=0do60i=1,n-jpj=pj+a(i)*a(i+j)/aj60conitnuepj=(pj-sl(3)*sl(3)/sl(5)uj=abs(pj*sqrt(aj)if(uj.gt.1.96)sl(9)=sl(9)+1.070continueend,programmaindim