A Levy Alpha Distribution Random Generator in Julia

LevyAlphaRand.jl

module LevyAlphaRand

#a module that provides a levy alpha-stable random number in julia
#uses "McCulloch's algorithm", as provided in:
#  Leccardi, M., Comparison of Three Algorithms for Levy Noise Generation
#
#  Note two errors in the McCulloch formula as provided in this paper
#  One:  The phi and w values are swapped
#  Two:  The sin term should have an exponent of one, a fractional exponent would result in complex values

export rande
export randla

function rande(d...)
	- log(ones(d) - rand(d))
end

function randla (d... ;alpha = 1.8, beta = 0)

	#typically levy-alpha random formulae allow for shifts and scales, but these are
	#simple affine transformations of the data so will be omitted.

	#check the range on alpha to make sure it's sensible
	if (alpha < 0.1 || alpha > 2)
		error("levy alpha parameter is outside of defined range")
	end

	res = zeros(d)

	if (alpha == 2)
		#gaussian case, beta has no effect
		res = randn(d)

	elseif (beta != 0)
 		#explicit Algorithm

		#generate a random variable that's uniformly-distributed from [-pi/2, pi/2]
		phi = pi * (rand(d) - 0.5 * ones(d))
		#generate a random variable that's exponentially-distributed
		w = rande(d)

		#generate a correction factor, which is going to be used several times
		c = atan(beta * tan(alpha * pi / 2))

		#levy random formula
		res = sin(alpha * phi + c) .* ((cos((1-alpha) * phi + c)) .^ (1/alpha - 1) ./ w) ./ ((cos(c) * cos(phi))^(1/alpha))

	elseif (alpha == 1)
		#cauchy case (beta == 0)
		res = tan(pi * (rand(d) + 0.5 * ones(d)))

	else
		#McCulloch Algorithm
		#generate a random variable that's uniformly-distributed from [-pi/2, pi/2]
		phi = pi * (rand(d) - 0.5 * ones(d))
		#generate a random variable that's exponentially-distributed
		w = - log(ones(d) - rand(d))

		#note that w and phi are matrices of dimensions specified by d
		res = ((cos((1-alpha) * phi) ./ w) .^ (1/alpha - 1)) .* sin(alpha * phi) ./ (cos(phi) .^ (1/alpha))
	end

	res

end #function

end #module