NSolve with NIntegrate returns empty list

Question

I have defined a function $p\left[x,V\right]$

p[x_,V_]:= (2/Sqrt[Pi V^2]) Exp[-x^2 / V^2]

just a Gaussian, normalised such that integrating over all positive x I get 1.

I have functions

cdf[V_,lim_]:= NIntegrate[p[x,V],{x,0,lim}]
median[V_]:=FindRoot[cdf[V,lim]==0.5,{lim,1.}]

which calculates the median of $p[x,V]$.

I would now like to write a function which takes a number $\alpha$ as input, and returns the variance $V$ for which $median[V] = \alpha$. I do this as follows

fun[α_]:=NSolve[median[V]==α,V]

This function however only returns empty lists $\{ \}$, even when I can verify that a solution exists (by e.g. plotting median[V], or "manually" finding the correct $V$).

How should I rewrite my functions to find the $V$ such that $median[V] = \alpha$ ?

Answer 1

The definition for the probability density function is:

p[x_, variance_] := 2/Sqrt[π variance^2] Exp[-x^2/variance^2]

Here is a plot with different variance parameters.

Plot[{p[x, 0.2], p[x, 0.4], p[x, 0.6], p[x, 0.8], p[x, 1.0]},
  {x, -2, 2}, PlotRange -> All]

Here is the definition of the cumulative distribution function. We force the input arguments to be numerical to facilitate subsequent processing in FindRoot (avoids error messages).

cdf[variance_ /; variance > 0, limit_ /; limit > 0] := 
 NIntegrate[p[x, variance], {x, 0, limit}]

We are interested in the correlation between the variance and the x value corresponding to cdf = 0.5.

Here is a graph of it with different variance parameters.

Show[
 Plot[{cdf[0.2, x], cdf[0.4, x], cdf[0.6, x], cdf[0.8, x], 
   cdf[1.0, x]}, {x, 0, 1}, PlotRange -> {0, 1}],
 Plot[0.5, {x, 0, 2}, PlotStyle -> Black]
 ]

Here is a table of variance values and the corresponding x value where the cdf is 0.5.

median[variance_] := FindRoot[cdf[variance, limit] == 0.5, {limit, 0.2}]

varianceData = 
 Table[{limit /. median[variance], variance}, {variance, 0.2, 1, 0.2}]
(* {{0.0953873, 0.2}, {0.190775, 0.4}, {0.286162, 
  0.6}, {0.381549, 0.8}, {0.476936, 1.}} *)

We superimpose it on our plot.

Show[
 Plot[{cdf[0.2, x], cdf[0.4, x], cdf[0.6, x], cdf[0.8, x], 
   cdf[1.0, x]}, {x, 0, 1}, PlotRange -> {0, 1}],
 Plot[0.5, {x, 0, 2}, PlotStyle -> Black],
 ListPlot[
  Transpose@Join[{varianceData[[All, 1]]}, {ConstantArray[0.5, 5]}],
  PlotStyle -> {PointSize[0.025], Red}]
 ]

The question is, given a particular x = α, what is the corresponding variance?

That can be solved using the same procedure as median but swapping the variable and parameter.

fun[α_] := FindRoot[cdf[variance, α] == 0.5, {variance, 0.2}]

Let's test it by using as inputs the limit value from varianceData and see if we get back the same variances.

limitData = 
 Table[{limit, variance /. fun[limit]},
   {limit, varianceData[[All, 1]]}]

(* {{0.0953873, 0.2}, {0.190775, 0.4}, {0.286162, 
  0.6}, {0.381549, -0.8}, {0.476936, -1.}} *)

Sure enough, seems to work fine.

Answer 2

Sorry,as a novice I'm no allowed to comment your question:

The definition of your median-function seems to be wrong, because the parameter lim is not defined!

Perhaps

median[V_?NumericQ] := FindRoot[cdf[V, lim] == 0.5, {lim, 1.}]

is the function you're looking for.

In this definition you find also the appendix _?NumericQ, which helps to handle purely numerical functions like the first one(NIntegrate)

Answer 3

Aside from _?NumericQ to quiet some unimportant error, the big error is the missing lim /. ... in median[]. FindRoot returns a Rule, which needs to be converted to a numeric value.

ClearAll[p, cdf, fun, median];
p[x_, V_] := (2/Sqrt[Pi V^2]) Exp[-x^2/V^2]
cdf[V_?NumericQ, lim_?NumericQ] := NIntegrate[p[x, V], {x, 0, lim}]
median[V_?NumericQ] := lim /. FindRoot[cdf[V, lim] == 0.5, {lim, 1.}]
fun[α_] := NSolve[median[V] == α, V]

Checks: Left inverse.

fun[0.5`]
median /@ (V /. %)

NSolve::ifun: Inverse functions are being used by NSolve, so some solutions may not be found; use Reduce for complete solution information.

(*  {{V -> 1.04836}}  <-- fun[0.5] *)

(*  {0.5}             <-- {median[V /. fun[0.5]]} *)

Right inverse.

median[1.]
fun[%]

(*  0.476936     <-- median[1.] *)

NSolve::ifun: Inverse functions are being used by NSolve, so some solutions may not be found; use Reduce for complete solution information.

(*  {{V -> 1.}}  <-- fun[median[1.]] *)

Note the use of InverseFunction. Basically, NSolve, which is principally for analytic equations and not ones involving numerical methods, calls InverseFunction[median, 1, 1][α]. So another way to define fun is

fun[α_] := InverseFunction[median, 1, 1][α]

Note the original fun returns a Rule, but this returns a value:

median[2.]
fun[%]
(*
  0.953873
  2.
*)

Answer 4

To define your own probability distribution and make the most use of Mathematica's built-in capabilities, use ProbabilityDistribution. Note the change from V^2 to V to be proportional to the actual Variance of your distribution.

dist = ProbabilityDistribution[2/Sqrt[π V] Exp[-x^2/V], {x, 0, Infinity}, 
   Assumptions -> V > 0];

The DistributionParameterAssumptions are the Assumptions specified in ProbabilityDistribution.

assume = DistributionParameterAssumptions[dist]

(* V > 0 *)

PDF[dist, x]

Assuming[assume, Integrate[PDF[dist, x], {x, 0, Infinity}]]

(* 1 *)

Or using the CDF

CDF[dist, x]

Assuming[assume, Limit[CDF[dist, x], x -> Infinity]]

(* 1 *)

Having used ProbabilityDistribution, built-in functions such as Mean, Median, and Variance will work with your distribution. Note that the Variance is proportional to, but not equal to V

#[dist] & /@ {Mean, Median, Variance}

To find V corresponding to a Median of α

solV = Solve[Median[dist] == α, V][[1]]

However if you want the Variance corresponding to a Median of α

Variance[dist] /. solV

Plot[Evaluate[{V, Variance[dist]} /. solV],
 {α, 0, 5}, AxesLabel -> {Style[α, Bold, 14], None},
 PlotLegends -> Placed[{V, Variance}, {.85, .4}],
 PlotLabel -> Style["V and Variance for Median[dist] == α", 14, Bold]]