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
- References:
- InverseFunction of a CDF
- From: "Christian Winzer" <cw420@cam.ac.uk>
- InverseFunction of a CDF