## Normal distribution mean unknown, variance known document.SUBSCRIPTION_OPTIONS = { "thing": "thread", "subscribed": false, "url": "subscribe", "icon": { "css": "fa fa-envelope-o" } };

Shravan
2013-10-15
2013-10-16
• Shravan - 2013-10-15

Here is a classical example that appears in various textbooks.

Let the data x be iid. x=20.1,22.0,20.1,23.1,22.4. Variance ($\sigma^2$) is known at 0.5. Let the prior be N(m=12,v=1). We can work out posterior distribution N(m,v) using the standard formula:

derive.post<-function(n,x.bar,sigma2,m,v){
v.star<- 1/( (1/v) + n/sigma2 )
m.star<-v.star*( m/v + (n*x.bar/sigma2) )
return( list(m.star,v.star) )
}

data<-list(y=c(20.1,22.0,20.1,23.1,22.4))
## let this data come from N(mu,sigma2=.5)
## Let prior be N(12,1)

(post<-derive.post(n=length(data$y), x.bar=mean(data$y),
sigma2=0.5,m=12,v=1))

[[1]]
[1] 20.673

[[2]]
[1] 0.090909


Now, I thought I could replicate this using JAGS, but the JAGS model returns a posterior distribution with mean 12 and sd=0.9996. I'm clearly doing something wrong, but I don't see the error in my code below.

cat("
model
{
for(i in 1:5){
y[i] ~ dnorm(mu[i],1/0.5)
mu[i] <- beta0
}
##prior
beta0 ~ dnorm(12,1)
}",
file="JAGSmodels/normalexercise1.jag" )

track.variables<-c("beta0")

library(rjags,quietly=T)

normalex1.mod <- jags.model(
file = "JAGSmodels/normalexercise1.jag",
n.chains = 4,

normalex1.res <- coda.samples( normalex1.mod,
var = track.variables,
n.iter = 100000,
thin = 50 )

summary(normalex1.res)
Iterations = 50:1e+05
Thinning interval = 50
Number of chains = 4
Sample size per chain = 2000

1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:

Mean             SD       Naive SE Time-series SE
12.0023         0.9996         0.0112         0.0112

2. Quantiles for each variable:

2.5%   25%   50%   75% 97.5%
10.0  11.3  12.0  12.7  13.9


Last edit: Martyn Plummer 2013-10-15
• Martyn Plummer - 2013-10-15

Yes, you forgot to supply the data argument to jags.model() and it grabbed a different y from your global environment.

Note that other modelling functions (e.g. glm()) behave the same way if a data argument is not given, but I think I will suppress this feature in rjags 4.x by making giving the data argument a default value of NULL.

• Shravan - 2013-10-16

Oops, sorry about that. I've run into trouble before with global variables and JAGS.

• Martyn Plummer - 2013-10-16

Yes it is dangerous behaviour for JAGS to do something like this without at least warning you. I have changed the default behaviour for the development version (future rjags_4-1) so that the default data is an empty list.