Integrating bivariant Gaussian Distribution by integral2 in Polar coordinates

5 views (last 30 days)
Yirui Wang
Yirui Wang on 28 Oct 2021
Commented: Paul on 30 Oct 2021
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

Accepted Answer

Paul
Paul on 28 Oct 2021
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.
  2 Comments
Paul
Paul on 30 Oct 2021
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!