R语言模拟中随机数异常导致实体活跃时长与预期不符的问题

核心误解：单期活跃概率 vs 累积存活概率与期望时长 #

你看到的模拟结果是符合理论预期的，问题出在混淆了「单期活跃概率」和「存活到第n年的概率」以及「期望活跃时长」的计算逻辑：

• 单期活跃概率是指在已经存活到前一年的前提下，当年保持活跃的概率；
• 存活到第t年的概率是前t年所有单期活跃概率的乘积（必须每年都满足存活条件）；
• 期望活跃时长是所有「存活到第t年的概率」的总和（即 ( E[T] = \sum_{t=1}^{60} P(T \geq t) )，其中 ( P(T \geq t) ) 是实体至少存活到第t年的概率）。

以你的数据计算：

• 存活到第10年的概率 = ( 0.9429 \times 0.9357 \times 0.9277 \times ... \times 0.8411 \approx 35% )
• 把前若干年的 ( P(T \geq t) ) 累加：( 0.9429 + 0.882 + 0.818 + 0.752 + 0.684 + 0.615 + 0.546 + 0.479 + ... )，总和恰好接近8年，和你的模拟结果一致。

关于随机数的疑问 #

runif(60) 生成的是均匀分布随机数，0-1区间内每个数值段的出现概率均等，0.9以上的数值出现概率为10%，前10个位置出现高值是正常的随机波动，并非随机数生成异常。你可以通过以下代码验证：

# 生成10000个随机数，统计0.9以上的比例
rand <- runif(10000)
mean(rand > 0.9) # 结果应接近0.1
hist(rand) # 直方图应呈现均匀分布

代码优化：避免内存崩溃 #

你之前存储全量随机数矩阵导致内存问题，建议改为每轮模拟仅计算活跃时长并存储结果，无需保存所有随机数。示例代码如下：

# 定义60年的活跃概率向量（替换为你实际的概率数据）
active_probs <- c(0.9429, 0.9357, 0.9277, 0.9188,
                  # 填充第5-9年的概率（按你的趋势补充）
                  0.909, 0.899, 0.888, 0.877, 0.865,
                  0.8411, # 第10年
                  # 第11-20年（示例，替换为你的实际值）
                  seq(0.83, 0.6, length.out = 10),
                  # 第21-60年（示例，替换为你的实际值）
                  seq(0.59, 0.1, length.out = 40))

# 单轮模拟函数
calc_lifetime <- function(probs) {
  lifetime <- 0
  is_active <- TRUE
  for (p in probs) {
    if (!is_active) break
    if (runif(1) < p) {
      lifetime <- lifetime + 1
    } else {
      is_active <- FALSE
    }
  }
  lifetime
}

# 批量模拟（设置种子保证结果可重复）
set.seed(123)
n_simulations <- 1000
lifetimes <- replicate(n_simulations, calc_lifetime(active_probs))

# 查看平均活跃时长
mean(lifetimes)

此代码仅存储每轮的活跃时长，内存占用极低，可轻松完成上万次模拟。

内容的提问来源于stack exchange，提问作者Mr.Anugar