1

我正在尝试对连续时间马尔可夫链系统进行建模，其中在不同的时间间隔内我有不同的速率。

我像这样为每个时间段建立一个速率矩阵

make.rate.matrix <- function(c1, c2, m12, m21) {
  matrix(
    c(# State 1: lineages in different populations
      -(m12+m21), m21, m12, 0,
      # State 2: both lineages in population 1
      2*m12, -(2*m12+c1), 0, c1,
      # State 3: both lineages in population 2
      2*m21, 0, -(2*m21+c2), c2,
      # State 4: coalesced (catches both populations; absorbing)
      0, 0, 0, 0),
    byrow = TRUE, nrow=4)
}

（如果您有兴趣，它正在模拟具有迁移的二元系统中的聚结密度）

cs 和 ms 的速率在不同的时间段内有所不同，因此我想为每个时间段构建一个速率矩阵，然后为每个时间段构建一个转换概率矩阵。

有两个时期，我可以指定这样的费率

rates <- data.frame(c1 = c(1,2), c2 = c(2,1), m12 = c(0.2, 0.3), m21 = c(0.4, 0.2))

我想使用从时间 0 到 t 的第一组费率和从时间 t 到 s 的第二组费率。

所以我想要一个第一和第二周期的速率矩阵表，以及从状态 a 移动到 b 通过第一和第二周期的概率转换矩阵。

mlply(rates, make.rate.matrix)

给了我两个比率矩阵的列表，如果我想要一个可以轻松查找比率矩阵的表格，我可以做类似的事情

> xx <- array(unlist(mlply(rates, make.rate.matrix)), dim=c(4,4,2))
> xx[,,1]
     [,1] [,2] [,3] [,4]
[1,] -0.6  0.4  0.2    0
[2,]  0.4 -1.4  0.0    1
[3,]  0.8  0.0 -2.8    2
[4,]  0.0  0.0  0.0    0
> xx[,,2]
     [,1] [,2] [,3] [,4]
[1,] -0.5  0.2  0.3    0
[2,]  0.6 -2.6  0.0    2
[3,]  0.4  0.0 -1.4    1
[4,]  0.0  0.0  0.0    0

然后我可以得到概率转换矩阵，如

> library(Matrix)
> t <- 1; s <- 2
> P1 <- expm(xx[,,1] * t)
> P2 <- expm(xx[,,2] * (s - t))

但我不知道如何获得这些表格，就像我可以为费率矩阵获得的表格一样。

我觉得我应该能够轻松到达那里，但我不知道如何到达那里。

如何获得一个表 P，其中 P[,,1] 是 P1，P[,,2] 是 P2？

4

0 回答 0