export Semicircle

"""
Wigner semicircle distribution

By default, creates a standardized distribution
with mean 0 and variance 1 (radius 2)
"""
immutable Semicircle{T<:Real} <: ContinuousUnivariateDistribution
    mean::T
    radius::T
    Semicircle(μ::T,r::T) = new(μ,r)
end
Semicircle{T<:Real}(μ::T=0.0, r::T=2.0) = r > 0 ? Semicircle{T}(μ, r) :
    throw(ArgumentError("radius r must be positive, got $r"))

# Distribution function methods
###############################

# cumulative distribution function
function cdf{T<:Real}(d::Semicircle{T}, x::T)
    r, a = d.mean, d.radius
    if insupport(d, x)
        return 0.5 + (x-a)/(π*r^2) * √(r^2 - (x-a)^2) + asin((x-a)/r)/π
    elseif x ≥ a
        return one(T)
    else
	return zero(T)
    end
end

# probability density function
function pdf{T<:Real}(d::Semicircle{T}, x::T)
    r, a = d.mean, d.radius
    if insupport(d, x)
        return 2/(π*r^2) * √(r^2 - (x-a)^2)
    else
	return zero(T)
    end
end

# predicate is x in the support of the distribution?
insupport{T<:Real}(d::Semicircle{T}, x::T)=abs(x)<d.radius

#Entropy methods
################

# entropy of distribution in nats
entropy(X::Semicircle)=log(π*X.radius) - 0.5

#Measures of central tendency methods
#####################################

# mean of distribution
mean(X::Semicircle)=X.mean

# median of distribution
median(X::Semicircle)=X.mean

# mode(s) of distribution as vector
modes{T}(X::Semicircle{T})=T[X.mean]

# kurtosis of the distribution
kurtosis{T}(X::Semicircle{T})=T(2)

# skewness of the distribution
skewness{T}(X::Semicircle{T})=T(0)

# standard deviation of distribution
std(X::Semicircle)=X.radius/2

# variance of distribution
var(X::Semicircle)=std(X)^2

# moment of distribution
function moment{T<:Real}(X::Semicircle{T}, order::Integer)
    a, r = X.mean, X.radius
    if X.mean != 0
        return a^n*hypergeom([(1-n)/2, -n/2], 2, (r/a)^2)
    elseif iseven(order)
        return (0.5*r)^(2order) * T(catalannum(order÷2))
    else
	return zero(T)
    end
end

# cumulant of distribution
function cumulant{T<:Real}(d::Semicircle{T}, order::Integer)
    if d.mean != 0
        throw(ArgumentError("Non-central Semicircle not supported"))
    end

    if order==0
	return one(T)
    elseif iseven(order)
        return (0.5*d.radius)^(2order) * T(lassallenum(order÷2))
    else
        return zero(T)
    end
end

# free cumulant of distribution
function freecumulant{T<:Real}(d::Semicircle{T}, order::Int)
    if order == 0
        return one(T)
    elseif order == 1
        return mean(X)
    elseif order == 2
        return var(X)
    else
        return zero(T)
    end
end

#Generating function methods
############################

# characteristic function
function cf(d::Semicircle, t::Real)
  r = t / d.mean
  2 * besselj(1, r)/r
end

# moment generating function
function mgf(d::Semicircle, t::Real)
  r = t * d.mean
  2 * besseli(1, r)/r
end

#Sampling methods
#################

# random sampler
# Use relationship with beta distribution
function rand(X::Semicircle)
    Y = rand(Beta(1.5, 1.5))
    X.mean + 2 * X.radius * Y - X.radius
end
