Integrating to get volume under bivariate normal distribution

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.

[supernoob]

Update! I found my mistake: I was squaring the variance in the denominators of the exponential but they are already squared. Now everything works as expected. I'll post this as the answer and leave this question up in case it is useful to somebody. If you think I should take it down, let me know.

Answer 1

[supernoob]

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