# Question on accuracy and precision of NIntegrate

As a Mathematica newbie, I was testing the accuracy/precision of NIntegrate (9.0.1.0 on Mac) and have obtained a very peculiar result.

f[x_] := (1/2) PDF[NormalDistribution[-100, 1], x]
+ (1/2) PDF[NormalDistribution[+100, 1], x]
g[n_] := NIntegrate[f[x] Log2[1/f[x]], {x, -Infinity, Infinity},
AccuracyGoal -> Infinity, PrecisionGoal -> Automatic,
MaxRecursion -> 1000, WorkingPrecision -> n]
Plot[g[n], {n, 50, 100}]


The correct value, which cannot be calculated analytically, should be around 3.0471. There was no error message for any value of n.

I am not interested in calculating this specific integral, but I am curious whether there is any way I can feel assured that the numerical value Mathematica returns is correct and by how much. In comparison, the GNU Scientific Library gives me an error bound (although it uses machine floating point numbers, which can be another source of trouble).

-
The issue seems to have nothing to do with precision or accuracy, but with sampling. – Michael E2 Jan 18 '14 at 19:58
@MichaelE2 What do you mean by sampling? – numerical-integrator Jan 18 '14 at 20:00
Numerical integration is done by evaluating the function at a few points in its domain. The points are called sampling points. See NIntegrate Integration Strategies. For instance, g[n_] := NIntegrate[f[x] Log2[1/f[x]], {x, -Infinity, Infinity}, Method -> "DoubleExponential", MaxRecursion -> 20, WorkingPrecision -> n] seems more reliable. – Michael E2 Jan 18 '14 at 20:05
@MichaelE2 Thanks for the explanation. What I am worried about, however, is that without any complaint the answer could be completely off. I clearly asked NIntegrate to give me an answer within certain accuracy/precision. Is there any way we can bound the error? – numerical-integrator Jan 18 '14 at 20:09
I agree that it is a sampling problem. Compare the results of NIntegrate[ f[x] Log2[1/f[x]], {x, -\[Infinity], -100, 100, \[Infinity]}] and NIntegrate[f[x] Log2[1/f[x]], {x, -\[Infinity], \[Infinity]}]. By using the {x, x1, x2, ..., xn} form for the range of integration (see the documentation "Details" for NIntegrate) to suggest the location of the "singularities", NIntegrate can then compute the correct result. – Stephen Luttrell Jan 18 '14 at 21:50

The integrand is nearly zero except for two maxima near -100 and 100. For some reason I do not understand, changing the precision causes the sampling to miss the rather small regions where the integrand is significantly greater than zero. (Changing the precision changes the calculations, and perhaps that is the whole reason.)

A fix would be to control the sampling to include the peaks.

## Illustration of the sampling

Here is a look at the function and value away from peaks -- the value is quite small!:

f[x_] := (1/2) PDF[NormalDistribution[-100, 1], x] +
(1/2) PDF[NormalDistribution[+100, 1], x]

Plot[f[x] Log2[1/f[x]], {x, -150, 150}, PlotRange -> All]


{f[10], f'[10]} // N
(*
{2.554043924441158*10^-1760, 2.298639531997042*10^-1758}
*)


Here is a look at the sampling for two settings of WorkingPrecision, using Sow to stored each value of x at which the integrand is evaluated:

g[n_] := NIntegrate[f[x] Log2[1/f[x]], {x, -Infinity, Infinity},
MaxRecursion -> 20, WorkingPrecision -> n,
EvaluationMonitor :> Sow[x]];

{int, samp} = Reap @ g[62];
int
ListPlot[samp, PlotRange -> {-400, 400}]
(*
3.0470955851806411027016019880549497362225614404580137317711359
*)


{int, samp} = Reap @ g[63];
int
ListPlot[samp, PlotRange -> {-400, 400}]
(*
1.52354779259032055135080099402747486811128072022900686588556600
*)


We can see above that the first integral found both maxima and the second found only one. Mathematica estimates the error from the function evaluations (perhaps from its derivatives, too -- I'm ignorant of the actual workings). If it misses a peak and all the values are extremely close to zero, it's likely that the estimate of the error would be nearly zero.

One thing to notice is that the sampling is biased. Initial sampling starts near zero and grows exponentially.

## Getting NIntegrate to include the peaks

This has to be done "by hand" in some sense. If you know about where the peaks are, you can use FindMaximum to get closer to them. If it isn't known where they all are, then it may be that not enough is known yet to solve the problem.

As Stephen Luttrell points out in a comment, you can specify extra points in the domain. These are tested to see if there are singularities, but even if this function passes the test, the side effect is that the peaks are included in the sampling.

g2[n_] :=
NIntegrate[f[x] Log2[1/f[x]], {x, -Infinity, -100, 100, Infinity},
PrecisionGoal -> Automatic, WorkingPrecision -> n,
EvaluationMonitor :> Sow[x]]

{int, samp} = Reap@g2[63];
int
ListPlot[samp, PlotRange -> {-400, 400}]
(*
3.04709558518064110270160198805494973622256144045801373177111896
*)


-
I really appreciate your detailed analysis of what's going wrong and suggested workaround. However, I am still wondering how many digits from the answer "3.04709558518064110270160198805494973622256144045801373177111896" are indeed correct. – numerical-integrator Jan 18 '14 at 23:22