当值随时间变化时求幂矩阵

Question

我正在使用带有 R 编程矩阵的马尔可夫模型来解决健康经济问题。

我有一个矩阵，我想对其进行 200 次取幂。矩阵中的值随时间变化。例如，位置 x12 上的值在第一个矩阵 0.1 中，第二个 0.2 和第三个 0.3 中。我不想写下 200 个矩阵，而是写下一个矩阵并将其乘以 200。对于 x12：我可以使用具有 200 个值的向量并且矩阵将使用该向量吗？

对于常量值，没问题：

############### Imatinib Base-line Strategy ##########################
trans_matrix_dasa_no2nd <- matrix (,nrow=3,ncol=3)

trans_matrix_dasa_no2nd[1,] <- c(0.9,0.05,0.05)
trans_matrix_dasa_no2nd[2,] <- c(0,0.8,0.2)
trans_matrix_dasa_no2nd[3,] <- c(0,0,1)

cycle_patient_distribution_dasa_no2nd <- matrix (,nrow=2,ncol=3)

for (i in 1:2){cycle_patient_distribution_dasa_no2nd[i,] <- 
           initial_patient_distribtion%*%(trans_matrix_dasa_no2nd%^%i)}

cycle_patient_distribution_dasa_no2nd
     [,1] [,2] [,3]
[1,]  900   50   50
[2,]  810   85  105

但是如果我使用向量，它就不再起作用了：

############### Imatinib Base-line Strategy ##########################
aaa <- c(0.9,0.5)
bbb <- c(0.05, 0.4)
ccc <- c(0.05, 0.1)

trans_matrix_dasa_no2nd <- matrix (,nrow=3,ncol=3)

trans_matrix_dasa_no2nd[1,] <- c(aaa,bbb,ccc)
Error in trans_matrix_dasa_no2nd[1, ] <- c(aaa, bbb, ccc) : 
  number of items to replace is not a multiple of replacement length
trans_matrix_dasa_no2nd[2,] <- c(0,0.8,0.2)
trans_matrix_dasa_no2nd[3,] <- c(0,0,1)

cycle_patient_distribution_dasa_no2nd <- matrix (,nrow=2,ncol=3)

for (i in 1:2){cycle_patient_distribution_dasa_no2nd[i,] <- 
         initial_patient_distribtion%*%(trans_matrix_dasa_no2nd%^%i)}

cycle_patient_distribution_dasa_no2nd
     [,1] [,2] [,3]
[1,]   NA   NA   NA
[2,]   NA   NA   NA

如果我将i引入矩阵，结果不正确：

############### Imatinib Base-line Strategy ##########################
aaa <- c(0.9,0.5)
bbb <- c(0.05, 0.4)
ccc <- c(0.05, 0.1)

trans_matrix_dasa_no2nd <- matrix (,nrow=3,ncol=3)

trans_matrix_dasa_no2nd[1,] <- c(aaa[i],bbb[i],ccc[i])
trans_matrix_dasa_no2nd[2,] <- c(0,0.8,0.2)
trans_matrix_dasa_no2nd[3,] <- c(0,0,1)

cycle_patient_distribution_dasa_no2nd <- matrix (,nrow=2,ncol=3)

for (i in 1:2){cycle_patient_distribution_dasa_no2nd[i,] <- 
      initial_patient_distribtion%*%(trans_matrix_dasa_no2nd%^%i)}

cycle_patient_distribution_dasa_no2nd
     [,1] [,2] [,3]
[1,]  500  400  100
[2,]  250  520  230

我该如何解决这个问题？ 非常感谢！

Answer 1

先在循环外注释掉这一行：

# trans_matrix_dasa_no2nd[1,] <- c(aaa[i],bbb[i],ccc[i])

然后将其放入循环中，以便它可以适当地访问值并更新先前的状态值。您需要以不同的方式处理第一种情况，以便 i-1 上的索引不会引用 0：

 cycle_patient_distribution_dasa_no2nd[1,] <-initial_patient_distribtion
 for (i in 2:n){ 
        trans_matrix_dasa_no2nd[1,] <- c(aaa[i],bbb[i],ccc[i])
        cycle_patient_distribution_dasa_no2nd[i,] <-
               cycle_patient_distribution_dasa_no2nd[i-1,]%*%    
                                      (trans_matrix_dasa_no2nd)
                 }

Answer 2

问题已解决：循环前插入矩阵

initial_patient_distribution <- c (1000,0,0)
aaa <- c(1,0.7,0.6,0.5,0.4)
bbb <- c(1, 0.2,0.3, 0.4, 0.5)
ccc <- c(1, 0.1,0.1,0.1,0.1)

cycle_patient_distribution_dasa_no2nd <- matrix (,nrow=5,ncol=3)

cycle_patient_distribution_dasa_no2nd[1,] <-initial_patient_distribution
 for (i in 2:5){
 trans_matrix_dasa_no2nd <- matrix (,nrow=3,ncol=3)
 trans_matrix_dasa_no2nd[1,] <- c(aaa[i],bbb[i],ccc[i])
 trans_matrix_dasa_no2nd[2,] <- c(0,0.5,0.5)
 trans_matrix_dasa_no2nd[3,] <- c(0,0,1)
cycle_patient_distribution_dasa_no2nd[i,] <- cycle_patient_distribution_dasa_no2nd[i-    1,]%*%(trans_matrix_dasa_no2nd)
}

cycle_patient_distribution_dasa_no2nd
     [,1]  [,2]  [,3]
[1,] 1000   0.0   0.0
[2,]  700 200.0 100.0
[3,]  420 310.0 270.0
[4,]  210 323.0 467.0
[5,]   84 266.5 649.5