Is there a way to obtain and store positions of matrix element in JAGS?

Question

I am developing a bayesian hierarchical model in R with BUGS code in JAGS.

In my model, I have two matrices that contain relevant information about each another in the same exact matrix position. My information is structured by rows. I apply a mathematical operation to the first matrix, Distmat, by row:

diffmat[i,j] <- abs(Distmat[birthterr[i],j] - Dist[i])

I am interested to record the column position of every minimum value in each row of diffmat in a new vector, to then apply this vector to the second matrix. This would be relatively easy in regular R code using functions "which" or "which.min":

a <- numeric()
for (i in 1:dim(diffmat)[1])
  for (j in 1:dim(diffmat)[2])
a[i] <- which.min(diffmat[i,])

And then apply vector "a" to the second matrix (terrmat) to obtain the values associated with Distmat positions:

b <- numeric(0)
for (i in 1:dim(diffmat)[1])
  for (j in 1:dim(diffmat)[2])
b[i] <- terrmat[i, a[i]]

However, apparently BUGS code does not recognize either which or which.min(), and I am struggling to find a way to store these matrix row positions in vectors. Perhaps there is a very simple solution to this, but I really got stuck there. Hope my info was enough clear.

Any suggestions would be very appreciated. Thanks for your time!

score 0 · Accepted Answer · answered Mar 16 '21 at 18:02

Here's a minimal working example. The analogs here are that x would be your diffmat. I'm drawing it at random here, but it should still work if you're defining it otherwise. Essentially, you rank the values of x in each row and make a new matrix e that is a dummy matrix that is coded 1 if x[i,j] is ranked 1 and 0 otherwise. Then you take the inner product of that and a vector of values from 1:ncol(terrmat) assuming terrmat and diffmat are of the same dimensions. Then that gives you the column index of the first ranked value for observation i. The ymat in the example below is where your terrmat would go. I think it will be pretty slow on any real-sized problem, but it appears to work from the output below.

dl <- list(
ymat = matrix(1:3, ncol=3, nrow=5, byrow=TRUE),
yinds = 1:3
)
mods <- "model{
  for(i in 1:5){
    for(j in 1:3){
      x[i,j] ~ dnorm(0,1)
      e[i,j] <- equals(rx[i,j], 1)
    }
    rx[i,1:3] <- rank(x[i,1:3])
    ind[i] <- inprod(e[i,], yinds)
    yval[i] <- ymat[i,ind[i]]
  }
  
}"

library(runjags)

out <- run.jags(mods, data=dl, monitor="yval")

out
# 
# JAGS model summary statistics from 20000 samples (chains = 2; adapt+burnin = 5000):
#   
#         Lower95 Median Upper95   Mean      SD Mode     MCerr MC%ofSD SSeff       AC.10    psrf
# yval[1]       1      2       3 1.9973  0.8139    2 0.0058421     0.7 19409  -0.0077146 0.99996
# yval[2]       1      2       3 2.0067 0.81605    3 0.0057704     0.7 20000  0.00049096  1.0003
# yval[3]       1      2       3 1.9895  0.8142    2 0.0057573     0.7 20000  0.00066309       1
# yval[4]       1      2       3 1.9973 0.81638    1 0.0057727     0.7 20000 -0.00040016 0.99998
# yval[5]       1      2       3  1.993 0.81611    1 0.0057708     0.7 20000  -0.0027988 0.99996
# 
# Total time taken: 0.7 seconds

Hi Dave, thanks for this suggestion. With this code, my model now runs smoothly and converges properly. As you said, computational time increases quite a lot, but at the moment I cannot think of a more efficient alternative. Thanks again! — JAB, Mar 26 '21 at 10:03

Is there a way to obtain and store positions of matrix element in JAGS?

1 Answers1