如何利用R中的Galton过程求出灭绝概率=1？

Question

我正在用几何分布模拟一个基本的Galton过程(GWP) .我用这个来找出每一代人的灭绝概率。我的问题是，我如何找到灭绝概率等于1的世代？

例如，我可以为GWP创建如下函数：

# Galton-Watson Process for geometric distribution
GWP <- function(n, p) {
  Sn <- c(1, rep(0, n))
  
  for (i in 2:(n + 1)) {
    Sn[i] <- sum(rgeom(Sn[i - 1], p))
  }
  return(Sn)
}

其中，n是代数。

如果我设定几何分布参数p=0.25.然后，为了计算第10代的灭绝概率，我只需要这样做：

N <- 10 # Number of elements in the initial population. 
GWn <- replicate(N, GWP(10, 0.25)[10])
probExtinction <- sum(GWn==0)/N 
probExtinction

这会给我第十代人灭绝的可能性..。为了找到每一代的灭绝概率，我必须在创建GWn时更改索引值(对应的代数).但我要做的是找出灭绝概率= 1的世代。

对于我如何解决这个问题，有什么建议吗？

Answer 1

原则上，我可以告诉你如何解决这个问题，但我建议你可能会遇到一些困难(如果你已经知道我要说的一切，就把它当作给下一个读者的建议吧……)

• theoretically，Galton-Watson过程的灭绝概率永远不会精确到1(除非prob==1，或者在无限时间内)，当然，对于任何给定的复制和随机数种子，您可以计算出所有血统灭绝的第一个时间点(如果有的话)。这将是高度可变的运行，取决于随机数种子.
• 的分布是极不平衡的；没有立即灭绝的血统将持续很长的时间.

我用两种方法修改了GWP函数，使其更有效：(1)在世系灭绝时停止模拟；(2)将几何偏差之和替换为单个负二项式偏差(见here)。

GWP <- function(n, p) {
  Sn <- c(1, rep(0, n))
  for (i in 2:(n + 1)) {
    Sn[i] <- rnbinom(1, size=Sn[i - 1], prob=p)
    if (Sn[i]==0) break ## extinct, bail out
  }
  return(Sn)
}

现在的基本策略是：(1)运行一段时间的模拟，保持整个轨迹；(2)计算每一代的灭绝概率；(3)找到第一代p==1。

set.seed(101)
N <- 10 # Number of elements in the initial population.
maxgen <- 100
GWn <- replicate(N, GWP(maxgen, 0.5), simplify="array")
probExtinction <- rowSums(GWn==0)/N
which(probExtinction==1)[1]

(如果要从第0代开始索引，则从最后的结果中减去1。)在这种情况下，答案是NA，因为有1/10的血统能够存活下来(而且确实变得非常大，因此它可能几乎永远存在)。

plot(0:maxgen, probExtinction, type="s")    ## plot extinction probability
matplot(1+GWn,type="l",lty=1,col=1,log="y") ## plot lineage sizes (log(1+x) scale)

​

​

​

​

## demonstration that (sum(rgeom(n,...)) is equiv to rnbinom(1,size=n,...)
nmax <- 70
plot(prop.table(table(replicate(10000, sum(rgeom(10, prob=0.3))))),
     xlim=c(0,nmax))
points(0:nmax,dnbinom(0:nmax, size=10, prob=0.3), col=2,pch=16)