本科生《统计计算》教材，采用R语言和Julia语言，包括误差(2)

当\(N\)很大时\(\hat I\)近似服从 \(\text{N}(I, [M(b-a)]^2 p(1-p) / N)\)分布， 称此近似分布的方差\([M(b-a)]^2 p(1-p) / N\)为\(\hat I\)的渐近方差。 计算渐近方差可以用\(\hat p\)代替\(p\)估计为 \([M(b-a)]^2 \hat p(1 - \hat p) / N\)。 式说明 \(\hat I\)的误差为\(O_p(\frac{1}{\sqrt{N}})\), 这样，计算\(\hat I\)的精度每增加一位小数， 计算量需要增加100倍。随机模拟积分一般都服从这样的规律。

为了计算\(I = \int_a^b h(x) \,dx\)， 上面的随机投点法用了类似于舍选法的做法， 在非随机问题中引入随机性时用了二维均匀分布和两点分布， 靠求两点分布概率来估计积分\(I\)。 随机投点法容易理解， 但是效率较低， 另一种效率更高的方法是利用期望值的估计。 取\(U \sim \text{U}(a,b)\)， 则 \[\begin{aligned} E [h(U)] =& \int_a^b h(u) \frac{1}{b-a} du = \frac{I}{b-a} \\ I =& (b-a) \cdot E h(U)\end{aligned}\] 若取\(\{ U_i, i=1,\ldots,N \}\) 独立同U(\(a,b\))分布, 则\(Y_i = h(U_i), i=1,2,\dots,N\)是iid随机变量列, 由强大数律， \[\begin{aligned} \bar Y = \frac{1}{N} \sum_{i=1}^N h(U_i) \to E h(U) = \frac{I}{b-a}, \quad \text{a.s.} \ (N \to \infty)\end{aligned}\] 于是\[\begin{align} \tilde I = \frac{b-a}{N} \sum_{i=1}^N h(U_i) \tag{11.4}\end{align}\]

是\(I\)的强相合估计。称这样计算定积分\(I\)的方法为平均值法。

由中心极限定理有 \[\begin{aligned} \sqrt{N}(\tilde I - I) \stackrel{\mbox{d}}{\longrightarrow} \text{N}\big(0, (b-a)^2 \text{Var}(h(U)) \big)\end{aligned}\] 其中\[\begin{align} \text{Var}[h(U)] = \int_a^b \big[h(u) - Eh(U) \big]^2 \frac{1}{b-a} du \tag{11.5}\end{align}\]仅与\(h\)有关， 仍有\(\tilde I - I = O_p(\frac{1}{\sqrt{N}})\), 但是的渐近方差小于的渐近方差(见下一节)。 \(\text{Var}[h(U)]\)可以用模拟样本\(\{ Y_i = h(U_i) \}\)估计为\[\begin{align*} \text{Var}(h(U)) \approx \frac{1}{N} \sum_{i=1}^N (Y_i - \bar Y)^2.\end{align*}\]

如果定积分区间是无穷区间， 比如\(\int_0^\infty h(x) \, dx\)， 为了使用均匀分布随机数以及平均值法计算积分可以做积分变换， 使积分区间变成有限区间。 例如，作变换\(t = 1/(x+1)\), 则 \[\begin{aligned} \int_0^\infty h(x) \, dx = \int_0^1 h(\frac{1}{t} - 1) \frac{1}{t^2} \,dt.\end{aligned}\]