蒙特卡洛方法(MC)
一、蒙特卡洛方法简介
蒙特卡洛方法得名于摩纳哥的蒙特卡洛赌场,是大数定律的经典应用,即用大样本数据计算出来的频率估计概率,采样越多,越近似最优解,但永远不是最优解。
蒙特卡罗方法可以粗略地分成两类:一类是所求解的问题本身具有内在的随机性,可以借助计算机进行模拟,例如在物理中,对光子传输的模拟、分子运动的模拟以及犹豫不决,量子力学;或者经济学上,模拟市场走势等。
另一类是所求解的问题可以转化为概率分布问题,通过统计实验结果来求解。例如求圆周率,求定积分。
二、利用蒙特卡罗方法计算圆周率
代码:
N = 1000
x = matrix(runif(N*2), N, 2)
plot(x)
L = x[,1]^2 + x[,2]^2 < 1
points(x[L,], col = 'red')
n = sum(L)
p = 4*n/N
没什么高深的内核 咱就是说暴力美学出奇迹(bushi)
三、利用蒙特卡洛方法求定积分
首先,需要了解连续随机变量X的期望定义为:
E
(
X
)
=
∫
−
∞
+
∞
x
f
(
x
)
d
x
E(X) = \int_{- \infty}^{+ \infty} xf(x) dx
E(X)=∫−∞+∞xf(x)dx
连续变量X的任意函数g(x)也是一个随机变量,于是下式成立:
E
(
X
)
=
∫
−
∞
+
∞
g
(
x
)
f
(
x
)
d
x
E(X) = \int_{- \infty}^{+ \infty} g(x)f(x) dx
E(X)=∫−∞+∞g(x)f(x)dx
利用蒙特卡罗方法求定积分的主要思想在于,将积分转换为期望的形式,再利用矩估计,用样本均值代替总体期望,样本均值也就近似为定积分结果:
∫
−
∞
+
∞
g
(
x
)
d
x
=
1
f
(
x
)
∫
−
∞
+
∞
g
(
x
)
f
(
x
)
d
x
=
1
f
(
x
)
E
(
X
)
\int_{- \infty}^{+ \infty} g(x)dx=\frac{1}{f(x)}\int_{- \infty}^{+ \infty} g(x)f(x) dx=\frac{1}{f(x)}E(X)
∫−∞+∞g(x)dx=f(x)1∫−∞+∞g(x)f(x)dx=f(x)1E(X)
可以通过以下几个例子,加深理解。
例1
计算下列定积分
∫
0
1
x
d
x
\int_{0}^{1} x dx
∫01xdx
手工推导:
代码展示:
# MC estimator
n = 3000
x = runif(n) # 生成n个0到1之间服从均匀分布的随机数
g = x # g(x)
mean(g) # 矩估计:样本均值代替总体期望
例2
计算下列定积分
∫
0
1
e
−
x
d
x
\int_{0}^{1} e^{-x} dx
∫01e−xdx
手工推导:
代码展示:
n = 3000
x = runif(n)
g = exp(1)^(-x)
mean(g)
1-1/exp(1) # 解析解
例3
计算下列定积分
∫
4
8
x
2
d
x
\int_{4}^{8}x^{2} dx
∫48x2dx
手工推导:
代码展示:
n = 3000
x = runif(n,4,8)
g = x**2
(8-4)*mean(g)
总结:
经过上面三个练习,可以将原理推广为:
∫
a
b
g
(
x
)
d
x
=
(
b
−
a
)
∫
−
∞
+
∞
g
(
x
)
f
(
x
)
d
x
=
(
b
−
a
)
E
(
X
)
\int_{a}^{b} g(x)dx=(b-a )\int_{- \infty}^{+ \infty} g(x)f(x) dx=(b-a)E(X)
∫abg(x)dx=(b−a)∫−∞+∞g(x)f(x)dx=(b−a)E(X)
但是该方法存在一个极大的缺点就是,计算效率低!!!但是耐不住它效果好,与解析解误差较小,所以应用依旧很广泛。