Bayesian Monte Carlo Integration

Introduction - bivariate normal Monte Carlo integration

In this tutorial, we will

  1. Draw from the multivariate normal distribution.
  2. Calculate sample mean and the variance-covariance matrix.
  3. Plot a histogram of the $\theta_1$ and $\theta_2$ samples.

Background

Monte Carlo integration is a simple but rarely feasible method for estimating parameters using an assumed posterior distribution. Monte Carlo integration requires that the posterior distribution can be directly drawn from, something that occurs only in the few special cases. In these cases, model parameters are computed by drawing a random sample from the posterior and averaging the appropriate functions of the draws. Given that $g(\theta)$ is a function of the parameter of interest,

$$\hat{g} = \frac{\sum_{r=1}^{R} g(\theta^{(r)})}{R}$$.

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$.

Compute draws from a bivariate random normal distribution

Prior to drawing from the distribution, we set the specifications of the Monte Carlo integration such that we will have 10,000 total draws. In addition, we will set $\rho$ equal to 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;

In order to induce the proper correlation from our standard normal 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 specified matrix.

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

Next, we draw our 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 normals
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 implicit 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. GAUSS utilizes a plotControl structure to programmatically format graph appearances. This structure is a convenient package that stores all of the information about how you would like a specific graph to be displayed. Using a plotControl structure requires two easy steps:

  1. Declare an instance of the structure.
  2. Initialize the structure instance by filling it with default values.
//Declare 'myMCPlot' to be an instance of a plotControl structure
struct plotControl myMCPlot;

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

Once these steps are completed you may change any of the plot settings in the structure using the suite of plotSet commands.

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

Finally, because we have two variables to graph we may want to combine two separate plots, one for $\theta_1$ and one for $\theta_2$. We do this by changing the graph layout using plotLayout. plotLayout divides the graph canvas into a grid with a specified number of rows and columns, and places the next drawn graph into a specified cell in this grid. plotLayout takes the following input arguments

r
Scalar, number of rows in the grid used to divide the graph canvas
c
Scalar, number of columns in the grid used to divide the graph canvas
location
Scalar, cell location in which to place the next created graph.
//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");

//Keep 2x1 grid. 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():

Have a Specific Question?

Get a real answer from a real person

Need Support?

Get help from our friendly experts.

REQUEST A FREE QUOTE

Thank you for your interest in the GAUSS family of products.

© Aptech Systems, Inc. All rights reserved.

Privacy Policy | Sitemap