MATLAB: Integrating to get volume under bivariate normal distribution

bivariate normal distributiondouble integralmultivariate normal distributionpdfpolar coordinates

Hi there, I have 2 gaussian random variables which together form a bivariate normal distribution. Let's call them A and B. Say I did this:
AB = [A,B]; %A and B are gaussian distributed RV's
AB_mean = mean(AB); 
AB_std = std(AB); 
AB_var = var(AB); 
AB_covar = cov(AB);
I can create and plot the multivariate normal distribution as follows:
xsteps = linspace(-5e-3,5e-3, 1000);
ysteps = xsteps;
[X,Y] = meshgrid(xsteps,ysteps);
F = mvnpdf([X(:) Y(:)],AB_mean,AB_covar);
F = reshape(F,length(ysteps),length(xsteps));
figure(53); s = surf(xsteps,ysteps,Fnorm, 'FaceAlpha',0.5); s.EdgeColor = 'none';
xlabel('x'); ylabel('y'); zlabel('Probability Density');
How do I integrate this distribution in 2d to get the volume under some portion of the surface? I haven't been able to find a way which works. Note what I really want is to integrate over the pdf in polar coordinates. I tried doing this (using the cartesian definition of the pdf from http://mathworld.wolfram.com/BivariateNormalDistribution.html):
AB_mean = mean(AB); mu1 = AB_mean(1); mu2 = AB_mean(2);
AB_std = std(AB); sigma1 = AB_std(1); sigma2 = AB_std(2);
AB_var = var(AB); var1 = AB_var(1); var2 = AB_var(2);
AB_covar = cov(AB);
rmax = 2e-3;
fun = @(x,y) (1/(2*pi*sigma1*sigma2*sqrt(1-rho^2)))...
    .*exp((-1/(2*(1-rho^2)))*(((x-mu1)/var1).^2+((y-mu2)/var2).^2-((2*rho*(x-mu1).*(y-mu2))/(sigma1*sigma2))));
polarfun = @(theta,r) fun(r.*cos(theta),r.*sin(theta)).*r;
q = integral2(polarfun,0,2*pi,0,rmax);
But this always gives me zero, even if I make r very large, and I can't figure out what is wrong. If things were working, I would expect to get 1 for very large r. Please help! I have attached my 2 RV's as a matrix AB.

Best Answer

The mistake is squaring the variance in the denominator of the exponential terms. The variance is already squared. The correct expression for the probability density function is the following:
fun = @(x,y) (1/(2*pi*sigma1*sigma2*sqrt(1-rho^2)))...
    .*exp((-1/(2*(1-rho^2)))*(((x-mu1).^2/var1)+((y-mu2).^2/var2)-((2*rho*(x-mu1).*(y-mu2))/(sigma1*sigma2))));
Related Question