from numarray import fft from numarray import * def levy(alpha, skew, width, n, center_on_mode=True): x = arange(1,n+1) ** -(alpha+1) x /= sum(x) y = fft.real_fft(x) i = width ** alpha y = (y**((1+skew)*i)) * (conjugate(y)**((1-skew)*i)) x = fft.inverse_real_fft(y) if center_on_mode: mode = argmax(x) x = concatenate((x[mode:],x[:mode])) return concatenate((x[n/2:], x[:n/2])) #Example import pylab pylab.plot(levy(0.5, 0.5, 1024, 65536)) pylab.show()

Use 0 < alpha <= 2, -1 < skew < 1. As FFTs wrap, the tails are inaccurate. If you need accuracy, choose n much larger than width and don't use the tips of the tails. To fit these distributions to data, I've been using numerical optimization to find the parameters for which the liklihood of the data is maximized (Maximum Likelihood estimation).

**Addendum:** Here's another way to calculate the distribution, based on an inverse Fourier transform of the characteristic function:

def levy(alpha, beta, n=10000): """ Note: becomes inaccurate for alpha close to 1.0, beta non-zero. """ factor = log(n) # <- somewhat arbitrary neglog_range = -log((arange(1,n)+0.5) / float(n)) t = neglog_range * factor if alpha == 1.0: F = -2/pi * log(t) else: F = tan(pi/2 * alpha) characteristic = -(t ** alpha) * (1-complex(0.0,beta)*F) mag = exp(characteristic.real + log(factor) + neglog_range) angle = characteristic.imag return lambda x: sum(mag * cos(asarray(x)[...,NewAxis]*t + angle), -1) / (pi*n) #Example import pylab pylab.plot(levy(1.1,1.0)( arange(-100,100)*0.1 )) pylab.show()

This trick is in choosing good sample points for the numerical inverse Fourier transform. The above seems to do pretty well. Some problems for alpha close to one and beta non-zero, and for very small alpha.