如何通过Python实现蒙特卡罗模拟算法
本文主要介绍蒙特卡罗模拟算法，以及如何通过Python来模拟问题。
什么是蒙特卡罗（Monte Carlo）方法？
蒙特卡罗（Monte Carlo）方法，又称随机抽样或统计试验方法，是通过使用随机数（或更常见的伪随机数）来解决很多计算问题的方法，将所求解的问题同一定的概率模型相联系，用计算机实现统计模拟或抽样，以获得问题的近似解。
蒙特卡罗解题归结为三个主要步骤：
- 构造或描述概率过程；
- 实现从已知概率分布抽样；
- 建立各种估计量。
接下来我们介绍3个简单的案例，看一下如何在实际问题中应用这3个步骤进行求解。
案例1:
的计算
如何使用蒙特卡罗方法计算圆周率
？
正方形内部有一个相切的圆，假设圆形的半径为r ，那么正方形的边长则为2r ，因此，圆形和正方形的面积之比为：
现在，在正方形内部，随机产生一系列的点，计算他们与中心点的距离，从而判断是否落在圆的内部，从而：
即可估算出圆周率
的大小。
Python模拟计算
Step1 构造或描述概率过程
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
r = 1 # 半径
a, b = (0.0, 0.0) # 圆心
# 正方形边界
xmin, xmax = a-1, a+1
ymin, ymax = b-1, b+1Step2 实现从已知概率分布抽样
由于每个点落入不同位置的概率是一致的，所以我们使用均匀分布进行抽样。
# 随机抽取n个点
n = 10000
# 在正方形范围内，随机抽取n个点
x = np.random.uniform(xmin, xmax, n)
y = np.random.uniform(ymin, ymax, n)Step3 建立各种估计量
# 计算每个点到圆心的距离
d = np.sqrt((x - a) ** 2 + (y - b) ** 2)
# 圆内点的个数
res = sum(np.where(d<r, 1, 0)) # 距离小于半径即为圆内
# pi
pi = res / n * 4模拟结果：
3.1572可以看到，模拟结果和
还是比较接近的，我们可以通过增大样本点的个数n ，使得这个估计结果更加精确。
把所有点画出来：
# 绘制正方形
fig = plt.figure(figsize=(10, 10))
axes = fig.add_subplot(1, 1, 1)
plt.plot(x, y, 'ro', markersize=1)
plt.axis('equal')
# 绘制圆形
circle = Circle(xy=(a, b), radius=r, alpha=0.5, color='gray')
axes.add_patch(circle)
plt.grid(True, linestyle='--', linewidth=0.8)
案例2: 定积分的计算
问题描述
计算
在[0, 1]区间的定积分，也就是下图中的阴影面积大小。
按照蒙特卡罗模拟的思想，我们可以计算有多少点落在积分范围内（判断条件高度
），落在阴影范围内的点数跟所有抽样点数的比值就是所要求的积分值。
Python模拟计算
Step1 构造或描述概率过程
# 正方形边界
xmin, xmax = 0.0, 1.0
ymin, ymax = 0.0, 1.0
# 函数曲线
def f(x):
return x ** 2Step2 实现从已知概率分布抽样
由于每个点落入不同位置的概率是一致的，所以我们使用均匀分布进行抽样。
# 随机抽取n个点
n = 10000
# 在正方形范围内，随机抽取n个点(均匀分布)
x = np.random.uniform(xmin, xmax, n)
y = np.random.uniform(ymin, ymax, n)Step3 建立各种估计量
# 统计落在函数图像下方点的数目
res = sum(np.where(y < f(x), 1, 0))
# 计算积分的近似值
integral = res / n模拟结果：
0.3302我们实际计算一下，跟模拟结果进行对比：
可以看到模拟结果还是相对准确的。
把所有点画出来：
# 绘制正方形
fig = plt.figure(figsize=(10, 10))
axes = fig.add_subplot(1, 1, 1)
plt.plot(x, y, 'ro', markersize=1)
plt.axis('equal')
# 绘制函数曲线
xi = np.linspace(0, 1, 100)
yi = xi ** 2
plt.plot(xi, yi, '--k')
plt.fill_between(xi, yi, 0, color='gray', alpha=0.5, label='area') # 填充曲线下的区域
案例3: 排队上厕所问题
问题描述
一场电影结束后会有20个人想要上厕所，假设只有一个厕所，同时：
- 这20个人会在0-10分钟内全部到达厕所；
- 每个人上厕所的时间在1-3分钟之内；
求开始时间、结束时间、等待时间、上厕所时长和厕所空闲时间分布情况。
Python模拟
Step1 构造或描述概率过程
starting_time = [0] * num # 每个人的开始时间
finish_time = [0] * num # 每个人的结束时间
waiting_time = [0] * num # 每个人的等待时间
empty_time = [0] * num # 厕所的空闲时间Step2 实现从已知概率分布抽样
num = 20 # 人数
arriving_time = np.random.uniform(0, 10, size=num) # 到达厕所时长分布
arriving_time.sort() # 按照到达厕所时间从小到大排序
working_time = np.random.uniform(1, 3, size=num) # 上厕所时长分布Step3 建立各种估计量
# 初始化第一个数据
starting_time[0] = arriving_time[0]
finish_time[0] = starting_time[0] + working_time[0]
waiting_time[0] = starting_time[0] - arriving_time[0] # 0
for i in range(1, num):
last_people_finish_time = finish_time[i-1]
# 到达时间早于上一个人的结束时间
if arriving_time[i] < last_people_finish_time:
starting_time[i] = last_people_finish_time
waiting_time[i] = last_people_finish_time - arriving_time[i]
empty_time[i] = 0
# 到达时间晚于等于上一个人的结束时间
else:
starting_time[i] = arriving_time[i]
waiting_time[i] = 0
empty_time[i] = arriving_time[i] - last_people_finish_time
finish_time[i] = starting_time[i] + working_time[i]模拟结果可视化：
# 等待曲线
plt.plot(starting_time, label='starting')
plt.plot(finish_time, label='finish')
plt.plot(waiting_time, label='waiting')
plt.plot(empty_time, label='empty')
plt.plot(working_time, label='working')
plt.legend()
小结
本文首先介绍蒙特卡罗模拟算法的思想，及其求解问题的3个步骤。
接着，通过3个简单的案例讲解了如何使用Python实现蒙特卡罗模拟算法。
说明：本文问题来源于网易云课堂的数据分析师（python）课程。