1 A generic function for metropolis sampling

met_gauss <- function (iters = 10000, log_f, stepsizes = 0.5, ini_value,
              iters_imc = 1,  ...)
{
    state <- ini_value
    no_var <- length (state)
    logf <- log_f (ini_value,...)
    rej <- 0

    if (!is.finite (logf)) stop ("Initial value has 0 probability")

    one_mc <- function ()
    {
        new_state <- rnorm (no_var, state, stepsizes)
        new_logf <- log_f (new_state,...)

        if (log (runif(1)) < new_logf - logf)
        {
            state <<- new_state
            logf <<- new_logf
        }
        else rej <<- rej + 1
    }

    one_sup <-  function ()
    {
        replicate (iters_imc, one_mc())
        state
    }

    mcsample <- replicate (iters, one_sup () )
    attr (mcsample, "rej.rate") <- rej / iters_imc / iters
    mcsample
}

2 Sample from truncated normal

log_tnorm <- function (x, mu, sigma, lb, ub) 
{
    if (x > lb & x < ub) dnorm (x, mu, sigma, log = T)
    else -Inf
}

par (mfrow=c(1,3))
mcsample <-  met_gauss (log_f = log_tnorm, stepsizes = 0.1, iters = 10000, 
                         mu = 0, sigma = 2, lb = 4, ub = 5, ini_value = 4.5)
attr (mcsample, "rej.rate")
## [1] 0.1127
plot(mcsample)
acf (mcsample)
hist(mcsample)

mcsample <-  met_gauss (log_f = log_tnorm, stepsizes = 1, iters = 10000, 
                         mu = 0, sigma = 2, lb = 4, ub = 5, ini_value = 4.5)
attr (mcsample, "rej.rate")
## [1] 0.6954
plot(mcsample)
acf (mcsample)
hist(mcsample)

mcsample <-  met_gauss (log_f = log_tnorm, stepsizes = 2, iters = 10000, 
                         mu = 0, sigma = 2, lb = 4, ub = 5, ini_value = 4.5)
attr (mcsample, "rej.rate")
## [1] 0.8371
plot(mcsample)
acf (mcsample)
hist(mcsample)

3 MH sampling for normal data

3.1 Sampling function

## this function performs MH sampling for normal data
## iters --- iterations of Gibbs sampling
## x --- data
## mu0 --- prior mean for mu parameter
## sigma0 -- prior variance for mu parameter
## alpha --- degree freedom for sigma parameter
## w --- prior mean for sigma parameter

mh_norm <- function (iters, stepsizes, x, mu0, sigma0, alpha, w)
{
  logmusigma <- function (mulogsigma)
  {
    mu <- mulogsigma[1]
    logsigma <- mulogsigma[2]
    sigma <- exp (logsigma)

    return (sum(dnorm (x, mu, sqrt(sigma), log = T)) +
    dnorm (mu, mu0, sqrt (sigma0), log = T) +
    dgamma (sigma, alpha/2, alpha * w/2, log = T) +
    logsigma )
  }

  mc_mulogsigma <- met_gauss (
  iters = iters, log_f = logmusigma, stepsizes = stepsizes, ini_value = c(0,0))
  cat ("Rejection rate is ", attr (mc_mulogsigma, "rej.rate"), "\n")

  list (mu = mc_mulogsigma[1,], sd =  exp (0.5*mc_mulogsigma[2,]))
}

3.2 Test 1

3.2.1 Generate A Dataset

x <- rnorm (50, 10, 2)

3.2.2 Run preliminary mcmc to determine the stepsize

mcsamples1 <- mh_norm (10000, c(20,20)/sqrt(length(x)),
                      x, 0, 1E10, 1E-5, 1E-5)
## Rejection rate is  0.983
mcsamples2 <- mh_norm (10000, c(10,10)/sqrt(length(x)),
              x, 0, 1E10, 1E-5, 1E-5)
## Rejection rate is  0.9364
## 
mcsamples3 <- mh_norm (10000, c(4,4)/sqrt(length(x)),
                      x, 0, 1E10, 1E-5, 1E-5)
## Rejection rate is  0.7417
## 
mcsamples4 <- mh_norm (10000, c(2,2)/sqrt(length(x)),
              x, 0, 1E10, 1E-5, 1E-5)
## Rejection rate is  0.4914
par (mfcol = c(2,4))
plot (mcsamples1$mu[1:1000], main = "MC trace of mu", type = "l")
acf (mcsamples1$mu[-(1:500)])
plot (mcsamples2$mu[1:1000], main = "MC trace of mu", type = "l")
acf (mcsamples2$mu[-(1:500)])
plot (mcsamples3$mu[1:1000], main = "MC trace of mu", type = "l")
acf (mcsamples3$mu[-(1:500)])
plot (mcsamples4$mu[1:1000], main = "MC trace of mu", type = "l")
acf (mcsamples4$mu[-(1:500)])

3.2.3 Run Long Chain with Appropriate Stepsize and Make Inference

mcsamples <- mh_norm (50000, c(4,4)/sqrt(length(x)),
              x, 0, 1E10, 1E-5, 1E-5)
## Rejection rate is  0.74854
par (mfrow = c(2,2))
plot (mcsamples$mu, main = "MC trace of mu", type = "l")
plot (mcsamples$sd, main = "MC trace of sd", type = "l")
acf (mcsamples$mu[-(1:200)])
plot (mcsamples$mu[-(1:200)], mcsamples$sd[-(1:200)])

## numerical summary
summary (mcsamples$mu[-(1:200)])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   8.201   9.445   9.655   9.653   9.868  11.013
summary (mcsamples$sd[-(1:200)])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.575   2.122   2.260   2.284   2.425   3.578
quantile (mcsamples$mu[-(1:200)], probs = c(0, 0.025, 0.5, 0.975, 1))
##        0%      2.5%       50%     97.5%      100% 
##  8.200504  9.018791  9.654869 10.275374 11.013467
quantile (mcsamples$sd[-(1:200)], probs = c(0, 0.025, 0.5, 0.975, 1))
##       0%     2.5%      50%    97.5%     100% 
## 1.574785 1.887410 2.259879 2.809227 3.577934