探秘离散时间更新过程：固定配额下的稳态年龄分布研究

1. 问题背景
考虑一个包含 $n$ 个位置的离散时间动态系统。在每一期，系统独立且均匀地选取 $m$ 个不同的位置（无放回抽样）。系统的演化规则如下：
年龄增长：将所有 $n$ 个位置的年龄加 1。
重置操作：将选中的 $m$ 个位置的年龄重置为 0。
经过长时间演化，系统达到稳态。我们关注的问题是：在稳态下，将 $n$ 个位置按年龄从小到大排序，年龄恰好为 $k$ 的位置个数 $N_k$ 服从什么分布？
2. 理论推导
2.1 基础性质与 $k=0$ 的特例
首先考察年龄为 0 的位置个数 $N_0$。根据操作规则，每期结束后，被选中的 $m$ 个位置年龄被置为 0。因此，在任何时刻，年龄为 0 的位置数量是确定性的。
$$ N_0 = m \quad \text{(概率为 1)} $$
这是一个常数分布，也是系统演化的“源”。
2.2 一般情形：链式超几何分布
对于 $k \ge 1$，年龄为 $k$ 的位置集合是上一期年龄为 $k-1$ 的位置集合的子集。具体来说，一个在第 $t$ 期年龄为 $k-1$ 的位置，只有在第 $t+1$ 期未被选中时，其年龄才会变为 $k$。
由于每期有 $m$ 个位置被选中，这意味着有 $n-m$ 个位置未被选中（年龄增加）。这构成了一个无放回抽样过程。已知 $N_{k-1}$ 的情况下，$N_k$ 的条件分布服从超几何分布。
设 $q = 1 - \frac{m}{n}$ 为单个位置在一期内未被选中的概率。
情形一：$k=1$ 的分布
$N_1$ 表示上一期年龄为 0 的 $m$ 个位置中，本期未被选中的个数。这是一个标准的超几何分布：
$$ N_1 \mid N_0 \sim \text{Hypergeometric}(n, N_0, n-m) $$
由于 $N_0=m$，其概率质量函数（PMF）为：
$$ P(N_1 = x) = \frac{\binom{m}{m-x}\binom{n-m}{x}}{\binom{n}{m}} $$
利用组合对称性 $\binom{m}{m-x} = \binom{m}{x}$，公式可简化为：
$$ P(N_1 = x) = \frac{\binom{m}{x}\binom{n-m}{x}}{\binom{n}{m}} $$
其中 $x$ 的取值范围为 $[\max(0, 2m-n), m]$。
情形二：$k > 1$ 的递推分布
对于 $k \ge 2$，分布由上一期的状态决定，形成一条马尔可夫链。我们需要计算边缘分布。
设 $P(N_{k-1} = y)$ 为上一期的分布，则 $N_k$ 的分布由全概率公式给出：
$$ P(N_k = x) = \sum_{y} P(N_{k-1} = y) \cdot P(N_k = x \mid N_{k-1} = y) $$
其中的条件概率核心项为：
$$ P(N_k = x \mid N_{k-1} = y) = \frac{\binom{y}{y-x}\binom{n-y}{m-y+x}}{\binom{n}{m}} $$
该公式的物理含义是：从 $n$ 个位置中选出 $m$ 个重置，其中从 $y$ 个“潜在候选”（上期年龄 $k-1$）中选出了 $y-x$ 个（使其归零，剩下 $x$ 个存活），从剩余 $n-y$ 个位置中选出了 $m-y+x$ 个。
2.3 期望与渐近性质
虽然精确分布需要递推计算，但其期望具有简洁的封闭形式。
利用全期望公式和超几何分布期望公式：
$$ E[N_k \mid N_{k-1}] = N_{k-1} \left(1 - \frac{m}{n}\right) = q N_{k-1} $$
迭代可得稳态期望：
$$ E[N_k] = m q^k = m \left(1 - \frac{m}{n}\right)^k $$
这印证了每个位置年龄服从几何分布的直觉。当 $n$ 较大时，无放回抽样近似于有放回抽样，$N_k$ 的分布近似于二项分布：
$$ N_k \stackrel{\text{approx}}{\sim} \text{Binomial}\left(m, q^k\right) $$
3. 实测验证 ($n=33, m=6$)
为了验证上述理论，我们选取 $n=33, m=6$ 的系统进行仿真测试。我们计算了理论值并与实测数据进行了对比。
3.1 $k=1$ 的分布验证
理论模型：$N_1 \sim \text{Hypergeometric}(33, 6, 27)$。
结论：理论值与实测值高度吻合，峰值位置（众数）均在 $x=5$。
3.2 $k > 1$ 的分布验证
我们利用 $k=1$ 的精确分布，通过卷积公式递推计算 $k=2, 3, 4$ 的理论分布。
$k=2$ 分布验证：
$k=4$ 分布验证（经过多次递推，误差潜在累积）：
结论：随着 $k$ 的增加，虽然方差效应累积导致微小偏差，但整体分布形态、峰值位置依然与理论预测保持一致。这证明了链式超几何分布模型的准确性。
4. 总结
对于固定配额 $m$ 的离散更新系统，其稳态年龄分布并非简单的独立同分布，而是一个具有强相关性的链式结构：
$N_0$ 是常数 $m$。
$N_1$ 服从精确的超几何分布。
$N_k (k>1)$ 服从由递推公式定义的边缘分布。
该模型精确地描述了“无放回抽样”带来的负相关性，为排队论、缓存替换策略等领域提供了有力的理论支撑。

