A4_MCMC.Rmd

---
title: "STAT 812 - Assignment 4"
author: "Kangjie Zhang"
date: "December 9, 2018"
output:
  md_document:
    variant: markdown_github
    df_print: paged
    toc: yes   
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Question 1
### Joint distribution
Joint distribution of $\mu,\ \sigma^2$ and latent variances $v_1,...,v_n$ given data $y_1,...,y_n$ is proportional to
$$
(\prod_{i=1}^nP(y_i|\mu,v_i)P(v_i|\sigma^2))P(\mu)P(\sigma^2)
$$

$$
=(\prod_{i=1}^n \frac{1}{\sqrt{ 2\pi v_i}}e^{-\frac{(y_i-\mu)^2}{2v_i}} \times 
\frac{1}{\Gamma(\frac{\nu}{2})}v_i^{-(\nu/2+1)}(\frac{\nu \sigma^2}{2})^{\nu/2} e^{-\frac{\nu \sigma^2}{2v_i}}) \times
\frac{1}{\sqrt{2\pi100}}e^{-\frac{\mu^2}{200}} \times
\frac{1}{\sigma^2}
$$

$$
\propto (\prod_{i=1}^n v_i)^{-\frac{\nu}{2}- \frac{3}{2}} 
(\sigma^2)^{\frac{\nu n}{2}-1}
e^{-\sum_{i=1}^{n}\frac{(y_i-\mu)^2}{2v_i}} 
e^{-\sum_{i=1}^n\frac{\nu \sigma^2}{2v_i}} 
e^{-\frac{\mu^2}{200}}
$$

### Conditional distribution of $\mu$
$$
P(\mu|\mathbf{v},\mathbf{y},\sigma^2) \propto e^{-\sum_{i=1}^{n}\frac{(y_i-\mu)^2}{2v_i}} e^{-\frac{\mu^2}{200}} \\
\propto \exp(-\frac{1}{2}[(\sum_{i=1}^n \frac{1}{v_i}+\frac{1}{100})\mu^2-2\mu \sum_{i=1}^n\frac{y_i}{v_i}])
$$

$$
\frac{1}{\sigma_1^2}=(\sum_{i=1}^n \frac{1}{v_i}+\frac{1}{100}) \Rightarrow \ \sigma_1^2=1/(\sum_{i=1}^n \frac{1}{v_i}+\frac{1}{100})
$$

$$
\frac{\mu_1}{\sigma_1^2}=\sum_{i=1}^n \frac{y_i}{v_i} \Rightarrow \ \mu_1=\sigma_1^2 \sum_{i=1}^n \frac{y_i}{v_i}
$$
$$
\mu|\mathbf{v},\mathbf{y},\sigma^2 \sim N(\mu_1,\sigma_1^2) 
$$

### Conditional distribution of $\sigma^2$
$$
P(\sigma^2|\mathbf{v},\mathbf{y},\mu) \propto  (\sigma^2)^{\frac{\nu n}{2}-1}
e^{-\sum_{i=1}^n\frac{\nu \sigma^2}{2v_i}} 
$$

$$
\sigma^2|\mathbf{v},\mathbf{y},\mu \sim Gamma(\frac{\nu n}{2},\ \sum_{i=1}^n\frac{\nu}{2v_i})
$$

### Conditional distribution of v
$$
P(\mathbf{v}|\sigma^2,\mathbf{y},\mu) \propto  
(\prod_{i=1}^n v_i)^{-\frac{\nu}{2}- \frac{3}{2}} 
e^{-\sum_{i=1}^{n}\frac{(y_i-\mu)^2}{2v_i}} 
e^{-\sum_{i=1}^n\frac{\nu \sigma^2}{2v_i}}  \\
\propto \prod_{i=1}^n[v_i^{-\frac{\nu}{2}-\frac{3}{2}}\exp(\frac{-(y_i-\mu)^2-\nu \sigma^2}{2v_i}) ]
$$
Above is the joint distribution of $\mathbf{v}$, so $v_i$ follows an inverse gamma distribution as following:
$$
v_i|\sigma^2,y_i,\mu \sim IG(\frac{\nu +1}{2},\ \frac{(y_i-\mu)^2+\nu \sigma^2}{2})
$$

