I have a probability density function p(az, el), where az is an angle that can take on values from -pi to pi, and el is an angle that can take on values from -pi/2 to pi/2. I am using the integral2 function to evaluate the probability obtained by integrating p(az, el) over its entire domain, so I expect the answer to be 1.0. In the code below I expect the result of azElProb_0 and azElProb_1 to both be 1.0, but azElProb_0 evaluates to 0.0. I believe that the issue with the azElProb_0 calculation is that I am integrating over the entire domain of the az (-pi to pi) and el (-pi/2 to pi/2) but the density function is only significant over a very small portion of this domain for the case of interest. The calculation of azElProb_1 restricts the domain of integration to the region of significant density, and returns the correct result (azElProb_1 = 1.0). My question is whether there is a better way to use integral2 for my problem so that I don't have to "artificially" restrict the domain of integration. I have tried different tolerances but have had little luck – integral2 often just fails completely and returns a NaN. While restricting the domain of integration is not an unreasonable solution the primary drawback is that I have to pick that domain to capture the region of significant density, which can be a tricky problem itself. It's far preferable to simply integrate over the entire domain of the density function. I realize that the computation I am doing in this example is silly – I know the answer is 1.0 without doing anything. But I need to make sure I can integrate over my density function using integral2 and get "valid" answers for other quantities I want to compute based on this density function (e.g., expectations and covariances).
Code:
r = [99923.8614955483;1744.17749028302;3489.9496702501];P = [ 10000 0 0; 0 22500 1125; 0 1125 5625];% Integrate over the density function to compute the probability. The az
% angle for this density function varies from -pi to pi. The el angle
% varies from -pi/2 to pi/2. The correct answer for azElProb_0 is 1.0. This
% first evaluation returns azElProb_0 = 0 (as wrong as the answer can get).
azElProb_0 = computeAzElProbability(r, P, -pi, pi, -pi/2, pi/2);% This time I restrict the region of integration to the area where I know the
% density is significant. In this case the resulting azElProb is ~1.0 as
% expected.
az_mean = 1;el_mean = -2;delta = pi / 180;azmin = az_mean * pi / 180 - delta;azmax = az_mean * pi / 180 + delta;elmin = el_mean * pi / 180 - delta;elmax = el_mean * pi / 180 + delta;azElProb_1 = computeAzElProbability(r, P, azmin, azmax, elmin, elmax);function azElProb = computeAzElProbability(r_bdy, P_bdy, azmin, azmax, elmin, elmax) azElProb = integral2(@integrand, azmin, azmax, elmin, elmax, ... 'RelTol', 1e-9, 'Method', 'iterated', 'AbsTol', 0); function p = integrand(az, el) %INTEGRAND Joint az/el density function
%
% P = INTEGRAND(AZ, EL) defines the integrand as the joint az/el density
% function.
% Initialize output
p = 0 * az; % Determinant of error covariance
det_P = det(P_bdy); % Define unit vector specified by az / el. The array rhat is n x 3.
rhat = [cos(az(:)) .* cos(el(:)) sin(az(:)) .* cos(el(:)) -sin(el(:))]; % P_bdy \ rhat' is 3 x n and a is n x 1.
a = 0.5 * dot(rhat, (P_bdy \ rhat')', 2); % rhat is n x 3 and P_bdy \ r_bdy(:) is 3 x 1 so b is n x 1
b = -0.5 * rhat * (P_bdy \ r_bdy(:)); % The value c is a scalar
c = 0.5 * r_bdy(:)' * (P_bdy \ r_bdy(:)); c1 = cos( el(:) ) / (2 * pi)^(3 / 2) / sqrt(det_P) / 4 ./ a.^(5 / 2); c2 = sqrt(pi) * ( 2 * b.^2 + a ) .* (1 - erf( b ./ sqrt(a) ) ); c3 = 2 * b .* sqrt(a) * exp(-c); p(:) = c1 .* (c2 .* exp(b.^2 ./ a - c) - c3); endend
Best Answer