1 Gibbs Sampling Normal Data

1.1 R function for Gibbs Sampling

## this function performs Gibbs 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
gibbs_norm <- function (iters, x, mu0, sigma0, alpha, w)
{
  sumx <- sum (x)
  n <- length (x)

  ## set and initial Markov chain state
  mu <- 0
  sigma <- 1

  one_gibbs <- function ()
  {
    ## update mu
    post_var_mu <- 1 / (n/sigma + 1/sigma0)
    mu <<- rnorm (1, (sumx / sigma + mu0 / sigma0) * post_var_mu,
                 sqrt (post_var_mu))
    sigma <<- 1/rgamma (1, (alpha + n)/2, (alpha * w + sum ((x-mu)^2))/2 )
    c(mu, sigma)
  }
  mc_musigma <- replicate (iters, one_gibbs ())
  list (mu = mc_musigma[1,], sd = sqrt (mc_musigma[2,]) )
}

1.2 Test 1

x <- rnorm (50)

mcsamples <- gibbs_norm (10000, x, 0, 1E10, 1E-5, 1E-5)
plot (mcsamples$mu, main = "MC trace of mu", type = "l")

plot (density (mcsamples$mu))

acf (mcsamples$mu)

quantile (mcsamples$mu, probs = c(0, 0.025, 0.5, 0.975, 1))
##          0%        2.5%         50%       97.5%        100% 
## -0.52059278 -0.23050144  0.08057914  0.38242030  0.64444340
mean (mcsamples$mu)
## [1] 0.07852107
sd (mcsamples$mu)
## [1] 0.1563053
plot (mcsamples$sd, main = "MC trace of sd", type = "l")

plot (density ((mcsamples$sd)))

quantile (mcsamples$sd, probs = c(0, 0.025, 0.5, 0.975, 1))
##        0%      2.5%       50%     97.5%      100% 
## 0.7572009 0.9101382 1.0982354 1.3572587 1.6285138
plot (mcsamples$mu, mcsamples$sd)

1.3 Test 2

x <- rnorm (50, 10, 2)

mcsamples <- gibbs_norm (10000, x, 0, 1E10, 1E-5, 1E-5)
plot (mcsamples$mu, main = "MC trace of mu", type = "l")

plot (density (mcsamples$mu))

acf (mcsamples$mu)

quantile (mcsamples$mu, probs = c(0, 0.025, 0.5, 0.975, 1))
##       0%     2.5%      50%    97.5%     100% 
##  9.57146 10.09034 10.63268 11.16682 11.97448
mean (mcsamples$mu)
## [1] 10.63223
sd (mcsamples$mu)
## [1] 0.2734884
plot (mcsamples$sd, main = "MC trace of sd", type = "l")

plot (density ((mcsamples$sd)))

quantile (mcsamples$sd, probs = c(0, 0.025, 0.5, 0.975, 1))
##       0%     2.5%      50%    97.5%     100% 
## 1.360298 1.597500 1.922926 2.385870 2.823925
plot (mcsamples$mu, mcsamples$sd)

1.4 Test 3

x <- rnorm (5000, 10, 2)


mcsamples <- gibbs_norm (10000, x, 0, 1E10, 1E-5, 1E-5)
plot (mcsamples$mu, main = "MC trace of mu", type = "l")

plot (density (mcsamples$mu))

acf (mcsamples$mu)

quantile (mcsamples$mu, probs = c(0, 0.025, 0.5, 0.975, 1))
##        0%      2.5%       50%     97.5%      100% 
##  9.861831  9.910068  9.965338 10.021130 10.073268
mean (mcsamples$mu)
## [1] 9.965285
sd (mcsamples$mu)
## [1] 0.02848199
plot (mcsamples$sd, main = "MC trace of sd", type = "l")

plot (density ((mcsamples$sd)))

quantile (mcsamples$sd, probs = c(0, 0.025, 0.5, 0.975, 1))
##       0%     2.5%      50%    97.5%     100% 
## 1.948127 1.975361 2.014816 2.054728 2.094786
plot (mcsamples$mu, mcsamples$sd)

2 Gibbs Sampling for Censored Normal Data

2.1 R function for Gibbs Sampling

library (truncnorm)
## Warning: package 'truncnorm' was built under R version 3.4.3
## this function performs Gibbs sampling for normal data
## iters --- iterations of Gibbs sampling
## y --- data, the observed lower limits 
## l the length of interval
## mu0 --- prior mean for mu parameter
## sigma0 -- prior variance for mu parameter
## alpha --- degree freedom for sigma parameter
## w --- prior mean for sigma parameter
gibbs_norm_interval <- function (iters, y, mu0, sigma0, alpha, w, l = 1)
{
    n <- length (y)
    
    ## set and initial Markov chain state
    mu <- 0
    sigma <- 1
    x <- rep (0, n)
    one_gibbs <- function ()
    {
        for (i in 1:n)
            x[i] <<- rtruncnorm (1, a = y[i], b = y[i] + l, 
                                 mean = mu, sd = sqrt(sigma))
        sumx <- sum (x)
        ## update mu
        post_var_mu <- 1 / (n/sigma + 1/sigma0)
        post_mean_mu <- (sumx / sigma + mu0 / sigma0) * post_var_mu
        mu <<- rnorm (1,post_mean_mu, sqrt (post_var_mu))
        sigma <<- 1/rgamma (1, (alpha + n)/2, (alpha * w + sum ((x-mu)^2))/2 )
        
        
        return(c(mu, sigma, x))
    }
    mc_musigma <- replicate (iters, one_gibbs ())
    
    list (mu = mc_musigma[1,], sd = sqrt (mc_musigma[2,]), 
          x = mc_musigma[-(1:2),] )
}

2.2 A Test with Normal Data Observed at Tens

x <- rnorm (200, 100, 40); 
y <- floor (x/10)*10; 

head(data.frame (x,y))
##           x   y
## 1  45.54441  40
## 2  80.54547  80
## 3  32.42080  30
## 4  32.21573  30
## 5 147.15844 140
## 6  92.72340  90
mean (y); sd (y)
## [1] 94.85
## [1] 43.83589
mcsamples <- gibbs_norm_interval (10000, y, 0, 1E10, 1E-5, 1E-5, 10)

plot (mcsamples$mu, main = "MC trace of mu", type = "l")

plot (density (mcsamples$mu))

acf (mcsamples$mu)

quantile (mcsamples$mu, probs = c(0, 0.025, 0.5, 0.975, 1))
##        0%      2.5%       50%     97.5%      100% 
##  88.09703  93.76378  99.86840 105.97811 113.10170
mean (mcsamples$mu)
## [1] 99.89053
sd (mcsamples$mu)
## [1] 3.125735
plot (mcsamples$sd, main = "MC trace of sd", type = "l")

plot (density ((mcsamples$sd)))

quantile (mcsamples$sd, probs = c(0, 0.025, 0.5, 0.975, 1))
##       0%     2.5%      50%    97.5%     100% 
## 36.68927 39.81848 43.82814 48.30569 52.92712
plot (mcsamples$mu, mcsamples$sd)

plot (mcsamples$mu[1:100], mcsamples$sd[1:100], type = "b")

plot (mcsamples$x[1, 1:100], type = "l")

hist (mcsamples$x[1,])