## Numerical Solution to Inverse Integral (Pseudo Random Number Generation)

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$

-
 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! – alext87 Sep 14 2010 at 19:24

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.

-
 Looks good, have copied it to the top. L – Lucas Sep 14 2010 at 22:41

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.

-

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.

-
 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 :) – Lucas Sep 15 2010 at 5:39

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

-

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

-
 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) – J. M. Sep 26 2010 at 23:40