7
$\begingroup$

Bug introduced in 10.2 or earlier and persisting through 11.0.1 or later


Here's a simplified example of what I'm trying to do:

RegionPlot[
  NIntegrate[PDF[NormalDistribution[0, 1], a], {a, 0, y}] >= 0.2,
                {x, -1, 1}, {y, 0.1, 0.7}]

(yes, it doesn't depend on x)

This returns an inscrutable error:

Throw::nocatch: "Uncaught \!\(Throw[\(-Holonomic`DifferentialRootReduceDump`y[NIntegrate`LevinRuleDump`x]\)
+ \*SuperscriptBox[\"Holonomic`DifferentialRootReduceDump`y\", \"\[Prime]\",
MultilineFunction->None][NIntegrate`LevinRuleDump`x],
NIntegrate`LevinRuleDump`FastLookupHolonomicDifferentialEquation]\) returned to top level."

I'm using Mathematica 10.2. Any ideas what's wrong?

EDIT: I'm not looking for algebraic simplifications or substituting Integrate for NIntegrate. The above code is just an example to reproduce the error. In the code I actually want to run, NIntegrate is the only option.

$\endgroup$
  • 1
    $\begingroup$ Have you tried to figure out where the error come from? To izolate the problem try divide and conquer approach first... $\endgroup$ – mmal Sep 10 '15 at 20:02
  • $\begingroup$ I can reproduce that error with 10.1, even defining a pure numeric function to pass to RegionPlot. The NIntegrate by itself works fine. B-G? $\endgroup$ – george2079 Sep 10 '15 at 20:34
10
$\begingroup$

This is a bug in RegionPlot. For a possible workaround, try the following undocumented option

RegionPlot[NIntegrate[PDF[NormalDistribution[0, 1], a], {a, 0, y}] >= 0.2, 
   {x, -1, 1}, {y, 0.1, 0.7}, "NumericalFunction" -> False]

Mathematica graphics

$\endgroup$
  • $\begingroup$ I cannot find documentation of the option "NumericalFunction" for RegionPlot nor Graphics $\endgroup$ – Bob Hanlon Sep 11 '15 at 16:22
  • $\begingroup$ It seems that 11.3 introduces additional error NIntegrate::nlim: a = y is not a valid limit of integration. and also the bug is still not solved, any idea that why the new error shows up? $\endgroup$ – luyuwuli Feb 23 at 16:21
5
$\begingroup$

EDIT : Changed for your edited question

$Version

(* "10.2.0 for Mac OS X x86 (64-bit) (July 7, 2015)" *)

Define a helper function that is defined only for numeric arguments

f[y_?NumericQ] :=
  NIntegrate[
   PDF[NormalDistribution[0, 1], a],
   {a, 0, y}];

rgn = ImplicitRegion[
   f[y] >= 0.2 && -1 <= x <= 1 && 0.1 <= y <= 0.7,
   {x, y}];

However, this is very sloo...oow

RegionPlot[rgn, PlotRange -> {{-1, 1}, {0.1, 0.7}}] // 
  AbsoluteTiming // Column

enter image description here

ContourPlot is much, much faster

f[0.62] >= 0.2

(* True *)

ContourPlot[f[y],
   {x, -1, 1}, {y, 0.1, 0.7},
   Contours -> {0.2},
   Epilog -> Text["f[y] \[GreaterEqual] 0.2", {0, 0.62}]] //
  AbsoluteTiming // Column

enter image description here


The largest that your integral can be is for y = 0.7

dist = NormalDistribution[0, 1];

Integrate[PDF[dist, a], {a, 0, 0.7}]

(* 0.258036 *)

This is equivalent to

CDF[dist, 0.7] - CDF[dist, 0]

(* 0.258036 *)

Even if you were to integrate from -Infinity, the largest that the integral could be is

Integrate[PDF[dist, a], {a, -Infinity, 0.7}]

(* 0.758036 *)

or equivalently,

CDF[dist, 0.7] - CDF[dist, -Infinity]

(* 0.758036 *)

or more simply

CDF[dist, 0.7]

(* 0.758036 *)

Consequently, since you are looking for the region for which the integral is greater than or equal to 0.95, your region is empty. Note the use of Integrate rather than NIntegrate

RegionPlot[
 Integrate[PDF[NormalDistribution[0, 1], a], {a, 0, y}] >= 0.95, {x, -1, 
  1}, {y, 0.1, 0.7}]

