calculation of SD seems off

Help
Andrew
2013-02-21
2013-02-25
  • Andrew

    Andrew - 2013-02-21

    Hello all
    Perhaps I'm being daft but the sd function in JAGS 3.3.0 is giving me something else to that in R 2.15.2. To illustrate this problem, please try running the following modification to the schools example. You will see that the means match but not the standard deviations. I trust that there is a very sensible explanation to this?

    Thanks!

    library(rjags)
    library(R2jags)

    //////data
    J <- 8.0
    y <- c(28.4,7.9,-2.8,6.8,-0.6,0.6,18.0,12.2)
    sd <- c(14.9,10.2,16.3,11.0,9.4,11.4,10.4,17.6)

    m1 <- function() {
    for (j in 1:J){ # J=8, the number of schools
    y[j] ~ dnorm (theta[j], tau.y[j]) # data model: the likelihood
    tau.y[j] <- pow(sd[j], -2) # tau = 1/sigma^2
    }
    for (j in 1:J){
    theta[j] ~ dnorm (mu, tau) # hierarchical model for theta
    }
    tau <- pow(sigma, -2) # tau = 1/sigma^2
    mu ~ dnorm (0.0, 1.0E-6) # noninformative prior on mu
    sigma ~ dunif (0, 1000) # noninformative prior on sigma

    m1 <- mean(theta)
    s1 <- sd(theta)
    }

    ////////this writes the model and sets up the data
    write.model(m1, con="testy.bug")
    jags.data <- list("y","sd","J")
    jags.params <- c("m1","s1","theta")
    jags.inits <- function(){
    list("mu"=rnorm(1),"sigma"=runif(1),"theta"=rnorm(J))
    }

    jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,
    n.iter=500, model.file="testy.bug")

    //////these match
    mean(jagsfit$BUGSoutput$sims.list$theta[1,])
    jagsfit$BUGSoutput$sims.list$m1[1]
    //////but not these
    sd(jagsfit$BUGSoutput$sims.list$theta[1,])
    jagsfit$BUGSoutput$sims.list$s1[1]

     
  • Martyn Plummer

    Martyn Plummer - 2013-02-25

    It's not making the 1 degree of freedom correction for estimating the mean, i.e. dividing by n instead of (n-1).

    This is now fixed in the sources and will be correct in the next release. In the mean time you can make the correction by hand. Sorry about this.

     
    • Andrew

      Andrew - 2013-02-25

      No need to apologize. Thanks for a great piece of software!

      Andrew

       

Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

JavaScript is required for this form.





No, thanks