Goals
This tutorial looks at one of the work horses of Bayesian estimation, the Gibbs sampler. In this tutorial, we will:
- Use the Gibbs sampler to generate bivariate normal draws.
- Create side-by-side plots of the parameter paths.
- Calculate the drawn distribution's mean and variance-covariance matrix.
Introduction
The Gibbs sampler draws iteratively from posterior conditional distributions rather than drawing directly from the joint posterior distribution. This makes the Gibbs Sampler particularly useful, as the joint posterior is not always easy to work with.
The Gibbs sampler works by restructuring the joint estimation problem as a series of smaller, easier estimation problems. For example, consider the case where the parameter vector can be broken into two blocks: $ \theta' = [\theta_{1}'\text{ }\theta_{2}'] $.
The Gibbs sampler steps
The bivariate general Gibbs Sampler can be broken down into simple steps:
- Set up sampler specifications including the number of iterations and the number of burn-ins draws.
- Choose a starting value $\:p(\theta_1|y,\:\theta_2^{(0)})$.
- Draw $\theta_2^{(r)}$ from $\:p(\theta_2|y,\:\theta_1^{(r-1)})$.
- Draw $\theta_1^{(r)}$ from $\:p(\theta_1|y,\:\theta_2^{(r)})$.
- Repeat 2-3 for desired number of iterations.
- Final calculations including means, standard deviations and plots.
Posterior distribution
In this tutorial, we consider a bivariate normal posterior distribution such that
$$ \begin{eqnarray*} \begin{pmatrix} \theta_1\\ \theta_2 \end{pmatrix} & \sim & N\left[\left(\begin{array}{c} 0\\ 0 \end{array}\right),\left(\begin{array}{ccc} 1 & \rho\\ \rho & 1 \end{array}\right)\right]\\ \end{eqnarray*} $$
where $\theta_1$ and $\theta_2$ are unknown parameters of the model, while $\rho$ is the known posterior correlation between $\theta_1$ and $\theta_2$.
Using the properties of the multivariate normal distribution
$$\theta_1|\theta_2,\: y \text{ ~ } N(\rho\theta_2,\: 1-\rho^2) \text{ ~ } \rho\theta_2 + \sqrt{1-\rho^2}N(0,\:1)$$
and
$$\theta_2|\theta_1,\: y \text{ ~ } N(\rho\theta_1,\: 1-\rho^2) \text{ ~ } \rho\theta_1 + \sqrt{1-\rho^2}N(0,\:1)$$
Set Gibbs sampler options
Our first step is to set the specifications of the Gibbs sampler such that:
- We have 10,000 total draws.
- We will make 1000 burn-in draws.
- The known $\rho$ is equal to 0.6.
/*
** Gibbs Sampler Specifications
** Known correlation parameter
*/
rho = 0.6;
// Burn-in for the Gibbs sampler
burn_in = 1000;
// Draws to keep from sampler
keep_draws = 10000; Pre-initialize vectors to hold the draws
Next, we need to initialize the storage matrices for our parameters, $\theta_1$ and $\theta_2$. The GAUSS function zeros allows us to build zero matrices with specified numbers of rows and columns.
// Initialize the variables
theta_1_do = zeros(burn_in + keep_draws, 1);
theta_2_do = zeros(burn_in + keep_draws, 1);Set initial values
The last step before starting our Gibbs sampler is to set our initial values for our parameters. In this example, we will start with $\theta_1^{(1)} = \theta_2^{(1)} = 0$. As we have already set our matrices to zero for all iterations, we can effectively set our initial values by starting our loop at r = 2.
Run Gibbs sampler using a for loop
The for loop offers a method for looping which, in general, is the preferred method for situations like the Gibbs sampler. The for loop often provides a fast, compact for iterative computations.
The for loop repeats a set of commands for a specified number of iterations based on a start, stop and step specification.
/*
** General anatomy of a 'for' loop
** 'i' is the loop counter
** Do not run this code snippet
*/
for i(start, stop, step);
// loop body work
endfor;Using this structure to run our Gibbs sampler:
// Initialize the variables for the loop to zero again
theta_1 = zeros(burn_in + keep_draws, 1);
theta_2 = zeros(burn_in + keep_draws, 1);
// Set random seed for repeatable random numbers
rndseed 45352;
// Start for loop
for i(2, burn_in + keep_draws, 1);
// Draw theta_2 given theta_1
theta_2[i, 1] = ((1 - rho^2)^0.5)*rndn(1, 1)
+ rho*theta_1[i-1, 1];
// Draw theta_1 given theta_2
theta_1[i, 1] = ((1 - rho^2)^0.5)*rndn(1, 1)
+ rho*theta_2[i, 1];
endfor;
print "The first 5 values of 'theta_1' and 'theta_2' are:";
print theta_1[1:5]~theta_2[1:5];The output reads:
The first 5 values of 'theta_1' and 'theta_2' are:
0.0000000 0.0000000
-0.65117043 -0.15106836
-0.32630327 -1.6688055
-0.98322799 -0.62352043
-0.45213979 -0.24999669
Plot the variable paths
To check the behavior of our Gibbs sampler, we can plot the path of our variables. We will first set the plot appearance using a GAUSS control structure.
// Declare 'myGibbsPlot` to be a plotControl structure
struct plotControl myGibbsPlot;
// Fill 'myGibbsPlot' with 'XY' defaults
myGibbsPLot = plotGetDefaults("xy");
// Set the text interpreter to LaTex for axis labels only
plotSetTextInterpreter(&myGibbsPlot, "latex", "axes");
// Set plot title
plotSetTitle(&myGibbsPlot, "Bivariate Normal Gibbs Sampler", "Arial", 18);
// Set axis labels
plotSetYLabel(&myGibbsPlot, "\\theta_1", "Arial", 14);We will again combine two separate plots, one for $\theta_1$ and one for $\theta_2$ using the plotLayout command. In our first cell, we will plot the pay of $\theta_1$ against the iteration number:
/*
** Change plot layout to 2 x 1 grid and place
** next plot in 1st cell
*/
plotLayout(2, 1, 1);
// Create a sequential vector counting iteration numbers
x_plot = seqa(1, 1, burn_in + keep_draws);
// Plot theta_1
plotXY(myGibbsPlot, x_plot, theta_1);Next, plot $\theta_2$ next to the graph of $\theta_1$:
// Turn off title for second plot
plotSetTitle(&myGibbsPlot, "");
// Set axis labels
plotSetYLabel(&myGibbsPlot, "\\theta_2", "Arial", 14);
/*
** Change plot layout to 2 x 1 grid and place
** next plot in 2nd cell
*/
plotLayout(2, 1, 2);
// Plot theta_1
plotXY(myGibbsPlot, x_plot, theta_2);
// Turn off plot grid for next plot
plotClearLayout();The above code should produce a plot similar to the following:
Finally, we finish our Gibbs sampler by examining the mean and variance-covariance matrix of the Gibbs sampler distribution.
// Output matrix concatenating theta_1 and theta_2
ygs = theta_1~theta_2;
ygs = ygs[burn_in+1:burn_in+keep_draws, .];
"---------------------------------------------------";
"Variance-Covariance Matrix from Gibbs Sampler";
"---------------------------------------------------";
"Sample Mean: ";; meanc(ygs);
"Sample Variance-Covariance Matrix: ";
vcx(ygs);The above code should print the following output:
---------------------------------------------------
Variance-Covariance Matrix from Gibbs Sampler
---------------------------------------------------
Sample Mean:
0.0060170250
0.0032377421
Sample Variance-Covariance Matrix:
0.97698091 0.58179239
0.58179239 0.98180547
Conclusion
Congratulations! You have:
- Used the Gibbs sampler to generate bivariate normal draws.
- Created side-by-side plots of the parameter paths.
- Calculated the drawn distribution's mean and variance-covariance matrix.
The next tutorial introduces the Metroplis Hastings sampler.
For convenience, the full program text is reproduced below.
/*
** Gibbs Sampler Specifications
** Known correlation parameter
*/
rho = 0.6;
// Burn-in for the Gibbs sampler
burn_in = 1000;
// Draws to keep from sampler
keep_draws = 10000;
// Initialize the variables for the loop to zero again
theta_1 = zeros(burn_in + keep_draws, 1);
theta_2 = zeros(burn_in + keep_draws, 1);
// Set random seed for repeatable random numbers
rndseed 45352;
// Start for loop
for i(2, burn_in + keep_draws, 1);
// Draw theta_2 given theta_1
theta_2[i, 1] = ((1 - rho^2)^0.5)*rndn(1, 1)
+ rho*theta_1[i-1, 1];
// Draw theta_1 given theta_2
theta_1[i, 1] = ((1 - rho^2)^0.5)*rndn(1, 1)
+ rho*theta_2[i, 1];
endfor;
print "The first 5 values of 'theta_1' and 'theta_2' are:";
print theta_1[1:5]~theta_2[1:5];
// Declare 'myGibbsPlot` to be a plotControl structure
struct plotControl myGibbsPlot;
// Fill 'myGibbsPlot' with 'XY' defaults
myGibbsPLot = plotGetDefaults("xy");
// Set the text interpreter to LaTex for axis labels only
plotSetTextInterpreter(&myGibbsPlot, "latex", "axes");
// Set plot title
plotSetTitle(&myGibbsPlot, "Bivariate Normal Gibbs Sampler", "Arial", 18);
// Set axis labels
plotSetYLabel(&myGibbsPlot, "\\theta_1", "Arial", 14);
/*
** Change plot layout to 2 x 1 grid and place
** next plot in 1st cell
*/
plotLayout(2, 1, 1);
// Create a sequential vector counting iteration numbers
x_plot = seqa(1, 1, burn_in + keep_draws);
// Plot theta_1
plotXY(myGibbsPlot, x_plot, theta_1);
/*
** Change plot layout to 2 x 1 grid and place
** next plot in 2nd cell
*/
plotLayout(2, 1, 2);
// Plot theta_1
plotXY(myGibbsPlot, x_plot, theta_2);
// Turn off plot grid for next plot
plotClearLayout();
// Output matrix concatenating theta_1 and theta_2
ygs = theta_1~theta_2;
ygs = ygs[burn_in+1:burn_in+keep_draws, .];
"---------------------------------------------------";
"Variance-Covariance Matrix from Gibbs Sampler";
"---------------------------------------------------";
"Sample Mean: ";; meanc(ygs);
"Sample Variance-Covariance Matrix: ";
vcx(ygs);