单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,*,R,语言中几个简单随机模拟的应用例子,赵明扬,2011.3.28,问题一车和羊的游戏,假设你在进行一个游戏节目。现给三扇门供你选择：一扇门后面是一辆轿车，另两扇门后面分别都是一头山羊。你的目的当然是要想得到比较值钱的轿车，但你却并不能看到门后面的真实情况。主持人先让你作第一次选择。在你选择了一扇门后，知道其余两扇门后面是什么的主持人，打开了另一扇门给你看，而且，当然，那里有一头山羊。现在主持人告诉你，你还有一次选择的机会。那么，请你考虑一下，你是坚持第一次的选择不变，还是改变第一次的选择，更有可能得到轿车？,广场杂志,刊登出这个题目后，竟引起全美大学生的举国辩论，许多大学的教授们也参与了进来。真可谓盛况空前。据,纽约时报,报道，这个问题也在中央情报局的办公室内和波斯湾飞机驾驶员的营房里引起了争论，它还被麻省理工学院的数学家们和新墨哥州洛斯阿拉莫斯实验室的计算机程序员们进行过分析。,问题分析,在一次实验中，如果第一次选择选中了轿车（概率为,1/3,），那么主持人打开一扇门后，如果坚持原来的选择，则能得到轿车，反之，改变第一次选择则不能得到轿车。如果第一次没有选中轿车（概率为,2/3,），那么其余两扇门后面必有一个是轿车，主持人只能打开有山羊有那扇门，则剩下的一扇门后面是轿车，此时坚持原来的选择不能得到轿车，改变第一次的选择必能得到轿车。因此，经过分析，坚持第一次的选择不变得到轿车的概率为,1/3,，改变第一次的选择得到轿车的概率为,2/3,。,实际上，在只有三扇门的情况下，那么改不改变选择效果并不明显。如果有,100,扇门，参与的嘉宾选择了其中的一扇，而主持人随后把剩下的,99,扇门中间的,98,扇门都打开，这,98,扇门后面都没有奖品，这时应该改变选择，毕竟最开始自己选择的那扇门中奖的概率只是,1%,而已。,需要注意的是，主持人是在知道其他两扇门后面都有什么的情况下选择一个门打开的。这种情况下三个门后是轿车的概率因为主持人知道结果并参与其中而关联在一起，而不是孤立等同的。如果打开门的不是主持人，而是另一个参与者，并且当他打开门时发现什么也没有，那么，剩下的两个门后是轿车的概率才是相等的。,计算机模拟,为了验证这一结果，我们就要比较不改变选择中奖的几率和改变选择中奖的几率。模拟方法是：我们从,0,，,1,，,2,这,3,个数中随机一个为轿车（即中奖号码），另随机一个数为你的选择。如果你的选择与中奖号相同，则计这次为不改变选择中奖；如果你的选择不对，则是改变选择中奖。分别累积出不改变选择中奖和改变选择中奖的次数，就可以得到不改变选择中奖的几率和改变选择中奖的几率了。,为了将结果表示的明显，我们可以假设有,100,扇门，参与的嘉宾选择了其中的一扇，而主持人随后把剩下的,99,扇门中间的,98,扇门都打开，这,98,扇门后面都没有奖品，然后模拟并比较不改变选择中奖的几率和改变选择中奖的几率。此时的情况也是相同的，只是每次随即都是从,0,到,99,中随机数而已。,程序：,f=function(n),stick=0;,alter=0;,x=floor(runif(n,min=0,max=3)+1);,y=floor(runif(n,min=0,max=3)+1);,for(i in 1:n),x=floor(runif(1,min=0,max=3)+1);,y=floor(runif(1,min=0,max=3)+1);,if(x=y),stick=stick+1,else,alter=alter+1;,stick=stick/n;alter=alter/n;,print(data.frame(trial=n,stick=stick,alter=alter);,for(i in 1:9)f(i*104),结果为：,trial stick alter,1 10000 0.3324 0.6676,trial stick alter,1 20000 0.3289 0.6711,trial stick alter,1 30000 0.3313333 0.6686667,trial stick alter,1 40000 0.331475 0.668525,trial stick alter,1 50000 0.33258 0.66742,trial stick alter,1 60000 0.3345167 0.6654833,trial stick alter,1 70000 0.3325143 0.6674857,trial stick alter,1 80000 0.3316125 0.6683875,trial stick alter,1 90000 0.3320444 0.6679556,对应的,100,个选择的程序为,f=function(n),stick=0;,alter=0;,x=floor(runif(n,min=0,max=3)+1);,y=floor(runif(n,min=0,max=3)+1);,for(i in 1:n),x=floor(runif(1,min=0,max=100)+1);,y=floor(runif(1,min=0,max=100)+1);,if(x=y),stick=stick+1,else,alter=alter+1;,stick=stick/n;alter=alter/n;,print(data.frame(trial=n,stick=stick,alter=alter);,for(i in 1:9)f(i*104),结果为：,trial stick alter,1 10000 0.0094 0.9906,trial stick alter,1 20000 0.01005 0.98995,trial stick alter,1 30000 0.009533333 0.9904667,trial stick alter,1 40000 0.0101 0.9899,trial stick alter,1 50000 0.01018 0.98982,trial stick alter,1 60000 0.0105 0.9895,trial stick alter,1 70000 0.009542857 0.9904571,trial stick alter,1 80000 0.0100875 0.9899125,trial stick alter,1 90000 0.009611111 0.9903889,当模拟的次数逐渐的增多时，其模拟值越接近理论值，这说明模拟的效果越好。在随机事件的大量重复出现中，往往呈现几乎必然的规律，这个规律就是大数定律。通俗地说，这个定理就是，在试验不变的条件下，重复试验多次，随机事件的频率近似于它的概率。,偶然必然中包含着必然。此次模拟试验也正好用实际的模拟例子说明了大数定理的正确性和应用性。,问题二蒲丰投针问题,1777,年法国科学家蒲丰提出的一种计算圆周率的方法,-,随机投针法，这种方法的具体步,骤是：,1,、取一张白纸，在上面画上许多条间距为,d,的平行直线。,2,、取一根长为,l(ld),的针，随机向画有平行直线的纸上掷,n,次，观察针与直线相交的次数，记为,m,。,3,、计算针与直线相交的概率,问题分析,针与平行线相交的充要条件是：,x l*sin/2,其中,0 xd/2,0,建立直角坐标系,(,，,x),，上述条件在坐标系下将是曲线所围成的曲边梯形区域：,由几何概率知：,程序为：,f=function(n),d=2;l=1.5;,m=0;,for(i in 1:n),x=runif(1,min=0,max=d/2);,theta=runif(1,min=0,max=pi);,if(x proc.time()-mm,用户 系统 流逝,11.19 0.03 14.23,由概率知识可知：,P=1-(4*5*6)/(6*6*6)=5/9=0.4445,。这与以上模拟出的值基本上是一样,的。,问题四无记忆性的例子,考虑有两个营业员的邮局。假设当,Smith,先生进入邮局的时候，他看到两个营业员分别在为,Jones,先生和,Brown,先生服务，并且被告知先服务完的营业员即将为他服务。再假设为每顾客服务的时间都服从参数为,lambda,的指数分布。,1).,求这三个顾客中,Smith,最后离开邮局的概率？,Jones,和,Brown,呢？,2).,求,Smith,在邮局花费的平均时间,ET,？,问题分析,当,Smith,接受服务的时候，,Jones,和,Brown,只有一个人离开，另一人仍在接受服务。然而，由指数分布的无记忆性，从现在起剩下的这个人被服务的时间是参数为,lambda,的指数分布，即与,Smith,是同分布的。再由对称性，他与,Smith,先离开的概率一样，都为,1/2,。第,1,问答案为,1/2,1/4,1/4,。,ET=Emin X;Y+S=1/(2*)+1/,其中,S,为,Smith,服务花费的时间,X;Y,分别为,Jones,和,Brown,服务所花费的时间。,程序为：,f=function(lambda,n=105),jones=0;,brown=0;,smith=0;,sumT=0;,for(i in 1:n),t1=-log(runif(1)/lambda;,t2=-log(runif(1)/lambda;,t3=-log(runif(1)/lambda;,if(t1+t3t2),brown=brown+1;,t=t1+t3;,else if(t2+t3T3,，假设,T1,T2,T3,均是随机变量，,且,T2N(30,2,2,),设,r1,r2,是,(0,1),区间上均匀分布的随机数，则,T1,和,T3,的分布律的模拟公式为,那么,t1,t2,可以看成,T1,T3,的一个观察值,.,令,t2,是服从正态分布,N(30,2,2,),的随机数，则将,t2,看成火车运行时间,T2,的一个观测值,.,在每次试验中，产生两个,U(0,1),的随机数,t1,t3,，一个,N(30,22),的随机数,t2,，当,t1+t2t3,，认为试验成功,(,能够赶上火车,).,若在了，次试验中，有,k,次成功，则用频率,k/n,作为此人赶上火车的概率,.,当,n,很大时，频率值与概率值近似相等,.,f=function(n),r1=runif(n);,r2=runif(n);,t2=rnorm(n,30,sd=2);,t1=0;t3=0;,num=0;,for(i in 1:n),if(r1i=0.7),t1i=0,else if(r1i=0.9),t1i=5,else t1i=10;,;,for(j in 1:n),if(r2j=0.3),t3j=28,else if(r2j=0.7),t3j=30,else if(r2jt2k),num=num+1;,p=num/n;,print(p);,for(i in 1:6)f(i*104),结果为：,1 0.661,1 0.65325,1 0.659,1 0.655225,1 0.65566,1 0.6569833,问题八计算积分问题,#,用,hit or miss,法计算积分值,#,要积分的函数,g(x),h=function(x)sqrt(1-x2);,#,给定一个比,g(x),在制定积分区间最大值大的一个常量,c,#,积分区间为,c(a,b),积分中使用的随机数个数为,n,默认为,10,f=function(g,a,b,c,n=10),x=runif(n,a,b);#,在区间,c(a,b),产生,n,个均匀随机数,y=c*runif(n,0,1);#,产生区间,c(0,c),之间的随机数,p=sum(g(x)=y)/n;#,计算落入制定区域内的频率,S=p*c*(b-a);SE=c*(b-a)*sqrt(