程序设计与计算思维¶

Computer Programming and Computational Thinking

第 11 讲：随机模拟的建模¶

2025—2026学年度春季学期

清华大学 韩文弢

问题引入¶

很多情况下容易知道单个个体的行为方式（微观视角，microscopic），例如：

  • 抛一个硬币会产生什么结果
  • 一个灯泡在某一天坏掉的可能性是多少
  • 疾病在人和人之间如何传染

从研究的角度，更想知道整个群体会呈现什么情况（宏观视角，macroscopic），例如：

  • 抛很多次硬币之后会有什么样的分布
  • 灯泡在哪天坏掉是什么样的分布
  • 人群中被感染的群体随着时间如何变化

如何用计算思维研究此类问题？

可以使用模拟（simulation）方法。

对成功时间（或失败时间）的建模¶

考虑一个简单的成功时间模型：假设玩一个游戏时每轮成功的概率为 $p$，问需要多少轮才能成功？

例如，掷骰子需要掷多少次才能得到6？掷两个骰子需要多少次才能得到双6？

此类模型可以用于多种场景：

  • 灯泡或机器的失效时间
  • 放射性原子核衰变时间
  • 感染恢复时间
  • ……

其基本问题是：

假设第0天有 $N$ 个正常工作的灯泡。

  • 如果每个灯泡每天出现故障的概率为 $p$，第 $t$ 天还有多少个灯泡在工作（没出故障）？
  • 平均而言，一个灯泡能使用多久？
  • 真实情况下灯泡会精确地在每天午夜失效吗？能否建立更真实的模型？

准备工作¶

In [1]:
from collections import Counter

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import patches
from ipywidgets import interact, FloatSlider, IntSlider

plt.rcParams['font.sans-serif'] = ['Noto Sans CJK SC']

可视化设备故障¶

模拟和可视化设备故障的过程。

In [2]:
def simulate_grid(M, prob):
    """模拟 M*M 个设备的故障过程，返回每个设备的故障时间。"""
    v = np.zeros((M, M), dtype=int)
    t = 0
    # 注意：下面的操作都在矩阵上进行
    while np.any(v == 0) and t < 100:
        t += 1
        # 创建一个掩码：值为 True 的位置是要被选中的位置
        mask = (np.random.rand(M, M) < prob) & (v == 0)
        v[mask] = t
    return v
In [3]:
def plot_grid(M, prob, tt):
    simulation = simulate_grid(M, prob)

    fig, axes = plt.subplots(1, 3, figsize=(16, 5))

    ax = axes[0]
    for i in range(M):
        for j in range(M):
            rect = patches.FancyBboxPatch((i, j), 0.9, 0.9,
                                          boxstyle="round,pad=0.05",
                                          facecolor="gray", alpha=0.5)
            ax.add_patch(rect)
            if simulation[i, j] >= tt:
                color = "lime"
            elif simulation[i, j] == tt - 1:
                color = "red"
            else:
                color = "purple"
            circle = plt.Circle((i + 0.45, j + 0.45), 0.3, color=color)
            ax.add_patch(circle)
            if simulation[i, j] < tt:
                ax.text(i + 0.45, j + 0.45, str(simulation[i, j]),
                       ha="center", va="center", color="white", fontsize=7)
    ax.set_xlim(-0.1, M + 0.1)
    ax.set_ylim(-0.1, M + 0.1)
    ax.set_aspect("equal")
    ax.set_title(f"设备状态：时间 {tt - 1}，故障个数 {np.sum(simulation < tt)}")
    ax.axis("off")

    ax2 = axes[1]
    cdf = [np.sum(simulation < i) for i in range(tt)]
    surviving = [np.sum(simulation > i) for i in range(tt)]
    x_range = range(tt)
    ax2.bar(x_range, cdf, color="purple", alpha=0.8)
    ax2.bar(x_range, surviving, color="lime", alpha=0.8)
    ax2.set_xlabel("时间")
    ax2.set_ylabel("个数")
    ax2.set_title("设备统计信息")

    ax3 = axes[2]
    flat = simulation.flatten()
    counts = Counter(flat)
    times = range(tt)
    vals = [counts.get(t, 0) for t in times]
    ax3.bar(times, vals, color="red", alpha=0.8)
    ax3.set_xlabel("时间")
    ax3.set_ylabel("个数")
    ax3.set_title("设备故障信息")

    plt.tight_layout()
    plt.show()

