随机数的生成及随机变量抽样 实验目的 实验内容 学习主要的随机变量抽样方法 1、均匀分布U(0,1)的随机数的产生 2、其他各种分布的随机数的产生方法 3、随机数生成实例 4、实验作业

随机数的生成 随机数的产生是实现MC计算的先决条件。而大多数概率分布的随机数的产生都是基于均匀分布U(0,1)的随机数。 其次，介绍服从其他各种分布的随机数的产生方法。以及服从正态分布的随机数的产生方法。 最后，关于随机数的几点注。

一、均匀分布U(0,1)的随机数的产生 产生均匀分布的标准算法在很多高级计算机语言的书都可以看到。算法简单，容易实现。使用者可以自己手动编程实现。Matlab 中也提供给我们用于产生均匀分布的各种函数。我们的重点是怎样通过均匀分布产生服从其他分布的随机数。因此，直接使用Matlab提供的可靠安全的标准函数，当然不用费事了。

RNUN/DRNUN: 产生均匀分布的随机数 IMSL库中的函数使用 RNSET: 种子的设定 CALL RNSET (ISEED) RNOPT: 产生器的类型的设定 CALL RNOPT (IOPT) RNUN/DRNUN: 产生均匀分布的随机数 CALL RNUN (NR, R)

生成1行1000列501—1000上离散均匀分布的随机数。 并画经验分布函数曲线。 例1生成1行1000列的1—10上离散均匀分布的随机数； 生成1行1000列21—30上离散均匀分布的随机数； 生成1行1000列501—1000上离散均匀分布的随机数。 并画经验分布函数曲线。 cdfplot(x) Randnum=unidrnd(10,1,10000)；cdfplot(Randnum)； pause Randnum=unidrnd(10,1,10000)+10；cdfplot(Randnum)； Randnum=unidrnd(500,1,10000)+500； cdfplot(Randnum)

例2设总体X的密度函数为 其中 >0, 生成 1行10000列的随机数. 并画经验分布函数曲线。 解:由密度函数知 其中 >0, 生成 1行10000列的随机数. 并画经验分布函数曲线。 解:由密度函数知 具有均值为 的指数分布 Randnum=exprnd(2,1,10000)+5 cdfplot(Randnum)

二、其他各种分布的随机数的产生 基本方法有如下三种： 逆变换法 合成法 筛选法

逆变换法 设随机变量 的分布函数为 ，定义 定理 设随机变量 服从 上的均匀分布，则 的分布函数为 。 设随机变量 的分布函数为 ，定义 定理 设随机变量 服从 上的均匀分布，则 的分布函数为 。 因此，要产生来自 的随机数，只要先产生来自 的随机数，然后计算 即可。 其步骤为

例3 设密度函数为 为常数 并画经验分布函数曲线。

例4 设X分布函数为F(X) 生成n=20的1行10000列随机数，并画经验分布函数曲线。

n=20 Randnum=1-(1-unifrnd(0,1,1,10000)).^(1/n); cdfplot(Randnum)

例5 设密度函数为 为常数 并画经验分布函数曲线。

合成法 合成法的应用最早见于Butlter 的书中。 构思如下： 如果 的密度函数 难于抽样，而 关于 的条件密度函数 以及 的密度函数 均易于抽样，则 的随机数可如下产生： 可以证明由此得到 的服从 。

筛选抽样 假设我们要从 抽样，如果可以将 表示成 ，其中 是一个密度函数且易于抽样，而 ， 是常数，则 的抽样可如下进行： 假设我们要从 抽样，如果可以将 表示成 ，其中 是一个密度函数且易于抽样，而 ， 是常数，则 的抽样可如下进行： 定理 设 的密度函数 ，且 ，其中 ， ， 是一个密度函数。令 和 分别服从 和 ，则在 的条件 下， 的条件密度为

三、生成标准正态分布的随机数 的随机数产生方法很多。简要介绍三种。 法1、 变换法（Box 和Muller 1958） 设 ， 是独立同分布的 变量，令 则 与 独立，均服从标准正态分布。 法2、 结合合成法与筛选法。（略） 法3、 近似方法（利用中心极限定理） 即用 个 变量产生一个 变量。 其中 是抽自 的随机数， 可近似为一 个 变量。

例6 生成单位圆上均匀分布的1行10000列随机数，并画经验分布函数曲线。 Randnum=unifrnd(0,2*pi,1,10000); xRandnum=cos(Randnum) [Y,II] =sort(xRandnum) yRandnum=sin(Randnum) plot(xRandnum(II),yRandnum(II),'.')

