用Python玩转蒙特卡罗：从骰子模拟到圆周率计算（附完整代码）

用Python玩转蒙特卡罗：从骰子模拟到圆周率计算（附完整代码）

如果你曾经好奇，为什么扔一万次骰子后，每个点数出现的次数会趋于平均？或者，如何用“撒豆子”这种看似儿戏的方法，去计算一个像圆周率π这样精确到万亿位的数学常数？这背后，都藏着一个强大而迷人的思想——蒙特卡罗方法。它不是什么高深莫测的数学魔法，而是一种用“随机性”去解决“确定性”问题的智慧。对于Python初学者和数据分析爱好者来说，蒙特卡罗就像一把万能钥匙，它能将复杂的积分、概率预测甚至金融风险评估，都转化为我们熟悉的随机模拟游戏。这篇文章，我将带你亲手用Python和几个核心库，从最生活化的骰子、硬币案例开始，一步步构建属于你自己的蒙特卡罗模拟世界，最终我们会一起“算”出圆周率，并理解其背后的统计原理。整个过程，代码就是我们的实验工具，可视化图表则是我们观察随机世界规律的窗口。

1. 蒙特卡罗初体验：从生活游戏理解“随机逼近”

很多人第一次听说蒙特卡罗，会觉得它和赌博有关。确实，它的名字来源于摩纳哥著名的赌城蒙特卡洛，因为其核心的随机抽样思想，与轮盘赌的不可预测性神似。但它的本质，是一种统计模拟方法。简单来说，当你遇到一个很难直接计算答案的问题时，可以设计一个随机实验，通过大量重复这个实验，用实验结果出现的频率去逼近我们想求的概率或数值。

想想我们小时候玩的抛硬币。理论上，正面朝上的概率是50%。但如果你只抛10次，很可能出现7次正面、3次反面的情况，这并不能很好地反映理论概率。但如果你有耐心抛上一万次、十万次，那么正面出现的比例会稳定在0.5附近。这个“频率稳定于概率”的现象，就是蒙特卡罗方法赖以生存的基石——大数定律。

让我们用Python来亲手验证一下。我们将模拟抛一枚均匀硬币的过程。

import random
import matplotlib.pyplot as plt

def coin_toss_experiment(num_tosses):
    """模拟抛硬币实验，返回正面朝上的频率列表"""
    heads_count = 0
    frequency_list = []
    
    for toss in range(1, num_tosses + 1):
        # random.choice 随机选择 'H'（正面）或 'T'（反面）
        if random.choice(['H', 'T']) == 'H':
            heads_count += 1
        current_frequency = heads_count / toss
        frequency_list.append(current_frequency)
    
    return frequency_list

# 进行10000次抛硬币模拟
num_simulations = 10000
frequencies = coin_toss_experiment(num_simulations)

# 绘制频率变化图
plt.figure(figsize=(10, 6))
plt.plot(range(1, num_simulations + 1), frequencies, linewidth=0.5, alpha=0.7, label='正面出现频率')
plt.axhline(y=0.5, color='r', linestyle='--', linewidth=2, label='理论概率 (0.5)')
plt.xlabel('抛硬币次数')
plt.ylabel('正面出现频率')
plt.title('大数定律的直观演示：抛硬币频率趋于稳定')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

运行这段代码，你会看到一幅生动的图表：起初，频率曲线剧烈波动，但随着实验次数增加，它像被一只无形的手安抚，越来越贴近那条红色的0.5理论线。这就是蒙特卡罗的“魔法”开始的地方——用大量随机样本的统计特性，去揭示系统的内在规律。

提示：在蒙特卡罗模拟中，random模块生成的其实是“伪随机数”。它们由确定的算法产生，但具有足够好的统计随机性，对于我们的模拟来说完全够用。在更严格的科学计算中，可能会用到numpy.random，它提供了更高效、更多样的分布生成器。

理解了基础思想后，我们可以把游戏升级。假设一个袋子里有10个球，编号1-10，我们进行有放回抽取。连续3次都抽到偶数号的概率是多少？用数学算，是(5/10)^3 = 0.125。用蒙特卡罗模拟，我们只需要重复这个“抽三次球”的实验很多遍，然后数一数其中“三次皆偶数”的次数占比。