interact(plot_grid,
         M=IntSlider(min=2, max=20, value=8, description="M", continuous_update=False),
         prob=FloatSlider(min=0.01, max=1.0, step=0.01, value=0.1, description="prob", continuous_update=False),
         tt=IntSlider(min=1, max=100, value=1, description="t", continuous_update=False));

interactive(children=(IntSlider(value=8, continuous_update=False, description='M', max=20, min=2), FloatSlider…

数学基础：伯努利随机变量¶

  • 概念：模拟不均匀的硬币。
    • 以概率 $p$ 取值 1（发生/成功）。
    • 以概率 $(1 - p)$ 取值 0（不发生/失败）。
In [4]:
def bernoulli(p):
    return int(np.random.rand() < p)
  • np.random.rand() < p 返回布尔值。
  • 利用 int() 转换为 1 或 0 完成单次采样。
In [5]:
flips = [bernoulli(0.25) for _ in range(100)]
print(f"前 20 次结果: {flips[:20]}")

前 20 次结果: [0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0]

期望与方差¶

In [6]:
print(f"均值 = {np.mean(flips):.4f}  (理论值 = 0.25)")

均值 = 0.2600  (理论值 = 0.25)
  • 均值（即 1 出现的比例）会趋近于概率 $p$。

思考练习：如何计算伯努利随机变量的方差？（提示：要做中心化）

代码设计：面向对象（OOP）封装¶

  • 将数学概念代码化，创建一个类来专门表示和处理“伯努利随机变量”这一实体。
In [7]:
class Bernoulli:
    def __init__(self, p):
        self.p = p

    def rand(self):
        return int(np.random.rand() < self.p)

    def mean(self):
        return self.p

    def __repr__(self):
        return f"Bernoulli(p={self.p})"
In [8]:
B = Bernoulli(0.25)
print(B)
print(f"一次采样: {B.rand()}")
print(f"理论均值: {B.mean()}")

Bernoulli(p=0.25)
一次采样: 0
理论均值: 0.25

注：在实际使用中，可直接调用现成库：scipy.stats.bernoulli。

进行随机模拟¶

In [9]:
def step(infectious, p):
    for i in range(len(infectious)):
        if infectious[i] and bernoulli(p):
            infectious[i] = False
    return infectious

def simulate_recovery(N, p, T):
    infectious = [True] * N
    num_infectious = [N]
    for t in range(1, T + 1):
        step(infectious, p)
        num_infectious.append(sum(infectious))
    return num_infectious
In [10]:
def plot_simulation(N, p, T):
    fig, ax = plt.subplots(figsize=(8, 4))
    for run in range(1, 3):
        data = simulate_recovery(N, p, T)
        ax.plot(data, "o-", alpha=0.5, lw=2, markersize=3, label=f"第 {run} 次模拟")
    ax.set_xlabel("时间")
    ax.set_ylabel("仍然正常工作的数量")
    ax.legend()
    plt.show()

interact(plot_simulation,
         N=IntSlider(min=1, max=1000, value=70, description="N", continuous_update=False),
         p=FloatSlider(min=0.0, max=1.0, step=0.01, value=0.25, description="p", continuous_update=False),
         T=IntSlider(min=10, max=200, value=100, description="T", continuous_update=False));

interactive(children=(IntSlider(value=70, continuous_update=False, description='N', max=1000, min=1), FloatSli…

多次运行与平均化处理¶

In [11]:
N_bulbs = 70
pp = 0.05
T = 100

all_data = [simulate_recovery(N_bulbs, pp, T) for _ in range(30)]

fig, ax = plt.subplots(figsize=(8, 5))
for d in all_data:
    ax.plot(d, "o-", alpha=0.1, markersize=1, color="steelblue")

mean_data = np.mean(all_data, axis=0)
ax.plot(mean_data, "o-", lw=3, color="red", alpha=0.7, markersize=3, label="均值")

ax.set_xlabel("时间")
ax.set_ylabel("仍然正常工作的数量")
ax.legend()
plt.title("多次随机模拟和均值")
plt.show()
No description has been provided for this image

对数尺度下的视角¶

In [12]:
fig, ax = plt.subplots(figsize=(8, 5))
for d in all_data:
    d_plot = [x if x > 0 else np.nan for x in d]
    ax.plot(d_plot, "o-", alpha=0.1, markersize=1, color="steelblue")

mean_plot = [x if x > 0 else np.nan for x in mean_data]
ax.plot(mean_plot, "o-", lw=3, color="red", alpha=0.7, markersize=3, label="均值")

ax.set_yscale("log")
ax.set_xlabel("时间")
ax.set_ylabel("仍然正常工作的数量")
ax.legend()
plt.title("对数尺度的模拟结果")
plt.show()
No description has been provided for this image

均值的时间演化：理论推导¶

设 $I_t$ 为时刻 $t$ 时正常工作的灯泡数量。

  • 故障平均数：$p \cdot I_t$
  • 差分递推关系： $$\Delta I_t = I_{t+1} - I_t = -p \, I_t$$
  • 推导解析解：$I_{t+1} = (1 - p) I_t$ $$I_t = (1-p)^t \, I_0$$

对比精确解与数值结果¶

In [13]:
exact = [N_bulbs * (1 - pp)**t for t in range(T + 1)]

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(mean_data, "o", alpha=0.5, markersize=4, label="随机模拟")
ax.plot(exact, lw=3, alpha=0.8, label="理论计算")
ax.set_xlabel("时间")
ax.set_ylabel("仍然正常工作的数量")
ax.set_title("实验与理论的对比")
ax.legend(loc="upper right")
plt.show()
No description has been provided for this image
  • 结论：理论模型与随机均值高度吻合。
  • 大数定理：群体规模越大，系统层面的随机涨落（波动）越小。

进阶分布：二项分布¶

  • 初始时存在 $N_0$ 个元件，第一步故障总数 $\Delta N_0$，是随机变量。
  • 它是由 $N_0$ 个独立的伯努利变量相加构成： $$\Delta N_0 = \sum_{i=1}^{N_0} B_0^{(i)}$$
  • 这一分布被称为二项分布 (binomial distribution)。
In [14]:
class Binomial:
    def __init__(self, N, p):
        self.N = N
        self.p = p

    def rand(self):
        return sum(Bernoulli(self.p).rand() for _ in range(self.N))

    def mean(self):
        return self.N * self.p

    def __repr__(self):
        return f"Binomial(N={self.N}, p={self.p})"
In [15]:
print(f"Binomial(10, 0.25) sample: {Binomial(10, 0.25).rand()}")

Binomial(10, 0.25) sample: 1
In [16]:
def plot_binomial(binomial_N, binomial_p):
    data = [Binomial(binomial_N, binomial_p).rand() for _ in range(10000)]
    counts = Counter(data)
    xs = sorted(counts.keys())
    ys = [counts[x] for x in xs]

    fig, ax = plt.subplots(figsize=(8, 4))
    ax.bar(xs, ys, alpha=0.5, color="steelblue")
    ax.set_xlabel("值")
    ax.set_ylabel("次数")
    ax.set_title(f"Binomial(N={binomial_N}, p={binomial_p}) - 10000 次取样")
    plt.show()

interact(plot_binomial,
         binomial_N=IntSlider(min=1, max=100, value=1, description="N", continuous_update=False),
         binomial_p=FloatSlider(min=0.0, max=1.0, step=0.01, value=0.0, description="p", continuous_update=False));

interactive(children=(IntSlider(value=1, continuous_update=False, description='N', min=1), FloatSlider(value=0…

二项系数¶

  • 设保留概率 $q := 1 - p$。
  • 从 $n$ 个对象中恰好有 $k$ 个失效的组合方式，由二项系数定义： $$\binom{n}{k} := \frac{n!}{k! \, (n-k)!}$$

连续时间演化¶

  • 时间步缩小：设时间增量为 $\delta t$。当 $\delta t \to 0$ 时，引入单位时间的失效速率 $\lambda$： $$p(\delta t) \simeq \lambda \, \delta t$$
  • 微分方程推导： $$I(t + \delta t) - I(t) \simeq -\lambda \,\delta t \, I(t)$$ $$\frac{dI(t)}{dt} = -\lambda \, I(t)$$
  • 精确解与复利：
    • 指数衰减：$I(t) = I_0 \exp(-\lambda \, t)$
    • 相当于复利公式：$I_{t} = (1 - \lambda \, \delta t)^{(t / \delta t)} I_0$

从离散到连续¶

In [17]:
def plot_cumulative(ax, p, N, delta=1, **kwargs):
    ps = [p * (1 - p)**(n - 1) for n in range(1, N + 1)]
    cumulative = np.cumsum(ps)
    ax.scatter([n * delta for n in range(1, N + 1)], cumulative, **kwargs)

fig, ax = plt.subplots(figsize=(8, 4))
plot_cumulative(ax, 0.1, 30, 1.0, s=10, c="red", alpha=1, label="dt=1.0")
plot_cumulative(ax, 0.05, 60, 0.5, s=10, c="limegreen", alpha=1, label="dt=0.5")
plot_cumulative(ax, 0.025, 120, 0.25, s=5, c="limegreen", alpha=0.7, label="dt=0.25")
plot_cumulative(ax, 0.0125, 240, 0.125, s=3, c="limegreen", alpha=0.5, label="dt=0.125")

lam = -np.log(1 - 0.1)
t_cont = np.linspace(0, 20, 500)
ax.plot(t_cont, 1 - np.exp(-lam * t_cont), lw=2, color="blue", label=f"exact: $\\lambda$={lam:.4f}")
ax.plot(t_cont, 1 - np.exp(-0.1 * t_cont), lw=2, color="orange", label="approx: $\\lambda$=0.1")

ax.set_xlabel("时间")
ax.set_ylabel("积累概率")
ax.set_title("从离散到连续")
ax.legend(fontsize=8)
plt.show()
No description has been provided for this image
  • 速率与概率：速率是单位时间框架内的概率测度。
  • 现实意义：宏观群体的衰减/故障发生在 $\delta t \to 0$ 的时间连续变化中。

SIR 模型¶

现在将方法扩展到更加完整的 SIR 模型，$S \to I \to R$，常用来对传染病建模。其中，

  • S：susceptible，易感者
  • I：infectious，感染者
  • R：recovered，康复者

设 $S_t$ 和 $I_t$ 分别为时刻 $t$ 的易感者和感染者数量，$N$ 为总人口。

平均而言，每一轮每个感染者有机会与另一个体接触。该个体从总人口 $N$ 中均匀随机选择。但只有当选中的个体是易感者（概率 $S_t/N$）且感染成功（概率 $b$）时，才会发生新感染。

因此感染者数量的变化为

$$\Delta I_t = I_{t+1} - I_t = b \, I_t \, \left(\frac{S_t}{N} \right)$$

$S_t$ 的减少也是 $\Delta I_t$。

简便起见，用 $N$ 来归一化，定义

$$s_t := \frac{S_t}{N}; \quad i_t := \frac{I_t}{N}; \quad r_t := \frac{R_t}{N}$$

加入恢复概率 $c$，得到离散时间 SIR 模型：

$$\begin{align} s_{t+1} &= s_t - b \, s_t \, i_t \\ i_{t+1} &= i_t + b \, s_t \, i_t - c \, i_t\\ r_{t+1} &= r_t + c \, i_t \end{align}$$

同样可以允许过程以步长 $\delta t$ 进行，取极限 $\delta t \to 0$。用速率 $\beta$ 和 $\gamma$ 得到标准的（连续时间）SIR 模型：

$$\begin{align} \frac{ds(t)}{dt} &= -\beta \, s(t) \, i(t) \\ \frac{di(t)}{dt} &= +\beta \, s(t) \, i(t) - \gamma \, i(t)\\ \frac{dr(t)}{dt} &= +\gamma \, i(t) \end{align}$$

由于该常微分方程（ordinary differential equation，ODE）没有解析解，可以使用离散时间模型的模拟来求数值解。

注意：求解 ODE 系统的最简单数值方法欧拉法，基本上就是求解离散时间模型。

更高级的 ODE 求解器工具可在 Python 的 scipy.integrate 中找到。

In [18]:
# 固定参数
NN = 100
SS = NN - 1
II = 1
RR = 0
In [19]:
ss, ii, rr = SS/NN, II/NN, RR/NN
In [20]:
p_infection, p_recovery = 0.1, 0.01
In [21]:
TT = 1000
In [22]:
def discrete_SIR(s0, i0, r0, T=1000, p_infection=0.1, p_recovery=0.01):
    """离散时间 SIR 模型"""
    s, i, r = s0, i0, r0
    results = [(s, i, r)]
    
    for t in range(1, T+1):
        delta_i = p_infection * s * i
        delta_r = p_recovery * i
        
        s_new = s - delta_i
        i_new = i + delta_i - delta_r
        r_new = r + delta_r
        
        results.append((s_new, i_new, r_new))
        s, i, r = s_new, i_new, r_new
    
    return results
In [23]:
SIR = discrete_SIR(ss, ii, rr)
In [24]:
ts = range(len(SIR))
plt.figure(figsize=(10, 6))
plt.plot(ts, [x[0] for x in SIR], 'o-', alpha=0.2, label='S', markersize=2)
plt.plot(ts, [x[1] for x in SIR], 'o-', alpha=0.2, label='I', markersize=2)
plt.plot(ts, [x[2] for x in SIR], 'o-', alpha=0.2, label='R', markersize=2)
plt.xlabel('时间')
plt.ylabel('人口比例')
plt.legend()
plt.xlim(0, 500)
plt.grid(True, alpha=0.3)
plt.show()
No description has been provided for this image
In [25]:
# 交互式 SIR 模型
@interact
def interactive_SIR(N=IntSlider(min=50, max=500, step=10, value=100, description='总人口'),
                  I0=IntSlider(min=1, max=10, step=1, value=1, description='初始感染者'),
                  p_infection=FloatSlider(min=0.01, max=0.5, step=0.01, value=0.1, description='感染概率'),
                  p_recovery=FloatSlider(min=0.001, max=0.1, step=0.001, value=0.01, description='恢复概率')):
    S0 = N - I0
    R0 = 0
    ss, ii, rr = S0/N, I0/N, R0/N
    
    SIR = discrete_SIR(ss, ii, rr, 500, p_infection, p_recovery)
    
    ts = range(len(SIR))
    plt.figure(figsize=(12, 5))
    
    plt.subplot(1, 2, 1)
    plt.plot(ts, [x[0] for x in SIR], '-', alpha=0.7, lw=2, label='易感者(S)')
    plt.plot(ts, [x[1] for x in SIR], '-', alpha=0.7, lw=2, label='感染者(I)')
    plt.plot(ts, [x[2] for x in SIR], '-', alpha=0.7, lw=2, label='恢复者(R)')
    plt.xlabel('时间')
    plt.ylabel('人口比例')
    R0_val = p_infection / p_recovery
    plt.title(f'SIR 模型 ($R_0$ = {R0_val:.2f})')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    plt.subplot(1, 2, 2)
    plt.stackplot(ts, [x[0] for x in SIR], [x[1] for x in SIR], [x[2] for x in SIR],
                  labels=['易感者', '感染者', '恢复者'], alpha=0.7, colors=['blue', 'red', 'green'])
    plt.xlabel('时间')
    plt.ylabel('人口比例')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()

interactive(children=(IntSlider(value=100, description='总人口', max=500, min=50, step=10), IntSlider(value=1, de…

本讲小结¶

随机模拟的建模过程：

  • 在个体模型中，明确指定每个个体的行为规则
  • 从个体模型出发，推导出确定性的宏观方程
    • 可以是离散时间的递推关系
    • 也可以取连续极限转换为常微分方程
  • 通过计算模拟结果，寻找规律