Monte Carlo方法由冯·诺伊曼于二战时提出,1777年法国人布丰用此思路估计pi值被认为是Monte Carlo的起源,这个方法简单又快速,通过产生随机数,将数值计算问题变成随机变量的某些特征值(比如期望值)。

积分运算,和估计pi值一样,用hit and miss法估计。

hit_miss <- function(fun, lower, upper, n=500) {
    # Monte Carlo integration using the hit and miss method
    x <- runif(n, lower, upper)
    f.value <- sapply(seq(lower, upper, 0.001), fun)
    f.min <- min(f.value)
    f.max <- max(f.value)
    y <- runif(n, f.min, f.max)
    hits <- sum(y <= fun(x))
    area <- (upper-lower) * f.min + hits/n * (upper-lower) * (f.max-f.min)
    return(area)
}

hit and miss方法收敛太慢,效率并不高,通常所说的MC积分是指下面这个方法。

按照Riemann对积分的定义: $$ \int{a}^{b}f(x)\,dx = \sum{i=0}^{n-1} f(x_i) \Delta{x} $$ 如果把f(x)在[a,b]区间上做N等分,那么

$$ \Delta{x} = (b-a)/N $$

这样子任一x值, 落入$ x_i + \Delta{x} $的概率就等于1/N,那么f(x)的期望值就是:

$$E(f(x)) = \sum_{i=0}^{n-1} f(x_i) P(X=xi) = \frac{\sum{i=0}^{n-1} f(x_i)}{N}$$

继而可以推导出, $$ \int_{a}^{b}f(x)dx = (b-a) \cdot E(f(x)) $$

底乘以平均高度,这跟算矩阵面积没啥区别,实现起来也非常容易。

mc.integrate <- function(fun, lower, upper, n=100) {
    x <- runif(n, lower, upper)
    y <- fun(x)
    fmean <- mean(y)
    area <- fmean * (upper-lower)
    return(area)
}

以正态分布为例,hit_miss精度较低,而上面这个方法,误差非常小。

> hit_miss(dnorm, 1,3)
[1] 0.1437858
> mc.integrate(dnorm, 1,3)
[1] 0.1578009
> integrate(dnorm, 1,3)
0.1573054 with absolute error < 1.7e-15

simpson法比起来,这个实现要简单得多。 当然在1维积分上,Monte Carlo的表现比不上simpson法,但是在高维积分上,Monte Carlo效率就很高。

比如求单位球体积,在xy平面的半球方程是 $$f(x,y) = \sqrt{r^2 - x^2 -y^2} $$

> r <- 1 
> n <- 10000
> set.seed=456
> x <- runif(n, -r, r) 
> y <- runif(n, -r, r) 
> 
> f2 <- r^2 - x^2 - y^2
> f=sqrt(f2[which(f2 > 0)])
> fmean = mean(f)
> I <- 2*pi*r^2*fmean
> I
[1] 4.182649
> 4*pi/3
[1] 4.18879

结果和真实值非常接近,而实现起来和一维积分一样简单。

参考资料:Monte Carlo Integration: http://beam.acclab.helsinki.fi/~knordlun/mc/mc5nc.pdf