[学习] 数值积分与蒙特卡洛方法
前言
这学期选了门期待已久的课，叫《计算机图形学研究进展》（神课，我在别的文章里提到过它），其中有一项作业是 paper reading：选择一篇近几年的 SIGGRAPH 论文，撰写阅读报告。
找来找去，最后找了一篇讲蒙特卡洛方法优化的论文，其它的论文要么满篇都是高深的数学推导，要么是神乎其神玄乎其玄的神经网络，把参数输进去大锅乱炖，最后产生一些莫名其妙的输出……（很遗憾，因为我对 AI 方面的东西基本一窍不通，所以得出的就是这样的结论）
话虽如此，我对于蒙特卡洛方法的认知仅仅局限于听说过这个名词，好像它可以用来估计圆周率。因此，我决定快速学习一下相关的知识。
如果你好奇为什么要 “快速” 学习，因为这门课的大作业更为硬核——复现一篇近几年的 SIGGRAPH 论文，如果想做出点有意思的东西，就得尽早把其它事做完。
我初步搜索了下，了解到蒙特卡洛方法一般用来计算数值积分。由于自己并非计算机专业，不要说蒙特卡洛方法了，连数值积分也没有接触过，所以打算先铺垫一点数值积分的基础，再学习蒙特卡洛方法。
以下内容是我的学习笔记。
数值积分
数值积分是利用黎曼积分等数学定义，用数值逼近的方法近似计算给定的定积分值。——《百度百科：数值积分》
为什么需要数值积分？
这里提供我的两种理解。
- 计算机是离散的
给你一个函数 \(f(x)=e^x\)，如果让你计算 \(\int{}_1^2f(x)dx\)，根据牛顿—莱布尼茨公式，你首先会写出其原函数 \(F(x)=e^x\)，然后就能得到答案为 \(F(2)-F(1)=e^2-e\)。问题在于，写出其原函数 对于人类来说可能非常轻松，但对于计算机则不然。准确地来说，对于给计算机编程的程序员来说，实现“计算任意给定函数的原函数” 并不容易用代码实现，因为除去最简单的多项式函数，许多函数可能包含指数、对数，甚至由一些非初等函数构成。
我们需要某种适合计算机的算法！
- 许多积分不具有解析表示
许多积分的原函数是非初等函数，举个经典例子：\(\int{}e^{x^2}dx\)，由于无法解析地表示 \(e^{x^2}\) 的原函数，因此不能使用牛顿—莱布尼茨公式。
梯形（trapezoid）法
人们已经发明出了许多数值积分方法，由于时间问题，我只学习了其中最简单的一种——梯形法。
回想积分的原始定义（黎曼和）：
如上。我们将积分区间 $[a,b] $划分为许多小区间 \([x_0,x_1], [x_1, x_2],\ldots{},[x_{n-1},x_n]\)，计算这些区间与函数值围成的矩形面积之和，当区间个数取得足够多，每个区间足够短时，上式就等于 \(\int_{a}^{b}f(x)dx\)。
因此，我们可以从定义出发，计算一个函数 \(f(x)\) 在 \([a,b]\) 上的定积分：
既然我们已经想到用矩形来估算定积分，为什么不试试梯形呢？直观上看，梯形更能贴合函数变化的趋势：
（图源：FNC 5.6. Numerical integration）
写成数学表达式就是：
为了证明它确实有效，我写了一个简单的程序估算 \(\int_{1}^{2}e^xdx\) 的值：
from math import *
def trapezoid(f, a, b, n):
res = 0
dt = (b - a) / n
offset = a
for i in range(n):
res += (f(offset) + f(offset + dt)) / 2 * dt
offset += dt
return res
print("True Value: {}".format(exp(2) - exp(1)))
for i in range(1, 9):
segNum = 10**i
print(
"Estimated Value ({} segments): {}".format(segNum, trapzoid(exp, 1, 2, segNum))
)
运行结果：
True Value: 4.670774270471606
Estimated Value (10 segments): 4.674665933799437
Estimated Value (100 segments): 4.670813193525656
Estimated Value (1000 segments): 4.670774659702487
Estimated Value (10000 segments): 4.670774274363632
Estimated Value (100000 segments): 4.670774270528357
Estimated Value (1000000 segments): 4.670774270248262
Estimated Value (10000000 segments): 4.6707742720582806
Estimated Value (100000000 segments): 4.670774253951506
注意区间并不一定需要等间隔分割，有时候根据待积函数的特性进行变间隔分割，或许可以得到更精确的估计结果，这里只是为了方便，采取等间隔分割的策略。
蒙特卡洛方法
蒙特卡洛方法是计算数值积分的强大工具，相比于梯形法，它的思想和前者截然不同——它是一种基于概率的方法。
我觉得概率真的是非常非常神奇的东西：第一，概率符合我们的直觉，因为生活中绝大多数的事件都是按概率发生的事件；第二，许多复杂的组合优化问题都可以使用概率的思想大大地简化；第三，概率和信息论直接相关，而后者有许多有趣的应用，等等。
有些扯远了。不过我想让你知道的是，我一直都很喜欢概率和概率论，无论是从知识的角度，还是在游戏中，抑或是对于单纯的概率练习题来说……还有，我的概率论是满绩，虽然现在已经忘记了大半。
PS 1：下面的内容会涉及一些概率论的术语，假定你已经了解过它们——本文只会进行简单的回顾。
PS 2：仅针对一维展开讨论，结论可以很容易地推广到高维
首先，我们需要一些前置知识：
概率密度函数与（累积）分布函数
对于连续型随机变量 \(X\sim{}p(x)\)，其中 \(p(x)\) 称为随机变量 \(X\) 的 概率密度函数（Probability Density Function, PDF），满足 \(\int_{-\infty}^{\infty}p(x)dx=1\) 且 \(p(x)\geq{}0\)。注意：\(p(x)\) 本身的值并不一定小于 \(1\)。
记 \(F(x)\) 满足 \(\frac{dF(x)}{dx}=p(x)\)，\(F(x)\) 称为 累积分布函数（Cumulative Distribution Function, CDF）。
根据累积分布函数的定义，有：
数学期望
定义连续型随机变量 \(X\) 的数学期望为 \(E(X)=\int_{-\infty}^{\infty}xp(x)dx\)，相应地可以得到 \(E[f(X)]=\int_{-\infty}^{\infty}f(x)p(x)dx\)，其中 \(f\) 为任意函数。
（辛钦）大数定律
设 \(X_1,X_2,\ldots{},X_n\) 是相互独立的随机变量序列，服从相同分布，且具有有限的数学期望 \(\mu\)，则 \(\forall\varepsilon>0\)，有
简单地讲就是频率趋近于概率，比如投掷硬币，虽然刚开始可能出现连续的正/反面，但当投掷足够多次之后，正/反面出现的频率趋近于概率 \(\frac{1}{2}\)
蒙特卡洛方法的定义
我们的目标是计算 \(\int_{a}^{b}h(x)dx,\ a<b\)，其中 \(h(x)\) 的表达式已知，为了简化问题，我们还要求 \(h(x)\) 在区间 \([a,b]\) 处处可积。
根据大数定律，如果我们能构造一个相互独立且同分布的随机变量序列 \(X_1,X_2,\ldots{},X_n\) 使得 \(E(X)=\int_{a}^{b}h(x)dx\)，那么只需要 \(n\) 足够大，\(\frac{1}{n}\sum_{i=1}^n X_i\) 就会依概率收敛到 \(\int_{a}^{b}h(x)dx\)。通过控制 \(n\) 的大小，我们可以调整估算的精度，或者说，控制误差的大小。
根据期望的定义：
其中 \(p(x)\) 为随机变量 \(X\) 的概率密度函数。
对比 \(E(X)=\int_{a}^{b}h(x)dx\)，如果令 \(f(x)=\frac{h(x)}{p(x)}\)，那么 \(E[\frac{h(X)}{p(X)}]=\int_{a}^{b}\frac{h(x)}{p(x)}p(x)dx=\int_{a}^{b}h(x)dx\)，不就是我们的目标吗？
这就是蒙特卡洛方法的基本思想。
定义随机变量 \(Y=\frac{h(X)}{p(X)}\)，其中 \(X\) 的概率密度函数为 \(p(x)\)，则 \(E(Y)=\int_{a}^{b}h(x)dx\)。我们可以使用 \(F_n=\frac{1}{n}\sum_{i=1}^n Y_i=\frac{1}{n}\sum_{i=1}^n\frac{h(X_i)}{p(X_i)}\) 来估计 \(\int_{a}^{b}h(x)dx\)，\(n\) 越大则结果越准确。
问题：如何选择 \(p(x)\) ？
接下来的问题是如何选取合适的 \(p(x)\)。注意，在上面的推导中，我们没有对 \(p(x)\) 给出任何限制，这意味着任何满足 \(\int_{a}^{b}p(x)=1\) 且 \(p(x)\geq{}0\) 的 \(p(x)\) 都是可以的。
当然，我们并不能随意选择 \(p(x)\)，否则本节也就没有存在的必要了。
主要有两个因素（应该不止两个，但我目前只学到这两个）制约着 \(p(x)\) 的选择。
因素 1：收敛率（convergence rate）
仅仅知道 \(F_n\) 能依概率收敛到 \(E(Y)\) 还不够，还需要知道收敛得有多快，因为计算资源非常珍贵，若能使用很小的 \(n\) 就达成很高的精度，那么就不需要增大 \(n\)。
方差（variance）可以用来衡量 \(F_n\) 的收敛程度：
因为方差的定义里面有平方，所以两边开方，得到误差 \(\varepsilon{}\propto\frac{1}{\sqrt{n}}\)，也就是收敛的程度与 \(\frac{1}{\sqrt{n}}\) 成正比。
由于 \(p(x)\) 的选择不同，\(V(Y)\) 也会不同。尽管从渐进的意义上看，不同的 \(p(x)\) 对应的收敛率都相同，但更小的 \(V(Y)\) 意味着能更快地收敛，所以选择使方差尽可能小的 \(p(x)\) 是有必要的。
重要性采样（Importance Sampling）
这里举一个通过选取合适的 \(p(x)\) 来降低方差的例子。
很容易想到，最简单的 \(p(x)\) 是均匀分布，即 \(p(x)=\frac{1}{b-a}\)，这里我写了一个程序，和前面一样估算 \(\int_{1}^{2}e^xdx\) 的值，不过用的是蒙特卡洛方法：
from math import *
import random
def mc(f, a, b, n):
res = 0
for i in range(n):
x = random.uniform(a, b)
res += f(x) / (1 / (b - a))
res /= n
return res
print("True Value: {}".format(exp(2) - exp(1)))
for i in range(1, 9):
segNum = 10**i
val = mc(exp, 1, 2, segNum)
print("Estimated Value ({} samples): {}".format(segNum, val))
运行结果：
True Value: 4.670774270471606
Estimated Value (10 samples): 4.773362441512867
Estimated Value (100 samples): 4.748707630229044
Estimated Value (1000 samples): 4.640803580143513
Estimated Value (10000 samples): 4.642255632409945
Estimated Value (100000 samples): 4.66787982737648
Estimated Value (1000000 samples): 4.6719327404508855
Estimated Value (10000000 samples): 4.671300362268291
Estimated Value (100000000 samples): 4.670712953211018
结果并不是非常令人满意。对比前面的梯形法，你会发现使用蒙特卡洛方法，在 \(n\) 相同的情况下误差要更大。
那么，如何在 \(n\) 相同的前提下提高估算的精度 / 减小误差 / 降低方差呢？
信息论告诉我们：不确定性和信息相挂钩。如果我们将蒙特卡洛方法视为对 待积函数 \(h(x)\) 的 \(n\) 次采样，那么提高精度的关键就在于增加每次采样所能获得的信息量。
考虑这样一个函数 \(h(x)\)：
（图源：pbr book）
如果我们想求 \(\int_0^1h(x)dx\)，此时采用 \([0,1]\) 上的均匀分布作为 \(p(x)\) 就不太恰当了，因为函数只在 \(0.5\) 附近有非零值，其它部分对积分没有贡献。如果我们能针对 \(0.5\) 附近的区间进行采样，可以想象，估算的精度会得到提高。
此方法称为重要性采样。在重要性采样中，我们一般选取和 \(h(x)\) 形状相似的 \(p(x)\)，这样就能保证每次采样平均获得的信息量比均匀分布采样要大。
依然以估算 \(\int_{1}^{2}e^xdx\) 的值为例，因为 \(h(x)\) 是一个上凹函数，所以我选择 \(p(x)=cx^2,\ 1\leq{}x\leq{}2\)（这里只是我随便想的，最终结果不一定会更好），\(c\) 是归一化系数，令 \(\int_{1}^{2}cx^2dx=1\) 可得 \(c=\frac{3}{7}\) 。
那么我们就可以写出程序……问题来了，这样的非均匀分布要怎么实现呢？
请继续看。
因素 2：用程序产生符合 \(p(x)\) 的分布并不容易
绝大部分的编程语言所提供的随机数发生器服从的都是 均匀分布，或者说它们的概率密度函数在给定区间上的取值是一个常数。（其实还有正态分布、泊松分布之类的，不过这些分布事实上都可以用均匀分布生成）
如果 \(p(x)\) 是非均匀分布呢？我们需要设计某种算法，将均匀分布映射到想要的非均匀分布。这里提供两种方法。
- 方法 1：逆累积分布函数
先考虑离散的情况。
定义随机变量 \(X\) 的分布列如下：
| X | P |
|---|---|
| 123 | 0.2 |
| 456 | 0.5 |
| 789 | 0.3 |
现在让你用程序生成符合该分布的随机数。
假设标准库提供了 rand() 函数用于生成 \(0\sim{}1\) 的随机数，不妨将这个随机数发生器抽象为随机变量 \(\xi\)，其概率密度函数 \(p(x)=1\)。
很容易可以想到下面这种方法：
\(\xi{}\) 落在哪个范围内，就输出对应的随机数（123、456 或 789）。
或者，可以写成累积分布函数的形式：
程序实现略。
现在考虑连续的情况。所谓连续，就是离散的极限情况，当这些分段的数量趋近于正无穷，每个分段对应的就不再是概率，而是概率密度，此时的累积分布函数可能形如下图：
依样画葫芦，如果我们想知道此时 \(\xi{}\) 对应的 \(X\) 取值为多少，就需要把累积分布函数 \(F(x)\) 的反函数 \(F^{-1}(x)\) 解出来，代入 \(F^{-1}(\xi{})\) 就得到了此时随机变量 \(X\) 的取值。
现在我们来实现前面没有实现的重要性采样程序，其中 \(p(x)=\frac{3}{7}x^2,\ 1\leq{}x\leq{}2\) 。
首先计算得到 \(F(x)=\int_{1}^xp(x)dx=\frac{1}{7}x^3-\frac{1}{7}\)，所以 \(F^{-1}(x)=\sqrt[3]{7x+1}\)
写出程序（重复的代码就不贴上来了）：
def mc(f, a, b, n):
res = 0
for i in range(n):
x = (7 * random.uniform(0, 1) + 1) ** (1 / 3)
res += f(x) / (3 / 7 * x**2)
res /= n
return res
运行结果：
True Value: 4.670774270471606
Estimated Value (10 samples): 4.741825153692039
Estimated Value (100 samples): 4.695969646551462
Estimated Value (1000 samples): 4.68496747191308
Estimated Value (10000 samples): 4.670736267923529
Estimated Value (100000 samples): 4.669402728657563
Estimated Value (1000000 samples): 4.671530413109727
Estimated Value (10000000 samples): 4.671033902580099
Estimated Value (100000000 samples): 4.670737814323587
变化不大，不过收敛值是正确的。
其实这里的例子用得不太好，刚开始写本文的时候还没考虑这么多。如果待积函数和前面那个仅在 \(0.5\) 附近有非零值的函数较为相似的话，那么重要性采样应该就能发挥很大作用了。请读者自行尝试。
- 方法 2：接受—拒绝采样
逆累积分布函数法需要知道累积分布函数的反函数，但有时候这个反函数非常难表示，或者能表示但计算代价非常大，此时就需要使用别的方法，其中一种方法称为接受—拒绝采样。
之前学编程的时候，你可能写过如下的代码：
while(不满足条件)
{
重新生成
}
程序会一直执行 while 循环里面的语句，直至条件被满足，这就是接受—拒绝采样的基本原理。
首先我们选取满足如下条件的函数 \(g(x)\) 和常数 \(C\)，满足 \(C\cdot{}g(x)\geq{}p(x)\) 在积分区间内恒成立:
（图源：如何理解蒙特卡洛方法的接受拒绝采样？）
接下来按照如下的伪代码进行采样：
do
{
生成一个随机数 x（在积分区间 [a, b] 上均匀分布）
计算 y = p(x) / (C * g(x))
生成一个随机数 u（在 [0, 1] 上均匀分布）
} while(u >= y);
循环会一直执行，直到 u < y，此时称为 “接受”，并将 \(x\) 作为结果；否则称为 “拒绝”。
如果你理解了逆累积分布函数的思想，应该也能很快理解接受—拒绝采样。这里并不论证该方法的正确性和可靠性，仅给出代码实现。
被积函数和前面相同，不再复述。取 \(Cg(x)=\frac{9}{7}x-\frac{6}{7}\) 在 \([1,2]\) 上恒大于等于 \(p(x)=\frac{3}{7}x^2\)，这也是我随便取的。
def mc(f, a, b, n):
res = 0
for i in range(n):
x = random.uniform(1, 2)
y = 3 / 7 * x**2 / (9 / 7 * x - 6 / 7)
u = random.uniform(0, 1)
while u >= y:
x = random.uniform(1, 2)
y = 3 / 7 * x**2 / (9 / 7 * x - 6 / 7)
u = random.uniform(0, 1)
res += f(x) / (3 / 7 * x**2)
res /= n
return res
运行结果：
True Value: 4.670774270471606
Estimated Value (10 segments): 4.586496885027307
Estimated Value (100 segments): 4.8672148582933765
Estimated Value (1000 segments): 4.8679642905803995
Estimated Value (10000 segments): 4.855200375221184
Estimated Value (100000 segments): 4.856935634120809
Estimated Value (1000000 segments): 4.856994731849525
Estimated Value (10000000 segments): 4.857045843023234
Estimated Value (100000000 segments): 4.856954213978123
这里的误差比前面还要大，主要有两个原因：
- 重要性采样选取的 \(p(x)\) 不合适（和待积函数 \(h(x)\) 的相似度不高）
- 建议分布函数 \(Cg(x)\) 不合适
一些没有讲到的东西（将来可能会补充）
-
梯形法的误差分析
-
无偏 vs. 有偏（到目前为止我们讨论的都是无偏估计。有些 \(F_n\) 的期望并不等于 \(\int_a^bh(x)dx\)，但它们有更高的收敛率）
-
其它提高收敛率的方法：参考 pbr-book
总结
本文简单地介绍了数值积分的概念，分析了蒙特卡洛方法计算定积分的原理、可行性、收敛率和提高收敛率的措施，提供了两种生成符合任意分布的随机数的方法。
参考资料
-
[1] Matt Pharr, Wenzel Jakob, Greg Humphreys. Physically Based Rendering: From Theory To Implementation[EB/OL]. https://www.pbr-book.org
-
[2] 李航. 统计学习方法[M]. 第 2 版. 北京：清华大学出版社，2019
（因为学得比较匆忙，所以可能会出现许多不恰当的表述和错误的地方。欢迎你在评论区提出意见，在此提前感谢！）
