蒲丰氏问题ppt课件

1、蒙特卡罗方法编程作业,用蒲丰投针法在计算机上计算值，取a=4、l=3。 分别用理论计算和计算机模拟计算，求连续掷两颗骰子，点数之和大于6且第一次掷出的点数大于第二次掷出点数的概率,数值求解过程,分析问题,建立模型,确立算法,程序设计,蒲丰氏问题,其中为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题,为了求得圆周率值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（ la）的平行线相交的频率代替概率P，再利用准确的关系式,求出值,蒲丰氏问题的求解模型,设针投到地面上的位置可以用一组参数（x,）来描述，x为针中心的坐标，为针

2、与平行线的夹角，如图所示。 任意投针，就是意味着x与都是任意取的，但x的范围限于0，a，夹角的范围限于0，。在此情况下，针与平行线相交的数学条件是,针在平行线间的位置,如何产生任意的（x,）？x在0，a上任意取值，表示x在0，a上是均匀分布的，其分布密度函数为： 类似地，的分布密度函数为： 因此，产生任意的（x,）的过程就变成了由f1(x)抽样x及由f2()抽样的过程了。由此得到： 其中1，2均为（0,1）上均匀分布的随机变量,蒲丰氏问题的算法,每次投针试验，实际上变成在计算机上从两个均匀分布的随机变量中抽样得到（x,），然后定义描述针与平行线相交状况的随机变量s(x,)，为 如果投针次，则

3、是针与平行线相交概率的估计值。事实上， 于是有,Matlab 算例,clear all; clc; aa=515; MM=248; 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,if isempty(a) a=4.0; end if isempty(l) l=3.0; end if isempty(N) N=100000; en

4、d,产生伪随机数的乘同余方法,为什么不用rand(1),参照： 在a，b上均匀分布的抽样,在a，b上均匀分布的分布函数为,则,其分布密度函数为,s=0; for n=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; end end,屏幕输出结果,蒙特卡罗方法解蒲丰氏问题 作者： 请输入平行线相间的半距离（默认值a=4.0） a=4 请输入针的半长度（默

5、认值l=3.0） l=3 请输入模拟试验次数（默认值N=100000） N=100000 真值pi = 3.141593e+000 计算pi = 3.141230e+000,p=s*1.0/N; result=2*l/(a*p); fprintf(真值pi = %d n,pi); fprintf(计算pi = %d n,result,randomx1=1; randomx2=1; while (randomx12+randomx22)1 x2=mod(aa*x1,MM); x1=x2; randomx1=x2*1.0/MM; x2=mod(aa*x1,MM); x1=x2; randomx2=

6、x2*1.0/MM; end sinphi=2*randomx1*randomx2/(randomx12+randomx22); if (x=l*sinphi) s=s+1; end,s=0; for n=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; end end,参考：散射方位角余弦分布的抽样,散射方位角在0,2上均匀分布，则其正弦和余弦sin

7、和cos服从如下分布： 直接抽样方法为,令=2，则在0,上均匀分布，作变换 其中01，0，则 (x,y) 表示上半个单位圆内的点。如果 (x,y) 在上半个单位圆内均匀分布，则在0,上均匀分布，由于,因此抽样sin和cos的问题就变成在上半个单位圆内均匀抽样 (x,y) 的问题。 为获得上半个单位圆内 的均匀点，采用挑选法，在 上半个单位圆的外切矩形内 均匀投点（如图）。 舍弃圆外的点，余下的就是所要求的点。 抽样方法为： 抽样效率 E=/40.785,randomx1=1; randomx2=1; while (randomx12+randomx22)1 x2=mod(aa*x1,MM);

8、x1=x2; randomx1=x2*1.0/MM; randomx1=2*randomx1-1; x2=mod(aa*x1,MM); x1=x2; randomx2=x2*1.0/MM; end,sinphi=randomx2/sqrt(randomx12+randomx22); if (x=l*sinphi) s=s+1; end,替换+挑选抽样1,randomx1=1; randomx2=1; while (randomx12+randomx22)1 x2=mod(aa*x1,MM); x1=x2; randomx1=x2*1.0/MM; x2=mod(aa*x1,MM); x1=x2;

9、 randomx2=x2*1.0/MM; end,sinphi=2*randomx1*randomx2/(randomx12+randomx22); if (x=l*sinphi) s=s+1; end,替换+挑选抽样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

10、=100000 真值pi = 3.141593e+000 计算pi = 3.149011e+000,C/C+算例,include #include #include #include #define PI 3.1415926535897932384626433832795 _int64 rdm = 5; _int64 a = pow(5,15);/30517578125 _int64 M = pow(2,48);/281474976710656 _int64 rdm_n = 0; /产生随机数 double random01() _int64 rand1; double rand2; rand

11、1 = (a*rdm)%M; if (rand10) rand1=M+rand1; rdm = rand1; rand2 = rand1*1.0/M; rdm_n+; return rand2;,产生伪随机数的乘同余方法,注意数据类型注意幂运算注意取模运算,int main() coutaa; coutll; coutN; / return 0;,s=0; for(n=1;n1) randomx1=random01(); randomx2=random01(); sinphi=2*randomx1*randomx2/ (randomx1*randomx1+randomx2*randomx2);

12、 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的概率为： 选取随机数，如 则 在等概率的情况下，可使用如下更简单的方法： 其中表示取整数,clear all;clc; aa=515; MM=248; x1=5; fprintf(蒙特卡罗方法解蒲丰氏问题 n); fprintf(作者：向东 2010年3月 n); N=input(请输入模拟试验次数（默认值N=100000） N=); if isempty(N