程序设计与计算思维¶

Computer Programming and Computational Thinking

第 10 讲：抽样与随机变量¶

2025—2026学年度春季学期

清华大学 韩文弢

准备工作¶

In [ ]:
import string
from collections import defaultdict, Counter

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from scipy import stats
from ipywidgets import interact, FloatSlider, IntSlider

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

伪随机数生成的基本原理¶

计算机中的"随机数"实际上是通过确定性算法生成的，因此称为伪随机数。

线性同余生成器¶

最简单的伪随机数生成器之一：线性同余生成器 (linear congruential generator, LCG)

公式：$X_{n+1} = (a \cdot X_n + c) \mod m$

其中：

  • $X_n$：当前状态（种子）
  • $a$：乘数
  • $c$：增量
  • $m$：模数
In [ ]:
def lcg(seed, a, c, m):
    """线性同余生成器：返回下一个随机数"""
    return (a * seed + c) % m
In [ ]:
# 测试 LCG
seed = 42
a, c, m = 1664525, 1013904223, 2**32

for _ in range(5):
    seed = lcg(seed, a, c, m)
    print(seed)
In [ ]:
# 产生 1-6 之间的整数
seed = 42
a, c, m = 1664525, 1013904223, 2**32

counter = defaultdict(int)
for _ in range(10000):
    seed = lcg(seed, a, c, m)
    x = seed % 6 + 1
    counter[x] += 1
for x, n in sorted(counter.items()):
    print(f'{x}: {n}')

方法一：使用全局变量保存状态¶

在调用随机函数时，简便起见一般不希望带着状态参数。可以使用全局变量来实现这个目的。

In [ ]:
_seed = 42
_a, _c, _m = 1664525, 1013904223, 2**32

def next_random():
    """使用全局变量保存状态"""
    global _seed  # 要修改全局变量，必须用 global 声明为全局
    _seed = (_a * _seed + _c) % _m  # 使用全局变量可以不用 global 声明
    return _seed
In [ ]:
for _ in range(5):
    print(next_random())

全局变量方法的问题：

  1. 难以管理多个随机数生成器：所有函数共享同一个状态
  2. 代码耦合度高：函数依赖外部变量，不易复用
  3. 容易出错：全局状态容易被其他部分的代码误改
  4. 并发不安全：多线程环境下会产生竞争条件
In [ ]:
# 问题示例：无法同时使用两个独立的随机数生成器
print("生成器1:")
for _ in range(3):
    print(next_random(), end=" ")

print("\n生成器2:")  # 无法创建独立的第二个生成器！
for _ in range(3):
    print(next_random(), end=" ")

方法二：面向对象方法¶

使用面向对象方法可以将状态封装在对象中，避免暴露在全局空间，从而解决上述问题。

In [ ]:
class RandomGenerator:
    """线性同余随机数生成器类"""
    
    def __init__(self, seed=42, a=1664525, c=1013904223, m=2**32):
        """初始化一个随机数生成器"""
        self._seed = seed
        self._a = a
        self._c = c
        self._m = m
    
    def next(self):
        """生成下一个随机数"""
        self._seed = (self._a * self._seed + self._c) % self._m
        return self._seed
    
    def random(self):
        """生成 [0, 1) 区间的随机浮点数"""
        return self.next() / self._m
In [ ]:
# 创建多个独立的随机数生成器
gen1 = RandomGenerator(seed=100)
gen2 = RandomGenerator(seed=200)

print("生成器1:")
for _ in range(5):
    print(gen1.next(), end=" ")

print("\n生成器2:")
for _ in range(5):
    print(gen2.next(), end=" ")
In [ ]:
# 生成 [0, 1) 区间的随机数
gen = RandomGenerator(seed=42)
[gen.random() for _ in range(5)]

面向对象方法的优势¶

  1. 封装性：状态（seed）被封装在对象内部，不会污染全局命名空间
  2. 独立性：可以创建多个独立的生成器实例，互不干扰
  3. 可扩展性：易于添加新方法（如 random(), randint() 等）
  4. 可重用性：可以传递对象引用，其他代码无需了解内部实现
  5. 可维护性：相关代码组织在一起，易于理解和维护
