# Metropolis-Hastings within Gibbs sampling

Suppose we have the following classical normal linear regression model:

$$y_i = \beta_1 x_{1i} + \beta_2x_{2i} + \beta_3x_{3i} + e_i$$

where $e_{i} \sim iid.N(0, \sigma^2)$ for all $i = 1, 2, \cdots, n$ and $x_{1i} = 1$ for all $i = 1, 2, \cdots, n$.

I artificially generate data values for both $x_{2i}$ and $x_{3i}$ such that $x_{2i}$ and $x_{3i}$ $\sim$ $iid N(0,1)$ for all $i = 1, 2, \cdots, n$. Define $\boldsymbol{\beta} = (\beta_1, \beta_2, \beta_3)'$ and assume a non-informative prior of the form $p(\boldsymbol{\beta}, \sigma) \propto \frac{1}{\sigma}$.

Part 1: (Already solved here but sets the context for the next part)

Use the conditional posterior pdfs $p(\boldsymbol{\beta}|\sigma, \mathbf{y})$ and $p(\sigma|\boldsymbol{\beta},\mathbf{y})$ to run a Gibbs sampler and estimate the posterior pdf for $\displaystyle \psi = \frac{\beta_2 + \beta_3}{\sigma^2}$ with $n=20$ as the sample size and $M=1000$ as the number of draws after a burn in period of $B = 100$.

Part 2: (Attempted but unsure regarding my solution).

