蒙特卡羅（Monte Carlo）方法——從數(shù)學原理到實際案例

根據(jù)實際問題的特點，構(gòu)造簡單而又便于實現(xiàn)的概率統(tǒng)計模型，使所求的解恰好是所求問題的概率分布或數(shù)學期望。對于本身就具有隨機性質(zhì)的問題，主要是描述和模擬其概率過程；對于本身不具備隨即性質(zhì)的確定性問題（比如計算定積分），就必須先構(gòu)造一個具有隨機性質(zhì)概率過程，并使它的參數(shù)正好是待求問題的解（比如通過投飛鏢試驗來求解不規(guī)則圖形的面積） ，即將不具有隨機性質(zhì)的問題轉(zhuǎn)化為隨機性質(zhì)的問題。

確定性數(shù)學問題：首先建立一個與待求解有關(guān)的概率模型，使待求解就是所建立模型的概率分布或數(shù)學期望。然后對這個模型進行隨機抽樣觀察，即產(chǎn)生隨機變量。最后用算術(shù)平均值作為待求解的近似估計值。比如，求積分、矩陣求逆、解方程（組）、偏微分方程邊值問題、計算微分方程特征值等。

隨機性問題：一般采用直接模擬，根據(jù)實際物理情況的概率法則，用計算機進行抽樣試驗。

②利用概率分布進行抽樣（進行隨機試驗）

通過計算機產(chǎn)生已知概率分布的隨機變量。常用的概率分布包括均勻分布、正態(tài)分布、指數(shù)分布、泊松分布、卡方分布等。這一步是Monte Carlo方法的核心，因為它涉及到如何準確地模擬隨機事件。

在上一步構(gòu)造了概率模型后，由于各種概率模型都是由各種概率分布構(gòu)成的，因此產(chǎn)生已知概率分布的隨機變量就成了實現(xiàn)Monte Carlo方法模擬試驗的基本手段，這也是Monte Carlo方法被稱為隨機抽樣的原因。實現(xiàn)Monte Carlo模擬的基本工具是隨機數(shù)，以均勻分布為例，隨機數(shù)可以看成具有均勻分布的隨機變量，隨機數(shù)序列就是具有均勻分布的總體的一個簡單子樣（服從均勻分布的相互獨立的隨機變量序列）。因此，在概率分布上抽樣的過程實際上就是產(chǎn)生隨機數(shù)的過程。

隨機數(shù)的產(chǎn)生有物理方法和數(shù)學方法。在計算機上通過物理方法產(chǎn)生的是真正的隨機數(shù)，但是代價較高。我們一般使用數(shù)學遞推公式來產(chǎn)生偽隨機數(shù)序列。大量試驗已經(jīng)證明，偽隨機數(shù)與隨機數(shù)具有相近的性質(zhì)，因此可以作為真正的隨機數(shù)來使用。

③建立估計量

構(gòu)造隨機概率模型并從中抽樣后，需要確定一個隨機變量作為所求問題的解，一般稱之為無偏估計。通常，使用多次隨機抽樣結(jié)果的算術(shù)平均值作為解的近似值。

④重復實驗分析結(jié)果

為減少估計方差并提高準確性，需要進行大量重復試驗。最后對試驗結(jié)果進行分析，計算估計值的置信區(qū)間和誤差，評估解的可靠性。

Monte Carlo方法的模擬手段

模擬是Monte Carlo方法的重要手段，一般通過物理或者數(shù)學模型來近似類比，模擬現(xiàn)實系統(tǒng)的演變過程，來尋找其規(guī)律。模擬的基本思想是建立一個試驗?zāi)Ｐ?，這個模型包含了所研究系統(tǒng)的主要特點，通過運行該模型來獲取需要的信息。模擬的實現(xiàn)方法有物理模擬和數(shù)學模擬。

