Comparing JAGS and WinBUGs posterior distribution

Shravan
2013-09-30
2013-09-30
  • Shravan
    Shravan
    2013-09-30

    Hi all,

    I'm trying to replicate a result I got with WinBUGS in JAGS (sorry that I write math in latex; I couldn't figure out how to get Markdown to display math correctly).

    This example is also from the BUGS book (p. 60 and 61 of the Kindle version). Suppose that we have data $y_i$ which we believe comes from a heavy-tailed t-distribution:

    y.sim<-rt(100,d=4)
    

    Assume that we don't know any of the parameters and set as priors for this data:

    • $\mu \sim N(\gamma,\omega^2)$
    • $1/\sigma^2 \sim Gamma(\alpha,\beta) $
    • $df \sim Unif(2,30)$

    Then, the joint posterior distribution of each parameter is:

    $p(\mu,\sigma^2,df\mid y) \propto \hbox{Likelihood} \times [\hbox{prior} \mu] \times [\hbox{prior} \sigma^2] \times [\hbox{prior} d]$

    Now, if we want the marginal distributions, $p(\mu\mid y), p(\sigma^2\mid y), p(df\mid y)$, we can use JAGS:

    n.samples<-100
    y.sim<-rt(n.samples,df=4)
    data<-list(n=n.samples,gamma=0,
               inv.omega.squared=(1/100)^2,
               alpha=10^(-3),beta=10^(-3),y=y.sim)
    
    cat("
    model
       {
        for(i in 1:n){ y[i] ~ dt(mu,r,d)}
        mu ~ dnorm(gamma,inv.omega.squared)
        r ~ dgamma(alpha,beta)
        d ~ dcat(p[])
        p[1]<-0
        for(i in 2:30){p[i]<-1/29}
       }",
         file="multimcmcexample2.jag" )
    
    track.variables<-c("mu","r","d")
    
    library(rjags)
    
    inits <- list (list(mu=0,r=1,d=10))
    
    system.time(
    mcmceg.mod <- jags.model( 
      file = "multimcmcexample2.jag",
                         data=data,
                         inits=inits,
                         n.chains = 1,
                          n.adapt =10000 )
                )
    summary(mcmceg.res)
    plot(mcmceg.res)
    

    The corresponding code in WinBUGS is:

    cat("model {
         for (i in 1:n){
           y[i] ~ dt (mu, r, d)}
           mu ~ dnorm(gamma,inv.omega.squared)
           r ~ dgamma(alpha,beta)
           d ~ dcat(p[])
           p[1]<-0
           for(j in 2:30){p[j]<-1/29}
       }",
         file="WinBUGsmodels/multimcmcexample2.bug" )
    
    library(arm)
    
    n.samp<-100
    y.sim<-rt(n.samp,df=4)
    
    data<-list(n=n.samp,gamma=0,inv.omega.squared=0.0001,
               alpha=0.001,beta=0.001,y=y.sim)
    
    inits <- list (list(mu=0,r=1,d=10))
    parameters <- c("mu", "r", "d")
    
    ## change this for your system:
    mybugsdir<-"C:/Users/Shravan/Desktop/winbugs14/WinBUGS14/"
    
    model.sim <- bugs (data, inits, parameters, "multiparammcmc.bug", 
                         n.chains=1, n.iter=10000,debug=TRUE,
                         bugs.directory=mybugsdir)
    

    print(model.sim)
    plot(model.sim)

    My question is: When I plot the posterior distributions of d in JAGS vs WinBUGS, I get quite different distributions. In JAGS, d is almost like a uniform distribution Unif(1,30), with a spike at x=1 or so, while in WinBUGS the distribution is somewhat like a gamma distribution (roughly, dgamma(x,shape=5) in R; see p. 61 of the BUGS book for the figure).

    Why is there such a big difference? I'm probably "mis-translating" WinBUGS to JAGS but I haven't yet managed to figure out what the mistake is.

    Many thanks!

     
  • Martyn Plummer
    Martyn Plummer
    2013-09-30

    To get the same results you need to use the same data. In your JAGS model, you have inv.omega.squared=(1/100)^2=1e-4 whereas in your WinBUGS model you have inv.omega.squared=0.001=1e-3. You also need to set the random seed in R before simulating the data to ensure that your outcome variables are the same:

    set.seed(1234)
    y.sim <- rt(n.samp, df=4)
    

    In addition, make sure you are using the right model file. In your WinBUGS model, you are writing to the file WinBUGSmodels/multimcmcmexample2.bug but reading from the file multiparammcmc.bug.

    With these changes, I get comparable results from the two models. Using OpenBUGS, the mixing of d appears to be worse for WinBUGS than for JAGS, so a relatively longer run is required for WinBUGS.

     
  • Shravan
    Shravan
    2013-09-30

    Ah, that was pretty careless of me. Thanks, Martyn!

    Regarding the data being different, I had just thought that the results would be roughly comparable even though the data are not the same.

     
  • Martyn Plummer
    Martyn Plummer
    2013-09-30

    Normally yes. But when simulating from long-tailed distributions you can get some highly influential observations.

     
  • Shravan
    Shravan
    2013-09-30

    Just a small comment: the original JAGS model I posted also produces the posterior distribution for d that WinBUGS does, it's just that the plot function in JAGS uses a different binning for the histogram. When I generated the histogram myself, it has the expected shape.

    It doesn't seem to matter whether inv.omega.squared is 0.001 or 0.0001. Here is the code that works for me, in case anyone else finds it helpful (I hope that I have finally stopped making stupid mistakes :).

    n.samples<-100
    y.sim<-rt(n.samples,df=4)

    the data:

    data<-list(n=n.samples,gamma=0,
    inv.omega.squared=0.0001,
    alpha=0.001,beta=0.001,y=y.sim)

    cat("
    model
    {
    for(i in 1:n){ y[i] ~ dt(mu,r,d)}
    mu ~ dnorm(gamma,inv.omega.squared)
    r ~ dgamma(alpha,beta)
    d ~ dcat(p[])
    p[1]<-0
    for(i in 2:30){p[i]<-1/29}
    }",
    file="JAGSmodels/multimcmcexample2.jag" )

    track.variables<-c("mu","r","d")
    library(rjags)

    inits <- list (list(mu=0,r=1,d=10))

    system.time(
    mcmceg.mod <- jags.model(
    file = "JAGSmodels/multimcmcexample2.jag",
    data=data,
    inits=inits,
    n.chains = 1,
    n.adapt =10000 )
    )

    system.time(
    mcmceg.res <- coda.samples( mcmceg.mod,
    var = track.variables,
    n.iter = 10000,
    thin = 20 ) )

    summary(mcmceg.res)

    post<-jags.samples(mcmceg.mod,var=track.variables,n.iter=10000)
    op<-par(mfrow=c(1,3),pty="s")
    hist(post$d)
    hist(post$mu)
    hist(post$r)