程序设计与计算思维¶

Computer Programming and Computational Thinking

第 13 讲：随机游走¶

2025—2026学年度春季学期

清华大学 韩文弢

问题引入¶

  • 布朗运动（Brownian motion，1827年）：观察到水中大颗粒（如花粉微粒）呈现随机游动。
  • 爱因斯坦在1905年解释布朗运动：水分子反复碰撞颗粒的宏观表现。
  • 微观确定与宏观随机：
    • 模拟相互碰撞的圆盘
    • 每个圆盘严格遵循牛顿定律
    • 但仅观察单一圆盘时，其轨迹呈现高度随机性

准备工作¶

In [1]:
import timeit  # 用于计时

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

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

可视化随机游走¶

随机游走 (Random Walk) 模型

  • 定义：模拟物体在时空中的随机运动过程。
  • 规则：在每个离散时间步，物体朝随机方向移动固定距离。
  • 目标：在二维空间中模拟并可视化其运动轨迹。
In [2]:
def walker2d_step():
    directions = [[1, 0], [0, 1], [-1, 0], [0, -1]]
    return directions[np.random.randint(4)]

def trajectory_2d(N):
    x, y = 0, 0
    positions = [(x, y)]
    for _ in range(N):
        dx, dy = walker2d_step()
        x += dx
        y += dy
        positions.append((x, y))
    return positions

_traj2d = trajectory_2d(10)
_N_val = 1

import math

@interact(N=IntSlider(min=1, max=6, step=1, value=1, description="N:"),
          t=IntSlider(min=1, max=100, step=1, value=1, description="t (%):"))
def plot_walk2d(N, t):
    global _traj2d, _N_val
    if N != _N_val:
        _N_val = N
        _traj2d = trajectory_2d(10**N)
    traj = _traj2d
    t = int(math.ceil(len(traj) * t / 100))

    xs = [p[0] for p in traj[:t]]
    ys = [p[1] for p in traj[:t]]

    fig, ax = plt.subplots(figsize=(8, 8))
    ax.plot(xs, ys, alpha=0.5, linewidth=2)
    ax.scatter([xs[0]], [ys[0]], c="red", s=100, zorder=5, label="起点")
    ax.scatter([xs[-1]], [ys[-1]], c="green", s=100, zorder=5, label="终点")

    all_xs = [p[0] for p in traj]
    all_ys = [p[1] for p in traj]
    ax.set_xlim(min(all_xs) - 1, max(all_xs) + 1)
    ax.set_ylim(min(all_ys) - 1, max(all_ys) + 1)
    ax.set_aspect("equal")
    ax.legend()
    plt.show()

