R语言中随机模拟的例子PPT资料.ppt

1、,问题分析 在一次实验中，如果第一次选择选中了轿车（概率为1/3），那么主持人打开一扇门后，如果坚持原来的选择，则能得到轿车，反之，改变第一次选择则不能得到轿车。如果第一次没有选中轿车（概率为2/3），那么其余两扇门后面必有一个是轿车，主持人只能打开有山羊有那扇门，则剩下的一扇门后面是轿车，此时坚持原来的选择不能得到轿车，改变第一次的选择必能得到轿车。因此，经过分析，坚持第一次的选择不变得到轿车的概率为1/3，改变第一次的选择得到轿车的概率为2/3。实际上，在只有三扇门的情况下，那么改不改变选择效果并不明显。如果有100扇门，参与的嘉宾选择了其中的一扇，而主持人随后把剩下的99扇门中间的98扇

2、门都打开，这98扇门后面都没有奖品，这时应该改变选择，毕竟最开始自己选择的那扇门中奖的概率只是1%而已。,需要注意的是，主持人是在知道其他两扇门后面都有什么的情况下选择一个门打开的。这种情况下三个门后是轿车的概率因为主持人知道结果并参与其中而关联在一起，而不是孤立等同的。如果打开门的不是主持人，而是另一个参与者，并且当他打开门时发现什么也没有，那么，剩下的两个门后是轿车的概率才是相等的。,计算机模拟 为了验证这一结果，我们就要比较不改变选择中奖的几率和改变选择中奖的几率。模拟方法是：我们从0，1，2这3个数中随机一个为轿车（即中奖号码），另随机一个数为你的选择。如果你的选择与中奖号相同，则计这

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,

4、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+1else alter=alter+1;stick=stick/n;alter=alter/n;print（data.frame（trial=n,stick=stick,alter=alter）;9）f（i*104）,结果为：trial stick alter1 10000 0.3324 0.6676 trial stick alter1 20000 0.3289 0.6711 tri

5、al stick alter1 30000 0.3313333 0.6686667 trial stick alter1 40000 0.331475 0.668525 trial stick alter1 50000 0.33258 0.66742 trial stick alter1 60000 0.3345167 0.6654833 trial stick alter1 70000 0.3325143 0.6674857 trial stick alter1 80000 0.3316125 0.6683875 trial stick alter1 90000 0.3320444 0.66

6、79556,对应的100个选择的程序为f=function（n）stick=0;n）x=floor（runif（1,min=0,max=100）+1）;y=floor（runif（1,min=0,max=100）+1）;trial stick alter1 10000 0.0094 0.9906 trial stick alter1 20000 0.01005 0.98995 trial stick alter1 30000 0.009533333 0.9904667 trial stick alter1 40000 0.0101 0.9899 trial stick alter1 50000

7、 0.01018 0.98982 trial stick alter1 60000 0.0105 0.9895 trial stick alter1 70000 0.009542857 0.9904571 trial stick alter1 80000 0.0100875 0.9899125 trial stick alter1 90000 0.009611111 0.9903889,当模拟的次数逐渐的增多时，其模拟值越接近理论值，这说明模拟的效果越好。在随机事件的大量重复出现中，往往呈现几乎必然的规律，这个规律就是大数定律。通俗地说，这个定理就是，在试验不变的条件下，重复试验多次，随机事件

8、的频率近似于它的概率。偶然必然中包含着必然。此次模拟试验也正好用实际的模拟例子说明了大数定理的正确性和应用性。,问题二蒲丰投针问题,1777年法国科学家蒲丰提出的一种计算圆周率的方法-随机投针法，这种方法的具体步 骤是：1、取一张白纸，在上面画上许多条间距为d的平行直线。2、取一根长为l（ld）的针，随机向画有平行直线的纸上掷n次，观察针与直线相交的次数，记为m。3、计算针与直线相交的概率,问题分析 针与平行线相交的充要条件是：x l*sin/2,其中0 xd/2,0 建立直角坐标系（，x），上述条件在坐标系下将是曲线所围成的曲边梯形区域：,由几何概率知：,程序为：f=function（n）d

9、=2;l=1.5;m=0;n）x=runif（1,min=0,max=d/2）;theta=runif（1,min=0,max=pi）;if（x=l*sin（theta）/2）m=m+1;estimate=2*l*n/（m*d）;result=data.frame（num=n,estimate=estimate）print（result）;for（i in 2:5）f（10i）,结果为：num estimate1 100 3.125 num estimate1 1000 2.946955 num estimate1 10000 3.111388 num estimate1 1e+05 3.14

10、327 可见，模拟的次数越多，所得的结果就越准确。这与理论所知的结果是一样的。,问题三掷骰子问题,掷三只相同的色子，求至少有两个点数相同的概率。问题分析 构造3个整数随机变量，使之服从1到6之间的均匀分布，比较每次模拟出的3个值，如果至少有两个相同，则将成功的次数加1，分析可知，成功的频率即为所求的概率。f=function（n）success=0;n）x=floor（runif（1）*6+1）;y=floor（runif（1）*6+1）;z=floor（runif（1）*6+1）;if（x=y|x=z|y=z）success=success+1;p=success/n;print（data.

11、frame（num=n,p=p）;mm=proc.time（）;9）f（i*104）;proc.time（）-mm,结果为：num p1 20000 0.4447 num p1 30000 0.4428333 num p1 40000 0.441525 num p1 50000 0.44162 num p1 60000 0.4456667 num p1 70000 0.4445857 num p1 80000 0.4463875 num p1 90000 0.4451111 proc.time（）-mm 用户 系统 流逝 11.19 0.03 14.23 由概率知识可知：P=1-（4*5*6）

12、/（6*6*6）=5/9=0.4445。这与以上模拟出的值基本上是一样的。,问题四无记忆性的例子,考虑有两个营业员的邮局。假设当Smith先生进入邮局的时候，他看到两个营业员分别在为Jones先生和Brown先生服务，并且被告知先服务完的营业员即将为他服务。再假设为每顾客服务的时间都服从参数为lambda的指数分布。1）.求这三个顾客中Smith最后离开邮局的概率？Jones和Brown呢？2）.求Smith在邮局花费的平均时间ET？问题分析 当Smith接受服务的时候，Jones和Brown只有一个人离开，另一人仍在接受服务。然而，由指数分布的无记忆性，从现在起剩下的这个人被服务的时间是参数

13、为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;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

14、;else if（t2+t3t1）jones=jones+1;t=t2+t3;else smith=smith+1;t=t3+min（t1,t2）;sumT=sumT+t;simth=smith/n;jones=jones/n;brown=brown/n;ET=sumT/n;ETreal=3/（2*lambda）;#真是的值print（data.frame（smith=smith,jones=jones,brown=brown,ET=ET,ETreal=ETreal）for（i in 1:9）f（3,i*103）,结果为：smith jones brown ET ETreal1 506 0.232 0.262 0.493107 0.5 smith jones brown ET ETreal1 979 0.264 0.2465 0.4938358 0.5 smith jones brown ET ETreal1 1459 0.25