[R] R2winbugs : vectorization

[Vitalie Spinu]

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.