蒲豐氏問(wèn)題課件

蒙特卡羅方法編程作業(yè)

用蒲豐投針?lè)ㄔ谟?jì)算機(jī)上計(jì)算π值，取a=4、l=3。分別用理論計(jì)算和計(jì)算機(jī)模擬計(jì)算，求連續(xù)擲兩顆骰子，點(diǎn)數(shù)之和大于6且第一次擲出的點(diǎn)數(shù)大于第二次擲出點(diǎn)數(shù)的概率。蒙特卡羅方法編程作業(yè)用蒲豐投針?lè)ㄔ谟?jì)算機(jī)上計(jì)算π值，取a=1數(shù)值求解過(guò)程分析問(wèn)題建立模型確立算法程序設(shè)計(jì)數(shù)值求解過(guò)程分析問(wèn)題建立模型確立算法程序設(shè)計(jì)2蒲豐氏問(wèn)題其中Ｎ為投計(jì)次數(shù)，n為針與平行線(xiàn)相交次數(shù)。這就是古典概率論中著名的蒲豐氏問(wèn)題。

為了求得圓周率π值，在十九世紀(jì)后期，有很多人作了這樣的試驗(yàn)：將長(zhǎng)為2l的一根針任意投到地面上，用針與一組相間距離為2a（l＜a）的平行線(xiàn)相交的頻率代替概率P，再利用準(zhǔn)確的關(guān)系式：求出π值蒲豐氏問(wèn)題其中Ｎ為投計(jì)次數(shù)，n為針與平3蒲豐氏問(wèn)題的求解模型設(shè)針投到地面上的位置可以用一組參數(shù)（x,θ）來(lái)描述，x為針中心的坐標(biāo)，θ為針與平行線(xiàn)的夾角，如圖所示。任意投針，就是意味著x與θ都是任意取的，但x的范圍限于［0，a］，夾角θ的范圍限于［0，π］。在此情況下，針與平行線(xiàn)相交的數(shù)學(xué)條件是針在平行線(xiàn)間的位置

蒲豐氏問(wèn)題的求解模型設(shè)針投到地面上的位置可以4如何產(chǎn)生任意的（x,θ）？x在［0，a］上任意取值，表示x在［0，a］上是均勻分布的，其分布密度函數(shù)為：類(lèi)似地，θ的分布密度函數(shù)為：因此，產(chǎn)生任意的（x,θ）的過(guò)程就變成了由f1(x)抽樣x及由f2(θ)抽樣θ的過(guò)程了。由此得到：

其中ξ1，ξ2均為（0,1）上均勻分布的隨機(jī)變量。

蒲豐氏問(wèn)題的算法如何產(chǎn)生任意的（x,θ）？x在［05

每次投針試驗(yàn)，實(shí)際上變成在計(jì)算機(jī)上從兩個(gè)均勻分布的隨機(jī)變量中抽樣得到（x,θ），然后定義描述針與平行線(xiàn)相交狀況的隨機(jī)變量s(x,θ)，為如果投針Ｎ次，則是針與平行線(xiàn)相交概率Ｐ的估計(jì)值。事實(shí)上，于是有每次投針試驗(yàn)，實(shí)際上變成在計(jì)算機(jī)上從6Matlab算例clearall;clc;aa=5^15;MM=2^48;x1=5;fprintf('蒙特卡羅方法解蒲豐氏問(wèn)題\n');fprintf('作者：向東2010年3月\n');a=input('請(qǐng)輸入平行線(xiàn)相間的半距離（默認(rèn)值a=4.0）a=');l=input('請(qǐng)輸入針的半長(zhǎng)度（默認(rèn)值l=3.0）l=');N=input('請(qǐng)輸入模擬試驗(yàn)次數(shù)（默認(rèn)值N=100000）N=');

ifisempty(a)a=4.0;endifisempty(l)l=3.0;endifisempty(N)N=100000;end產(chǎn)生偽隨機(jī)數(shù)的乘同余方法為什么不用rand(1)?Matlab算例clearall;clc;ifise7參照：在［a，b］上均勻分布的抽樣在［a，b］上均勻分布的分布函數(shù)為：則其分布密度函數(shù)為：參照：在［a，b］上均勻分布的抽樣8s=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;endends=0;9屏幕輸出結(jié)果：蒙特卡羅方法解蒲豐氏問(wèn)題作者：請(qǐng)輸入平行線(xiàn)相間的半距離（默認(rèn)值a=4.0）a=4請(qǐng)輸入針的半長(zhǎng)度（默認(rèn)值l=3.0）l=3請(qǐng)輸入模擬試驗(yàn)次數(shù)（默認(rèn)值N=100000）N=100000真值pi=3.141593e+000計(jì)算pi=3.141230e+000

p=s*1.0/N;result=2*l/(a*p);fprintf('真值pi=%d\n',pi);fprintf('計(jì)算pi=%d\n',result);屏幕輸出結(jié)果：蒙特卡羅方法解蒲豐氏問(wèn)題p=s*1.0/10randomx1=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;endendrandomx1=1;randomx2=1;s=11參考：散射方位角余弦分布的抽樣