def even_ball_experiment(num_trials):
    """模拟从10个球中抽取3次，记录三次均为偶数的频率"""
    success_count = 0
    
    for _ in range(num_trials):
        # 模拟三次有放回抽取
        draws = [random.randint(1, 10) for _ in range(3)]
        # 检查是否全部为偶数
        if all(draw % 2 == 0 for draw in draws):
            success_count += 1
    
    return success_count / num_trials

# 尝试不同的模拟次数
trial_counts = [100, 1000, 5000, 10000, 50000]
estimates = []
for n in trial_counts:
    est = even_ball_experiment(n)
    estimates.append(est)
    print(f"模拟次数 {n:6d} | 估计概率: {est:.6f} | 与理论值(0.125)的绝对误差: {abs(est - 0.125):.6f}")

你会发现，模拟次数越多，估计值通常越接近0.125。但请注意，由于随机性，50000次模拟的结果不一定比10000次更接近理论值，但它的波动范围会更小，结果更可靠。这引出了蒙特卡罗方法中关于误差和效率的重要话题，我们稍后会深入探讨。

2. 核心武器库：NumPy的随机引擎与向量化加速

当我们从简单的单个实验转向需要成千上万次、甚至百万次模拟时，纯Python的循环会变得异常缓慢。这时，NumPy就该登场了。NumPy的核心优势在于向量化操作和高性能的随机数生成器。它允许我们一次性生成海量的随机样本，并用数组运算替代循环，效率提升可达数十甚至上百倍。

让我们用NumPy重写抛硬币实验，感受一下速度的飞跃。

import numpy as np
import time

def coin_toss_numpy(num_tosses):
    """使用NumPy向量化操作模拟抛硬币"""
    # 一次性生成所有抛掷结果：1代表正面，0代表反面
    # np.random.choice 比 random.choice 快得多
    tosses = np.random.choice([1, 0], size=num_tosses, p=[0.5, 0.5])
    # 计算累积正面数，然后除以累积次数得到频率
    cumulative_heads = np.cumsum(tosses)
    trial_numbers = np.arange(1, num_tosses + 1)
    frequencies = cumulative_heads / trial_numbers
    return frequencies

# 对比性能
num_tosses = 1000000

start = time.time()
frequencies_py = coin_toss_experiment(num_tosses)  # 使用之前纯Python的函数
end = time.time()
print(f"纯Python循环耗时: {end - start:.4f} 秒")

start = time.time()
frequencies_np = coin_toss_numpy(num_tosses)
end = time.time()
print(f"NumPy向量化耗时: {end - start:.4f} 秒")

在我的测试中，NumPy版本通常比纯Python循环快50倍以上。对于蒙特卡罗模拟这种计算密集型任务，这意味着一小时的计算可以缩短到一分钟内完成。

NumPy的random模块提供了丰富的概率分布，这是我们构建更复杂模型的工具箱。下面这个表格列出了最常用的几种分布及其在蒙特卡罗中的应用场景：

分布类型	NumPy生成函数	主要参数	典型应用场景
均匀分布	np.random.uniform(low, high, size)	下限low, 上限high	等概率随机投点、生成任意区间内的随机数
正态分布	np.random.normal(loc, scale, size)	均值loc, 标准差scale	模拟测量误差、股票收益率、自然现象
二项分布	np.random.binomial(n, p, size)	试验次数n, 成功概率p	模拟固定次数独立试验的成功次数（如抛硬币正面次数）
泊松分布	np.random.poisson(lam, size)	发生率lam	模拟单位时间内随机事件发生的次数（如客服接电话数）
指数分布	np.random.exponential(scale, size)	尺度参数scale (均值=scale)	模拟随机事件的时间间隔（如设备故障间隔）

例如，模拟一个骰子投掷6000次，并统计各点数的出现频率，用NumPy可以一行代码搞定：

# 模拟一个公平的六面骰子投掷6000次
num_rolls = 6000
rolls = np.random.randint(1, 7, size=num_rolls)  # 生成1到6之间的随机整数

# 统计频率
values, counts = np.unique(rolls, return_counts=True)
frequencies = counts / num_rolls

print("点数 : 频率")
for v, f in zip(values, frequencies):
    print(f"  {v}   : {f:.4f}")

理论上，每个点数的频率都应接近1/6 ≈ 0.1667。运行代码，你会发现结果非常接近。这种从大量随机样本中快速提取统计信息的能力，是进行高效蒙特卡罗分析的基础。

3. 可视化艺术：用Matplotlib洞察随机过程

