Add cdf(::TracyWidom)

Author: jiahao

src/RandomMatrices.jl

@@ -6,14 +6,15 @@ using GSL
 using ODE
 
 import Base: isinf, rand
-import Distributions: ContinuousMatrixDistribution, Chi, pdf
+import Distributions: ContinuousMatrixDistribution, Chi, cdf, 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
-       pdf
+       pdf,
+       cdf
 
 typealias Dim2 Tuple{Int, Int} #Dimensions of a rectangular matrix
 

src/densities/TracyWidom.jl

@@ -9,21 +9,38 @@ 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
+function pdf{S<:Real}(d::TracyWidom, t::S, t0::S = convert(S, -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]
+    ts, y = _solve_painleve_ii(t0, t)
+
+    ΔF2=exp(-y[end-1][3]) - exp(-y[end][3]) # the cumulative distribution
+    f2=ΔF2/(ts[end]-ts[end-1])              # the density at t
+end
+pdf(d::Type{TracyWidom}, t::Real) = pdf(d(), t)
 
+"""
+Computes the Tracy-Widom distribution by directly solving the
+Painlevé II equation by numerical integration
+"""
+function cdf{S<:Real}(d::TracyWidom, t::S, t0::S = convert(S, -8.0))
+    t≤t0 && return 0.0
+    t≥5  && return 1.0
+
+    ts, y = _solve_painleve_ii(t0, t)
+
+    F2=exp(-y[end][3])
+end
+cdf(d::Type{TracyWidom}, t::Real) = cdf(d(), t)
+
+function _solve_painleve_ii{S<:Real}(t0::S, t::S)
+    deq(t, y) = [y[2], t*y[1]+2y[1]^3, y[4], y[1]^2]
     a0 = airy(t0)
     T = typeof(a0)
     y0=T[a0, airy(1, t0), 0, airy(t0)^2]    # Initial conditions
     (ts, y)=ode23(deq, y0, [t0, t])         # Solve the ODE
-    ΔF2=exp(-y[end-1][3]) - exp(-y[end][3]) # the cumulative distribution
-    f2=ΔF2/(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

test/densities/TracyWidom.jl

@@ -3,11 +3,14 @@ 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
+@test pdf(TracyWidom, -10) == pdf(TracyWidom, 10) == 0
+@test cdf(TracyWidom, -10) == 0
+@test cdf(TracyWidom, 10) == 1
 
 if isdefined(:ODE) && isa(ODE, Module)
     t = rand()
-    @test RandomMatrices.pdf(TracyWidom, t) > 0
+    @test pdf(TracyWidom, t) > 0
+    @test 0 < cdf(TracyWidom, t) < 1
 end
 
 @test isfinite(rand(TracyWidom, 10))