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!
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