In [ ]:
# 扩展示例：添加更多功能
class AdvancedRandomGenerator(RandomGenerator):
    def randint(self, a, b):
        """生成 [a, b] 区间的随机整数"""
        return a + self.next() % (b - a + 1)
    
    def choice(self, seq):
        """从序列中随机选择一个元素"""
        return seq[self.randint(0, len(seq) - 1)]
In [ ]:
adv_gen = AdvancedRandomGenerator(seed=42)

print("随机整数 (1-10):", [adv_gen.randint(1, 10) for _ in range(5)])
print("随机选择:", [adv_gen.choice(['A', 'B', 'C']) for _ in range(5)])

伪随机数生成小结¶

  • 伪随机数是通过确定性算法生成的数列，看起来随机但实际上是确定的
  • 线性同余生成器是最简单的伪随机数生成算法
  • 状态管理是随机数生成的核心问题
  • 面向对象方法相比全局变量具有封装性、独立性、可扩展性等优势

随机采样基础¶

  • 使用 NumPy 的 random 模块从各类对象中生成随机数据：
In [ ]:
# 从 1 到 6 的整数中随机采样（模拟掷骰子）
np.random.randint(1, 7)
In [ ]:
# 从列表中随机选择一个元素
np.random.choice([2, 3, 5, 7, 11])
In [ ]:
# 从字符串中随机选择一个字符
np.random.choice(list("THU"))
In [ ]:
# 从 a-z 中随机选择一个字母
np.random.choice(list(string.ascii_lowercase))
In [ ]:
# 生成 0 到 1 之间的随机浮点数
np.random.random()
  • 从任意列表集合（如颜色）中随机选取：
In [ ]:
colors = ['red', 'green', 'blue', 'orange', 'purple',
          'cyan', 'magenta', 'yellow', 'brown', 'pink']
print(f"随机选取的颜色: {np.random.choice(colors)}")

批量生成随机对象¶

  • 方法一：列表推导式
In [ ]:
[np.random.randint(1, 7) for _ in range(10)]
  • 方法二：利用 size 参数（推荐做法，性能更高）
In [ ]:
np.random.randint(1, 7, size=10)
  • 扩展：生成随机矩阵
In [ ]:
np.random.randint(1, 7, size=(10, 12))
In [ ]:
# 生成随机颜色矩阵并可视化
color_list = ['red', 'green', 'blue', 'orange', 'purple']
color_indices = np.random.randint(0, len(color_list), size=(10, 10))

fig, ax = plt.subplots(figsize=(6, 6))
cmap = ListedColormap(color_list)
ax.imshow(color_indices, cmap=cmap)
ax.set_title("随机颜色矩阵")
ax.set_xticks([])
ax.set_yticks([])
plt.show()

均匀采样 (Uniform Sampling)¶

定义：每个可能的对象被选中的概率相同。

示例：模拟抛硬币并用 Counter 统计结果：

In [ ]:
tosses = np.random.choice(["正面", "反面"], size=10000)
toss_counts = Counter(tosses)
toss_counts
In [ ]:
toss_counts["反面"]

Counter 返回类字典对象：键 $\to$ 频率/出现次数

In [ ]:
prob_tail = toss_counts["反面"] / len(tosses)
print(f"反面的概率 \u2248 {prob_tail:.4f}")

大数定律：随试验次数增加，实际频率无限趋近于期望概率 $1/2$。

非均匀采样：加权硬币¶

目标：模拟非对称概率 (例如：正面 $p = 0.7$，反面 $q = 0.3$)

思路：离散整数区间映射 (1-7 映射为正面，8-10 为反面)

In [ ]:
def simple_weighted_coin():
    if np.random.randint(1, 11) <= 7:
        return "正面"
    else:
        return "反面"

simple_weighted_coin()

推广：推广至任意概率 $p \in [0, 1]$

伯努利试验 (Bernoulli Trial)：生成 $[0, 1)$ 的均匀浮点数，判断其是否小于 $p$。

In [ ]:
np.random.random() < 0.314159

比较运算输出布尔值 (True/False)，将其封装为通用函数：

In [ ]:
def bernoulli(p):
    return np.random.random() < p

bernoulli(0.7)
In [ ]:
@interact(p=FloatSlider(min=0.0, max=1.0, step=0.01, value=0.7,
                        description='p:', readout_format='.2f'))
