# Data for Tracy-Widom distribution

I have data for a random variable and I wish to test whether it conforms to the Tracy-Widom distribution. However, the T-W distribution is hard to compute. Is there a readly available table of its values I can just download?

It's actually not all that complicated to calculate the Tracy-Widom CDF just from its definition: see On The Numerical Evaluation Of Fredholm Determinants by Folkmar Bornemann. The Wikipedia page gives the definition as $$F_2(s) = \det(I - A_s),$$ where $A_s$ is the integral operator on $[s,\infty)$ with the kernel $$\DeclareMathOperator{\Ai}{Ai} K(x,y) = \frac{\Ai(x)\Ai'(y)-\Ai'(x)\Ai(y)}{x-y},$$ $$A_s[f](x) = \int_s^\infty K(x,y)f(y)\,\mathrm{d}y.$$

Just going from this, you can take the weights $w_i$ and nodes $x_i$ of any good quadrature rule on the interval $[s,\infty)$, such as the excellent tanh-sinh quadrature, and compute the determinant of the matrix $$\det\big( \delta_{ij} - w_i^{1/2} K(x_i, x_j) w_j^{1/2}\big)_{i,j}.$$ This works because the kernel is smooth.

To make sure this actually works, I wrote the code below (using Julia 0.6): there is a little bit of inaccuracy at $F_2(-8)$, but it's fine otherwise. Most of the time is spent calculating Airy functions, so it should be okay for generating tables of values for future use.

module TracyWidom

using SpecialFunctions
using Base.Test

function K(x, y)
try
if x == y
airyaiprime(x)^2 - x * airyai(x)^2
else
(airyai(x)*airyaiprime(y) - airyaiprime(x)*airyai(y)) / (x - y)
end
catch ex # underflow is not an error here
if isa(ex, SpecialFunctions.AmosException) && ex.id == 4
return 0.0
else
rethrow(ex)
end
end
end

function tanhSinhNodes(s, n=64)
tf = 3.158587875883174 # fzero(t -> log1p(-tanh(pi/2*sinh(t))) - log(eps(typeof(s))), 3.1)
t = linspace(-tf, tf, 2n+1)
u = tanh.(pi/2*sinh.(t))
w = pi/2*cosh.(t)./cosh.(pi/2*sinh.(t)).^2
w * step(t) .* 2./(u-1).^2, s + (u+1) ./ (1-u)
end

function tanhSinh(f, s, n=64)
w, x = tanhSinhNodes(s, n)
dot(w, f.(x))
end

function F2(s, n=128)
w, x = tanhSinhNodes(s, n)
m = length(w)
det([Int(i==j) - sqrt(w[i])*sqrt(w[j])*K(x[i], x[j]) for i=1:m, j=1:m])
end

function test()
@testset "TracyWidom" begin
@testset "tanh-sinh" begin
@test tanhSinh(x -> 1/(1+x^2), 0., 32) ≈ π/2
@test tanhSinh(x -> 1/(1+x^2), 1., 32) ≈ π/2 - atan(1.)
@test tanhSinh(x -> (1+x^2)^-1.25, 2.0, 32) ≈ 0.20885622523501382
end

@testset "F2" begin
@test F2(0.1) ≈ 0.9754704606594619
@test F2(-0.1) ≈ 0.9620142937272265
@test F2(-1.0) ≈ 0.8072142419992853
@test F2(-2.0) ≈ 0.41322414250512257
@test F2(-3.0) ≈ 0.08031955293933454
@test F2(-4.0) ≈ 0.0035445535955092003
@test F2(-5.0) ≈ 2.135996984741116e-5
@test F2(-6.0) ≈ 1.062254674124451e-8
@test F2(-7.0) ≈ 2.639614767246062e-13
@test F2(-8.0) ≈ 1.9859004257636574e-19

@test F2(1.0) ≈ 0.9975054381493893
@test F2(2.0) ≈ 0.9998875536983092
@test F2(3.0) ≈ 0.9999970059566077
@test F2(4.0) ≈ 0.9999999504208784
@test F2(5.0) ≈ 0.9999999994682207
@test F2(6.0) ≈ 0.9999999999961827
@test F2(7.0) ≈ 0.9999999999999811
@test F2(8.0) ≈ 0.9999999999999999
end
end
nothing
end

end


Mathematica has the TW distributions:

http://reference.wolfram.com/language/ref/TracyWidomDistribution.html