应用随机过程第4章随机模拟

第4章随机模拟2016-2017学年第2学期统计与信息学院张建新2017/4/19第4章随机模拟StochasticSimulation›4.1随机模拟›4.2随机数的抽样›4.3随机游走模拟›4.4泊松过程模拟›4.5马尔可夫链蒙特卡罗方法›参考文献4.1随机模拟›随机模拟是运用计算机对随机系统进行的一种仿真研究方法。因其简单实用、适用面广，已经成为科学技术各领域的有力研究手段。›随机模拟方法也称为“蒙特卡罗（MonteCarlo）”方法，其灵魂在于由计算机生成随机数序列，从而使人们可以由此模拟出各种随机现象（随机事件）。›注：蒙特卡罗是摩纳哥国的世界闻名赌城›随机模拟方法要求研究者具有概率论与随机过程的基本知识以及最基本的编程能力。4.1随机模拟›蒙特卡罗（MonteCarlo）方法定义：theuseofrandomsamplingtechniquesandoftentheuseofcomputersimulationtoobtainapproximatesolutionstomathematicalorphysicalproblemsespeciallyintermsofarangeofvalueseachofwhichhasacalculatedprobabilityofbeingthesolution(Merriam-Webster,Inc.,1994,pp.754-755)›蒙特卡罗模拟的基本思想：当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。4.1随机模拟›蒙特卡罗模拟的试验过程：计算机模拟试验过程，就是将试验过程转化为数学问题，在计算机上实现。在解决实际问题的时候应用蒙特卡罗模拟主要有两部分工作：用蒙特卡罗模拟某一过程时，需要产生各种概率分布的随机变量，然后用统计方法把模型的数字特征估计出来，从而得到实际问题的数值解。›蒙特卡罗模拟主要适用于哪些情况？当我们手上的数据并不能满足统计理论的假设时，蒙特卡罗模拟就是理论预测的一个很好的替代。蒙特卡罗模拟的一个主要应用是评估估计结果是否违反最初的假设，另一个应用是当一个样本没有理论分布时，我们可以用蒙特卡罗模拟来决定样本分布。4.1随机模拟›用蒙特卡罗模拟掷骰子（SAS实现）›利用SAS中的随机数函数与数学函数，模拟掷骰子游戏。每次游戏独立地掷骰子2次，分别用变量N1,N2记录2次抛掷出现的点数，用变量Sum记录2次抛掷出现的点数之和，用变量Id记录游戏的序号。4.1随机模拟›用蒙特卡罗模拟掷骰子（SAS实现）/*建立数据集Dice*/DATADICE;DOId=1TO10000;***掷10000次骰子;N1=1+INT(6*RANUNI(123));***第一次掷骰子的结果;N2=1+INT(6*RANUNI(123));***第二次掷骰子的结果;SUM=N1+N2;***两次掷骰子的结果的和;OUTPUT;***保存两次掷骰子的结果及两次掷骰子的结果的和;END;RUN;4.1随机模拟›用蒙特卡罗模拟掷骰子（SAS实现）这里我们可以看到两次投出的点数之和sum=7的Percent在1/6=16.67附近，也即验证了概率论计算的结果。4.1随机模拟›用蒙特卡罗模拟掷骰子（SAS实现）›变量Sum的理论分布律Sum23456789101112P1/362/363/364/365/366/365/364/363/362/361/36P0.0280.0560.0830.1110.1390.1670.1390.1110.0830.0560.0284.1随机模拟›在随机过程课程中学习随机模拟的意义–随机过程是典型的随机系统，表现出复杂的不确定性，不容易被理解和掌握。通过计算机仿真，帮助同学们进一步认识与感知随机过程。–利用“随机模拟”方法研究随机过程。–除理论研究、科学实验之外，在认识论层面拓宽同学们的思路，补充“随机模拟”这种研究方法。搭建“随机过程”课程理论与实际应用之间的桥梁。4.2随机数的抽样›标准分布随机数生成–伪随机数生成算法›非标准分布随机数生成–逆变换法–接受拒绝法›多维联合分布抽样4.2随机数的抽样›标准分布随机数生成–伪随机数生成算法–蒙特卡罗模拟的核心是随机数据的生成器，而各种随机数据的生成器的核心是随机均匀分布数据的生成器，因为各种特定分布的随机数据都可以通过随机均匀分布数据得到。–产生随机数，可以通过物理方法取得，但当今最为普遍的乃是在计算机上利用数学方法产生随机数。这种随机数根据特定的迭代公式计算出来，初值确定后，序列就可以预测出来，所以不能算是真正的随机数（就成为“伪随机数”）。4.2随机数的抽样›标准分布随机数生成–伪随机数生成算法–RANUNI生成的随机均匀分布数据有着广泛的应用，它主要是通过以下递归公式来产生随机数据，–R(i+1)=(a*R(i)+c)(modm)i=0,1,2…，–这里R(i)为第i个随机数，a为乘子，c为增量，(modm)是指将(a*R(i)+c)得到的数据取模运算，使得得到的数据都在(0,m)这个区间内。–要产生的随机数序列R(i)从第一个随机数R(0)开始和决定，然后得到在(0,m)这个区间的均匀分布的随机数，SAS会通过除以m得到(0,1)区间的随机数。4.2随机数的抽样›标准分布随机数生成–伪随机数生成算法–在SAS系统中，–常数a=397,204,094–m=2^31-1=2,147,483,647（是一个素数）–c=0–种子R(0)必须是一个整数并且其值介于1到m-1之间。–这里c=0的数据生成器被称为multiplicativecongruentialgenerator，被广泛地应用。4.2随机数的抽样›标准分布随机数生成–SAS随机数函数4.2随机数的抽样›标准分布随机数生成›利用SAS生成标准分布随机数一般有两种方法–由随机数函数产生随机数序列其语法为：var=name(seed,arg)–CALL子程序产生随机数序列其语法为：callname(seed,arg,var)。ｰ两种方法的主要区别在于：ｰ随机数函数产生随机数序列时，其序列的值只由第一个随机数种子的值决定，而用CALL子程序时，每一次调用随机函数，都会重新产生新的随机数种子。4.2随机数的抽样›标准分布随机数生成，利用SAS生成标准分布随机数›Seed为零时的随机数序列：–随机数据集由随机数种子(seed)决定，它的值可以是（-（2^31-2），2^31-2）中的任何一个数。其中，当随机数种子(seed)的值大于0时，每次产生的R(0)的值都是一样的，但当的seed值小于或等于零时，每次产生的R(0)的值是不同的，它会以系统的时间为种子产生随机数值。–当Seed为零时，每一次执行CALL子程序所产生的随机数序列结果都会不一致，我们可以将SEED=0看作是随机种子。4.2随机数的抽样›标准分布随机数生成，利用SAS生成标准分布随机数›生成大量不重复的seed序列–在实际的应用中，我们经常会遇到需要大量随机数序列的情况，这时候我们就不能靠手工输入随机数种子。–当SEED=0时，我们可以用这个随机种子产生大量的随机数序列，但是这里产生的随机数序列并不一定能保证这些随机数序列不重复。–这里介绍一个产生不重复的随机数种子的宏4.2随机数的抽样›生成大量不重复的seed序列产生随机数种子的原理，是要产生多少个随机数种子，就按一定步长递增多少次，然后得到一个随机数作为种子。这个宏有个缺点，就是当步长*随机数种子数量2**31-1时，可能得不到要求得到的随机数种子数量。4.2随机数的抽样›由SAS生成50000个服从正态分布的随机数，并对生成的样本作描述性统计分析和拟合分布检验。4.2随机数的抽样›由SAS生成50000个服从正态分布的随机数，并对生成的样本作描述性统计分析和拟合分布检验。分布拟合检验结果：模拟样本来自标准正态总体。4.2随机数的抽样›由SAS生成50000个服从正态分布的随机数，并对生成的样本作描述性统计分析和拟合分布检验。模拟样本数据直方图对标准正态分布的密度函数拟合良好。4.2随机数的抽样›由SAS生成1000个服从指定离散分布（1/0.12/0.23/0.44/0.25/0.1）的随机数，对生成的样本作频数统计，并作分布拟合检验。4.2随机数的抽样›由SAS生成1000个服从指定离散分布（1/0.12/0.23/0.44/0.25/0.1）的随机数，对生成的样本作频数统计，并作分布拟合检验。样本容量为1000时，样本分布与期望分布略有差异。分布拟合检验结果接受原假设，即认为样本分布与期望分布相同。4.2随机数的抽样›由SAS生成50000个服从指定离散分布（1/0.12/0.23/0.44/0.25/0.1）的随机数，对生成的样本作频数统计，并作分布拟合检验。样本容量为50000时，样本分布与期望分布的差异减小，（检验统计量的值由7.0625减小为2.1371）。分布拟合检验结果接受原假设，即认为样本分布与期望分布相同。4.2随机数的抽样›非标准分布随机数生成：逆变换法4.2随机数的抽样›非标准分布随机数生成：接受-拒绝法›接受拒绝法（AcceptanceRejectionMethod）：假设希望生成的随机数的概率密度函数（PDF）为f，则首先找到一个PDF为g的随机数发生器与常数M，使得f伪随机数发生器f（x）≤Mg（x），然后根据接收拒绝算法求解。›由于算法平均运算M次才能得到一个希望生成的随机数，因此c的取值必须尽可能小。显然，该算法的缺点是较难确定g与M。4.2随机数的抽样›非标准分布随机数生成：接受-拒绝法›接受拒绝法（AcceptanceRejectionMethod）算法：›1.选择一个容易抽样的分布g作为建议分布，确定常数M1使得在f的定义域上f（x）≤Mg（x）成立；›2.生成服从分布g的建议随机数y；›3.生成一个服从均匀分布U(0,1)的随机数u；›4.计算接受比函数h(x)=f(x)/(Mg(x)),如果uh(x)，则接受建议随机数y，否则拒绝y，返回2。›由此生成的随机数序列即为服从目标分布f的随机数。证明略。4.2随机数的抽样›非标准分布随机数生成：接受-拒绝法举例›目标分布（概率密度函数）f(x)=6(x-0.5)2/7，x∈[0,2]›建议分布为均匀分布U(0,2),由f在区间[0,2]上的最大值1.929得到M=3.858。›SAS程序：***Histogram***;goptionsftext=SWISSctext=BLACKhtext=1cells;symbolv=SQUAREc=BLUEh=1cells;procunivariatedata=arm_expnoprint;varX;histogram/caxes=BLACKcframe=CXF7E1C2waxis=1cbarline=BLACKcfill=BLUEpfill=SOLIDvscale=percenthminor=0vminor=0name='HIST';run;dataARM_exp;N=50000;x=0;doi=1toN;y=2*ranuni(0);u=ranuni(0);h=6*((y-0.5)**2)/7/1.929;ifuhthenx=y;output;end;run;4.2随机数的抽样接受-拒绝法举例›目标分布（概率密度函数）f(x)=6(x-0.5)2/7,x∈[0,2]›SAS生成随机数的直方图4.3随机游走模拟›随机游走模拟›随机游走应用--赌徒破产4.3随机游走模拟›随机游走模拟x-100102030405060708090100110120130140150160170180190n010020030040050060070080090010004.3随机游走模拟›随机游走模拟datasimulation_randomwalk_2