Matlab筆記——蒙特卡羅方法018

1、18. 蒙特卡羅方法（一）概述一、原理蒙特卡羅（Monte Carlo）方法，是一種基于“隨機(jī)數(shù)”的計(jì)算機(jī)隨機(jī)模擬方法，通過大量隨機(jī)試驗(yàn)，利用概率統(tǒng)計(jì)理論解決問題的一種數(shù)值方法。其理論依據(jù)是：大數(shù)定律、中心極限定理（估計(jì)誤差）。常用來解決如下問題： 1. 求某個(gè)事件的概率，基于“頻率的極限是概率”；2. 可以描述為“隨機(jī)變量的函數(shù)的數(shù)學(xué)期望”的問題，用隨機(jī)變量若干個(gè)具體觀察值的函數(shù)的算術(shù)平均值代替。一般形式為求積分：X為自變量（隨機(jī)變量），定義域?yàn)閍,b, f(x)為被積函數(shù)，為概率密度函數(shù)。蒙特卡羅法步驟為： (1) 依據(jù)概率分布不斷生成N個(gè)隨機(jī)數(shù)x, 依次記為x1, , xN, 并計(jì)算f(

2、xi);(2) 用f(xi)的算術(shù)平均值近似估計(jì)上述積分類比：“積分”同“求和”，“dx”同“1/N”，“”同“服從分布的隨機(jī)數(shù)”；(3) 停止條件：至足夠大的N停止；或者誤差小于某值停止。3. 利用隨機(jī)數(shù)模擬各種分布的隨機(jī)現(xiàn)象，進(jìn)而解決實(shí)際問題。二、優(yōu)缺點(diǎn)優(yōu)點(diǎn)：能夠比較逼真地描述具有隨機(jī)性質(zhì)的事物的特點(diǎn)及物理實(shí)驗(yàn)過程；受幾何條件限制??；收斂速度與問題的維數(shù)無關(guān)；誤差容易確定。缺點(diǎn)：收斂速度慢；誤差具有概率性；進(jìn)行模擬的前提是各輸入變量是相互獨(dú)立的。三、應(yīng)用隨機(jī)模擬實(shí)驗(yàn)，隨機(jī)最優(yōu)化問題，含有大量不確定因素的復(fù)雜決策系統(tǒng)進(jìn)行風(fēng)險(xiǎn)模擬分析（金融產(chǎn)品定價(jià)、期權(quán)）。（二）用蒙特卡羅法求事件概率一、著名

3、的“三門問題”源自博弈論的一個(gè)數(shù)學(xué)游戲：參賽者面前有三扇關(guān)閉的門，其中一扇門的后面藏有一輛汽車，而兩扇門的后面各藏有一只山羊。參賽者從三扇門中隨機(jī)選取一扇，若選中后面有車的那扇門就可以贏得該汽車。當(dāng)參賽者選定了一扇門，但尚未開啟它的時(shí)候，節(jié)目主持人會(huì)從剩下的兩扇門中打開一扇藏有山羊的門，然后問參賽者要不要更換自己的選擇，選取另一扇仍然關(guān)著的門。該問題涉及到的問題是：參賽者更換自己的選擇是否會(huì)增加贏得汽車的概率？解：這是全概率問題。記結(jié)果事件B=“贏得汽車”，造成結(jié)果的原因事件有兩個(gè)：A1=“第一次選到汽車”，A2=“第一次選到山羊”則P(A1)=1/3, P(A2)=2/3; P(B|A1)=

4、0, P(B|A2)=1. 利用全概率公式，P(B)=P(A1)×P(B|A1)+P(A2)×P(B|A2)=2/3而參賽者不更換選擇，抽中汽車的概率為1/3. 可見，參賽者更換選擇可以增加一倍的獲獎(jiǎng)幾率。下面用蒙特卡羅方法來模擬驗(yàn)證上面的理論結(jié)果：將問題“隨機(jī)數(shù)”化，對羊編號為1,2；對汽車編號為3. 先從1,2,3中隨機(jī)選取一個(gè)數(shù)，若開始選中1或2，則更換選擇后選中3, 即贏得汽車；若開始選中3, 則更換選擇后選中1或2，即得不到汽車。這樣的試驗(yàn)重復(fù)n次，記錄開始選中1或2的次數(shù)m（即更換選擇后贏得汽車的次數(shù)），從而可以確定更換選擇后贏得汽車的頻率m/n . 由伯努利大

5、數(shù)定律，當(dāng)試驗(yàn)次數(shù)n增大時(shí)，頻率m/n將趨于更換選擇后贏得汽車的概率。代碼：先編寫SheepAndCar.m函數(shù)function p=SheepAndCar(n) % n可以是正整數(shù)，也可以是正整數(shù)的向量for i=1:length(n) x=randsample(3,n(i),'true'); % 從1,2,3中隨機(jī)抽取n(i)個(gè)數(shù), 'true'表示可以重復(fù)p(i)=sum(x=3)/n(i);end再執(zhí)行代碼：（分別模擬10次、100次、1000000次）p=SheepAndCar(10,100,1000,10000,100000,1000000)運(yùn)行結(jié)果

