MATLAB概率论试题与数理统计程序设计MATLAB概率论试题与数理统计程序设计概率论试题部分胡尧概率论试题部分胡尧贵州大学理學院数学系贵州大学理学院数学系EilihdEmail:sciyhugzueducnQQ:贵州师范学院主讲内容主讲内容一、概率基础与随机试验二、MonteCarlo仿真实验三、MonteCarlo积分四、随机数的产生五、建模案例分析一、概率分布与随机试验随机事件与概率基本概念、概率分布与随机试验基本概念：样本点、样本空间、随机事件、概率、樣本点、样本空间、随机事件、概率、独立性、测度、代数、概率测度空间??基本公式：加法公式、乘法公式、全概率公式、Bayes公式概率模型：古典概型、几何概型、Bernoulli概型概率模型：古典概型、几何概型、Bernoulli概型样本空间随样本空间随机样本随机事件试验本点验???对立完備事件组,,,,互不相容,概率条件概率概率的定义及性质非负性可列可加性非负性完整性可列可加性其它运算公式可测空间、代数??古典几何概测空率度间Buffon可测空间、代数古典几何概型率度间Buffon问题概率公理化抽样概率公理化定义及计算抽样模型Bertrand悖论问题分房生日DeMere问题分赌本问题問题问题乘法公式样本空间划分条全概率公式条件概事件概率的独立性率Bayes公式独立性Bayes公式Bernoulli概型Bernoulli概型Bernoulli分布二几Pascal分布Bernoulli分布二项分几何分Pascal分布分咘分布多项分布随机游动布布随机游动经典例题赌徒破产模型SBanach火柴盒问题经典例题产模型火柴盒问题随机向量及概率分布随机向量及概率汾布基本概念：随机变量ＲV、分布律、概率密度pdf分布函数CDF随机变量函数的分布常用随机变量分布函数CDF、随机变量函数的分布常用随机变量：Bernoulli分布、Poisson分布、Geometric分布Bernoulli分布、Poisson分布、Geometric分布Uniform分布、Exponential分布、Guass分布RVDiscreteRV分布律几类常用的DiscreteRVDiscreteRV分布Ｂinomial分布Ｂinomial分布Poisson分布几何、超几何、负二项等RVContinuousRVRV常用的ContinuousRVpdfＥxponential分布Uniform汾布Normal分布p分布Nomemory性P{c<X<d}两个参数的意义查表随机变量的分布函数随机变量的分布函数单调不减性归性右连续性归一性F(x)…f(x)非负性非负性P{a<X<b}连续型随机變量的概率密度的概率密度随机变量的数字特征基本概念与定义：基本概念与定义：数学期望、方差、协方差、相关系数及矩阵形式数学期望、方差、协方差、相关系数及矩阵形式条件期望、GeneratingFunction、矩Moment特征函数、EntropyandInfermation、中心矩期望方差性质、常用随机向量的数学期望及方差数学期望嘚定义DiscreteContinuous数学期望的性质随机变量函数常用随机变量的期望常用随机变量的期望),(~pBernoulliXEXp?~(,)XBinomialnpEXnp?)(~?PoissonXEX??)(~pGeometricXEX?VarXpq?VarXnpq?VarX??qVarX?)(?PoissonXEX?)(pGeometricXEXp),(~baUniformXabEX??),(~??NormalXEX??VarX??VarXp?()baVarX??VarX??)(~?lExponentiaXEX??),(~??GammaXEX???VarX??VarX???理论与试验并重概率论试题与数理统计模拟仿真试验投色子试验投币试验离散型与连续型随机变量分布图形試验Poisson分布试验Poisson定理试验二项分布试验维与维正态分布试验CLT试验一维与二维正态分布试验CLT试验经验分布试验点估计的相合性与无偏性试验经驗分布试验点估计的相合性与无偏性试验假设检验中第第二类错误演示试验假设检验中第一、第二类错误演示试验二MtCl仿真实验二、MonteCarlo仿真实驗MCl仿真原理MonteCarlo方法的的基本思想是首先建立一个概率模型使所MonteCarlo仿真原理求问题的解正好是该模型参数或其他有关特征量然后通过模拟统计试驗即多次随机抽样试验（确定和）统计拟统计试验即多次随机抽样试验（确定m和n）,统计出某事件发生的百分比只要试验次数很大该百分比僦近拟于事件发生的概率。这实际上就是事件发生概率的统计定义利用建立的概率模型求出要估计的参数MonteCarlo属于试验数学建立的概率模型求出要估计的参数。MonteCarlo属于试验数学的一个分支Buffonproblem平面上画着一些间隔距离相等MonteCarlo仿真案例Buffonproblem平面上画着些间隔距离相等的平行线,它们之间的距離等于a（a＞）,在某一高度向此平面任投一长度为l（l＜a）的针,求此针与平行线相交的概率长度为l（l＜a）的针,求此针与平行线相交的概率解:以x表示针的中点到最***行线的距离,?表示该平行线与针的交角由于有ax??????故角,由于有,x??????故(,)|,aGxx??????????????(,)|,Gxx???????????xaalasinl?Glx???Og?事件A｛针与平行线相交｝发生il??即事件A=｛针与平行线相交｝发生?sinx??,即()|i()l???????(,)|sin,(,)lgxxxG????????????于是得于是得()()AAmSPA??sinldgl????的面积()()PAmS????aGa?????的面积注(Remark):①已知l、a则可求()PA①已知l、a,则可求()PA②已知()PA,l、a则可求?③计算机模拟:MonteCarlo方法McMc:MarkovchainMontecarlo方法③计算机模拟:MonteCarlo方法McMc:MarkovchainMontecarlo方法MATLAB实现Buffon问题仿真求解程序程序clearallL=针的长度平行线间的距离（）程序d=平行线间的距离（d>L）m=统计满足針与线相交条件的次数并赋初值n=投针试验次数fork=:n迭代次数fork=:n迭代次数x=unifrnd(,d)随机产生数的长度即投针之后针中点与平行线的距离p=unifrnd(,pi)随机产生的针与线相茭的角度ifx<=L*sin(p)针与线相交的条件m=m针与线相交则记数elseendendp=vpa(mn,)n次中与平行线相交的次数的频率比即相交的概率,vpa()以任意精度（位小数点默认值为位）显示出來pim=vpa((*L*n)(m*d),)利用投针频率估计圆周率pi,vpa()以任意精度（位小数点默认值为位）显示出来位小数点默认值为位）显示出来p=运行结果ppim=clearallN=循环迭代次数赋初值每佽循环迭之后的针与线相交的概率的记录值程序P=zeros(,N)赋初值每次循环迭之后的针与线相交的概率p的记录值Pim=zeros(,N)赋初值每次循环迭之后的圆周率pim的记錄值fori=:NL=针的长度d=平行线间的距离（d>L）d=平行线间的距离（d>L）m=统计满足针与线相交条件的次数并赋初值n=投针试验次数fork=:n迭代次数x=unifrnd(,d)随机产生数的长度即投针之后针中点与平行线的距离p=unifrnd(,pi)随机产生的针与线相交的角度ifx<=L*sin(p)针与线相交的充要条件m=m针与线相交则记数elseendendp=mnn次中与平行线相交的次数的频率仳即相交的概率pim=(*L*n)(m*d)利用投针频率估计圆周率piP(i)=p记录第i次循环之后的相交概率值P(,i)=p记录第i次循环之后的相交概率值Pim(,i)=(pim)记录第i次循环之后的圆周率pi值i=i进叺下次循环迭代endP=P无“”则显示每次的相交概率值Pim=Pim无“”则显示每次的圆周率pi值Pmean=mean(P)显示N次迭代之后的相交概率均值Pimmean=mean(Pim)显示N次迭代之后的圆周率pi均徝P运行结果Pmean=Pimmean=赌徒输光问题MonteCarlo仿真案例两个赌徒甲、乙两人将进行一系列赌博。在每一局中甲获胜的概率为p而乙获胜的概率为q(pq??)在每一局後失败者都要支付一元线给而乙获胜的概率为q(pq?)。在每局后失败者都要支付元线给胜利者在开始时甲拥有赌本a元而乙拥有赌本b元两个赌徒直到甲输光或乙输光为止。求甲输光的概率通过理论分析可知甲输光的概率是：b?()bbpqabP????????(),()babpqpqpq????????模拟赌博过程思路：在每一次模拟中,随机产生一个数,如果该数小于p,说时赌徒甲获胜,相应甲得到一元钱,而乙付出一元钱反之甲拿出一元钱给乙这里对甲的賭本a、乙的赌应甲得到元钱,而乙付出元钱反之甲拿出元钱给乙这里对甲的赌本a、乙的赌本b、甲赢的概率p取不同的数值进行次赌博过程模拟楿应程序如下：clcclearclosealla=甲的赌本b=乙的赌本p=甲赢的概率计数设置为S=计数设置为N=模拟次数m=设定随机数状态值（）改变这个值可以进行不同的实验rand('state',m)设置隨机数状态fork=:Nfork=:Nat=a初始化甲的赌本bt=b初始化乙的赌本whileat>bt>模拟整个赌博过程r=(rand<p)*算输赢at=atr交换赌本bt=btr交换赌本endS=S(at<)如甲输累加甲输的次娄endP=SN计算甲输的概率值g=ppPo=g^bg^(ab)返回甲输咣的概率理论值P=PPo=Binomial(二项分布)的使用MC案例Galton板实验一个级Galton板实验系统如下图(A)所示(源程序见后)和(B)所示(源程序见后)(源程序见后)和(B)所示(源程序见后)（A）（B）在图（A）中当小球从顶部向下降落时遇到第一层竖隔板此时小球分别向左右下落的概率各占一半（）当小球继续下落遇到第二层竖隔板时小球仍以左右下落的概率各占半（）当小球继续下落遇到第二层竖隔板时小球仍以左右相同的概率往下落以后每层均如此（如图（B）所示）。最后到了第层底部小球将落入底部个槽中的一个但是小球落入每个槽内的概率是不一样的如查将球将落入底部个槽中的一个但昰小球落入每个槽内的概率是不一样的。如查将这个个槽编号为：、、、、、、、、计算小球落入第个槽的概率理论分析这是个经典的②项分布概率模型考虑在第K层小球运动方向有两种理论分析：这是一个经典的二项分布概率模型。考虑在第K层小球运动方向有两种可能用kX表示则kX服从两点分布这里用kX?表示向右侧竖隔板方向运动用kX?表示向左侧竖隔板方向运动它们发生的概率均为。最终位置X由kkXX??决定即②项分布决定上述第层即有~(,)XBinomial上述该二k项分布的概率密度可用下列程序获得：bpdf=binopdf(:,,)运行如下：clearall源程序bpdf=binopdf(:,,)运行结果bpdf=运行结果注：第到第个槽内的概率徝即落入第个槽内的概率为下面用Matlab模拟小球下落到底部个槽内的概率分布在各层中用和分别表示向左和向右运动。在小球下落到底层的槽内时一个位的二进制数就完全确定了而所落入槽的编号应该是这个位二进制数各位数按十进制数相加的结果。重复模拟小球下落过程佽可以统计出小球落入各槽内的次数（即频球下落过程次可以统计出小球落入各槽内的次数（即频数）画出个频数数据的直方图。如下圖所示（源程序见后）注：灰色线条表示模拟值而红色线条表示理论值（其中把概率理论值的最大值调整为模拟值的最大值以便于对比）嘚最大值调整为模拟值的最大值以便于对比）(A)图源程序clearclccloseallfigurexlim(,)ylim(,)axisequalholdonL=No=levelnumberforN=:NomN=NoNfork=:Nfork=:Nplot(mNk*,,,NoN,'k','linewidth',)text(,NoN,numstr(N),'fontsize',,'fontname','timesnewroman')endendarg=linspace(,pi*,)col=,,fill(No*cos(arg),No*sin(arg),col,'Edgecolor',col)plot(No,No,,,No,No,,,'k')plot(No,No,No,No,No,No,,,'k')plot(,No,,,'k','linewidth',)fork=:Notext(knumstr(k)'fontsize'text(k,,numstr(k),'fontsize',,'fontname','timesnewroman')plot(k*cos(arg),*sin(arg),'k')endset(gcf,'color','w')axisoff（B）图源程序Galton板实验clearclccloseallfigurexlim(,)ylim()ylim(,)axisequalholdonL=No=levelnumberforN=:NomN=NoNfork=:Nplot(mNk*,,,NoN,'k','linewidth',)text(,NoN,numstr(N),'fontsize',,'fontname','timesnewroman')endendendarg=linspace(,pi*,)col=,,fill(No*cos(arg),No*sin(arg),col,'Edgecolor',col)plot(No,No,,,No,No,,,'k')plot(No,No,No,No,No,No,,,'k')plot(,No,,,'k','linewidth',)fkNfork=:Notext(k,,numstr(k),'fontsize',,'fontname','timesnewroman')plot(k*cos(arg),*sin(arg),'k')endset(gcf,'color','w')axisoff频数直方图模拟源程序clcclccloseallcleard('tt')固定随机数产生状态rand('state',)固定随机数产生状態R=unidrnd(,,)产生个~的随机数test=sum(R)对随机数求和统计各数出现的频数h=hist(test,)统计各数出现的频数bpdf=binopdf(:,,)计算二项分布出现各取值（即第至第槽内）的理论值bpdf=bpdfmax(bpdf)*max(h)使模拟数据與理论数据保持量级一致b=bar(:,hbpdf')绘制直方图set(b(),'facecolor','r')设置右侧的直线条填充颜色为红色set(b(),'facecolor',,,)把左侧直第填充颜色设置为灰色三、MonteCarlo积分随机投点法用MonteCarlo方法计算定積分用方法计算定积分（随机投点法）用MonteCarlo方法计算定积分（随机投点法）设()fx??求()Jfxdx??设(,)XY服从正方形{,}XY????上的Uniform分布(,){,}则~(,)XUnifrom~(,)YUnifrom,且相互独立记倳件{()}AYfX??{()}f则A的概率为()(){()}()fxpPAPYfXJdydxfxdxJ??????????(){()}()pPAPYfXJdydxfxdxJ????即定积分的值J就是事件A发生的概率p由BernoulliLLN我们可以用重复试验中由BernoulliLLN我们可以用重复试验中({()})AYfX??出现的频率作为概率p的估计值这({()})f出现频率作为概率p计种求定积分的方法也称为随机投点法,即将(,)XY看成是向正方形{,}XY????内随机投点,用隨机点落在区域{()}yfx?中的频率作为定积分的近拟值。MtCl方法来得到A出现的频率MonteCarlo方法来得到A出现的频率：()先用计算机产生(,)上Uniform分布的n个随机数:()先用計算机产生(,)上Uniform分布的n个随机数:,(,,,)iixyin??,这里n可以很大,如n?,甚至n?等()对n对数据(,),,,,iixyin??,记录满足不等式()iiyfx?的次数,这就是事件A发生的频数n?由此可得事件A发生的频率nn?则(())nJfxdxn????例：计算定积分()Jxxdx???源程序如下：源程序如下：clearallcloseallN=,,,,,输入参数N是随机投点的个数S=int('sqrt(x)x',,)面积的理论值（解析解）S=quad((x)sqrt(x)x,,)面积的悝论值（数值解）计算阴影区域的面积的MonteCarlo模拟值fori=:length(N)x=rand(N(i),)点的横坐标y=rand(N(i),)点的纵坐标yad((),)点的纵坐标n=sum(sqrt(x)>=yy>=x)落在阴影区域内点的频数Sn(i)=nN(i)落到阴影区域内点的频率即概率的模拟值endendSS运行结果S=Sn=运行结果例：计算定积分xJedx????源程序如下：clearallcloseallclearallcloseallN=,,,,,,输入参数N是随机投点的个数symsxS=int('exp(x^)sqrt(*pi)',,)面积的理论值（解析解）S(S)函数得出确定小數点位数后的确切解S=vpa(S,)vpa函数得出确定小数点位数后的确切解S=quad((x)exp(x^)sqrt(*pi),,)面积的理论值（数值解）计算阴影区域的面积的MonteCarlo模拟值filh(N)fori=:length(N)x=unifrnd(,,N(i),)点的横坐标y=unifrnd(,,N(i),)点的纵坐标n=sum(exp(x(i)^)sqrt(*pi)>=y)落茬阴影区域内点的频数Sn(i)=vpa(nN(i),)落到阴影区域内点的频率即概率的模拟值end运行结果SSnS运行结果S=Sn=源程序如下：源程序如下：clearallcloseallN=,,,,,,输入参数N是随机投点的个数個数symsxS=int('exp(x^)sqrt(*pi)',,)面积的理论值（解析解）S=vpa(S,)vpa函数得出确定小数点位数后的确切解S=quad((x)exp(x^)sqrt(*pi),,)面积的理论值（数值解）计算阴影区域的面积的MonteCarlo模拟值fori=:length(N)x=unifrnd(,,N(i),)点的横坐标y=unifrnd(,,N(i),)点的縱坐标n=sum(exp(x^)sqrt(*pi)>=y)落在阴影区域内点的频数Sn(i)=vpa(nN(i),)落到阴影区域内点的频率即概率的模拟值endS运行结果SnSS=Sn=,,,,,,对于一般区间,ab上的定积分()baJgxdx???a?作线性变换xayba???即鈳化成,区间上的积分进一步若可令()cgxd??,可令(){()}fygabaycdc?????则()fy??此时有()()()()()baJgxdxbadcfydycba?????????用上述MonteCarlo方法计算上式中()fydy?代入上式即可求出出()baJgxdx???例计算定积分xJedx???xax?作线性变换xaxyba?????即化成,区间上的积分但()xcegxeed???故令的积分但()cegxeed?????故令(){()}fygabaycd????dc?{exp()}yeee??????ee即()fy??此时有(){exp()}xJedxeeyedyeee??????????????按上述Montecarlo方示计算出上述积分项再代入即可编写程序实现J的随机投点法计算J的随机投点法计算clearallcloseallN=,,,,,,输叺参数N是随机投点的个数面积的理论值（解析解）symsxS=int('exp(x)',,)面积的理论值（解析解）S=vpa(S,)vpa函数得出确定小数点位数后的确切解symsxS=quad((x)(exp(x)),,)面积的理论值（数值解）S=vpa(S,)vpa函数得出确定小数点位数后的确切解计算积分项阴影区域的面积的MonteCarlo模拟值fori=:length(N)x=unifrnd(,,N(i),)点的横坐标y=unifrnd(N(i))点的纵坐标y=unifrnd(,,N(i),)点的纵坐标n=sum((exp(*x^)exp())(exp()exp())>=y)计算上述式子积分项中被积函數落在阴影区域内点的频数Sn(i)=vpa(nN(i),)上述积分项落到阴影区域内点的频率即概率的模拟值endSSn输出积分项值S=vpa(Sn^**(exp()exp())*exp(),)将上述式子中的积分项代回JSpa(S(ep()ep())ep(),)将述式子中的积汾项代回J中并取小数点后位值S=运行结果S=Sn=,,,,,,S=,,,,,,用MonteCarlo方法计算重积分见文献：谢中华MATLAB统计分析与应用个案例分析北京航空航天大学出版社P页北京航空航天大学出版社P页三、MonteCarlo积分平均值法用MonteCarlo方法计算定积分用MonteCarlo方法计算定积分（平均值法）用MonteCarlo方法计算定积分（平均值法）为计算定积分为计算定积分()Jfxdx??设则的数学期望设~(,)XUniform则()YfX?的数学期望()()EfXfxdxJ???()()ff?所以估计J的值就是估计()fX的数学期望的值由辛钦大数定律可以用()fX的观察值的平均去估计()fX钦大数定律,可以用()fX的观察值的平均去估计()fX的数学期望的值具体方法如下用计算机随机产生n个(,)上Uniform分布的随机数：,(,,,)ixin??,然后对每个ix计算()ifx,最後得J的估计值为()niiJfxn???i计算定积分xJedx??????编写程序实现clearallcloseallclearallcloseallN=,,,,,,输入参数N是随机产生的样本点的个数smsSint('ep(^)sqrt(*pi)')面积的理论值（解析解）symsxS=int('exp(x^)sqrt(*pi)',,)面积的理论值（解析解）S=vpa(S,)vpa函数得出确定小数点位数后的确切解S=quad((x)exp(x^)sqrt(*pi),,)面积的理论值（数值解）计算了随机数对应的函数的平均值计算了随机数对应的函数的平均值fori=:length(N)x=unifrnd(,,N(i),)随机产生点的横坐标应参数的数值求的初值J(i)=zeros()对应参数N的函数值求和的初值forj=:N(i)J(i)=J(i)vpa(exp(x(j)^)sqrt(*pi),)计算对应的函数值endSn(i)=vpa(J(i)N(i),)对应随机产生数N的函数平均值end运行结果SSn上述程序运行时间(h左右)较长运行结果S=Sn=,,,,,,例计算定积分xeJdxe????()xxxJedxeedx???????a．用MonteCarlo方法计算积分（随机投点法）J的随机投点法计算clearallcloseallN=,,,,,,输入参数N是隨机投点的个数symsxS=int('(exp(x))(exp())',,)面积的理论值（解析解）S=quad((x)(exp(x))(exp()),,)面积的理论值（数值解）S=vpa(S,)vpa函数得出确定小数点位数后的确切解计算阴影区域的面积的MonteCarlo模拟值filth(N)fori=:length(N)x=unifrnd(,,N(i),)点的橫坐标y=unifrnd(,,N(i),)点的纵坐标n=sum((exp(x))(exp())>=y)落在阴影区域内点的频数nsum((exp(x))(exp())>y)落在阴影区域内点的频数Sn(i)=vpa(nN(i),)落到阴影区域内点的频率即概率的模拟值endS运行结果Sn运行结果S=Sn=,,,,,,J的随机投點法计算clearallcloseall编写程序实现clearallcloseallN=,,,,,,输入参数N是随机投点个数symsxS=int('exp(x)',,)面积的理论值（解析解）S=vpa(S,)vpa函数得出确定小数点位数后的确切解symsxS=quad((x)(exp(x)),,)面积的理论值（数值解）S=vpa(S,)vpa函數得出确定小数点位数后的确切解计算阴影区域的面积的MonteCarlo模拟值fori=:length(N)fori:length(N)x=unifrnd(,,N(i),)点的横坐标y=unifrnd(,,N(i),)点的纵坐标n=sum((exp(*x^)exp())(exp()exp())>=y)落在阴影区域内点的频数Sn(i)=vpa(nN(i),)落到阴影区域内点的频率即概率的模拟值endendSSnS=vpa(Sn^**(exp()exp())*exp(),)S=Sn=,,,,,,运行结果S=,,,,,,b．用MonteCarlo方法计算积分（平均值法）J的平均值法计算J的平均值法计算clearallcloseallN=,,,,,,输入参数N是随机产生的样本点的个数编写程序实現本点的个数symsxS=int('(exp(x))(exp())',,)面积的理论值（解析解）S=vpa(S,)vpa函数得出确定小数点位数后的确切解d(()(())(()))面积的理论值（数值解）S=quad((x)(exp(x))(exp()),,)面积的理论值（数值解）S=vpa(S,)vpa函数得出确萣小数点位数后的确切解计算了随机数对应的函数的平均值fori=:length(N)x=unifrnd(,,N(i),)随机产生点的横坐标J(i)=zeros()对应参数N的函数值求和的初值forj=:N(i)J(i)=J(i)(exp(x(j)))(exp())计算对应的函数值endSn(i)=vpa(J(i)N(i),)对应随机產生数N的函数平均值endS运行结果SnS=Sn=,,,,,,运行结果J的平均值法计算llllll编写程序实现clearallcloseallN=,,,,,,输入参数N是随机产生的样本点的个数symsxS=int('exp(x)',,)面积的理论值（解析解）S=vpa(S,)vpa函数得絀确定小数点位数后的确切解S=quad((x)exp(x),,)面积的理论值（数值解）S=vpa(S,)vpa函数得出确定小数点位数后的确切解计算了随机数对应的函数的平均值fori=:length(N)x=unifrnd(,,N(i),)随机产生点嘚横坐标J(i)=zeros()对应参数N的函数值求和的初值forj=:N(i)J(i)=J(i)exp(x(j))exp(x(j))计算对应的函数值endSn(i)=vpa(J(i)N(i),)对应随机产生数N的函数平均值Sn(i)vpa(J(i)N(i),)对应随机产生数N的函数平均值endSSn运行结果SnS=Sn=,,,,,,四、随机数嘚产生四、随机数的产生四、随机数的产生随机数的定义和特性随机数的定义和特性四、随机数的产生什么是随机数?单个的数字不是隨机数?是指个数列其中的每个体称为随机数其值与数列中?是指一个数列其中的每一个体称为随机数其值与数列中的其它数无关?在一個均匀分布的随机数中每一个体出现的概率是均等的?例在,区间上均匀分布的随机数序列中与出现的机会均等出现的机会均等随机数应具囿的基本特性?考虑一个对高能粒子反应过程的模拟：需用随机数确定：出射粒的属性能量方向?出射粒子的属性：能量、方向、…?粒孓与介质的相互作用?对这一过程的模拟应满足以下要求（相空间产生过程）：?出射粒子的属性应是互不相关的即每一粒子的属?出射粒子的属性应是互不相关的即每一粒子的属性的确定独立于其它的粒子的属性的确定?粒子的属性的分布应满足物理所要求的理论分布?粒子的属性的分布应满足物理所要求的理论分布?所模拟的物理过程要求随机数应具有下列特性：?随机数序列应是独立的、互不相关的(uncorrelated)：即序列中的任一子序列应与其它的子序列无关?长的周期(longperiod)：gp实际应用中随机数都是用数学方法计算出来的这些算法具有周期性即当序列達到定长度后会重复算法具有周期性即当序列达到一定长度后会重复?均匀分布的随机数应满足均匀性(Uniformity)：?均匀分布的随机数应满足均匀性(Uniformity)：随机数序列应是均匀的、无偏的即：如果两个子区间的“面积”相等则落于这两个子区间内的随机数的个数应相等。例如对)区间均匀汾布的随机数如果产生了足够多的随机例如：对,)区间均匀分布的随机数如果产生了足够多的随机数而有一半的随机数落于区间,?不满足均勻性如果均匀性不满足则会出现序列中的多组随机数相关的情况?均匀性与互不相关的特性是有联系的?有效性（Efficiency):模拟结果可靠?模拟产苼的样本容量大?模拟产生的样本容量大?所需的随机数的数量大?随机数的产生必须快速、有效最好能够进行并行计算随机数的产生隨机数的产生?,区间上均匀分布的随机数是MonteCarlo模拟的基础:?服从任意分布的随机数序列可以用,区间均匀分布的随机数序列作适当的变换或舍選后求得?,均匀分布的随机数的产生方法：随机数序列作适当的变换或舍选后求得?利用一些具有内在的随机性过程：?放射性衰变过程（radioactivedecay)?放射性衰变过程（radioactivedecay)?热噪声(thermalnoise)?宇宙线的到达时间（cosmicrayarrival)?宇宙线的到达时间（cosmicrayarrival)?…?缺点：模拟的结果不可再现使得模拟程序的找错困难缺点模拟的可再现使得模拟程的找错难?利用事先制订好的随机数表:?缺点：表的容量有限不适合需要大量随机数的应用?利用数学递推公式在计算机中产生随机数(伪随机数)),,,(?????knnnknrrrTr其中：T为某个函数给定初值r,r,…,rk,可按上式确定rn,n=,,…?随机数序列算法：产生M区间上的整数In然后利用公式rn=InM返回,区间上的实数nn返回,区间上的实数优点：占用计算机的内存少–占用计算机的内存少–产生速度快复前次的模拟结便程序的找–可以重复前次的模拟结果便于程序的找错缺点：?不满足随机数之间相互独立的要求：公式和初值确定后序列就唯一地确定了?不满足均匀性：计算机能表示的,区间内的数是有限的（由字长确定）?伪随机数（PdRdNb)?递推到一定次数后出现周期性的重复现象?伪随机数（PseudoRandomNumber)线性塖同余方法线性乘同余方法（（LinearLinearCongruentialCongruentialMethod)Method)线性乘同余方法线性乘同余方法（（LinearLinearCongruentialCongruentialMethod)Method)年Lehmer提出的一种产生伪随机数的方法,是最常用的方法。A、递推公式：mcaIInnmod)(???A、递推公式：其中：I：初始值（种子seed)caImca,??：初始值（种子)a：乘法器（multiplier)c：增值（additiveconstant)：模数（caIm,,?m：模数（modulus)mod：取模运算：(aInc)除以m后的余数a,c和m皆为整数?产生整型的随机数序列,随机性来源于取模运算产生整型的随机数序列,随机性来源于取模运算B、实型随机数序列：如果c=?乘同余法：速度更快也可产生长的随机数序列),)(??fltIrnn?mInB、实型随机数序列：,)()(???mfloatIrmfloatnn??mInn)(mfloatC、特点：“摇号”抽签等个中抽取小于等个以）最大容量为m：mIn??等于个可以但大于则重号！）独立性和均匀性取决于参数a和c的选择例：a=c=I=,m=?,,,,,,,,…重复、不连续D、模数m的选择：?m应尽可能地大因为序列的周期不可能大于m?通常将m取为计算机所能表示的最大的整型量在位计通常将m取为计算机所能表示的最大的整型量在位计算机上m==xE、乘数因子a的選择：Model关键E、乘数因子a的选择：年MGreenberger证明：用线性乘同余方法产生的随机数序列具有周期的条件是：机数序列具有周期m的条件是：c和m为互质數c和m为互质数a是质数p的倍数,其中p是a和m的共约数如果m是的倍数,a也是的倍数随机分配例：a=,c=,m=,I=?周期=m=分配编号,,,,,,,,,,,,,,,,,,,,,,伪随机数的产生“摇号”或“分配编號”的难度在于总数即M变种子数即I初值在变核发指标即拟中签数事先可确定均匀性、随机性与独立性入学(幼儿园、小学、中学等)摇号买房、购车的摇号┇中国进入了摇号时代！！┇中国进入了摇号时代伪随机数的产生国内已开发的摇号系统伪随机数的产生国内已开发的摇號系统伪随机数的产生国内已开发的摇号系统伪随机数的产生国内已开发的摇号系统伪随机数的产生国内已开发的摇号系统我们的摇号系統摇号编号分配系统我们的摇号系统摇号系统五、建模案例分析五、建模案例分析摇号编码分配设计．摇号编号分配系统设计算法第步从報名系统中读有效申请编并统总数第一步：从报名系统中读取有效申请编码并统计总数M申编位位申请编码位数位：XXXXXXXXX后九位数据为中的任意隨机数第二步：GBT标准第二步：GBT标准乘同余法MultiplicativeCongruentialMethod迭代()(mod)nnXaXCM??????(,)nnXRUniformM???????,,,nM??XSeednumber???随机种子数取值随机种子数取值第三步：摇号编號生成通过上述国家标准模型随机产生申请编码对应的唯一随机摇号编号其摇号编号是一组连续数即为~M。通过摇号编号分配系统随机生成嘚编号该编号将进入摇号系统的摇号数即摇号系统的摇号总数即为M其摇号系统摇出的编号数个数即为该期专段号牌核发小客汽车指标数个數即为该期专段号牌核发小客汽车指标数抽样方式：无放回抽样抽样方式无放回抽样需克服的关键技术：需克服的关键技术：?随机模型生成的摇号编号不重号、不漏号?保证摇号编号的随机性、独立性与重现性?避免报名先后与随机编号即摇号编号无相关性。建模过程伍、建模案例分析五、建模案例分析摇号系统中签设计．摇号系统设计算法第一步：基本参数输入该期核发号牌有效编码总数M该期核发该期核发号牌有效编码总数M、该期核发小客汽车专段号牌指标数K与随机种子数X小客汽车专段号牌指标数K与随机种子数X第二步：迭代GBT标准模型塖同余法MultiplicativeCongruentialMethod()(mod)nnXaXCMX??????(,)nnXRUniformMXSeednumber???????,,,nk??XSeednumber???第三步摇号系统生成中签第三步：摇号系统生成中签通过上述国家标准模型随机生成中嘚通过上述国家标准模型随机生成~M中的k个摇号编号即为小客汽车专段号牌中签编号k个摇号编号即为小客汽车专段号牌中签编号抽样方式：无放回简单随机抽样抽样方式：无放回简单随机抽样需克服的关键技术需克服的关键技术：随机模型生成的中签编号具有均匀性随机性．随机模型生成的中签编号具有均匀性、随机性．保证中签编号的不可重复性．每期核发指标的中签编号具有重现性无“特权号”与“不圉号”．无特权号与不幸号．每个有效编号被等概率抽中每个种子数后的组合被等组合概率抽中。建模过程