JAGS＆R：将矩阵乘法乘积存储在3维数组中的结果

编程问答 2022-06-12

问题描述

我正在R中的JAGS中运行分层贝叶斯模型，并通过定义四个不同的矩阵来编写我的似然函数，然后使用矩阵乘法来获得每种情况的完全似然。更具体地说（如果有帮助的话），我有一个状态空间多事件标记捕获模型，每个矩阵都是一个状态转换，但是有4个转换需要考虑。我有一个运行良好的模型，但是现在我想给模型增加时间依赖性，因此我希望四个矩阵中的每个矩阵都具有第三维（时间）。我的问题是，如何在我的JAGS模型规范中对这些3-D数组进行矩阵乘法，并将所得矩阵集存储在3-D数组中。

在我的原始模型中，我这样定义每个矩阵（不是我的真实矩阵或维）：

  model {  #begin JAGS model specification

    px1[1,1:4]<-c(phiA,1-phiA)
    px1[2,1:4]<-c(0,phiA,1-phiA)
    px1[3,1-phiA)
    px1[4,1)
    
    px2[1,1:5]<-c(...)
    ##more code here
    px2[4,1:5]<-c(...)
    
    px3[1,1:4]<-c(...)
    ##more code here
    px3[5,1:4]<-c(...)
    
    px <- px1 %*% px2 %*% px3  ##define final matrix using matrix multiplication

    #more JAGS model specification code here
  } #JAGS model

一切正常，但是现在，我希望每个时间步都有一个不同的矩阵。所以我像这样设置矩阵：

model {     #begin JAGS model specification
    
    for (t in 1:yrs) {
       px1[1,1:4,t]<-c(phiA[t],1-phiA[t])
       px1[2,t]<-c(0,phiA[t],1-phiA[t])
       ##more code here to finish defining px1 matrix...
    
       px2[1,1:5,t]<-c(...)
       ###more code here to finish px2 matrix...
    
       px3[1,t]<-c(...)
       ##more code here to finish px3...
    } #t

    # more JAGS code here
  } #JAGS model

但是如何在JAGS模型规范内进行矩阵乘法，以使每个时间步长都得到最终的矩阵(px1 %*% px2 %*% px3)？我尝试了各种选择，但是无法在JAGS中使用它们。这是我尝试过的一些事情：

第一次尝试：

  model {
    #...JAGS code
    for (t in 1:yrs) {
      #[define 3-D matrices here as above]
    
      px[,t]<-px1[,t] %*% px2[,t] %*% px3[,t]
     } #t

    # more JAGS code...
  } #JAGS model

第二次尝试：

    for (t in 1:yrs) {
      #[define 3-D matrices here as above]
    
      px[1:4,t]<-px1[1:4,t] %*% px2[1:4,t] %*% px3[1:5,t]
    } #t

第三次尝试：

    for (t in 1:yrs) {
      #[define matrix 3-D px1 matrix as above]
      px1.tmp[1:4,1:4]<-px1[,t]  #store first 2 dimensions in temporary matrix
    
      #[define matrix 2 as above]
      px2.tmp[1:4,1:5]<-px2[,t]
    
      #[define matrix 3 as above]
      px3.tmp[1:5,1:4]<-px3[,t]
    
      px.tmp<-px1.tmp %*% px2.tmp %*% px3.tmp
    
      px[,t]<-px.tmp
    }

我已经尝试了这些基本主题的其他变体，包括在JAGS模型定义之外设置最终的px数组。

    model {
      #JAGS model code here...
    } #end JAGS model code

    px<-array(dim=c(4,4,Years)

根据我的方法，我会得到略有不同的错误，但是它们都围绕矩阵乘法运算。第三次尝试+在模型定义之外定义px数组时，出现以下错误

>RUNTIME ERROR:
>Compilation error on line 1021.
>Cannot evaluate subset expression for px

1021行是jags文件中的矩阵乘法行。

在此方面，我将不胜感激，谢谢！

解决方法

假设calc是您计算矩阵乘积的函数，您可以lapply将其乘以年份，然后使用abind::abind函数将结果列表项绑定到3-D数组中

library(abind)


phiA <- seq(0,1,.01)
calc <- function(t) {
    px1<-array(dim=c(4,4))
    px2<-array(dim=c(4,4))
    px3<-array(dim=c(4,4))
    px1[1,1:4]<-c(phiA[t],1-phiA[t])
    px1[2,1:4]<-c(0,phiA[t],1-phiA[t])
    px1[3,1-phiA[t])
    px1[4,1)
    
    px2[4,1-phiA[t])
    px2[1,1-phiA[t])
    px2[2,1-phiA[t])
    px2[3,1)
    
    px3[2,1-phiA[t])
    px3[1,1-phiA[t])
    px3[4,1-phiA[t])
    px3[3,1)
    px1 %*%  px2 %*%  px3
}

abind(lapply(1:length(phiA),calc),along=3)
arraysarraysjagsmatrix-multiplicationrr