異なる時間間隔で異なるレートを持つ連続時間マルコフ連鎖のシステムをモデル化しようとしています。
このように各期間のレートマトリックスを作成します
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)
}
（興味がある場合は、移行を伴う2分割システムでの合体密度のモデル化です）
レート、cs および ms は期間によって異なるため、各期間のレート マトリックスを作成し、次に各期間の遷移確率マトリックスを作成します。
2 つの期間を使用して、このようにレートを指定できます
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 までの 2 番目のレート セットを使用したいとします。
そこで、第 1 期と第 2 期の率行列と、第 1 期と第 2 期を通じて状態 a から b に移行する確率遷移行列の表が必要です。
mlply(rates, make.rate.matrix)
2 つのレート マトリックスのリストが表示されます。レート マトリックスを簡単に検索できるテーブルが必要な場合は、次のようにします。
> 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[,,1] が P1 で、P[,,2] が P2 であるテーブル P を取得するにはどうすればよいですか?