物理模擬：通過與實際系統(tǒng)相近的實物系統(tǒng)模擬。比如，導彈發(fā)射試驗、拋擲硬幣試驗等。物理模擬通常代價較大，難以重復進行，甚至無法實現(xiàn)。

數(shù)學模擬：數(shù)學模擬一般通過計算機實現(xiàn)，也叫做計算機模擬。計算機模擬可以反復進行，但是在實際問題中，建模需要進行種種簡化和假設(shè)。

Monte Carlo方法的收斂性與誤差

Monte Carlo方法作為一種計算方法，有必要討論其收斂性以及誤差。

收斂性

蒙特卡羅方法的收斂性可以由大數(shù)定律給出。

蒙特卡羅方法是由隨機變量X的簡單子樣X1,X2,...,Xn的算術(shù)平均值來作為所求解的近似值。

由大數(shù)定律可知，如果X1,X2,...,Xn獨立同分布，且具有有限期望值（E(X)<∞），那么

也就是說，當子樣數(shù)N足夠大時，隨機變量X的算術(shù)平均值以概率1收斂到它的數(shù)學期望值。蒙特卡羅方法的收斂是概率意義上的收斂，雖然不能斷言蒙特卡羅方法的誤差不會超過某個值，但是能指出它的誤差以接近1的概率不超過某個界限。

誤差

蒙特卡羅方法的近似值與真值的誤差可以由中心極限定理給出。

中心極限定理指出，如果X1,X2,...,Xn獨立同分布，并且具有有限非零的方差

其中，f(x)是X的密度函數(shù)。那么

當N充分大時，有如下的近似表達

我們把α叫做置信度，1-α叫做置信水平。該近似表明，不等式

近似地以概率1-α成立，且誤差收斂速度的階為

一般將蒙特克羅方法的誤差ε定義為

誤差ε表達式中的正態(tài)差λ與置信度α是一一對應(yīng)的，根據(jù)問題要求確定出置信水平之后，通過查閱正態(tài)分布表，就可以確定出λ的值。下面給出常用的λ與α的數(shù)值對照表。

當α=0.5時，誤差超過ε的概率和誤差小于ε的概率相等，即α=1-α=0.5，此時稱ε為概然誤差。

蒙特卡羅方法的誤差是一種概率誤差，這和其他數(shù)值計算方法是有區(qū)別的。并且誤差表達式中的σ是未知的，必須使用它的估計值來代替，估計值計算方法如下

效率

通過誤差表達式可以看出，在給定置信度α后，誤差ε由σ和N決定。想要減小誤差ε，要么增大試驗次數(shù)N，要么減小方差σ2。并且，在標準差σ（標準差也叫均方差）固定的情況下，要想把精度提升一個數(shù)量級，試驗次數(shù)要增加兩個數(shù)量級（因為N和ε是平方根倒數(shù)的關(guān)系）。如果通過減小均方差σ，例如σ減小一半，那么誤差ε也能減小一半，這相當于試驗次數(shù)N會增大4倍。

一般來說，降低方差往往會使得觀察一個子樣的時間增加。在固定時間內(nèi)，使觀察的樣本數(shù)減少。所以，該方法的優(yōu)劣，需要由方差和觀察一個子樣的費用（計算機耗費的時間）綜合來衡量。這就是蒙特卡羅方法中效率的概念，效率定義為 σ2·c ，其中c是觀察一個子樣的平均費用， σ2·c 越小，效率越高。假設(shè)在固定誤差ε，和產(chǎn)生一個隨機變量X的平均費用c不變的條件下，如果均方差σ縮小為原來的1/10，那么試驗次數(shù)N可以減小為原來的1/100。我們將誤差表達式變形可得，N=(λσ/ε)2 ，那么總的費用 N·c=(λσ/ε)2·c=(λ/ε)2·σ2c。由此可見要想提高蒙特卡羅方法的效率，應(yīng)該盡可能的降低 σ2c 的值。

