I am using the Chambers-Mallows-Stuck method for simulation of skewed stable random variables http://www.sciencedirect.com/science/article/pii/0167715295001131.
There are two equations, one for $\alpha=1$ and one for leftover values of $\alpha$.
I would expect that both equation will return similar output for similar values of $\alpha$, but it is not.
I double checked the code against the paper. And I have not found any problem.
I attached the resulting histograms and the Python code at bottom.
According to the: https://en.wikipedia.org/wiki/Stable_distribution#/media/File:Levy_distributionPDF.png
It seems that my implementation (or algorithm by itself) is really wrong. Can somebody verify if the described equations in paper makes sense? (eq 3.8 and 3.9)
def levy_noise(n, alpha=2., beta=1., sigma=1., position=0.):
# correct the inputs or throw error
alpha = float(alpha)
beta = float(beta)
check_type_or_raise(n, int, "n")
if not 0. <= alpha <= 2.:
raise ValueError("Alpha must be between 0 and 2")
if not -1. <= beta <= 1.:
raise ValueError("Beta must be between -1 and 1")
# generate random series
v = np.random.uniform(-0.5*np.pi, 0.5*np.pi, n)
w = np.random.exponential(1, n)
if alpha == 1:
arg1 = (0.5 * np.pi) + (beta * v)
arg2 = w * np.cos(v)
arg3 = 0.5 * np.pi * beta * sigma + np.log10(sigma)
x = 0.5 * np.pi * (arg1 * np.tan(v) - beta * np.log10(arg2 / arg1))
return (sigma * x) + arg3 + position
else:
arg1 = 0.5*np.pi*alpha
b_ab = np.arctan(beta*np.tan(arg1)) / alpha
s_ab = (1 + (beta**2) * np.tan(arg1)**2)**(1/(2.*alpha))
arg2 = alpha * (v + b_ab)
n1 = np.sin(arg2)
d1 = np.cos(v)**(1/alpha)
n2 = np.cos(v - arg2)
x = s_ab * (n1/d1) * (n2/w)**((1-alpha)/alpha)
return (sigma * x) + position