生成单位正方形上均匀分布的1行10000列随机数，并画散点图。 例7 生成单位正方形上均匀分布的1行10000列随机数，并画散点图。 mm=10000;Randnum=unifrnd(0,4,1,mm); xRandnum=zeros(1,mm);yRandnum=zeros(1,mm); for ii=1:mm if Randnum(1,ii)<=1 xRandnum(1,ii)=0; yRandnum(1,ii)=Randnum(1,ii); else if Randnum(1,ii)<=2 xRandnum(1,ii)=Randnum(1,ii)-1; yRandnum(1,ii)=1; else if Randnum(1,ii)<=3 xRandnum(1,ii)=1; yRandnum(1,ii)=1-(Randnum(1,ii)-2); else xRandnum(1,ii)=1-(Randnum(1,ii)-3); yRandnum(1,ii)=0; end end end end [Y,JJ] =sort(xRandnum);plot(xRandnum(JJ),yRandnum(JJ),'.')

离散型随机变量的生成 离散型随机变量X,它的取值是非光滑连续的值，它只能间断地即离散地取值x1,x2,x3，…，xn，且规定x1＜x2＜x3＜…<xn。其概率密度函数为 p(xi)=p{X= xi} 概率分布函数为 例10 对某车间每天需求某种零件的数量历史数据中统计获得表1的结果。生成1行1000列零件需求的随机数。并画经验分布函数曲线。 表1 某零件每天需求量 X

③若u＞F(xi)，转回到第②步，否则转至④； ④输出得 X＝xi。 概率P(x) 累积概率F(x) 可分配的随机数范围 X1=10 0.10 F(X1)=0.10 （ .00 -- .10] X2=20 0.20 F(X2)=0.30 （ .10 -- .30] X3=30 0.40 F(X3)=0.70 （ .30 -- .70] X4=40 0.25 F(X4)=0.95 （ .70 -- .95] X5=50 0.05 F(X5)=1.00 （ .95 – 1) 随机变量生成的算法为 ①产生一个u(0，1)，并令i=0； ②令i=i+1； ③若u＞F(xi)，转回到第②步，否则转至④； ④输出得 X＝xi。

mm=10000;Randnum=unifrnd(0,1,1,mm);xRandnum=zeros(1,mm); for ii=1:mm if Randnum(1,ii)<=0.1 xRandnum(1,ii)=10; else if Randnum(1,ii)<=0.3 xRandnum(1,ii)=20; else if Randnum(1,ii)<=0.7 xRandnum(1,ii)=30; else if Randnum(1,ii)<=0.95 xRandnum(1,ii)=40; else xRandnum(1,ii)=50; end cdfplot(xRandnum)

三角分布(a,m,b)的随机变量其密度函数为 其分布函数为

随机向量的抽样方法 在用Monte Carlo等方法解应用问题时,随机向量的抽样也是经常用到的. 若随机向量各分量相互独立,则它等价于多个一元随机变量的抽样。

例8 生成单位正方形内均匀分布的1行10000列随机数，并画散点图。 mm=10000 xRandnum=unifrnd(0,1,1,mm); yRandnum=unifrnd(0,1,1,mm); plot(xRandnum,yRandnum,'.')

mm=100000 xRandnum=unifrnd(0,1,1,mm); yRandnum=unifrnd(0,1,1,mm); [Y,JJ] =sort(xRandnum) plot(xRandnum(JJ),yRandnum(JJ),'.')

生成单位圆内均匀分布的1行10000列随机数，并画散点图。 例9 生成单位圆内均匀分布的1行10000列随机数，并画散点图。 mm=10000； Randnum1=unifrnd(-1,1,1,2*mm); Randnum2=unifrnd(-1,1,1,2*mm); xRandnum=zeros(1,mm);yRandnum=zeros(1,mm); s=Randnum1.^2+Randnum2.^2；ii=1;jj=1; while ii<mm if s(1,jj)<=1; xRandnum(1,ii)=Randnum1(1,jj); yRandnum(1,ii)=Randnum2(1,jj); ii=ii+1; end jj=jj+1; plot(xRandnum,yRandnum,'.')

关于随机数的几点注 注1 由于均匀分布的随机数的产生总是采用某个确定的模型进行的，从理论上讲，总会有周期现象出现的。初值确定后，所有随机数也随之确定，并不满足真正随机数的要求。因此通常把由数学方法产生的随机数成为伪随机数。 但其周期又相当长，在实际应用中几乎不可能出现。因此，这种由计算机产生的伪随机数可以当作真正的随机数来处理。 注2 应对所产生的伪随机数作各种统计检验，如独立性检验，分布检验，功率谱检验等等。

