How can I get integral2 to work for my problem?

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

    end

end

Answer 1

[nsunjaya]

The "integral2" function in MATLAB was designed to have a more generic interface in order to house more than one method of integration. Unfortunately, it was not designed to offer a lot of custom tuning such as turning off the singularity-weakening transform (which makes the 'tiled' method slightly different to the "TwoD" function, and ultimately causes it to return unexpected result when integrating your PDF).

Getting an adaptive quadrature algorithm successfully started is a ubiquitous problem. Different implementations will exhibit different behaviors on problems that are challenging to start integrating, challenging because the integrand is zero or almost zero almost everywhere. Generally speaking, it is always possible to engineer a failure on such problems by modifying the limits of integration and/or the intensity of the spike. There will be different “tipping points”, and the success of one implementation or the other in a given circumstance usually doesn’t mean much; a lucky sample point here, a more conservative start there, these things matter but on the whole don’t make one method clearly better than the other in general. With 1D integration we provide a “waypoints” feature that can help, but you would need to know where those points are to supply them. Since we don’t have a waypoints feature for 2D integration, we would want to locate those spikes, split the region of integration, and put those spikes in small regions, perhaps on boundaries. That will make your "computeAzElProbability" function considerably more complicated as it figures out where to split the region and then performs and adds up integrations over subregions, but ultimately that will make it much more reliable.

As a workaround, you might want to use the "quad2d" function in MATLAB instead:

   >> azElProb = quad2d(@integrand, azmin, azmax, elmin, elmax,  'RelTol', 1e-9, 'AbsTol', 0, 'singular', false);

You can refer to the following blog post to learn more about some of the differences when using "integral2" and "quad2d": https://blogs.mathworks.com/loren/2014/02/12/double-integration-in-matlab-methods-and-handling-discontinuities-singularities-and-more/

Answer 2

[John]

Interestingly, I downloaded the TwoD.m code written by L.F. Shampine referenced in the integral2 documentation and it has no problem returning the result of 1.0 when integrating over the entire domain of my az and el angles. I just replaced the integral2 function call with TwoD(@integrand, azmin, azmax, elmin, elmax, 'reltol', 1e-9, 'abstol', 0, 'vectorized', true) and everything worked fine. I can still find cases where TwoD has trouble with the integrals I need to evaluate for my problems but TwoD still seems far more robust than integral2.