MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: InverseFunction of a CDF

  • To: mathgroup at smc.vnet.net
  • Subject: [mg102798] Re: [mg102767] InverseFunction of a CDF
  • From: DrMajorBob <btreat1 at austin.rr.com>
  • Date: Thu, 27 Aug 2009 06:36:49 -0400 (EDT)
  • References: <200908261155.HAA11857@smc.vnet.net>
  • Reply-to: drmajorbob at yahoo.com

Inputs 78 and 80 fail because Mathematica simply doesn't know the  
InverseFunction or Quantile. It never computes them; they're built in as  
definitions, for standard functions and distributions.

Input 82 fails, in part, because Mathematica attempts a symbolic solution,  
but x must be numeric in order to evaluate xCDF.

Eliminating THAT problem gets us back to the InverseFunction problem:

Clear[xCDF]
xCDF[x_?NumericQ] :=
  NIntegrate[xCDFm[x, m]*mPDF[m], {m, -\[Infinity], \[Infinity]}]

NSolve[xCDF[x] == 0.5, x]

InverseFunction::ifun: Inverse functions are being used. Values may be  
lost for multivalued inverses. >>

InverseFunction::ifun: Inverse functions are being used. Values may be  
lost for multivalued inverses. >>

{{x -> \!\(\*
TagBox[
SuperscriptBox["xCDF",
RowBox[{"(",
RowBox[{"-", "1"}], ")"}]],
{InverseFunction, 1, 1},
Editable->False]\)[0.5]}}

The plot helps a great deal, however:

p = Plot[xCDF[x], {x, -3, 3}];
lines = Cases[p, Line[a_List] :> a, Infinity];
Length@lines
pts = Reverse /@ lines[[1]];
Dimensions@pts

1

{306, 2}

I just verified that the plot contains a single line with 306 data points,  
and I reversed x and y to form "pts". Now to interpolate an inverse:

Checking with a plot, we get a very nice Identity function:

Plot[inv1@xCDF@x, {x, -3, 3}]

If more precision is needed, FindRoot can be used with inv1 as a first  
guess:

Clear[inv2]
inv2[y_] /; y <= 0 := -3
inv2[y_] /; y >= 1 := 3
inv2[y_?NumericQ] :=
  Module[{x}, x /. FindRoot[xCDF[x] - y, {x, inv1@y, -3, 3}]]

Plot[inv2@xCDF@x, {x, -3, 3}]

You could use different limits, of course, in place of -3 and +3.

Bobby

On Wed, 26 Aug 2009 06:55:06 -0500, Christian Winzer <cw420 at cam.ac.uk>  
wrote:

> Hello,
>
> I am calculating the cumulated density function of a variable x based on  
> the
> variables m and z using the following code:
>
> (* Distribution of the stochastic variable m*)
> M = NormalDistribution[0, 1];
> mPDF = PDF[M];
> mCDF = CDF[M];
>
> (* Weight of m in x*)
> w = 0.5;
>
> (* Distribution of the stochastic variable z*)
> Z = NormalDistribution[0, 1];
> zPDF = PDF[Z];
> zCDF = CDF[Z];
>
> (*This section calculates the distribution of x = w * m + z *)
> (*Conditional cumulated density xCDFm[x,m] of the stochastic variable x
> conditional on m *)
> xCDFm = Evaluate[zCDF[(#1 - w * #2)/(Sqrt[1 - w^2])]] &;
>
> (*Unconditional cumulated density xCDF[x] of the stochastic variable x *)
> xCDF = NIntegrate[xCDFm[#1, m]* mPDF[m], {m, -\[Infinity], \[Infinity]}]  
> &;
>
>
> The calculation seems to work, because I can plot xCDF correctly.
> But for some reason, I can not get the inverse function. I have tried 3
> different ways:
>
>
> In[78]:=
> x = InverseFunction[xCDF][#1] &;
> Evaluate[x[0.5]]
>
> Out[79]=
> InverseFunction[NIntegrate[xCDFm[#1, m] mPDF[m], {m, -\[Infinity],
> \[Infinity]}] &][0.5]
>
>
> In[80]:=
> Evaluate[Quantile[xCDF, 0.5]]
>
> Out[81]=
> Quantile[NIntegrate[xCDFm[#1, m] mPDF[m], {m, -\[Infinity],  
> \[Infinity]}] &,
> 0.5]
>
>
> In[82]:=
> Solve[xCDF[x] == 0.5, x]
>
> During evaluation of In[73]:= NIntegrate::inumr: The integrand  
> (E^(-(m^2/2))
> (1+Erf[0.816497 (-0.5 m+((xCDF^(-1))[Slot[<<1>>]]&))]))/(2 Sqrt[2 \[Pi]])
> has evaluated to non-numerical values for all sampling points in the  
> region
> with boundaries {{-\[Infinity],0.}}. >>
>
> During evaluation of In[73]:= NIntegrate::inumr: The integrand  
> (E^(-(m^2/2))
> (1+Erf[0.816497 (-0.5 m+((xCDF^(-1))[Slot[<<1>>]]&))]))/(2 Sqrt[2 \[Pi]])
> has evaluated to non-numerical values for all sampling points in the  
> region
> with boundaries {{-\[Infinity],0.}}. >>
>
> During evaluation of In[73]:= NIntegrate::inumr: The integrand  
> (E^(-(m^2/2))
> (1+Erf[0.816497 (-0.5 m+((xCDF^(-1))[Slot[<<1>>]]&))]))/(2 Sqrt[2 \[Pi]])
> has evaluated to non-numerical values for all sampling points in the  
> region
> with boundaries {{-\[Infinity],0.}}. >>
>
> During evaluation of In[73]:= General::stop: Further output of
> NIntegrate::inumr will be suppressed during this calculation. >>
>
> During evaluation of In[73]:= Solve::ifun: Inverse functions are being  
> used
> by Solve, so some solutions may not be found; use Reduce for complete
> solution information. >>
>
> Out[83]= {}
>
> In this particular example I know without calculation, that xCDF is the
> CDF[NormalDistribution[0,Sqrt[1.25]]], however I want to use different
> distributions for m and z in the future, and then I can't always tell the
> distribution of x, so I will need to calculate it. Any ideas how I could  
> do
> that?
>
> Thank you very much in advance and best wishes,
> Christian
>
>
>
>
>



-- 
DrMajorBob at yahoo.com


  • Prev by Date: Re: Re: Viewing packages in mathematica
  • Next by Date: Re: Re: Viewing packages in mathematica
  • Previous by thread: InverseFunction of a CDF
  • Next by thread: Re: InverseFunction of a CDF