## R function - Gibbs sampling
```{r}
## Gibbs sampling
## no_gibbs --- iterations of Gibbs sampling
## sigma.2 --- sigma^2
## y --- data, generated from function data_y
gibbs <- function(no_gibbs,y,nu)
{
  n <- length (y)
  ## set and initial Markov chain state
  mu <- 0
  sigma.2 <- 1
  v <- rep (1, n)

  one_gibbs <- function ()
  {
    for (i in 1:n){
      v[i] <<- 1/rgamma(1,(nu+1)/2,((y[i]-mu)^2+nu*sigma.2)/2)
    }
    ## update mu and sigma.2
    post_var_mu <- 1 / (sum(1/v) + 1/100)
    post_mean_mu <- sum(y/v) * post_var_mu    
    mu <<- rnorm (1,post_mean_mu, sqrt (post_var_mu))
    sigma.2 <<- rgamma (1, (nu* n)/2, sum (nu/(2*v)) )
    return(c(mu, sigma.2, v))
  }
  mc_musigma <- replicate (no_gibbs, one_gibbs ())

  list (mu = mc_musigma[1,], sigma.2 = mc_musigma[2,], 
        v = mc_musigma[-(1:2),] )
}

## data generation
data_y<-function(n,nu,mu,sigma.2){
  v<-rep(0,n)
  y<-rep(0,n)
  for (i in 1:n) {
    v[i] <- 1/rgamma(1,nu/2,nu*sigma.2/2)
    y[i] <- rnorm(1,mu,sqrt(v[i]))
  }
  y
}

```

## Test 1 ($\nu=2,\mu=-3,\sigma^2=1$) with different number of n 
```{r}
# n=100, no_gibbs=1000
y<-data_y(n=100,nu=2,mu=-3,sigma.2 = 1)
samples<-gibbs(1000,y,nu=2)
## posterior mean, median, 95% credible intervals of mu and sigma.2
mean (samples$mu)
quantile (samples$mu, probs = c(0.025, 0.5, 0.975))
mean (samples$sigma.2)
quantile (samples$sigma.2, probs = c(0.025, 0.5, 0.975))
## plots
par (mfrow=c(2,3))
plot (samples$mu, main = "MC trace of mu", type = "l");abline(h=-3)
plot (density (samples$mu),,main = "density of mu")
acf (samples$mu)
plot (samples$sigma.2, main = "MC trace of sigma.2", type = "l");abline(h=1)
plot (density (samples$sigma.2),main = "density of sigma.2")
acf (samples$sigma.2)

# n=1000, no_gibbs=1000
y<-data_y(n=1000,nu=2,mu=-3,sigma.2 = 1)
samples<-gibbs(1000,y,nu=2)
## posterior mean, median, 95% credible intervals of mu and sigma.2
mean (samples$mu)
quantile (samples$mu, probs = c(0.025, 0.5, 0.975))
mean (samples$sigma.2)
quantile (samples$sigma.2, probs = c(0.025, 0.5, 0.975))
## plots
par (mfrow=c(2,3))
plot (samples$mu, main = "MC trace of mu", type = "l");abline(h=-3)
plot (density (samples$mu),,main = "density of mu")
acf (samples$mu)
plot (samples$sigma.2, main = "MC trace of sigma.2", type = "l");abline(h=1)
plot (density (samples$sigma.2),main = "density of sigma.2")
acf (samples$sigma.2)

# n=10000, no_gibbs=1000
y<-data_y(n=10000,nu=2,mu=-3,sigma.2 = 1)
samples<-gibbs(1000,y,nu=2)
## posterior mean, median, 95% credible intervals of mu and sigma.2
mean (samples$mu)
quantile (samples$mu, probs = c(0.025, 0.5, 0.975))
mean (samples$sigma.2)
quantile (samples$sigma.2, probs = c(0.025, 0.5, 0.975))
## plots
par (mfrow=c(2,3))
plot (samples$mu, main = "MC trace of mu", type = "l");abline(h=-3)
plot (density (samples$mu),,main = "density of mu")
acf (samples$mu)
plot (samples$sigma.2, main = "MC trace of sigma.2", type = "l");abline(h=1)
plot (density (samples$sigma.2),main = "density of sigma.2")
acf (samples$sigma.2)
```