蒙特卡罗模拟生成的是数据，而可视化是将这些数据转化为洞察力的关键。Matplotlib库能帮助我们直观地看到随机过程如何收敛，以及结果的分布形态。一个好的图表，有时比千言万语更能说明问题。

让我们结合前面学到的NumPy，创建一个更综合的可视化示例：模拟投掷两个骰子，并观察其点数之和的分布。

import matplotlib.pyplot as plt

def dice_sum_simulation(num_rolls):
    """模拟投掷两个骰子，记录点数之和"""
    # 一次性生成所有投掷结果
    die1 = np.random.randint(1, 7, size=num_rolls)
    die2 = np.random.randint(1, 7, size=num_rolls)
    sums = die1 + die2
    return sums

# 模拟10万次投掷
num_rolls = 100000
sums = dice_sum_simulation(num_rolls)

# 计算理论概率分布
# 两个骰子点数之和从2到12，共11种可能
possible_sums = np.arange(2, 13)
theoretical_probs = np.array([1,2,3,4,5,6,5,4,3,2,1]) / 36.0

# 计算模拟得到的频率
unique_sums, counts = np.unique(sums, return_counts=True)
simulated_probs = counts / num_rolls

# 创建对比图表
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 子图1：模拟结果的直方图
ax1.hist(sums, bins=np.arange(1.5, 13.5, 1), density=True, edgecolor='black', alpha=0.7, label='模拟频率')
ax1.plot(possible_sums, theoretical_probs, 'ro-', linewidth=2, markersize=8, label='理论概率')
ax1.set_xlabel('两点数之和')
ax1.set_ylabel('概率/频率')
ax1.set_title(f'两个骰子点数之和的分布 (模拟次数: {num_rolls:,})')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 子图2：模拟与理论的偏差
ax2.bar(possible_sums, simulated_probs - theoretical_probs, color='orange', edgecolor='black')
ax2.axhline(y=0, color='r', linestyle='--')
ax2.set_xlabel('两点数之和')
ax2.set_ylabel('频率 - 理论概率')
ax2.set_title('模拟频率与理论概率的偏差')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

这段代码会生成两张并排的图表。左边直方图清晰地展示了经典的概率分布形态——中间（和为7）的概率最高，向两边递减。红色的理论概率点完美地落在直方图柱子上方，表明模拟结果与理论高度吻合。右边的偏差图则量化了这种吻合程度，所有偏差值都紧密分布在0线附近，这证明了我们模拟的准确性。

注意：在绘制直方图时，density=True参数非常重要，它会让直方图归一化，使得总面积和为1，这样才能与概率密度进行直接比较。bins参数我们设置为从1.5到12.5，步长为1，这样每个柱子的中心正好对应整数和（2,3,...,12）。

可视化不仅能用于展示结果，还能动态展示收敛过程。还记得我们估算圆周率的想法吗？在进入具体代码前，我们可以先构思一下如何用动画来展示：随机点如何在正方形中逐渐撒开，圆内点的比例如何一步步逼近π/4。虽然本文不展开动画代码，但你可以尝试用matplotlib.animation模块来实现，它会让你对蒙特卡罗的“逐渐逼近”有更深刻的体会。

4. 经典项目实战：蒙特卡罗方法估算圆周率

现在，让我们挑战蒙特卡罗最著名的演示项目：估算圆周率π。这个项目的优雅之处在于，它将一个纯粹的几何问题，转化为了一个可以通过“撒点”解决的随机实验。

原理回顾： 在一个边长为2的正方形中，内切一个半径为1的圆。它们的面积分别为：

• 正方形面积：( S_{square} = (2)^2 = 4 )
• 圆形面积：( S_{circle} = \pi * (1)^2 = \pi )

如果我们在这个正方形内均匀随机地投点，那么一个点落在圆内的概率，就等于圆的面积与正方形的面积之比： ( P = \frac{S_{circle}}{S_{square}} = \frac{\pi}{4} )

因此，如果我们投了N个点，其中有M个落在圆内，那么频率 ( \frac{M}{N} ) 应该近似等于 ( \frac{\pi}{4} )。于是，π的估计值就是： ( \pi_{estimate} \approx 4 * \frac{M}{N} )

下面，我们用Python实现这个想法，并深入分析其误差和收敛速度。

