Efficient calculation of inverse CDF function for non-analytical CDF
Roman Werpachowski
I am tackling a problem with inverting a cumulative distribution function
(CDF) of a random variable which has an analytical probability density
function (PDF) on (-infty, infty), but no analytical CDF, which hence has
to be calculated numerically. The naive approach in Matlab goes like that
(symbolically)
cdf = quad( pdf );
icdf(p) = fzero( cdf(x) - p )
It is obviously very slow. Any suggestions at how this could be implemented
faster? Tabularizing the inverse CDF is not an option, because it is
defined by several free parameters and I need to be able to invert it for
any of them.
Regards,
Roman Werpachowski
Herman Rubin
>I am tackling a problem with inverting a cumulative distribution function
>(CDF) of a random variable which has an analytical probability density
>function (PDF) on (-infty, infty), but no analytical CDF, which hence has
>to be calculated numerically. The naive approach in Matlab goes like that
>(symbolically)
Your terminology is bad. If the density is real analytical,
so is the cdf. If you mean no closed form, say it.
>cdf = quad( pdf );
>icdf(p) = fzero( cdf(x) - p )
>It is obviously very slow. Any suggestions at how this could be implemented
>faster? Tabularizing the inverse CDF is not an option, because it is
>defined by several free parameters and I need to be able to invert it for
>any of them.
One can use higher order iteration methods. For example,
solving x+q*exp(x) = b can be done using a fifth-order
method with only three starting cases, and one can restrict
q to be +- 1.
--
This address is for information only. I do not claim that these views
are those of the Statistics Department or of Purdue University.
Herman Rubin, Department of Statistics, Purdue University
hru...@stat.purdue.edu Phone: (765)494-6054 FAX: (765)494-0558
Roman Werpachowski
news:h2rj27$3m...@odds.stat.purdue.edu:
> In article <h2qstg$mef$1...@inews.gazeta.pl>,
> Roman Werpachowski <roman_dot_werpachowski@gmail_dot_com.org> wrote:
>>Hello,
>
>>I am tackling a problem with inverting a cumulative distribution
>>function (CDF) of a random variable which has an analytical
>>probability density function (PDF) on (-infty, infty), but no
>>analytical CDF, which hence has to be calculated numerically. The
>>naive approach in Matlab goes like that (symbolically)
>
> Your terminology is bad. If the density is real analytical,
> so is the cdf. If you mean no closed form, say it.
Sorry. Now I see you're a genuine mathematician ;-)
>
> One can use higher order iteration methods. For example,
> solving x+q*exp(x) = b can be done using a fifth-order
> method with only three starting cases, and one can restrict
> q to be +- 1.
That would require the use of higher-order derivatives of the inverted CDF,
right (as in going from Newton's to Halley's method)? What about numerical
stability, is it going to be a problem with higher-order methods?
Regards,
Roman Werpachowski
Axel Vogt
You need the derivatives of CDF (not of the inverse). For the tails
one may try to switch to logarithms (and/or use asymptotics for CDF
to be inverted)
Herman Rubin
Roman Werpachowski <roman_dot_werpachowski@gmail_dot_com.org> wrote:
>hru...@odds.stat.purdue.edu (Herman Rubin) wrote in
>news:h2rj27$3m...@odds.stat.purdue.edu:
>> In article <h2qstg$mef$1...@inews.gazeta.pl>,
>> Roman Werpachowski <roman_dot_werpachowski@gmail_dot_com.org> wrote:
>>>Hello,
>>>I am tackling a problem with inverting a cumulative distribution
>>>function (CDF) of a random variable which has an analytical
>>>probability density function (PDF) on (-infty, infty), but no
>>>analytical CDF, which hence has to be calculated numerically. The
>>>naive approach in Matlab goes like that (symbolically)
>> Your terminology is bad. If the density is real analytical,
>> so is the cdf. If you mean no closed form, say it.
>Sorry. Now I see you're a genuine mathematician ;-)
>> One can use higher order iteration methods. For example,
>> solving x+q*exp(x) = b can be done using a fifth-order
>> method with only three starting cases, and one can restrict
>> q to be +- 1.
>That would require the use of higher-order derivatives of the inverted CDF,
>right (as in going from Newton's to Halley's method)? What about numerical
>stability, is it going to be a problem with higher-order methods?
I was thinking of methods using rational approximations,
and using square root as a primitive operator. This
changes the game quite a bit.
There is little problem with numerical stability, unless
both the first and second derivatives are small, or the
derivatives are approximated numerically. If the density
has a closed form, so do its derivatives.
Roman Werpachowski
news:7bdpdoF...@mid.individual.net:
Yes, sorry for the mistake. The problem with the tails is that I may find
myself inverting the CDF for very low probabilities, and I want to control
the error there as well.
RW
Roman Werpachowski
> In article <h2s5ht$rt7$1...@inews.gazeta.pl>,
> Roman Werpachowski <roman_dot_werpachowski@gmail_dot_com.org> wrote:
>>That would require the use of higher-order derivatives of the inverted
>>CDF, right (as in going from Newton's to Halley's method)? What about
>>numerical stability, is it going to be a problem with higher-order
>>methods?
>
> I was thinking of methods using rational approximations,
> and using square root as a primitive operator. This
> changes the game quite a bit.
Do you mean something like this:
http://en.wikipedia.org/wiki/Simple_rational_approximation
?
>
> There is little problem with numerical stability, unless
> both the first and second derivatives are small, or the
> derivatives are approximated numerically. If the density
> has a closed form, so do its derivatives.
The PDF is composed of exp and Bessel function, so the computation of its
derivatives may be more expensive and depends on how accurately the Bessel
functions are calculated.
Regards,
Roman Werpachowski
Axel Vogt
May be you want to state it concrete.