Insert a Metropolis Hastings sub-step into the program in part 1 at iteration $j-1$ of the Gibbs algorithm to draw from the conditional posterior of $\boldsymbol{\beta}$. Use the following candidate pdf: $$q(\boldsymbol{\beta}|\sigma^{j-1}, \mathbf{y})=(2\pi)^{-\frac{3}{2}}\left|\sigma^{(j-1)2}\boldsymbol{\Omega}^{-1}\right|^{-\frac{1}{2}}\exp\left[-\frac{1}{2\sigma^{(j-1)2}}(\boldsymbol{\beta}-\boldsymbol{b})'\boldsymbol{\Omega}(\boldsymbol{\beta}-\boldsymbol{b})\right]$$

where $\sigma^{(j-1)}$ is the draw of $\sigma$ at iteration $j-1$ of the Gibbs algorithm. In specifying $\boldsymbol{\Omega}^{-1}$, mimic the basic magnitudes of the values in $(\boldsymbol{X}'\boldsymbol{X})^{-1}$ but making the numbers different.

The following is my attempt at part 2 and coded in R.

Since the MH is applied within an outer Gibbs algorithm, the MH chain does not need to be iterated as convergence occurs immediately and only one value needs to be draw from the MH algorithm. The steps of the MH (within the Gibbs) algorithm is as follows:

1. Specify the candidate function $q$ as required in the question. I multiplied each value in $(\boldsymbol{X}'\boldsymbol{X})^{-1}$ by $0.9$ for my specification of $\boldsymbol{\Omega}^{-1}$

2. Draw a candidate draw, $\boldsymbol{\beta}^c$, from $q(.)$

3. Calculate the probability $\alpha = \min\left(\frac{p(\boldsymbol{\beta}^c|\sigma^{(j-1)},\boldsymbol{y})/q(\boldsymbol{\beta}^c|\sigma^{(j-1)},\boldsymbol{y})}{p(\boldsymbol{\beta}^{(j-1)}|\sigma^{(j-1)},\boldsymbol{y})/q(\boldsymbol{\beta}^{(j-1)}|\sigma^{(j-1)},\boldsymbol{y})},1\right)$

4. With probability $\alpha$, set $\boldsymbol{\beta}^{(j)}=\boldsymbol{\beta}^c$, otherwise set $\boldsymbol{\beta}^{(j)}=\boldsymbol{\beta}^{(j-1)}$

5. Use $\boldsymbol{\beta}^{(j)}$ to define the conditional posterior of $\sigma$ from which $\sigma^{(j)}$ is drawn and so on.

My code in R is:

####################################################################################
###########     Functions to be used in the MH algorithm      ######################
####################################################################################

## Specifying the conditional posterior pdf for beta (vector)

condbetavec = function(beta,sigma){

intcont = (2*pi)^(-3/2)*sigma^(-3)*(det(invxx))^(-1/2)

cond_pdf = intcont*exp(-(1/(2*sigma^2))*t(beta - betahat)%*%solve(invxx)%*%(beta - betahat))

return(cond_pdf)

}

## Specifying the candidate pdf

cand = function(cand_beta,cand_sigma){

intcont = (2*pi)^(-3/2)*(cand_sigma)^(-3)*(det(cand_invxx))^(-1/2)

cand_pdf = intcont*exp(-(1/(2*(cand_sigma)^2))*t(cand_beta - betahat)%*%solve(cand_invxx)%*%(cand_beta - betahat))

return(cand_pdf)

}

###########################################################
########       Generation of data = y               #######
###########################################################

true_beta = matrix(c(2,0.5,0.7),3,1)      ## True value of beta (vector) used in the data generating process (dgp)

true_sig = 0.3                            ## True value of sig used in the data generating process (dgp)

## Setting the random number seed

set.seed(123456, kind = NULL, normal.kind = NULL)

## Number of observations

nobs=20

## Generating the values for x1, x2, x3 and y. Using Gaussian distribution here for convenience.

x1 = matrix(1,nobs,1)
x2 = matrix(c(rnorm(nobs,0,1)),nobs,1)
x3 = matrix(c(rnorm(nobs,0,1)),nobs,1)
u  = matrix(c(rnorm(nobs,0,1)),nobs,1)%*%true_sig

x = cbind(x1,x2,x3)

y = x%*%true_beta + u

###########################################################
### Specification of sample statistics
###########################################################

invxx = solve(t(x)%*%x)
betahat = invxx%*%t(x)%*%y
cand_invxx = invxx*0.9  ## Specifying the var-covar matrix for the candidate pdf

df = nobs - 3
sighat = sqrt((t(y-x%*%betahat)%*%(y-x%*%betahat))/df)

###############################################################
###
###  Generation of draws of beta/sigma/psi via Gibbs+MH sampling
###
###############################################################

B = 100+1        ## Extra +1 since the Gibbs sampler coded below starts at j=2, hence burn-in period is from j=2 to j=101
M = 1000         ## Number of draws after burn in period.
repl = B+M       ## Total number of Gibbs iterations including burn in and number of draws after burn in period.

## Dimensioning the matrix of beta/sigma/psi draws ##

betav = matrix(0,repl,3)
sigv = matrix(0,repl,1)
psiv = matrix(0,repl,1)

## Generating underlying normal random draws to be used in Gibbs+MH algorithm

stnbeta = matrix(c(rnorm((repl*3),0,1)),repl,3)

## Specifying starting value for Gibbs+MH chain

sigv[1] = true_sig
betav[1,] = true_beta

## Outer Gibbs Chain

for(j in 2:repl){

## MH sub-step within the Gibbs to generate beta (vector) via candidate pdf

meanbeta = betahat
varbeta = (sigv[j-1])^2*cand_invxx

betagen = chol(varbeta)%*%stnbeta[j,] + meanbeta

if (runif(1,0,1)<min(((condbetavec(betagen,sigv[j-1]))/(cand(betagen,sigv[j-1])))/((condbetavec(betav[j-1,],sigv[j-1]))/(cand(betav[j-1,],sigv[j-1]))),1)){
betav[j,] = betagen
}else{
betav[j,] = betav[j-1,]
}

## Generation of sigma via its known inverted gamma conditional

a = t(y-x%*%betav[j,])%*%(y-x%*%betav[j,])

zvec = matrix(c(rnorm(nobs,0,1)),nobs,1)

chisq = t(zvec)%*%zvec

sigv[j] = sqrt(a/chisq)

## Draws of psi from its marginal posterior pdf

psiv[j] = (betav[j,2]+betav[j,3])/(sigv[j])^2

}

#################################################################################
## Estimation of the marginal posterior pdf of psi using kernel density smoothing
#################################################################################

## Remove burn in period of psiv which is from j=1 to j=101 to create the vector psi which conists of 1000 draws after the burn in period

psi = psiv[102:repl]

## Histogram of draws of psi

hist(psi,30,ylab="Percentage frequencies",xlab=expression(psi),
main=expression(paste("Histogram of draws of ",psi)))

## kernel density plot of the draws of psi

d=density(psi)
plot(d,main=expression(paste("Kernel density estimate of marginal of ",psi)),
xlab=expression(psi),ylab=expression(paste("p(",psi,"|y)")))


Just wondering if anyone could verify whether the basic structure of my code is correct. Thanks.

• If you are still interested in the question, check my answer to a similar question. – Xi'an Dec 1 '15 at 9:08
• i am sorry but could you please explain why "Since the MH is applied within an outer Gibbs algorithm, the MH chain does not need to be iterated as convergence occurs immediately and only one value needs to be draw from the MH algorithm"? – Z宇飛 Jan 30 at 15:23