Estimating Normal Density Parameters via Metropolis

# library(R2WinBUGS)
rm(list = ls())
verbose <- FALSE
# generate data
mu.true <- 3; sig.true <- 5
set.seed(1234L); y <- rnorm(1000L,mu.true,sig.true)
hist(y)

# initial vector setting and parameter values
nbatch <- 1e2L; B <- nbatch%/%10; B1 <- B+1
mu <- sig <- numeric(nbatch)

# initial parameter values
mu[1] <- 0; sig[1] <- 1
u.mu <- u.sig <- runif(nbatch)

# rejection counter
REJmu <- 0; REJsig <- 0

# log posterior density (up to a constant)
logpost_factory <- function(y)
function(mu,sig){
  loglike <- sum(dnorm(y,mu,sig,log=TRUE))
  return(loglike - log(sig))
}

logpost <- logpost_factory(y)

Sampling loop

kappa <- 0.5

t <- 2
for (t in 2:nbatch) {
  #  print(t)
  mut <- mu[t-1]; sigt <- sig[t-1]
  
  # uniform proposals with kappa=kappa
  mucand <- mut + runif(1,-kappa,kappa)
  if(verbose) cat(sprintf("mucand=%.2f\n", mucand))
  
  alph.mu <- logpost(mucand,sigt)-logpost(mut,sigt)
  if(verbose) {
  sprintf("old loglike=%.2f new loglike=%.2f\n",logpost(mut,sigt),logpost(mucand,sigt))
  sprintf("comparing: log(u.mu[t])=%f alpha.mu=%f\n", log(u.mu[t]),alph.mu)}
  if (log(u.mu[t]) <= alph.mu) { 
    mu[t] <- mucand;cat(sprintf("updating: mu[%d]=%.2f\n", t,mucand))
  } else { mu[t] <- mut; REJmu <- REJmu+1; cat("REJmu=", REJmu, "\n") }
  
  sigcand <- abs(sigt + runif(1,-kappa,kappa)); cat(sprintf("sigcand=%.2f\n", sigcand))
  if(verbose) cat(sprintf("old loglike=%.2f new loglike=%.2f\n",logpost(mu[t],sigt),logpost(mu[t],sigcand)))
  alph.sig <- logpost(mu[t],sigcand)-logpost(mu[t],sigt)
   if(verbose) cat(sprintf("comparing: log(u.sig[t])=%f with gain of loglike alpha.sig=%f\n", log(u.sig[t]),alph.sig))
  if (log(u.sig[t]) <= alph.sig) {
    if(verbose) cat(sprintf("updating: sig[%d]=%.2f\n", t,sigcand))
    sig[t] <- sigcand
  } else { sig[t] <- sigt; REJsig <- REJsig+1;  cat("REJsig=", REJsig, "\n") }
}
REJmu= 1 
sigcand=1.43
updating: mu[3]=0.19
sigcand=1.76
updating: mu[4]=0.21
sigcand=1.88
updating: mu[5]=0.23
sigcand=2.06
updating: mu[6]=0.70
sigcand=2.08
updating: mu[7]=0.85
sigcand=2.08
REJmu= 2 
sigcand=1.60
REJsig= 1 
REJmu= 3 
sigcand=2.47
updating: mu[10]=0.89
sigcand=2.07
REJsig= 2 
REJmu= 4 
sigcand=2.39
REJsig= 3 
updating: mu[12]=1.31
sigcand=2.42
REJsig= 4 
REJmu= 5 
sigcand=2.77
REJmu= 6 
sigcand=2.81
REJmu= 7 
sigcand=3.04
updating: mu[16]=1.63
sigcand=2.90
REJsig= 5 
REJmu= 8 
sigcand=2.85
REJsig= 6 
REJmu= 9 
sigcand=3.12
updating: mu[19]=1.75
sigcand=3.36
REJmu= 10 
sigcand=3.44
updating: mu[21]=1.85
sigcand=3.51
REJmu= 11 
sigcand=3.72
REJmu= 12 
sigcand=3.64
REJsig= 7 
REJmu= 13 
sigcand=3.39
REJsig= 8 
REJmu= 14 
sigcand=3.95
updating: mu[26]=2.23
sigcand=3.69
REJsig= 9 
updating: mu[27]=2.33
sigcand=4.27
updating: mu[28]=2.36
sigcand=4.67
REJmu= 15 
sigcand=5.00
REJmu= 16 
sigcand=4.89
REJmu= 17 
sigcand=4.64
REJsig= 10 
updating: mu[32]=2.38
sigcand=4.75
REJsig= 11 
updating: mu[33]=2.66
sigcand=4.53
REJsig= 12 
updating: mu[34]=2.78
sigcand=5.00
updating: mu[35]=2.55
sigcand=5.11
updating: mu[36]=2.76
sigcand=5.17
updating: mu[37]=2.91
sigcand=5.33
REJsig= 13 
REJmu= 18 
sigcand=4.76
updating: mu[39]=2.74
sigcand=4.28
REJsig= 14 
REJmu= 19 
sigcand=4.45
REJsig= 15 
updating: mu[41]=2.70
sigcand=4.81
updating: mu[42]=2.82
sigcand=4.54
REJsig= 16 
REJmu= 20 
sigcand=5.09
REJmu= 21 
sigcand=5.50
REJsig= 17 
updating: mu[45]=2.99
sigcand=4.87
updating: mu[46]=2.95
sigcand=4.97
REJmu= 22 
sigcand=5.35
REJsig= 18 
REJmu= 23 
sigcand=5.34
REJsig= 19 
updating: mu[49]=2.86
sigcand=4.73
REJsig= 20 
updating: mu[50]=2.81
sigcand=4.71
REJsig= 21 
updating: mu[51]=2.80
sigcand=4.75
REJsig= 22 
updating: mu[52]=2.85
sigcand=4.90
REJmu= 24 
sigcand=4.87
REJsig= 23 
updating: mu[54]=2.88
sigcand=5.18
REJsig= 24 
REJmu= 25 
sigcand=4.49
REJsig= 25 
REJmu= 26 
sigcand=5.09
REJsig= 26 
REJmu= 27 
sigcand=4.88
REJmu= 28 
sigcand=4.49
REJsig= 27 
updating: mu[59]=2.66
sigcand=4.62
REJsig= 28 
updating: mu[60]=3.06
sigcand=4.63
REJsig= 29 
updating: mu[61]=2.99
sigcand=5.11
REJsig= 30 
updating: mu[62]=2.88
sigcand=5.22
REJsig= 31 
updating: mu[63]=2.55
sigcand=4.82
updating: mu[64]=2.54
sigcand=4.75
REJsig= 32 
updating: mu[65]=2.48
sigcand=5.27
updating: mu[66]=2.51
sigcand=5.59
REJsig= 33 
updating: mu[67]=2.68
sigcand=5.37
REJsig= 34 
updating: mu[68]=3.00
sigcand=5.50
REJsig= 35 
updating: mu[69]=3.13
sigcand=5.74
REJsig= 36 
REJmu= 29 
sigcand=5.33
REJmu= 30 
sigcand=5.39
updating: mu[72]=2.63
sigcand=5.81
REJsig= 37 
REJmu= 31 
sigcand=5.58
REJsig= 38 
updating: mu[74]=2.84
sigcand=4.90
REJmu= 32 
sigcand=5.30
REJsig= 39 
REJmu= 33 
sigcand=4.89
updating: mu[77]=2.95
sigcand=4.88
updating: mu[78]=2.78
sigcand=4.81
REJsig= 40 
REJmu= 34 
sigcand=5.13
REJsig= 41 
updating: mu[80]=3.10
sigcand=4.55
REJsig= 42 
updating: mu[81]=2.80
sigcand=4.63
REJsig= 43 
updating: mu[82]=2.95
sigcand=5.17
updating: mu[83]=3.11
sigcand=5.51
REJsig= 44 
REJmu= 35 
sigcand=4.93
updating: mu[85]=3.19
sigcand=5.31
REJsig= 45 
REJmu= 36 
sigcand=4.69
REJsig= 46 
REJmu= 37 
sigcand=5.18
REJsig= 47 
updating: mu[88]=3.03
sigcand=4.89
REJmu= 38 
sigcand=5.14
REJsig= 48 
updating: mu[90]=3.09
sigcand=5.04
REJmu= 39 
sigcand=5.49
REJsig= 49 
updating: mu[92]=2.67
sigcand=5.07
updating: mu[93]=3.01
sigcand=5.57
REJsig= 50 
updating: mu[94]=2.99
sigcand=5.19
updating: mu[95]=3.03
sigcand=5.09
updating: mu[96]=2.89
sigcand=5.12
updating: mu[97]=3.24
sigcand=5.62
REJsig= 51 
REJmu= 40 
sigcand=4.83
REJsig= 52 
REJmu= 41 
sigcand=4.72
REJsig= 53 
updating: mu[100]=2.91
sigcand=5.58
REJsig= 54