由誤差定義可知，在給定置信水平的情況下，蒙特卡羅方法的收斂速度為 O(N^-0.5)，這個值和問題本身的維數(shù)無關(guān)，問題本身維數(shù)的變化，只會引起抽樣時間以及估計量計算時間的變化，不會對誤差造成影響。也就是說，使用蒙特卡羅方法時，抽取的子樣總數(shù)N與維數(shù)s無關(guān)，維數(shù)的增加只會引起計算量的增加，不會影響誤差。這個特點使得蒙特卡羅方法對解決多維問題有很好的適應(yīng)性。而一般的數(shù)值計算方法，比如計算定積分，計算時間會隨著維數(shù)的冪次方而增加，并且由于分點數(shù)與維數(shù)的冪次方成正比，需要占用一定的存儲空間。

大數(shù)定律和中心極限定理

在完美世界中，我們可以準確的回答任何問題。但是在現(xiàn)實世界中，要精確的計算出某些概率是非常困難甚至不可能的。不過，對于許多現(xiàn)實問題，一般不需要知道確切答案，只要得到一定范圍內(nèi)的答案就足夠了。

極限定理是概率論中最重要的理論結(jié)果，而極限定理中最重要的便是大數(shù)定律和中心極限定理。通常，當隨機變量序列的平均值在某種條件下收斂到數(shù)學期望（均值）的時候，就是我們所說的大數(shù)定律。當大量隨機變量之和的分布在某種條件下逼近于正態(tài)分布時，就稱之為中心極限定理。

大數(shù)定律

大數(shù)定律也叫做大數(shù)法則，通俗來講，就是在樣本數(shù)量很大的時候，樣本均值和真實值充分接近。大數(shù)定律與中心極限定理結(jié)合在一起將成為現(xiàn)代概率統(tǒng)計科學的基石。

切比雪夫不等式

首先介紹兩個不等式

我們知道，如果知道了隨機變量的概率分布，我們可以直接計算出概率值，但是有時候概率分布是未知的。此時，馬爾可夫不等式和切比雪夫不等式表明：當我們只知道隨機變量的期望，或者同時知道期望和方差時，可以導出概率的上界。

收斂類型

大數(shù)定律有強大數(shù)定律和弱大數(shù)定律兩種類型，二者主要考察的是獨立隨機變量之和的性質(zhì)，而在很多情況下，這些性質(zhì)都與極限有關(guān)。所以，有必要給出以下幾種收斂類型。

在后面，我們會經(jīng)常用到依概率收斂的概念，實際上依概率收斂的含義是：隨機序列距X的偏差不小于任一給定值的可能性會隨著n的增大而越來越小。

弱大數(shù)定律（辛欽大數(shù)定律）

如果在一個固定分布中進行獨立抽樣，那么樣本均值任意接近于總體均值的概率趨近于1。簡言之，樣本均值依概率收斂于期望（弱大數(shù)定律的收斂是一種依概率收斂）。所以弱大數(shù)定律也可以表示為下面形式。

強大數(shù)定律

強大數(shù)定律表明，獨立同分布的隨機變量序列，前n個觀察值的平均值以概率為1地收斂到分布的平均值。強大數(shù)定律中的收斂是指幾乎必然收斂。

強大數(shù)定律和弱大數(shù)定律的區(qū)別：

弱大數(shù)律只能保證對充分大的 n * ，隨機變量(X1+…+Xn * )/n * 靠近u。但它不能保證對一切n>n*，(X1+·+Xn)/n也一定在u的附近。這樣，|(X1+·+Xn)/n-u|可以無限多次離開0（盡管出現(xiàn)較大偏離的頻率不會很高)。而強大數(shù)律能保證這種情況不會出現(xiàn)，強大數(shù)律能夠以概率為1地保證，

只可能出現(xiàn)有限次。

伯努利大數(shù)定律

