蒲丰氏问题

蒙特卡罗方法编程作业1)用蒲丰投针法在计算机上计算π值，取a=4、l=3。2)分别用理论计算和计算机模拟计算，求连续掷两颗骰子，点数之和大于6且第一次掷出的点数大于第二次掷出点数的概率。数值求解过程分析问题建立模型确立算法程序设计蒲丰氏问题其中Ｎ为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。alP2)(22nNalaPl为了求得圆周率π值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（l＜a）的平行线相交的频率代替概率P，再利用准确的关系式：求出π值蒲丰氏问题的求解模型设针投到地面上的位置可以用一组参数（x,θ）来描述，x为针中心的坐标，θ为针与平行线的夹角，如图所示。任意投针，就是意味着x与θ都是任意取的，但x的范围限于［0，a］，夹角θ的范围限于［0，π］。在此情况下，针与平行线相交的数学条件是针在平行线间的位置sinlx如何产生任意的（x,θ）？x在［0，a］上任意取值，表示x在［0，a］上是均匀分布的，其分布密度函数为：类似地，θ的分布密度函数为：因此，产生任意的（x,θ）的过程就变成了由f1(x)抽样x及由f2(θ)抽样θ的过程了。由此得到：其中ξ1，ξ2均为（0,1）上均匀分布的随机变量。其他,00,/1)(1axaxf其他,00,/1)(2f21ax蒲丰氏问题的算法每次投针试验，实际上变成在计算机上从两个均匀分布的随机变量中抽样得到（x,θ），然后定义描述针与平行线相交状况的随机变量s(x,θ)，为如果投针Ｎ次，则是针与平行线相交概率Ｐ的估计值。事实上，于是有其他当,0sin,1),(lxxsNiiiNxsNs1),(1aladxddxdfxfxsPl2)()(),(sin0021NsalaPl22Matlab算例clearall;clc;aa=5^15;MM=2^48;x1=5;fprintf('蒙特卡罗方法解蒲丰氏问题\n');fprintf('作者：向东2010年3月\n');a=input('请输入平行线相间的半距离（默认值a=4.0）a=');l=input('请输入针的半长度（默认值l=3.0）l=');N=input('请输入模拟试验次数（默认值N=100000）N=');ifisempty(a)a=4.0;endifisempty(l)l=3.0;endifisempty(N)N=100000;end)(mod,1Mxaxii,2,1,11iMxii产生伪随机数的乘同余方法为什么不用rand(1)?参照：在［a，b］上均匀分布的抽样在［a，b］上均匀分布的分布函数为：bxbxaabaxaxxF当当当10)()(abaXF则其他,0),/(1)(bxaabxf其分布密度函数为：s=0;forn=1:1:N;x2=mod(aa*x1,MM);x1=x2;randomx=x2*1.0/MM;x=a*randomx;x2=mod(aa*x1,MM);x1=x2;randomx=x2*1.0/MM;phi=pi*randomx;if(x=l*sin(phi))s=s+1;endend其他,00,/1)(1axaxf1ax其他,00,/1)(2f其他当,0sin,1),(lxxs2屏幕输出结果：蒙特卡罗方法解蒲丰氏问题作者：请输入平行线相间的半距离（默认值a=4.0）a=4请输入针的半长度（默认值l=3.0）l=3请输入模拟试验次数（默认值N=100000）N=100000真值pi=3.141593e+000计算pi=3.141230e+000p=s*1.0/N;result=2*l/(a*p);fprintf('真值pi=%d\n',pi);fprintf('计算pi=%d\n',result);NsalaPl22randomx1=1;randomx2=1;while(randomx1^2+randomx2^2)1x2=mod(aa*x1,MM);x1=x2;randomx1=x2*1.0/MM;x2=mod(aa*x1,MM);x1=x2;randomx2=x2*1.0/MM;endsinphi=2*randomx1*randomx2/(randomx1^2+randomx2^2);if(x=l*sinphi)s=s+1;ends=0;forn=1:1:N;x2=mod(aa*x1,MM);x1=x2;randomx=x2*1.0/MM;x=a*randomx;x2=mod(aa*x1,MM);x1=x2;randomx=x2*1.0/MM;phi=pi*randomx;if(x=l*sin(phi))s=s+1;endend参考：散射方位角余弦分布的抽样散射方位角φ在［0,2π］上均匀分布，则其正弦和余弦sinφ和cosφ服从如下分布：直接抽样方法为：其它当011111)(2xxxf2coscos2sinsin令φ=2θ，则θ在［0,π］上均匀分布，作变换其中0≤ρ≤1，0≤ρ≤π，则(x,y)表示上半个单位圆内的点。如果(x,y)在上半个单位圆内均匀分布，则θ在［0,π］上均匀分布，由于sincosyx2222sincosyxyyxx222222222cossin22sinsinsincos2coscosyxxyyxyx因此抽样sinφ和cosφ的问题就变成在上半个单位圆内均匀抽样(x,y)的问题。为获得上半个单位圆内的均匀点，采用挑选法，在上半个单位圆的外切矩形内均匀投点（如图）。舍弃圆外的点，余下的就是所要求的点。抽样方法为：抽样效率E=π/4≈0.78521yx2221212221222122212sin,cos1randomx1=1;randomx2=1;while(randomx1^2+randomx2^2)1x2=mod(aa*x1,MM);x1=x2;randomx1=x2*1.0/MM;randomx1=2*randomx1-1;x2=mod(aa*x1,MM);x1=x2;randomx2=x2*1.0/MM;endsinphi=randomx2/sqrt(randomx1^2+randomx2^2);if(x=l*sinphi)s=s+1;end替换+挑选抽样1：222122221)12(sin1)12(randomx1=1;randomx2=1;while(randomx1^2+randomx2^2)1x2=mod(aa*x1,MM);x1=x2;randomx1=x2*1.0/MM;x2=mod(aa*x1,MM);x1=x2;randomx2=x2*1.0/MM;endsinphi=2*randomx1*randomx2/(randomx1^2+randomx2^2);if(x=l*sinphi)s=s+1;end22212122212sin1替换+挑选抽样2：屏幕输出结果：p=s*1.0/N;result=2*l/(a*p);fprintf('真值pi=%d\n',pi);fprintf('计算pi=%d\n',result);蒙特卡罗方法解蒲丰氏问题作者：向东2010年3月请输入平行线相间的半距离（默认值a=4.0）a=4请输入针的半长度（默认值l=3.0）l=3请输入模拟试验次数（默认值N=100000）N=100000真值pi=3.141593e+000计算pi=3.149011e+000NsalaPl22C/C++算例#includestdio.h#includestdlib.h#includemath.h#includeiostream.h#definePI3.1415926535897932384626433832795_int64rdm=5;_int64a=pow(5,15);//30517578125_int64M=pow(2,48);//281474976710656_int64rdm_n=0;//产生随机数doublerandom01(){_int64rand1;doublerand2;rand1=(a*rdm)%M;if(rand10)rand1=M+rand1;rdm=rand1;rand2=rand1*1.0/M;rdm_n++;returnrand2;})(mod,1Mxaxii,2,1,11iMxii产生伪随机数的乘同余方法注意数据类型注意幂运算注意取模运算intmain(){cout蒙特卡罗方法解蒲丰氏问题\n;cout作者：向东2010年3月\n;doubleaa,ll;doublex,p,phi,sinphi,randomx1,randomx2,result;longN,s,n;cout请输入平行线相间的半距离（默认值a=4.0）a=;cinaa;cout请输入针的半长度（默认值l=3.0）l=;cinll;cout请输入模拟试验次数（默认值N=100000）N=;cinN;//……return0;}s=0;for(n=1;nN;n++){x=aa*random01();randomx1=1;randomx2=1;while(randomx1*randomx1+randomx2*randomx21){randomx1=random01();randomx2=random01();}sinphi=2*randomx1*randomx2/(randomx1*randomx1+randomx2*randomx2);if(x=ll*sinphi)s++;}p=s*1.0/N;result=2*ll/(aa*p);printf(真值pi=%f\n,PI);printf(计算pi=%f\n,result);作业2（选做）分别用理论计算和计算机模拟计算，求连续掷两颗骰子，点数之和大于6且第一次掷出的点数大于第二次掷出点数的概率。蒙特卡罗方法解问题作者：请输入模拟试验次数（默认值N=100000）N=100000点数之和大于6且第一次掷出的点数大于第二次掷出点数的概率=2.503600e-001参考：掷骰子点数的抽样掷骰子点数X=n的概率为：选取随机数ξ，如则在等概率的情况下，可使用如下更简单的方法：其中［］表示取整数。61)(nXP661nnnXF1]6[FX1]6[FXclearall;clc;aa=5^15;MM=2^48;x1=5;fprintf('蒙特卡罗方法解蒲丰氏问题\n');fprintf('作者：向东2010年3月\n');N=input('请输入模拟试验次数（默认值N=100000）N=');ifisempty(N)N=100000;ends=0;forn=1:1:N;x2=mod(aa*x1,MM);x1=x2;randomx=x2*1.0/MM;test1=fix(6*randomx)+1;x2=mod(aa*x1,MM);x1=x2;randomx=x2*1.0/MM;test2=fix(6*randomx)+1;if((test1+test2)6)&&(test1test2)s=s+1;endendp=s*1.0/N;fprintf('点数之和大于6且第一次掷出的点数大于第二次掷出点数的概率=%d\n',p);Matlab算例