sequence of sampled values and ACF plots

Path of sampled values

Kernel density plots

Monte Carlo Standard Errors

library(mcmcse, quietly = TRUE)
Warning in loadNamespace(package, lib.loc): NAs introduced by coercion
Warning in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]): NAs introduced by coercion

Warning in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]): NAs introduced by coercion

Warning in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]): NAs introduced by coercion
Warning in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck = vI[[j]]): NAs introduced by coercion

Warning in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck = vI[[j]]): NAs introduced by coercion

Warning in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck = vI[[j]]): NAs introduced by coercion
Warning in methods::cacheMetaData(ns, TRUE, ns): NAs introduced by coercion
D <- data.frame(mu[B1:nbatch],sig[B1:nbatch])
mcse.mat(D)
                    est        se
mu.B1.nbatch.  2.631311 0.2147759
sig.B1.nbatch. 4.679381 0.2984221

Acceptance Rates

ACCmu <- 1-REJmu/nbatch
ACCsig <- 1-REJsig/nbatch
cat("Acceptance Rate mu =",ACCmu,"\n")
Acceptance Rate mu = 0.59 
cat("Acceptance Rate sigma = ",ACCsig, "\n")
Acceptance Rate sigma =  0.46

kernel countour plots

Using MASS

Warning in loadNamespace(package, lib.loc): NAs introduced by coercion

Estimates of Effective Sample Sizes

Sample size adjusted for autocorrelation using coda

library(coda)
Warning in loadNamespace(package, lib.loc): NAs introduced by coercion
Warning in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck = vI[[j]]): NAs introduced by coercion
?coda::effectiveSize
starting httpd help server ... done
coda::effectiveSize(mu[B1:nbatch])
    var1 
5.121767 
effectiveSize(sig[B1:nbatch])
    var1 
3.433514

Univariate effective sample size (ESS) as described in Gong and Flgal (2015)

library(mcmcse)
?mcmcse::ess
ess(D)
 mu.B1.nbatch. sig.B1.nbatch. 
      6.152171       5.723003 
multiESS(D)
[1] 7.77002

Posterior probability on hypothesis mu < 3

sum(mu[B1:nbatch] < 3)/(nbatch-B)
[1] 0.7888889