散射方位角φ在［0,2π］上均勻分布，則其正弦和余弦sinφ和cosφ服從如下分布： 直接抽樣方法為：參考：散射方位角余弦分布的抽樣 散射方位角12 令φ=2θ，則θ在［0,π］上均勻分布，作變換 其中0≤ρ≤1，0≤ρ≤π，則 (x,y)表示上半個(gè)單位圓內(nèi)的點(diǎn)。如果(x,y)在上半個(gè)單位圓內(nèi)均勻分布，則θ在［0,π］上均勻分布，由于 令φ=2θ，則θ在［0,π］上均勻分布，作變換13 因此抽樣sinφ和cosφ的問(wèn)題就變成在上半個(gè)單位圓內(nèi)均勻抽樣(x,y)的問(wèn)題。 為獲得上半個(gè)單位圓內(nèi) 的均勻點(diǎn)，采用挑選法，在 上半個(gè)單位圓的外切矩形內(nèi) 均勻投點(diǎn)（如圖）。 舍棄圓外的點(diǎn)，余下的就是所要求的點(diǎn)。 抽樣方法為： 抽樣效率

E=π/4≈0.785> 因此抽樣sinφ和cosφ的問(wèn)題就變成在上14randomx1=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;

end

sinphi=randomx2/sqrt(randomx1^2+randomx2^2);

if(x<=l*sinphi)s=s+1;end替換+挑選抽樣1：>randomx1=1;randomx2=1;15randomx1=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;

end

sinphi=2*randomx1*randomx2/(randomx1^2+randomx2^2);

if(x<=l*sinphi)s=s+1;end>替換+挑選抽樣2：randomx1=1;randomx2=1;16屏幕輸出結(jié)果：p=s*1.0/N;result=2*l/(a*p);fprintf('真值pi=%d\n',pi);fprintf('計(jì)算pi=%d\n',result);蒙特卡羅方法解蒲豐氏問(wèn)題作者：向東2010年3月請(qǐng)輸入平行線(xiàn)相間的半距離（默認(rèn)值a=4.0）a=4請(qǐng)輸入針的半長(zhǎng)度（默認(rèn)值l=3.0）l=3請(qǐng)輸入模擬試驗(yàn)次數(shù)（默認(rèn)值N=100000）N=100000真值pi=3.141593e+000計(jì)算pi=3.149011e+000

屏幕輸出結(jié)果：p=s*1.0/N;蒙特卡羅方法解蒲豐氏問(wèn)17C/C++算例#include<stdio.h>#include<stdlib.h>#include<math.h>#include<iostream.h>#definePI3.1415926535897932384626433832795_int64rdm=5;_int64a=pow(5,15);//30517578125_int64M=pow(2,48);//281474976710656_int64rdm_n=0;//產(chǎn)生隨機(jī)數(shù)doublerandom01(){ _int64rand1; doublerand2;

rand1=(a*rdm)%M; if(rand1<0)rand1=M+rand1; rdm=rand1; rand2=rand1*1.0/M; rdm_n++; returnrand2;}產(chǎn)生偽隨機(jī)數(shù)的乘同余方法注意數(shù)據(jù)類(lèi)型

注意冪運(yùn)算

注意取模運(yùn)算C/C++算例#include<stdio.h>產(chǎn)生偽隨機(jī)18intmain(){ cout<<"蒙特卡羅方法解蒲豐氏問(wèn)題\n"; cout<<"作者：向東2010年3月\n"; doubleaa,ll; doublex,p,phi,sinphi,randomx1,randomx2,result; longN,s,n; cout<<"請(qǐng)輸入平行線(xiàn)相間的半距離（默認(rèn)值a=4.0）a="; cin>>aa; cout<<"請(qǐng)輸入針的半長(zhǎng)度（默認(rèn)值l=3.0）l="; cin>>ll; cout<<"請(qǐng)輸入模擬試驗(yàn)次數(shù)（默認(rèn)值N=100000）N="; cin>>N;

//……

return0;}intmain()19 s=0; for(n=1;n<N;n++) {

x=aa*random01();

randomx1=1; randomx2=1;

while(randomx1*randomx1+randomx2*randomx2>1) { 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("計(jì)算pi=%f\n",result); s=0;20作業(yè)2（選做）分別用理論計(jì)算和計(jì)算機(jī)模擬計(jì)算，求連續(xù)擲兩顆骰子，點(diǎn)數(shù)之和大于6且第一次擲出的點(diǎn)數(shù)大于第二次擲出點(diǎn)數(shù)的概率。蒙特卡羅方法解問(wèn)題作者：請(qǐng)輸入模擬試驗(yàn)次數(shù)（默認(rèn)值N=100000）N=100000點(diǎn)數(shù)之和大于6且第一次擲出的點(diǎn)數(shù)大于第二次擲出點(diǎn)數(shù)的概率=2.503600e-001

作業(yè)2（選做）蒙特卡羅方法解問(wèn)題21參考：擲骰子點(diǎn)數(shù)的抽樣擲骰子點(diǎn)數(shù)X=n的概率為：選取隨機(jī)數(shù)ξ，如