interactive(children=(IntSlider(value=1, description='N:', max=6, min=1), IntSlider(value=1, description='t (%…

观察结论：这种宏观的数学动力学，在定性上极其相似于物理中硬盘碰撞的杂乱运动。

为什么要使用随机游走建模？¶

随机过程的优势：

  1. 规避数据爆炸：无需获取海量微观数据来精确描述复杂系统。
  2. 本质洞察：通过简化的概率假设，更容易提取系统的宏观物理法则。

典型应用场景：

  • 金融：股票价格的涨跌波动
  • 环境：污染物在空气中的扩散
  • 生物：中性基因在种群中的传播

简单随机游走¶

  • 最简模型：在一维整数坐标上，每步向左或向右跳跃。
  • 数学抽象：
    • 每一步是独立的随机变量。
    • 取值为 $\pm 1$，概率均为 $1/2$。
  • 本质：抛硬币实验（伯努利随机变量）的线性缩放映射。

Python：基准测试¶

寻找最快的 $\pm 1$ 生成方式：

  • 模拟需要运行数百万次的核心循环。
  • 微小的性能差异将被无限放大。
  • 工具：Python 内置 timeit 模块。
  • 原理：重复运行代码片段并统计耗时分布。
In [3]:
def step1():
    return np.random.choice([-1, +1])

def step2():
    return 2 * int(np.random.random() < 0.5) - 1

def step3():
    return 2 * np.random.randint(0, 2) - 1

def step4():
    return int(np.sign(np.random.randn()))
In [4]:
n = 1000000

print("step1() (np.random.choice):")
result1 = timeit.timeit("step1()", globals=globals(), number=n)
print(f"  平均时间: {result1/n*1e6:.3f} 微秒")

print("\nstep2() (np.random.random < 0.5):")
result2 = timeit.timeit("step2()", globals=globals(), number=n)
print(f"  平均时间: {result2/n*1e6:.3f} 微秒")

print("\nstep3() (np.random.randint):")
result3 = timeit.timeit("step3()", globals=globals(), number=n)
print(f"  平均时间: {result3/n*1e6:.3f} 微秒")

print("\nstep4() (np.sign(np.random.randn)):")
result4 = timeit.timeit("step4()", globals=globals(), number=n)
print(f"  平均时间: {result4/n*1e6:.3f} 微秒")

step1() (np.random.choice):
  平均时间: 3.193 微秒

step2() (np.random.random < 0.5):
  平均时间: 0.188 微秒

step3() (np.random.randint):
  平均时间: 1.004 微秒

step4() (np.sign(np.random.randn)):
  平均时间: 0.322 微秒

计算随机游走轨迹¶

  • 轨迹记录：从初始位置 $x=0$ 开始，记录连续迈步后的累积位置序列。
In [5]:
def walk1D(N):
    x = 0
    xs = [x]

    for i in range(N):
        x += step2()
        xs.append(x)

    return xs
In [6]:
fig, ax = plt.subplots(figsize=(10, 6))

for i in range(10):
    ax.plot(walk1D(100), alpha=0.5, linewidth=2)

plt.show()
No description has been provided for this image

使用类实现的随机游走¶

需求拓展：支持 1D、2D 甚至多维的各种游走规则。

最佳实践：使用面向对象编程中的类与多态。

  • 建立抽象基类 Walker。
  • 统一行为接口：获取位置、生成步长、更新状态。
In [7]:
from abc import ABC, abstractmethod

class Walker(ABC):
    @abstractmethod
    def position(self):
        pass

    @abstractmethod
    def step(self):
        pass

    @abstractmethod
    def update(self, step_val):
        pass
In [8]:
class Walker1D(Walker):
    def __init__(self, pos=0):
        self.pos = pos

    def position(self):
        return self.pos

    def step(self):
        return step2()

    def update(self, step_val):
        return Walker1D(self.pos + step_val)
In [9]:
class Walker2D(Walker):
    def __init__(self, x=0, y=0):
        self.x = x
        self.y = y

    def position(self):
        return (self.x, self.y)

    def step(self):
        directions = [[1, 0], [0, 1], [-1, 0], [0, -1]]
        return directions[np.random.randint(4)]

    def update(self, step_val):
        return Walker2D(self.x + step_val[0], self.y + step_val[1])
In [10]:
def trajectory(walker, N):
    positions = [walker.position()]

    for i in range(N):
        step_val = walker.step()
        walker = walker.update(step_val)
        positions.append(walker.position())

    return positions
In [11]:
trajectory(Walker2D(0, 0), 10)
Out[11]:
[(0, 0),
 (-1, 0),
 (-2, 0),
 (-2, 1),
 (-2, 2),
 (-2, 3),
 (-2, 2),
 (-1, 2),
 (-1, 1),
 (0, 1),
 (0, 2)]

随机游走的数学原理：随机变量之和¶

  • 在时间 $n$ 的位置 $S_n$，是 $n$ 个随机步长之和： $$S_n = X_1 + \cdots + X_n$$
  • $X_i$：第 $i$ 步的独立同分布 (independent and identically distributed, i.i.d., iid, 或 IID) 随机变量。
  • 核心难点：虽然每步 $X_i$ 独立，但中间状态 $S_n$ 与 $S_{n-1}$ 有状态依赖。

快速计算：累积和¶

首先，生成一组独立的随机步长向量：

In [12]:
steps = np.random.choice([-1, +1], size=10)
steps
Out[12]:
array([-1, -1, -1, -1, -1, -1, -1,  1,  1, -1])

轨迹等价于部分和 (Partial Sums)： $$\begin{align} S_1 & = X_1 \\ S_2 &= X_1 + X_2 \\ S_3 &= X_1 + X_2 + X_3 \end{align}$$

  • NumPy 的 cumsum (前缀和/扫描) 函数可直接、高效地计算整条轨迹：
In [13]:
np.cumsum(steps)
Out[13]:
array([-1, -2, -3, -4, -5, -6, -7, -6, -5, -6])

可视化验证：

In [14]:
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(np.cumsum(steps), "o-", linewidth=2, markersize=8)
plt.show()
No description has been provided for this image

从轨迹到概率分布的演化¶

视角切换：

  • 运行成千上万次模拟
  • 观察空间中的概率分布

主方程 (Master Equation)：

  • 设 $p_i^{(t)}$ 为时间 $t$ 位于位置 $i$ 的概率。
  • $t+1$ 时刻位于 $i$ 的概率，必然来自上一时刻其相邻位置的流入： $$p_i^{(t+1)} = \frac{1}{2} (p_{i-1}^{(t)} + p_{i+1}^{(t)})$$

设整条一维空间的概率分布为向量 $\mathbf{p}^{(t)}$，定义演化函数：

In [15]:
def evolve(p):
    p_new = np.zeros_like(p)

    for i in range(1, len(p) - 1):
        p_new[i] = 0.5 * (p[i-1] + p[i+1])

    p_new[0] = 0
    p_new[-1] = 0

    return p_new

数学观察：卷积

  • 这个邻域平均操作，本质上是一维离散卷积。
  • 类似图像处理中的“模糊/平滑滤波”。

两项基本设定：

  1. 边界条件 (Boundary Conditions)

    • 强制设两端为 0。
    • 物理意义：吸收边界（粒子游走到边界即消失/逃逸）。

  2. 初始条件 (Initial Condition)

    • $t=0$ 时的向量 $\mathbf{p}^{(0)}$。
    • 设定：所有粒子从系统中心点出发（中心概率为 1，其余为 0）。
In [16]:
def initial_condition(n):
    p0 = np.zeros(n)
    p0[n // 2 + 1] = 1
    return p0

概率演化的可视化呈现¶

In [17]:
def time_evolution(p0, N):
    ps = [p0.copy()]
    p = p0.copy()

    for i in range(N):
        p = evolve(p)
        ps.append(p.copy())

    return ps
In [18]:
p0 = initial_condition(101)
ps = time_evolution(p0, 100)
In [19]:
ps[2]
Out[19]:
array([0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ,
       0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ,
       0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ,
       0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ,
       0.  , 0.  , 0.  , 0.  , 0.  , 0.25, 0.  , 0.5 , 0.  , 0.25, 0.  ,
       0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ,
       0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ,
       0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ,
       0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ,
       0.  , 0.  ])
In [20]:
@interact(t=IntSlider(min=0, max=100, step=1, value=0, description="时间步 t:"))
def plot_distribution(t):
    fig, ax = plt.subplots(figsize=(10, 5))
    ax.plot(ps[t], linewidth=2)
    ax.set_ylim(0, 1)
    plt.show()

interactive(children=(IntSlider(value=0, description='时间步 t:'), Output()), _dom_classes=('widget-interact',))
In [21]:
M = np.array(ps)
fig, ax = plt.subplots(figsize=(10, 8))
im = ax.imshow(M, aspect="auto", origin="upper", cmap="hot")
plt.colorbar(im)
plt.show()
No description has been provided for this image

本讲小结¶

  • 随机游走模型：从布朗运动的物理背景出发，建立了随机游走的数学模型——每步朝随机方向移动固定距离，宏观呈现随机性。
  • 建模动机：随机建模可规避数据爆炸，以简化的概率假设提取系统的宏观法则，广泛应用于金融、环境、生物等领域。
  • 面向对象设计：利用抽象基类 Walker 与多态，统一了 1D/2D 游走的接口（position、step、update），使轨迹计算函数 trajectory 适用于任意维度的游走者。
  • 概率分布演化：从单轨迹转向概率分布视角，主方程 $p_i^{t+1} = \frac{1}{2}(p_{i-1}^t + p_{i+1}^t)$ 描述了概率随时间的演化，其本质是一维离散卷积；通过吸收边界条件与中心初始条件，可视化了概率分布从尖峰逐步扩散的过程。