作业: 1.生成单位球内均匀分布的1行10000列随机数，并画散点图。 2.设密度函数为 为常数 并画经验分布函数曲线。

作业: 3.生成三角分布(0,1,2)的1行10000列随机数，并画散点图。 并画经验分布函数曲线。

§5.2随机数与随机变量的生成 ui的数学期望和方差分别为 § 5.2.1随机数的生成  § 5.2.1随机数的生成 在系统模拟中只要有随机变量，则在模拟运行的每一步中都要对随机变量确定一个具体的值。我们将会遇到各种概率分布的随机变量，但其中最简单或最基本的随机变量是在(0， 1)区间上均匀分布的随机变量。服从某一分布的随机变量都可以通过对(0， 1)均匀分布的随机变量进行适当转换而得到。(0， 1)均匀分布的随机变量的取值也是在(0，1)区间上均匀分布的随机数ui序列(流)的独立采样，其密度函数是 ui的数学期望和方差分别为

因此，若能获得(0，1)均匀分布的随机数，也就能通过对其适当的转换而获得某一规定分布的随机变量的取值，这就是随机变量的生成。为此，首先要掌握(0，1)区间上均匀分布随机数的生成方法。 均匀分布随机数必须具备均匀性和独立性的要求；要生成符合上述要求的随机数流，现在多用数学算法来产生，一般是采用递推算法，确定一个初始值(种子数)以后，逐次递推算得随机数流。数学算法获得的随机数、常称之为伪随机数(Pseudo Random Number)序列。 数学方法计算产生的随机数流必须满足下列要求： (1)尽可能在(0，1)区间均匀分布； (2)具有统计上的独立性； (3)产生的随机数流能够重复出现，即给以相同的初值(种子数)能获得相同的随机数流； (4)有足够长的周期，即在出现周期性重复之前，能生成足够多个的随机数； (5)算法占用计算机内存较少而计算生成速度较快。 目前广泛应用的算法是线性同余法( Linear congruential Method)，其中又分为： 1．混合线性同余法。它是由Lehmer于1951年提出的，其算式为

xi+1=(axi+c) mod m ui+1=xi+1/m 式中 a——乘数(常数)； C——增量(常数)； x0——种子数； m——模数。 a,c,m和x0的选取对随机数流的统计特性和周期长度有极大影响。 上述第一式的含义是 式中[]表示取整数， a，c，m皆为整常数。 2、  乘法线性同余法。若混合线性同余法中c=0，则为乘法线性同余法，其算式为 x i+1=ax imod m u i+1=xi+1/m (5.7)

可参考选用的数据有： （1） a=16807, m=2147483647, x0=123457; （2）    a=655393, m=33554432。 § 5.2.2随机数流的检验 一、均匀分布性检验 1．参数检验。检验ui的数字特征，如均值、方差的估计值和其理论值的差异是否显著。 设有u1, u2,…，un随机数流，则它们的 若 ui序列在(0，1)上均匀分布，可假设：u的期望和方差分别为

2的期望和方差分别为 则上列假设(8．10)与(8．11)应该成立。据此，可对n个ui计算下列统计量 若取显著性水平a=0．05， 当|V1|≤1．96时。则可认为假设(8．10)式成立； 当|V2|≤1．96时，则可认为假设(8．11)式成立。因而可以接受此假设，检验通过;否则拒绝接受。 2．均匀性检验。它是检验所生成的随机数落在(0， 1)各子区间的频率的均匀程度，是否与理论上的均匀分布频率有显著性差异。此处介绍常用方法之一，x2检验方法如下： 将(0，1)区间划分为相等的k个子区间，假如落在第i个( i＝l，2，3，…，k)子区间的随机数有ni个；而在理论上第i个子区间的随机数个数为mi = N／k，其中 N为随机数流总个数(拟检验的)。由此，可计算 x2统计量

再按k-1为自由度、显著性水平取0．05，查得2(a)表值。当算得统计量x2≤2(a)时，可认为在显著性水平a下能接受为均匀分布假设。 二、独立性检验 独立性检验是检验随机数流中前后各数之间是否存在相关性。常用的方法是进行自相关检验。此外，还有Poker Test和Run检验，一般应用较少。   § 5.2.3随机变量的生成 一、离散型随机变量的生成 离散型随机变量 X,它的取值是非光滑连续的值，它只能间断地即离散地取值x1,x2,x3，…，xn，且规定x1＜x2＜x3＜…<xn。其概率密度函数为 p(xi)=p{X= xi} 概率分布函数为

