Mathematica failing to compute function to calculate integral over a region bounded by straight line

Question

Statement of the problem

Consider the following situation:

• You have a model which employs a bivariate distribution with known parameters.
• You have a random realization from the distribution consisting of a pair of random variables.

• You wish to calculate a probability associated with those two random variables using the following methodology:

  1. Calculate the value V of the PDF at the value of the random realization.
  2. Determine the equation of the hyperplane that is tangential to the PDF at that value.
  3. Set the z value of the hyperplane equal to V, resulting in a linear equation only involving the x and y variables.
  4. The probability to be calculated is the volume of the PDF over the region bounded by the straight line, where that probability is, broadly speaking, away from the mean value of the distribution, so that we choose the value that is less than 0.5. For convenience I will call this the tangential probability.

The above is not terribly difficult and is what I asked for help with in here: Integrating a bivariate distribution over a region bounded by a straight line.

The problem that I am asking for help with in this post involves finding all points (that would have been random realisations in the formulation above) such that the tangential probability equals alpha, where alpha is typically a smallish value such as 1%.

I envisage that the solution would involve finding the region defined by the tangents as illustrated in the figure below:

Attempted solution

The following code (using Mathematica v8) attempts to find the tangential probability associated with an arbitrary point (ax,by) for a Mulitvariate normal distribution:

bivModel = MultinormalDistribution[{-1, 1}, {{1, 0.5 2}, {0.5 2, 4}}]
myF[x_Real, y_Real] := PDF[bivModel, {x, y}];
    NewF[ax_Real, by_Real] := 
     NIntegrate[
      PDF[bivModel, {x, y}]*
       Boole[y >= -(D[myF[x, y], x] /. {x -> ax, y -> by})/(D[
               myF[x, y], y] /. {{x -> ax, y -> by}})*(x - ax) + 
          by], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}, 
      WorkingPrecision -> 7]

(The WorkingPrecision specification was suggested in the solution to the previous post.) Then region of interest for alpha = 0.01 should be possible to visualize using:

RegionPlot[NewF[ax, by] >= 0.01, {ax, -7, 3}, {by, -2, 4}]

Unfortunately I can't even get NewF to work properly. For example, when I calculate NewF[2.7, 1.6] I expect to get 0.01550654 but instead obtain an error:

NIntegrate::inumr: The integrand 0.0918881 E^(1/2 (-(1+x) (1.33333 Plus[<<2>>]-0.333333 Plus[<<2>>])-(-0.333333 Plus[<<2>>]+0.333333 Plus[<<2>>]) (-1+y))) Boole[y>={1.6 +4.58065 (-2.7+x)}]
 has evaluated to non-numerical values for all sampling points in the region with boundaries {{-\[Infinity],0},{-\[Infinity],0}}. >>

Do you know why the above code is not working?

The following gives further detail. It shows step by step derivation of the assertion that the result should give 0.01550654, including a visualization of the area over which the integration occurs:

contourval = myF[-2.7, 1.6];
myyVec = D[myF[xz, yz], yz] /. {xz -> -2.7, yz -> 1.6};
myxVec = D[myF[xz, yz], xz] /. {xz -> -2.7, yz -> 1.6};
Show[RegionPlot[
  PDF[bivModel, {x, y}] >= contourval, {x, -7, 3}, {y, -4, 7}], 
 Plot[-myxVec/myyVec*(x + 2.7) + 1.6, {x, -7, 3}], 
 Graphics[Point[{-2.7, 1.6}]]];
RegionPlot[
 y >= -myxVec/myyVec*(x + 2.0) + 1.0, {x, -7, 3}, {y, -2, 4}];
NIntegrate[
 PDF[bivModel, {x, y}]*
  Boole[y >= -myxVec/myyVec*(x + 2.7) + 1.6], {x, -Infinity, 
  Infinity}, {y, -Infinity, Infinity}, WorkingPrecision -> 7]

This gives the following output:

The figure immediately below shows the region bounded by the contour line corresponding with the random realization at (2.7,1.6).

The region over which the integral is to be calculated is given in the following figure:

The output value of the integration is 0.01550654.

Answer 1

I corrected the problem code, so that the integral NewF[-2.7, 1.6]is calculated even more accurately than at WorkingPrecision -> 7. But when calculating a domain, division by 0 occurs, since the derivative Derivative[0, 1][myF][ax, by] vanishes at some points in the desired domain.

bivModel = MultinormalDistribution[{-1, 1}, {{1, 1}, {1, 4}}];
myF[x_, y_] := PDF[bivModel, {x, y}];
NewF[ax_, by_] := 
 NIntegrate[
  myF[x, y]*
   Boole[y >= -Derivative[1, 0][myF][ax, by]/
        Derivative[0, 1][myF][ax, by]*(x - ax) + by], {x, -Infinity, 
   Infinity}, {y, -Infinity, Infinity}]


NewF[-2.7, 1.6]

Out[]= 0.0154967