伯努利大數(shù)定律表明， 事件發(fā)生的頻率依概率收斂于事件發(fā)生的概率。該定理以嚴格的數(shù)學形式表達了頻率的穩(wěn)定性，也就是說當n很大時，事件發(fā)生的頻率與總體概率有較大偏差的可能性很小。（概率是頻率的穩(wěn)定值）。

比如最簡單的伯努利試驗拋硬幣，理論上正面朝上和反面朝上的概率都是0.5。依據(jù)伯努利大數(shù)定律，我們重復的拋擲一枚均勻的硬幣，當拋擲次數(shù)n足夠大時，正面和反面出現(xiàn)的次數(shù)因該是一樣的，都是0.5n，并且次數(shù)n越大，這個值越接近。

中心極限定理

中心極限定理用粗略的語言來說，大量的獨立隨機變量之和的分布近似為正態(tài)分布。中心極限定理為計算獨立隨機變量和的有關(guān)概率提供了理論依據(jù)，同時也解釋了現(xiàn)實世界中許多實際分布的頻率曲線為鐘形曲線（正態(tài)概率密度）的原因。中心極限定理揭示了測量誤差近似正態(tài)分布這一重要的統(tǒng)計規(guī)律，因此也被稱為誤差頻率定律。

中心極限定理的陳述

首先來看什么是正態(tài)分布。正態(tài)分布是概率中最重要的分布，正態(tài)分布也稱為高斯分布和鐘形曲線。

中心極限定理有很多版本，各版本的區(qū)別在于假設(shè)條件和最終收斂類型的不同。下面給出幾種常見形式的中心極限定理。

中心極限定理為統(tǒng)計學提供了正態(tài)分布的廣泛應(yīng)用可能性，因為許多統(tǒng)計方法都是基于正態(tài)分布的性質(zhì)來推導的。在實際應(yīng)用中，中心極限定理使得我們能夠通過對樣本數(shù)據(jù)的分析來推斷總體的參數(shù)，即使我們不知道總體的具體分布形態(tài)。這對于數(shù)據(jù)分析和決策制定有著重要的實際意義。中心極限定理還解釋了為什么在實際中可以使用樣本均值來估計總體均值，并且可以通過樣本均值的分布來估計總體均值的可靠性。中心極限定理表明，在一定條件下，大量獨立隨機變量的和或平均值會趨向于正態(tài)分布，無論這些隨機變量本身服從何種分布。這一發(fā)現(xiàn)極大地簡化了概率計算和統(tǒng)計分析的過程。

中心極限定理與高斯白噪聲

在使用KF、EKF、UKF進行濾波估計時，我們一般會假設(shè)傳感器的觀測噪聲為高斯白噪聲，并且在很多場合都會假設(shè)噪聲是服從正態(tài)分布的。那么為什么要假設(shè)噪聲為正態(tài)分布，這種假設(shè)的理論依據(jù)又是什么呢？

首先，高斯白噪聲在不同時刻下的觀測值相互獨立。在實際系統(tǒng)中，一般存在眾多噪聲源（比如通過傳感器測量角度時，影響測量結(jié)果的因素有很多），實際上這些噪聲源大多滿足相互獨立的假設(shè)（比如聲吶每次測量角度時，測量結(jié)果都是相互獨立的），當噪聲源的數(shù)量足夠多且每個噪聲源對于總體的貢獻可忽略不計時，根據(jù)中心極限定理可知，這些噪聲源的累加結(jié)果服從正態(tài)分布（無論這些噪聲本身服從什么分布，只要他們相互獨立，那么他們的累加和也就是最終的噪聲就是服從正態(tài)分布的）。 簡單來說，根據(jù)中心極限定理，多個獨立隨機變量的和趨近于高斯分布，無論這些變量本身的分布如何。在實際環(huán)境中，噪聲往往是多種不同來源的隨機變量的疊加，因此其整體分布趨向于高斯分布。

中心極限定理的應(yīng)用

