# Calculate median of pdf using NSolve [closed]

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.

• 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].
– JimB
Oct 15, 2021 at 22:28
• @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. Oct 15, 2021 at 23:34
• 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. Oct 16, 2021 at 0:08
• @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$.
– JimB
Oct 16, 2021 at 5:00
• @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).
– JimB
Oct 16, 2021 at 17:10

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

• 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." Oct 16, 2021 at 4:01
• 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? Oct 16, 2021 at 10:54