Consider the following code (Mathematica 8):

mu1 = -1.0; mu2 = 1.0; ss1 = 1; ss2 = 4; cov = 1;
bivModel = 
 MultinormalDistribution[{mu1, mu2}, {{ss1, cov}, {cov, ss2}}]
myFun[x_, y_] := PDF[bivModel, {x, y}];
a = -2.7; b = 1.6;
yVec = D[myFun[x, y], y];
xVec = D[myFun[x, y], x]
m = (mu2 - b)/(mu1 - a)
c = mu2 - m*mu1
expr3[z_] := (-xVec //. {x -> z, y -> (m*z + c)})/(yVec //. {x -> z, 
   y -> (m*z + c)})*(x - z) + m*z + c

Using the straight line defined by expr, I want to evaluate a double integral which I do in two steps under Method 1:

Part1F[x_Real, z_Real] := 
 Integrate[PDF[bivModel, {x, y}], {y, expr3[z], Infinity}]

   Part2F[z_Real] := NIntegrate[Part1F[x, z], {x, -Infinity, Infinity}]

So for example when I type in

Part2F[1.0]

I get a series of error messages and a result of 0.805675 - 1.96891*10^-18 I.

The errors are:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {-1.69485}. NIntegrate obtained 0.805675 -1.96891*10^-18 I
 and 0.0005750462582540412` for the integral and error estimates. >>

But under Method 2 when I type in:

Part1K[x_Real, z_Real] := 
 Integrate[
  PDF[bivModel, {x, y}], {y, 
   11/17 - (6 z)/17 - 
    1/(2/3 (6/17 + (6 z)/17) + (2 (1 + z))/3) E^(
     1/2 (-(1/3 (-1 - z) + 1/3 (-(6/17) - (6 z)/17)) (-(6/17) - (6 z)/
            17) - (1 + z) (1/3 (6/17 + (6 z)/17) + (4 (1 + z))/3)) + 
      1/2 ((1/3 (-1 - z) + 1/3 (-(6/17) - (6 z)/17)) (-(6/17) - (6 z)/
            17) + (1 + z) (1/3 (6/17 + (6 z)/17) + (4 (1 + z))/
            3))) (x - z) (2/3 (-(6/17) - (6 z)/17) - (8 (1 + z))/3), 
   Infinity}]

(obtained by pasting in the expansion of exp3[z] into the expression for Part1F) and

Part2K[z_Real] := NIntegrate[Part1K[x, z], {x, -Infinity, Infinity}]
Part2K[1.0]

I get no errors and a result of 0.994423 - 8.13629*10^-21 I.

If under Method 3 I simply start off with a value of z of 1.0 and use the following code:

Part1H[x_Real] := 
 Integrate[
  PDF[bivModel, {x, y}], {y, 
   0.29411764705882343` + 3.2173913043478257` (-1.` + x), Infinity}]NIntegrate[Part1H[x], {x, -Infinity, Infinity}]

I get:

0.994423 - 8.13629*10^-21 I

Questions:

  1. Why am I getting different results between the first method and the second and third?
  2. I want to be able to vary the value of z and indeed the value of m. Is there a way of keeping this functionality? The formulae for Part1F and Part2F are the closest I can get to this but I have more confidence in the results from Part1G and Part2G and Part 1H and its integral.
  3. I am calculating an integral of the bivariate Normal distribution. Why on earth am I getting imaginary numbers (albeit very small)?
  • 2
    The small imaginary part in your results come from using machine arithmetic in your computations. Look up Chop in documentation to see how eliminate them. Or rework your computations to use exact or arbitrary precision arithmetic. – m_goldberg Dec 4 at 22:38
  • Thanks, will give that a go! – GerardF123 Dec 5 at 9:29
up vote 1 down vote accepted
$Version

(* "11.3.0 for Mac OS X x86 (64-bit) (March 7, 2018)" *)

Clear["Global`*"]

To avoid imaginary artifacts resulting from calculating with machine precision, use exact numbers when possible and Simplify intermediate results. (EDIT: Unless told otherwise, Mathematica operates in the complex plane. Calculations which would cancel the imaginary components if done exactly, can result in imaginary artifacts when done with machine precision.)

mu1 = -1; mu2 = 1; ss1 = 1; ss2 = 4; cov = 1;
bivModel = MultinormalDistribution[{mu1, mu2}, {{ss1, cov}, {cov, ss2}}];
myFun[x_, y_] := Evaluate@Simplify@PDF[bivModel, {x, y}]

a = -27/10; b = 8/5;
{xVec, yVec} = D[myFun[x, y], {{x, y}}] // Simplify;
m = (mu2 - b)/(mu1 - a);
c = mu2 - m*mu1;
expr3[z_] := 
 Evaluate@Simplify[(-xVec /. {x -> z, y -> (m*z + c)})/(yVec /. {x -> z, 
         y -> (m*z + c)})*(x - z) + m*z + c]

Note the simplified form of expr3

expr3[z]

(* 1/391 (253 + 1258 x - 1396 z) *)

Since you are using SetDelayed in defining Part1F, use Evaluate so that the integration is only done once rather than for each evaluation of Part1F

Part1F[x_, z_] := 
 Evaluate@Integrate[PDF[bivModel, {x, y}], {y, expr3[z], Infinity}]

Part1F[x, z]

(* (E^(-(1/2) (1 + x)^2) (1 + 
   Erf[(529 - 867 x + 1396 z)/(391 Sqrt[6])]))/(2 Sqrt[2 π]) *)

Part2F[z_Real] := NIntegrate[Part1F[x, z], {x, -Infinity, Infinity}]

sol1 = Part2F[1.0]

(* 0.994423 *)

Part1K[x_, z_] := 
 Evaluate@Integrate[
   PDF[bivModel, {x, y}], {y, 
    11/17 - (6 z)/17 - 
     1/(2/3 (6/17 + (6 z)/17) + (2 (1 + z))/
          3) E^(1/2 (-(1/3 (-1 - z) + 
                1/3 (-(6/17) - (6 z)/17)) (-(6/17) - (6 z)/17) - (1 + 
               z) (1/3 (6/17 + (6 z)/17) + (4 (1 + z))/3)) + 
         1/2 ((1/3 (-1 - z) + 
               1/3 (-(6/17) - (6 z)/17)) (-(6/17) - (6 z)/17) + (1 + 
               z) (1/3 (6/17 + (6 z)/17) + (4 (1 + z))/3))) (x - 
        z) (2/3 (-(6/17) - (6 z)/17) - (8 (1 + z))/3), Infinity}]

The integrals are identical

Part1K[x, z] === Part1F[x, z]

(* True *)

Part2K[z_Real] := NIntegrate[Part1K[x, z], {x, -Infinity, Infinity}]

sol2 = Part2K[1.0]

(* 0.994423 *)

The results are identical

sol1 === sol2

(* True *)
  • Thank you - this is extremely useful. – GerardF123 2 days ago

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.