蒙特卡洛模拟 – 极客岛

冯·诺伊曼 (John von Neumann)等三名科学家在 20 世纪 40 年代发明了蒙特卡洛模拟。他们以摩纳哥著名的赌城——蒙特卡洛 (Monte Carlo)——为其命名。

蒙特卡洛模拟 (Monte Carlo simulation)，也称统计模拟方法，是以概率统计理论为核心的数值计算方法。蒙特卡洛模拟可以提供多种可能的结果以及通过大量随机数据样本得出的每种结果的概率。

蒙特卡洛模拟是一种利用随机抽样来解决复杂问题的计算方法。它不依赖于精确的数学公式，而是通过模拟大量可能的结果，来揭示问题的本质。这种方法就像是在不确定性的海洋中，用随机的船桨划出一条通往答案的航道。

蒙特卡洛模拟的七个步骤

1. 定义问题：在开始之前，我们需要明确我们要探索的问题。它可能是金融市场的波动，工程设计的挑战，或是物理学中的未解之谜。

2. 构建模型：为问题构建一个数学模型，确定影响结果的关键变量和它们的概率分布。

3. 随机抽样：使用随机数生成器，为每个变量生成符合其概率分布的随机样本。

4. 运行模拟：将这些随机样本作为输入，运行模拟模型，观察输出结果。

5. 重复实验：为了提高结果的准确性，需要重复模拟过程多次。

6. 收集数据：记录每次模拟的结果，积累大量的数据。

7. 分析结果：通过统计分析，从数据中提取有用的信息，如均值、方差、置信区间等。

蒙特卡洛模拟的优势在于其灵活性和直观性，它适用于各种复杂系统和问题，无需解析解即可提供解决方案。然而，它也有局限性，如计算成本高、结果的不确定性以及对随机数生成器质量的依赖。

蒙特卡洛模拟应用举例—估计圆周率π的值

1. 定义问题：

我们想要估计圆周率π的值。圆的面积是π乘以半径的平方（A = πr²），而圆的外接正方形的面积是边长的平方（B = a² = 4r²）。根据蒙特卡洛模拟，单位圆面积与正方形面积的比是π/4，所以可以通过计算在单位圆内的点数与总点数的比例来估计π的值。

2. 构建模型：

定义一个模拟范围，在一个边长为2的正方形内，内切一个半径为1的圆。

3. 随机抽样：

在这个正方形内随机生成大量点（x, y），其中x和y的值都在[-1, 1]范围内。

4. 运行模拟：

对于每个点，检查它是否在圆内。如果点的坐标满足x²+𝑦²≤1，则认为它在圆内。

计算所有生成的点中有多少个在圆内。

使用圆内点数与总点数的比例来估计π的值。具体公式为：

π≈4×(圆内点数/总点数)

5. 重复实验或增大随机数数量：

为了提高估计的准确性，可以重复上述过程多次，然后取平均值。也可以不断增大随机数数量，提高估算精确度。

6. 收集数据：

在蒙特卡洛模拟中，数据收集主要是记录每次模拟的结果，即落在单位圆内的点的数量。这些数据将用于计算π的估计值。随着模拟次数的增加，我们可以收集到更多的数据点，这些数据点将帮助我们更好地理解模拟结果的分布和稳定性。通过收集这些数据，我们可以分析模拟结果的收敛性，并评估估计值的准确性。

7. 分析结果：

包括几个方面：

估计值的计算：根据落在圆内的点数与总点数的比例，使用公式 π≈4×(圆内点数/总点数) 来计算π的估计值。

收敛性分析：观察随着随机点数量的增加，π的估计值是否趋于稳定。通常，随着点数的增加，估计值应该越来越接近真实的π值。

误差评估：计算估计值与真实值之间的差异，即误差。可以通过多次模拟取平均值来减少随机误差。

置信区间：根据收集到的数据，可以构建估计值的置信区间，以评估结果的可靠性。

可视化：通过图表来可视化模拟结果，比如绘制估计值随点数增加的变化趋势图，或者显示随机点在单位圆内外的分布图。

Python示例代码：

# 导入了Numpy和Matplotlib库，分别用于数学计算和绘图。

import numpy as np
import matplotlib.pyplot as plt

# 生成一个500×2的矩阵X，其中的元素是随机分布在[-1,1]区间内的。然后将这个矩阵的列分别赋值给x和y，这样x和y就是500个随机点的横纵坐标。
X = np.random.uniform(-1, 1, size = (500,2))
x = X[:,0]
y = X[:,1]

#计算每个点(x,y)到原点的距离，如果这个距离小于1，则这个点在单位圆内，对应的masks值为True。
masks = np.sqrt(x**2 + y**2) < 1

#根据蒙特卡洛方法，单位圆面积与正方形面积的比是π/4，所以可以通过计算在单位圆内的点数与总点数的比例来估计π的值。
pi_est = 4 * sum(masks)/len(x)

#创建一个绘图窗口和坐标系。
fig, ax = plt.subplots()

# 在坐标系中添加一个单位圆。

circ = plt.Circle((0, 0), radius=1, edgecolor=‘k’, facecolor=‘None’)
ax.add_patch(circ)

#将落在单位圆内的点用蓝色“x”标记，落在单位圆外的点用红色“.”标记。设置坐标轴比例相等，并添加标题显示估计的π值。

# plot data inside the circle
plt.scatter(x[masks], y[masks], marker=“x”, alpha=0.5, color = ‘b’)

# plot data outside the circle
plt.scatter(x[~masks], y[~masks], marker=“.”, alpha=0.5, color = ‘r’)
plt.axis(‘scaled’)
plt.title(‘Estimated $pi$ = %1.3f’ %(pi_est))
plt.xlim(-1, 1)
plt.ylim(-1, 1)

# 定义一个函数est_pi，用于根据输入的点数n估计π的值。
# define a function of estimating pi
def est_pi(n):
    X = np.random.uniform(-1, 1, size = (int(n),2))
    x = X[:,0]
    y = X[:,1]

masks = np.sqrt(x**2 + y**2) < 1

pi_est = 4 * sum(masks)/len(x)

return pi_est

# 创建一个数组n_array，包含从1000到100000的100个点数。创建一个空数组est_pi_array，用于存储对应的π估计值。

n_array = np.linspace(1000,1000*100,100)
est_pi_array = np.empty(len(n_array))

# 遍历n_array中的每个点数，调用est_pi函数估计π的值，并将结果存储在est_pi_array中。

# convergence of estimated pi
i = 0
for n in n_array:
    pi_est = est_pi(n)
    est_pi_array[i] = pi_est
    i = i + 1

# 创建一个新的绘图窗口，并以对数坐标轴绘制点数与π估计值的关系。添加一条红色水平线表示真实的π值。显示网格，并最终展示图形。
fig, ax = plt.subplots()

plt.semilogx(n_array, est_pi_array)
plt.xlabel(“Number of random number”)
plt.ylabel(“Estimated $pi$”)
plt.axhline(np.pi, color=“r”);

plt.grid(True, which=“both”, ls=“–“)

plt.show()

下图所示为一次随机数数量为 500 条件下，圆周率估算的结果。

下图所示为不断增大随机数数量，圆周率估算精确度不断提高。

发表评论 取消回复

发表评论取消回复