蒙特卡洛估算定积分平均值法Python实现
由think创建，最终由think 被浏览 7 用户
如何估算任意函数积分
对于给定的函数，我们可以通过公司计算积分。但对于计算机，对于任何复杂的函数，我们如何运行算力优势，快速估算积分。附python代码实现完整notebook，可以直接克隆运行。
原理
直观理解如上图所示，采样后求平均。写成求和公式：
$ \overline E = (b - a) \frac{1}{N} \sum\limits_{i=1}^{N}f(x_{i})$
原理来自于“强大数定理”，可证明得到以下公式：
\begin{equation*} Pr(\lim\limits_{N\rightarrow\infty}\frac{1}{N}\sum_{i=1}^{N}g(X_{i})=\int_{a}^{b}f(x)dx )=1 \end{equation*}
其中，$g(X_{i})=\frac{f(X_{i})}{p(X_{i})}$，$p(x)$是概率密度函数pdf (probability density function)。 也就是说在$p(x)$的分布下，抽取$x$样本，当N足够大时，可以采用均值来近似$f(x)$的积分：
\begin{equation*} \int _{a}^{b}{f(x)dx} \approx \frac{g(x_1) + g(x_2) + ... + g(x_N)}{N} \end{equation*} 以上图片的例子是选取了均匀分布函数$p(x)=\frac{1}{b-a}$的情况。
一般来说均匀分布计算的效率不会高，需要另外选取一个合理的分布提高计算效率，例如高斯函数在峰值附近积分贡献比较多，应该赋予更多的采样。这个技术叫做“重要性采样” (importance sampling)。均匀分布、高斯分布、Gamma分布等都属于简单分布采样（或者概率分布采样）。
如果是采用比较复杂的分布$p(X_{i})$，即使已知该分布，很难直接获取符合该分布的样本${x_{1}, x_{2}, ..., x_{N}}$，这时候需要使用“接受-拒绝采样”、 “马尔科夫链-蒙特卡洛 (MCMC) 采样”、“ Metropolis-Hastings (M-H)采样”、“Gibbs采样”等方法。
代码实现
定义要求积分的函数代码
# 函数 f(x) = x^2 求积分
def f(x):
return x ** 2
模特卡洛求积分代码
# 模特卡洛求解积分
def monte_carlo(func, a, b, sample_count):
"""
args:
func: 待求积分的函数
a: 下界
b: 上界
sample_count: 抽样数量，越大结果越准确，耗时越多
"""
# 从 [a, b) 间随机抽样 sample_count 的个数，e.g. x = [0.123, 0.234, ..]
x = np.random.uniform(a, b, sample_count)
# 计算在这些点上的 f(x) 的值, e.g. y = [0.123**2, ..]
y = func(x)
# 计算y的平均值 * (b - a)作为积分的值
return y.mean() * (b - a)
测试代码
测试从10轮到1000000轮不同采样次数的情况，可以看到随着轮次的增加，误差在快速下降
def test_monte_carlo(func, a, b):
"""测试一下"""
# 用 scipy 的 integrate 计算一下，作为实际值
from scipy.integrate import quad
expected = quad(func, a, b)[0]
for sample_count in [10, 100, 1000, 10000, 100000, 1000000]:
start_time = time.time()
actual = monte_carlo(func, a, b, sample_count)
T.log.info(f"--- 模特卡洛采样 ({sample_count}次)")
T.log.info(f"实际值: {expected}")
T.log.info(f"蒙特卡洛估计值: {actual}")
T.log.info(f"误差(比例): {abs(expected - actual)}({abs((expected - actual) / expected) * 100}%)")
T.log.info(f"耗时: {time.time() - start_time}s")
另一个方法：投点法
随机生成N个点，在$f(x)$曲线个数比例，乘以方形面积，就是函数f(x)的积分值。这个方法比较直观。参考代码如下
def MC_1(): # 蒙特卡洛求定积分1：投点法
n = 10**7
x_min, x_max = 0.0, 1.0
y_min, y_max = 0.0, 1.0
count = 0
for i in range(0, n):
x = random.uniform(x_min, x_max)
y = random.uniform(y_min, y_max)
# x*x > y，表示该点位于曲线的下面。所求的积分值即为曲线下方的面积与正方形面积的比。
if x * x > y:
count += 1
integral_value = count / n
return integral_value
代码Fork
点击克隆代码，可以在线运行。
https://bigquant.com/experimentshare/5b01b23e50574015b99dbbf6ea545964
\