enter image description here

If you reverse the inequality then

RegionPlot[
 Integrate[PDF[dist, a], {a, 0, y}] < 0.95, {x, -1, 1}, {y, 0.1, 0.7}]

enter image description here

Or the same result with

RegionPlot[CDF[dist, y] - CDF[dist, 0] < 0.95, {x, -1, 1}, {y, 0.1, 0.7}]
$\endgroup$
  • $\begingroup$ be that as it may, we shouldn't expect an "inscrutable" error message. $\endgroup$ – george2079 Sep 11 '15 at 2:30
  • 1
    $\begingroup$ Yes, it has nothing to do with NormalDistribution, CDF or PDF. Try RegionPlot[NIntegrate[Sin[a], {a, 0, y}] >= 0.95, {x, -1, 1}, {y, 0.1, 0.7}] vs RegionPlot[NIntegrate[a, {a, 0, y}] >= 0.95, {x, -1, 1}, {y, 0.1, 0.7}]. $\endgroup$ – mmal Sep 11 '15 at 7:04
  • $\begingroup$ Thanks for this, though substituting Integrate for NIntegrate isn't an option for me. See my edited post. $\endgroup$ – mlubin Sep 11 '15 at 14:45
2
$\begingroup$

Try this small variation over your original request

   RegionPlot[
Integrate[PDF[NormalDistribution[0, 1], a], {a, 0, x}] > 0.3 // 
  Evaluate, {x, -10, 10}, {y, -1, 1}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ doesn't work on v9 :( $\endgroup$ – Dr. belisarius Sep 10 '15 at 20:53
  • $\begingroup$ Release .. had to look that up, superseded in version 2 (1991!). (Evaluate works as well ) $\endgroup$ – george2079 Sep 11 '15 at 2:39
  • $\begingroup$ @george2079 oops… I am old and I forget :-) Then against does the job! $\endgroup$ – chris Sep 11 '15 at 5:34
1
$\begingroup$

While the behavior has already been confirmed a bug -- an uncaught Throw from an internal function must always be one, right? -- here are a couple more workarounds.

How I analyze the problem of finding a workaround: From the context NIntegrate`LevinRuleDump` in the error message, one might infer that NIntegrate is trying to determine whether to use (or even using) "LevinRule" for the integration method. Since it's for oscillatory integrands, it doesn't seem that appropriate. Probably it is occurring during the "SymbolicProcessing" phase of NIntegrate. And if not, we could pick an integration rule manually.

RegionPlot[
  NIntegrate[PDF[NormalDistribution[0, 1], a], {a, 0, y}, 
    Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}] >= 
   0.2, {x, -1, 1}, {y, 0.1, 0.7}]; // AbsoluteTiming

Mathematica graphics

RegionPlot[
   NIntegrate[PDF[NormalDistribution[0, 1], a], {a, 0, y}, 
     Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}] >= 
    0.2, {x, -1, 1}, {y, 0.1, 0.7}]; // AbsoluteTiming

Mathematica graphics

Turning off symbolic processing is fastest, but you lose automatic method selection for certain types of integral such as when the integrand is oscillatory. By comparison, ilian's mystery option "NumericalFunction" -> False takes about as long as picking the "GaussKronrodRule" explicitly. For some reason, this method does not produce a spurious error message (which comes from RegionPlot evaluating the argument symbolically).

RegionPlot[
   NIntegrate[PDF[NormalDistribution[0, 1], a], {a, 0, y}] >= 0.2,
   {x, -1, 1}, {y, 0.1, 0.7}, 
   "NumericalFunction" -> False]; // AbsoluteTiming

Mathematica graphics

All produce plots that look the same:

Mathematica graphics

$\endgroup$
  • 1
    $\begingroup$ This is a good approach to analyze the problem (I upvoted). I would like to point out that singularity handling for one dimensional integrals is not switched off by "SymbolicProcessing"->0. That is done with the option setting "SingularityHandler" -> None in the integration strategy methods. $\endgroup$ – Anton Antonov Oct 3 '15 at 14:35
  • 1
    $\begingroup$ @AntonAntonov Thanks for the correction. I hope the update is more accurate. (And thanks for the upvote.) $\endgroup$ – Michael E2 Oct 3 '15 at 23:45

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.