Calculating Lévy stable distributions

There seems rather a lack of code to do this on the web. Here is a very hackish first attempt on my part. It works on the simple principle that, by definition, repeatedly convolving a distribution with itself will eventually give you a Lévy stable distribution.

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.