Bayesian approach systematically overestimates sigma (SD)

sigma-bias

library(dclone)
set.seed(1)

## Bayesian estimate of sigma
N <- 9
intcept <- 1
sigma <- 1
Y <- rnorm(N, intcept, sigma)

dat <- list(Y=Y, N=N)
mod <- custommodel("model {
    intcept ~ dnorm(0, 0.01)
    sigma ~ dunif(0, 10)
    for (i in 1:N) {
        Y_exp[i] <- intcept
        Y[i] ~ dnorm(Y_exp[i], 1/sigma^2)
    }
}")

m <- jags.fit(dat, c("intcept", "sigma"), mod)
c(Post=coef(m), Data.intcept=mean(Y), Data.sigma=sd(Y))


fun1 <- function() {
    Y <- rnorm(N, intcept, sigma)
    dat <- list(Y=Y, N=N)
    m <- jags.fit(dat, c("intcept", "sigma"), mod)
    c(Post=coef(m), Data.intcept=mean(Y), Data.sigma=sd(Y))
}

x1 <- replicate(100, fun1())
summary(t(x1))
#  Post.intcept      Post.sigma      Data.intcept      Data.sigma
# Min.   :0.3388   Min.   :0.7067   Min.   :0.3398   Min.   :0.5821
# 1st Qu.:0.7811   1st Qu.:1.0060   1st Qu.:0.7853   1st Qu.:0.8403
# Median :1.0148   Median :1.2109   Median :1.0199   Median :1.0069
# Mean   :0.9820   Mean   :1.2291   Mean   :0.9837   Mean   :1.0208
# 3rd Qu.:1.1751   3rd Qu.:1.4323   3rd Qu.:1.1778   3rd Qu.:1.1918
# Max.   :1.6046   Max.   :2.1949   Max.   :1.6099   Max.   :1.8278

## MLE

dc <- dc.fit(dat, c("intcept", "sigma"), mod,
    n.clones=c(1, 10, 100), multiply="N")
summary(dc)
(dct <- dctable(dc))
c(Post.intcept=dct$intcept[1,"mean"], Post.sigma=dct$sigma[1,"mean"],
    MLE=coef(dc), Data.intcept=mean(Y), Data.sigma=sd(Y))


fun2 <- function() {
    Y <- rnorm(N, intcept, sigma)
    dat <- list(Y=Y, N=N)
    dc <- dc.fit(dat, c("intcept", "sigma"), mod,
        n.clones=c(1, 100), multiply="N")
    dct <- dctable(dc)
    c(Post.intcept=dct$intcept[1,"mean"], Post.sigma=dct$sigma[1,"mean"],
          MLE=coef(dc), Data.intcept=mean(Y), Data.sigma=sd(Y))
}

x2 <- replicate(100, fun2())
summary(t(x2))
#  Post.intcept      Post.sigma      MLE.intcept       MLE.sigma       Data.intcept
# Min.   :0.2056   Min.   :0.5033   Min.   :0.2071   Min.   :0.3940   Min.   :0.2069
# 1st Qu.:0.7295   1st Qu.:0.9569   1st Qu.:0.7311   1st Qu.:0.7504   1st Qu.:0.7310
# Median :0.9612   Median :1.1910   Median :0.9694   Median :0.9364   Median :0.9695
# Mean   :0.9942   Mean   :1.1782   Mean   :0.9960   Mean   :0.9259   Mean   :0.9960
# 3rd Qu.:1.2649   3rd Qu.:1.3977   3rd Qu.:1.2649   3rd Qu.:1.1071   3rd Qu.:1.2652
# Max.   :1.7428   Max.   :1.8451   Max.   :1.7435   Max.   :1.4376   Max.   :1.7433
#   Data.sigma
# Min.   :0.4171
# 1st Qu.:0.7944
# Median :0.9912
# Mean   :0.9802
# 3rd Qu.:1.1717
# Max.   :1.5212