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: θ’=[θ’1 θ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(θ1|y, θ2(0)).
  • Draw θ2(r) from pθ1|y, θ1(r-1).
  • Draw θ1(r) from p(θ1|y, θ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

where θ1 and θ2 are unknown parameters of the model, while p is the known posterior correlation between θ1 and θ2.

 

Using the properties of the multivariate normal distribution

and

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 p 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, θ1 and θ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 θ1(1) = θ2(1). 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 θ1 and one for θ2 using the plotLayout command. In our first cell, we will plot the pay of θ1 against the iteration number: 

** Change plot layout to 2 x 1 grid and place
 ** next plot in 1st cell
 */
 plotLayout(2, 1, 1);

© 2024 Aptech Systems, Inc. All rights reserved.

 

 /*
 
// 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 θ2 next to the graph of θ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:

Plot of bivariate normal gibbs sampler variable paths.

 

 

 

 

 

 

 

 

 

 

 

 

 

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);