WOLFRAM

Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 2009

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

Search the Archive

Re: InverseFunction of a CDF

  • To: mathgroup at smc.vnet.net
  • Subject: [mg102807] Re: [mg102767] InverseFunction of a CDF
  • From: "Christian Winzer" <cw420 at cam.ac.uk>
  • Date: Thu, 27 Aug 2009 06:38:33 -0400 (EDT)
  • References: <200908261155.HAA11857@smc.vnet.net> <op.uzahb5ootgfoz2@bobbys-imac-2.local>
Thank you for the quick response!
I have replaced the corresponding lines in the example with the code you sent.
To be on the safe side I have also restarted the kernel before the evaluation.
However, it still does not work :
NSolve returns the rule {{x->xCDF^(-1)[0.5]}} but I can't get a numerical
value for this x.
If I replace y=x/.% that does not do the job and neither if I ask for the
numerical value thereof via 
N[y[[1]]]. 
The plots are also empty..

I have attached the adjusted code including the output below.
Thank you for your help and best wishes!
Christian
____________________________________________________________________________
_____________________________

In[1]:= (* 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])]] &;

(*Calculate unconditional cumulated density xCDF[x,m] of the \
stochastic variable x *)
Clear[xCDF]
xCDF[x_?NumericQ] := 
 NIntegrate[xCDFm[x, m]*mPDF[m], {m, -\[Infinity], \[Infinity]}]
NSolve[xCDF[x] == .5, x]

During evaluation of In[1]:= InverseFunction::ifun: Inverse functions are
being used. Values may be lost for multivalued inverses. >>

During evaluation of In[1]:= InverseFunction::ifun: Inverse functions are
being used. Values may be lost for multivalued inverses. >>

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

In[12]:= y = x /. %

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

In[13]:= N[y[[1]]]

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

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

Out[16]= 1

Out[18]= {306, 2}

In[19]:= Plot[inv1@xCDF@x, {x, -3, 3}]

Out[19]= \!\(\*
GraphicsBox[GraphicsComplexBox[CompressedData["
1:eJxTTMoPSmViYGAwAmIQnbx+wvv//9kPMEDBrNnKk87yIfjz73fkB3Kywfli
F4L9NUVZ4fwTAcEup1VZ4HwPTdUSZSVmON9qXcyy8/pMcP6vrjRhI01GOP/r
9yVT49QZ4PwjSpNM9uX+2w/jW7PWb9ZO+QPnT3EwN3es/QXnJ8+JePpq/g84
//xyOU7GWd/gfDmeDYb7t32B87kUnp1/ueoTnJ8T0jpjxYoPcH7C4fun5Sa9
hfOVbP8FHm18CeevzNj8cPe2p3B+1ueNMd2zHsL5enWskWzKt+F84cor3NNm
X4bzr6fkaK1/fAzOL6/S6js5dR2cv6rfwmGJwVp7GL8vs/FaQtNxOJ/3i4/1
msDLcP6DkmMmZS234fx3O30PlZ94COcva1+YkNL5FM5X/bZauGL2Szj/3ib5
vdYZb+H80OkPOl5P+ADnz1vas9dn2ic4nyU69y1r5Rc4/1zCyuVBLd/g/K28
UrHfc37A+WqzFF+w+/yC83+VtfxdHfwHzn/9S5pvhcU/OH9BevsVaWUGBxh/
guOZHklVRjifL3+102FRJjh/n8uOvdukmOH8v7nxZmF8LHA+I/ulmdp/EXwz
j6JDXgxscP6vO2teiH5C8KH5Ac4HAPdUymg=
"], {}],
AspectRatio->NCache[GoldenRatio^(-1), 0.6180339887498948],
Axes->True,
AxesOrigin->{0, 0},
PlotRange->{{-3, 3}, {0., 1.}},
PlotRangeClipping->True,
PlotRangePadding->{
Scaled[0.02], 
Scaled[0.02]}]\)

In[20]:= 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}]


