Add test of Tracy-Widom distribution

Make GSL a required dependency, and get rid of conditional testing code

Drop obsolete dependency on Docile

Author: jiahao

REQUIRE

@@ -1,5 +1,6 @@
 julia 0.4
 Combinatorics 0.2
-Docile
 Distributions
 ODE
+GSL
+

src/Haar.jl

@@ -31,8 +31,6 @@ function permutations_in_Sn(n::Integer)
 end
 
 
-if _HAVE_GSL
-
 function compose(P::Ptr{gsl_permutation}, Q::Ptr{gsl_permutation})
     #Compose the permutations
     n=convert(Int64, permutation_size(P))
@@ -66,8 +64,6 @@ end
 data(P::Ptr{gsl_permutation}) = [convert(Int64, x)+1 for x in
     pointer_to_array(permutation_data(P), (convert(Int64, permutation_size(P)) ,))]
 
-end #_HAVE_GSL
-
 
 type Haar <: ContinuousMatrixDistribution
     beta::Real
@@ -268,16 +264,12 @@ function expectation(X::Expr)
     eval(X)
 end
 
-if _HAVE_GSL
-
 #Computes the Weingarten function for permutations
 function WeingartenUnitary(P::Ptr{gsl_permutation})
     C = cycle_structure(P)
     WeingartenUnitary(C)
 end
 
-end #_HAVE_GSL
-
 #Computes the Weingarten function for partitions
 function WeingartenUnitary(P::partition)
     n = sum(P)

src/RandomMatrices.jl

@@ -1,25 +1,19 @@
 module RandomMatrices
-import Distributions
-using Combinatorics,Compat
 
-import Distributions: ContinuousMatrixDistribution, Chi
-
-
-#If the GNU Scientific Library is present, turn on additional functionality.
-_HAVE_GSL = try
-  using GSL
-  true
-catch
-   false
-end
+using Combinatorics
+using Compat
+using GSL
+using ODE
 
 import Base: isinf, rand
+import Distributions: ContinuousMatrixDistribution, Chi, pdf
 
 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
+       eigvaljpdf, #Eigenvalue joint probability density
+       pdf
 
 typealias Dim2 Tuple{Int, Int} #Dimensions of a rectangular matrix
 

src/densities/Semicircle.jl

@@ -66,19 +66,13 @@ std(X::Semicircle)=X.radius/2
 var(X::Semicircle)=std(X)^2
 
 # moment of distribution
-if _HAVE_GSL
-  function moment(X::Semicircle, order::Integer)
+function moment(X::Semicircle, order::Integer)
     a, r = X.mean, X.radius
     if X.mean != 0
-      a^n*hypergeom([(1-n)/2, -n/2], 2, (r/a)^2)
+        a^n*hypergeom([(1-n)/2, -n/2], 2, (r/a)^2)
     else
-      order%2 ? (0.5*r)^(2n) * catalan(div(order,2)) : 0
+        order%2 ? (0.5*r)^(2n) * catalan(div(order,2)) : 0
     end
-  end
-else
-  function moment(X::Semicircle, order::Integer)
-    error("The moment of a semicircle distribution requires the GSL module to be loaded.")
-  end
 end
 
 # cumulant of distribution

src/densities/TracyWidom.jl

@@ -1,32 +1,40 @@
-#Computes the Tracy-Widom distribution by directly solving the
-#Painlevé II equation
-import Distributions
 export TracyWidom
-using ODE
 
 include("gradient.jl")
 
-immutable TracyWidom <: ContinuousUnivariateDistribution
-end
+immutable TracyWidom <: ContinuousUnivariateDistribution end
+
 
+"""
+Computes the Tracy-Widom distribution by directly solving the
+Painlevé II equation by numerical integration
+"""
 function pdf(d::TracyWidom, t::Real)
-  t0 = -8.0
-  t>t0 ? (return 0.0) : nothing
-  t>5  ? (return 0.0) : nothing
+    t0 = -8.0
+    t≤t0 && return 0.0
+    t≥5  && return 0.0
 
-  deq(t, y) = [y[2]; t*y[1]+2y[1]^3; y[4]; y[1]^2]
+    deq(t, y) = [y[2]; t*y[1]+2y[1]^3; y[4]; y[1]^2]
 
-  y0=[airy(t0); airy(1, t0); 0; airy(t0)^2] # Initial conditions
-  (ts, y)=ode23(deq, [t0, t], y0)           # Solve the ODE
-  F2=exp(-y[:,3])                           # the cumulative distribution
-  f2=(F2[end]-F2[end-1])/(ts[end]-ts[end-1])# the density at t
+    y0=[airy(t0); airy(1, t0); 0; airy(t0)^2] # Initial conditions
+    (ts, y)=ode23(deq, [t0, t], y0)           # Solve the ODE
+    F2=exp(-y[:,3])                           # the cumulative distribution
+    f2=(F2[end]-F2[end-1])/(ts[end]-ts[end-1])# the density at t
 end
+pdf(d::Type{TracyWidom}, t::Real) = pdf(d(), t)
 
-
-#Samples the largest eigenvalue of the n x n GUE matrix
+"""
+Samples the largest eigenvalue of the n x n GUE matrix
+"""
 function rand(d::TracyWidom, n::Integer)
-  k=int(n-10*n^(1/3)-1)
-  A=[chi(i) for i=(n-1):-1:k]
-  B=randn(n-k+1)
-  v=eigmax(SymTridiagonal(B, A))
+    n > 1 || error("Cannot construct $n × $n matrix")
+    if n < 100
+        k = n
+    else #Exploit the fact that the eigenvector is concentrated in the top left corner
+	k = round(Int, n-10*n^(1/3)-1)
+    end	
+    a=randn(n-k+1)
+    b=[χ(i) for i=(n-1):-1:k]
+    v=eigmax(SymTridiagonal(a, b))
 end
+rand(d::Type{TracyWidom}, t::Integer) = rand(d(), t)

test/densities/TracyWidom.jl

@@ -0,0 +1,14 @@
+using RandomMatrices
+using Base.Test
+
+#Test far outside support
+#Tracy-Widom has support on all x>0, but the integration won't pick it up
+@test RandomMatrices.pdf(TracyWidom, -10) == RandomMatrices.pdf(TracyWidom, 10) == 0
+
+if isdefined(:ODE) && isa(ODE, Module)
+    t = rand()
+    @test RandomMatrices.pdf(TracyWidom, t) > 0
+end
+
+@test rand(TracyWidom, 10) > 0
+@test rand(TracyWidom, 100) > 0

test/runtests.jl

@@ -6,3 +6,6 @@ include("GaussianEnsembles.jl")
 include("FormalPowerSeries.jl")
 include("Haar.jl")
 include("StochasticProcess.jl")
+
+include("densities/TracyWidom.jl")
+