其他
MC积分
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对积分的定义：
如果把f(x)在[a,b]区间上做N等分，那么
这样子任一x值， 落入
的概率就等于1/N，那么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平面的半球方程是
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
往期精彩