抛硬币问题

我们考虑一枚硬币（正面概率为 $p$，反面概率为 $1−p$）。不断投掷，直到首次出现连续 $k$ 个正面为止。我们要计算总投掷次数的数学期望。

下面我们使用马尔科夫链吸收理论讨论这个问题。

数学建模

定义状态：

  1. $S_j$：当前结尾连续 $j$ 个正面（$0 \leqslant j \leqslant k$）.
  2. $S_k$ 为吸收态，表示目标事件（连续 $k$ 个正面）已发生。

转移规律：

  • 在 $S_j$（$0 \leqslant j \leqslant k - 1$）：
    1. 掷正面（概率 $p$）：进入 $S_{j+1}$.
    2. 掷反面（概率 $1-p$）：回到 $S_{0}$.

因此我们得到一个吸收马尔科夫链，非吸收态为 $S_0, \cdots,S_{k-1}$。

矩阵推导

将状态划分为非吸收态 $S_0, \cdots,S_{k-1}$ 和吸收态 $S_{k}$.

转移矩阵写成标准分块形式：

$$
\begin{pmatrix}
Q & R \
0 & I
\end{pmatrix}
$$

其中 $Q$ 为 $k \times k$ 的非吸收转移矩阵，$R$ 为进入吸收态的部分。

构造矩阵 Q

例如 $k=3$ 时，非吸收态为 $(S_0, S_1, S_2)$，对应转移概率矩阵：

$$
Q =
\begin{pmatrix}
1-p & p & 0 \
1-p & 0 & p \
1-p & 0 & 0
\end{pmatrix}
$$

一般的 $k$ 时，$Q$ 的结构为：

  1. 每一行的第 0 列为 $1-p$（回到 $S_0$）。
  2. 对角线右上方元素为 $p$ （从 $S_j \rightarrow S_{j+1}$）
  3. 其余为0.

基本矩阵 N

定义：
$$
N = (I - Q)^{-1} = \sum_{t=0}^\infty Q^t
$$

其含义：$N_{ij}$ 是从状态 $S_i$ 出发，在到达吸收前平均会访问 $S_j$ 的次数。

期望停时向量

期望投掷次数向量：
$$
t = N \mathbb{1}
$$

其中 $\mathbb{1}$ 为全 1 列向量。于是
$$
t_i = \sum_{j=0}^{k-1} N_{ij},\ \text{即从状态} S_i ​\text{出发的期望投掷次数}.
$$

经过代数化简（可利用矩阵几何级数展开或递推关系），得到：
$$
t_j = E_j = \frac{1 - p^{k-j}}{(1-p)p^k}
$$

总结

通过阵法，我们得到解：
$$
\mathbb{E}[T] = \frac{1 - p^{k-j}}{(1-p)p^k}
$$
​

对公平硬币（$p=1/2$），该式化简为
$$
\mathbb{E}[T] = 2^{k+1} - 2
$$

python 代码模拟

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
import random
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams

# 设置字体以支持中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
plt.rcParams['axes.unicode_minus'] = False   # 正常显示负号

# 单次模拟
def simulate_once():
    count = 0
    consecutive_heads = 0
    while consecutive_heads < 3:
        flip = random.choice(['H', 'T'])  # 正面 H，反面 T
        count += 1
        if flip == 'H':
            consecutive_heads += 1
        else:
            consecutive_heads = 0
    return count

# 多次模拟
trials = 1000000
results = [simulate_once() for _ in range(trials)]

# 计算期望值
expected_value = np.mean(results)
print(f"模拟得到的投掷总次数的数学期望约为：{expected_value:.4f}")

# 绘制频率图
plt.figure(figsize=(10, 6))
plt.hist(results, bins=range(min(results), max(results) + 1), density=True, color='skyblue', edgecolor='black')
plt.axvline(expected_value, color='red', linestyle='dashed', linewidth=1.5, label=f'期望 ≈ {expected_value:.2f}')
plt.title("投掷至连续三次正面所需次数的频率分布图")
plt.xlabel("投掷次数")
plt.ylabel("相对频率")
plt.legend()
plt.grid(True)
plt.show()