用蒙特卡洛法求解三门问题和Π
一、三门问题
三门问题是一个经典的概率谜题：假设你参加一个游戏，有三扇关闭的门，其中一扇门后有汽车（中奖），另外两扇门后是山羊（不中奖）。规则如下：
-
你先随机选择一扇门（例如门 1）。
-
主持人知道哪扇门后有汽车，且一定会打开另一扇有山羊的门（例如打开门 2）。
-
此时你可以选择保持最初的选择，或切换到剩下的那扇未打开的门（例如门 3）。
问题：哪种策略（保持不变 / 切换选择）的中奖概率更高？
二、思路
-
模拟游戏过程：随机设定汽车所在的门、玩家初始选择、主持人打开的门。
-
两种策略对比：
-
策略 1：不换门，坚持最初的选择。
-
策略 2：换门，选择主持人未打开的另一扇门。
- 重复试验：通过大量模拟（如 100 万次），计算两种策略的中奖概率。
三、Python 实现代码
import numpy as np# 设置试验次数（越大结果越稳定）
n_trials = 1000000# 初始化两种策略的中奖次数
stick_wins = 0 # 不换门中奖次数
switch_wins = 0 # 换门中奖次数for _ in range(n_trials):# 1. 随机设定汽车所在的门（0,1,2分别代表三扇门）car_door = np.random.randint(0, 3)# 2. 玩家随机选择一扇门player_choice = np.random.randint(0, 3)# 3. 主持人打开一扇有山羊的门（既不是汽车门，也不是玩家选择的门）# 找出所有可能打开的门（排除汽车门和玩家选择的门）possible_doors = [door for door in [0,1,2] if door != car_door and door != player_choice]# 主持人随机打开其中一扇（当玩家选错时只有1种可能，选对时有2种可能）host_opens = np.random.choice(possible_doors)# 4. 计算两种策略的中奖情况# 策略1：不换门if player_choice == car_door:stick_wins += 1# 策略2：换门（选择剩下的未打开的门）switched_choice = [door for door in [0,1,2] if door != player_choice and door != host_opens][0]if switched_choice == car_door:switch_wins += 1# 计算概率
stick_prob = stick_wins / n_trials
switch_prob = switch_wins / n_trialsprint(f"试验次数：{n_trials}")
print(f"不换门中奖概率：{stick_prob:.4f}（约1/3）")
print(f"换门中奖概率：{switch_prob:.4f}（约2/3）")
试验次数：1000000不换门中奖概率：0.3335（约1/3）换门中奖概率：0.6665（约2/3）
数学结论：
-
初始选择正确的概率为 1/3，此时换门会失败。
-
初始选择错误的概率为 2/3，此时主持人打开另一扇错误的门后，换门一定能选中汽车，因此换门成功的概率为 2/3。
Π
import numpy as np
import matplotlib.pyplot as plt# 1. 设置参数
n_points = 100000 # 随机点数量（数量越多，结果越精确）# 2. 生成随机点（x和y的范围均为[-1,1]）
x = np.random.uniform(-1, 1, n_points) # 生成n_points个[-1,1)之间的随机数
y = np.random.uniform(-1, 1, n_points)# 3. 判断点是否落在单位圆内（x² + y² ≤ 1）
inside_circle = (x**2 + y**2) <= 1
m = np.sum(inside_circle) # 落在圆内的点数# 4. 计算Π的近似值
pi_approx = 4 * m / n_points
print(f"随机点数量：{n_points}")
print(f"落在圆内的点数：{m}")
print(f"Π的近似值：{pi_approx}")# 5. 可视化结果
plt.figure(figsize=(6, 6)) # 创建正方形画布
# 绘制落在圆内的点（蓝色）和圆外的点（红色）
plt.scatter(x[inside_circle], y[inside_circle], c='blue', s=1, alpha=0.5, label='圆内点')
plt.scatter(x[~inside_circle], y[~inside_circle], c='red', s=1, alpha=0.5, label='圆外点')
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.axis('equal') # 保证x和y轴比例一致
plt.title(f'蒙特卡洛法求Π（近似值：{pi_approx:.4f}）')
plt.legend()
plt.show()