2
$\begingroup$

I am trying to calculate the area under a constrained part of a function built from the empirical CDFs of two distributions. When I calculate the area of either CDF using NIntegrate it works fine, but when I combine them using (for example) the Min of the two functions, it takes much longer and throws an error warning.

Here's a minimal working example.

TestData1 = RandomVariate[NormalDistribution[0.5, 0.1], 100];
TestCDF1 = CDF[EmpiricalDistribution[TestData1], x];
TestCDF2 = Min[CDF[EmpiricalDistribution[TestData1], x], 1 - CDF[EmpiricalDistribution[TestData1], x]];
Data1Area = NIntegrate[TestCDF1, {x, Min[Flatten[TestData1]], Max[Flatten[TestData1]]}, Method -> "TrapezoidalRule"]
Data1Area = NIntegrate[TestCDF2, {x, Min[Flatten[TestData1]], Max[Flatten[TestData1]]}, Method -> "TrapezoidalRule"]

It generates data from a normal distribution and uses it to create the empirical CDF. The integrals of the ECDF and 1-ECDF are correct and no error message, but the Min[ECDF, 1-ECDF] integration produces the following error warning:

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} = {0.595167}. 
NIntegrate obtained 0.08174208369697061` and 7.858621548929746`*^-6 
for the integral and error estimates. >>

None of the suspected reasons listed sounds appropriate for this data. This error occurs no matter which method I choose, so I chose the Trapezoidal method because it should yield an exactly correct answer for CDF data. The value of the answer is the same for repeated integrations of the same data, and the x-value of the error is the same too. But that x-value isn't the peak where ECDF and 1-ECDF intersect.

If anybody knows a way to compute this (seemingly simple) integral without it throwing an error I would love to find out. Due to the step-like behavior of empirical CDFs it would be possible to write my own micro-area summing function by going from point to point and calculating the trapezoid areas, but I'm pretty sure this functionality is built into Mathematica and I just need to learn to use it properly.

$\endgroup$
  • $\begingroup$ I don't know if this is something you want or not, but it seems to solve the problem of your function being somewhat sparse: fun1[y_]:=TestCDF1/.x->y and then Data1Area = NIntegrate[fun1[y], {y, Min[Flatten[TestData1]], Max[Flatten[TestData1]]}, Method -> "TrapezoidalRule"] $\endgroup$ – MathX Apr 16 '16 at 16:00
3
$\begingroup$

TestCDF2 has a number of singularities, each of which might cause a convergence warning message:

(* Data prepared with SeedRandom[0] for reproducibility *)
Plot[TestCDF2, {x, Min[Flatten[TestData1]], Max[Flatten[TestData1]]}, 
 PlotPoints -> 200]

Mathematica graphics

There is probably a limit in both time and number of singularities that NIntegrate will detect on its own. We can pass the singularities explicitly in the integration interval argument. (Note the Automatic method works faster than the trapezoidal rule.)

NIntegrate[TestCDF2, 
 Evaluate@{x, Sequence @@ Sort[Flatten[TestData1]]}]
(*  0.0759967  *)
$\endgroup$
  • $\begingroup$ That's a great answer, and basically the kind of thing I was thinking needed to be done: somehow tell Mathematica where the segments are for piecewise integration. I'm still confused about why Mathematica can handle both the CDFs, but not that (poorly named) function which has all the same data points. But as long as this trick works for all discrete step-like data I'm in business! $\endgroup$ – Aaron Bramson Apr 17 '16 at 12:08
  • $\begingroup$ Also, though I didn't check the times, I did check that Automatic gives the same answer as the trapezoidal rule, which it does. And that answer is slightly (but not insignificantly) different from the output when it encounters an error, so correcting this was more important than just avoiding annoying warning messages. $\endgroup$ – Aaron Bramson Apr 17 '16 at 12:12
  • $\begingroup$ @AaronBramson Try PiecewiseExpand[TestCDF1] (fast) vs. PiecewiseExpand[TestCDF2] (I aborted it after a minute). While simplifying TestCDF2 looks easy to me (as a human), Mathematica might be checking all combinations of the cases or some fraction of them. $\endgroup$ – Michael E2 Apr 17 '16 at 12:21
  • 1
    $\begingroup$ @AaronBramson I didn't check exactly what's going on, but the usual default method is Gauss-Kronrod. It would be applied to each subinterval defined by TestData1. It has "open" sampling, which means it does not include the end points, so that the CDF is not evaluated at the jump discontinuities. The trapezoidal rule is "closed." That means in estimating the error, Gauss-Kronrod will estimate the error to be zero on the first pass, and do no recursive subdivision. The trapezoidal rule will subdivide near the discontinuities, until the precision goal is reached. $\endgroup$ – Michael E2 Apr 17 '16 at 12:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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