3
$\begingroup$

If I have an arbitrary positive monotonically decreasing function $f(x), x \in [0,\infty]$, is there an 'efficient' method for finding $y$ in:

$r = \int\limits_0^y f(x) dx $

for a known $r \in [0, \int\limits_0^\infty f(x) dx]$. By efficient I guess I mean more efficient than doing numerical integration until one lands in within a given distance from $r$.

I particularly care about the cases where the integral of $f$ has no closed form.

Best answer so far (rewrite of below), doubt it gets much better than this:

Start with $s_0 = r$, $y_0 = 0$, then

$y_{i+1} = y_i + \frac{s_i}{f(y_i)} $

$s_{i+1} = s_i - \int\limits_{y_i}^{y_{i+1}}f(x)dx$

Improve this question
$\endgroup$
1
  • $\begingroup$ If you let me know which software packages you have available I may be able to explain how to implement root finding with that software or if you are keen write one! $\endgroup$
    – alext87
    Commented Sep 14, 2010 at 19:24

5 Answers 5

Reset to default
3
$\begingroup$

I guess the only non-trivial thing about the problem is that: $$ x f(0) \geq \int_0^{x} f(t) dt \geq x f(x). $$ So you start by computing the integral $$ r_1 = \int_0^{y_1} f(t) dt,\quad y_1 = \frac{r}{f(0)}. $$ Then replace $r$ by $r - r_1$. I think under reasonable assumptions this should converge pretty quickly (and always lower bound). It should be noted that I only use one of the inqualities, one can probably optimize it by using the other one.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Looks good, have copied it to the top. L $\endgroup$
    – Lucas
    Commented Sep 14, 2010 at 22:41
1
$\begingroup$

There's always the option of expanding $f$ as a series, integrating that, and performing Lagrange inversion to obtain a series for the inverse function, from which you can determine suitable function approximations.

But for the purpose of generating non-uniform pseudorandom numbers, transformation isn't always the best choice; you might have better luck with the rejection approach. I suppose it depends on what distribution you're dealing with.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks JM. I am actually implementing the ziggurat algorithm, the problem with it being that when you get to the tail you have to use something else (or it again). It's looking like the wrong approach :) $\endgroup$
    – Lucas
    Commented Sep 15, 2010 at 5:39
1
$\begingroup$

Luc Devroye has written some timeago a superb book "Non-uniform random variate generation". The whole book is available for free on his webpage.

Improve this answer
$\endgroup$
0
$\begingroup$

Firstly this problem is not well-defined as stated. Let $f(x)=e^{-100x}$ and let $r=f(0)=1$. Since, $\int_0^{\infty}f(x)dx=\frac{1}{100}<1$ there is no $y$ that solves your problem. I am suspecting that you have $f$ so you should convince yourself that the problem actually has a solution.

I always solve this sort of problem by root finding:

Let $H(y)=\int_0^y f(x)dx-r$ then the problem becomes finding a root of $H$. Now there are lots of software packages to solve this, these will give you much better accuracy than trial and error.

In Matlab, the function fzero(fun,x0) will solve this problem.

Improve this answer
$\endgroup$
0
0
$\begingroup$

Newton's method? The derivative should be fairly straightforward to compute...

Improve this answer
$\endgroup$
1
  • $\begingroup$ That would be a good addition to the method I proposed above (using the approximation derived from Lagrange inversion to start up a Newton-Raphson iteration); but sometimes, you have to assess if the increased accuracy in evaluation is worth the computational burden (for generating a whole lot of random variates, you'd want to minimize transcendental function evaluations because they're so slow) $\endgroup$
    – J. M. isn't a mathematician
    Commented Sep 26, 2010 at 23:40

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .