随机数产生原理.ppt

第四章 随机数产生原理一、引言 二、伪随机数产生原理 三、[0，1]均匀分布随机数的算法 四、其他分布随机数的产生 五、正态分布随机数的产生 六、MATLAB统计库中的随机数发生器 七、随机数的检验 八、案例3 九、习题一、引言n以随机数产生为基础的Monte-Carlo方法已成为现代科研的重要手段之一其意义早以超出了概率论与数理统计的范畴广泛应用于计算方法、随机数规划、管理科学、物理化学中的高分子结构的研究等领域我们来看一些例子1、数值计算的研究数值计算的研究可以说是一切Monte-Carlo应用的基础，在计算数学领域我们遇到了很多的复杂计算，一个典型的例子是计算积分对于一维、二维的问题早已获得解决但是当遇到高维积分问题时，我们传统的数值方法都由于计算量太大而陷于了困境但是高维积分问题又偏偏是物理、高分子化学、运筹学和最优化问题迫切而必须解决的问题我们来看一个例子这里这是一个众所周知的积分公式，我们当然也可以把它一般的看为是一个高维积分，如果从传统的数值计算方法来看待，则高维问题是随着维数的增加计算量成指数增加的，计算很快就失去控制但是如果我们换一个角度来看待这个问题，从概率论的角度，它实际是：即是f（x）的均值，对于均值我们有一个很好的估计，即【例4.1.1】 用Monte-Carlo 对 积分解：将积分区域和值域看成是一个边长为一的正方形。

利用均匀分 布随机数将点撒在正方形中，计算小于函数的个数并除全部点数 这就是积分的近似值 利用Monte-Carlo方法计算定积分 x=rand(1,1000); x_2=x.^2; JF=mean(x_2) % 作Monte-Carlo积分示意图 for i=1:1000xx=rand(1,100);yy=rand(1,100); end x1=linspace(0,1,50); y1=x1.^2; plot(xx,yy,'.',x1,y1,'linewidth',2) axis equal h = legend('x-y','x^2'); JF = 0.3346面积计算结果为： s = 0.3482【例4.1.2】 利用Monte-Carlo方法计算定积分解：抽两组随机数，求每组元素的平方代入给定的函数，然后求平 均值即得积分的近似值 % 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.5613 JF_Sin_xy = 0.5612 当抽样数很大时结果很接近。

我们可以从Monte-Carlo方法中看出 ，如果维数增加实际计算难度并不增加，因此是处理高维问题的有 效方法这里 x 是积分定义域中的均匀分布随机数，这是革命性的一个视角从这个视角，我们把繁难的积分计算变成了简单的算术平均，而大量的抽样计算又是计算机的拿手好戏，更重要的是当维数增加时并不增加计算难度，从而用 Monte-Carlo 方法研究高维积分问题已是当今计算数学界的热门课题2、管理科学的系统仿真研究 管理科学中的系统仿真研究，如服务系统、库存系统等其共性就是研究的对象是随机数，如顾客到达时间一般是一个服从指数分布的随机数，而服务时间也可以看成是服从某种分布的随机数，当一个系统是多队多服务体系时，问题就变的相当复杂了我们很难用数学的解析式来表达这时Monte-Carlo方法也是有利的武器对于这个领域的已有各种比较成熟的专用软件如GPSS、SIMULATION等可以使用3. 物理化学中的分子领域50年代科学家已经在高分子领域使用Monte-Carlo方法了这一领域所研究的问题更加复杂，计算量非常之大高分子材料是由几乎是无穷的高分子链组成，而每一个链的长度又是10的好几次方而链的构象又是千差万别，而且是随机游动的。

