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)} \stackrel{d}{\longrightarrow}&
\text{N}(0,1), \\
\sqrt{N}(\hat I_1 - I) \stackrel{d}{\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%。
※※※※※
