Hello Bugs users,
I have noticed that the inprod function in Winbugs is extremely slow compared to multiplying individual elements one by one and adding them up. Has anyone else encountered this or know why this is the case?
I am trying to model random effects of different plant species on a response variable. Not all plant species are present in all samples, so I have a matrix of presence/absence (1/0) with samples as rows and species as columns. My random species effect (beta) is therefore a vector, with each item in the vector corresponding to a different species effect. It seems logical to me to program this as an inner product of the beta vector and the row from the predictor matrix for a given sample:
for (i in 1:N){ #for each sample
Y[i] ~ dnorm(mu[i], tau)
mu[i] <- beta0 + inprod(beta1[], X1[i,])
}
However this model runs extremely slowly. To make sure the model was doing what I wanted it to, I also programmed this the "brute force" way:
for (i in 1:N){
Y[i] ~ dnorm(mu[i], tau)
mu[i] <- beta0 + beta1[1]*X1[i,1] + beta1[2]*X1[i,2] + beta1[3]*X1[i,3]
+ beta1[4]*X1[i,4] + beta1[5]*X1[i,5] + beta1[6]*X1[i,6]
}
This gives the exact same results and runs 17 times faster. Looking through the archives, I haven't seen anything directly addressing this, although I have found a few posts about the slow speed of models that happen to contain the inprod function. For example:
https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=ind05&L=BUGS&P=R16625&1=BUGS&9=A&I=-3&J=on&X=53E83E639C8444807D&Y=elgersma%40eden.rutgers.edu&d=No+Match%3BMatch%3BMatches&z=4
Below I'm pasting full code for anyone who's interested. This is R code that calls Winbugs using R2WinBUGS. I've included a model using the inprod function as well as one doing the element-by-element multiplication. For this simple model, the added time due to the inprod function is not so severe, but I am also running much more complicated models, and the added time becomes more prohibitive with the complex models, and programming the element-by-element multiplication is also a lot of extra work. So any comments on this issue are much appreciated.
Thanks,
Kenneth Elgersma
library(R2WinBUGS)
#X1 is a matrix with samples as rows, species presence/absence as columns
X1 <- matrix(c(
0,0,0,1,0,1,1,0,0,1,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,1,1,1,
0,0,0,1,1,0,1,0,0,0,1,0,1,0,0,1,0,0,0,0,0,0,1,0,1,1,0,1,0,0,0,0,0,1,0,0,1,1,1,
1,0,1,1,1,1,1,0,0,0,0,0,0,0,1,0,1,0,1,1,1,0,1,0,1,1,0,1,1,1,1,1,1,1,1,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,1,0,0,0,0,
0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,1,1,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,
0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,1,1,
0,0,1,0,0,0,1,1,1,0,0,1,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,1,0,1,1,1,1,1,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,1,0,1,
1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,
0,1,0,0,1,1,0,0,1,1,0,1,0,0,0,1,0,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,
0,0,0,0,0,0,0,0,0), nrow=99, ncol=6)
#X1 is a different matrix with samples as rows, species presence/absence as columns
X2 <- matrix(c(
1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,0,1,1,1,0,1,0,0,1,1,0,1,1,1,1,1,1,1,1,1,1,
1,1,1,0,0,0,0,1,0,1,0,1,1,1,1,1,0,1,1,1,1,1,1,1,0,1,1,1,0,1,1,1,0,0,1,0,0,0,1,
1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,0,0,0,0,0,1,0,1,1,1,1,1,1,0,1,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,1,1,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0,0,0,
0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,0,0,0,
1,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,1,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,0,0,0,1,0,0,1,0,0,0,0,1,1,0,
0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0), nrow=99, ncol=6)
#Y is the response variable
Y <- c(9.8930187,-1.490307,-0.4936855,0.5393929,-1.57666,-0.728428,-1.5674914,
0.7069833,-2.8345468,-0.6530553,-0.747831,-4.0941493,-3.2709876,1.4481509,
-1.6980442,0.7288624,-4.4118578,-1.1038497,-4.238852,2.1894815,3.3573135,
-0.6171794,-2.7569383,-0.2444025,-3.5896847,-3.0592437,-2.240281,-0.9253375,
-3.0509816,0.7698062,-3.0803311,-5.8194451,-3.1037302,2.489518,-1.0438584,
-3.6437553,-2.2660199,-4.6090409,-4.4674625,-2.1653731,2.6709991,-2.9697435,
-1.3261689,-3.3756437,-3.9558114,1.2298576,8.3432455,0.8881382,-3.8493501,
1.8933769,-1.5663625,-3.0006472,-0.1109214,5.8355305,15.1267489,-1.2012711,
7.2557701,1.0296247,0.8402131,1.1968549,3.9400009,1.4664672,1.7834119,
-2.7578203,-0.7995056,0.6273709,-1.3193063,-0.6902929,-2.0584048,-0.2920117,
5.9116684,-3.3724693,0.2567529,-2.0187592,-0.3093973,1.093452,6.2728027,
16.8341751,0.3946647,-2.4653587,1.7957032,0.2632702,0.3318346,-2.1720797,
3.7329973,1.9678246,1.0105771,1.158189,0.9274523,0.3195322,-2.7850753,4.0820113,
-0.7489046,-1.5531917,-2.7992032,0.7910987,-0.7379023,1.2350259,-0.8007561)
N <- length(Y) #N is the number of samples
mymodel <- function() {
#data model
for (i in 1:N){
Y[i] ~ dnorm(mu[i], tau)
mu[i] <- beta0 + inprod(beta1[], X1[i,])
}
#priors
beta0 ~ dnorm(0,1.0E-6)
beta1[1] <- 0
for (i in 2:6){beta1[i] ~dnorm(0,tau1)}
tau ~ dgamma(0.001,0.001)
tau1 ~ dgamma(0.001,0.001)
sigma <- 1/sqrt(tau)
sigma1 <- 1/sqrt(tau1)
}
filename <- file.path(tempdir(), "model.bug")
write.model(mymodel, filename)
data <- list("X1", "N", "Y")
parameters <- c("beta0", "beta1", "sigma", "sigma1")
inits <- list(
list(beta1=c(NA,rep(0,times=5)), tau1=1, beta0=0, tau=1),
list(beta1=c(NA,rep(3,times=5)), tau1=3, beta0=3, tau=3)
)
results <- bugs(data,inits,parameters,filename,n.chains=2, n.iter=50000, n.burnin=1000,n.thin=5, debug=T)
#updates took 102 seconds
newmodel <- function() {
#data model
for (i in 1:N){
Y[i] ~ dnorm(mu[i], tau)
mu[i] <- beta0 +
beta1[1]*X1[i,1] +
beta1[2]*X1[i,2] +
beta1[3]*X1[i,3] +
beta1[4]*X1[i,4] +
beta1[5]*X1[i,5] +
beta1[6]*X1[i,6]
}
#priors
beta0 ~ dnorm(0,1.0E-6)
beta1[1] <- 0
for (i in 2:6){beta1[i] ~dnorm(0,tau1)}
tau ~ dgamma(0.001,0.001)
tau1 ~ dgamma(0.001,0.001)
sigma <- 1/sqrt(tau)
sigma1 <- 1/sqrt(tau1)
}
filename <- file.path(tempdir(), "model.bug")
write.model(newmodel, filename)
data <- list("X1", "N", "Y")
parameters <- c("beta0", "beta1", "sigma", "sigma1")
inits <- list(
list(beta1=c(NA,rep(0,times=5)), tau1=1, beta0=0, tau=1),
list(beta1=c(NA,rep(3,times=5)), tau1=3, beta0=3, tau=3)
)
results <- bugs(data,inits,parameters,filename,n.chains=2, n.iter=50000, n.burnin=1000,n.thin=5, debug=T)
#updates took 6 seconds
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|
