数学建模：蒙特卡罗法

1.蒙特卡罗法

蒙特卡罗法又称统计模拟法，是一种随机模拟方法，以概率和统计理论为基础的一种计算方法。通过将所求解的问题同一定的概率模型相联系，以获得问题的近似解。

2.模型原理

由大数定理可知，当样本容量足够大时，事件的发生频率即为其概率。（蒙特卡罗法准确来说是一种思想，没有通用代码）。

3.案例演示

3.1 求圆周率
即用落在圆内的点除以落在正方形里的点再乘以4来得到圆周率。

import numpy as np
import matplotlib.pyplot as plt

# 参赛初始化，投放10000个点，圆半径为1，圆心坐标为（1,1）
p=10000
r=1
x0,y0=1,1
n=0  # 初始时还未投放点，有0个点在圆内

# 设置绘图窗口
plt.figure()  # 配置一个绘图窗口
plt.title('Monte Carlo Simulation for Estimating Pi')
plt.xlabel('x')
plt.ylabel('y')

# 保持绘图窗口，多次绘图
for i in range(p):
    # p是在正方形内的随机点
    px=np.random.rand()*2  # 随机生成该点的横坐标
    # 生成的随机数乘以2，就可以把范围从0-1扩展到0-2
    py=np.random.rand()*2  # 随机生成该点的纵坐标
    # 判断点是否在圆内
    if(px-x0)**2+(py-y0)**2<r**2:  # 横纵坐标的平方和小于半径的平方，则在圆内
        plt.plot(px,py,'.',color='b')  # 绘制圆内的点，用蓝色表示
        n+=1
    else:
        plt.plot(px,py,'.',color='r')  # 绘制圆外的点，用红色表示
plt.axis('equal')  # 绘制时横纵坐标单位长度相同，便于观察圆
plt.show()

# 计算n的估计值
# n/p近似等于圆面积与正方形面积的比例，即 π/4
# 为了得到圆面积占正方形面积的比例的近似值，需将n/p乘以4
s=(n/p)*4  # 计算圆面积占正方形面积的比例
pi_estimate=s
print(f'π的估计值为：{pi_estimate}')

输出如下，可知 π 的估计值为 3.1304，与 3.1415 具有0.0111的误差，若投放的随机点越多，则误差越小直到可忽略不计：

3.2 三门问题
接下来我们来看一个很有意思的题：三门问题。
题目是：你参加一档电视节目，节目组提供了A、B、C三扇门，主持人告诉你，其中一扇门后面有辆汽车，其他两扇门后面是一头山羊，你可以选择一扇门打开获得门后的东西。假如你选择了B门，这时，主持人打开了C门，让你看到C门后面是一只山羊，然后问你要不要改选A门？（你想要汽车）

一般人的直觉：
在三个门，参赛者选一个，主持人排除掉一个错误的，还剩一个。那么此时的命题变为：

• 条件：两个门里，一个是汽车，一个是山羊，任选一个。
• 问题：交换选择是否能提高胜率？
• 答案：不能，二选一的问题，选哪个门都是50%的胜率。

然而事实上并不是这样。接下来我们来通过代码计算交换选择是否能提高胜率：

import numpy as np

# 在成功的条件下的概率
n=100000  # n代表蒙特卡罗模拟重复次数
a=0  # a表示不改变主意时能赢得汽车的次数
b=0  # b表示改变主意时能赢得汽车的次数
for i in range(n):  # 开始模拟n次
    x=np.random.randint(1,4)  # 随机生成一个1-3之间的整数x表示汽车出现在第x扇门后
    y=np.random.randint(1,4)  # 随机生成1-3之间的整数y表示自己选的门
    # 下面分为两种情况讨论，x=y和x!=y
    if x==y:  # 如果x和y相同，那么我们只有不改变主意时才能赢
        a+=1
    else:  # 如果x和y不同，那么我们只有改变主意时才能赢
        b+=1

print('蒙特卡罗方法得到的不改变主意时的获奖概率为：',a/n)
print('蒙特卡罗方法得到的改变主意时的获奖概率为：',b/n)

由输出结果可知改变主意的获奖概率竟大于不改变主意的获奖概率，甚至为不改主意时的2倍，这是为什么呢？ 其实这可以理解为一个条件概率问题。

① 若参赛者开始选中的是后面有汽车的门B（p=1/3），主持人打开后面有山羊的门C。这时参赛者如果换门，那么选中汽车的概率为1/3*0=0；如果不换门，那么选中汽车的概率为1/3*1=1/3 。

② 若参赛者开始选中的是后面有山羊的门A（p=1/3），主持人打开后面有山羊的门C。这时参赛者如果换门，那么选中汽车的概率为1/3*1=1；如果不换门，那么选中汽车的概率为1/3*0=0 。

③ 若参赛者开始选中的是后面有山羊的门C（p=1/3），主持人打开后面有山羊的门A。这时参赛者如果换门，那么选中汽车的概率为1/3*1=1；如果不换门，那么选中汽车的概率为1/3*0=0 。

于是可用下图表示为：

 综合上述三种情况，参赛者如果选择换门，那么选中汽车的概率为1/3+1/3=2/3；如果选择不换门，选中汽车的概率为1/3 。

而会出现这两种结果的关键在于：只有当参赛者开始选中后面有汽车的门时，在主持人打开后面有山羊的一扇门后，参赛者选择不换门，才会选中汽车；而参赛者在开始选中后面有汽车的门的概率仅为1/3，开始选中后面有山羊的门的概率为2/3，换门后，则把两种结果出现的概率扭转了过来。

代码实现：

import numpy as np
# 考虑失败情况的代码（无条件概率）
n=100000  # n代表蒙特卡罗模拟重复次数
a=0  # a表示不改变主意时能赢得汽车的次数
b=0  # b表示改变主意时能赢得汽车的次数
c=0  # c表示没有获奖的次数
for i in range(n):  # 开始模拟n次
    x=np.random.randint(1,4)  # 随机生成一个1-3之间的整数x表示汽车出现在第x扇门后
    y=np.random.randint(1,4)  # 随机生成1-3之间的整数y表示自己选的门
    change=np.random.randint(0,2)  # 0不改变主意，1改变主意
    # 下面分为两种情况讨论，x=y和x!=y
    if x==y:  # 如果x和y相同，那么我们只有不改变主意时才能赢
        if change==0:
            a+=1
        else:
            c+=1
    else:  # 如果x和y不同，那么我们只有改变主意时才能赢
        if change==1:
            b+=1
        else:
            c+=1  

print('蒙特卡罗方法得到的不改变注意时的获奖概率为：',a/n)
print('蒙特卡罗方法得到的改变注意时的获奖概率为：',b/n)
print('蒙特卡罗方法得到的没有获奖的概率为：',c/n)

所以，综上所述，答案是换门比不换门能提高胜率 。