概率實(shí)驗(yàn)報(bào)告蒙特卡洛積分

1、本 科 實(shí) 驗(yàn) 報(bào) 告實(shí)驗(yàn)名稱： 概率與統(tǒng)計(jì)隨機(jī)模擬實(shí)驗(yàn) 隨機(jī)模擬實(shí)驗(yàn)實(shí)驗(yàn)一設(shè)隨機(jī)變量X的分布律為試產(chǎn)生該分部的隨機(jī)數(shù)1000個(gè)，并作出頻率直方圖。一、實(shí)驗(yàn)原理采用直接抽樣法：定理：設(shè)U是服從0,1上的均勻分布的隨機(jī)變量，則隨機(jī)變量 與X有相同的分布函數(shù)（為F(x)的逆函數(shù)），即的分部函數(shù)為.二、題目分析易得題中X的分布函數(shù)為 若用ceil表示對(duì)小數(shù)向正無(wú)窮方向取整，則F(x)的反函數(shù)為產(chǎn)生服從0,1上的均勻分布的隨機(jī)變量a，則m=F-1(a)則為題中需要產(chǎn)生的隨機(jī)數(shù)。三、MATLAB實(shí)現(xiàn)f=;i=1;while i<=1000 a=unifrnd(0,1); %產(chǎn)生隨機(jī)數(shù)a，服從【0,

2、1】上的均勻分布 m=log(1-a)/log(1/2); b=ceil(m); %對(duì)m向正無(wú)窮取整 f=f,b; i=i+1;enddisplay(f);n,xout=hist(f);bar(xout,n/1000,1)產(chǎn)生的隨機(jī)數(shù)（取1000個(gè)中的20個(gè)）如下：i12345678910f(i)2233124133i11121314151617181920f(i)1221111112頻率分布直方圖實(shí)驗(yàn)二設(shè)隨機(jī)變量X的密度函數(shù)為試產(chǎn)生該分布的隨機(jī)數(shù)1000個(gè)，并作出頻率直方圖一、實(shí)驗(yàn)原理取舍抽樣方法，當(dāng)分布函數(shù)的逆函數(shù)難以求出時(shí)，可采用此方法。取舍抽樣算法的流程為：（1） 選取一個(gè)參考分布，其

3、選取原則，一是該分布的隨機(jī)樣本容易產(chǎn)生；二是存在常數(shù)，使得。（2） 產(chǎn)生參考分布的隨機(jī)樣本；（3） 獨(dú)立產(chǎn)生0,1上的均勻分布隨機(jī)數(shù)；（4） 若，則保留x0,作為所需的隨機(jī)樣本；否則舍棄。二、題目分析選取參考分布則有，密度函數(shù)為g(x)隨機(jī)變量的產(chǎn)生采用題1中的方法即可三、MATLAB實(shí)現(xiàn)f=;i=1;while i<=1000 a=rand(1); %產(chǎn)生0,1上的隨機(jī)分布的數(shù)a b=-log(1-a); %利用直接抽樣法產(chǎn)生密度函數(shù)為g(x)的隨機(jī)數(shù)b c=rand(1); %產(chǎn)生0,1上的隨機(jī)分布的數(shù)a if c*2*(1-a)<=4*b*exp(-2*b) %判斷c*2*g

4、(b)是否小于f(b),若滿足，則將b作 為所需的隨機(jī)樣本保留 f=f,b; i=i+1; endenddisplay(f);n,xout=hist(f);bar(xout,n/sum(n),1)產(chǎn)生的部分隨機(jī)數(shù)（取1000個(gè)中的20個(gè)）如下：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頻率分布直方圖實(shí)驗(yàn)三兩汽車司機(jī)每天

5、早上8:008:01駕車到達(dá)某十字路口。若到達(dá)的時(shí)間間隔小于15秒，他們必須停車相互避讓。問(wèn)他們過(guò)此十字路口需要停車避讓的概率？如果是三輛、四輛汽車呢？此概率與汽車數(shù)的關(guān)系怎樣？有多少汽車時(shí)有必要安裝紅綠燈呢？分析求解：首先我們考慮樣本空間。如果把n臺(tái)車的到達(dá)時(shí)間看為一個(gè)序列Xn?？紤]兩臺(tái)車到達(dá)時(shí)間分別為在同一個(gè)時(shí)間點(diǎn)內(nèi)到達(dá)的概率則可認(rèn)為序列幾乎兩兩不相等。若考慮2臺(tái)車在60秒內(nèi)，間隔時(shí)間超過(guò)15秒算不擁堵，設(shè)為不會(huì)擁堵,則：則由此可知2臺(tái)車在60秒內(nèi)，間隔時(shí)間超過(guò)15秒算不擁堵，不擁堵的概率為。若考慮三臺(tái)車，四臺(tái)車只需將2臺(tái)車的積分公式進(jìn)行拓展便可求出答案?，F(xiàn)在考慮n臺(tái)車在t秒內(nèi)，間隔時(shí)間超