如何在中微观上几乎是无规律的现象中去判断其宏观的性质？用数学的解析式来说明这样的现象是苍白无力的Monte-Carlo方法是一个很好的工具，它使得科学家用Monte-Carlo方法去探索高分子运动的规律一个典型的例子是：对于高分子多链体的研究这是一个很复杂的问题，直到标度理论和重正化群理论方法的引入，才使得单链构象统计问题得到了较好的解决例：用计算机模拟高分子链链的末端距n末端距：空间一链的末端与 始端的距离为末端距，由于 我们将始端放在坐标原点， 所以末端距的 计算公式为 ：n末端距=（X2+Y2+Z2)1/2n这里X，Y，Z为链的末端点 的坐标n显然末端距随链的不同而不同，即为随机变量三根链的起点 (0,0,0)蓝链的末端距绿链的末端距红链的末端距二、伪随机数产生原理 n前面Monte-Carlo方法的例子是以高质量的随机数为 基础的通过完全的随机抽样或调查可以产生随机序 列n当我们用Monte-Carlo方法研究一个实际问题时，我 们需要快速地获得大量的随机数用计算机产生这样 的随机数是非常方便的，用数学方法在计算机上产生 的随机数称为伪随机数 基本定理：如果随机变量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以 表示随机变量 R 的分布函数，则有(4.2.1)=证毕图4.2.1(4.2.2)基本定理给出了任一随机变量和均匀分布R之间的关系而有些随机变量可以通过分布函数的逆变换来获得，因此如果我们可以产生高质量的均匀分布，我们就可以通过变换获得高质量 的其他分布见公式 (4.2.3)(4.2.3)例4.2.1 求指数分布的随机数令 从而我们用服从[0，1]上的随机数R，通过上面的公式就可以得到指数分布的随机数了例4.2.1 产生1000个均匀分布随机数，利用变换产生λ=6的指数 分布并进行拟合优度检验clc,clear x=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, [ex expcdf(ex, 6)]) % 拟合优度检验结果为：H=0, 接受原假设，变换后的确为λ=6的指数分布三、 （0，1）均匀分布随机数的产生n1、算法要求n（1） 产生的数值序列要具有均匀总体简 单子样的一些概率统计特性，通常包括分布 的均匀性，抽样的随机性，试验的独立性和 前后的一致性。

n（2） 产生的随机数要有足够长的周期n（3） 产生随机数速度快，占用内存小为了达到快速的要求，一般采用递推公式（4.3.1)目前最常用的方法是上述方法的一个特例： 混合同余法(4.3.2)其中a，b，M 以及初值 y 都是正整数，容易看出 x 满足0≤x≤1其中 mod M 运算定义为：任一整数 y 可唯一表示为公式则乘同余法当 b = 0 时，有 (4.3.4)加同余法 以下形式的同余法称为加同余法 (3.4.5)例4.3.1 历史上比较有名的称为“菲波那西”数列为加同余法的特例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）有 ，从而混合同余发生器达到最大周期的算法为：（3.4.8)其中c，d为任意整数混合同余发生器是否达到最大周期M与初始值无关对于乘同余发生器，由同余运算的定义，知其由如下性质(1) 如果则有：（2）如果则其中（c，M）是c，M的最大公约数利用这些性质可得到以下定理定理4.3.2 对乘同余发生器，若 ，则使成立得最小正整数 V 就是此发生器得周期。

在数论中称 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四、其他分布随机数的产生 n1、 直接抽样法n 由基本定理我们知道，对于有些随机变量可以建立 与R的函数关系，因此只需对R进行抽样，利用函数的 映射关系我们就可以方便地得到该随机变量的抽样了 如前面的指数分布随机数n2、变换抽样n 产生随机变量的变换抽样方法，是讨论均匀分布的 不同函数分布，为随机变量抽样提供一些简单可行的算 法在概率论中，从不同的角度出发，对随机变量函数 进行了讨论，以下列出一些结果设随机变量 X 具有密度函数 ，是对随机变量 X 的变换，且 的逆函数存在，记为 有一阶连续导数，则随机变量的密度函数，（4.4.1)例4.4.1 R，1-R均为（0，1）上的均匀分布随机变量，设随机向量（X，Y）具有二维联合密度。

对于随机变量 X，Y 进行函数变换其中，函数 ， 的逆变换存在，记为且存在一阶偏导数，设 J 为 Jacobi 变换则随机变量的二维联合密度为(4.4.2)例4.4.2 用变换抽样产生标准正态分布的随机变量 U 随机变量U 的密度函数为:取随机数，则相互独立、服从N（0，1）分布解上面的两个方程，得逆变 换公式由(4.4.2)，随机变。