在微積分中，我們可以對一個函數(shù)在某個區(qū)間上積分，并得到積分值。但實際上，求積分的前提是，可以找到一個給定函數(shù)的原函數(shù)的完美解析表達式。在很多時候，這是一件非常困難甚至不可能實現(xiàn)的事。在很多時候為了計算積分，需要借助無窮級數(shù)展開（誤差函數(shù)）來實現(xiàn)。于是，就有了中心極限定理最重要的應(yīng)用之一：蒙特卡羅積分。

假設(shè)有一個面積為R的正方形，在正方形內(nèi)有一個不規(guī)則的區(qū)域A，我們現(xiàn)在要求這個區(qū)域A的面積。如果我們可以用函數(shù)來較好的描述區(qū)域A，那么很自然的我們會想到通過對函數(shù)積分來求A的面積，此時需要知道原函數(shù)，即：若A可以描述為

那么A的面積可以通過以下積分來求

如果函數(shù)f和g的原函數(shù)難以得到，那么積分法就很難實現(xiàn)了。這時就可以采用蒙特卡羅法，其可以描述為：在正方形R內(nèi)隨機取點，那么A的面積應(yīng)該接近于落在A中的點數(shù)占總共抽取的點數(shù)的比例乘以R。以投飛鏢為例，向R中投擲N次飛鏢，那么

這就是蒙特卡羅法。再比如下面圖形求陰影面積，使用蒙特卡羅將是一個很好的選擇。

隨機數(shù)的產(chǎn)生和隨機數(shù)的意義

在前面說過，隨機數(shù)（RNs）是Monte Carlo的基石，隨機試驗是Monte Carlo的實現(xiàn)核心。最初，隨機抽樣試驗的形式是投擲飛鏢、擲硬幣、投骰子等，并觀察理論概率和頻率。隨著計算機技術(shù)的出現(xiàn)和發(fā)展，這種多次重復進行的隨機抽樣試驗特別適合計算機來執(zhí)行。而隨機數(shù)的產(chǎn)生也經(jīng)歷了大量的研究，比如各種隨機抽樣數(shù)字表，但是這種隨機數(shù)表中的隨機數(shù)總是有限的，而有時候我們要做的抽樣試驗次數(shù)要遠大于隨機數(shù)表中給出的隨機數(shù)個數(shù)。利用物理現(xiàn)象產(chǎn)生的隨機數(shù)是完全隨機的，比如RAND公司曾經(jīng)以隨機脈沖源為信息源，用電子旋轉(zhuǎn)輪產(chǎn)生隨機數(shù)表，但是物理方法的代價太過昂貴。

通過數(shù)學方法，不斷迭代，可以產(chǎn)生一系列偽隨機數(shù)，但是顯然偽隨機數(shù)并不是真正隨機的。不過，在迭代過程開始之前，每一項都是不可預測的，對產(chǎn)生的偽隨機數(shù)，只要可以通過一系列局部隨機性檢驗（比如均勻性、獨立性檢驗等）就可以作為隨機數(shù)使用，因此通過數(shù)學方法產(chǎn)生的隨機數(shù)稱為偽隨機數(shù)。我們所熟知的matlab中的rand()函數(shù)就是用于產(chǎn)生0~1之間服從均勻分布的偽隨機數(shù)。

隨機數(shù)的意義

下面依然通過蒙特卡羅積分的例子來說明隨機數(shù)的重要意義。

假設(shè)我們現(xiàn)在要估計一個積分值

首先，我們引入一個在(0,1)上服從均勻分布的隨機變量x并構(gòu)造一個新的隨機變量 y=g(x) 。進一步，通過

我們把待求的未知量 I （積分值）表示成了一個隨機變量y的期望值。假設(shè)隨機變量x代表某個具體實驗中的物理量，那么，我們便可以用期望的頻率解釋去估計 I ，通過多次重復試驗，然后觀察每次試驗x的取值xi，并計算相應(yīng)的yi值 yi=g(xi) ，然后再求平均值，以此來求出 I 的一個近似值。

