## calculation of SD seems off document.SUBSCRIPTION_OPTIONS = { "thing": "thread", "subscribed": false, "url": "subscribe", "icon": { "css": "fa fa-envelope-o" } };

Help
Andrew
2013-02-21
2013-02-25
• 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 - 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 - 2013-02-25

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

Andrew