@author：https://blog.csdn.net/qq_39559641/article/details/81050845

白袍小夫轉，分類在Graphy,後面再放一個C++And HLSL or GLSL的。

? ?

? ?

本文講述一下蒙特卡洛模擬方法與定積分計算，首先從一個題目開始：設 0≤f(x)≤10≤f(x)≤1 ，用蒙特卡洛模擬法求定積分 J=∫10f(x)dxJ=∫01f(x)dx 的值。

隨機投點法

設 (X,Y)(X,Y) 服從正方形 {0≤x≤1,0≤y≤1}{0≤x≤1,0≤y≤1} 上的均勻分布，則可知 X,YX,Y 分別服從 [0,1] 上的均勻分布，且 X,YX,Y 相互獨立。記事件 A={Y≤f(X)}A={Y≤f(X)} ，則 AA 的概率為

P(A)=P(Y≤f(X))=∫10∫f(x)0dydx=∫10f(x)dx=JP(A)=P(Y≤f(X))=∫01∫0f(x)dydx=∫01f(x)dx=J

即定積分 JJ 的值 就是事件 AA 出現的頻率。同時，由伯努利大數定律，我們可以用重復試驗中 AA 出現的頻率作為 pp 的估計值。即將 (X,Y)(X,Y) 看成是正方形 {0≤x≤1,0≤y≤1}{0≤x≤1,0≤y≤1} 內的隨機投點，用隨機點落在區域 y≤f(x)y≤f(x) 中的頻率作為定積分的近似值。這種方法就叫隨機投點法，具體做法如下：

1、首先產生服從 [0,1][0,1] 上的均勻分布的 2n2n 個隨機數（ nn 為隨機投點個數，可以取很大，如 n=104n=104 ）並將其配對。

2、對這 nn 對數據 (xi,yi),i=1,2,...;,n(xi,yi),i=1,2,...;,n ，記錄滿足不等式 yi≤f(xi)yi≤f(xi) 的個數，這就是事件 AA 發生的頻數 μnμn ，由此可得事件 AA 發生的頻率 μnnμnn ，則 J≈μnnJ≈μnn 。

舉一實例，譬如要計算 ∫10e?x2/2/√2πdx∫01e?x2/2/2πdx ，模擬次數 n=104n=104 時，R 代碼如下：

n=10^4;
x=runif(n);
y=runif(n);
f=function(x)
{
exp(-x^2/2)/sqrt(2*pi)
}
mu_n=sum(y<f(x));

模擬次數 n=105n=105 時，令 n=105n=105 , 其余不變。

定積分 ∫10e?x2/2/√2πdx∫01e?x2/2/2πdx 的精確值和模擬值如下：

表 1

註：精確值用 integrate(f,0,1) 求得

擴展

如果你很細心，你會發現這個方法目前只適用於積分區間 [0,1][0,1] ，且積分函數 f(x)f(x) 在區間 [0,1][0,1] 上的取值也位於 [0,1][0,1] 內的情況。那麽，對於一般區間 [a,b][a,b] 上的定積分 J′=∫bag(x)dxJ′=∫abg(x)dx 呢？一個很明顯的思路，如果我們可以將 J′J′ 與 JJ 建立代數關系就可以了。

首先，做線性變換，令 y=(x?a)/(b?a)y=(x?a)/(b?a) ，此時，

x=(b?a)y+ax=(b?a)y+a , J′=(b?a)∫10g[(b?a)y+a]dyJ′=(b?a)∫01g[(b?a)y+a]dy 。

進一步如果在區間 [a,b][a,b] 上有 c≤g(x)≤dc≤g(x)≤d ，令

f(y)=1d?cg(x)?c=1d?cg[a+(b?a)y]?cf(y)=1d?cg(x)?c=1d?cg[a+(b?a)y]?c ，

則 0≤f(y)≤10≤f(y)≤1 。此時，可以得到\J‘=\int_{a}^{b}g(x)dx=S_0J+c(b-a) \)。

其中，S0=(b?a)(d?c)S0=(b?a)(d?c) , J=∫10f(y)dyJ=∫01f(y)dy 。

這說明，用隨機投點法計算定積分方法具有普遍意義。

舉一個實例，求定積分 J′=∫52e?x2/2/√2πdxJ′=∫25e?x2/2/2πdx 。

顯然 a=2,b=5a=2,b=5 ，c=min{g(x)|2≤x≤5},d=max{g(x)|2≤x≤5}c=min{g(x)|2≤x≤5},d=max{g(x)|2≤x≤5} ，由於 g(x)=e?x2/2/√2πg(x)=e?x2/2/2π 在 [2,5][2,5] 上時單調減函數，所以 c=g(5),d=g(2)c=g(5),d=g(2) ，S0=(b?a)(d?c)S0=(b?a)(d?c) 。R 中代碼為

a=2;
b=5;
g=function(x)
{
exp(-x^2/2)/sqrt(2*pi);
}
f=function(y)
{
(g(a+(b-a)*y)-c)/(d-c);
}
c=g(5);d=g(2);s_0=(b-a)*(d-c);
n=10^4;
x=runif(n);y=runif(n);
mu_n=sum(y<=f(x));
J=mu_n/n;
J_0=s_0*J+c*(b-a);

定積分 J′=∫52e?x2/2/√2πdxJ′=∫25e?x2/2/2πdx 的精確值和模擬值如下：

表 2

註：精確值用 integrate(g,2,5) 求得

平均值法

蒙特卡洛模擬法計算定積分時還有另一種方法，叫平均值法。這個原理也很簡單：設隨機變量 XX 服從 [0,1][0,1] 上的均勻分布，則 Y=f(X)Y=f(X) 的數學期望為

E(f(X))=∫10f(x)dx=JE(f(X))=∫01f(x)dx=J

所以估計 JJ 的值就是估計 f(X)f(X) 的數學期望值。由辛欽大數定律，可以用 f(X)f(X) 的觀察值的均值取估計 f(X)f(X) 的數學期望。具體做法：

先用計算機產生 nn 個服從 [0,1][0,1] 上均勻分布的隨機數：xi,i=1,2,...;,nxi,i=1,2,...;,n 。

對每一個 xixi ，計算 f(xi)f(xi) 。

計算 J≈1n∑ni=1f(xi)J≈1n∑i=1nf(xi) 。

譬如，計算 J=∫10e?x2/2/√2πdxJ=∫01e?x2/2/2πdx ，R 中的代碼為

n=10^4;

x=runif(n);
f=function(x)
{
exp(-x^2/2)/sqrt(2*pi)
}
J=mean(f(x));

其精確值和模擬值如下：

表 3

平均值法與隨機投點法類似可以進行擴展，這裏不再贅述。

結論

用蒙特卡洛模擬法計算定積分具有普遍意義。蒙特卡洛模擬方法為我們提供了一個看待世界的新思路，即一個不具隨機性的事件可以通過一定的方法用隨機事件來模擬或逼近。

? ?

來自 <https://cosx.org/2010/03/monte-carlo-method-to-compute-integration>

（轉）蒙特卡洛方法與定積分計算