# Integrating bivariant Gaussian Distribution by integral2 in Polar coordinates

6 views (last 30 days)

Show older comments

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

##### 0 Comments

### Accepted Answer

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]

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
on 30 Oct 2021

### More Answers (0)

### See Also

### Community Treasure Hunt

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

Start Hunting!