我們將具有某些性質(zhì)的數(shù)據(jù)xi稱為隨機數(shù)，這樣，只要我們能夠產(chǎn)生這些隨機數(shù)，就可以確定出待求積分值 I 。這就是隨機數(shù)的意義。

但是，實際上“什么是隨機數(shù)？隨機數(shù)如何產(chǎn)生？能否得到真正的隨機數(shù)？”這一系列問題并沒有一個標準的答案。因為隨機數(shù)就像概率一樣，在理論上和經(jīng)驗上的意義是完全不同的。在理論上，隨機數(shù)是按照一個抽象模型定義的抽象概念抽象概念。在經(jīng)驗上，隨機數(shù)是一個實數(shù)序列，它要么來自一個具體實驗的物理數(shù)據(jù)，要么是一個確定的計算機程序的輸出。

富蘭克林（Franklin）對隨機數(shù)的定義：如果一個序列具有一個均勻分布的隨機變量獨立樣本構(gòu)成的所有無窮序列所具有的所有性質(zhì)，稱這個序列是隨機的。

萊莫（Lehmer）對隨機數(shù)的定義：隨機序列是一個模糊的概念，包含了這樣一些思想，即序列的每一項對非初始項來說是不可預測的，同時其數(shù)字通過了一定數(shù)目的檢驗，這一概念與統(tǒng)計學家的習慣一致，并在某種程度上取決于這些序列的用途。

下面是從實際用途考慮給出的隨機數(shù)的兩種定義。

①概念上的定義：如果一個序列 xi 是在重復試驗空間上獨立同分布的隨機變量序列 Xi 的樣本，即 xi=Xi(ξ) ，則被稱為是隨機的。

這個定義和富蘭克林的定義類似。區(qū)別在于，在富蘭克林定義中，數(shù)列 xi 具有獨立同分布隨機變量所具有的所有性質(zhì)，而在本定義中，xi 等于獨立隨機變量序列 Xi 的樣本。

②經(jīng)驗上的定義：如果一個數(shù)列 xi 的統(tǒng)計性質(zhì)與從隨機試驗中得到的隨機數(shù)據(jù)是相同的，則被稱為是隨機的。

隨機數(shù)的產(chǎn)生

在蒙特卡羅計算中，隨機數(shù)主要通過計算機程序產(chǎn)生，當然，隨機數(shù)也可以是隨機試驗產(chǎn)生的觀測數(shù)據(jù)。比如，投擲一個骰子可以產(chǎn)生1-6的隨機序列，放射性物質(zhì)粒子發(fā)射之間的時間間隔生成的指數(shù)分布的觀測樣本。

目前計算機并不能直接產(chǎn)生任意分布的隨機數(shù)序列，已有的算法只能生成服從 (0,m) 上均勻分布的隨機整數(shù)序列，（但是任意分布的序列可以在此基礎(chǔ)上間接得到）。

下面是一個最簡單的產(chǎn)生隨機數(shù)序列的算法

其中，f 是一個由最近的已得到的 r 個數(shù)確定的函數(shù)，z 是 f 除以 m 的余數(shù)。序列 z 的生成是非線性的，遞歸關(guān)系由整數(shù) m ，函數(shù) f ，以及初值 z1確定。這樣設(shè)計的隨機數(shù)發(fā)生器所生成的隨機數(shù)的質(zhì)量依賴于函數(shù) f 的形式，f 越復雜，產(chǎn)生的隨機數(shù)質(zhì)量越高。

萊莫算法：一種最簡單的遞歸產(chǎn)生隨機數(shù)的算法

其中，m 是一個很大的素數(shù)（萊莫建議大小為2^31-1），a 是一個整數(shù)。

隨機數(shù)的產(chǎn)生算法很多，我們不需要去深究，因為現(xiàn)在的編程語言大都已經(jīng)封裝好了產(chǎn)生隨機數(shù)的函數(shù)，我們只需要去調(diào)用即可。并且也不建議自己去生成隨機數(shù)序列（確切的說是偽隨機數(shù)）。

