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
vito...@osu.edu
______________________________________________
R-help@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.
