复杂积分的计算机模拟数值计算方法

本篇文章讲讲Monte Carlo方法在数值计算上的应用，这一问题是MCMC（Markov Chain Monte Carlo）采样的子问题。
假定需求积分为$I=\int_{\theta }P(\theta )d\theta$，其中$P(\theta )$的积分式很难求得，比如$P(\theta )=\frac{sin\theta }{\theta }$。

其问题可以转化为求一个随机变量的函数的期望的问题。设$\theta\sim Q(\theta)$，$\theta$服从概率密度函数为$Q(\theta)$的分布。将$\frac{P(\theta)}{Q(\theta)}$看作是关于随机变量$\theta$的函数，上述积分可转化为$I=\int_{\theta }\frac{P(\theta)}{Q(\theta)}\cdot Q(\theta)d\theta=E[\frac{P(\theta)}{Q(\theta)}]$关于随机变量$\theta$的函数$\frac{P(\theta)}{Q(\theta)}$的数学期望。

Monte Carlo方法的核心是对于随机变量$\theta$大规模采样，可以依$\widehat{E}[\frac{P(\theta)}{Q(\theta)}]=\frac{1}{n}\sum_{i=1}^{n}\frac{P(\theta_{i})}{Q(\theta_{i})}$得到期望的估计值$\widehat{E}$.

依大数定律：

$$lim_{n \to \infty }\widehat{E}=E$$
当采样的规模足够大时,有：
$$\widehat{E}\approx E$$

需要注意的是若积分限为[a,b]则所选取的$\theta\sim Q(\theta)$中的概率密度函数的定义域也是[a,b]。

举个例子，还是以$P(\theta )=\frac{sin\theta }{\theta }$为例，设所求积分为$I_{1}=\int_{0 }^{1}\frac{sin\theta }{\theta} d\theta$，设$\theta\sim U(0，1)$服从[0,1]上的均匀分布，概率密度函数为$f(\theta)=1, (0\leqslant \theta \leqslant 1)$。

$$I_{1}==\int_{\theta }\frac{sin\theta }{\theta}\cdot 1 d\theta=E[\frac{sin\theta }{\theta }]$$
积分$I_{1}$的近似估计值为:
$$I_{1}\approx \frac{1}{n}\sum_{i=1}^{n}\frac{sin\theta_{i}}{\theta_{i}}$$

Python代码：

#coding:utf-8
import numpy as np
import random

n = 10000 #采样数

def func(x):
    x = random.uniform(0, 1)
    return np.sin(x)/x

print sum(map(func, range(n)))/n

#输出为0.945611582984

计算$sin\theta$在[0,1]内积分，符号结果为$cos0 - cos1=1 - cos1$

#将上述函数改为np.sin(x) 后输出为0.460831340906
>>> import numpy as np
>>> 1 - np.cos(1)
0.45969769413186035

可见采样数为10000时结果相差不大。