Monte Carlo方法實例

投硬幣模擬試驗

硬幣投擲試驗是最經(jīng)典的隨機試驗之一，下面通過計算機程序來模擬硬幣投擲試驗。隨機投擲一枚均勻的硬幣，正面朝上的次數(shù)應(yīng)該是服從參數(shù)為 (1,p) 的二項分布，p 為正面朝上的概率。

matlab中產(chǎn)生二項分布偽隨機數(shù)的函數(shù)為 r=binornd(n,p) ，參數(shù)與返回值如下：

n - Number of trials? ? ? ? ? ? positive integer | array of positive integersp - Probability of success for each trial? ? ? ? ? ? scalar value | array of scalar values

r - Random numbers from binomial distribution? ? ? ? ? ? scalar value | array of scalar values

matlab試驗代碼如下

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %?file:	硬幣投擲模擬試驗% %?author:	mindtechnist% %% %?log:% %?% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%end-of-header%%%%%%%%%%%%%%%%%%%%%%%%%%%%close all; clear; clc;N =?2000; % 仿真試驗次數(shù)count =?0; % 結(jié)果為正面的次數(shù)freq = zeros(1,N); % 頻率prob =?0.5*ones(1,N); % 概率for?i =?1:N? ? count = count + binornd(1,?0.5); % 均勻硬幣，正反面出現(xiàn)的概率都是0.5? ? freq(i) = count/i;endfigurehold on;plot(1:N, freq,?'-r');plot(1:N, prob,?'-k');legend("頻率",?"概率");

仿真結(jié)果如下

該試驗證明了，頻率穩(wěn)定于概率。

古典概率試驗

在一個袋子中有10個完全相同的球，分別標號為1~10。每次任取一球，記錄號碼后將球放回，然后繼續(xù)取球。即有放回的抽取?，F(xiàn)在有放回的抽取3個球，則3個球的號碼均為偶數(shù)的概率是多少？

通過數(shù)學計算可以得到理論概率為0.125。下面做500次重復的仿真試驗，并取其均值。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % file:? ?古典概率模擬試驗
% % author:?mindtechnist
% %
% % log:
% %?
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%end-of-header%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close?all;?clear;?clc;

N?=?500; ?

M?=?1000;

freq?=?zeros(1,N);

for?k?=?1:N
? ??freq(k) =?subfun(M);
end

freq_mean?=?mean(freq);

figure
hold?on;?box?on;
plot(freq,?'-r');
line([0,N],[freq_mean,freq_mean],'LineWidth',2,'Color','k');
legend("估計概率",?"均值");

function?p?=?subfun(M)

? ??count?=?0;?% 全為偶數(shù)的次數(shù)
? ? ?
? ??for?k=1:M
? ? ?
? ? ? ??ball1?=?unidrnd(10);?% r = unidrnd(n) 從由最大值 n 指定的離散均勻分布中生成隨機數(shù)。
? ? ? ??ball2?=?unidrnd(10);
? ? ? ??ball3?=?unidrnd(10);
? ? ?
? ? ? ??if( (mod(ball1,2)==0)?&&?(mod(ball2,2)==0)?&&?(mod(ball3,2)==0) )
? ? ? ? ? ??count?=?count?+?1;
? ? ? ??end
? ??end
? ? ?
? ??p?=?count/M;?% 計算頻率

end

仿真結(jié)果如下

蒙特卡羅方法求圓周率

如下圖所示，假設(shè)圓的半徑為1，正方形邊長為2，均勻的產(chǎn)生隨機變量點 (x,y) ，x和y的取值都在[-1,1]內(nèi)，并以圓心為原點。共產(chǎn)生N個隨機點，計算距離原點不大于1的點的個數(shù)（即落在圓內(nèi)的點的個數(shù)）n，那么圓周率的值近似為4n/N。( πr2/4 = n/N )

