數(shù)學建?！队嬎銠C隨機模擬》

計算機隨機模擬(MonteCarlo)（建模培訓）張洪波主講有些問題，由于隨機因素很多，用概率論的方法進行求解可能很難和復雜，這時就需要借助隨機模擬方法得到近似解答。隨機模擬法也叫蒙特卡羅(MonteCarlo)方法，也稱為計算機隨機模擬方法。由于MonteCarlo法計算量大，精度不高，因而需要借助計算機，并僅適合一些用解析方法或常規(guī)數(shù)值方法難以解決問題的低精度求解和驗證一、簡介命令說明rand(m,n)[0,1]上均勻分布unifrnd(a,b,m,n)[a,b]上均勻分布unifrnd(N,m,n)1，…,N的等概率分布randn(m,n)標準正態(tài)分布N(0,1)exprnd(λ,m,n)均值為λ的指數(shù)分布poissrnd(λ,m,n)均值為λ的泊松分布normrnd(μ,σ,m,n)正態(tài)分布N(μ,σ2)二、Matlab中的部分隨機數(shù)產(chǎn)生命令注：以下都是產(chǎn)生不同分布m×n階隨機矩陣正態(tài)分布變量X的數(shù)學期望，方差

2

，密度函數(shù)

計算密度值函數(shù):normpdf(x，mu，sigma)

累積分布函數(shù)，即積分上限函數(shù)

計算概率值函數(shù):normcdf(x，mu，sigma)

逆累積分布函數(shù)值，即已知概率值p，求z使得

計算函數(shù):z=norminv(p，mu，sigma)

有關(guān)正態(tài)分布matlab計算函數(shù)直方圖繪圖函數(shù):hist(data，n)

其中,data是需要處理的數(shù)據(jù)塊,

繪圖原理:利用data中最小數(shù)和最大數(shù)構(gòu)成一區(qū)間,將區(qū)間等分為n個小區(qū)間，統(tǒng)計落入每個小區(qū)間的數(shù)據(jù)量,以數(shù)據(jù)量為高度繪小矩形，形成直方圖。如果省略參數(shù)n，MATLAB將n的默認值取為10。直方圖統(tǒng)計計算函數(shù)：N=hist(data，n)

計算結(jié)果N是n個數(shù)的一維數(shù)組，分別表示data中各個小區(qū)間的數(shù)據(jù)量。這種方式只計算而不繪圖。其他用法：hist(data,x)或N＝hist(data,x)統(tǒng)計直方圖matlab函數(shù)hist三、一些隨機模擬的例子[例1](相遇問題)甲、乙兩船在24小時內(nèi)獨立地隨機到達碼頭.設兩船到達碼頭時刻都服從[0,24]上的均勻分布，如果甲船到達碼頭后停留2小時,乙船到達碼頭后停留1小時.問兩船相遇的概率有多大？（可用幾何概率，此處略）分析：如果知道甲、乙兩船到達的時刻x和y，兩船能相遇的條件就是：兩船到達的時刻x和y可以隨機生成，生成一組到達時刻，可以確定是否能相遇。如果重復很多次，統(tǒng)計相遇的比例就可近似為相遇的概率。這也是MonteCarlo思想.兩船到達碼頭時刻服從[0,24]上的均勻分布，甲船停留2小時,乙船停留1小時,相遇概率？方法一利用軟件產(chǎn)生一組[0，24]上的隨機數(shù)xi:

10.45.30.39.020.16.87.8……5.4

分別代表甲船到達碼頭時間；再產(chǎn)生一組[0，24]上的隨機數(shù)yi:

6.03.48.117.50.813.4……14.0

分別代表乙船到達碼頭時間；把（xi，yi）(i＝1，2…)當作一天中甲船乙船到達的時間，統(tǒng)計出能相遇的頻數(shù)，計算出相遇的頻率做為相遇的概率?？梢允褂密浖╡xcel、C等實現(xiàn)）兩船到達碼頭時刻服從[0,24]上的均勻分布，甲船停留2小時,乙船停留1小時,相遇概率？方法二clc;N=1000;c=0;

%模擬次數(shù)N，相遇次數(shù)c清零fori=1:N

%重復N次到達時間

x=unifrnd(0,24,1,1);

%甲船到達時間x(隨機數(shù))

y=unifrnd(0,24,1,1);

%乙船到達時間y(隨機數(shù))

if((x<=y&y<=x+2)|(y<=x&x<=y+1))c=c+1;

%如果能相遇，則計數(shù)器加1

endendP=c/N

