I'm new to bugs, so please bear with me. Can someone tell me if the
following two models are doing the same thing? The reason I ask is
that with the same data, the first (based on 4 separate coeffs
a1--a4) takes about 50 secs, while the second (based on a vectorized
form, a[]) takes about 300. The means are about the same, though
R-hat's in the second version are quite a bit better.
(Also, and completely unrelated: is there any way to get more than 2
decimal places in the display of the means?)
Thanks!!
Here are the two models: (these are censored regressions, the first
is essentially a copy of code in Gelman+Hill):
===================== model 1 : 4 separate a's
model{
for (i in 1:n){
z.lo[i]<- C * equals(y[i],C)
z[i]~dnorm(z.hat[i],tau.y)I(z.lo[i],)
z.hat[i]<-a1*x[i,1]+a2*x[i,2]+a3*x[i,3]+a4*x[i,4]
}
a1~dunif(0,100)
a2~dunif(0,100)
a3~dunif(0,100)
a4~dunif(0,100)
tau.y<-pow(sigma.y,-2)
sigma.y~dunif(0,100)
}
============== model 2 : vector of a's
model{
for (i in 1:n){
z.lo[i]<- C * equals(y[i],C)
z[i]~dnorm(z.hat[i],tau.y)I(z.lo[i],)
z.hat[i]<-inprod(a[],x[i,])
}
for (j in 1:k){
a[j]~dunif(0,100)
}
tau.y<-pow(sigma.y,-2)
sigma.y~dunif(0,100)
}
and here, for reference, is the R calling code:
x<-as.matrix(iv)
y<-dv
C<-cens
z<-ifelse(y==C,NA,y)
n<-length(z)
data1<-list(x=x,y=y,z=z,n=n,C=C)
inits1<-function(){
list(a1=runif(1),a2=runif(1),a3=runif(1),a4=runif(1),sigma.y=runif(1))}
params1<-c("a1","a2","a3","a4","sigma.y")
## now the bugs call for model 1
proc.time()
aasho.1<-bugs(data1,inits1,params1,"aasho1.bug",n.iter=10000,debug=FALSE)
proc.time()
print(aasho.1,digits=4)
now we try a vector approach
k<-4 # niv
data2<-list(x=x,y=y,z=z,n=n,C=C,k=k)
inits2<-function(){
list(a=runif(k),sigma.y=runif(1))}
params2<-c("a","sigma.y")
## now the bugs call for model 2
proc.time()
aasho.2<-bugs(data2,inits2,params2,"aasho2.bug",n.iter=10000,debug=FALSE)
proc.time()
print(aasho.2,digits=6)
------------------------
Philip A. Viton
City Planning, Ohio State University
275 West Woodruff Avenue, Columbus OH 43210
viton.1 at osu.edu
Philip A. Viton wrote:> > > I'm new to bugs, so please bear with me. Can someone tell me if the > following two models are doing the same thing? The reason I ask is > that with the same data, the first (based on 4 separate coeffs > a1--a4) takes about 50 secs, while the second (based on a vectorized > form, a[]) takes about 300. The means are about the same, though > R-hat's in the second version are quite a bit better. >The model definitions do seem to be formally equivalent, but BUGS is a pretty big-ugly-black-box thing, and your timing (and Rhats) seem to be pretty good evidence that WinBUGS is *not* constructing the same sampler in both cases. (Do the posterior distributions agree once you get to convergence?) You could try the same model in JAGS (the R2jags package lets you use the same R code to call it, just substituting "jags()" for "bugs()" if you wanted to pursue this further. In my limited experience JAGS is usually slower than WinBUGS but it's a good cross-check. Philip A. Viton wrote:> > (Also, and completely unrelated: is there any way to get more than 2 > decimal places in the display of the means?) > > [snip snip snip] >print(bugs.object,digits=5) Ben Bolker -- View this message in context: http://www.nabble.com/R2winbugs-%3A-vectorization-tp21019193p21021187.html Sent from the R help mailing list archive at Nabble.com.
I remember having similar problem with inprod function. As far as I could remember a sole deference in my models was that I used inprod instead of explicit sum (exactly as you did). In my case the inprod version was faster but result were completely aberrant. So I abandoned the inprod as unreliable. I did use OpenBugs (it's newer version of WinBugs) and BRugs interface from R. On Mon, 15 Dec 2008 18:23:39 +0100, Philip A. Viton <viton.1 at osu.edu> wrote:> > I'm new to bugs, so please bear with me. Can someone tell me if the > following two models are doing the same thing? The reason I ask is that > with the same data, the first (based on 4 separate coeffs a1--a4) takes > about 50 secs, while the second (based on a vectorized form, a[]) takes > about 300. The means are about the same, though R-hat's in the second > version are quite a bit better. > > > (Also, and completely unrelated: is there any way to get more than 2 > decimal places in the display of the means?) > > > Thanks!! > > > > Here are the two models: (these are censored regressions, the first is > essentially a copy of code in Gelman+Hill): > > ===================== model 1 : 4 separate a's > model{ > for (i in 1:n){ > z.lo[i]<- C * equals(y[i],C) > z[i]~dnorm(z.hat[i],tau.y)I(z.lo[i],) > z.hat[i]<-a1*x[i,1]+a2*x[i,2]+a3*x[i,3]+a4*x[i,4] > } > a1~dunif(0,100) > a2~dunif(0,100) > a3~dunif(0,100) > a4~dunif(0,100) > tau.y<-pow(sigma.y,-2) > sigma.y~dunif(0,100) > } > > > ============== model 2 : vector of a's > model{ > for (i in 1:n){ > z.lo[i]<- C * equals(y[i],C) > z[i]~dnorm(z.hat[i],tau.y)I(z.lo[i],) > z.hat[i]<-inprod(a[],x[i,]) > } > for (j in 1:k){ > a[j]~dunif(0,100) > } > tau.y<-pow(sigma.y,-2) > sigma.y~dunif(0,100) > } > > > and here, for reference, is the R calling code: > > x<-as.matrix(iv) > y<-dv > C<-cens > z<-ifelse(y==C,NA,y) > n<-length(z) > data1<-list(x=x,y=y,z=z,n=n,C=C) > inits1<-function(){ > list(a1=runif(1),a2=runif(1),a3=runif(1),a4=runif(1),sigma.y=runif(1))} > params1<-c("a1","a2","a3","a4","sigma.y") > > ## now the bugs call for model 1 > proc.time() > aasho.1<-bugs(data1,inits1,params1,"aasho1.bug",n.iter=10000,debug=FALSE) > proc.time() > print(aasho.1,digits=4) > > now we try a vector approach > k<-4 # niv > data2<-list(x=x,y=y,z=z,n=n,C=C,k=k) > inits2<-function(){ > list(a=runif(k),sigma.y=runif(1))} > params2<-c("a","sigma.y") > > ## now the bugs call for model 2 > proc.time() > aasho.2<-bugs(data2,inits2,params2,"aasho2.bug",n.iter=10000,debug=FALSE) > proc.time() > print(aasho.2,digits=6) > > ------------------------ > Philip A. Viton > City Planning, Ohio State University > 275 West Woodruff Avenue, Columbus OH 43210 > viton.1 at osu.edu > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.
Many thanks to all responders. It turns out that there is a winbugs
update available, which defines a new version of inprod, inprod2.
Using inprod2 the vectorized code runs in about the same time as the
"scalar" version.
On the print problem: the issue here turns out to be that the arm
package (or one of the packages it loads) resets digits to 2. So the
solution is to reset digits ( options("digits"=5) or whatever) after
loading your packages.
Thanks again!
------------------------
Philip A. Viton
City Planning, Ohio State University
275 West Woodruff Avenue, Columbus OH 43210
viton.1 at osu.edu