## Test 2 ($\nu=2,\mu=5,\sigma^2=2$) with different number of n 
```{r}
# n=100, no_gibbs=1000
y<-data_y(n=100,nu=2,mu=5,sigma.2 = 2)
samples<-gibbs(1000,y,nu=2)
## posterior mean, median, 95% credible intervals of mu and sigma.2
mean (samples$mu)
quantile (samples$mu, probs = c(0.025, 0.5, 0.975))
mean (samples$sigma.2)
quantile (samples$sigma.2, probs = c(0.025, 0.5, 0.975))
## plots
par (mfrow=c(2,3))
plot (samples$mu, main = "MC trace of mu", type = "l");abline(h=5)
plot (density (samples$mu),,main = "density of mu")
acf (samples$mu)
plot (samples$sigma.2, main = "MC trace of sigma.2", type = "l");abline(h=2)
plot (density (samples$sigma.2),main = "density of sigma.2")
acf (samples$sigma.2)

# n=1000, no_gibbs=1000
y<-data_y(n=1000,nu=2,mu=5,sigma.2 = 2)
samples<-gibbs(1000,y,nu=2)
## posterior mean, median, 95% credible intervals of mu and sigma.2
mean (samples$mu)
quantile (samples$mu, probs = c(0.025, 0.5, 0.975))
mean (samples$sigma.2)
quantile (samples$sigma.2, probs = c(0.025, 0.5, 0.975))
## plots
par (mfrow=c(2,3))
plot (samples$mu, main = "MC trace of mu", type = "l");abline(h=5)
plot (density (samples$mu),,main = "density of mu")
acf (samples$mu)
plot (samples$sigma.2, main = "MC trace of sigma.2", type = "l");abline(h=2)
plot (density (samples$sigma.2),main = "density of sigma.2")
acf (samples$sigma.2)

# n=10000, no_gibbs=1000
y<-data_y(n=10000,nu=2,mu=5,sigma.2 = 2)
samples<-gibbs(1000,y,nu=2)
## posterior mean, median, 95% credible intervals of mu and sigma.2
mean (samples$mu)
quantile (samples$mu, probs = c(0.025, 0.5, 0.975))
mean (samples$sigma.2)
quantile (samples$sigma.2, probs = c(0.025, 0.5, 0.975))
## plots
par (mfrow=c(2,3))
plot (samples$mu, main = "MC trace of mu", type = "l");abline(h=5)
plot (density (samples$mu),,main = "density of mu")
acf (samples$mu)
plot (samples$sigma.2, main = "MC trace of sigma.2", type = "l");abline(h=2)
plot (density (samples$sigma.2),main = "density of sigma.2")
acf (samples$sigma.2)
```

## Compare with MLE found by nlm
```{r}
## MLE estimates found by nlm
neg_loglike <- function(b,y){
  n<-length(y)
  (2+1)/2*sum(log(1+(y-b[1])^2/(b[2]*2)))+n/2*log(b[2])
}

mle_nlm <- nlm (neg_loglike, c(3, 1.7) , y=y, hessian = T); mle_nlm
# hessian matrix 
mle_nlm$hessian
# find MLE standard errors using hessian matrix of negative log likelihood from nlm
se.nlm <- sqrt (diag (solve(mle_nlm$hessian))); se.nlm
# standard deviation of samples from gibbs sampling
sd (samples$mu);sd (samples$sigma.2)
# Compare se.nlm and standard deviation(sd) from gibbs sampling
compare <- data.frame(rbind(c(mean (samples$mu),mean (samples$sigma.2)),
                      mle_nlm$estimate,c(sd (samples$mu),sd(samples$sigma.2)),se.nlm),
                       row.names = c("gibbs","mle_nlm","gibbs_sd","se_nlm"))
colnames(compare) <- c("mu","sigma.2")
# for gibbs sampling, generate data from: y<-data_y(n=10000,nu=2,mu=5,sigma.2 = 2)
compare
```

# Question 2
## rstan 
```{r}
## rstan
# generate data from: y<-data_y(n=10000,nu=2,mu=5,sigma.2 = 2)
library (rstan)
code <- "
data {
  int<lower=1> n;
  real nu;
  vector[n] y;
  real<lower=0> alpha;
  real<lower=0> w;
}

parameters {
  vector<lower=0> [n] v;          
  real mu;
  real<lower = 0> sigma2;
}

model {
  for (i in 1:n){
    v[i]~inv_gamma(nu/2,nu*sigma2/2);
    y[i]~normal(mu,sqrt(v[i]));
  }
  mu~normal(0,10);
  sigma2~inv_gamma(alpha/2,alpha*w/2);
}
"

## either read external file .stan or read code
inits1<-list(v=rep(0.5,length(y)),mu=1,sigma2=0.1)
inits2<-list(v=rep(0.8,length(y)),mu=10,sigma2=5)
stanfit = stan(model_code=code, data=list(n=length(y),nu=2,y=y,alpha=1E-100,w=1),
                chains=2, warmup=500, iter=1000, save_warmup=T)
    
print(stanfit, digits=2,pars=c("mu","sigma2"),probs=c(0.025, 0.5, 0.975)) 
traceplot(stanfit,pars=c("mu","sigma2"), inc_warmup = T)
pairs(stanfit,pars=c("mu","sigma2"))
```

## JAGS
```{r}
## rjags
# generate data from: y<-data_y(n=10000,nu=2,mu=5,sigma.2 = 2)
library ("R2jags")
model_jags <- "
model{
  mu ~ dnorm(0,0.01)
  isigma2 ~ dgamma(alpha/2,alpha*w/2)
  sigma2 <-1/isigma2
  for (i in 1:n){
    iv[i] ~ dgamma(nu/2,nu*sigma2/2)
    v[i] <- 1/iv[i]
    y[i] ~ dnorm(mu,1/v[i])
  }
}
"
write (model_jags, file = "model.bug")
data <- list (y=y,n=length(y), nu = 2,alpha=1E-100,w=1)
# define a function to generate initial values for parameters
inits <- function ()
{
  list (mu = rnorm (1), isigma2 = runif(1,1,10))
}
# call jags to simulate MCMC
jagfit <- jags (model.file = "model.bug",
              data = data, inits = inits, 
              parameters = c("mu", "sigma2"),
              n.chains = 2, n.thin = 1, n.burnin = 500, n.iter = 1000)
jagfit
traceplot (jagfit)

# convert to MCMC object
jagfit.mcmc <- as.mcmc (jagfit) 
summary (jagfit.mcmc)
# the first chain
chain1 <- jagfit.mcmc[[1]]
# plot
plot(chain1[,2:3])

# all chains
jagfit.matrix <- as.matrix (jagfit.mcmc)
# scatterplot
plot(jagfit.matrix[,2:3])
```


## Comments for both Q1 and Q2
From test 1 and test 2, R function - gibbs got close estimates as true values. As n (sample size) increases, the estimates getting closer to true values and with smaller variances for samples. Compared with MLE estimates found by nlm with a defined objective function (used log-t here), the estimates are close but nlm has samller standard error. 

Using blackbox like jags and rstan, it will be more efficient, and the starndard deviation for posterior samples are smaller. Also, rstan got smaller standard error than jags but more time consuming.