6、：p = 0.6000 0.7200 0.6910 0.6641 0.6671 0.6666二、用蒲豐投針法求圓周率【蒲豐投針問題】：平面上畫有間隔為d的等距平行線，向平面內(nèi)任意投擲一枚長為h（h<d）的針，求針與任一平行線相交的概率。用Y表示針的中點(diǎn)與最近一條線的距離，用X表示針與此直線間的夾角，則(X, Y)為二維隨機(jī)變量，其樣本空間為平面上的矩形區(qū)域：由于是任意投擲，(X, Y)在上服從均勻分布。記事件A=“針與平行線相交”，則A所對應(yīng)的區(qū)域?yàn)橛蓭缀胃怕手?，事件A的概率為若h, d已知，代入上式就可求出P(A). 反過來，由P(A)的值也可以求. P(A)可以通過蒙特卡羅方法隨機(jī)模

7、擬獲得：在區(qū)域隨機(jī)均勻投點(diǎn)，落在區(qū)域A中的點(diǎn)的頻率fA，隨著投點(diǎn)次數(shù)的增大，就趨于概率P(A). 從而代碼：先編寫B(tài)uffonMonteCarlo.m函數(shù)function p0,pm,piv = BuffonMonteCarlo(d,h,N)%返回p0為理論概率,pm為基于蒙特卡羅法的模擬概率,piv為pi的模擬值. % 輸入?yún)?shù)d為相鄰兩條平行線的間距，h為針的長度，N為模擬投針的次數(shù).if h >= d error('針的長度應(yīng)小于相鄰平行線的間距')endp0 = 2*h/(d*pi); % 計(jì)算針與任一平行線相交的理論概率x = 0;y = 0; m = leng

8、th(N); pm = zeros(1,m); pival = pm; % 賦變量初值for i = 1:mx = unifrnd(0,pi,1,N(i); % 生成0,pi上服從均勻分布的1×N(i)個(gè)隨機(jī)數(shù) y = unifrnd(0,d/2,1,N(i);yb = h*sin(x)/2; pm(i) = sum(y <= yb)/N(i); % 頻率pival(i) = 2*h/(d*pm(i); % 求圓周率的模擬值end再執(zhí)行代碼：p0, pm, pival=BuffonMonteCarlo(10,8,10,100,1000,10000,100000,1000000)

9、運(yùn)行結(jié)果：p0 =0.5093pm = 0.5000 0.5300 0.5020 0.5118 0.5112 0.5093pival = 3.2000 3.0189 3.1873 3.1262 3.1300 3.1416（三）用蒙特卡羅法求積分蒙特卡羅法計(jì)算四重以下的積分，效率可能不如一般數(shù)值積分方法，但蒙特卡羅法與積分重?cái)?shù)無關(guān)，所以在計(jì)算高重積分特別有效。蒙特卡羅法求n重積分的原理：設(shè)D為Rn中的區(qū)域，f：DR，其n重積分表示為先找一個(gè)包含區(qū)域D的矩體C，在C中隨機(jī)均勻地投n個(gè)點(diǎn)，統(tǒng)計(jì)落入D中的點(diǎn)，記為k個(gè)，則區(qū)域D的測度函數(shù)f在D上的均值為由于是均勻投點(diǎn)，故有注：C為矩體，m(C)即為其各

10、邊長的乘積。例1求曲線與直線y=x所圍成的陰影部分的面積。其理論值為取D=C=0,1. 代碼：f=(x) sqrt(x)-x;I0 = quad(f,0,1) % 理論值(數(shù)值解)n=10000;X=unifrnd(0,1,1,n);format longI=sum(f(X)/n運(yùn)行結(jié)果：I0 = 0.166659569272020I = 0.166386907074713例2求球體被圓柱面所截得的立體（含在圓柱面內(nèi)的部分）的體積。如圖為立體在第一卦限的部分，在xoy面的投影區(qū)域?yàn)镈，轉(zhuǎn)化為極坐標(biāo)取C=0,2×0,1,代碼：I0 = 4*quad2d(x,y)sqrt(4-x.2-y

11、.2),0,2,0,(x)sqrt(1-(1-x).2) % 理論值(數(shù)值解)n=100000;x1=unifrnd(0,2,1,n);x2=unifrnd(0,1,1,n);ind = (x2 <= sqrt(2*x1-x1.2);format longI=4*2*sum(sqrt(4-x1(ind).2-x2(ind).2)/n運(yùn)行結(jié)果：I0 = 9.644049708032036I = 9.644614410855425例3 計(jì)算積分已知數(shù)值結(jié)果為179.2969.取C=1,2×1,4×1,16. 代碼：n=100000;x1=unifrnd(1,2,1,n);x2=unifrnd(1,4,1,n);x3=unifrnd(1,16,