Why is mathematica giving wrong and incomparable results for the integral?

Question

1) Integration of Gaussian Distribution with $(x,y,z)$ ranging from $-\infty$ to $\infty$ gives 1 as expected using this command in mathematica. (Total Probability = 1)

$\sigma = 200000$ and $\mu=1000000$

Integrate[(E^(-(x - 1000000)^2/(2*200000^2))*E^(-(y - 1000000)^2/(2*200000^2))*
          E^(-(z -1000000)^2/(2*200000^2))/(200000*Sqrt[2*Pi])^3) , {x, -Infinity, Infinity}, 
          {y, -Infinity,    Infinity}, {z, -Infinity, Infinity}]

2) Now I want to integrate over $(x,y,z)$ such that $$\frac{x}{1+0.1}+\frac{y}{(1+0.1)^2}+\frac{z}{(1+0.1)^3}\geq 2000000 $$ given that $$x,y,z\geq0$$ So I wrote this in mathematica:

Integrate[(E^(-(x - 1000000)^2/(2*200000^2))*E^(-(y - 1000000)^2/(2*200000^2))
       *E^(-(z -1000000)^2/(2*200000^2))/(200000*Sqrt[2*Pi])^3) Boole[x/(1 + 0.1) + y/(1 + 0.1)^2
         + z/(1+ 0.1)^3 - 2000000 >= 0], {x, 0,   Infinity}, {y, 0, Infinity}, {z, 0, Infinity}]

But the answer I am getting is some very big number ($1.7\times10^{17}$)

Can anyone point out the mistake?

Answer 1

There are two ways come to my mind to go.

1. Truncate the upper limits:

Since OP has an exponential decay term like what it reads, truncating the integral limits at 2000000 should give a reasonably precise result:

NIntegrate[(
            E^(-(x - 1000000)^2/(2*200000^2))
             *E^(-(y - 1000000)^2/(2*200000^2))
             *E^(-(z -1000000)^2/(2*200000^2))
             /(200000*Sqrt[2*Pi])^3
           )
           Boole[x/(1 + 0.1) + y/(1 + 0.1)^2 + z/(1+ 0.1)^3 - 2000000 >= 0],
           {x, 0, 2000000}, {y, 0, 2000000}, {z, 0, 2000000}
          ]

0.9545159272869793`

2. Normalize the integral:

This is my preferred way.

As anderstood commented, because of the magnitude of the coefficients in the exponential term, there are high chance we encounter numerical problems like overflow, precision losing etc. Increasing WorkingPrecision is expensive, the more preferred way would be to normalize the integral.

By introducing variables like $a = \frac{x-1000000}{200000\sqrt{2}}$ etc., the original integral will become:

1/π^(3/2) NIntegrate[
                     E^(-a^2 - b^2 - c^2)
                     Boole[121 a + 110 b + 100 c + 324/Sqrt[2] >= 0], 
                     {a, -(5/Sqrt[2]), ∞}, {b, -(5/Sqrt[2]), ∞}, {c, -(5/Sqrt[2]), ∞}
                    ]

0.9545167785163403`

Answer 2

Below are a couple of ways to add to Silvia's two. First a couple of remarks.

• Using 1/10 instead of 0.1 in the specification of the region allows Mathematica to apply exact methods. This may help, but it can also add significantly to computation time sometimes. What is important is to realize that there is a difference and to become familiar with the advantages and disadvantages of each.
• NIntegrate tends to begin its analysis of the function's relation to the domain by scaling out when the domain is infinite. This has a tendency to make the integration more sensitive to the endpoints/boundaries, if any, and to the origin. In your case, the center of the interesting behavior of integrand is relatively far away from the boundaries and origin. This is why Silvia's second method works. A translation without scaling would also work.

3. Compute the complement

To do this, one has to be able to compute the following, which can be done exactly:

Integrate[(E^(-(x - 1000000)^2/(2*200000^2))*
   E^(-(y - 1000000)^2/(2*200000^2))*
   E^(-(z - 1000000)^2/(2*200000^2))/(200000*Sqrt[2*Pi])^3), {x, 0, 
  Infinity}, {y, 0, Infinity}, {z, 0, Infinity}]
(* 1/8 (1 + Erf[5/Sqrt[2]])^3 *)

1/8 (1 + Erf[5/Sqrt[2]])^3 - 
 NIntegrate[(E^(-(x - 1000000)^2/(2*200000^2))*
     E^(-(y - 1000000)^2/(2*200000^2))*
     E^(-(z - 1000000)^2/(2*200000^2))/(200000*Sqrt[2*Pi])^3) *
    Boole[x/(1 + 0.1) + y/(1 + 0.1)^2 + z/(1 + 0.1)^3 - 2000000 <= 0],
  {x, 0, Infinity}, {y, 0, Infinity}, {z, 0, Infinity}]
