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