def estimate_pi(num_points):
    """
    使用蒙特卡罗方法估计圆周率。
    
    参数:
        num_points (int): 随机生成的点数。
    
    返回:
        float: π的估计值。
        numpy.ndarray: 所有点的x坐标。
        numpy.ndarray: 所有点的y坐标。
        numpy.ndarray: 布尔数组，True表示点在圆内。
    """
    # 在[-1, 1]区间内均匀生成随机点
    x = np.random.uniform(-1, 1, num_points)
    y = np.random.uniform(-1, 1, num_points)
    
    # 计算每个点到原点(0,0)的距离
    distances = np.sqrt(x**2 + y**2)
    
    # 判断点是否在单位圆内 (距离 <= 1)
    inside_circle = distances <= 1.0
    
    # 计算圆内点的数量
    num_inside = np.sum(inside_circle)
    
    # 估计π值
    pi_estimate = 4 * num_inside / num_points
    
    return pi_estimate, x, y, inside_circle

# 进行不同规模的点数实验
point_counts = [100, 1000, 5000, 10000, 50000, 100000, 500000]
estimates = []
errors = []

print("点数\t\t估计的π值\t\t绝对误差\t\t相对误差(%)")
print("-" * 70)
for n in point_counts:
    pi_est, _, _, _ = estimate_pi(n)
    abs_error = abs(pi_est - np.pi)
    rel_error = (abs_error / np.pi) * 100
    estimates.append(pi_est)
    errors.append(abs_error)
    print(f"{n:8,}\t{pi_est:.6f}\t\t{abs_error:.6f}\t\t{rel_error:.4f}%")

运行这段代码，你会看到一个清晰的表格，展示随着点数增加，估计值如何变化。通常，10万个点就能将误差控制在千分之一以内。但请注意，由于随机性，50万个点的结果不一定比10万个点更精确，但它的误差波动范围会更小。

为了更直观地理解这个过程，我们为一次中等规模的模拟（比如1万个点）绘制散点图。

def plot_pi_estimation(x, y, inside, pi_estimate, num_points):
    """绘制蒙特卡罗估算π的可视化图"""
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
    
    # 子图1：散点图
    ax1.scatter(x[inside], y[inside], color='blue', s=1, alpha=0.6, label='圆内点')
    ax1.scatter(x[~inside], y[~inside], color='red', s=1, alpha=0.6, label='圆外点')
    # 绘制圆形和正方形边框
    circle = plt.Circle((0, 0), 1, color='green', fill=False, linewidth=2)
    ax1.add_patch(circle)
    ax1.plot([-1, 1, 1, -1, -1], [-1, -1, 1, 1, -1], 'k-', linewidth=2) # 正方形
    ax1.set_xlim(-1.1, 1.1)
    ax1.set_ylim(-1.1, 1.1)
    ax1.set_aspect('equal', 'box')
    ax1.set_title(f'蒙特卡罗法估算π (点数: {num_points:,})')
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')
    ax1.legend(loc='upper right', markerscale=6)
    ax1.grid(True, alpha=0.3)
    
    # 在图中添加估计值文本
    textstr = f'$N = {num_points:,}$\n$N_{{in}} = {np.sum(inside):,}$\n$\pi \\approx 4 * ({np.sum(inside):,}/{num_points:,})$\n$\\mathbf{{\pi \\approx {pi_estimate:.5f}}}$'
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.8)
    ax1.text(0.05, 0.95, textstr, transform=ax1.transAxes, fontsize=12,
             verticalalignment='top', bbox=props)
    
    # 子图2：收敛过程（假设我们记录了每增加1000个点时的估计值）
    # 为了展示收敛，我们重新模拟并记录中间过程
    cumulative_inside = np.cumsum(inside)
    indices = np.arange(1, num_points + 1)
    running_estimates = 4 * cumulative_inside / indices
    
    ax2.plot(indices, running_estimates, 'b-', alpha=0.7, linewidth=0.5)
    ax2.axhline(y=np.pi, color='r', linestyle='--', linewidth=2, label='真实π值')
    ax2.set_xlabel('已投点数')
    ax2.set_ylabel('π的估计值')
    ax2.set_title('估计值随点数增加的收敛过程')
    ax2.set_xscale('log')  # 使用对数坐标能更清晰地展示初期波动
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()

# 运行一次模拟并绘图
num_points = 10000
pi_est, x_vals, y_vals, inside_flags = estimate_pi(num_points)
plot_pi_estimation(x_vals, y_vals, inside_flags, pi_est, num_points)