%顯示相遇的概率近似值注意：每次運行的結(jié)果一般都不一樣%計算機模擬程序（報童訣竅的簡化版）報童每天清晨從報社購進報紙零售，晚上退回沒有賣掉的報紙.若每份報紙的購進價為b=0.75元，售出價為a=1元，退回價為c=0.6元.每天需求量X是離散型隨機變量，其分布為[例2]X500510520P0.340.360.30問：如果報童每天購進報紙為n=510份，每天的平均利潤是多少？方法二：如果我們知道每天的需求量，可直接計算利潤。而每天需求量可以按分布生成（隨機模擬思想）方法一：概率方法（略）報童訣竅的簡化版售出價a=1，購進價b=0.75，退回價c=0.6，購進數(shù)量n=510份需求量X500510520概率P0.340.360.30按照需求量的分布規(guī)律，隨機生成N=20個數(shù)據(jù):510520500510520500500500500510500500500520510520510510代表20天的需求量，計算出報童在這20天的總利潤和平均利潤，用平均利潤來近似報童的平均收入。這也是MonteCarlo方法.問題如何按分布規(guī)律產(chǎn)生隨機數(shù)據(jù)？隨機數(shù)據(jù)很多時，如何編程？報童訣竅的簡化版需求量X500510520概率P0.340.360.30問題1：如何產(chǎn)生以下分布規(guī)律的隨機數(shù)據(jù)？注:rand(m,n)可以生成[0,1]上均勻分布隨機數(shù)把[0,1]分成長度為0.34、0.36、0.30的三個區(qū)間[0,0.34]、(0.34,0.70]、(0.70,1]用rand(1,1)產(chǎn)生1個[0,1]上均勻分布隨機數(shù)，

如該數(shù)在[0,0.34]、(0.34,0.70]或(0.70,1]內(nèi)，相當于該天的需求量相應為500、510和520重復多次就可以若干天的需求量了生成N=20天的需求量的matlab代碼可以為報童訣竅的簡化版生成N=20天的需求量的matlab代碼（定義函數(shù)）functiony=randfun1(N)y=zeros(1,N);fori=1:Nt=rand(1,1);

ift<=0.34

X=500;

elseift<=0.70

X=510;

else

X=520;

endy(i)=X;end根據(jù)隨機數(shù)t的范圍，確定需求量X值，并保存到數(shù)組的相應位置中（關(guān)鍵部分）產(chǎn)生1行N列的全0矩陣，目的分配好矩陣大小，可以省略文件名為randfun1.m報童訣竅的簡化版問題2：如何從需求量計算利潤？售出價a=1，購進價b=0.75，退回價c=0.6，購進數(shù)n=510函數(shù)文件名fun2.m

functiony=fun2(x)a=1;%售出價

b=0.75;%購進價

c=0.6;%退回價

n=510;%購進數(shù)量if(x>n)

y=n*(a-b);

else

y=x*(a-b)-(n-x)*(b-c)

end模擬程序代碼N=1000;x=randfun1(N);y=0;fori=1:Ny=y+fun2(x(i));endy/N報童訣竅的簡化版a=1;b=0.75;c=0.6;n=510;

n=510;N=1000;y=0;fori=1:Nt=rand(1,1);

if(t<=0.34)x=500;

elseif(t<=0.70)x=510;

elsex=520;

end

if(x>n)y=y+n*(a-b);

elsey=y+x*(a-b)-(n-x)*(b-c);

endEndfprintf(1,'平均利潤=%.3f',y/N);也可完整模擬程序：根據(jù)隨機數(shù)t，計算需求量x值根據(jù)需求量x，計算利潤并累加到y(tǒng)中。顯示平均利潤售出價a

購進價b退回價c

購進數(shù)n總利潤y模擬天數(shù)N報童訣竅的簡化版粒子分離器的某關(guān)鍵參數(shù)（記作y）由7個零件的參數(shù)（記作x1，x2，…，x7）決定，[例3]均值=標定值標準差=容差等級*標定值÷3關(guān)鍵參數(shù)y的目標值是1.50，當偏離為±0.1時，產(chǎn)品為次品，損失為1000元；當偏離為±0.3時，產(chǎn)品為廢品，損失為9000元。由于工藝原因，7個零件參數(shù)可以看著是正態(tài)隨機變量，在后面的標定值及容差等級情況下，求產(chǎn)品的平均損失？容差等級A=1%B=5%C=15%零件參數(shù)均值=標定值標準差=容差等級*標定值÷3容差等級A=1%B=5%C=15%零件損失Q是一個隨機變量，平均損失就是期望E(Q)由于y的表達式很復雜，要想計算y的分布和上述概率很困難，我們必須尋找較為有效的近似方法。模擬：通過產(chǎn)生指定分布的隨機數(shù)，來代表7個零件的參數(shù)值，計算y值，確定損失大小。多做幾次可得均值零件參數(shù)%編寫計算y值的函數(shù)：y=funY(x1,x2,x3,x4,x5,x6,x7)

%編寫計算損失費函數(shù):Q=funQ(y)D1=[0.1,0.3,0.1,0.1,1.5,16,0.75]%標定值

D2=[5,5,5,15,15,5,5,]./100%容差等級a=D1

b=D1.*D2./3x1=normrnd(a(1),b(1),1,1);

……

x7=normrnd(a(7),b(7),1,1)y=funY(x1,x2,x3,x4,x5,x6,x7)

Q=fu