6、過(guò)s秒算不擁堵，設(shè)為事件“不會(huì)擁堵“,則：上述積分式可用數(shù)學(xué)歸納法求得：則有：為了驗(yàn)證上述公式的正確性，我們用計(jì)算機(jī)進(jìn)行了模擬（模擬107組）。模擬的思路是產(chǎn)生服從給定時(shí)間范圍內(nèi)的均勻分布的大量隨機(jī)數(shù)，通過(guò)判斷條件來(lái)計(jì)數(shù)需要停車避讓的點(diǎn)，求得其出現(xiàn)頻率，由大數(shù)定理可得需要停車避讓的概率。隨機(jī)模擬C+代碼如下：設(shè)p0為上述公式計(jì)算值，p1為計(jì)算機(jī)模擬值。ntsP0P1260150.56250.562617360150.1250.1249852460150.003906250.0039277101000100.389416120.3893441301000100.000034490.0000378

7、0501000050.289310160.2893139由上表觀察到P0與p1符合程度較好?？苫菊J(rèn)為推導(dǎo)公式正確。兩汽車司機(jī)需要停車避讓的概率為0.5625,三輛汽車時(shí)概率為0.125，四輛汽車時(shí)0.0039。關(guān)于是否安裝紅綠燈，可根據(jù)若要使不堵車的概率大于a,則有：故當(dāng)時(shí)，需安裝紅綠燈進(jìn)行調(diào)控。實(shí)驗(yàn)四試采用兩種方法，模擬計(jì)算下定積分并分析其誤差。分析求解：要計(jì)算的定積分是在無(wú)窮區(qū)間上?？紤]到被積函數(shù)收斂速度很快（在x為40時(shí)已經(jīng)近似收斂到0），我們用自適應(yīng)辛普森公式（目前精度最高的一種差值積分算法）計(jì)算，求得其近似實(shí)際值作為隨機(jī)模擬誤差分析的依據(jù),C+程序如下：#include<st

8、dio.h>#include<math.h>#include<stdlib.h>#include<algorithm>#include<iostream>using namespace std;double x2100000;double f(double x) return exp(log(x)*0.9-x);int main() int n=1000000; for(int i=1;i<=n;+i) xi=i*2000.0/n; double ans=0; for(int i=1;i<=n/2;+i) ans+=2000.0

9、/(3*n)*(f(x2*i-2)+4*f(x2*i-1)+f(x2*i); printf("%.8lfn",ans); return 0;執(zhí)行后得到積分近似值0.961766.下面采用蒙特卡洛法的思想進(jìn)行隨機(jī)模擬計(jì)算定積分：（一）平均值法：原理：若，則有，因此要計(jì)算，只需求隨機(jī)變量的平均值，再乘以區(qū)間長(zhǎng)度即可。求解：Matlab程序如下：x=unifrnd(a,b,1 n);for i=1:n m(i)=fun(x(i);endmean(m);p=mean(m)*(b-a)1.當(dāng)時(shí)p=0，故可認(rèn)為在0到100上取x計(jì)算得到的積分值即為0到上的積分值;2.當(dāng)時(shí)重復(fù)五次實(shí)驗(yàn)得

10、序號(hào)12345P0.96740.96670.96820.96970.9643平均值,百分誤差.3.當(dāng)時(shí)重復(fù)五次實(shí)驗(yàn)得序號(hào)12345p0.95820.96410.96380.96190.9621平均值,百分誤差.分析：可以看出，所投點(diǎn)越多，所求值越接近實(shí)際值，誤差越?。ǘ╇S機(jī)投點(diǎn)法原理：設(shè)服從矩形上的均勻分布,則，且, 相互獨(dú)立。記事件則右側(cè)積分值u即為A出現(xiàn)的頻率。由伯努利大數(shù)定律，用重復(fù)試驗(yàn)中A出現(xiàn)的頻率作為u的估計(jì)值。將看成矩形內(nèi)的隨機(jī)投點(diǎn)，用隨機(jī)點(diǎn)落在區(qū)域中的頻率作為定積分的近似值。求解：（1）為了確定y的取值范圍使得投點(diǎn)范圍更小，先求在0到100上的最大值：syms x;y=(x)(-x0.9*exp(-x)>> c=fminbnd(y,0,100)c = 0.9000>> 0.90.9*exp(-0.9)ans = 0.3698得求0到100上的最大值為0.3698（2）在上投107個(gè)點(diǎn)，計(jì)算被包含的點(diǎn)出現(xiàn)的頻率N = 107;m=0;x=unifrnd(100,200,1,N);y=unifrnd(0,0.3698,1,N);m=y<(x.0.9).*exp(-x);n=sum(m);p=n/N*100*0.3698得p = 0故可認(rèn)為在0到100上取x計(jì)算得到的積分值即為0到上的積分值.（3）在上投107個(gè)點(diǎn)，計(jì)算被包含的點(diǎn)出現(xiàn)的