譬如，对某车间每天需求某种零件的数量历史数据中统计获得表5．2的结果。 表5．2某零件每天需求量 X 概率P(x) 累积概率F(x) 可分配的随机数范围 X1=10 0.10 F(X1)=0.10 .00 -- .09 X2=20 0.20 F(X2)=0.30 .10 -- .29 X3=30 0.40 F(X3)=0.70 .30 -- .69 X4=40 0.25 F(X4)=0.95 .70 -- .94 X5=50 0.05 F(X5)=1.00 .95 -- .99 随机变量生成的算法为 ①产生一个u(0， 1)，并令i=0； ②令i=i+1； ③若u＞F(xi)，转回到第②步，否则转至④； ④输出得 X＝xi。

1．逆变换法。任何概率分布的随机变量 X，其分布函数 F(x)的值域是(0,l)；因此，可令 二、连续型随机变量的生成 以下介绍常用几种概率分布的随机变量生成方法。 1．逆变换法。任何概率分布的随机变量 X，其分布函数 F(x)的值域是(0,l)；因此，可令 F(x)=u (5．17) 如果对上式能解出显式的逆函数 x=F-1(u) (5．18) 则调用随机数生成算法产生一个随机数u，将其代入显式(8．18)，即可获得随机变量 X的一次取值 x。可惜能求得显式逆函数的随机变量，只有如下几种分布。 ( l)在(a，b)区间均匀分布的随机变量 X。 其概率密度函数为 其分布函数为

均值为(a十 b)／2；方差为(b— a)2／12。 两函数的图形如图8．1。 令u=(x- a)／( b- a)，则可解得 x = a + (b – a).u (5．21) (2)三角分布( a， m，b)的随机变量 其密度函数为 其分布函数为

k U P 结果 0.443 P>0.0025 1 0.172 0.076 2 0.423 0.032 3 0.819 0.026 4 0.255 0.0066 5 0.749 0.0049 0.225 0.0011 P<0.0025,置X=6 k U P 结果 0.443 P>0.0025 1 0.172 0.076 2 0.423 0.032 3 0.819 0.026 4 0.255 0.0066 5 0.749 0.0049 0.225 0.0011 P<0.0025,置X=6 3．近似计算法。正态分布随机变量也是一种常见的随机变量。由于其分布函数无法积分得到显式，故不能应用逆变换法。可采用如下近似法，也可采用函数变换法。

若要产生均值为u，方差为2的正态分布随机变量 X，则 Z是均值u=0，方差为2 =1的近似正态分布；随着n值的增大，近似程度则愈大。实践表明，当取n=12时，近似性已能满意，则有 若要产生均值为u，方差为2的正态分布随机变量 X，则 例题：要产生 u = 7．3分钟， 2=11.7分2的正态分布随机变量 x，则 ①先相继产生12个u，比如为 0．390，0．977，0．492，0．323，0．434，0．994 0．886，0．919，0．153，0．522，0．014，0．303 ②按(8．34)式，

5.3随机系统模拟常规操作方法 时间步长法的模拟操作过程如图9．1所示。 从图9．1提示，对时间步长法应注意下述几个要点：  § 5.3.1 时间步长法 时间步长法的模拟操作过程如图9．1所示。 从图9．1提示，对时间步长法应注意下述几个要点： (l)首先应进行被模拟系统的系统分析，明确模拟目的，确定系统状态变量和决策变量，包括随机变量，建立模拟数学模型，搜集和确定有关参数、常数等数据，以及随机变量生成算法。 (2)建立模拟时钟，即确定一个称作“时钟”的变量T，其初值T。一般定为T0=0，以便对模拟过程进行计时。 (3)设定模拟总时间长度，它与时间步长t的确定有关。 (4)确定固定的时间步长t值。时间步长t愈小，则愈能较真实地考察、记录到系统的变化过程，但模拟工作量较大，占用机时较多；时间步长t愈大，则可能会遗漏一些系统的演变环节，造成模拟过程和结果的失真。 (5)确定模拟开始的系统初始状态及其有关参数值。 (6)建立随机数生成算法的子程序，确定其种子数。 (7)设计输出模拟结果的要求和方式。 (8)用一般高级语言编写模拟程序并调试。