############################################################### # # Calculation of bivariate normal distribution using # monte-carlo integration mu<-c(0,0) sigma<-c(1,1) rho<-.5 n<-1000000 lower.x<-(-.5) upper.x<-(.5) lower.y<-(-.5) upper.y<-(.5) # Function to draw n pairs of (x,y) from a bivariate # normal distribution. bvsim<-function(n,mu,sigma,rho) { x<-rnorm(n, mu[1], sigma[2]) y<-rnorm(n,mu[2]+(rho*sigma[2]/sigma[1]*(x-mu[1])), sigma[2]*sqrt(1-rho^2)) cbind(x,y) } # Call our function to simulate n pairs of (x,y) from # a bivariate normal distribution defined by # mu, sigma and rho x.y.dat<-bvsim(n, mu, sigma, rho) # Plot first 1000 points in scatterplot # with best fit least squares prediction line plot(x.y.dat[1:1000,]) x.y.fit<-lsfit(x.y.dat[1:1000,2], x.y.dat[1:1000,1]) abline(x.y.fit) prob.monte.carlo<-sum(x.y.dat[,1]<=upper.x & x.y.dat[,1]>=lower.x & x.y.dat[,2]<=upper.y & x.y.dat[,2]>=lower.y)/n prob.monte.carlo ################################################################### # # Here we are checking our monte-carlo integration # results against another method for integrating # the bivariate normal distribution. # # Calculation of multivariate normal # using function "pmvnorm". This function is used # to calculate the probabilities of the bivariate normal dist. # # This function computes probabilities based on methods discussed # in: # # Genz, A. (1992). Numerical computation of multivariate normal # probabilities. Journal of Computational and Graphical # Statistics, 1, 141–150 # Genz, A. (1993). Comparison of methods for the computation of # multivariate normal probabilities. Computing Science and # Statistics, 25, 400-405 # mean<-c(0,0) corr<-matrix(c(1, .5, .5, 1), ncol=2) corr lower<-c(-.5, -.5) upper<-c(.5, .5) library(mvtnorm) prob.pmvnorm<-pmvnorm(lower, upper, mean, corr) prob.pmvnorm