Why does this integral have a complex component?

Question

I wanted to find the probability of my normally-distributed random variable being at least 15, so I set up this integral:

Integrate[PDF[NormalDistribution[14, 3.7], x], {x, 15, Infinity}]

Imagine my surprise when I got this result:

0.393476 - 1.75334*10^-15 I

Granted, the angle is miniscule. But why is Mathematica producing this result? This happens even if I type in the PDF manually:

0.107822 E^(-0.036523 (-14 + x)^2)

My TI in rectangular mode returns a result equal to the real component of what Mathematica produces. No complex component.

It doesn't make sense for a real-valued function to have complex area anyway. What's going on here? Whenever I see e's and i's in the same room I tend to blame Euler, but I'm having trouble figuring out what Mathematica is doing.

Using the "proper" Mathematica command for computing this probability returns a value that agrees with my calculator.

Probability[15 <= x, x \[Distributed] NormalDistribution[14, 3.7]]

I just want to know why Mathematica thinks the integral has complex area.

Answer 1

I looked at 5 cases. Conclusion at bottom.

ClearAll[x, a, b];
pdf = PDF[NormalDistribution[14, 3.7], x]

a0 = pdf[[1]]

b0 = pdf[[2, 2, 1]]

Case 1

int1 = Integrate[ a Exp[b (-14 + x)^2], x]

int1 /. {a -> a0, b -> b0}

Simplify[(int1 /. x -> Infinity) - (int1 /. x -> 15)]

Case 2

 int2 = Integrate[ a Exp[b (-14 + x)^2], {x, 15, Infinity}]

 int2 /. {a -> a0, b -> b0}

Case 3

 int3 = Integrate[ a Exp[b (-14 + x)^2] /. {a -> a0, b -> b0}, x]

 Simplify[(int3 /. x -> Infinity) - (int3 /. x -> 15)]

Case 4

 int4 = Integrate[a Exp[b (-14 + x)^2] /. {a -> a0, b -> b0}, {x, 15, k}]

 Limit[int4, k -> Infinity]

Case 5

 int5 = Integrate[a Exp[b (-14 + x)^2] /. {a -> a0, b -> b0}, {x, 15, Infinity}]

 Chop[%]

Conclusion

The error can be Chopped. Complex number show up due the function Erfi which comes due to the term exp[x^2] in the integrand as you can see from above

 Expand[(-14 + x)^2]  (* 196 - 28 x + x^2 *)
 Integrate[Exp[x^2], x]

So this answers the question as to where does the complex number comes from. For example, a complex number can show up like this:

  Erfi[Sqrt[b] ] /. b -> b0

As a general point, it is best to delay as long as possible using floating point values in expressions, and substitute these at the very end to get the "most" accurate result. As you can see from the cases above, when the numbers are substituted at the end, less "noise" was present in the intermediate result.

Answer 2

Just to add a couple of more observations to Nasser's.

Case 6

As Daniel Lichtblau hints at in a comment, if we use an exact 37/10 in place of the approximate 3.7, we get an exact result with a zero imaginary component:

Integrate[PDF[NormalDistribution[14, 37/10], x], {x, 15, Infinity}]
N@%
(*
  1/2 Erfc[(5 Sqrt[2])/37]
  0.393476
*)

Case 7

Such a small imaginary component suggests it is due to round-off error. If it is, then using arbitrary precision numbers, say 3.7`16, instead of MachinePrecision numbers, e.g. 3.7, usually produces a zero component.

Integrate[PDF[NormalDistribution[14, 3.7`16], x], {x, 15, Infinity}]
(*
  0.3934761677149 + 0.*10^-15 I
*)

The zero imaginary component may be removed with Chop as Nasser points out.