Normal distribution mean unknown, variance known

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

    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,
        n.adapt =2000 ,quiet=T)
    
    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

    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

    Shravan - 2013-10-16

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

     
  • Martyn Plummer

    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.

     

Log in to post a comment.

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

Sign up for the SourceForge newsletter:





No, thanks