Numerical approximation of Levy Flight

I'm trying to produce a computer simulation of a Levy Flight in 2-dimensions; an approximation would be ok. Please excuse the simplistic level of this question: my maths is very rusty.

My proposed method of plotting the Levy Flight is as follows: Start at the origin, and plot a connected series of points, each being in a random direction $\theta$ and a distance f from the previous point, where the size of f conforms to a power law distribution. Wikipedia suggests this should be of the form $y=x^{-a}$ where $1<a<3$.

My problem is that I can't see how to generate values of f using a random number generator: I understand that I can use for example x = rand(1) which will with equal probability (allowing for the limitations of numerical accuracy and pseudo-randomness) take any value between 0 and 1. But here's where my limited math skills start to let me down: to transform this into a suitable value of f, I think I would need to use a 'cumulative distribution function' cdf(x), which maps my equally probable values of x onto some corresponding values of f. I also suspect that cdf(x) would basically be related the integral of $y=x^{-a}$, however I can't see exactly how to do this:

• Do I need to somehow scale $y=x^{-a}$ so the area under it sums to 1 (in order to be a valid probability density function)? I can't see how this is possible, unless I approximate using discrete values of $x$ and/or choosing a limited interval for values of $x$.

• Alternatively, could I do an approximation to this by the following method:(1) select a nominal start value of f, say 1. (2) if rand(1)>t (where t is an arbitrary threshold value) then set f=f*k (say k=1.2) and repeat step 2, otherwise stop and use the current value of f.

If someone could offer some pointers, I would be very grateful. Thank you.

EDIT:

Thank you to Chris for the prompt answer which works perfectly! So for $\alpha = 3$ and $x_{\mathrm{min}} = 1$, when I choose 100,000 random values of u in [0,1] I get this distribution of values of $F^{-1}(u) = x_{\mathrm{min}} (1-u)^{-1/\alpha}$:

... which is just what I'm after (there's a little blip at the end for truncated values of $x>8$).

In fact, since u is evenly distributed I can get exactly the same distribution using $F^{-1}(u) = u^{-1/\alpha}$ which is shorter.

Here's a typical flight plot (for $\alpha=1.5$, 1000 steps):

EDIT: Matlab code for the above example:

alpha=1.5;
s=1000;
x=zeros(1,s);
y=zeros(1,s);

for n=2:s;
theta=rand*2*pi;
f=rand^(-1/alpha);
x(n)=x(n-1)+f*cos(theta);
y(n)=y(n-1)+f*sin(theta);
end;

figure('color', 'white');
axis equal off;
line(x,y);

-
I was trying to figure out the generation and plotting of Lévy Flights. I stumbled upon your question. Looking at the answer of Chris Taylor I am now able to generate the random increments using the inverse cdf method. However, I am still unable to understand that $U[0, 1]$ will produce only positive values that's the Markov chain so produced will be going ahead only and will not step back. I will be thankful to you could explain it a bit. –  Asad Ali Feb 12 '13 at 6:22
Also, I am wondering as to how are you producing a two dimensional plot of Lévy flights because, as I understand it, there is only one sequence of Lévy flights. I will be really, thankful if you could give me the source code (in R or Matlab) of your above example. Cheers, Asad PS: Of course, I am quite new to Lévy Flights. –  Asad Ali Feb 12 '13 at 6:22
Hi @Asad Ali: the power law distribution is used to control the length of each step (always positive), but each step is in a random direction (see second paragraph in my question above), which produces a 2-d path. I will try to find some matlab code for you (but it might take a day or so!). –  Richard Inglis Feb 12 '13 at 10:44
I wanted to thank Richard for his prompt reply and I am keenly waiting for his matlab code of the above example. –  Asad Ali Feb 19 '13 at 4:50

To transform a uniformly distributed random variable into another distribution, you need to use the inverse cumulative distribution function. That is, if $F$ is a cumulative distribution function corresponding to the probability density $f$, and $u$ is a uniform random variable in [0,1] then

$$x = F^{-1}(u)$$

is distributed according to $F$. The cumulative distribution function for a pure power law distribution is

$$F(x) = 1 - \left( \frac{x}{x_{\mathrm{min}}} \right)^{-\alpha}$$

where $x_{\mathrm{min}}$ is the minimum value that your random variable can take, and therefore the inverse distribution function is

$$F^{-1}(u) = x_{\mathrm{min}} (1-u)^{-1/\alpha}$$

If you don't want to introduce an artificial minimum, you can consider two-sided power law distributions. These are not pure power laws (you have to fudge the distribution around zero so that the density is not infinite) but they are asymptotically power laws in the tails.

One such distribution is the Student's t distribution with $\nu$ degrees of freedom, which behaves like a power law with exponent $\nu$ as its values tend to $\pm\infty$.

-
Perfect - thank you! –  Richard Inglis Jul 21 '11 at 13:12
No problem at all. –  Chris Taylor Jul 21 '11 at 13:13
Transformation with the quantile function (i.e. inverse CDF) is good if you have a quantile function available. Otherwise, there are other methods to transform uniform variates appropriately. Devroye's book is a good reference for these things. –  Guess who it is. Jul 21 '11 at 14:34
Not to mention that even if you have the quantile function it may not be fast or stable to compute numerically, and other methods may give better results - a simple result that comes to mind is computing normal random variables by the Box-Muller transormation rather than using the inverse CDF. Thanks for the link J.M., it looks very interesting. –  Chris Taylor Jul 21 '11 at 14:43

$F_{\mathrm{inv}}(r) = x_{\min}\cdot(1 − r)^{−1/(α−1)}$ !

-