《R语言中随机模拟的例子》由会员分享,可在线阅读,更多相关《R语言中随机模拟的例子(40页珍藏版)》请在金锄头文库上搜索。
1、R语言中几个简单随机模拟的 应用例子赵明扬 2011.3.28问题一车和羊的游戏 假设你在进行一个游戏节目。现给三扇门供你选择:一扇门后面 是一辆轿车,另两扇门后面分别都是一头山羊。你的目的当然是要想 得到比较值钱的轿车,但你却并不能看到门后面的真实情况。主持人 先让你作第一次选择。在你选择了一扇门后,知道其余两扇门后面是 什么的主持人,打开了另一扇门给你看,而且,当然,那里有一头山 羊。现在主持人告诉你,你还有一次选择的机会。那么,请你考虑一 下,你是坚持第一次的选择不变,还是改变第一次的选择,更有可能 得到轿车? 广场杂志刊登出这个题目后,竟引起全美大学生的举国辩论, 许多大学的教授们也参
2、与了进来。真可谓盛况空前。据纽约时报 报道,这个问题也在中央情报局的办公室内和波斯湾飞机驾驶员的营 房里引起了争论,它还被麻省理工学院的数学家们和新墨哥州洛斯阿 拉莫斯实验室的计算机程序员们进行过分析。问题分析 在一次实验中,如果第一次选择选中了轿车(概率为1/3),那么 主持人打开一扇门后,如果坚持原来的选择,则能得到轿车,反之, 改变第一次选择则不能得到轿车。如果第一次没有选中轿车(概率为 2/3),那么其余两扇门后面必有一个是轿车,主持人只能打开有山 羊有那扇门,则剩下的一扇门后面是轿车,此时坚持原来的选择不能 得到轿车,改变第一次的选择必能得到轿车。因此,经过分析,坚持 第一次的选择不
3、变得到轿车的概率为1/3,改变第一次的选择得到轿 车的概率为2/3。 实际上,在只有三扇门的情况下,那么改不改变选择效果并不明 显。如果有100扇门,参与的嘉宾选择了其中的一扇,而主持人随后 把剩下的99扇门中间的98扇门都打开,这98扇门后面都没有奖品, 这时应该改变选择,毕竟最开始自己选择的那扇门中奖的概率只是 1%而已。 需要注意的是,主持人是在知道其他两 扇门后面都有什么的情况下选择一个门打 开的。这种情况下三个门后是轿车的概率 因为主持人知道结果并参与其中而关联在 一起,而不是孤立等同的。如果打开门的 不是主持人,而是另一个参与者,并且当 他打开门时发现什么也没有,那么,剩下 的两个
4、门后是轿车的概率才是相等的。 计算机模拟 为了验证这一结果,我们就要比较不改变选择中奖的 几率和改变选择中奖的几率。模拟方法是:我们从0,1, 2这3个数中随机一个为轿车(即中奖号码),另随机一个 数为你的选择。如果你的选择与中奖号相同,则计这次为 不改变选择中奖;如果你的选择不对,则是改变选择中奖 。分别累积出不改变选择中奖和改变选择中奖的次数,就 可以得到不改变选择中奖的几率和改变选择中奖的几率了 。 为了将结果表示的明显,我们可以假设有100扇门,参 与的嘉宾选择了其中的一扇,而主持人随后把剩下的99扇 门中间的98扇门都打开,这98扇门后面都没有奖品,然后 模拟并比较不改变选择中奖的几
5、率和改变选择中奖的几率 。此时的情况也是相同的,只是每次随即都是从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; pr
6、int(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
7、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=f
8、loor(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 tr
9、ial 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
10、0.9903889 当模拟的次数逐渐的增多时,其模拟值越接近理 论值,这说明模拟的效果越好。在随机事件的大 量重复出现中,往往呈现几乎必然的规律,这个 规律就是大数定律。通俗地说,这个定理就是, 在试验不变的条件下,重复试验多次,随机事件 的频率近似于它的概率。 偶然必然中包含着必然。此次模拟试验也正好 用实际的模拟例子说明了大数定理的正确性和应 用性。问题二蒲丰投针问题 1777年法国科学家蒲丰提出的一种计 算圆周率的方法-随机投针法,这种方法 的具体步 骤是:1、取一张白纸,在上面画上许 多条间距为d的平行直线。 2、取一根长为l(l proc.time()-mm 用户 系统 流逝 11.
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只有一个人离开,另一人
12、 仍在接受服务。然而,由指数分布的无记忆性,从现在起剩下的这个 人被服务的时间是参数为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(run
13、if(1)/lambda; t3=-log(runif(1)/lambda; if(t1+t3T3 ,假设T1,T2,T3均是随机变量 ,且T2N(30,22),设r1,r2是(0,1)区间上均匀 分布的随机数,则T1和T3的分布律的模拟 公式为 那么t1,t2可以看成T1,T3的一个观察值. 令t2是服从正态分布N(30,22)的随机数,则 将t2看成火车运行时间T2的一个观测值. 在每次试验中,产生两个U(0, 1)的随机数 t1,t3,一个N(30, 22)的随机数t2,当 t1+t2t3,认为试验成功(能够赶上火车).若 在了,次试验中,有k次成功,则用频率k/n 作为此人赶上火车的概
14、率.当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(r1it2k) 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,默认为10f=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; #计算落入制