概率论实战：用Python模拟经典概率问题（附代码示例）

概率论不仅是数学的重要分支，更是数据科学、人工智能和金融工程等领域的基石。但很多初学者常陷入纯理论学习的困境，难以将抽象概念与实际应用联系起来。本文将通过Python代码实现五个经典概率问题，从蒙特卡洛模拟到贝叶斯推理，带你感受"用代码思考概率"的独特魅力。

1. 蒙特卡洛方法：估算圆周率π

想象一下，如果给你无数个随机点，你能算出圆周率吗？18世纪法国数学家布丰提出的"投针问题"启发了蒙特卡洛方法。我们用Python再现这个思想实验：

import random
import matplotlib.pyplot as plt

def estimate_pi(num_points):
    points_inside = 0
    points_history = []
    
    for i in range(1, num_points+1):
        x, y = random.random(), random.random()
        if x**2 + y**2 <= 1:
            points_inside += 1
        if i % 1000 == 0:
            points_history.append(4 * points_inside / i)
    
    plt.plot(range(1000, num_points+1, 1000), points_history)
    plt.axhline(y=3.14159, color='r', linestyle='--')
    plt.xlabel('Number of points')
    plt.ylabel('Estimated π value')
    plt.show()
    return 4 * points_inside / num_points

print(estimate_pi(1000000))

运行这段代码，你会看到随着采样点增加，估计值逐渐收敛到真实π值。这揭示了蒙特卡洛方法的核心特点：

• 随机采样：通过大量随机实验逼近理论值
• 收敛速度：精度与√N成正比，N为样本量
• 可视化验证：动态观察估计过程增强理解

注意：实际应用中会根据问题特性采用重要性采样、马尔可夫链等方法提高效率

2. 贝叶斯定理：疾病检测中的反直觉

假设某种疾病在人群中的患病率为1%，检测准确率为99%。若某人检测为阳性，实际患病的概率是多少？这个经典问题常让人大跌眼镜——答案不是99%，而是约50%！

def bayes_theorem(p_disease, p_positive_given_disease, p_positive_given_no_disease):
    # 计算检测为阳性的总概率
    p_positive = (p_disease * p_positive_given_disease + 
                 (1 - p_disease) * p_positive_given_no_disease)
    # 应用贝叶斯公式
    p_disease_given_positive = (p_disease * p_positive_given_disease) / p_positive
    return p_disease_given_positive

result = bayes_theorem(p_disease=0.01,
                      p_positive_given_disease=0.99,
                      p_positive_given_no_disease=0.01)
print(f"检测阳性时实际患病概率: {result:.2%}")

这个案例展示了贝叶斯推理的关键要点：

概念	数值	说明
先验概率	1%	检测前的患病概率
似然度	99%	患者检测为阳性的概率
边际似然	1.98%	任意人检测为阳性的总概率
后验概率	50%	检测阳性后的患病概率

3. 赌徒谬误：独立事件的认知陷阱

"已经连续开了5次红色，下次肯定是黑色！"——这是赌场最常见的概率误解。我们用代码模拟轮盘赌的独立事件特性：

import numpy as np
import seaborn as sns

def roulette_simulation(num_spins):
    results = np.random.choice(['红','黑','绿'], 
                              size=num_spins,
                              p=[18/37, 18/37, 1/37])
    
    streaks = []
    current_streak = 1
    
    for i in range(1, num_spins):
        if results[i] == results[i-1] and results[i] != '绿':
            current_streak += 1
        else:
            if current_streak > 1:
                streaks.append(current_streak)
            current_streak = 1
    
    sns.histplot(streaks, bins=range(2,12))
    plt.xlabel('连续出现次数')
    plt.ylabel('发生频率')
    plt.show()

roulette_simulation(100000)

通过分析10万次旋转数据，我们会发现：

• 独立事件无记忆性：每次旋转都是独立的18/37概率
• 长连串并非不可能：在大量试验中，7-8次连续很常见
• 条件概率误区：P(红|连续5红) = P(红) = 18/37

4. 生日问题：惊人的概率直觉

"23个人中至少两人生日相同的概率有多大？"多数人猜测不到50%，实际概率却超过50%。这个反直觉现象用组合数学解释更清晰：

def birthday_paradox(max_people):
    probabilities = []
    for n in range(1, max_people+1):
        prob_unique = 1.0
        for i in range(1, n):
            prob_unique *= (365 - i) / 365
        probabilities.append(1 - prob_unique)
    
    plt.plot(range(1, max_people+1), probabilities)
    plt.axhline(y=0.5, color='r', linestyle='--')
    plt.xlabel('人数')
    plt.ylabel('至少两人同生日概率')
    plt.show()

birthday_paradox(60)

关键发现：

• 23人时概率达50.7%
• 57人时概率达99%
• 计算公式：1 - (365!/((365-n)! * 365^n))

5. 马尔可夫链：天气预测模型

最后我们构建一个简单的天气状态转移模型。假设天气只有"晴"、"雨"两种状态，转移矩阵如下：

import numpy as np

transition_matrix = np.array([
    [0.8, 0.2],  # 晴天→晴天80%，晴天→雨天20%
    [0.3, 0.7]   # 雨天→晴天30%，雨天→雨天70%
])

def forecast(initial_state, days):
    states = ['晴', '雨']
    current_state = initial_state
    sequence = [states[current_state]]
    
    for _ in range(days):
        current_state = np.random.choice(
            [0,1], 
            p=transition_matrix[current_state]
        )
        sequence.append(states[current_state])
    
    return sequence

print(forecast(0, 30))  # 从晴天开始预测30天

这个简单模型揭示了：

• 稳态分布：通过矩阵幂运算可求得长期天气分布
• 无后效性：下一状态仅取决于当前状态
• 应用场景：从PageRank算法到金融风险建模

进阶思考：概率编程实践

掌握了这些基础模型后，可以尝试用PyMC3等概率编程库解决更复杂问题。比如建立贝叶斯线性回归模型：

import pymc3 as pm
import numpy as np

# 生成模拟数据
np.random.seed(42)
X = np.linspace(0, 10, 100)
true_slope = 2.5
y = true_slope * X + np.random.normal(0, 2, size=100)

with pm.Model() as model:
    # 定义先验分布
    slope = pm.Normal('slope', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=1)
    
    # 定义似然函数
    likelihood = pm.Normal('y', mu=slope*X, 
                          sigma=sigma, observed=y)
    
    # 采样
    trace = pm.sample(2000, return_inferencedata=False)

print(pm.summary(trace))

这种建模方式允许我们：

1. 明确表达关于参数的先验信念
2. 基于观测数据更新认知
3. 量化估计的不确定性
4. 灵活扩展复杂模型结构

从这些案例可以看出，编程不仅是验证概率理论的工具，更是培养概率思维的催化剂。当你下次遇到概率难题时，不妨先写几行代码做个实验——这往往比纯数学推导更能揭示问题的本质。