1
$\begingroup$

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$ ?

$\endgroup$
2
  • $\begingroup$ Suggest that you add the code for p[x] $\endgroup$ Nov 28, 2017 at 11:57
  • $\begingroup$ Done, thanks. This is the simplest case, but still exhibits the problem. $\endgroup$
    – Mashy
    Nov 28, 2017 at 12:01

4 Answers 4

1
$\begingroup$

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]

Mathematica graphics

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]
 ]

Mathematica graphics

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}]
 ]

Mathematica graphics

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.

$\endgroup$
1
$\begingroup$

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)

$\endgroup$
1
$\begingroup$

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.
*)
$\endgroup$
1
$\begingroup$

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]

enter image description here

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

(* 1 *)

Or using the CDF

CDF[dist, x]

enter image description here

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}

enter image description here

To find V corresponding to a Median of α

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

enter image description here

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

Variance[dist] /. solV

enter image description here

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]]

enter image description here

$\endgroup$

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.