Integrating bivariant Gaussian Distribution by integral2 in Polar coordinates

The following codes' function is to calculate the probability of a point falling in a circle of radius R. When I running the following codes, it will give out the results of P = 1.86e-20, which is very small. This is incorrect, because the reasonble results should be very close to 1. Where is the problem?

R = 6378;
Sigma_bpara = [54.1622939975184	-497.142586464551; -497.142586464551	8986.94991970858];
Mu_bpara = [-0.423193291368079; -934.725277662277];
inputs = [Sigma_bpara(1,1), Sigma_bpara(2,2), Sigma_bpara(1,2), Mu_bpara(1), Mu_bpara(2)];
P = Pc_integral(inputs,R);

function Pc = Pc_integral(inputs,R)
sigma2_x = inputs(1);		
sigma2_y = inputs(2);	
sigma_xy = inputs(3);						
x_m 	 = inputs(4);						
y_m 	 = inputs(5);						

Mu    = [x_m; y_m];
Sigma = zeros(2,2);
Sigma(1,1) = sigma2_x;
Sigma(1,2) = sigma_xy;
Sigma(2,1) = Sigma(1,2);
Sigma(2,2) = sigma2_y;

Pc_function = @(x,y) (1/(2*pi*sqrt(det(Sigma)))) * exp(-0.5*transpose([x;y]-Mu)*inv(Sigma)*([x;y]-Mu));
pcoll=@(theta,r) r.*Pc_function(r.*cos(theta),r.*sin(theta));
Pc=integral2(@(theta,r) arrayfun(pcoll,theta,r), 0,2*pi,0,R);
end

Answer 1

[Paul]

It appears that the default settings for integral2 don't "realize" that the critical region that dominates the integral reisdes in relatively small area far from the origin. This particular problem can be solved by either changing the integration method or doing some work up front to narrow down the integration limits to only cover the region of non-trivial probability density.

Updated code to illustrate both options:

R = 6378;
Sigma_bpara = [54.1622939975184	-497.142586464551; -497.142586464551	8986.94991970858];
Mu_bpara = [-0.423193291368079; -934.725277662277];
inputs = [Sigma_bpara(1,1), Sigma_bpara(2,2), Sigma_bpara(1,2), Mu_bpara(1), Mu_bpara(2)];
P1 = Pc_integral(inputs,R,1);
P2 = Pc_integral(inputs,R,2);
[P1 P2]

ans = 1×2

1.0000 1.0000

function Pc = Pc_integral(inputs,R,option)
sigma2_x = inputs(1);		
sigma2_y = inputs(2);	
sigma_xy = inputs(3);						
x_m 	 = inputs(4);						
y_m 	 = inputs(5);						
Mu    = [x_m; y_m];
Sigma = zeros(2,2);
Sigma(1,1) = sigma2_x;
Sigma(1,2) = sigma_xy;
Sigma(2,1) = Sigma(1,2);
Sigma(2,2) = sigma2_y;
Pc_function = @(x,y) (1/(2*pi*sqrt(det(Sigma)))) * exp(-0.5*transpose([x;y]-Mu)*inv(Sigma)*([x;y]-Mu));
pcoll=@(theta,r) r.*Pc_function(r.*cos(theta),r.*sin(theta));
if option == 1
Pc=integral2(@(theta,r) arrayfun(pcoll,theta,r), 0,2*pi,0,R,'Method','iterated');
else
Pc=integral2(@(theta,r) arrayfun(pcoll,theta,r), 3*pi/2-atan(.1),3*pi/2+atan(.1),500,1500);
end
end

For option 2 I hard-wired in the limits of integration based a plot of the pdf. In reality, one would have to come up with a procedure to pick those limits based on the mean and covariance.

I don't know enough about the integral2() methods to have an informed opinion on whether or 'iterated' is alwasy satisfactory. My preference is to be more insightful in choosing the limits of integration.

[Yirui Wang]

Thank you for your answer!

I try another 'inputs', and found the method 'option == 1' may still doesn't work :(

inputs = [45.5751992476038	396.118890570531	-119.136767607612	-1.39494701212971	-2484.16322726011];

ans = 2.1298e-30

Maybe I can try the method 'option == 2'. Thank you!

[Paul]

This second set of inputs yields a probability density with a critical region that is much, much, smaller in area than the first. So I guess that option 1 fails because the integration steps can't find that very small region. But I was able to get option 2 to work pretty quickly, again just by estimating reasonable limits of integration by eye. I think the more robust approach is to determine the region of high density and then set the limits of integration appropriately for whatever probabilty you're trying to compute, i.e., option 2.