-
-
Save gragusa/53a316a0c0e94a879046 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Utilities needed in rand() | |
# (exp(x) - 1) / x | |
function d2(z::Real) | |
const p1 = 0.840066852536483239e3 | |
const p2 = 0.200011141589964569e2 | |
const q1 = 0.168013370507926648e4 | |
const q2 = 0.18013370407390023e3 | |
const q3 = 1.0 | |
if abs(z) > 0.1 | |
return (exp(z) - 1) / z | |
else | |
zz = z * z | |
pv = p1 + zz * p2 | |
return 2 * pv / (q1 + zz * (q2 + zz * q3) - z * pv) | |
end | |
end | |
# tan(x) / x | |
function tan2(xarg::Real) | |
const p0 = 0.129221035e3 | |
const p1 = -0.887662377e1 | |
const p2 = 0.528644456e-1 | |
const q0 = 0.164529332e3 | |
const q1 = -0.451320561e2 | |
const q2 = 1.0 | |
const piby4 = π / 4 | |
x = abs(xarg) | |
if x > piby4 | |
return tan(xarg) / xarg | |
else | |
x /= piby4 | |
xx = x * x | |
return (p0 + xx * (p1 + xx * p2)) / | |
(piby4 * (q0 + xx * (q1 + xx * q2))) | |
end | |
end | |
immutable Stable <: ContinuousUnivariateDistribution | |
alpha::Float64 # Characteristic exponent | |
beta::Float64 # Skewness in standard parameterization | |
bprime::Float64 # Skewness in revised parameterization | |
gamma::Float64 # Scale | |
delta::Float64 # Location | |
function Stable(a::Real, b::Real, g::Real, d::Real) | |
0 < a <= 2 || throw(ArgumentError("alpha must be in (0, 2]")) | |
abs(b) <= 1 || throw(ArgumentError("beta must be in [-1, 1]")) | |
g > 0 || throw(ArgumentError("gamma must be non-negative")) | |
phi0 = pi / 2 * b * ((1 - abs(1 - a)) / a) # Always zero for a = 2 | |
bprime = a == 1.0 ? b : -tan(pi / 2 * (1 - a)) * tan(a * phi0) | |
new(float64(a), float64(b), float64(bprime), float64(g), float64(d)) | |
end | |
end | |
# 2-parameter stable distribution | |
Stable(a::Real, b::Real) = Stable(a, b, 1.0, 0.0) | |
# Algorithm from Mallows | |
function Distributions.rand(d::Stable) | |
const piby2 = π / 2 | |
α, β, γ, δ = d.alpha, d.bprime, d.gamma, d.delta | |
u = rand() | |
w = rand(Exponential()) | |
ε = 1 - α | |
# Compute some tangents | |
phiby2 = piby2 * (u - 0.5) | |
a = phiby2 * tan2(phiby2) | |
bb = tan2(ε * phiby2) | |
b = ε * phiby2 * bb | |
if ε > -0.99 | |
τ = β / (tan2(ε * piby2) * piby2) | |
else | |
τ = β * piby2 * ε * (1 - ε) * tan2((1 - ε) * piby2) | |
end | |
# Compute some necessary subexpressions | |
da = a^2 | |
db = b^2 | |
a2 = 1.0 - da | |
a2p = 1.0 + da | |
b2 = 1.0 - db | |
b2p = 1.0 + db | |
# compute coefficient | |
z = a2p * (b2 + 2 * phiby2 * bb * τ) / (w * a2 * b2p) | |
# compute the exponential-type expression | |
logz = log(z) | |
d = d2(ε * logz / (1 - ε)) * (logz / (1 - ε)) | |
# compute "standard" stable with std = sqrt(2) | |
s = (1 + ε * d) * 2 * | |
((a - b) * (1 + a * b) - phiby2 * τ * bb * (b * a2 - 2 * a)) | |
s /= (a2 * b2p) | |
s += τ * d | |
# shift scale and location | |
s = s / sqrt(2) * γ + δ | |
return s | |
end |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment