Math 534 Final

i. Gibbs-Sampling Conditional Distributions

The goal is to generate values from the posterior distribution of \((\alpha_{1},\dots,\alpha_{30},\beta_{1},...,\beta{30,\tau, \alpha,\beta})\) through Gibbs Sampling by using the conditional distributions of these parameters.

The posterior distribution, in question, is proportional to the joint distribution of \((\alpha_{1},\dots,\alpha_{30},\beta_{1},...,\beta{30,\tau, \alpha,\beta})\), from which we could derive the required conditional distributions. Now that we know how the parameters above are distributed, and we know the likelihood, we could write the joint distribution as follows,

\[ \prod_{i = 1}^{30}\prod_{j = 1}^{5} \left[f(\mathbf{y_{ij}}|\alpha_{i},\beta_{i},\sigma^{2}) \right]\prod_{i = 1}^{30} \left[ f(\alpha_{i}|\alpha) f(\beta_{i}|\beta) \right] f(\alpha)f(\beta)f(\tau) \]

Conditional of \(\tau\) :

\[ f(\tau|\mathbf{y_{ij},\alpha_{i},\beta_{i}}, \alpha, \beta) \propto (\frac{1}{\sigma^2})^{10^{-3} + 75 -1} \exp \left( -\frac{1}{\sigma^2} (10^{-3} + \sum_{i = 1}^{30} \sum_{j = 1}^{5} (y_{ij} - (\alpha_{i} + \beta_{i}(x_{j}-\bar{x}) ))^2\right) \\ \] The above is a gamma kernel which implies the following for \(\tau\),

\[ \tau \sim gamma(10^{-3} + 75, (10^{-3} + \sum_{i = 1}^{30} \sum_{j = 1}^{5} (y_{ij} - (\alpha_{i} + \beta_{i}(x_{j}-\bar{x}) ))^2) \]

Conditional for \(\beta\),

\[ f(\beta | \mathbf{y_{ij},\alpha_{i},\beta_{i}}, \alpha, \tau) \propto e^{-\frac{30\beta^2 -2\beta \sum_{i = 1}^{30}\beta_i}{2}}e^{-\frac{\beta^{2}}{2000}} \]

Simplifying and completing the square yields,

\[ f(\beta | \mathbf{y_{ij},\alpha_{i},\beta_{i}}, \alpha, \tau) \propto e^{-\frac{30000 + 1}{2 \cdot 1000}(\beta - \frac{1000\sum_{i = 1}^{30}\beta_{i}}{30000 + 1})^2} \]

This is normal kernel, which implies the following for \(\beta\),

\[ \beta \sim N\left(\frac{1000\sum_{i = 1}^{30}\beta_{i}}{30000 + 1}, \frac{1000}{(30000 + 1)}\right), \]

Conditional for \(\alpha\):

\[ f(\alpha|\mathbf{y_{ij},\alpha_{i},\beta_{i}}, \beta, \tau) \propto e^{-\frac{\left(30\alpha^2 -2\alpha \sum_{i = 1}^{30}\alpha_i)\right)}{2\cdot 100}}e^{-\frac{\alpha^{2}}{2000}} \]

Simplifying and completing the square yields the following,

\[ f(\alpha|\mathbf{y_{ij},\alpha_{i},\beta_{i}}, \beta, \tau) \propto e^{\frac{300 + 1}{2\cdot 1000}\left(\alpha - \frac{10\sum_{i = 1}^{30}\alpha_{i}}{300 + 1}\right)^2} \]

This implies the following for \(\alpha\),

\[ \alpha \sim N\left( \frac{10\sum_{i = 1}^{30}\alpha_{i}}{300 + 1},\frac{1000}{(300+1)}\right) \]

The conditional for \(\beta_{i}\), for \(i = 1, \dots,30\),

\[ f(\beta_{i}|\mathbf{y_{ij},\alpha_{i},\beta_{-i}}, \alpha,\beta, \tau) \propto e^{-\frac{(\beta^{2}\sum_{j = 1}^{5}(x_{j} - \bar{x})^{2} -2\beta_{i}\sum_{j = 1}^{5} y_{ij}(x_{j}-\bar{x}))}{2\sigma^2}}e^{-\frac{(\beta_{i}^2 - 2\beta_{i}\beta)}{2}} \]

Simplifying and completing the square yields the following,

\[ f(\beta_{i}|\mathbf{y_{ij},\alpha_{i},\beta_{-i}}, \alpha,\beta, \tau) \propto e^{\frac{\sum_{j =1}^{5}(x_{j}-\bar{x})^2}{2\sigma^2}\left(\beta_{i} - \frac{\sum_{j = 1}^{5}y_{ij}(x_{j}-\bar{x}) + \beta\sigma^2}{\sum_{j = 1}^{5}(x_{j}-\bar{x})^2 + \sigma^2}\right)^2} \]

This implies the following for \(\beta_{i}\),

\[ \beta_{i} \sim N\left( \frac{\sum_{j = 1}^{5}y_{ij}(x_{j}-\bar{x}) + \beta\sigma^2}{\sum_{j = 1}^{5}(x_{j}-\bar{x})^2 + \sigma^2}, \frac{\sigma^2}{\sum_{j = 1}^{5}(x_{j} - \bar{x})^2 + \sigma^2}\right) \]

The conditional for \(\alpha_{i}\), for \(i = 1,\dots, 30\),

\[ f(\alpha_{i}|\mathbf{y_{ij},\alpha_{-i},\beta_{i}}, \alpha,\beta, \tau) \propto e^{-\frac{(5\alpha_{i}^2 - 2\alpha_{i}\sum_{j = 1}^{5} y_{ij})}{2\sigma^2}}e^{-\frac{\alpha_{i}^{2} -2\alpha_{i}\alpha}{200}} \]

Simplifying and completing the square yields the following,

\[ f(\alpha_{i}|\mathbf{y_{ij},\alpha_{-i},\beta_{i}}, \alpha,\beta, \tau) \propto e^{-\frac{500 + \sigma^2}{2\sigma^2 \cdot 100}\left(\alpha_{i} - \frac{100\sum_{j = 1}^{5} y_{ij} + \alpha\sigma^2}{500 + \sigma^2}\right)^2} \]

This implies the following for \(\alpha_{i}\)

\[ \alpha_{i} \sim N\left(\frac{100\sum_{j = 1}^{5} y_{ij} + \alpha\sigma^2}{500 + \sigma^2}, \frac{100\sigma^2}{500 + \sigma^2}\right) \]

ii. Gibbs-Sampling Algorithm for problem 1.

0. Set t = 0
1. Set \(\alpha^{(t)}\), \(\beta^{(t)}\), and \(\tau^{(t)}\) equal to their respective means of their distributions. In this case, set \(\alpha^{(t)} =0\), \(\beta^{(t)}=0\), and \(\tau^{(t)}=1\). Furthermore, set \((\sigma^2)^{(t)}\) = \(\frac{1}{\tau^{(t)}}\)
2. Given \(\alpha^{(t)}\), \(\beta^{(t)}\), and \(\tau^{(t)}\), NOTE: \((\sigma^2)^{(t)} = \frac{1}{\tau^{(t)}}\), draw \(\alpha_{i}^{(t)}\) and \(\beta_{i}^{(t)}\) using the following conditional distributions:

\[ \alpha_{i}^{(t)} \sim N\left(\frac{100\sum_{j = 1}^{5} y_{ij} + \alpha^{(t)}(\sigma^2)^{(t)}}{500 + (\sigma^2)^{(t)}}, \frac{100(\sigma^2)^{(t)}}{500 + (\sigma^2)^{(t)}}\right) \\ \]

\[ \beta_{i}^{(t)} \sim N\left( \frac{\sum_{j = 1}^{5}y_{ij}(x_{j}-\bar{x}) + \beta^{(t)}(\sigma^2)^{(t)}}{\sum_{j = 1}^{5}(x_{j}-\bar{x})^2 + (\sigma^2)^{(t)}}, \frac{(\sigma^2)^{(t)}}{\sum_{j = 1}^{5}(x_{j} - \bar{x})^2 + (\sigma^2)^{(t)}}\right) \]

for \(i = 1,\dots, 30\) 3. Now given \(\alpha_{i}^{(t)}\) and \(\beta_{i}^{(t)}\), for \(i = 1,\dots,30\), draw \(\alpha^{(t+1)}\), \(\beta^{(t+1)}\), \(\tau^{(t+1)}\) using the following conditional distributions,

\[ \alpha^{(t+1)} \sim N\left( \frac{10\sum_{i = 1}^{30}\alpha_{i}^{(t)}}{300 + 1},\frac{1000}{(300+1)}\right) \\ \]

\[ \beta^{(t+1)} \sim N\left(\frac{1000\sum_{i = 1}^{30}\beta_{i}^{(t)}}{30000 + 1}, \frac{1000}{(30000 + 1)}\right), \\ \]

\[ \tau^{(t+1)} \sim gamma(10^{-3} + 75, (10^{-3} + \sum_{i = 1}^{30} \sum_{j = 1}^{5} (y_{ij} - (\alpha_{i}^{(t)} + \beta_{i}^{(t)}(x_{j}-\bar{x}) ))^2) \]

4. Let \((\sigma^2)^{(t + 1)} = \frac{1}{\tau^{(t+1)}}\), increment t by 1, and go to step 2.

iii.

gibbs = function(nsim,burn){
gibb = matrix(NA,nrow = nsim, ncol = 63)
colnames(gibb) = c(rep("",60), "alpha", "beta", "tau")
for(i in 1:30){
  colnames(gibb)[i] = paste("a",i,sep = "")
  colnames(gibb)[30+i] = paste("b",i,sep= "")
}
# gibb[,1:30] = alpha_i's
# gibb[,31:60] = beta_i's
# gibb[,61:63] = alpha, beta, and tau respectively 
# starting values
gibb[1,61] = 0 # alpha expected value of N(0,10000)
gibb[1,62] = 0 # beta expected value of N(0,10000)
gibb[1,63] = 1 # set tau as expected values of gamma(10^-3, 10^-3)

# set values equal to variables for computation:
sig2 = 1/gibb[1,63]; alpha = gibb[1,61]; beta = gibb[1,62]

# Calculate first set of alpha_i's and beta_i's
gibb[1,c(1:30)] = rnorm(30, (100*rowSums(y) +  sig2*alpha)/(500 + sig2),
                        sqrt((sig2*100)/(500+sig2))) # alpha_i
gibb[1,c(31:60)] = rnorm(30,(y%*%xj_xb + sig2*beta)/(sum(xj_xb^2) + sig2),
                         sqrt((sig2)/(sum(xj_xb^2)+sig2))) #beta_i
# we have all initial values now we must generate next row priors based on 
# previous row alpha_i beta_i
for(i in 2:nsim){
  gibb[i,61] = rnorm(1, 10*sum(gibb[i-1,1:30])/(300+1), sqrt(1000/(300 + 1))) #alpha
  gibb[i,62] = rnorm(1, (1000*sum(gibb[i-1,31:60]))/(30000 + 1), sqrt(1000/(30000 + 1))) #beta
  # computation to avoid longer for loop in sig_sum computation
  bxj = c()
  for(w in 1:length(xj_xb)){
    temp = xj_xb[w]*gibb[i-1,31:60]
    bxj = cbind(bxj,temp)
  }
  sig_sum = sum(rowSums(((y- gibb[i-1,1:30]) - bxj)^2))/2
  gibb[i,63] = rgamma(1, 10^{-3} + 75, 10^{-3} + sig_sum) #tau
  
  sig2 = 1/gibb[i,63]; alpha = gibb[i,61]; beta = gibb[i,62] # reassign values
  
  
  gibb[i,c(1:30)] = rnorm(30, (100*rowSums(y) +  sig2*alpha)/(500 + sig2),
                    sqrt((sig2*100)/(500+sig2))) # alpha_i
  
  gibb[i,c(31:60)] = rnorm(30,(y%*%xj_xb + sig2*beta)/(sum(xj_xb^2) + sig2),
                         sqrt((sig2)/(sum(xj_xb^2)+sig2))) #beta_i

}
gibb[-c(1:burn),]
}

iv.

nsim = 50000
burn = 10000
gs.rat = gibbs(nsim,burn)

The mixing of, \(\alpha_{1}, \beta_{1}, \alpha_{15}, \beta_{15}\), and \(\tau\), as could be seen above, look ideal. None of the traces above, show a parameter getting stuck in a particular area of the parameter space. That is none, of the parameters had spurts of concentration in one particular area of the parameter space, only to jump to a very different area of the space. Although, \(\tau\) may be very close to look as a bad mix, as it appears to have holes near the 22,000, 27,000, 44,000 iteration marks, where it appears that \(\tau\) may have gotten stuck. However, the mixing looks good overall.

v.

The following are point estimates for \(\alpha_{1}, \beta_{1}, \alpha_{15}, \beta_{15},\) and \(\tau\).

sum_df = data.frame("Mean" = apply(gs.rat[,c("a1","b1","a15","b15","tau")], 2,mean)) # point estimates
rownames(sum_df) = c("Alpha 1", "Beta 1", "Alpha 15", "Beta 15", "Tau")
kable(t(sum_df))%>%
  kable_styling(full_width = F)

	Alpha 1	Beta 1	Alpha 15	Beta 15	Tau
Mean	239.9342	6.03862	242.7326	5.256386	0.0274962

nsim = 10000
x0 = c(0,0)
mu = c(1,2)
sig = matrix(c(1,.9,.9,1),2,2)
s = c(.001,.001)
bi.mh = randomwalk(nsim,c(0,0),mu, sig,s)

i.

bi.mh$rate

## [1] 0.9974 0.9967

mh.draws = mcmc(bi.mh$rw) # mcmc object
acceptance.rate = 1 - rejectionRate(mh.draws)
acceptance.rate

##      var1      var2 
## 0.9974997 0.9967997

The acceptance rate is unreasonably high for both variables. A rule of thumb states that a good acceptance rate for M-H with few variables is around 50%, and a range between 25% - 50% preferable in general. If many points are being accepted it suggest that the algorithm is most likely stuck in a particular area, taking tiny steps and not moving away from the intial value. This is caused by the extremely small choice in \(\sigma_{1}, \sigma_{2}\), as the algorithm will only generate points in the vicinity of \(x_{0}\), given such small standard deviation.

ii.

mh.draws = mcmc(bi.mh$rw, start = 1001, end = nsim) # burn-in 1000 points

The trace for both variables, are eratic and shows no clear sign of convergence. It could be seen that the variables vary at a neighborhood around 0, which were the initial values for both variables, however, it is clear that the algorithm is getting stuck at particular points and then jumping to another point only to get stuck again. What we would like to see instead of this behavior, is a couple of initial jumps (before burning), and then a random walk about a very small interval of a point, rather than jumping from point to point. Furthermore, the densities for both variables have many bumps almost resembeling a mixture of densities, rather than the normal densities it should be representing, about mean 1 and 2 respectively. Again, this is a representation of the many different places where the algorithm got stuck and produced several points before jumping to another point.

iii.

bvn = gen(9000,2,mu,sig)

The overlay of our 9000 Metropolis-Hastings draws over actual draws from a bivariate normal makes it visually clear that the M-H draws do not resemble the bivariate normal. As expected, the tiny size of our standard deviations limit the algorithm to a small neighborhood around the initial point, \(x_{0} = (0,0)\).