计算机模拟第16课

1、5.4 随机数与伪随机数一、随机数与伪随机数 真随机数：不可预计性，不可重复 伪随机数：利用算法或公式产生的随机数。 1. 真随机数的发生源 晶体管噪声发生器、放射粒子计数器等2.伪随机数：例如0，1区间均匀分布的随机数1二、 伪随机数产生方法1、平方取中法 从某个初始的2k位整数开始，求出这个数的平方，去头截尾取其中间2k位作为一个新的随机数，重复以上过程，则得到一列随机数。例1 已知种子x0=3187， 利用平方取中法产生四位数的随机数序列。 (3187)2=10156969 x1=1569 (1569)2=02461761 x2=4617 . .2例2、 已知x0=44，试 产生两位数的

2、随机数序列。 (44)2 =1936 x1=93 (93)2 =8649 x2=64 (64)2 =4096 x3=09 (09)2 =0081 x4=08 (08)2 =0064 x5=06 (06)2 =0036 x6=03 (03)2 =0009 x7=00 (00)2 =0000 x8=00例3、若利用平方取中法得到xi的一个中间值为4500，继续类推，则有以下结果： (4500)2=20250000 xi+1=2500 (2500)2=06250000 xi+2=2500缺点： 重复周期短； 可能产生例外； 较长的随机数序列可能无法通过统计性检验。3平方取中法的函数 r=pfqz(k

3、,x0,n)其中，k：随机数种子位数的一半 x0：随机数种子 n：产生的随机数个数 r: 产生的随机数序列42. 线性同余法 若两整数A，B之差是m的整数倍， 则称A和B按m同余。 A-B=m*k 记为B=A mod m xn+1=(axn+c) mod m线性同余式(0 x m-1)例: 设a=5,c=3,m=16,取x0=7，用线性同余法产生随机数序列。 x0=7 x1=(57+3) mod 16=6 依次有 x2=1 x3=8 x4=11 x5=10 x6=5 x7=12 x8=15 x9=14 x10=9 x11=0 x12=3 x13=2 x14=13 x15=4 x16=7 令 R

4、i=xi/m 即可得到0，1区间分布的随机数。5介绍一种利用素数模乘同余法产生随机数的程序。 xn+1=(axn) mod p axn-xn+1=kp xn+1=axn-kp =axn+k1p-k1p-kp 素数p=ab+c(0bp,0ca)67素数模乘同余法的函数 r=primod(x0,n)其中，x0：随机数种子 n：产生的随机数个数 r: 产生的随机数序列85.5 任意分布的伪随机变量的抽样一、直接抽样法1. 连续随机变量的直接抽样 设连续随机变量 的分布密度函数为f(x),则 分布函数F(x)= ，若F(x)存在反函数 F-1(x)，令为（0，1）区间的一个均匀分布的随机数，再令F()

5、，则F-1() 证明：随机变量x的概率为: 9例1:设某随机变量的分布函数由下式给出, 且产生的均匀分布的随机数为0.1021, 0.2162, 0.7621,现将它们转化为上述分布下的随 机变量。 10设为(0,1)区间的均匀随机数,令F(),则11例2：产生（a,b）区间均匀分布的随机数，已知12 例3：产生指数分布的随机数，已知其分布函数：13142.离散随机变量的直接抽样设离散随机变量X的可能取值为x1,x2, xk 密度函数为Pk=P(X=xk), k=1,2,3 其分布函数为 (1)取为(0,1)区间的均匀随机数(2)求非负整数k，使得满足 F(xk-1)F(xk)(3)令=xk,

6、即为所求随机数15 例4：产生取值可能为0,1,2,3,4的离散分布随机函数，其概率函数为Pk=(k+1)/15, k=0,1,2,3,416例5：产生几何分布的随机数，其概率函数为：其分布函数为：17 当x=0, qx=1; x=,qx=0 18 function jihe(k,q) r=rand(k) ix=alog(r)/alog(q) jihe=ix end二、函数变换法 利用随机变量函数的概率分布特性，通过函数变换关系，可以从均匀分布的随机数产生非均匀分布的随机数，或者从一种分布的随机数产生另外一种分布的随机数。1、随机变量19 设连续型随机变量x的密度函数为P(x),其形式为P(x

7、)=f(x),则其分布函数 不存在反函数形式 作变量替换，令x=g(y)，且 存在 其分布函数20若G(y)存在反函数则令=G(y) y=G-1() x=g(y)=g(G-1()2、二维随机向量 设(1, 2)为一个二维随机向量，其联合密度函数为f (1, 2),对(1, 2)作函数变换 其逆变换存在21 若1, 2存在连续的一阶导数，则函数变换的雅可比行列式为： 例6：用函数变换法产生标准正态分布的随机数 22 解：设 r1,r2 (0,1)区间的两个相互独立的均匀随机数，作函数变换 其逆变换2324 function tran(k1,k2,y1,y2) r1=rand(k1) r2=ran

8、d(k2) p=sqrt(-2.0*alog(r1) q=2*3.1415926*r2 y1=p*cos(q) y2=p*sin(q) tran=y1 end25例2：用函数变换法产生柯西分布的随机数 其中 为实数 , 特别, 当设x1,x2是两个相互独立的标准正态分布的随机 数,作如下变换 ,其逆变换为26 所以关于(y1,y2)的联合密度函数为 27关于y2的边缘密度函数：28subroutine cch(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

9、=y1/y2 z=z*a+u end 293.随机数的检验（经验检验）1)参数检验：均值，方差2)均匀性检验3)独立性检验：Pjj阶自相关系数(j=0不相关) subroutine test(n, k, a, b, ia, sl) dimension a(n),b(k),ia(k),sl(9) sl(1)=0.5 sl(2)=0.5 do 10 I=3,9 sl(i)=0.010 continue an=float(n)30 ak=float(k) do 20 i=1,k ai=float(i) ia(i)=0 b(i)=ai/ak20 continue do 30 j=1,n yf1=a(j

10、) if(yf1.gt.sl(1) sl(1)=yf1 if(yf1.lt.sl(2) sl(2)=yf1 if(yf1-0.5).gt.0.0) sl(6)=sl(6)+1.0 sl(3)=sl(3)+yf1/an 31 sl(4)=sl(4)+yf1*2/an do 40 i=1,k if(yf1.gt.b(i) goto 40 ia(i)=ia(i)+1 goto 3040 continue30 continue sl(7)=sl(6)/an sl(5)=sl(4)-sl(3)*sl(3) do 50 i=1,k b(i)=ia(i)-an/ak x2=b(i)*b(i)32 sl(8)=sl(8)+x2*ak/an50 continue do 70 j=1,20 aj=float(n-j) pj=0 do 60 i=1,n-j pj=pj+a(i)*a(i+j)/aj60 conitnue pj=(pj-sl(3)*sl(3)/sl(5) uj=abs(pj*sqrt(aj) if(uj.gt.1.96) sl(9)=sl(9)+1.0 70 continue end 33 program main dime