4
$\begingroup$

I was trying to compute the median of a pdf using NSolve, defined such that $$ \int_{0}^{\rm median}p(x)dx = \frac{1}{2} $$ but NSolve seems unhappy about solving for an upper limit of an integral. A simple example:

F[x_] := 2*Pi^(-1/2)*NIntegrate[Exp[-y^2], {y, 0, x}]
NSolve[F[x] == 0.5, x]

This gives the warning

NIntegrate::nlim: y = x is not a valid limit of integration.

I find this surprising - what is causing NSolve to not be able to handle this and what is a workaround? Obviously in this example Mathematica can symbolically integrate to give ${\rm erf}$ and then reports

NSolve[Erf[x] == 0.5]

During evaluation of In[30]:= NSolve::ifun: Inverse functions are being used by NSolve, so some solutions may not be found; use Reduce for complete solution information.

{{x -> 0.476936}}

but symbolically integrating and then using inverse functions seems huge overkill for what could in this case easily be achieved by a few function calls along with interval bisection.

$\endgroup$
9
  • $\begingroup$ What about something like the following: dist = ProbabilityDistribution[2 \[Pi]^(-1/2) Exp[-y^2], {y, 0, \[Infinity]}]; NSolve[CDF[dist, y] == 1/2, y]. $\endgroup$
    – JimB
    Oct 15, 2021 at 22:28
  • $\begingroup$ @JimB yes that works in this case. I'd still really like to know what NSolve is doing since I can't seem to get it to play nicely with NIntegrate in a variety of contexts. $\endgroup$
    – jacob1729
    Oct 15, 2021 at 23:34
  • 2
    $\begingroup$ Since the definition of F uses a numeric technique (i.e., NIntegrate), the argument to F should be restricted to numeric values: F[x_?NumericQ] := ... You will then get the expected results. $\endgroup$
    – Bob Hanlon
    Oct 16, 2021 at 0:08
  • 1
    $\begingroup$ @user64494 Not true. The pdf integrates to 1 over 0 to $\infty$. It is not a normal distribution but rather a half-normal distribution with parameter $\sqrt\pi$. $\endgroup$
    – JimB
    Oct 16, 2021 at 5:00
  • 1
    $\begingroup$ @user64494 I get {{y -> 0.225312}} (Mathematica 12.3.1.0, Windows 10). And CDF[dist, 0.2253120550121781] gets me 0.25 (the correct value). $\endgroup$
    – JimB
    Oct 16, 2021 at 17:10

1 Answer 1

5
$\begingroup$

We can replace NIntegrate by Integrate and use FindRoot.

F[x_] := 2*Pi^(-1/2)*Integrate[Exp[-y^2], {y, 0, x}]
FindRoot[F[x] == 0.5, {x, 0}]

{x -> 0.476936}

Another way is again use NDSolve.

NDSolve[{F'[x] == 2*Pi^(-1/2)*Exp[-x^2], F[0] == 0, 
  WhenEvent[F[x] == 0.5, {Print[x], "StopIntegration"}]}, 
 F[x], {x, 0, ∞}]

0.476936

$\endgroup$
2
  • $\begingroup$ F[x_?NumericQ] := 2*Pi^(-1/2)*NIntegrate[Exp[-y^2], {y, 0, x}]; FindRoot[F[x] == 0.5, {x, 1}] also works in 12.3.1, producing {x -> 0.476936} and "NIntegrate::nlim: y = x is not a valid limit of integration." $\endgroup$
    – user64494
    Oct 16, 2021 at 4:01
  • $\begingroup$ I agree with some other users that replacing F[x_?NumericQ] in my original suggestion (ie using NIntegrate and NSolve) seems to work in this case. Is there any difference between using NSolve and FindRoot? $\endgroup$
    – jacob1729
    Oct 16, 2021 at 10:54

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