Hi, does anyone know if there is an implementation in any package of the skew normal distribution? There doesn’t appear to be one in Distributions.jl. Thank you!
This article discusses an implementation of Owen’s T function, needed for the skew normal, and provides Fortran code.
Unfortunately, I see no mention of a license agreement. So actually making any use of the code would be difficult.
Are there open source implementations of Owen’s T somewhere?
EDIT:
Good news:
Code distributed with JSS articles uses the GNU General Public License (GPL) version 2 or version 3 or a GPL-compatible license. JSS does NOT consider software distributed under other licenses.
Information on Mission | Journal of Statistical Software
While most of Julia is MIT, the statistical functions from Distributions.jl are GPL anyway.
That would cover the cdf, and SpecialFunctions.jl would cover what’s needed for the pdf.
What features do you need?
Thanks for looking into it! Basically I’m just looking to draw from the distribution, having specified the relevant parameters. For example:
mydist = Normal(0, 1)
rand(mydist)
Ideally I would love to be able to just plug in SkewNormal for Normal. But I guess all I really need is the cdf and the already-existing ability to draw from Uniform(0,1), correct?
Since computing the CDF is tricky, I would suggest using rejection sampling. Something like the following should work:
using StatsFuns
function rand_skewnormal(alpha)
while true
z = randn()
u = rand()
if u < StatsFuns.normcdf(z*alpha)
return z
end
end
end
2 Likes
You could also try the sample() function from ApproxFun, click here.
1 Like
Do you need this distribution specifically, or do you just want to model skewness?
What methods do you need, (log)pdf, cdf, quantiles, or random draws?
Depending on the answer to this, if you don’t need skew normal, there may be an easier solution (using various transformations of the uniform or the standard normal).
To sample from a skew normal, proposition 2.1 in Castillo et al. 2011 might be of some help.
Here’s a function I wrote a while ago (in julia 0.5, so it might need dusting):
immutable SkewNormalDP{T <: Real} <: ContinuousUnivariateDistribution
ξ::T
ω::T
α::T
SkewNormalDP(ξ, ω, α) = new(ξ, ω, α)#(@check_args(SkewNormalDP, ω > zero(ω)); new(ξ, ω, α))
end
#### Outer constructors
SkewNormalDP{T<:Real}(ξ::T, ω::T, α::T) = SkewNormalDP{T}(ξ, ω, α)
SkewNormalDP{T<:Real}(α::T) = SkewNormalDP{T}(zero(α), one(α), α)
# #### Conversions
#=convert{T <: Real, S <: Real}(::Type{Normal{T}}, μ::S, σ::S) = Normal(T(μ), T(σ))=#
#=convert{T <: Real, S <: Real}(::Type{Normal{T}}, d::Normal{S}) = Normal(T(d.μ), T(d.σ))=#
@distr_support SkewNormalDP -Inf Inf
#### Parameters
params(d::SkewNormalDP) = (d.ξ, d.ω, d.α)
#### Statistics
mean_z(d::SkewNormalDP) = √(2/π) * delta(d)
std_z(d::SkewNormalDP) = 1 - 2/π * delta(d)^2
delta(d::SkewNormalDP) = d.α/√(1+d.α^2)
mean(d::SkewNormalDP) = d.ξ + d.ω * mean_z(d)
var(d::SkewNormalDP) = abs2(d.ω)*(1-mean_z(d)^2)
std(d::SkewNormalDP) = √var(d)
skewness(d::SkewNormalDP) = (4-π)/2 * mean_z(d)^3 / (1-mean_z(d)^2)^(3/2)
pdf(d::SkewNormalDP, x::Real) = 2/d.ω*normpdf((x-d.ξ)/d.ω)*normcdf(d.α*(x-d.ξ)/d.ω)
logpdf(d::SkewNormalDP, x::Real) = log(2)-log(d.ω)+normlogpdf((x-d.ξ)/d.ω)+normlogcdf(d.α*(x-d.ξ)/d.ω)
function rand(d::SkewNormalDP)
u0 = randn()
v = randn()
δ = delta(d)
u1 = δ * u0 + √(1-δ^2) * v
return d.ξ + d.ω * sign(u0) * u1
end
immutable SkewNormalCP{T <: Real} <: ContinuousUnivariateDistribution
μ::T
σ::T
γ1::T
SkewNormalCP(μ, σ, γ1) = new(μ, σ, γ1)#(@check_args(SkewNormalCP, σ > zero(σ)); new(μ, σ, γ1))
end
#### Outer constructors
SkewNormalCP{T<:Real}(ξ::T, ω::T, α::T) = SkewNormalCP{T}(ξ, ω, α)
SkewNormalCP{T<:Real}(α::T) = SkewNormalCP{T}(zero(α), one(α), α)
#### Parameters
params(d::SkewNormalCP) = (d.μ, d.σ, d.γ1)
mean(d::SkewNormalCP) = d.μ
var(d::SkewNormalCP) = abs2(d.σ)
std(d::SkewNormalCP) = d.σ
skewness(d::SkewNormalCP) = d.γ1
function convert{T<:Real}(::Type{SkewNormalCP{T}}, d::SkewNormalDP{T})
μ = mean(d)
σ = std(d)
γ1 = skewness(d)
return SkewNormalCP{T}(μ, σ, γ1)
end
function convert{T<:Real}(::Type{SkewNormalDP{T}}, d::SkewNormalCP{T})
# ξ = d.μ - - d.σ/std_z(d)*mean_z(d)
# ω = d.σ / std_z(d)
c = cbrt(2.0*d.γ1 / (4.0-π))
μ_z = c / √(1+c^2)
α = √(π/2.) * c / √(1+(1.0-π/2.)*c^2.)
δ = α / √(1. +α^2.)
σ_z = 1.0-2.0/π * δ^2.
ω = sqrt(d.σ^2.0 / (1.0-μ_z^2.0))
ξ = d.μ - ω*μ_z
return SkewNormalDP{T}(ξ, ω, α)
end
pdf{T<:Real}(d::SkewNormalCP{T}, x::Real) = pdf(SkewNormalDP{T}(d),x)
logpdf{T<:Real}(d::SkewNormalCP{T}, x::Real) = logpdf(SkewNormalDP{T}(d),x)
function rand{T<:Real}(d::SkewNormalCP{T})
return rand(SkewNormalDP{T}(d))
end
2 Likes
@bsprung I submitted a PR to include skew-normal in Distributions.jl: https://github.com/JuliaStats/Distributions.jl/pull/1104 (still haven’t heard back)
Let me know if you guys have any feedback
1 Like
On a related note it would be nice for the Julia community to have the remaining missing distributions in Distributions.jl.
See how fragmented the distributions ecosystem is in R.
This could be an awesome task for a student.
Probably not fancy enough for GSoC or JSoC but definitely useful.
@logankilpatrick if you’re taking suggestions for possible projects.
The projects are already set for this year but I suggest opening an issue with the “help wanted” label.