During evaluation of In[20]:= FindRoot::srect: Value inv1[0.00135044] in
search specification {x$19329,inv1[0.00135044],-3,3} is not a number or
array of numbers. >>

During evaluation of In[20]:= ReplaceAll::reps:
{FindRoot[xCDF[x$19329]-0.00135044,{x$19329,inv1[0.00135044],-3,3}]} is
neither a list of replacement rules nor a valid dispatch table, and so
cannot be used for replacing. >>

During evaluation of In[20]:= FindRoot::srect: Value inv1[0.00135044] in
search specification {x$19329,inv1[0.00135044],-3,3} is not a number or
array of numbers. >>

During evaluation of In[20]:= ReplaceAll::reps:
{FindRoot[xCDF[x$19329]-0.00135044,{x$19329,inv1[0.00135044],-3,3}]} is
neither a list of replacement rules nor a valid dispatch table, and so
cannot be used for replacing. >>

During evaluation of In[20]:= FindRoot::srect: Value inv1[0.00135044] in
search specification {x$19329,inv1[0.00135044],-3.,3.} is not a number or
array of numbers. >>

During evaluation of In[20]:= General::stop: Further output of
FindRoot::srect will be suppressed during this calculation. >>

During evaluation of In[20]:= ReplaceAll::reps: {FindRoot[xCDF[x$19329]-1.
0.00135044,{x$19329,inv1[0.00135044],-3.,3.}]} is neither a list of
replacement rules nor a valid dispatch table, and so cannot be used for
replacing. >>

During evaluation of In[20]:= General::stop: Further output of
ReplaceAll::reps will be suppressed during this calculation. >>

Out[24]= \!\(\*
GraphicsBox[GraphicsComplexBox[CompressedData["
1:eJxTTMoPSmViYGAwAmIQnbx+wvv//9kPMEDBrNnKk87yIfjz73fkB3Kywfli
F4L9NUVZ4fwTAcEup1VZ4HwPTdUSZSVmON9qXcyy8/pMcP6vrjRhI01GOP/r
9yVT49QZ4PwjSpNM9uX+2w/jW7PWb9ZO+QPnT3EwN3es/QXnJ8+JePpq/g84
//xyOU7GWd/gfDmeDYb7t32B87kUnp1/ueoTnJ8T0jpjxYoPcH7C4fun5Sa9
hfOVbP8FHm18CeevzNj8cPe2p3B+1ueNMd2zHsL5enWskWzKt+F84cor3NNm
X4bzr6fkaK1/fAzOL6/S6js5dR2cv6rfwmGJwVp7GL8vs/FaQtNxOJ/3i4/1
msDLcP6DkmMmZS234fx3O30PlZ94COcva1+YkNL5FM5X/bZauGL2Szj/3ib5
vdYZb+H80OkPOl5P+ADnz1vas9dn2ic4nyU69y1r5Rc4/1zCyuVBLd/g/K28
UrHfc37A+WqzFF+w+/yC83+VtfxdHfwHzn/9S5pvhcU/OH9BevsVaWUGBxh/
guOZHklVRjifL3+102FRJjh/n8uOvdukmOH8v7nxZmF8LHA+I/ulmdp/EXwz
j6JDXgxscP6vO2teiH5C8KH5Ac4HAPdUymg=
"], {}],
AspectRatio->NCache[GoldenRatio^(-1), 0.6180339887498948],
Axes->True,
AxesOrigin->{0, 0},
PlotRange->{{-3, 3}, {0., 1.}},
PlotRangeClipping->True,
PlotRangePadding->{
Scaled[0.02], 
Scaled[0.02]}]\)

> -----Original Message-----
> From: DrMajorBob [mailto:btreat1 at austin.rr.com]
> Sent: 26 August 2009 22:44
> To: Christian Winzer; mathgroup at smc.vnet.net
> Subject: Re: [mg102767] InverseFunction of a CDF
> 
> 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: accumulate coefficients of a polynomial
  • Previous by thread: Re: InverseFunction of a CDF
  • Next by thread: Re: Re: InverseFunction of a CDF