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