def show_bernoulli(p):
    results = Counter(bernoulli(p) for _ in range(1000))
    print(f"1000 次伯努利试验结果 (p={p:.2f}):")
    print(f"  True  (成功): {results.get(True, 0)}")
    print(f"  False (失败): {results.get(False, 0)}")

数据可视化：条形图¶

  • 目标：统计并可视化离散类别的出现频次
  • 示例：掷公平骰子 100,000 次
In [ ]:
rolls = np.random.randint(1, 7, size=100000)

朴素统计法：使用列表推导式或其他统计工具

In [ ]:
# 统计每个面出现的次数
counts = [np.sum(rolls == i) for i in range(1, 7)]
counts

条形图：适用于展示分类/离散数据

In [ ]:
fig, ax = plt.subplots(figsize=(8, 4))
ax.bar(range(1, 7), counts, alpha=0.5)
ax.axhline(y=len(rolls) / 6, linestyle='--', color='orange', label='期望值')
ax.set_title(f"掷骰子次数 = {len(rolls)}")
ax.set_xlabel("点数")
ax.set_ylabel("次数")
ax.set_ylim(0, len(rolls) / 3)
ax.legend()
plt.show()

概率密度与中心极限定理¶

场景：统计抛掷多个骰子的点数之和¶

In [ ]:
def roll_dice(n, sides=12):
    """掷 n 个 sides 面的骰子并返回点数之和"""
    return np.sum(np.random.randint(1, sides + 1, size=n))

trials = 10**6

分布形状的收敛¶

交互观察（改变骰子数量 $n$）：

  1. 原始直方图：随 $n$ 增大，由于聚合效应，形状趋向钟形曲线。
  2. 归一化 $y$ 轴：转换为密度，使得阴影区域总面积为 1。
  3. 标准化 $x$ 轴：数据去均值并除以标准差，使其收敛至标准正态分布。
In [ ]:
@interact(n=IntSlider(min=1, max=50, step=1, value=1, description='n:'))
def show_dice_analysis(n):
    # 向量化生成：一次性掷 trials 组骰子
    data = np.sum(np.random.randint(1, 13, size=(trials, n)), axis=1)
    sigma = np.std(data)
    mu = np.mean(data)
    normalised_data = (data - mu) / sigma

    fig, axes = plt.subplots(1, 3, figsize=(18, 4))

    # 图1: 原始直方图
    axes[0].hist(data, bins=200, alpha=0.5, color='lightsalmon')
    axes[0].set_title(f"原始直方图 (n = {n})")
    axes[0].set_xlabel("点数之和")
    axes[0].set_ylabel("频次")

    # 图2: 归一化 y 轴
    axes[1].hist(data, bins=50, density=True, alpha=0.5, color='lightsalmon')
    axes[1].set_title(f"归一化 y 轴 (n = {n})")
    axes[1].set_xlabel("点数之和")
    axes[1].set_ylabel("概率密度")
    axes[1].set_ylim(0, 0.05)

    # 图3: 归一化 x 轴 + 标准正态分布
    axes[2].hist(normalised_data,
                 bins=np.arange(-5 - 1/(2*sigma), 5, 1/sigma),
                 density=True, alpha=0.5, color='lightsalmon')
    x = np.linspace(-5, 5, 200)
    axes[2].plot(x, np.exp(-x**2 / 2) / np.sqrt(2 * np.pi),
                 'r-', linewidth=2, alpha=0.5,
                 label=r'$\frac{1}{\sqrt{2\pi}} e^{-x^2/2}$')
    axes[2].set_ylim(0, 0.41)
    axes[2].set_title(f"归一化 x+y 轴 (n = {n})")
    axes[2].set_xlabel("标准化值")
    axes[2].set_ylabel("概率密度")
    axes[2].legend()

    plt.tight_layout()
    plt.show()
    print(f"均值 = {mu:.2f}, 标准差 = {sigma:.2f}")

钟形曲线¶

image.png

从离散分布到连续分布¶

  • 直方图 (Histogram)：统计落入数据"箱" (bin) 中的频数，表现连续或密集离散数据。
  • 缩放漂移对策：
    • 归一化 ($y$ 轴)：总面积固定为 1。
    • 标准化 ($x$ 轴)：对齐中心 ($0$)，统一宽度 (除以标准差)。
  • 概率密度函数 (PDF)：在连续极限下，单点发生概率为 0；在区间 $[a, b]$ 内的积分面积代表落入该区间的概率。

