蒙特卡洛方法以及python实现
- 1.什么是蒙特卡洛方法(Monte Carlo method)
- 2.蒙特卡洛方法的基本思想
- 3.应用： 蒙特卡洛求定积分常见方法
- 3.1投点法：
- 3.2期望法：
- 3.3蒙特卡洛求定积分
- 4.蒙特卡洛方法python实例
1.什么是蒙特卡洛方法(Monte Carlo method)
蒙特卡罗方法也称统计模拟方法，是1940年代中期由于科学技术的发展和电子计算机的发明，而提出的一种以概率统计理论为指导的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。
20世纪40年代，在冯·诺伊曼，斯塔尼斯拉夫·乌拉姆和尼古拉斯·梅特罗波利斯在洛斯阿拉莫斯国家实验室为核武器计划工作时，发明了蒙特卡罗方法。因为乌拉姆的叔叔经常在摩纳哥的蒙特卡洛赌场输钱得名，而蒙特卡罗方法正是以概率为基础的方法。
2.蒙特卡洛方法的基本思想
百度百科:
蒙特卡罗法也称统计模拟法、统计试验法。是把概率现象作为研究对象的数值模拟方法。是按抽样调查法求取统计值来推定未知特性量的计算方法。蒙特卡罗是摩纳哥的著名赌城，该法为表明其随机抽样的本质而命名。故适用于对离散系统进行计算仿真试验。在计算仿真中，通过构造一个和系统性能相近似的概率模型，并在数字计算机上进行随机试验，可以模拟系统的随机特性。
通常蒙特卡罗方法可以粗略地分成两类：一类是所求解的问题本身具有内在的随机性，借助计算机的运算能力可以直接模拟这种随机的过程。例如在核物理研究中，分析中子在反应堆中的传输过程。中子与原子核作用受到量子力学规律的制约，人们只能知道它们相互作用发生的概率，却无法准确获得中子与原子核作用时的位置以及裂变产生的新中子的行进速率和方向。科学家依据其概率进行随机抽样得到裂变位置、速度和方向，这样模拟大量中子的行为后，经过统计就能获得中子传输的范围，作为反应堆设计的依据。
另一种类型是所求解问题可以转化为某种随机分布的特征数，比如随机事件出现的概率，或者随机变量的期望值。通过随机抽样的方法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解。这种方法多用于求解复杂的多维积分问题
3.应用： 蒙特卡洛求定积分常见方法
3.1投点法：
这个方法也常常被用来求
π π
π
值。现在我们用它来求函数的定积分。如下图所示，有一个函数
f ( x ) {f(x)}
f
(
x
)
，若要求它从a到b的定积分，其实就是求曲线下方的面积。
这时我们可以用一个比较容易算得面积的矩型罩在函数的积分区间上（假设其面积为
A r e a Area
A
r
e
a
）。然后随机地向这个矩形框里面投点，其中落在函数
f ( x ) {f(x)}
f
(
x
)
下方的点为绿色，其它点为红色。然后统计绿色点的数量占所有点（红色+绿色）数量的比例为
r r
r
，那么就可以据此估算出函数
f ( x ) {f(x)}
f
(
x
)
从aa到bb的定积分为
A r e a × r Area×r
A
r
e
a
×
r
。
注意由蒙特卡洛法得出的值并不是一个精确值，而是一个近似值。而且当投点的数量越来越大时，这个近似值也越接近真实值。
3.2期望法：
下面我们来重点介绍一下利用蒙特卡洛法求定积分的第二种方法——期望法，有时也称为平均值法。详情移步该博主文章：https://blog.csdn.net/baimafujinji/article/details/53869358
3.3蒙特卡洛求定积分
蒙特卡洛方法的一个重要应用就是求定积分。来看下面的一个例子(参考文献1)
当我们在[a,b]之间随机取一点x时，它对应的函数值就是f(x)。接下来我们就可以用f(x) * (b - a)来粗略估计曲线下方的面积，也就是我们需要求的积分值，当然这种估计（或近似）是非常粗略的,如下图所示：
进行四次随机采样如下：
进行四次进行四次随机采样如下：
在此图中，做了四次随机采样，得到了四个随机样本
x 1 , x 2 , x 3 , x 4 {x_1},{x_2},{x_3},{x_4}
x
1
,
x
2
,
x
3
,
x
4
，并且得到了这四个样本的f(xi)f(xi)的值分别为
f ( x 1 ) , f ( x 2 ) , f ( x 3 ) , f ( x 4 ) f(x_1), f(x_2), f(x_3), f(x_4)
f
(
x
1
)
,
f
(
x
2
)
,
f
(
x
3
)
,
f
(
x
4
)
对于这四个样本，每个样本能求一个近似的面积值，大小为
f ( x i ) ∗ ( b − a ) f({x_i})*(b - a)
f
(
x
i
)
∗
(
b
−
a
)
。
为什么能这么干么？对照图下面那部分很容易理解，每个样本都是对原函数
f ( x ) f({x})
f
(
x
)
的近似，所以我们认为
f ( x ) f({x})
f
(
x
)
的值一直都等于
f ( x i ) f({x_i})
f
(
x
i
)
。
按照图中的提示，求出上述面积的数学期望，就完成了蒙特卡洛积分。
如果用数学公式表达上述过程：
S = 1 4 ( f ( x 1 ) ( b − a ) + f ( x 2 ) ( b − a ) + f ( x 3 ) ( b − a ) + f ( x 4 ) ( b − a ) ) = 1 4 ( b − a ) ( f ( x 1 ) + f ( x 2 ) + f ( x 3 ) + f ( x 4 ) ) = 1 4 ( b − a ) ∑ i = 1 4 f ( x i ) \begin{aligned} S & = \frac{1}{4}(f(x_1)(b-a) + f(x_2)(b-a) + f(x_3)(b-a) + f(x_4)(b-a)) \\ & = \frac{1}{4}(b-a)(f(x_1) + f(x_2) + f(x_3) + f(x_4)) \\ & = \frac{1}{4}(b-a) \sum_{i=1}^4 f(x_i) \end{aligned}
S
=
4
1
(
f
(
x
1
)
(
b
−
a
)
+
f
(
x
2
)
(
b
−
a
)
+
f
(
x
3
)
(
b
−
a
)
+
f
(
x
4
)
(
b
−
a
)
)
=
4
1
(
b
−
a
)
(
f
(
x
1
)
+
f
(
x
2
)
+
f
(
x
3
)
+
f
(
x
4
)
)
=
4
1
(
b
−
a
)
i
=
1
∑
4
f
(
x
i
)
对于更一般的情况，假设要计算的积分如下：
I = ∫ a b f ( x ) d x I = \int _a ^ b f(x) dx
I
=
∫
a
b
f
(
x
)
d
x
取N趋于无穷大，则有：
I ‾ = b − a N ∑ i = 1 N f ( X i ) \overline I = \frac{b-a}{N}\sum_{i=1}^Nf(X_i) I = N b − a i = 1 ∑ N f ( X i )
此时可以认为： I ‾ ≈ I \overline I \approx I I ≈ I
如果从直观上理解这个式子也非常简洁明了：
在[a,b]区间上按均匀分布取N个随机样本
x i {x_i}
x
i
,计算
f ( x 1 ) f(x_1)
f
(
x
1
)
并取均值，得到的相当于Y坐标值，然后乘以
( b − a ) (b-a)
(
b
−
a
)
为X坐标长度，得到的即为对应矩形的面积，即积分值。
4.蒙特卡洛方法python实例
首先看一个经典的用蒙特卡洛方法求π值。
简单介绍下思路：
正方形内部有一个相切的圆，它们的面积之比是π/4。现在，在这个正方形内部，随机产生n个点，计算它们与中心点的距离，并且判断是否落在圆的内部。若这些点均匀分布，则圆周率 pi=4 * count/n, 其中count表示落到圆内投点数 n:表示总的投点数。
//求π值
import
random
n
=
1000000
r
=
1.0
a
,
b
=
(
0.0
,
0.0
)
x_neg
,
x_pos
=
a
-
r
,
a
+
r
y_neg
,
y_pos
=
b
-
r
,
b
+
r
count
=
0.0
for
i
in
range
(
0
,
n
)
:
x
=
random
.
uniform
(
x_neg
,
x_pos
)
y
=
random
.
uniform
(
y_neg
,
y_pos
)
if
x
*
x
+
y
*
y
<=
1.0
:
count
+=
1
result
=
(
count
/
float
(
n
)
)
*
4
print
(
"result is "
,
result
)
运行结果：
result is 3.14036
再举个栗子
看一个求定积分的例子。
假设我们想求
∫ 0 1 x 2 d x \int_0^1 x^2 dx
∫
0
1
x
2
d
x
的值。具体代码如下:
（x*x > y，表示该点位于曲线的下面。所求的积分值即为曲线下方的面积与正方形面积的比。）
import
random
def
zero2onetox2
(
)
:
n
=
1000000
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
)
if
x
*
x
>=
y
:
count
+=
1
result
=
count
/
float
(
n
)
print
(
"result is "
,
result
)
然后运行：
zero2onetox2
(
)
运行结果：
result is 0.333084
由此可见，基于蒙特卡洛原理的求解值与解析值的结果误差还是比较小的。
参考文献：
此部分代码已上传github：https://github.com/ShuaiWang-Code/Monte_Carlo
1.https://blog.csdn.net/baimafujinji/article/details/53869358
2.http://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods-in-practice/monte-carlo-integration
3.https://blog.csdn.net/bitcarmanlee/article/details/82716641 参考1中的期望法
