數值積分

Monte Carlo Method

先前我們處理的數值積分,函數 \(f(x)\) 從 \(R^1\) 映到 \(R^1\),如果我們的定義域不是一維的而是高維的,例如 \(R^3\),我們會使用重積分的概念來切割我們的定義域,這裡出現一個有趣的問題。若定義域的維度很高,例如10,假若我們在一個維度上將定義域切割成100個片段,在整個定義域上就是 \(100^10\) 個小片段。這個量會隨著step size 的減少而快速的增加,這樣會大幅地降低計算的速度。Monte Carlo Method 是一個透過隨機方式來估計積分值的方法,這個方法的收斂速度很慢,但是在高維度時,卻展現他的好處。

什麼是Monte Carlo Method 呢?假設我們要計算函數 \(f(x)\) 在 \([0,1]\)區間上的積分,而我們使用 \([0,1]\) 區間上的均勻分配取亂數的方式,在 \([0,1]\) 區間上取 \(N\) 個點 \(x_i\),然後我們把對應的 \(f(x_i)\) 取其平均。這個平均值,其實就是 \(\int_0^1 f(x) dx\) 的一個估計值。真的嗎?我們看看以下程式。

 

>>> import random as rn
>>> x = [rn.uniform(0,1) for i in xrang(100)]
>>> x = np.array(x)
>>> def myfun(x):
return x**2

>>> x = np.array(x)
>>> y = myfun(x)
>>> sum(y)/len(x)
0.29306117602929921
>>>
>>> x = [rn.uniform(0,1) for i in xrange(1000)]
>>> x = np.array(x)
>>> y = myfun(x)
>>> sum(y)/len(x)
0.33414800736424988

我們知道 \(\int_0^1 x^2 dx = 1/3\),我們透過上述的亂數法估計值為 0.293,若我們將點數增加到1000個點,估計值為0.334148。這個方法很暴力,但是也不賴對吧。簡單說,當我們會因為維度過高而造成計算速度過慢時,我們就在定義域裡均勻選取取值點,然後計算其平均值就對了(別忘了要調整定義域的寬度)。