R语言中随机模拟的例子.ppt
- 文档编号:832067
- 上传时间:2023-04-30
- 格式:PPT
- 页数:40
- 大小:238KB
R语言中随机模拟的例子.ppt
《R语言中随机模拟的例子.ppt》由会员分享,可在线阅读,更多相关《R语言中随机模拟的例子.ppt(40页珍藏版)》请在冰点文库上搜索。
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(iin1:
n)x=floor(runif(1,min=0,max=3)+1);y=floor(runif(1,min=0,max=3)+1);if(x=y)stick=stick+1elsealter=alter+1;stick=stick/n;alter=alter/n;print(data.frame(trial=n,stick=stick,alter=alter);for(iin1:
9)f(i*104),结果为:
trialstickalter1100000.33240.6676trialstickalter1200000.32890.6711trialstickalter1300000.33133330.6686667trialstickalter1400000.3314750.668525trialstickalter1500000.332580.66742trialstickalter1600000.33451670.6654833trialstickalter1700000.33251430.6674857trialstickalter1800000.33161250.6683875trialstickalter1900000.33204440.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(iin1:
n)x=floor(runif(1,min=0,max=100)+1);y=floor(runif(1,min=0,max=100)+1);if(x=y)stick=stick+1elsealter=alter+1;stick=stick/n;alter=alter/n;print(data.frame(trial=n,stick=stick,alter=alter);for(iin1:
9)f(i*104),结果为:
trialstickalter1100000.00940.9906trialstickalter1200000.010050.98995trialstickalter1300000.0095333330.9904667trialstickalter1400000.01010.9899trialstickalter1500000.010180.98982trialstickalter1600000.01050.9895trialstickalter1700000.0095428570.9904571trialstickalter1800000.01008750.9899125trialstickalter1900000.0096111110.9903889,当模拟的次数逐渐的增多时,其模拟值越接近理论值,这说明模拟的效果越好。
在随机事件的大量重复出现中,往往呈现几乎必然的规律,这个规律就是大数定律。
通俗地说,这个定理就是,在试验不变的条件下,重复试验多次,随机事件的频率近似于它的概率。
偶然必然中包含着必然。
此次模拟试验也正好用实际的模拟例子说明了大数定理的正确性和应用性。
问题二蒲丰投针问题,1777年法国科学家蒲丰提出的一种计算圆周率的方法-随机投针法,这种方法的具体步骤是:
1、取一张白纸,在上面画上许多条间距为d的平行直线。
2、取一根长为l(ld)的针,随机向画有平行直线的纸上掷n次,观察针与直线相交的次数,记为m。
3、计算针与直线相交的概率,问题分析针与平行线相交的充要条件是:
xl*sin/2,其中0xd/2,0建立直角坐标系(,x),上述条件在坐标系下将是曲线所围成的曲边梯形区域:
由几何概率知:
程序为:
f=function(n)d=2;l=1.5;m=0;for(iin1:
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(iin2:
5)f(10i),结果为:
numestimate11003.125numestimate110002.946955numestimate1100003.111388numestimate11e+053.14327可见,模拟的次数越多,所得的结果就越准确。
这与理论所知的结果是一样的。
问题三掷骰子问题,掷三只相同的色子,求至少有两个点数相同的概率。
问题分析构造3个整数随机变量,使之服从1到6之间的均匀分布,比较每次模拟出的3个值,如果至少有两个相同,则将成功的次数加1,分析可知,成功的频率即为所求的概率。
程序为:
f=function(n)success=0;for(iin1:
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.frame(num=n,p=p);mm=proc.time();for(iin2:
9)f(i*104);proc.time()-mm,结果为:
nump1200000.4447nump1300000.4428333nump1400000.441525nump1500000.44162nump1600000.4456667nump1700000.4445857nump1800000.4463875nump1900000.4451111proc.time()-mm用户系统流逝11.190.0314.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=EminX;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(iin1:
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;elseif(t2+t3t1)jones=jones+1;t=t2+t3;elsesmith=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(iin1:
9)f(3,i*103),结果为:
smithjonesbrownETETreal15060.2320.2620.4931070.5smithjonesbrownETETreal19790.2640.24650.49383580.5smithjonesbrownETETreal114590.25533330.25833330.49128720.5smithjonesbrownETETreal120350.246250.2450.50535070.5smithjonesbrownETETreal125580.24780.24060.50286690.5smithjonesbrownETETreal129830.24433330.25850.49651260.5smithjonesbrownETETreal135380.24757140.2470.50549180.5smithjonesbrownETETreal140190.24950.2481250.49952210.5smithjonesbrownETETreal145380.24166670.25411110.50181040.5,问题五生日问题,有n个人的班级中至少有两个人生日是同一天的概率,模拟之。
问题分析假设:
所有人生日是相互独立的。
设一年有y天,为求至少两个人生日在同一天的概率,我们可以先求n个人生日全不同的概率:
n个人生日所有可能的组合(不考虑次序)种数有:
而n个人生日都不一样的可能的组合方式种数有:
故任意两人生日不同的理论概论是:
(这样写可以方便编程计算)从而有至少两人生日同一天的概率是p=1-p0注:
显然当学生人数n比一年的天数y大时,有抽屉原理易知:
p=1.,程序为:
f=function(n,tn=100,y=365)#n:
班级人数#tn:
实验次数#y:
一年的天数k=0;#记录至少有两人同一天生日的次数for(iin1:
tn)#重复实验次数a=0;#记录n个学生生日的数组for(jin1:
n)#给每个学生一个生日aj=floor(y*runif
(1)+1);j=j+1;b=unique(a);#抽取a中不在同一天的生日if(length(a)!
=length(b)#判断是否有至少两个人在同一天生日k=k+1;p=k/tn;print(data.frame(stu_num=n,tn=tn,p=p);for(iin1:
9)f(50,i*100),结果为:
stu_numtnp1501000.96stu_numtnp1502000.99stu_numtnp1503000.95stu_numtnp1504000.975stu_numtnp1505000.96stu_numtnp1506000.9716667stu_numtnp1507000.9685714stu_numtnp1508000.98375stu_numtnp1509000.9677778,问题六抽奖模型,假设底部的每个桶中装有不同价值的奖品,当你向每个桶中投入一个小球时,需要付出1元代价,一般抽奖者总是希望奖品回报的期望大于所付出的代价,而设置这个抽奖游戏的人很明显是希望赚钱的,即要使奖品回报的期望小于每次投球所付出的代价,游戏设定者会如何设定奖品的价值呢?
假设有一种奖品设定方案如下:
很明显要设定合适的球比例数才能使收益大,同时也要让顾客满意,假使现在每个桶中的比例相同程序为:
f=function(n)x=rbinom(n=6,size=n,prob=rep(1/6,5);price=c(5,1,0.2,0.2,1,5);benefit=sum(x*price);benefit_mean=benefit/n;benefit_real=1;print(data.frame(times=n,benefit_mean=benefit_mean,benefit_real);for(iin1:
9)f(i*104),结果为:
timesbenefit_meanbenefit_real1100002.01511timesbenefit_meanbenefit_real1200002.069491timesbenefit_meanbenefit_real1300002.0662931timesbenefit_meanbenefit_real1400002.0899351timesbenefit_meanbenefit_real1500002.0775961timesbenefit_meanbenefit_real1600002.064611timesbenefit_meanbenefit_real1700002.0735231timesbenefit_meanbenefit_real1800002.0524051timesbenefit_meanbenefit_real1900002.0891911显然这样做是会亏损的,先调整桶中球的比率为:
f=function(n)p=1/100;x=rbinom(n=6,size=n,prob=c(5*p,10*p,30*p,40*p,10*p,5*p);price=c(5,1,0.2,0.2,1,5);benefit=sum(x*price);benefit_mean=benefit/n;benefit_real=1;print(data.frame(times=n,benefit_mean=benefit_mean,benefit_real);for(iin1:
9)f(i*104),结果为:
timesbenefit_meanbenefit_real1100000.848081timesbenefit_meanbenefit_real1200000.837881timesbenefit_meanbenefit_real1300000.83024671timesbenefit_meanbenefit_real1400000.8497951timesbenefit_meanbenefit_real1500000.8432041timesbenefit_meanbenefit_real1600000.838861timesbenefit_meanbenefit_real1700000.84865711timesbenefit_meanbenefit_real1800000.836681timesbenefit_meanbenefit_real1900000.83865331,问题七赶火车问题,一列火车从A站经过B站开往C站,某人每天赶往B站乘这趟火车。
已知火车从A站到B站运行时间为均值30分钟、标准差为2分钟的正态随机变量.火车大约在下午1点离开A站。
离开时刻的频率分布为,该人到达B站时的频率分布为:
用计算机仿真火车开出、火车到达B站、这个人到达B站情况,并给出能赶上火车的仿真结果。
引入以下变量:
T1火车从A站开出的时刻;T2火车从A站运行到B站所需要的时间;T3此人到达B站的时刻;该人能赶上火车的的充分必要条件是T1+T2T3,假设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作为此人赶上火车的概率.当n很大时,频率值与概率值近似相等.,f=function(n)r1=runif(n);r2=runif(n);t2=rnorm(n,30,sd=2);t1=0;t3=0;num=0;for(iin1:
n)if(r1it2k)num=num+1;p=num/n;print(p);for(iin1:
6)f(i*104),结果为:
10.66110.6532510.65910.65522510.6556610.6569833,问题八计算积分问题,#用hitormiss法计算积分值#要积分的函数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;#计算落入制定区域内的频率S=p*c*(b-a);SE=c*(b-a)*sqrt(p*(1-p)/n);#计算相应的积分值和计算估计误差print(S);print(SE);#用样本均值法计算积分值h=function(x)sqrt(1-x2);f=function(h,a,b,n=10)x=runif(n,a,b);y=h(x);S=mean(y)*(b-a);SE=sd(y)/sqrt(n)*(b-a);c(S,SE);for(iin10(1:
5)cat(i,f(h,0,1,i)*4,n);,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 语言 随机 模拟 例子