(* 0.954516788746158` *)

Or in V10:

region = ImplicitRegion[
   x/(1 + 1/10) + y/(1 + 1/10)^2 + z/(1 + 1/10)^3 - 2000000 <= 0 && 
    x >= 0 && y >= 0 && z >= 0, {x, y, z}];

1/8 (1 + Erf[5/Sqrt[2]])^3 - 
 NIntegrate[(E^(-(x - 1000000)^2/(2*200000^2))*
    E^(-(y - 1000000)^2/(2*200000^2))*
    E^(-(z - 1000000)^2/(2*200000^2))/(200000*Sqrt[2*Pi])^3),
   {x, y, z} ∈ region]

4. Break down the region

First, the region of integration is fairly simple to analyze algebraically. Even Integrate makes some progress:

prob = Integrate[(E^(-(x - 1000000)^2/(2*200000^2)) *
     E^(-(y - 1000000)^2/(2*200000^2)) *
     E^(-(z - 1000000)^2/(2*200000^2)) / (200000*Sqrt[2*Pi])^3) *
   Boole[x/(1 + 1/10) + y/(1 + 1/10)^2 + z/(1 + 1/10)^3 - 2000000 >= 0],
  {x, 0, Infinity}, {y, 0, Infinity}, {z, 0, Infinity}]
(*
  (1/8)*(1 + Erf[5/Sqrt[2]])^2*Erfc[3*Sqrt[2]] + 
   Integrate[<>, {x, 0, 2200000}] + 
   Integrate[<>, {x, 0, 2200000}, {y, 0, (1/10)*(24200000 - 11*x)}]/(160000000000*Pi)
*)

Then we can get the numerical result (see note N1 below):

N@prob
(* 0.9545168023458417` *)

Second, we can also break the region down ourselves:

redreg = Reduce[
  x/(1 + 1/10) + y/(1 + 1/10)^2 + z/(1 + 1/10)^3 - 2000000 > 0 && 
   x > 0 && y > 0 && z > 0, {x, y, z}]
(*
  (0 < x <= 2200000 &&
      ((0 < y <= 1/10 (24200000 - 11 x) && z > 1/100 (266200000 - 121 x - 110 y)) ||
       (y > 1/10 (24200000 - 11 x) && z > 0))) ||
  (x > 2200000 && y > 0 && z > 0)
*)

LogicalExpand will distribute the And over the Or, so that we can break the region into three pieces:

pieces = Map[
 Reduce[#, {x, y, z}] /. {And -> List, 
    HoldPattern@Inequality[a_, _, v_, _, b_] :> {v, a, b}, 
    v_ > a_ :> {v, a, Infinity}} &,
 (LogicalExpand[redreg] /. Or -> List)
 ]
(*
  {{{x, 2200000, ∞}, {y, 0, ∞}, {z, 0, ∞}},
   {{x, 0, 2200000}, {y, 1/10 (24200000 - 11 x), ∞}, {z, 0, ∞}},
   {{x, 0, 2200000}, {y, 0, 1/10 (24200000 - 11 x)}, {z, 1/100 (266200000 - 121 x - 110 y), ∞}}}
*)

We can then integrate over each of the three pieces and add up the results

ParallelMap[
  NIntegrate[(E^(-(x - 1000000)^2/(2*200000^2))*
        E^(-(y - 1000000)^2/(2*200000^2))*
        E^(-(z - 1000000)^2/(2*200000^2))/(200000*
            Sqrt[2*Pi])^3), ##] & @@ # &,
  pieces
  ] // Total
(* warnings about slow convergence *)
(* 0.9545167832232275` *)

Notes

N1. N@prob takes about 1/20 sec. the first time it is evaluated and almost 16 sec. the second time. (After the first time Mathematica tries to do the unevaluated Integrates again.) One way around the delay is the following which substitutes NIntegrate for Integrate before prob is evaluated:

Hold[prob] /. (OwnValues[prob] /. Integrate -> NIntegrate) // ReleaseHold

N2. The built-in singularity handling of NIntegrate helps, but it is not really up to the task in this case. Adding numbers to the intervals of integration affects the sampling around these numbers. We can add the peak (x = 10^6 etc.) and a point after the peak (x = 2 * 10^6 etc.), or even more. I did not pursue this line of investigation far enough to say much more. By increasing MinRecursion, one can achieve four accurate digits instead of the two in the example below.

Example use:

NIntegrate[(E^(-(x - 1000000)^2/(2*200000^2))*
    E^(-(y - 1000000)^2/(2*200000^2))*
    E^(-(z - 1000000)^2/(2*200000^2))/(200000*Sqrt[2*Pi])^3) *
  Boole[x/(1 + 0.1) + y/(1 + 0.1)^2 + z/(1 + 0.1)^3 - 2000000 >= 0],
 {x, 0, 10^6, 2 10^6, ∞}, {y, 0, 10^6, 2 10^6, ∞}, {z, 0, 10^6, 2 10^6, ∞}]
(* warnings about slow convergence and error not converging *)
(* 0.9596152356949542` *)