随机投点法计算定积分java_11 随机模拟积分 | 统计计算
11.4 高維定積分
上面的兩種計算一元函數定積分的方法可以很容易地推廣到多元函數定積分，
或稱高維定積分。
設\(d\)元函數\(h(x_1, x_2, \dots, x_d)\)定義于超矩形
\[\begin{aligned}
C = \{(x_1, x_2, \ldots, x_d): a_i \leq x_i \leq b_i, i=1,2,\ldots,d \}
\end{aligned}\]
且
\[\begin{aligned}
0 \leq h(x_1, \ldots, x_d) \leq M,
\ \forall x \in C.
\end{aligned}\]
令
\[\begin{aligned}
D =& \{(x_1, x_2, \ldots, x_d, y): (x_1, x_2, \ldots, x_d) \in C,
\ 0 \leq y \leq h(x_1, x_2, \ldots, x_d) \}, \\
G =& \{(x_1, x_2, \ldots, x_d, y): (x_1, x_2, \ldots, x_d) \in C,
\ 0 \leq y \leq M \}
\end{aligned}\]
為計算\(d\)維定積分
\[\begin{align}
I = \int_{a_d}^{b_d} \cdots \int_{a_2}^{b_2} \int_{a_1}^{b_1}
h(x_1, x_2, \ldots,x_d) \, dx_1 d x_2 \cdots dx_d,
\tag{11.18}
\end{align}\]
產生服從\(d+1\)維空間中的超矩形\(G\)內的均勻分布的獨立抽樣
\(\boldsymbol Z_1, \boldsymbol Z_2, \ldots, \boldsymbol Z_N\), 令
\[\begin{aligned}
\xi_i = \begin{cases}
1, & \boldsymbol Z_i \in D \\
0, & \boldsymbol Z_i \in G-D
\end{cases}, \quad i=1,2,\ldots,N
\end{aligned}\]
則\(\xi_i\) iid
b(1,\(p\)),
\[\begin{aligned}
p = P(\boldsymbol Z_i \in D) = \frac{V(D)}{V(G)}
= \frac{I}{M V(C)}
= \frac{I}{M \prod_{j=1}^d (b_j-a_j)}
\end{aligned}\]
其中\(V(\cdot)\)表示區域體積。
令\(\hat p\)為\(N\)個隨機點中落入\(D\)的百分比，則
\[\begin{aligned}
\hat p =& \frac{\sum \xi_i}{N} \to p,
\ \text{a.s.} (N \to \infty),
\end{aligned}\]
用
\[\begin{align}
\hat I_1 = \hat p V(G) = \hat p \cdot M \, V(C)
= \hat p \cdot M \prod_{j=1}^d (b_j-a_j)
\tag{11.19}
\end{align}\]
估計積分\(I\)，
則\(\hat I_1\)是\(I\)的無偏估計和強相合估計。
稱用式(11.19)計算高維定積分\(I\)的方法為隨機投點法。
由中心極限定理知
\[\begin{aligned}
\sqrt{N}(\hat p - p)/\sqrt{p(1-p)} \stackrelozvdkddzhkzd{\longrightarrow}&
\text{N}(0,1), \\
\sqrt{N}(\hat I_1 - I) \stackrelozvdkddzhkzd{\longrightarrow}&
\text{N}\left(0, \left(M \prod_{j=1}^n (b_j - a_j) \right)^2 p(1-p)\right),
\end{aligned}
\]
\(\hat I_1\)的漸近方差為
\[\begin{align}
\frac{\left(M \prod_{j=1}^n (b_j - a_j)\right)^2 p(1-p)}{N}
\tag{11.20}
\end{align}\]
所以\(\hat I_1\)的隨機誤差仍為\(O_p\left(\frac{1}{\sqrt{N}}\right)\),
\(N\to\infty\)時的誤差階不受維數\(d\)的影響，
這是隨機模擬方法與其它數值計算方法相比一個重大優勢。
在計算高維積分(11.18)時，
仍可以通過估計\(E h(\boldsymbol U)\)來獲得,
其中\(\boldsymbol U\)服從\(R^d\)中的超矩形\(C\)上的均勻分布。 設\(\boldsymbol U_i \sim\) iid
U(\(C\))，\(i=1,2,\ldots,N\), 則\(h(\boldsymbol U_i), i=1,2,\ldots,N\)是iid隨機變量列,
\[\begin{aligned}
E h(\boldsymbol U_i) = \int_C h(\boldsymbol u) \frac{1}{\prod_{j=1}^d (b_j - a_j)} d \boldsymbol u
= \frac{I}{\prod_{j=1}^d (b_j - a_j)},
\end{aligned}\]
估計\(I\)為
\[\begin{align}
\hat I_2 = \prod_{j=1}^d (b_j - a_j)
\cdot \frac{1}{N} \sum_{i=1}^N h(U_i),
\tag{11.21}
\end{align}\]
稱用式(11.21)計算高維定積分\(I\)的方法為平均值法。
由強大數律
\[\begin{aligned}
\hat I_2 \to \prod_{j=1}^d (b_j - a_j) E h(U) = I,
\ \text{a.s.} (N \to \infty),
\end{aligned}\]
\(\hat I_2\)的漸近方差為
\[\begin{align}
\frac{(V(C) )^2 \text{Var}(h(\boldsymbol U))}{N}
= \frac{\left(\prod_{j=1}^d (b_j - a_j) \right)^2 \text{Var}(h(\boldsymbol U))}{N}.
\tag{11.22}
\end{align}\]
\(N \to \infty\)時的誤差階也不受維數\(d\)的影響。
我們來比較隨機投點法(11.19)與平均值法(11.21)的精度，
只要比較其漸近方差。
對\(I = \int_C h(\boldsymbol x) d \boldsymbol x\)，
設\(\hat I_1\)為隨機投點法的估計，
\(\hat I_2\)為平均值法的估計。
因設\(0 \leq h(\boldsymbol x) \leq M\)，
不妨設\(0 \leq h(\boldsymbol x) \leq 1\), 取\(h(\boldsymbol x)\)的上界\(M=1\)。
令\(\boldsymbol X_i \sim\) iid U(\(C\)),
\(\eta_i = h(\boldsymbol X_i)\), \(Y_i \sim\) iid
U(0,1)與\(\{\boldsymbol X_i\}\)獨立，
\[\begin{aligned}
\xi_i = \begin{cases}
1, & \text{as} Y_i \leq h(\boldsymbol X_i) \\
0, & \text{as} Y_i > h(\boldsymbol X_i)
\end{cases}
\ i=1,2,\ldots,N,
\end{aligned}\]
這時有
\[\begin{aligned}
\hat I_1 &= V(C) \frac{1}{N} \sum_{i=1}^N \xi_i, \qquad
\hat I_2 = V(C) \frac{1}{N} \sum_{i=1}^N \eta_i \\
\text{Var}(\hat I_1) &= \frac{1}{N} V^2(C) \cdot \frac{I}{V(C)}
\left(1 - \frac{I}{V(C)} \right) \\
\text{Var}(\hat I_2) &= \frac{1}{N} V^2(C) \cdot \left(
\frac{1}{V(C)} \int_C h^2(\boldsymbol x)d \boldsymbol x
- \left( \frac{I}{V(C)} \right)^2 \right) \\
\text{Var}(\hat I_1) - \text{Var}(\hat I_2) &=
\frac{V(C)}{N} \int_C \left\{ h(\boldsymbol x) - h^2(\boldsymbol x) \right \} \,d\boldsymbol x \geq 0
\end{aligned}\]
可見平均值法精度更高。
事實上，隨機投點法多用了隨機數\(Y_i\)，必然會增加抽樣隨機誤差。
在計算高維積分(11.18)時， 如果用網格法作數值積分，
把超矩形\(C = [a_1, b_1] \times [a_2, b_2] \times \cdots \times [a_d, b_d]\)
的每個邊分成\(n\)個格子，就需要\(N=n^d\)個格子點，
如果用每個小超矩形的中心作為代表點， 可以達到\(O(n^{-2})\)精度，
即\(O(N^{-2/d})\)， 當維數增加時為了提高一倍精度需要\(2^{d/2}\)倍的代表點。
比如\(d=8\)，精度只有\(O(N^{-1/4})\)。
高維的問題當維數增加時計算量會迅猛增加， 以至于無法計算，
這個問題稱為維數詛咒。
如果用Monte Carlo積分，則精度為\(O_p(N^{-1/2})\)， 與\(d\)關系不大。
所以Monte Carlo方法在高維積分中有重要應用。
為了提高積分計算精度，需要減小\(O_\text{p}(N^{-1/2})\)中的常數項，
即減小\(\hat I\)的漸近方差。
例11.3考慮
\[\begin{aligned}
& f(x_1, x_2, \ldots, x_d) = x_1^2 x_2^2 \cdots x_d^2,
\ 0 \leq x_1 \leq 1, 0 \leq x_2 \leq 1, \dots, 0 \leq x_d \leq 1
\end{aligned}\]
的積分
\[\begin{aligned}
I = \int_0^1 \cdots \int_0^1 \int_0^1
x_1^2 x_2^2 \cdots x_d^2
\, dx_1 dx_2 \cdots dx_d
\end{aligned}\]
當然，這個積分是有精確解\(I = (1/3)^d\)的，
對估計\(I\)的結果我們可以直接計算誤差。
以\(d=8\)為例比較網格點法、隨機投點法、平均值法的精度，
這時真值\(I = (1/3)^8 \approx 1.5241587 \times 10^{-4}\)。
用\(n\)網格點法，\(N=n^d\)，公式為
\[\begin{aligned}
I_0 = \frac{1}{N} \sum_{i_1=1}^n \sum_{i_2=1}^n \cdots \sum_{i_d=1}^n
f(\frac{2i_1-1}{2n}, \frac{2i_2-1}{2n}, \ldots, \frac{2i_d-1}{2n})
\end{aligned}\]
誤差絕對值為\(e_0 = | I_0 - I |\)。
如果取\(n=5, d=8\)，需要計算\(N=390625\)個點的函數值，
計算量相當大。
計算程序：
積分真值：
網格法:
網格法誤差和相對誤差:
相對誤差約8%，誤差較大。
用隨機投點法求\(I\)，
先在\((0,1)^d \times (0,1)\)均勻抽樣
\((\xi_i^{(1)}, \xi_i^{(2)}, \ldots, \xi_i^{(d)}, y_i), i=1,\ldots, N\)，
令\(\hat I_1\)為\(y_i \leq f(\xi_i^{(1)}, \xi_i^{(2)}, \ldots, \xi_i^{(d)})\)成立的百分比。
一次估計的程序和結果為
估計RMSE、MAE和MRE來度量誤差大小：
估計的相對誤差約為10個百分點，
誤差較大。
為了確認這樣的誤差估計，
將隨機投點法重復模擬\(B=100\)次，得到100個\(I\)的估計，程序為：
用\(B=100\)次重復模擬估計平均相對誤差：
這確認了從一次模擬給出的平均相對誤差估計，兩者基本相同。
用平均值方法，公式為
\[\begin{aligned}
\hat I_2 = \frac{1}{N} \sum_{i=1}^N f(\xi^{(1)}_i, \xi^{(2)}_i, \cdots, \xi^{(d)}_i)
\end{aligned}\]
其中\(\xi^{(j)}_i, i=1, \ldots, N, j=1,\ldots,d\) 獨立同U(0,1)分布。
一次計算的程序為:
平均值法的各種誤差度量為：
相對誤差約為1.3個百分點。
類似地， 重復\(B=100\)次，估計MRE：
用\(B=100\)次重復模擬估計平均相對誤差：
相對誤差約為1.4個百分點，
與從一次模擬估算的平均相對誤差相同。
注意實際編程時，
隨機投點法和平均值法可以使用同一組隨機數，
可以大大減少運算量。
隨機模擬程序是消耗計算資源特別嚴重的計算任務之一，
設計更高效的算法、程序進行優化、考慮并行計算，
都是隨機模擬問題的實現需要考慮的。
對隨機模擬問題，
高效的算法主要指估計的方差更小，
這樣誤差小，
達到一定精度需要的計算量就小。
最后，三種方法計算量相同(都計算\(N=5^8\)次函數值)的情況下，
平均相對誤差分別為：
隨機投點法的誤差與網格法相近；
平均值法的誤差只有隨機投點法的13%。
※※※※※
總結
以上是生活随笔為你收集整理的随机投点法计算定积分java_11 随机模拟积分 | 统计计算的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 卖房包过户还是不包过户好
- 下一篇： 大学数学建模大赛是用计算机,北京大学第十