直方图绘图选项解释¶

  • 核心参数：
    • density=True：归一化为密度图
    • bins=50：指定箱数（或输入边界数组）
    • histtype='step'：绘制无填充的阶梯线条
  • 通用样式：
    • alpha：透明度 (0 不可见 ~ 1 不透明)
    • color：支持名称 ('red')、十六进制 ('#FF5733') 或 RGB
    • linewidth (lw)：线宽设定
In [ ]:
data_example = np.sum(np.random.randint(1, 13, size=(trials, 10)), axis=1)
normalised_example = (data_example - np.mean(data_example)) / np.std(data_example)

fig, ax = plt.subplots(figsize=(8, 4))
ax.hist(normalised_example, bins=50, density=True, histtype='step', linewidth=2)
ax.set_title("阶梯直方图示例")
plt.show()

卡方分布¶

  • 卡方分布属性验证：自由度为 $k$ 的卡方分布，等价于 $k$ 个标准正态随机变量的平方和。
  • 交互验证（调节自由度 dof）：
    1. 左图：使用 np.random.chisquare 直接采样。
    2. 右图：手动生成 dof 个标准正态变量，求平方和。
  • 结论：两者分布的概率密度表现完全一致。
In [ ]:
@interact(dof=IntSlider(min=1, max=50, step=1, value=1, description='自由度:'))
def show_chisq(dof):
    chisq_data = np.random.chisquare(dof, size=100000)
    manual_chisq = np.sum(np.random.randn(100000, dof)**2, axis=1)

    fig, axes = plt.subplots(1, 2, figsize=(14, 4))

    # 左图: numpy 内置卡方分布
    axes[0].hist(chisq_data, bins=100, density=True, alpha=0.5)
    axes[0].set_xlim(0, 10 * np.sqrt(dof))
    axes[0].set_title(f"卡方分布 (自由度 = {dof})")
    axes[0].set_xlabel("值")
    axes[0].set_ylabel("概率密度")

    # 右图: 手动构造 (标准正态平方和)
    axes[1].hist(manual_chisq, bins=100, density=True, alpha=0.5)
    axes[1].set_xlim(0, 10 * np.sqrt(dof))
    axes[1].set_title(f"标准正态平方和 (自由度 = {dof})")
    axes[1].set_xlabel("值")
    axes[1].set_ylabel("概率密度")

    plt.tight_layout()
    plt.show()

统计学三大分布对比表¶

分布类型 定义 (Definition) 主要性质 (Properties) 典型例子 (Examples)
均匀分布
(Uniform)		 在给定区间 $[a, b]$ 内，随机变量取到任何值的概率密度都是相等的。记作 $X \sim U(a, b)$。 1. 等可能性：概率密度函数在区间内是平坦的常数。
2. 均值中心化：期望值位于区间正中间，即 $(a+b)/2$。
3. 独立性：常用于描述没有任何先验偏好的随机状态。		 等待每10分钟一班的公交车所需的时间、理想状态下掷骰子的结果（离散型）、计算机随机数生成器。
正态分布
(Normal)		 描述数据围绕一个中心值（均值）对称分布的连续型概率分布。记作 $X \sim N(\mu, \sigma^2)$。 1. 对称性：关于均值 $\mu$ 对称，呈钟形。
2. 中心极限定理：大量独立随机变量之和趋向于正态分布。
3. 68-95-99.7法则：约99.7%的数据落在 $\mu \pm 3\sigma$ 内。		 人的身高、成年人的体重、标准化考试（如SAT）分数、工业产品的尺寸误差。
卡方分布
(Chi-square)		 $k$ 个独立的标准正态分布随机变量的平方和所服从的分布。记作 $X \sim \chi^2(k)$。 1. 非负性：取值范围始终 $\ge 0$。
2. 偏态性：曲线右偏（长尾在右），随自由度 $k$ 增大而趋于对称。
3. 加法性：两个独立卡方变量之和仍服从卡方分布。		 拟合优度检验（验证骰子是否均匀）、独立性检验（研究吸烟与肺癌是否相关）、估计总体的方差。

本讲小结¶

  • 伪随机数生成的原理
    • 面向对象方法的使用场景
  • 随机采样与可视化
  • 三种重要的分布
    • 均匀分布
    • 正态分布
    • 卡方分布