以下是计算该稳态年龄分布的 Python 代码。代码中定义了一个通用的 calculate_distribution 函数，可以计算任意 kk 值的精确概率分布，并包含了针对 n=33,m=6n=33,m=6 案例的验证部分。

import math
from scipy.special import comb
import numpy as np

def calculate_distribution(n, m, k_max):
    """
    计算稳态下年龄为 k (0 到 k_max) 的位置个数分布。
    
    参数:
    n: 总位置数
    m: 每期重置的位置数
    k_max: 需要计算的最大 k 值
    
    返回:
    一个列表，包含每个 k 对应的概率分布数组。
    """
    
    # 结果存储列表：dists[k] 表示 N_k 的概率分布向量
    # 向量索引 i 代表个数 x=i 的概率
    dists = []
    
    # 1. 计算 k=0 的分布
    # N_0 恒等于 m (确定性分布)
    # 向量长度需足以容纳最大可能值。对于 N_k，最大值始终不超过 m。
    vec_len = m + 1
    prob_0 = np.zeros(vec_len)
    prob_0[m] = 1.0
    dists.append(prob_0)
    
    # 2. 递推计算 k=1 到 k_max 的分布
    for k in range(1, k_max + 1):
        # 上一期 k-1 的分布
        prob_prev = dists[k-1]
        # 本期 k 的分布初始化
        prob_curr = np.zeros(vec_len)
        
        # 遍历上一期的所有可能状态 y
        for y in range(vec_len):
            p_y = prob_prev[y]
            
            # 如果该状态概率为0，跳过
            if p_y == 0:
                continue
            
            # 遍历本期可能的个数 x
            # N_k = x 意味着上一期 y 个位置中有 x 个未被选中 (存活)
            # 未被选中的个数 x 范围: [max(0, y-m), min(y, n-m)]
            # 但实际上 N_k 不可能超过 m (因为 N_0=m 是源头)
            
            min_x = max(0, y - m)
            max_x = min(y, n - m)
            
            # 计算条件概率 P(N_k=x | N_{k-1}=y)
            # 这是一个超几何分布：
            # 总数 n, "坏"品 y (会被重置), "好"品 n-y
            # 抽取 n-m 个 (未被选中的个数)
            # 我们需要其中包含 x 个 "坏"品 (即 y 中有 x 个未被重置)
            for x in range(min_x, max_x + 1):
                # 超几何分布公式:
                # C(坏品总数, 抽中的坏品数) * C(好品总数, 抽中的好品数) / C(总数, 抽取总数)
                # = C(y, x) * C(n-y, (n-m)-x) / C(n, n-m)
                
                numerator = comb(y, x, exact=True) * comb(n - y, (n - m) - x, exact=True)
                denominator = comb(n, n - m, exact=True)
                
                p_cond = numerator / denominator
                
                # 累加概率 (全概率公式)
                prob_curr[x] += p_y * p_cond
                
        dists.append(prob_curr)
        
    return dists

# --- 验证部分：n=33, m=6, k=1至4 ---

n = 33
m = 6
k_max = 4

# 计算分布
distributions = calculate_distribution(n, m, k_max)

print(f"参数设置: n={n}, m={m}")
print("-" * 50)

# 实测数据 (用户提供的)
observed_data = {
    1: [0, 0.03, 0.5, 4.73, 23.9, 43.66, 27.16],
    2: [0.12, 0.86, 6.6, 23.1, 36.24, 26.515, 6.59],
    3: [0.71, 4.94, 18.8, 34.45, 27.8, 11.56, 1.74],
    4: [1.98, 12.74, 29.71, 32.19, 18.38, 4.79, 0.21]
}

# 格式化输出对比
for k in range(1, k_max + 1):
    probs = distributions[k]
    print(f"年龄 k={k} 的分布对比:")
    print(f"{'个数 x':<8} | {'理论值 (%)':<12} | {'实测值 (%)':<12} | {'偏差 (%)':<10}")
    print("-" * 50)
    
    obs_vals = observed_data[k]
    for x in range(m + 1):
        theory = probs[x] * 100
        obs = obs_vals[x]
        diff = theory - obs
        print(f"{x:<8} | {theory:<12.2f} | {obs:<12.2f} | {diff:<10.2f}")
    print("-" * 50)
    print()