蒙特卡洛方法的理解、推导和应用

最新推荐文章于 2024-08-21 19:59:48 发布

氢键H-H

最新推荐文章于 2024-08-21 19:59:48 发布

阅读量2.6w

点赞数 28

分类专栏：算法文章标签：蒙特卡洛方法

本文为氢键H-H原创文章，未经允许请勿用于商业用途，转载请注明出处

本文链接：https://blog.csdn.net/qq_32618327/article/details/90452846

版权

算法专栏收录该内容

18 篇文章 43 订阅

订阅专栏

蒙特卡洛方法的理解、推导和应用

1. 简介
2. 基本思想
3. 蒙特卡洛法求定积分
- 3.1. 随机投点法
- 3.2. 平均值法
4. 理解平均值法
5. 应用
- 5.1 蒙特卡洛方法求 π 值
- 5.2 蒙特卡洛方法求定积分

1. 简介

蒙特卡罗方法也称 统计模拟 方法
1940年代中期由于科学技术的发展和电子计算机的发明
而提出的一种以概率统计理论为指导的数值计算方法
是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法

20世纪40年代，在冯·诺伊曼，斯塔尼斯拉夫·乌拉姆和尼古拉斯·梅特罗波利斯在洛斯阿拉莫斯国家实验室为核武器计划工作时，发明了蒙特卡罗方法
因为乌拉姆的叔叔经常在摩纳哥的蒙特卡洛赌场输钱得名，而蒙特卡罗方法正是以概率为基础的方法

2. 基本思想

蒙特卡罗法也称统计模拟法、统计试验法，是把概率现象作为研究对象的数值模拟方法
是按抽样调查法求取统计值来推定未知特性量的计算方法

蒙特卡罗是摩纳哥的著名赌城，该法为表明其随机抽样的本质而命名
故适用于对离散系统进行计算仿真试验
在计算仿真中，通过构造一个和系统性能相近似的概率模型
并在数字计算机上进行随机试验，可以模拟系统的随机特性

通常蒙特卡罗方法可以粗略地分成两类：

所求解的问题本身具有内在的随机性，借助计算机的运算能力可以直接模拟随机的过程
例如在核物理研究中，分析中子在反应堆中的传输过程
中子与原子核作用受到量子力学规律的制约，人们只能知道它们相互作用发生的概率，却无法准确获得中子与原子核作用时的位置以及裂变产生的新中子的行进速率和方向
科学家依据其概率进行随机抽样得到裂变位置、速度和方向，这样模拟大量中子的行为后，经过统计就能获得中子传输的范围，作为反应堆设计的依据
所求解问题可以转化为某种随机分布的特征数
比如随机事件出现的概率，或者随机变量的期望值
通过随机抽样的方法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解
这种方法多用于求解复杂的多维积分问题

3. 蒙特卡洛法求定积分

3.1. 随机投点法

这个方法也常常被用来求 $π$ 值，现在用它来求函数的定积分
如下图所示，有一个函数 $f (x)$ ，若要求它从 a 到 b 的定积分，其实就是求曲线下方的面积

在这里插入图片描述
这时可以用一个比较容易算得面积的矩型罩在函数的积分区间上（假设其面积为 $A r e a$ ）
然后随机地向这个矩形框里面投点，其中落在函数 $f (x)$ 下方的点为绿色，其它点为红色
然后统计绿色点的数量占所有点（红色+绿色）数量的比例为 $r$
那么就可以据此估算出函数 $f (x)$ 从 a 到 b 的定积分为 $A r e a$ × r

注意由蒙特卡洛法得出的值并不是一个精确值，而是一个近似值
而且当投点的数量越来越大时，这个近似值也越接近真实值

3.2. 平均值法

利用蒙特卡洛法求定积分的第二种方法——平均值法，有时也称为期望法
详情移步博主：白马负金羁的文章：蒙特卡洛（Monte Carlo）法求定积分

4. 理解平均值法

来看下面的一个例子
在这里插入图片描述

当在[a,b]之间随机取一点x时，它对应的函数值就是 $f (x)$
接下来就可以用 $f (x) * (b - a)$ 来粗略估计曲线下方的面积，也就是需要求的积分值
当然这种估计（或近似）是非常粗略的,如下图所示：
在这里插入图片描述

所以，在[a,b]之间随机取一系列点 $x$ _i时（ $x$ _i 满足均匀分布）
把估算出来的面积取平均来作为积分估计的一个更好的近似值
如果这样的采样点越来越多，那么对于这个积分的估计也就越来越接近
进行四次随机采样（满足均匀分布）如下：
在这里插入图片描述
在此图中，做了四次随机采样，得到了四个随机样本 $x 1$ ， $x 2$ ， $x 3$ ， $x 4$
并且得到了这四个样本的 $f$ (x_i) 值分别为 $f$ (x₁)， $f$ (x₂)， $f$ (x₃)， $f$ (x₄)
对于这四个样本，每个样本能求一个近似的面积值，大小为 $f$ (x_i) $* (b - a)$

对照图下面那部分很容易理解，每个样本都是对原函数 $f (x)$ 的近似
所以认为 $f (x)$ 的值一直都等于 $f$ (x_i)

按照图中的提示，求出上述面积的数学期望，就完成了蒙特卡洛积分
如果用数学公式表达上述过程：
在这里插入图片描述
对于更一般的情况，假设要计算的积分如下：
取N趋于无穷大，则有：
就可以认为：

如果从直观上理解这个式子也非常简洁明了：
在 [a , b] 区间上按均匀分布取N个随机样本 $x$ _i ，计算 $f$ (x_i) 并取均值，得到的相当于Y坐标值
然后乘以 $(b - a)$ 为X坐标长度，得到的即为对应矩形的面积，即积分值

5. 应用

5.1 蒙特卡洛方法求 π 值

正方形内部有一个相切的圆，它们的面积之比是 $π$ / 4
在这里插入图片描述
现在，在这个正方形内部，随机产生n个点
计算它们与中心点的距离，并且判断是否落在圆的内部

若这些点均匀分布，令 count 表示落到圆内投点数，n 表示总投点数
则圆周率 pi = 4 * count / n

通过随机模拟300000个点， $π$ 的估算值与真实值相差0.01%

import random

r = 1.0  # 圆半径，假设为1
n = 300000  # 总投点数
count = 0  # 落到圆内投点数

# 投点x,y的范围
x_min, x_max = -r, r
y_min, y_max = -r, r

for i in range(0, n):

    # 在 [min, max] 范围内随机生成实数
    x = random.uniform(x_min, x_max)
    y = random.uniform(y_min, y_max)

    # 落到圆内投点数+1
    if x * x + y * y <= r:
        count += 1

pi = (count / float(n)) * 4
print("pi is ", pi)

# pi is 3.142053333333333

5.2 蒙特卡洛方法求定积分

假设求 $\int_{1}^{2} x^3\, dx$ 的值，其估算值与真实值相差0.7%
在这里插入图片描述

import random

n = 300000  # 总投点数
count = 0  # 位于曲线之下的投点数

# 投点x,y的范围
x_min, x_max = 1.0, 2.0
y_min, y_max = 0.0, 8.0

for i in range(0, n):

    # 在[min, max]范围内随机生成实数
    x = random.uniform(x_min, x_max)
    y = random.uniform(y_min, y_max)

    # 位于曲线之下的投点数+1
    if x * x * x >= y:
        count += 1

# 所求的积分值即为曲线下方的面积与正方形面积8的比
result = count / float(n) * 8
print("result is ", result)

# result is  3.7480533333333335