# Skew normal distribution?

**URL:** <https://discourse.julialang.org/t/skew-normal-distribution/21549>\
**Category:** Statistics\
**Created:** [March 6, 2019, 9:03pm UTC](https://discourse.julialang.org/t/skew-normal-distribution/21549 "2019-03-06T21:03:05Z")\
**Posts on this page:** 1\
**Showing post:** 7

<div class="post-metadata">

**Author:** ![maximerischard](https://avatars.discourse-cdn.com/v4/letter/m/c77e96/32.png) [@maximerischard](https://discourse.julialang.org/u/maximerischard)\
**Post date:** [March 16, 2019, 11:04pm UTC](https://discourse.julialang.org/t/skew-normal-distribution/21549/7 "2019-03-16T23:04:20Z")

</div>

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

```julia
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

```

---

_[View the full topic](https://discourse.julialang.org/t/skew-normal-distribution/21549)._