python代碼如下

import?numpy?as?np

print('投點次數(shù)：')
n?=?eval(input())

x?=?np.random.uniform(-1,?1,?n)?# 均勻分布
y?=?np.random.uniform(-1,?1,?n)

d?=?np.sqrt(x**2?+?y**2)
res?=?sum(np.where(d?<=?1.0,?1,?0))

pi?=?4.0?*?res?/?n

print("圓周率的估計值為：")
print(round(pi,?6))

運行結(jié)果

蒙特卡羅法計算定積分

蒙特卡羅方法在計算多重積分的問題中，有著良好的應(yīng)用。因為在多重積分運算中，蒙特卡羅方法的誤差與積分重數(shù)無關(guān)。盡管蒙特卡羅方法的精度較低，但是它能夠很快的給出一個低精度的模擬結(jié)果。

隨機投點法

假設(shè)待求的定積分表達式，以及函數(shù)圖像如下所示

求函數(shù) f(x) 的定積分實際上就是求圖中陰影部分的面積。假設(shè)矩形面積為 S=M(b-a) ，待求陰影部分面積為 θ ，現(xiàn)在我們向矩形區(qū)域隨機投點，那么落在陰影部分的概率可以通過下面的數(shù)學公式計算出來

這就是隨機投點法求定積分的原理。

回到上面的定積分計算問題，我們可以通過數(shù)學方法求出這個定積分的解為7.2432.下面我們編寫matlab程序，使用隨機投點法求解該定積分。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % file:? ?隨機投點法計算定積分
% % author:?mindtechnist
% %
% % log:
% %?
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%end-of-header%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close?all;?clear;?clc;

N?= [1000,?10000,?100000,?1000000];?% 投點總數(shù)
value?=?zeros(1,length(N));

for?i?=?1:length(N)
? ??value(i) =?simulation(N(i));
end

% figure
% hold on;
% plot(value);
value

function?theta?=?simulation(N)

? ??% 積分上下限
? ??a?=?0;?b?=?4;?
? ??
? ??M?=?4;?% f(x)的上界（矩形的上界）
? ??
? ??count?=?0;?% 落在陰影部分的點數(shù)
? ??
? ??for?i?=?1:N
? ? ??
? ? ? ??% 隨機投點
? ? ? ??x?=?unifrnd(a,b);
? ? ? ??y?=?unifrnd(0,M);
? ? ?
? ? ? ??if?(cos(x)+2.0)?>=?y?% y <= f(x) 說明落入陰影區(qū)
? ? ? ? ? ??count?=?count+1;?
? ? ? ??end
? ??end
? ? ?
? ??freq?=?count/N;
? ? ?
? ??theta?=?freq*(b-a)*M;

end

通過仿真結(jié)果可以看到，隨著仿真次數(shù)的增加，結(jié)果越來越接近數(shù)學上算出來精確解。

平均值法

假設(shè)要計算的積分形式為

matlab代碼如下

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % file:? ?平均值法計算定積分
% % author:?mindtechnist
% %
% % log:
% %?
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%end-of-header%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close?all;?clear;?clc;

% 積分上下限
a?=?0;?b?=?1;

N?=?100000;?

value?=?zeros(1,10);

for?k?=?1:10
? ??value(k) =?experiment(a,b,N);
end

value

function?result?=?experiment(a,?b,?N)

? ??count?=?0;
? ??xrandnum?=?unifrnd(a,b,1,N);

? ??for?k?=?1:N
? ? ? ??count?=?count?+?exp(xrandnum(1,k));
? ??end? ??

? ??result?=?count/N;

end

運行結(jié)果如下

計算定積分是蒙特卡羅方法的重要應(yīng)用之一，特別是在復雜積分運算中，通過蒙特卡羅方法模擬是一個有效的解決手段。

關(guān)注公眾號即可獲取python項目實戰(zhàn)視頻資源，快來關(guān)注試試吧！