概率实验报告-蒙特卡洛积分

本科实验报告实验名称:《概率与统计》随机模拟实验随机模拟实验实验设随机变量X的分布律为P{X=i}=2-i,i=1,2,3 试产生该分部的随机数1000个，并作出频率直方图。一、实验原理采用直接抽样法：定理：设U是服从[0,1]上的均匀分布的随机变量，则随机变量y=F-1(U)与X有相同的分布函数Y=F-1(U)(为F(x)的逆函数)，即Y=F-1(U)的分部函数为F(x).二、题目分析易得题中X的分布函数为F(x)=1--,i<x<i+1,i=0,1,2,3,……2i若用ceil表示对小数向正无穷方向取整，则F(x)的反函数为F-1(〃)=ceilI叱")I

[产生服从[0,1]上的均匀分布的随机变量a，则m=F-1(a)则为题中需要产生的随机数。三、MATLAB实现f=[];i=1;whilei<=1000a=unifrnd(0,1); %产生随机数a，服从【0,1】上的均匀分布m=log(1-a)/log(1/2);b=ceil(m);%对m向正无穷取整f=[f,b];i=i+1;enddisplay(f);[n,xout]=hist(f);bar(xout,n/1000,1)产生的随机数(取1000个中的20个)如下：

i12345678910f(i)2233124133i11121314151617181920f(i)1221111112频率分布直方图0.60.70.60.60.40.60.70.60.60.4□.3.0.20.10 .2. 46 8 1012实验二设随机变量X的密度函数为，/、f4xe-2x,x>0,

f(x)=<0,x<0试产生该分布的随机数1000个，并作出频率直方图一、实验原理取舍抽样方法，当分布函数的逆函数难以求出时，可采用此方法。取舍抽样算法的流程为：(1)选取一个参考分布，其选取原则，一是该分布的随机样本容易产生；二是存在常数C，使得f(x)<Cg(x)。

产生参考分布g(x)的随机样本x；0独立产生[0,1]产生参考分布g(x)的随机样本x；0独立产生[0,1]上的均匀分布随机数u；0若uCg(x)<f(x)00二、题目分析选取参考分布则保留x0,作为所需的随机样本；否则舍弃。Ie-x,x>0,g(x)=LcI0,x<0则有f(x)<2g(x)，密度函数为g(x)随机变量的产生采用题1中的方法即可三、MATLAB实现f=[];i=1;whilei<=1000a=rand(1); %产生[0,1]上的随机分布的数ab=-log(1-a); %利用直接抽样法产生密度函数为g(x)的随机数bc=rand(1); %产生[0,1]上的随机分布的数aifc*2*(1-a)<=4*b*exp(-2*b)%判断c*2*g(b)是否小于f(b),若满足，则将b作为所需的随机样本保留f=[f,b];i=i+1;endenddisplay(f);[n,xout]=hist(f);bar(xout,n/sum(n),1)产生的部分随机数(取1000个中的20个)如下：i12345678910f(i)0.76420.40641.55131.08330.07680.41271.34141.44290.33720.7056i11121314151617181920f(i)0.56730.17380.97900.93140.21520.19140.24050.26000.88501.2865频率分布直方图

实验三两汽车司机每天早上8:00—8:01驾车到达某十字路口。若到达的时间间隔小于15秒，他们必须停车相互避让。问他们过此十字路口需要停车避让的概率？如果是三辆、四辆汽车呢？此概率与汽车数的关系怎样？有多少汽车时有必要安装红绿灯呢？分析求解：首先我们考虑样本空间。如果把n台车的到达时间看为一个序列｛Xn｝。考虑两台车到达时间分别为了,%+At在同一个时间点内到达的概率limJtAtdx=limtAt=0At-00 At-0则可认为序列｛X｝几乎两两不相等。n若考虑2台车在60秒内，间隔时间超过15秒算不拥堵，设A(2,60,15)为不会拥堵,则：J45dxJ60dx9P(A(2,60,15))=升—1产2=-J60dxJ60dx 16x112x1则由此可知2台车在60秒内，间隔时间超过15秒算不拥堵，不拥堵的概率为2。16若考虑三台车，四台车只需将2台车的积分公式进行拓展便可求出答案。现在考虑n台车在t秒内，间隔时间超过s秒算不拥堵，设A(n,t,s)为事件“不会拥堵“,则：上述积分式Jt—(n—1)sdxJt—(n—2)sd上述积分式Jt—(n—1)sdxJt—(n—2)sdx...Jt—s20 s . x2+stdx...Jtdx1 2 n—10x x1 —20,elsedxn—1dx1+s n,t>(n—1)sx—1dxnJtdxJtdx...JtdxJt1 2 n—10x x x1 n—2 n—1dxnJt—(n—1)sdxJt—(n—2)sdx...Jt—s dx Jt dx1x1+sxn—2+sn—1x +sn—1可用数学归纳法求得：JtdxJt0 1x1dx...Jt2dxJtn—1dxntnJt—(n—1)s Jt—(n—2)s Jt—sJt (t—(n—1)s)nx C/v,^xJ C/v,^x...J C/v,^x J C/v,^xTOC\o"1-5"\h\z1 2 n—1 n n!0 x+s x+s x+s n1 n—2 n—1则有：t>(n—1)selse(1(nt>(n—1)selseP(A(n,t,s))=<1—t—/°，为了验证上述公式的正确性，我们用计算机进行了模拟(模拟107组)。模拟的思路是产生服从给定时间范围内的均匀分布的大量随机数，通过判断条件来计数需要停车避让的点，求得其出现频率，由大数定理可得需要停车避让的概率。随机模拟C++代码如下：#include<cmath>#include<dlgoritlmi>using-namespacestd;iLitTLfk;doubles,t;doublea二口口口口口二，；doublegetrand(doiibley)return(irand()-EtAND_MAX-iand(J)/(double)(EtAHD_MAX-RAND_MAX))*yintmain()Brand(tiir.e(0));cout<<,rn=,r;cin>>n;aout<<,rt=,r;cin>>t;cont<<"s=*';cin>>s;k=_ooooaoo；intent(0);for(intL=L;i<=k;--i)foriintj=L;j<=n;——j)a'j'=getrand.(t);与口r七(己一二ra-L-n);boolf口;for(intj=L;j<n.;——j)i£(a；j-1；-a；j；<5)■-：£=Q;I>reak;；if(fi-一ent：printf「pf二i=W. pow((L.0-(ti-L)"b"L.0/t>,n));printfi,rpl=^.WXX'n",cnt"L.(?/k);return0;设p0为上述公式计算值，p1为计算机模拟值。ntsP0P1260150.56250.562617360150.1250.1249852460150.003906250.0039277101000100.389416120.3893441301000100.000034490.00003780501000050.289310160.2893139由上表观察到P0与p1符合程度较好。可基本认为推导公式正确。两汽车司机需要停车避让的概率为0.5625,三辆汽车时概率为0.125，四辆汽车时0.0039。关于是否安装红绿灯，可根据

1一(n^s、Vt70,else若要使不堵车的概率大于a,则有：1-（n^Vt7故当V1故当V1-7）n<a时，需安装红绿灯进行调控。实验四试采用两种方法，模拟计算下定积分卜sx0.9e-xdx0并分析其误差。要计算的定积分是在无穷区间上。考虑到被积函数收敛速度很快（在x为40时已经近似收敛到0），我们用自适应辛普森公式（目前精度最高的一种差值积分算法）计算J2000x0.9e-xdx,求得其近似实际值作为随机模拟误差分析的依0据,C++程序如下：#include<stdio.h>#include<math.h>#include<stdlib.h>#include<algorithm>#include<iostream>usingnamespacestd;doublex[2100000];doublef(doublex){returnexp(log(x)*0.9-x);}intmain(){intn=1000000;for(inti=1;i<=n;++i)x[i]=i*2000.0/n;doubleans=0;for(inti=1;i<=n/2;++i)ans+=2000.0/(3*n)*(f(x[2*i-2])+4*f(x[2*i-1])+f(x[2*i]));printf("%.8lf\n",ans);return0;}执行后得到积分近似值0.961766.下面采用蒙特卡洛法的思想进行随机模拟计算定积分：(一)平均值法：原理：若％~U(a,b)，则有E(g(%))=fbg(x)f(x)dx=一--fbg(x)dx，因此要a b-aa计算fbg(x)dx，只需求随机变量y=g(x)的平均值，再乘以区间长度即可。a求解：Matlab程序如下：x=unifrnd(a,b,[1n]);fori=1:nm(i)=fun(x(i));endmean(m);p=mean(m)*(b-a).当a=100,b=200,n=106时p=0，故可认为在0到100上取x计算得到的积分值即为0到+8上的积分值；.当a=0,b=100,n=106时重复五次实验得序号12345P0.96740.96670.96820.96970.9643平均值P11-0.96726,|0.96726一0.961766|百分误差8 x100%“0.571241%.11 0.9617663.当a=0,b=100,n=107时重复五次实验得序号12345P0.95820.96410.96380.96190.9621平均值P12-0.96202,|0.96202一0.961766|百分误差8 x100%h0.026410%.12 0.961766分析：可以看出，所投点越多，所求值越接近实际值，误差越小(二)随机投点法原理：设r.v.(XY服从矩形他<%<b,c<y<d｝上的均匀分布，则％~U(a,b),y~U(c,d),且X,Y相互独立。记事件A=｛Y<f(x)｝,则P(A)=P(Y<f(x))=Jbff(x)dydx-Jb(f(%)-c)dx-Jbf(%)dx一c(b一a)-uac a a右侧积分值u即为A出现的频率。由伯努利大数定律，用重复试验中A出现的频率作为u的估计值。将(X,Y)看成矩形内的随机投点，用随机点落在区域Y<f(%)中的频率作为定积分的近似值。求解：(1)为了确定y的取值范围使得投点范围更小，先求f(%)-%0.9e-%在0到100上的最大值：symsx;y=@(x)(-xA0.9*exp(-x))>>c=fminbnd(y,0,100)c=0.9000>>0.9A0.9*exp(-0.9)ans=0.3698得求0到100上的最大值为0.3698(2)在｛100<%<OQ0 <03嵌｝ 上投107个点，计算被f(%)包含的点出现的频率N=10A7;m=0;x=unifrnd(100,200,[1,N]);

y=unifrnd(0,0.3698,[1,N]);m=y<((x.A0.9).*exp(-x));n=sum(m);p=n/N*100*0.3698得p=0故可认为在0到100上取x计算得到的积分值即为0到上的积分值.(3)在｛0<%<100,0<y<0.3698｝上投