GOALS

In this tutorial we will look at using Monte Carlo integration to draw from a bivariate normal distribution. By the end of this tutorial we will:

  • Draw from the multivariate normal distribution.
  • Calculate the sample mean and the variance-covariance matrix of our distributionn.
  • Plot a histogram of the θ1 and θ2 samples.

 

INTRODUCTION 

Monte Carlo integration is a simple but rarely feasible method for estimating parameters using an assumed posterior distribution. The difficulty of Monte Carlo integration is that it requires that the posterior distribution can be directly drawn from. This is something that occurs only in a few special cases.

 

In the feasible Monte Carlo cases, model parameters are computed by drawing a random sample from the posterior and averaging the appropriate functions of the draws.

 

In this particular tutorial, we consider two parameters, θ1 and θ2, drawn directly from a bivariate normal distribution such that

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

 

SET SAMPLER PARAMETERS 

Prior to drawing from the distribution, we set the parameters of the Monte Carlo integration such that:

  • There are 10,000 total draws.
  • The correlation parameter p = 0.6.
 
 /*
 ** Monte Carlo integration specifications
 */
 // Known correlation parameter
 rho = 0.6;

// Draws to keep from sampler
 keep_draws = 10000;

// Variance-covariance matrix
 v = 1~rho|rho~1;

SETTING UP CORRELATION MATRIX 

In order to induce the proper correlation for our draws, we need to find the matrix p such that

v = p’p

 

We can do this in GAUSS using the function chol which finds the Cholesky decomposition of a specified matrix.

 
 // Cholesky decomposition s.t. p'p = v
 p = chol(v);

 

IMPLEMENTING THE SAMPLER 

Next, we draw out the samples of our two parameters, θ1 in the first column and θ2 in the second column, from the standard normal distribution using the GAUSS rndn command. The proper correlation is induced using the matrix ‘p’ which we derived from the variance-covariance matrix.

 
 // Set random seed for repeatable random numbers
 rndseed 5234;

// Generate the standard normal draws
 theta = rndn(keep_draws, 2);

// Generate the correlation by appropriate rotation
 theta_corr = theta*p;

 

Note: For increased efficiency in real-world projects, you can use the built-in function rndMVn to compute random draws from a multivariate normal distribution.

 

ESTIMATE THE MEAN AND VARIANCE-COVARIANCE MATRIX 

Finally, we are able to estimate our parameters from the mean of our samples using the GAUSS meanc command and the sample variance-covariance using the GAUSS command vcx.

 
 "---------------------------------------------------";
 "Variance-Covariance Matrix from Cholesky Simulation";
 "---------------------------------------------------";
 "Sample Mean: ";; meanc(theta_corr);
 "Sample Variance-Covariance Matrix: ";
 vcx(theta_corr);

 

The print statements above produce the following output.

---------------------------------------------------
Variance-Covariance Matrix from Cholesky Simulation
---------------------------------------------------
Sample Mean:
 -0.0215
 -0.0255
Sample Variance-Covariance Matrix:

  1.0128   0.6063
  0.6063   0.9998

 

PLOT THE POSTERIOR DISTRIBUTION

To better understand the behavior of our Monte Carlo integration, we can plot the posterior distribution of our parameters. We will use two, side-by-side, histogram plots to visualize our random variables.

 

// Declare 'myMCPlot' to be an instance of a plotControl structure
 struct plotControl myMCPlot;

// Fill myMCPlot with 'histogram' defaults
 myMCPLot = plotGetDefaults("hist");

// Set the text interpreter to LaTex for the axis labels
 plotSetTextInterpreter(&myMCPlot, "latex", "axes");

// Set plot title
 plotSetTitle(&myMCPlot, "Bivariate Monte Carlo Integration", "Arial", 18);

 

© 2024 Aptech Systems, Inc. All rights reserved.

 

// Set axis labels
 plotSetYLabel(&myMCPlot, "\\theta_1", "Arial", 14);

// Change plot layout to 2 x 1 cells
 plotLayout(2, 1, 1);

 

In our first graph cell, we will plot a histogram of our draws of θ1.

 // Plot theta_1
 plotHist(myMCPlot, theta_corr[., 1], 20);

 

Now plot θ2 below the graph of θ1.

// Turn off title for second plot
 plotSetTitle(&myMCPlot, "");

/*
 ** Set axis labels. Font family and size will remain unchanged
 ** from previous call to 'plotSetYLabel' above
 */
 plotSetYLabel(&myMCPlot, "\\theta_2");

// Place next plot in second cell.
 plotLayout(2, 1, 2);

// Plot theta_2
 plotHist(myMCPlot, theta_corr[., 2], 20);

// Turn off layout grid for future plots
 plotClearLayout();

 

 

 

 

 

 

 

 

 

 

CONCLUSION 

Congratulations! You have successfully implemented Monte Carlo integration including:

  • Drawing from the multivariate normal distribution.
  • Calculating the sample mean and the variance-covariance matrix.
  • Plotting a histogram of the θ1 and θ2 samples.
  • The next tutorial introduces the Importance Sampling.

 

For your convenience, the entire code is below.

/*
 ** Monte Carlo integration specifications
 */
 // Known correlation parameter
 rho = 0.6;

// Draws to keep from sampler
 keep_draws = 10000;

// Variance-covariance matrix
 v = 1~rho|rho~1;

// Cholesky decomposition s.t. p'p = v
 p = chol(v);

// Set random seed for repeatable random numbers
 rndseed 5234;

// Generate the standard normal draws
 theta = rndn(keep_draws, 2);

// Generate the correlation by appropriate rotation
 theta_corr = theta*p;

"---------------------------------------------------";
 "Variance-Covariance Matrix from Cholesky Simulation";
 "---------------------------------------------------";
 "Sample Mean: ";; meanc(theta_corr);
 "Sample Variance-Covariance Matrix: ";
 vcx(theta_corr);

// Declare 'myMCPlot' to be an instance of a plotControl structure
 struct plotControl myMCPlot;

// Fill myMCPlot with 'histogram' defaults
 myMCPLot = plotGetDefaults("hist");

// Set the text interpreter to LaTex for the axis labels
 plotSetTextInterpreter(&myMCPlot, "latex", "axes");

// Set plot title
 plotSetTitle(&myMCPlot, "Bivariate Monte Carlo Integration", "Arial", 18);

// Set axis labels
 plotSetYLabel(&myMCPlot, "\\theta_1", "Arial", 14);

// Change plot layout to 2 x 1 cells
 plotLayout(2, 1, 1);

// Plot theta_1
 plotHist(myMCPlot, theta_corr[., 1], 20);

// Turn off title for second plot
 plotSetTitle(&myMCPlot, "");

/*
 ** Set axis labels. Font family and size will remain unchanged
 ** from previous call to 'plotSetYLabel' above
 */
 plotSetYLabel(&myMCPlot, "\\theta_2");

// Place next plot in second cell.
 plotLayout(2, 1, 2);

// Plot theta_2
 plotHist(myMCPlot, theta_corr[., 2], 20);

// Turn off layout grid for future plots
 plotClearLayout();