在R中创建分支进程

oprakyz7 于 2023-10-13 发布在其他

关注(0)|答案(1)|浏览(114)

我试图模拟一个具有负二项后代分布的分支过程（或Galton-Watson过程）。分支过程的一个例子如下所示，从Cui, et al中任意选取。
简而言之，在传染病的背景下，每个分支过程都从一个传染性个体（第0代）开始，它们感染了由某些后代分布定义的随机数量的继发病例（我使用的是具有参数的负二项后代（平均值= r 0，离散度=k））。后续世代中的每个病例感染随机数量的病例，依此类推（i.i.d.）。

我想象数据布局是一个矩阵，类似于这个，它代表了图中的前三代：

#Example of data layout from first three generations in figure
x<-matrix(c(0,1,1,2,2,2,1,1,2,1,2,3,2,2,1,1,0,3),ncol = 3,nrow = 6)
colnames(x)<-c("Generation","Case Number","Secondary Cases (Z)")

我不确定如何创建一个函数来梳理单个Z值的数量。我写了一个函数来模拟一个负二项分支过程，但是只对每一代的案例总数进行求和（即：在下图中，G（0）=1，G（1）=2，G（2）=3，G（3）=4）。此外，在这个函数中，我必须定义代数（“n”），更准确的代码将允许它自然地继续下去，直到一代中有0个案例。
下面的函数为我尝试做的事情提供了一些好处，但功能有限。

#Single branching process
nbbp<-function(n, r0, k){
  # initialize return vector
  Z<-integer(n)
  #Initiate branching process with 1 index case in generation 1
  Z[1]<-1
  for (i in seq_len(n)[-1]){
    if(Z[i-1]>0) {
         x<-rnbinom(Z[i-1], size=k, mu=r0)
         Z[i] <- sum(x)
    }
  }
  return(Z)
}
#Simulate multiple BP 
nbbp.ind<-function(num_sim,n,r0,k){
    x<-replicate(num_sim, nbbp(n,r0,k))
    }

r

来源：https://stackoverflow.com/questions/58434856/creating-a-branching-process-in-r

1条答案

按热度按时间

a64a0gku1#

一个广义分支过程函数，可以接受任何分布（distr），具有任意数量的参数，可以写为：

bp <- function(distr, gens = 20, init.size = 1,  ...){  
  Z <- list() 
  Z[[1]] <- init.size 
  i <- 1 
  while(sum(Z[[i]]) > 0 && i <= gens) { 
    Z[[i+1]] <- distr(sum(Z[[i]]), ...) 
    i <- i+1 
  } 
  return(Z)
}

注意gens设置了防止溢出的最大代数（如果需要，可以设置为Inf），init.size允许您更改第0代中的初始事件数。
这将返回一个向量列表，其中列表中的每个位置位置代表一个世代，而向量中的每个位置代表一个后代：

set.seed(12345)

# Negative binomial distribution
bp(distr = rnbinom, mu = 0.9, size = 0.25, gen = Inf, init.size = 1)
# [[1]]
# [1] 1
# 
# [[2]]
# [1] 1
# 
# [[3]]
# [1] 0

# Poisson distiribution
bp(distr = rpois, lambda = 0.9, gen = Inf, init.size = 1)

# [[1]]
# [1] 1
# 
# [[2]]
# [1] 1
# 
# [[3]]
# [1] 4
# 
# [[4]]
# [1] 0 0 1 0
# 
# [[5]]
# [1] 0

赞(0）回复(0）举报 2023-10-13

我来回答

在R中创建分支进程

1条答案

相关问题

热门标签

最新问答