我們可以使用蒙特卡羅模擬來近似積分的值。\int_a^b f(x) dx 可以解釋成曲邊梯形的面積（ 假設f(x)>=0 for all x \in [a,b] ), 蒙特卡羅方法適用于使用使用解析方法無法找到的的近似情形，尤其適用于高維積分的近似。

我們以一個簡單的例子介紹一下蒙特卡羅積分方法。
我們編寫蒙特卡羅模擬的北太天元腳本，以近似由y=x^2，y=0，x=0，x=1限定的曲邊三角形S的面積。這實際上是 \int_0^1 x^2 dx .

當然，對于這種情況，我們可以使用解析方法找到確切的值:
??????? \int_0^1 x^2 dx=1/3。
使用蒙特卡羅積分，我們可以使用隨機數(shù)來近似這個積分值。一種方法是從已知面積的矩形中重復選取隨機點，并確定哪些點屬于曲邊三角形S。請注意，這個曲邊三角形S完全位于單位正方形D=[0,1]×[0,1]中。我們可以通過取兩個均勻的隨機數(shù)x和y來選擇D中的隨機點。然后我們跟蹤y<=x^2的次數(shù)，以及位于曲邊三角形S中的隨機點數(shù)。

%北太天元用Monte Carlo(蒙特卡洛)方法計算 int_0^1 x^2 dx
n = 1000;
x = rand(n,1);
y = rand(n,1);?? ?
和 = sum(? double(x.^2 > y) );
sprintf("x^2 從0積到1的積分值是 %3.4f\n", 和/n)
ind1 = find(? y <= x.^2 );
sh1 = scatter( x(ind1), y(ind1), 'filled');
set(sh1, 'SizeData', 500);
hold on ?
ind2 = find( y > x.^2 );
sh2 = scatter( x(ind2), y(ind2), 'filled');
set(sh2, 'SizeData', 500);
legend('y=x^2曲線下方的點', 'y=x^2曲線上方的點','FontSize',20)
title("蒙特卡羅方法計算 x^2 的積分")
hold off

用1000個樣本點進行蒙特卡羅積分，以近似S的面積。藍點為
曲邊三角形S中的隨機點，紅點是不在曲邊三角形S中的任意點。
紅點的比例就是曲邊三角形的面積:
上面北太天元代碼你運行后可以能看到輸出:
"x^2 從0積到1的積分值是 0.3580\n"

上面使用的是n=400, 如果把 n 增加到2000，我們可以看到

"x^2 從0積到1的積分值是 0.3330\n"

這距離精確值1/3 更加接近了。

此時藍紅點的分布如下圖

還有另一種方法來構建蒙特卡羅積分技術。假設我們想要近似曲邊三角形R的面積

\int_a^b f(x) dx.

我們不是在該區(qū)域的矩形中選取成對的點而是在區(qū)間[a，b]中的隨機選取x_i, 一次選取就稱為一個實驗。如果我們重復這個實驗，f(x_i)的平均值應該接近 1/(b-a)\int_a^b f(x) dx 的平均值。在適當?shù)臈l件下，可以證明當n逼近∞時，f(x_i)*(b-a)逼近積分值 \int_a^b f(x) dx.

在北太天元學習系列中的，我們介紹了可以使用北太天元的內置函數(shù)rand
得到服從[0，1]上的均勻分布的隨機數(shù)，然后做一個伸縮平移變換后得到
y=（b?a）x+a 是一個服從[a,b]區(qū)間上的均勻分布的隨機數(shù)。因此，上面我們介紹的
[0,1]向上的積分方法可以用這個方式遷移到[a,b]上的積分。

與任何數(shù)值積分方法一樣, 我們不能期望一定完全精確地得到積分值。使用蒙特卡羅方法，我們可以觀察到誤差隨著隨機樣本數(shù)量的增加而減小。在蒙特卡羅積分法中，
通過計算平均值A_n= ( f(x_1)+...+f(x_n) )/n 作為近似積分值。有定理表明
A_n的方差接近 σ^2/n。 雖然我們可能不知道σ，但我們知道增加n會導致A_n的方差
以1/n成比例的速率減小。從而，標準偏差以與1/sqrt(n) 成比例的速率遞減.

北太天元可以用下面的命令來計算 \int_0^{1} x^2 dx 的數(shù)值積分
大家可以觀察到誤差大致還是隨著n增加而減小。
for k=1:8
?n = 100*2^k;
?disp(sprintf("n = %d , 誤差 %f \n",n ,? abs( mean( ( rand(n,1) ).^2? ) -1/3) ));
end

? 1x1 string
??? "n = 200 , 誤差 0.002196 \n"
? 1x1 string
??? "n = 400 , 誤差 0.001769 \n"
? 1x1 string
??? "n = 800 , 誤差 0.006589 \n"
? 1x1 string
??? "n = 1600 , 誤差 0.001157 \n"
? 1x1 string
??? "n = 3200 , 誤差 0.001773 \n"
? 1x1 string
??? "n = 6400 , 誤差 0.004136 \n"
? 1x1 string
??? "n = 12800 , 誤差 0.002946 \n"
? 1x1 string
??? "n = 25600 , 誤差 0.001589 \n"