R语言基于几何分布Galton-Watson过程如何找到灭绝概率为1的世代

核心逻辑 #

你要找的是蒙特卡洛模拟下，首次出现全部抽样种群均已灭绝的世代位置。原有写法存在样本量过小、逐代统计效率低、易误判结果的问题，可以按下面的流程优化落地。

具体实现步骤 #

• 理清原有代码的索引对应关系：你写的GWP函数返回的向量中，Sn[1]对应初始第0代（种群大小固定为1），Sn[k]对应第k-1代的种群规模，之前取索引10得到的实际是第9代的结果。
• 扩大模拟样本量：原来设置N=10的模拟规模太小，随机波动极强，很容易因为单次抽样的巧合，误判灭绝概率达到1，建议至少把样本量提到1000-2000的量级，让经验概率的结果足够稳定。
• 批量存储全世代模拟结果：每次运行GWP函数本来就会返回所有世代的种群大小，不需要逐次修改索引取单个世代的值，把所有模拟的全长度结果存成矩阵，就能一次性算出所有世代的灭绝概率。
• 自动定位目标世代：从第1代开始逐代扫描，找到第一个满足“该代所有模拟样本的种群值为0”的位置，就是你要找的、经验灭绝概率等于1的世代。

优化后的可运行代码 #

首先对原有GWP函数做提速优化，对已经灭绝的种群跳过后续无意义的随机数计算：

# 优化版GWP函数
GWP <- function(n, p) {
  Sn <- c(1, rep(0, n))
  for (i in 2:(n + 1)) {
    # 上一代种群已经灭绝的话，后续世代直接保持0，无需生成随机数
    if (Sn[i-1] == 0) {
      Sn[i] <- 0
      next
    }
    Sn[i] <- sum(rgeom(Sn[i - 1], p))
  }
  return(Sn)
}

执行批量模拟、逐代计算灭绝概率、自动定位目标世代：

p <- 0.25          # 几何分布参数
max_gen <- 100     # 预先设置足够大的最大观测世代数
N <- 2000          # 模拟样本量

# 生成N次模拟的全世代结果，转为矩阵后每一行对应一次完整模拟
sim_matrix <- t(replicate(N, GWP(max_gen, p)))
# 逐列计算每一代的灭绝概率（种群值为0的样本占比）
extinction_prob <- colMeans(sim_matrix == 0)
# 跳过第0代，找到第一个灭绝概率为1的位置，转换为实际世代数
first_hit_gen <- which(extinction_prob == 1)[1] - 1

# 输出结果
cat(sprintf("首次经验灭绝概率达到1的世代为第%d代，对应该代灭绝概率值为%.2f", 
            first_hit_gen, extinction_prob[first_hit_gen + 1]))

关键注意点 #

• 先校验参数适用场景：当你用p=0.25的参数时，几何分布后代数的期望为(1-p)/p = 3，属于超临界GWP，理论上存在固定的正概率让种群永久存续，长期灭绝概率会收敛到一个小于1的常数，永远不会严格等于1。如果你的目标是观测灭绝概率收敛到1的过程，需要把p调整到大于0.5的区间（次临界，后代期望小于1），才会出现最终全种群渐近灭绝的结果。
• 灵活调整预设世代数：如果运行代码后发现max_gen设置太小，也就是最后一代的灭绝概率还小于1，只需要把max_gen的数值调大，重新运行直到找到第一个全灭绝的世代即可。
• 明确结果的统计属性：样本量越大，首次摸到灭绝概率为1的世代会越靠后——因为样本量越大，越容易抽到小概率的未存续样本，经验值的结果越贴近理论上的渐近过程，不存在一个绝对固定的“真实世代值”。

内容的提问来源于stack exchange，提问作者Electrino