这张图是理解蒙特卡罗精髓的绝佳材料。左图，你能看到随机点如何填满正方形，蓝点（圆内）和红点（圆外）的比例直观地决定了π的估计值。右图则展示了收敛路径：初期估计值剧烈震荡，随着点数增加，震荡幅度减小，逐渐向红色虚线（真实π值）靠拢。这种“震荡衰减”的模式，是蒙特卡罗方法收敛的典型特征。

误差分析与效率： 蒙特卡罗方法的误差有一个非常重要的性质：误差与点数的平方根成反比。用公式表示，标准差大约为 ( \sigma \approx \frac{C}{\sqrt{N}} )，其中C是一个常数。这意味着，要想将误差减小为原来的十分之一，你需要将模拟次数（点数）增加一百倍。这就是为什么我们有时需要模拟数百万甚至数十亿个点来获得高精度结果。

我们可以通过实验来验证这个规律。多次运行不同点数的模拟，计算估计值的标准差。

def analyze_pi_error(max_points=100000, step=5000, trials_per_point=50):
    """分析不同点数下π估计值的误差分布"""
    point_range = np.arange(step, max_points + step, step)
    std_errors = []
    
    for n in point_range:
        estimates = []
        for _ in range(trials_per_point):
            pi_est, _, _, _ = estimate_pi(n)
            estimates.append(pi_est)
        std_errors.append(np.std(estimates))
    
    # 拟合误差曲线： std_error = A / sqrt(n)
    # 取对数： log(std_error) = log(A) - 0.5 * log(n)
    log_n = np.log(point_range)
    log_std = np.log(std_errors)
    coeffs = np.polyfit(log_n, log_std, 1)  # 一次多项式拟合
    A_fit = np.exp(coeffs[1])  # 截距的指数就是A
    slope_fit = coeffs[0]  # 斜率，理论值应为-0.5
    
    # 绘图
    plt.figure(figsize=(10, 6))
    plt.plot(point_range, std_errors, 'bo-', label='模拟标准差')
    # 绘制拟合曲线
    fit_curve = A_fit / np.sqrt(point_range)
    plt.plot(point_range, fit_curve, 'r--', linewidth=2, 
             label=f'拟合: σ ≈ {A_fit:.3f}/√N\n斜率: {slope_fit:.3f} (理论: -0.5)')
    plt.xlabel('点数 N')
    plt.ylabel('估计值的标准差 σ')
    plt.title(f'蒙特卡罗π估计的误差分析 (每个点数进行{trials_per_point}次试验)')
    plt.xscale('log')
    plt.yscale('log')
    plt.legend()
    plt.grid(True, which="both", ls="-", alpha=0.3)
    plt.show()
    
    return slope_fit

# 运行分析（这需要一些计算时间）
slope = analyze_pi_error(max_points=50000, step=2500, trials_per_point=30)
print(f"拟合得到的斜率（对数坐标下）为: {slope:.4f}")

如果一切正常，拟合得到的斜率应该非常接近-0.5，这有力地验证了误差与1/√N成正比的理论。理解这一点至关重要，它能帮助你合理规划计算资源：不要指望通过线性增加模拟次数来线性提升精度，那是指数级的代价。

5. 超越圆周率：蒙特卡罗的广阔应用天地

掌握了估算π的核心技能后，你已经拥有了蒙特卡罗方法的“元能力”。这种能力可以迁移到无数更复杂、更实用的场景中。蒙特卡罗不再是象牙塔里的数学游戏，而是工程师、金融分析师、数据科学家手中的利器。

应用一：计算复杂积分 很多函数的积分没有解析解，或者在高维情况下计算极其复杂（“维数灾难”）。蒙特卡罗积分却能优雅地处理。其思想与估算π如出一辙：在积分区域内随机采样，用函数值的平均值乘以区域体积来近似积分值。

例如，计算函数 ( f(x) = e^{-x^2} ) 在区间[0, 2]上的积分。解析解涉及误差函数，但用蒙特卡罗很简单：

def mc_integration(func, a, b, num_points):
    """使用蒙特卡罗方法计算一维定积分 [a, b]"""
    # 在[a, b]区间均匀采样
    x_samples = np.random.uniform(a, b, num_points)
    # 计算函数值
    f_values = func(x_samples)
    # 积分估计值 = 区间宽度 * 函数值的平均值
    integral_estimate = (b - a) * np.mean(f_values)
    return integral_estimate

