Skew normal distribution?


#1

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!


#2

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.
https://www.jstatsoft.org/pages/view/mission

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?


#3

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?


#4

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

#5

You could also try the sample() function from ApproxFun, click here.


#6

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).


#7

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