《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; #计算落入制