# 定义被积函数
def my_function(x):
    return np.exp(-x**2)

# 计算积分
a, b = 0, 2
num_points = 100000
integral_est = mc_integration(my_function, a, b, num_points)

# 与通过SciPy数值积分得到的高精度结果进行比较（仅作参考）
from scipy import integrate
integral_exact, _ = integrate.quad(my_function, a, b)

print(f"蒙特卡罗估计积分值: {integral_est:.6f}")
print(f"SciPy高精度积分值: {integral_exact:.6f}")
print(f"绝对误差: {abs(integral_est - integral_exact):.6f}")

应用二：项目管理与风险分析 在项目管理中，任务时长、成本往往具有不确定性。蒙特卡罗模拟可以用来预测项目总工期或总成本的概率分布。假设一个项目有三个任务，每个任务的工期估计如下（单位：天）：

任务	乐观时间	最可能时间	悲观时间
A	5	7	10
B	8	10	15
C	3	4	6

我们可以为每个任务假设一个三角分布，然后通过大量模拟来预测项目总工期。

def triangular_random(low, mode, high, size):
    """生成三角分布的随机数（简易实现）"""
    u = np.random.uniform(0, 1, size)
    condition = u < (mode - low) / (high - low)
    return np.where(condition,
                    low + np.sqrt(u * (high - low) * (mode - low)),
                    high - np.sqrt((1 - u) * (high - low) * (high - mode)))

def project_duration_simulation(num_simulations):
    """模拟项目总工期"""
    # 定义每个任务的三角分布参数
    task_A = (5, 7, 10)
    task_B = (8, 10, 15)
    task_C = (3, 4, 6)
    
    # 生成随机工期（假设任务顺序进行）
    duration_A = triangular_random(*task_A, num_simulations)
    duration_B = triangular_random(*task_B, num_simulations)
    duration_C = triangular_random(*task_C, num_simulations)
    
    total_duration = duration_A + duration_B + duration_C
    return total_duration

# 运行模拟
num_sims = 50000
durations = project_duration_simulation(num_sims)

# 分析结果
mean_duration = np.mean(durations)
median_duration = np.median(durations)
p90 = np.percentile(durations, 90)  # 90%分位数
p95 = np.percentile(durations, 95)

print(f"项目总工期模拟结果 (基于{num_sims:,}次模拟):")
print(f"  平均工期: {mean_duration:.2f} 天")
print(f"  中位数工期: {median_duration:.2f} 天")
print(f"  90%的可能工期在 {p90:.2f} 天以内")
print(f"  95%的可能工期在 {p95:.2f} 天以内")

# 绘制工期分布直方图
plt.figure(figsize=(10, 6))
plt.hist(durations, bins=50, density=True, edgecolor='black', alpha=0.7)
plt.axvline(mean_duration, color='r', linestyle='--', linewidth=2, label=f'均值 = {mean_duration:.2f}')
plt.axvline(p90, color='orange', linestyle=':', linewidth=2, label=f'90%分位点 = {p90:.2f}')
plt.xlabel('项目总工期 (天)')
plt.ylabel('概率密度')
plt.title('蒙特卡罗模拟：项目总工期不确定性分析')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

通过这样的模拟，项目经理不再只是给出一个单一的总工期估计（比如21天），而是可以自信地说：“根据模拟，我们有90%的把握在23.5天内完成项目。” 这种基于概率的决策，远比一个脆弱的点估计更有价值。

应用三：金融期权定价 在量化金融中，蒙特卡罗方法是给复杂金融衍生品（如亚式期权、障碍期权）定价的核心工具。其基本思路是模拟标的资产（如股票）未来可能的价格路径，然后在这些路径上计算期权的平均收益，最后折现回现值。虽然完整的期权定价模型涉及随机微分方程（如几何布朗运动），但其蒙特卡罗模拟的骨架与我们之前的例子一脉相承。

从估算圆周率这个具体的项目出发，我们实际上解锁了一套解决广泛问题的思维框架：定义问题空间 -> 构建随机模型 -> 大量抽样 -> 统计聚合 -> 分析误差。无论你未来是进入数据科学、量化金融、工程仿真还是人工智能领域，这套用随机性驾驭复杂性的能力，都将是你工具箱里一件趁手的兵器。