Bayesian Monte Carlo Integration

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:

1. Draw from the multivariate normal distribution.
2. Calculate the sample mean and the variance-covariance matrix of our distributionn.
3. Plot a histogram of the $\theta_1$ and $\theta_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, $\theta_1$ and $\theta_2$, drawn directly from a bivariate normal 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$.

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 $\rho = 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, $\theta_1$ in the first column and $\theta_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; 

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

// 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 $\theta_1$.

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

Now plot $\theta_2$ below the graph of $\theta_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 $\theta_1$ and $\theta_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();