Test Semicircle law

Add missing imports and exports

Make a lot of the Semicircle methods more type stable

Also modernize code - 4 space indents, use of iseven(),

Author: jiahao

src/RandomMatrices.jl

@@ -6,15 +6,18 @@ using GSL
 using ODE
 
 import Base: isinf, rand
-import Distributions: ContinuousMatrixDistribution, Chi, cdf, pdf
+import Distributions: ContinuousUnivariateDistribution,
+    ContinuousMatrixDistribution,
+    Beta, Chi,
+    cdf, pdf, entropy, insupport, mean, median, modes, kurtosis, skewness, std, var, moment
 
 export bidrand,    #Generate random bidiagonal matrix
        tridrand,   #Generate random tridiagonal matrix
        sprand,     #Generate random sparse matrix
        eigvalrand, #Generate random set of eigenvalues
        eigvaljpdf, #Eigenvalue joint probability density
-       pdf,
-       cdf
+       cumulant, freecumulant, cf, mgf,
+       cdf, pdf, entropy, insupport, mean, median, modes, kurtosis, skewness, std, var, moment
 
 typealias Dim2 Tuple{Int, Int} #Dimensions of a rectangular matrix
 

src/densities/Semicircle.jl

@@ -1,45 +1,51 @@
-#Wigner semicircle distribution
-
-import Distributions
-import Distributions: ContinuousUnivariateDistribution
-
 export Semicircle
 
-immutable Semicircle <: ContinuousUnivariateDistribution
-    mean::Float64
-    radius::Float64
-    Semicircle(mu, r) = r > 0 ? new(float64(mu), float64(r)) :
-        error("radius r must be positive") #Constructor
-end
+"""
+Wigner semicircle distribution
 
-#Standardized distribution - mean 0, variance 1
-Semicircle(mu) = Semicircle(mu, 2.0)
-Semicircle() = Semicircle(0.0, 2.0)
+By default, creates a standardized distribution
+with mean 0 and variance 1 (radius 2)
+"""
+immutable Semicircle{T<:Real} <: ContinuousUnivariateDistribution
+    mean::T
+    radius::T
+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(X::Semicircle, x::Real)
-    r, a = X.mean, X.radius
-    insupport(x) ? 0.5 + (x-a)/(pi*r^2) * sqrt(r^2 - (x-a)^2) + 1/pi * asin((x-a)/r) : (x>a ? 1.0 : 0.0)
+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(X::Semicircle, x::Real)
-    r, a = X.mean, X.radius
-    insupport(x) ? 2/(pi*r^2) * sqrt(r^2 - (x-a)^2) : 0.0
+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(X::Semicircle, x::Real)=abs(x)<=X.radius
-
+insupport{T<:Real}(d::Semicircle{T}, x::T)=abs(x)<d.radius
 
 #Entropy methods
 ################
 
 # entropy of distribution in nats
-entropy(X::Semicircle)=log(pi*X.radius) - 0.5
+entropy(X::Semicircle)=log(π*X.radius) - 0.5
 
 #Measures of central tendency methods
 #####################################
@@ -51,13 +57,13 @@ mean(X::Semicircle)=X.mean
 median(X::Semicircle)=X.mean
 
 # mode(s) of distribution as vector
-modes(X::Semicircle)=[X.mean]
+modes{T}(X::Semicircle{T})=T[X.mean]
 
 # kurtosis of the distribution
-kurtosis(X::Semicircle)=2
+kurtosis{T}(X::Semicircle{T})=T(2)
 
 # skewness of the distribution
-skewness(X::Semicircle)=0
+skewness{T}(X::Semicircle{T})=T(0)
 
 # standard deviation of distribution
 std(X::Semicircle)=X.radius/2
@@ -66,51 +72,57 @@ std(X::Semicircle)=X.radius/2
 var(X::Semicircle)=std(X)^2
 
 # moment of distribution
-function moment(X::Semicircle, order::Integer)
+function moment{T<:Real}(X::Semicircle{T}, order::Integer)
     a, r = X.mean, X.radius
     if X.mean != 0
-        a^n*hypergeom([(1-n)/2, -n/2], 2, (r/a)^2)
+        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
-        order%2 ? (0.5*r)^(2n) * catalan(div(order,2)) : 0
+	return zero(T)
     end
 end
 
 # cumulant of distribution
-function cumulant(X::Semicircle, order::Integer)
-  if X.mean != 0 error("not supported") end
-  if order%2
-    order==0 ? 1 : (0.5*r)^(2n) * lassalle(order/2)
-  else
-    0
-  end
+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(X::Semicircle, order::Integer)
-  if order == 0
-    return 1
-  elseif order == 1
-    return mean(X)
-  elseif order == 2
-    return var(X)
-  else
-    return 0
-  end
+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(X::Semicircle, t::Real)
-  r = t * X.mean
+function cf(d::Semicircle, t::Real)
+  r = t / d.mean
   2 * besselj(1, r)/r
 end
 
 # moment generating function
-function mgf(X::Semicircle, t::Real)
-  r = t * X.mean
+function mgf(d::Semicircle, t::Real)
+  r = t * d.mean
   2 * besseli(1, r)/r
 end
 
@@ -120,6 +132,6 @@ end
 # 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
+    Y = rand(Beta(1.5, 1.5))
+    X.mean + 2 * X.radius * Y - X.radius
 end

test/densities/Semicircle.jl

@@ -0,0 +1,30 @@
+using RandomMatrices
+using Base.Test
+
+let
+    d = Semicircle()
+    @test cdf(d, 2.0) == 1.0
+    @test cdf(d, -2.0) == 0.0
+    @test pdf(d, 2.0) == 0.0
+    @test pdf(d, -2.0) == 0.0
+
+    @test entropy(d) == log(2π) - 0.5
+    @test mean(d) == median(d) == skewness(d) == 0
+    @test modes(d) == [0]
+    @test std(d) == var(d) == 1
+    @test kurtosis(d) == 2
+
+    @test isfinite(rand(d))
+    @test moment(d, 4) == 2
+    @test cumulant(d, 4) == 1
+    @test moment(d, 5) == cumulant(d, 5) == freecumulant(d, 5) == 0
+
+    @test cf(d, 1.0) == 0.0
+
+    #XXX Throws AmosException - ???
+    #@show mgf(d, 1.0)
+
+    x = rand(d)
+    @test insupport(d, x)
+end
+

test/runtests.jl

@@ -7,5 +7,6 @@ include("FormalPowerSeries.jl")
 include("Haar.jl")
 include("StochasticProcess.jl")
 
+include("densities/Semicircle.jl")
 include("densities/TracyWidom.jl")