1、第四章第四章 离散事件系统仿真基础离散事件系统仿真基础主要内容主要内容n基本概念基本概念n随机变量模型的确定随机变量模型的确定n随机数的产生随机数的产生n随机变量的产生随机变量的产生第一节第一节 基本概念基本概念n离散事件系统离散事件系统n状态仅在离散时间点上变化，且离散时间点状态仅在离散时间点上变化，且离散时间点一般不确定一般不确定n面向事件；反映系统各部分相互作用的一些面向事件；反映系统各部分相互作用的一些事件，模型为反映事件状态的数集，仿真结事件，模型为反映事件状态的数集，仿真结果是产生处理这些事件的时间历程果是产生处理这些事件的时间历程n连续系统：时间常为均匀间隔计时；系统动连续系统：
2、时间常为均匀间隔计时；系统动力学模型由表征系统变量间关系的方程描写，力学模型由表征系统变量间关系的方程描写，结果常为变量随时间的变化历程结果常为变量随时间的变化历程例：单服务台排队系统例：单服务台排队系统n系统工作时间长度固定系统工作时间长度固定n顾客到达时间随机顾客到达时间随机n服务员服务时间随机服务员服务时间随机n要求通过仿真估计系统要求通过仿真估计系统工作情况，以决定是否工作情况，以决定是否增加服务台增加服务台n显然，离散事件系统一般有固有的随机显然，离散事件系统一般有固有的随机性（注意：连续系统也有随机性，如白性（注意：连续系统也有随机性，如白噪声，但两者行为不同）噪声，但两者行为不同
3、）n研究的理论基础：经典的概率及数理统研究的理论基础：经典的概率及数理统计理论、随机过程理论计理论、随机过程理论n简单系统可能有理论解析解，但对实际简单系统可能有理论解析解，但对实际系统，只有靠计算机仿真计算才有可能系统，只有靠计算机仿真计算才有可能提供较完整的结果提供较完整的结果常用概念常用概念n实体实体n永久实体：永久驻留在系统中，是系统处于永久实体：永久驻留在系统中，是系统处于活动的必要条件，如服务员活动的必要条件，如服务员n临时实体：仅在系统中存在一段时间，按一临时实体：仅在系统中存在一段时间，按一定规律到达，如顾客定规律到达，如顾客n关系：临时实体按一定规律不断产生，在永关系：临时实
4、体按一定规律不断产生，在永久实体作用下通过系统，最后离开系统久实体作用下通过系统，最后离开系统n事件事件n引起系统状态发生变化的行为引起系统状态发生变化的行为n离散事件系统本质是由事件驱动的离散事件系统本质是由事件驱动的n例：顾客到达事件使服务员状态由闲到忙，例：顾客到达事件使服务员状态由闲到忙，或使队列长度加或使队列长度加1n事件的发生一般与某一类实体相联系，放在事件的发生一般与某一类实体相联系，放在事件表中管理，事件表通常记录事件类型、事件表中管理，事件表通常记录事件类型、发生条件、时间及相关实体的有关属性发生条件、时间及相关实体的有关属性n活动活动n导致系统状态变化的一个过程为活动导致系
5、统状态变化的一个过程为活动n活动表示两个可区分事件之间的过程，标志活动表示两个可区分事件之间的过程，标志着系统状态的转移着系统状态的转移n如顾客到达事件与顾客开始接受服务事件之如顾客到达事件与顾客开始接受服务事件之间为一活动，使服务员忙及队列长度减间为一活动，使服务员忙及队列长度减1n进程进程n相当于系统的子集或子系统，包含若干个事相当于系统的子集或子系统，包含若干个事件及活动，并且描述了其所包含事件及活动件及活动，并且描述了其所包含事件及活动间的逻辑关系和时序关系间的逻辑关系和时序关系n如某一顾客在系统中的全部活动为一进程如某一顾客在系统中的全部活动为一进程事件、活动、进程的关系图事件、活动
6、、进程的关系图n仿真时钟仿真时钟n离散事件动态系统的状态本来就只在离散时间点上离散事件动态系统的状态本来就只在离散时间点上发生变化，因而不需要进行离散化处理。发生变化，因而不需要进行离散化处理。n离散事件系统一般不以时间推动，但事件间有时序离散事件系统一般不以时间推动，但事件间有时序关系，仿真中仍必须有控制时间的部件关系，仿真中仍必须有控制时间的部件n由于引起状态变化的事件发生时间的随机性，仿真由于引起状态变化的事件发生时间的随机性，仿真钟的推进步长则完全是随机的钟的推进步长则完全是随机的n两个相邻发生的事件之间系统状态不会发生任何变两个相邻发生的事件之间系统状态不会发生任何变化，因而仿真钟可
7、以跨过这些化，因而仿真钟可以跨过这些“不活动不活动”周期周期,n仿真钟的推进呈现跳跃性，推进速度具有随机性。仿真钟的推进呈现跳跃性，推进速度具有随机性。n统计计数器统计计数器n因固有的随机性，某一次仿真运行得到的状因固有的随机性，某一次仿真运行得到的状态变化过程只不过是随机过程的一次取样，态变化过程只不过是随机过程的一次取样，离散事件系统的仿真结果只有在统计意义下离散事件系统的仿真结果只有在统计意义下才有参考价值才有参考价值n在仿真模型中在仿真模型中,需要有一个统计计数部件，需要有一个统计计数部件，以便统计系统中的有关变量，如排队系统中以便统计系统中的有关变量，如排队系统中的顾客等待时间、队列
8、长度等的顾客等待时间、队列长度等离散事件系统仿真离散事件系统仿真的一般步骤的一般步骤n系统建模：系统建模：n一般用流程图描述，一般用流程图描述，反映临时实体在系统反映临时实体在系统内部历经的过程、永内部历经的过程、永久实体对临时实体的久实体对临时实体的作用及相互间逻辑关作用及相互间逻辑关系系n关键：确定随机变量关键：确定随机变量的模型的模型n确定仿真算法确定仿真算法n产生随机变量产生随机变量n确定仿真建模策略确定仿真建模策略n事件调度法：面向事件建立仿真模型事件调度法：面向事件建立仿真模型n活动扫描法：面向活动建模活动扫描法：面向活动建模n进程交互法：面向进程建模进程交互法：面向进程建模n三阶
9、段法：结合活动扫描与事件调度三阶段法：结合活动扫描与事件调度n图形仿真方法：图形仿真方法：Petri网网n建立仿真模型建立仿真模型n定义状态变量、定义系统事件及有关属性、定义状态变量、定义系统事件及有关属性、活动及进程、设计仿真钟的推进方法等活动及进程、设计仿真钟的推进方法等n仿真程序设计及运行仿真程序设计及运行n仿真语言或高级语言仿真语言或高级语言n长期运行或多次运行长期运行或多次运行n仿真结果分析仿真结果分析n统计结果、可信度分析等统计结果、可信度分析等第二节第二节 随机变量模型的确定随机变量模型的确定n无序中蕴含着有序，随机过程也有数学无序中蕴含着有序，随机过程也有数学描述形式，可近似归
10、纳总结为几种变量描述形式，可近似归纳总结为几种变量分布模式，使定量研究成为可能分布模式，使定量研究成为可能n没有绝对的无序和有序，如没有绝对的无序和有序，如混沌混沌n以单服务台排队系统中顾客到达时刻为以单服务台排队系统中顾客到达时刻为例，总可以找到一种接近的随机变量分例，总可以找到一种接近的随机变量分布布n通常需要从观测数据中寻找规律通常需要从观测数据中寻找规律n在寻找分布形式时，根据对随机变量在寻找分布形式时，根据对随机变量（Random variable,r.v.）的特性了解程）的特性了解程度，一般会遇到三种情况度，一般会遇到三种情况nr.v.分布类型已知，需要由观测数据确定分分布类型已知
11、，需要由观测数据确定分布参数布参数n需要由观测数据确定概率分布类型及参数需要由观测数据确定概率分布类型及参数n难以由观测数据确定理论分布形式，需要定难以由观测数据确定理论分布形式，需要定义实验分布义实验分布一、分布参数的确定一、分布参数的确定n分布参数的类型分布参数的类型 定义分布所采用的大多数参数，由物理定义分布所采用的大多数参数，由物理或几何解释，可分为三个基本类型或几何解释，可分为三个基本类型n位置参数位置参数n比例参数比例参数n形状参数形状参数n位置参数位置参数n确定了一个分布函数确定了一个分布函数取值范围的横坐标取值范围的横坐标n当改变时，分布函数当改变时，分布函数仅平移而无其它变化
12、，仅平移而无其它变化，又称位移参数又称位移参数n例均匀分布函数例均匀分布函数U(a,bU(a,b)，密度函数，密度函数1,()0,axbf xba其它其中其中a，b均可定均可定义为位置参数义为位置参数n比例参数比例参数n决定分布函数在其取决定分布函数在其取值范围内取值的比例值范围内取值的比例尺尺n的改变只压缩或扩张的改变只压缩或扩张分布函数，不改变基分布函数，不改变基本形状本形状n例：指数分布函数例：指数分布函数EXPO()，密度函，密度函数数1,0()0,0 xexf xxn形状参数形状参数n确定分布函数的形状，确定分布函数的形状，从而改变分布函数的从而改变分布函数的性质性质n例：韦伯分布例
13、：韦伯分布Weibull(,)，密，密度函数：度函数：()1,0()0,0 xxexf xx分布参数的估计分布参数的估计n常用方法常用方法n最大似然估计最大似然估计(maximum likelihood estimation)n在已经得到试验结果的情况下，应该寻找使这个在已经得到试验结果的情况下，应该寻找使这个结果出现的可能性最大的那个结果出现的可能性最大的那个参数值参数值作为真值的作为真值的估计估计n最小二乘估计最小二乘估计(least-square estimation)n无偏估计无偏估计(unbiased estimation)最大似然估计法估计分布参数最大似然估计法估计分布参数讨论一个
14、未知参数讨论一个未知参数的情形，设观测数据为的情形，设观测数据为n离散分布情形：可令离散分布情形：可令 为该分布的概率质量为该分布的概率质量函数，定义似然函数函数，定义似然函数L(L()为：为：的最大似然估计值的最大似然估计值 使使L(L()取最大值取最大值n连续分布情形：令连续分布情形：令 为概率密度函数，定义为概率密度函数，定义似然函数为似然函数为12,nx xx()P x12()()()()nLP x P xP x()fx12()()()()nLfxfxfx连续分布例：指数分布连续分布例：指数分布/1()xfxe(0)12/11111()()()()exp()nnxxxniiLeeex(
15、)Lln()L11()ln()lnniiRLnx 21110/()nniiiidRnxxnx nd 222312niid Rnxd()x nix220d Rd密度函数为密度函数为，被估计的参数，被估计的参数由由为求使为求使取最大值的取最大值的，取对数，取对数因因严格递增，严格递增，取最大值等价于取最大值等价于最大最大求极值求极值考虑到考虑到当当时时，由于，由于为正为正，显然，显然为最大值为最大值()Lln()L离散分布例：泊松分布离散分布例：泊松分布被估计参数被估计参数质量函数质量函数由由令令故为最大似然估计值故为最大似然估计值(0)()/!(0,1,2,)xxP xexx1121()()()
16、()/(!)niinxnniiLP x P xP xex11()ln()()lnln(!)nniiiiRLnxx 1110/()nniiiidRnxxnx nd 222110niid Rxd n对多个参数的估计问题，上述方法原则上适用，对多个参数的估计问题，上述方法原则上适用，只是计算更繁琐，有时不能解析计算，必须用只是计算更繁琐，有时不能解析计算，必须用数值计算方法数值计算方法n常见分布的参数的最大似然估计常见分布的参数的最大似然估计n均匀分布均匀分布 密度函数密度函数 参数参数(,)U a b1()0axbf xbaother1minii nax 1maxii nbx n正态分布正态分布
17、密度函数密度函数 位置参数位置参数 比例参数比例参数n二项分布二项分布 质量函数质量函数t为正整数，为正整数，t已知时已知时2(,)N 22()/221()2xf xe()x n21/21()nSnn(,)bin t p(1),0,1,()0,xt xtppxtP xxother (0,1)p()/px nt二、分布类型的假设二、分布类型的假设n对观测数据进行预处理后，才能对分布对观测数据进行预处理后，才能对分布类型进行假设类型进行假设n连续分布数据的预处理方法连续分布数据的预处理方法n点统计法、直方图法、概率图法点统计法、直方图法、概率图法n离散分布数据的预处理方法离散分布数据的预处理方法n
18、点统计法、线图法点统计法、线图法n定义实验分布定义实验分布连续分布类型的假设（连续分布类型的假设（1）n点统计法：基于连续分布的偏差系数特点统计法：基于连续分布的偏差系数特征进行分布类型假设征进行分布类型假设n偏差系数偏差系数 （方差、均值）（方差、均值）n常见分布的偏差系数取值范围常见分布的偏差系数取值范围n指数分布指数分布n均匀分布均匀分布 ，n取值取值 ，除，除0以外以外n正态分布正态分布 ，取值同上，取值同上()/()Var xE x()EXPO1(,)U a b()/3()baab(,)2(,)N/n预处理预处理 ，n则则n显然不能唯一确定分布类型，且计算不一定显然不能唯一确定分布类
19、型，且计算不一定无偏，但对确定理论分布具有指导作用无偏，但对确定理论分布具有指导作用1()()niiE xx nx221()()()/(1)niiVar xSnxx nn2()/()Snx n连续分布类型的假设（连续分布类型的假设（2）n直方图法直方图法n将将 的取值范围等分为的取值范围等分为k个区间个区间 ，n令令 表示第表示第j个区间上观测数据个数个区间上观测数据个数 占总数占总数的比例，即的比例，即n定义函数定义函数1nxx01121,),),)kkb bb bbb1,(1,2,)jjbbbjk jgjn/,(1,2,)jjgnnjk010()0ijjkxbh xgbxbxbn作出作出h
20、(x)的直方图，并与理论分布的密度函的直方图，并与理论分布的密度函数图形比较，选择接近的分布数图形比较，选择接近的分布n困难：区间长度困难：区间长度b的取值的取值n太大：丢失很多信息太大：丢失很多信息n太小：不能抑制观测噪声太小：不能抑制观测噪声n一般应多选几个一般应多选几个b，取使直方图包络线较光滑，取使直方图包络线较光滑的的区间长度区间长度连续分布类型的假设（连续分布类型的假设（3）n概率图法概率图法n通过比较分布函数确定理论分布通过比较分布函数确定理论分布n设设 有有m个取值，记为个取值，记为x(1)x(m)n定义实验分布定义实验分布 其中其中 表示小于或等于表示小于或等于x(i)的观测
21、数据的个的观测数据的个数，显然数，显然n为避免由有限个观测数据得到的实验分布函为避免由有限个观测数据得到的实验分布函数值数值1，修正为，修正为1nxx()/,(1,2,)iF x inn iminmnn()(0.5)/iF x innn概率图不能直接比较（基本都是概率图不能直接比较（基本都是S型），而型），而是采用是采用“分位点分位点”比较法比较法n定义：设定义：设0g0时称混合乘同余法，时称混合乘同余法，c=0称乘同称乘同余法余法素数取模乘同余法（素数取模乘同余法（PMMLCG）nm是小于是小于2b的最大素数，而的最大素数，而a的选择满足的选择满足al-1被被m整除的最小整数整除的最小整数l
22、=m-1，也就是说，也就是说能被能被m整除的整除的al-1的最小整数为的最小整数为am-1-1，那，那么得到的么得到的zj的周期为的周期为m-1，且在每个周期，且在每个周期内，内，1,2,m-1这些整数严格地只出现一这些整数严格地只出现一次。次。两个经过检验性能较好的两个经过检验性能较好的PMMLCG ZZii551(mod 235 31)ZZii851(mod 231 1)二、组合发生器二、组合发生器n为提高性能，可用一个发生器控制另一个发生为提高性能，可用一个发生器控制另一个发生器产生随机数，两种方式器产生随机数，两种方式n发生器发生器1产生产生 ，发生器，发生器2产生产生1,k上上均匀分
23、布的随机整数均匀分布的随机整数I，取出并重新产生，取出并重新产生UIn发生器发生器1，2产生产生 和和 ，在二进制下循在二进制下循环移位环移位 次得次得 ，再次与，再次与 “异或异或”后得后得n思路：思路：n减少递推公式的自相关性，提高独立性减少递推公式的自相关性，提高独立性n加长发生器周期，提高随机数密度和均匀性加长发生器周期，提高随机数密度和均匀性1(,)kUU(1)iZ(2)iZ(2)iZ(1)iZ(2)iZ(2)iZiZ三、随机数发生器的测试三、随机数发生器的测试n1.均匀性检验均匀性检验n常用频率检验：将随机数发生器取值范围常用频率检验：将随机数发生器取值范围0,1分为分为k个等长区
24、间，由该发生器产生个等长区间，由该发生器产生N个随机数，则落在每个子区间上随机数个数个随机数，则落在每个子区间上随机数个数理论值为理论值为n=N/k，称为理论频率，实际第，称为理论频率，实际第j个个区间上的数据个数区间上的数据个数nj 总会有偏差，由此可总会有偏差，由此可得得 ，其大小反映均匀性程度，其大小反映均匀性程度221()kjjnnnn2.独立性检验独立性检验n计算该随机数序列相邻一定间隔的随机数之计算该随机数序列相邻一定间隔的随机数之间的相关系数，判断相关程度间的相关系数，判断相关程度n先由随机数发生器产生先由随机数发生器产生N个随机数个随机数Ui，则前，则前后相隔为后相隔为j个数的
25、相关系数的均值为个数的相关系数的均值为 其中其中S2为随机数方差的估计值：为随机数方差的估计值：再根据统计理论处理再根据统计理论处理22111()/2njjiijiUUSNj22111()12NiiSUNC语言中与随机数有关的函数nrandomize(void);初始化nx=random(int M);产生0M 之间的随机数nx=rand(void);产生0215-1之间的随机数Matlab中与随机数有关的函数nRand(n)：n个0,1之间均匀分布随机数nRand(m,n)：m*n 个0,1之间均匀分布随机数nrandn(n):n个N(0,1)标准正态分布随机数nrandn(m,n)：m*n
26、个N(0,1)标准正态分布随机数例例 射击命中率射击命中率 在我方某前沿防守地域，敌人以一个炮排在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点经过长期观察发现，装并经常变换射击地点经过长期观察发现，我方指挥所对敌方目标的指示有我方指挥所对敌方目标的指示有50是准确的，是准确的，而我方火力单位，在指示正确时，有而我方火力单位，在指示正确时，有1/3的射的射击效果能毁伤敌人一门火炮，有击效果能毁伤敌人一门火炮，有1/6的射击效的射击效果能
27、全部消灭敌人果能全部消灭敌人 现在希望能用某种方式把我方将要对敌人现在希望能用某种方式把我方将要对敌人实施的实施的20次打击结果显现出来，确定有效射击次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值的比率及毁伤敌方火炮的平均值说明说明n这是一个概率问题，只涉及到简单的随机数的这是一个概率问题，只涉及到简单的随机数的产生和处理，可以通过理论计算得到相应的概产生和处理，可以通过理论计算得到相应的概率和期望值。但这样只能给出作战行动的最终率和期望值。但这样只能给出作战行动的最终静态结果，而显示不出作战行动的动态过程静态结果，而显示不出作战行动的动态过程n为了能显示我方为了能显示我方20次
28、射击的过程，考虑采用计次射击的过程，考虑采用计算机模拟算机模拟1.问题分析问题分析需要模拟出以下两件事：需要模拟出以下两件事：n(1)观察所对目标的指示正确与否观察所对目标的指示正确与否n模拟试验有两种结果，每一种结果出现的概率都是模拟试验有两种结果，每一种结果出现的概率都是1/2n因此，可以用因此，可以用投掷一枚硬币投掷一枚硬币的方式予以确定，当硬的方式予以确定，当硬币出现正面时为指示正确，反之为不正确币出现正面时为指示正确，反之为不正确n(2)当指示正确时，我方火力的射击结果当指示正确时，我方火力的射击结果模拟试验有三种结果：模拟试验有三种结果：n毁伤一门火炮的可能性为毁伤一门火炮的可能性
29、为1/3(即即2/6)，n毁伤两门的可能性为毁伤两门的可能性为1/6，n没能毁伤敌火炮的可能性为没能毁伤敌火炮的可能性为1/2(即即3/6)这时这时可用投掷骰子的方法来确定：可用投掷骰子的方法来确定：n出现、点：则认为没击中敌人；出现、点：则认为没击中敌人；n出现、点：出现、点：则认为击毁敌人一门火炮；则认为击毁敌人一门火炮；n出现点：出现点：则认为击毁敌人两门火炮则认为击毁敌人两门火炮2.符号假设符号假设i：要模拟的打击次数；：要模拟的打击次数；k1：没击中敌人火炮的射击总数；：没击中敌人火炮的射击总数；k2：击中敌人一门火炮的射击总数；：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的
30、射击总数：击中敌人两门火炮的射击总数E：有效射击比率；：有效射击比率；E1：20次射击平均每次毁伤敌人的火炮数次射击平均每次毁伤敌人的火炮数3.程序流图程序流图初始化初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子骰子点数点数?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i20?E=(k2+k3)/20 E1=0*k1/20+1*k2/20+2*k3/20停止停止硬币正面硬币正面?YNNY1，2，34，564.模拟结果模拟结果消灭敌人火炮数消灭敌人火炮数 试验试验 序号序号 投硬币投硬币 结结 果果 指示指示 正确正确 指指 示示 不正确不正确 掷骰子掷骰子 结结 果果
31、正 正 反 正 正 反 正 正 反 反 消灭敌人火炮数消灭敌人火炮数 试验试验 序号序号 投硬币投硬币 结结 果果 指示指示 正确正确 指指 示示 不正确不正确 掷骰子掷骰子 结结 果果 正 反 正 反 正 正 正 正 反 正 n从以上模拟结果可以计算出：从以上模拟结果可以计算出：有效射击有效射击 E=7/20=0.35平均值平均值113430120.5202020E 5.理论计算理论计算则由全概率公式：则由全概率公式：000()(0)(|0)(1)(|1)11100.25222EP AP jP AjP jP Aj 111120.33212E 111()(0)(|0)(1)(|1)111102
32、236P AP jP AjP jP Aj 222()(0)(|0)(1)(|1)1111022612P AP jP AjP jP Aj 6.结果比较结果比较n模拟结果与理论计算虽然结果不完全一致，但仍反映模拟结果与理论计算虽然结果不完全一致，但仍反映了事件发生的随机性。只要多次实验求平均值，模拟了事件发生的随机性。只要多次实验求平均值，模拟值就会很接近理论值值就会很接近理论值 理论计算和模拟结果的比较 分类 项目 无效射击 有效射击 平均值 模 拟 理 论 例例：赶火车的仿真：赶火车的仿真n一列火车从一列火车从A站经过站经过B站开往站开往C站，某人每天赶站，某人每天赶往往B站乘这趟火车。已知火
33、车从站乘这趟火车。已知火车从A站到站到B站运行站运行时间为均值时间为均值30分钟、标准差为分钟、标准差为2分钟的正态随分钟的正态随机变量。火车大约在下午机变量。火车大约在下午1点离开点离开A站，离开站，离开时刻的频率分布为时刻的频率分布为出发时刻（出发时刻（T）1：00 1：051：10 频率频率 0.7 0.2 0.1n此人到达此人到达B站时的频率分布为：站时的频率分布为：模拟火车开出、火车到达模拟火车开出、火车到达B站、此人到达站、此人到达B站的情况，站的情况，并给出能赶上火车的仿真结果。并给出能赶上火车的仿真结果。n引入以下变量：引入以下变量：nT1 火车从火车从A站开出的时刻；站开出的
34、时刻；nT2 火车从火车从A站运行到站运行到B站所需要的时间；站所需要的时间；nT3 此人到达此人到达B站的时刻站的时刻；0.10.1 0.20.2 0.40.4 0.3 0.3 频率频率1 1：3434 1 1:32:32 1 1:30:30 1 1:28:28到达时刻（到达时刻（T T）nT1,T2,T3 是随机变量，其概率分布为是随机变量，其概率分布为 x=0.7,x2=0.9,y1=0.3,y2=0.7,y3=0.9 0.10.1 0.20.2 0.70.7 P P（频率）（频率）10 10 5 5 0 0 T1(T1(分分)T3（分）（分）28 30 32 34 P(频率频率)0.3
35、 0.4 0.2 0.1开车时间的仿真测试开车时间的仿真测试 s1=0;s2=0;s3=0;x=rand(10000,1);for i=1:10000 if x(i)0.9 s3=s3+1;end ends1/10000,1-s1/10000-s3/10000,s3/10000 s1=0;s2=0;s3=0;s4=0;x=rand(10000,1);for i=1:10000 if x(i)0.3 s1=s1+1;elseif x(i)0.7 s2=s2+1;else if x(i)0.9 s3=s3+1;else s4=s4+1;end endends1/10000,s2/10000,s3/
36、10000,s4/10000人到达时刻人到达时刻仿真测试仿真测试火车运行时间的仿真测试火车运行时间的仿真测试x=randn(10000,1);for i=1:10000 y(i)=30+2*x(i);end赶上火车的仿真结果赶上火车的仿真结果s=0;x1=rand(10000,1);x2=rand(10000,1);x3=randn(10000,1);for i=1:10000if x1(i)0.7 T1=0;elseif x1(i)0.9 T1=5;else T1=10;endT2=30+2*x3(i);if x2(i)0.3 T3=28;elseif x2(i)0.7 T3=30;else
37、 if x2(i)0.9 T3=32;else T3=34;endendif T3T1+T2 s=s+1;endends/10000第四节第四节 随机变量的产生原理随机变量的产生原理n对产生随机变量的方法的性能要求对产生随机变量的方法的性能要求n准确性：即由这种方法产生的随机变量应准准确性：即由这种方法产生的随机变量应准确地具有所要求的分布确地具有所要求的分布 n快速性要求，在离散事件仿真中，一次运行快速性要求，在离散事件仿真中，一次运行往往需要产生几万甚至几十万个随机变量，往往需要产生几万甚至几十万个随机变量，产生随机变量的速度将极大地影响着仿真执产生随机变量的速度将极大地影响着仿真执行的效
38、率。行的效率。n产生随机变量的常用方法产生随机变量的常用方法n反变换法、组合法、卷积法及取舍法反变换法、组合法、卷积法及取舍法 一、反变换法一、反变换法n常用且直观，以概率积分变换定理为基础常用且直观，以概率积分变换定理为基础n设随机变量设随机变量x的分布函数为的分布函数为F(x)，则，则x的抽样的抽样值可按如下方法求：值可按如下方法求：产生在产生在0,1上的均匀分布的独立随机变量上的均匀分布的独立随机变量u，由反分布函数由反分布函数 得到的值即为所需要的得到的值即为所需要的随机变量随机变量xn该方法对分布函数进行反变换，故此得名该方法对分布函数进行反变换，故此得名1()xFu例：连续随机变量
39、例：连续随机变量n对指数分布对指数分布 产生产生x抽样值抽样值 解：用随机数发生器产生解：用随机数发生器产生1,0()0,0 xexF xx(0,1)uU/()1xuF xe 1()ln(1)xFuu 1(0,1)uU lnxu 令令则则因因(0,1)uU同分布同分布则有则有n由上面的例子可以看到，用反变换法产生随机由上面的例子可以看到，用反变换法产生随机变量时首先必须用随机数发生器产生在变量时首先必须用随机数发生器产生在0,1上均匀分布的独立的上均匀分布的独立的u，以此为基础得到的随，以此为基础得到的随机变量机变量x才能保证分布的正确性，可见，选择才能保证分布的正确性，可见，选择一个均匀性及
40、独立性较好的随机数发生器在产一个均匀性及独立性较好的随机数发生器在产生随机变量中是十分重要的生随机变量中是十分重要的n当当x是离散随机变量时是离散随机变量时，其反变换法的形式略，其反变换法的形式略有不同，原因在于离散随机变量的分布函数也有不同，原因在于离散随机变量的分布函数也是离散的，因而不能直接利用反函数来获得随是离散的，因而不能直接利用反函数来获得随机变量的抽样值机变量的抽样值n对有对有离散随机变量离散随机变量 x，分别以概率，分别以概率 取值取值 ，其中其中 ，由此可划分，由此可划分n个子区间：个子区间：用随机数发生器产生用随机数发生器产生 ，若，若 则取则取x=xi即可即可n显然，区间
41、搜索策略直接影响抽样值产生速度显然，区间搜索策略直接影响抽样值产生速度12(),(),()np xp xp x12,nx xx1()1niip x12nxxx1111211(0,(),(),()(),(),()nniiiip xp xp xp xp xp x100()()iijjjjp xup x(0,1)uU例：离散随机变量例：离散随机变量n质量函数及累积分布函数如下质量函数及累积分布函数如下 用反变换法产生随机变量用反变换法产生随机变量x 先由随机数发生器产生先由随机数发生器产生0,1区间上均匀分布区间上均匀分布的随机变量的随机变量u，设，设u=0.72，按反变换法，逐，按反变换法，逐区间
42、搜索比较，得区间搜索比较，得x=x3=3二、组合法二、组合法n若分布函数可表示成若干其它分布函数之和，若分布函数可表示成若干其它分布函数之和，且这些分布函数更易取样时，选择组合法且这些分布函数更易取样时，选择组合法n设随机变量设随机变量x的分布函数的分布函数 或密度函数或密度函数 其中其中 ，且，且 ，易于取样易于取样 则组合法步骤：则组合法步骤：1()()jjjF xp F x1()()jjjf xp fx0jp 01jjp(),()jjF xfxn(1)产生一个随机整数产生一个随机整数 J，满足，满足 即决定采用哪一个分布函数进行取样即决定采用哪一个分布函数进行取样n(2)产生具有分布函数
43、产生具有分布函数Fi(x)或密度函数或密度函数fi(x)的随机变量，则的随机变量，则x=xin第一步可用离散反变换法实现第一步可用离散反变换法实现 第二步可用反变换法或以下介绍的其它方法第二步可用反变换法或以下介绍的其它方法n显然，两步都需要产生相互独立的显然，两步都需要产生相互独立的,1,2,jP Jjpj(0,1)uU例：设密度函数为例：设密度函数为产生服从该分布的随机变量产生服从该分布的随机变量 x 直接用反变换法产生随机变量直接用反变换法产生随机变量 x很困难，分析很困难，分析f(x)的特点，可知它以纵轴为对称轴分为两部分的特点，可知它以纵轴为对称轴分为两部分f(x)可写成如下形式可写
44、成如下形式|()0.5xf xe(,0)0,)()0.5()0.5()xxf xe I xeI x(,0)1,(,0)()0,xI xother 0,)1,0,)()0,xI xother其中其中n用组合法产生随机变量用组合法产生随机变量x的方法如下：的方法如下：n由由U(0,1)产生随机变量产生随机变量u1及及u2n若若u10.5，则由，则由f1(x)的分布函数产生，可由反变换的分布函数产生，可由反变换法得到法得到n若若 ，则由，则由f2(x)的分布函数产生，即的分布函数产生，即令令从而从而11(,0)()()(0.5)xf xe I xp220,)()()(0.5)xfxeI xp21()
45、()iiif xp f x2lnxu2lnxu 10.5u 三、卷积法三、卷积法n若若 ，为独立同分布的随机为独立同分布的随机变量，则变量，则 x的分布函数与的分布函数与 的分布函数相的分布函数相同，称同，称 x的分布为的分布为Yi分布的分布的m重卷积重卷积n为产生为产生 x，可先独立地产生，可先独立地产生 ，并简单，并简单相加即得相加即得 x，称为卷积法，称为卷积法12mxYYY12,mY YY1miiY12,mY YY例：产生均值为例：产生均值为的的m维维Erlang分布随机变量分布随机变量xnm维维Erlang分布分布En()的概率密度为的概率密度为 难以直接产生难以直接产生xnErla
46、ng(m,)可表示为可表示为m个均值为个均值为 的独立的指的独立的指数随机变量之和，即数随机变量之和，即 以及以及解解:(1)独立产生独立产生m个个U(0,1)随机数随机数Ui (2)用反变换法分别产生用反变换法分别产生 (3)令令 即可即可m1()immximf xe1()(1)immxiF xeln,(1,)iiYuimm 1miixY1()(),(0;0)(1)!mxxf xexm n改进算法改进算法n实际仿真中，为提高算法速度，考虑到对数实际仿真中，为提高算法速度，考虑到对数运算速度较慢，可将运算速度较慢，可将(2)、(3)步改进为步改进为(2)计算计算(3)则则121mimiuu u
47、u1lnmiixum 四、取舍法四、取舍法n反变换法、组合法以及卷积法有一个共同的特反变换法、组合法以及卷积法有一个共同的特点，即直接面向分布函数，因而称为直接法，点，即直接面向分布函数，因而称为直接法，它们以反变换法为基础它们以反变换法为基础n当反变换法难于使用时当反变换法难于使用时(例如随机变量的分布例如随机变量的分布函数不存在封闭形式函数不存在封闭形式)，取舍法是主要方法之，取舍法是主要方法之一，其基本思想是根据概率密度函数一，其基本思想是根据概率密度函数(质量函质量函数数)的几何特征，拟合一个简单的分布形式，的几何特征，拟合一个简单的分布形式，以后者产生抽样值并按照某种判据进行取舍以后
48、者产生抽样值并按照某种判据进行取舍基本思想基本思想n设随机变量设随机变量 x密度函数为密度函数为f(x)，f(x)最大值为最大值为C，x取值范围取值范围0,1n若独立产生两个若独立产生两个0,1区间均匀分布的随机变量区间均匀分布的随机变量u1，u2，则，则Cu1是在是在0,C内均匀分布的随机变内均匀分布的随机变量，且满足量，且满足 的概率为的概率为n若上式成立，则认为若上式成立，则认为 x=u2，否则拒绝，否则拒绝u212()Cuf u1()00121()1f xdxdyP Cuf uCCn显然，抽样成功的概率显然，抽样成功的概率为为f(x)下的面积除以总面下的面积除以总面积积C(见图见图)，
49、直观上看成直观上看成功抽样的点符合所需要功抽样的点符合所需要的分布的分布n一般情形，是根据一般情形，是根据f(x)的的特征规定一个函数特征规定一个函数t(x)，其曲线处于其曲线处于f(x)上方且形上方且形状接近状接近f(x)，再按照判，再按照判据取舍据取舍n对对t(x)的要求是的要求是n n n易于从易于从t(x)进行反变换进行反变换n如果令如果令 ，则，则 从而可将从而可将r(x)看作是一个密度函数，并用看作是一个密度函数，并用r(x)代替代替f(x)取样，以得到所需要的随机变量取样，以得到所需要的随机变量()()t xf x()t x dxC 1()()r xt xC1()()1r x d
50、xt x dxCn由于由于r(x)不是要求的不是要求的f(x)，故此产生选取与舍弃，故此产生选取与舍弃问题，取舍算法如下问题，取舍算法如下n产生产生n由由r(x)独立地产生随机变量独立地产生随机变量u2n检验如下不等式检验如下不等式 若不等式成立，则令若不等式成立，则令x=u2，否则返回第一步，否则返回第一步nr(x)通常可以选择均匀、指数、正态分布等通常可以选择均匀、指数、正态分布等1(0,1)uU122()/()uf ut u五、典型随机变量的产生五、典型随机变量的产生n正态分布正态分布n正态分布正态分布 密度函数密度函数n因分布函数在极坐标下有封闭形式，可采用因分布函数在极坐标下有封闭形
