1
$\begingroup$

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?

$\endgroup$
5
$\begingroup$

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
$\endgroup$
3
$\begingroup$

Mathematica has the TW distributions:

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.