Re: InverseFunction of a CDF
- To: mathgroup at smc.vnet.net
- Subject: [mg102816] Re: [mg102767] InverseFunction of a CDF
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Fri, 28 Aug 2009 00:44:54 -0400 (EDT)
- References: <200908261155.HAA11857@smc.vnet.net>
- Reply-to: drmajorbob at yahoo.com
Is I noted already, the code
NSolve[xCDF[x] == .5, x]
does not work. It gets you to an undefined InverseFunction. Adding
y = x /. %
afterward doesn't help, and neither does
N[y[[1]]]
Then you plotted inv1 and defined inv2, which depends on inv1, without
DEFINING inv1.
And OOPS, my bad, I left out that definition! So here's the code that
works:
M = NormalDistribution[0, 1];
mPDF = PDF[M];
mCDF = CDF[M];
w = 0.5;
Z = NormalDistribution[0, 1];
zPDF = PDF[Z];
zCDF = CDF[Z];
xCDFm = Evaluate[zCDF[(#1 - w*#2)/(Sqrt[1 - w^2])]] &;
Clear[xCDF]
xCDF[x_?NumericQ] :=
NIntegrate[xCDFm[x, m]*mPDF[m], {m, -\[Infinity], \[Infinity]}]
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}
inv1 = Interpolation@pts;
Plot[inv1@xCDF@x, {x, -3, 3}]
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}]
inv1 and inv2 are both decent inverses to xCDF.
Bobby
On Thu, 27 Aug 2009 03:34:42 -0500, Christian Winzer <cw420 at cam.ac.uk>
wrote:
> 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
>
--
DrMajorBob at yahoo.com
- References:
- InverseFunction of a CDF
- From: "Christian Winzer" <cw420@cam.ac.uk>
- InverseFunction of a CDF