Numerical approximation of Levy Flight

Question

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);

Answer 1

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$.

Answer 2

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