Take the 2-minute tour ×
MathOverflow is a question and answer site for professional mathematicians. It's 100% free, no registration required.

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$

share|improve this question
    
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 '10 at 19:24
add comment

5 Answers 5

up vote 3 down vote accepted

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.

share|improve this answer
    
Looks good, have copied it to the top. L –  Lucas Sep 14 '10 at 22:41
add comment

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.

share|improve this answer
    
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 '10 at 5:39
add comment

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

share|improve this answer
add comment

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.

share|improve this answer
add comment

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

share|improve this answer
    
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 '10 at 23:40
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.