More code cleanup

- Fixed normalization bug in rand(GaussianHermite, .)
- renamed UniformHaar to Haar and added optional GSL wrappers to Haar code.
- moved LargestEigenvalue.jl code to densities/TracyWidom.jl
- added new abstract StochasticProcess type

Author: jiahao

src/GaussianEnsembles.jl

@@ -13,7 +13,8 @@
 # Alan Edelman, Per-Olof Persson and Brian D Sutton, "The fourfold way"
 # http://www-math.mit.edu/~edelman/homepage/papers/ffw.pdf
 
-export Wigner, GaussianHermite, GaussianLaguerre, GaussianJacobi 
+export GaussianHermite, GaussianLaguerre, GaussianJacobi,
+       Wigner, Wishart, MANOVA #Synonyms for Hermite, Laguerre and Jacobi
 
 #A convenience function to define a chi scalar random variable
 chi(df::Real) = sqrt(rand(Chisq(df)))
@@ -25,7 +26,7 @@ chi(df::Real) = sqrt(rand(Chisq(df)))
 type GaussianHermite <: ContinuousMatrixDistribution
   beta::Real
 end
-typealias Wigner GaussianHermite
+typealias Wigner  GaussianHermite
 
 # Generates a NxN symmetric Wigner matrix
 function rand(d::GaussianHermite, n::Integer)
@@ -41,7 +42,7 @@ function rand(d::GaussianHermite, n::Integer)
   else
     error(@sprintf("beta = %d is not implemented", d.beta))
   end
-  normalization = sqrt(2*d.beta)
+  normalization = sqrt(2*d.beta*n)
   return (A + A') / normalization
 end
 
@@ -63,7 +64,6 @@ tridrand(d::GaussianHermite, dims::Dim2) = dims[1]==dims[2] ? tridrand(d, dims[1
 #Return n eigenvalues distributed according to the Hermite ensemble
 eigvalrand(d::GaussianHermite, n::Integer) = eigvals(tridrand(d, b))
 
-
 #Calculate Vandermonde determinant term
 function VandermondeDeterminant{Eigenvalue<:Number}(lambda::Vector{Eigenvalue}, beta::Real)
   n = length(lambda)
@@ -83,10 +83,7 @@ function eigvaljpdf{Eigenvalue<:Number}(d::GaussianHermite, lambda::Vector{Eigen
   for j=1:n
     c *= gamma(1 + beta/2)/gamma(1 + beta*j/2)
   end
-
-  #Calculate argument of exponential
-  Energy = sum(lambda.^2/2)
-
+  Energy = sum(lambda.^2/2) #Calculate argument of exponential
   VandermondeDeterminant(lambda, beta) * exp(-Energy)
 end
 
@@ -100,6 +97,8 @@ type GaussianLaguerre <: ContinuousMatrixDistribution
   beta::Real
   a::Real
 end
+typealias Wishart GaussianLaguerre
+
 
 #Generates a NxN Hermitian Wishart matrix
 #TODO Check - the eigenvalue distribution looks funky
@@ -172,6 +171,7 @@ type GaussianJacobi <: ContinuousMatrixDistribution
   a::Real
   b::Real
 end
+typealias MANOVA  GaussianJacobi
 
 function rand(d::GaussianJacobi, m::Integer)
   w1 = Wishart(m, int(2.0*d.a/d.beta), d.beta)

src/Haar.jl

@@ -1,9 +1,10 @@
-using GSL
 using Catalan
 
-export permutations_in_Sn, compose, cycle_structure, data, part #Functions working with partitions and permutations
-    UniformHaar, expectation, WeingartenUnitary
+export permutations_in_Sn, compose, cycle_structure, data, part, #Functions working with partitions and permutations
+    partition, Haar, expectation, WeingartenUnitary
 
+
+typealias partition Vector{Int}
 #Functions working with partitions and permutations
 # TODO Migrate these functions to Catalan
 
@@ -31,6 +32,9 @@ function permutations_in_Sn(n::Integer)
     end
 end
 
+
+if _HAVE_GSL
+
 function compose(P::Ptr{gsl_permutation}, Q::Ptr{gsl_permutation})
     #Compose the permutations
     n=convert(Int64, permutation_size(P))
@@ -64,8 +68,10 @@ 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 UniformHaar <: AbstractMatrix
+type Haar <: ContinuousMatrixDistribution
     beta::Real
 end
 
@@ -105,20 +111,20 @@ function expectation(X::Expr)
     for i=1:n
         thingy=X.args[i+1]
         if isa(thingy, Symbol)
-            if isa(eval(thingy), UniformHaar)
+            if isa(eval(thingy), Haar)
                 if MyQ==None MyQ=thingy end
                 if MyQ == thingy
                     Qidx=[Qidx; i]
                 else
-                    warning("only one instance of UniformHaar supported, skipping the other guy ", thingy) end
+                    warning("only one instance of Haar supported, skipping the other guy ", thingy) end
             else
                 Others = [Others; (thingy, i, i+1)]
             end
             println(i, ' ', thingy, "::", typeof(eval(thingy)))
         elseif isa(thingy, Expr)
             println(i, ' ', thingy, "::Expr")
             if thingy.head==symbol('\'') && length(thingy.args)>=1 #Maybe this is a Q'
-                if isa(thingy.args[1], Symbol) && isa(eval(thingy.args[1]), UniformHaar)
+                if isa(thingy.args[1], Symbol) && isa(eval(thingy.args[1]), Haar)
                     println("Here is a Qtranspose")
                     Qpidx=[Qpidx; i]
                 end
@@ -264,11 +270,16 @@ 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/HaarMeasure.jl

@@ -21,7 +21,7 @@ import Base.rand
 #    Edelman and Rao, 2005
 #    Mezzadri, 2006, math-ph/0609050
 #TODO implement O(n^2) method
-function rand(W::UniformHaar, doCorrection::Integer)
+function rand(W::Haar, doCorrection::Integer)
     n, beta = W.N, W.beta
     M=rand(Ginibre(n, beta))
     q,r=qr(M)
@@ -54,7 +54,7 @@ end
 #By default, always do piecewise correction
 #For most applications where you use the HaarMatrix as a similarity transform
 #it doesn't matter, but better safe than sorry... let the user choose else
-rand(W::UniformHaar) = rand(W, 1)
+rand(W::Haar) = rand(W, 1)
 
 #A utility method to check if the piecewise correction is needed
 #This checks the R part of the QR factorization; if correctly done,

src/LargestEigenvalue.jl

@@ -26,14 +26,27 @@ typealias Dim2 (Int, Int) #Dimensions of a rectangular matrix
 #Classical Gaussian matrix ensembles
 include("GaussianEnsembles.jl")
 
+
 # Classical univariate distributions
 include("densities/Semicircle.jl")
 include("densities/TracyWidom.jl")
 
-#Fast histogrammer for matrix eigenvalues
+#Generating matrices of Haar measure
+include("Haar.jl")
+include("HaarMeasure.jl")
+include("HaarSymbolic.jl")
+
+#Fast histogrammer for matrix eigenvalues - hist_eig
 include("FastHistogram.jl")
 
 #Formal power series
 include("FormalPowerSeries.jl")
+
+#Statistical tests based on random matrix theory
+include("StatisticalTests.jl")
+
+#Stochastic processes
+include("StochasticProcess.jl")
+
 end
 

src/RandomMatrices.jl

@@ -1,15 +1,31 @@
+#Matrices related to stochastic processes
+
+abstract StochasticProcess
+export evolve
+
+immutable BrownianProcess <: StochasticProcess
+end
+
 # Generates Brownian motion
-BrownianProcess(dx::Real, xend::Real) = cumsum(WhiteNoise(dx, xend))
+evolve(p::BrownianProcess, dx::Real, xend::Real) = cumsum(WhiteNoiseProcess(dx, xend))
 
 # Generates a white noise process
 function WhiteNoiseProcess(dx::Real, xend::Real)
-    x=[0:dx:xend]
-    dW=randn(length(x),1)*sqrt(dx)
+  x=[0:dx:xend]
+  randn(length(x),1)*sqrt(dx)
+end
+
+################
+# Airy process #
+################
+
+immutable AiryProcess <: StochasticProcess
+  beta::Real
 end
 
 # Calculates the largest eigenvalue of a stochastic Airy process
 # with Brownian noise
-function StochasticAiryProcess(dx::Real, xend::Real, beta::Real)
+function evolve(p::AiryProcess, dx::Real, xend::Real)
     x=[0:dx:xend]
     N=length(x)
 
@@ -18,14 +34,14 @@ function StochasticAiryProcess(dx::Real, xend::Real, beta::Real)
     b=+(1/dx^2)*ones(N-1)
     #Plus noise
     dW=WhiteNoiseProcess(dx, xend)
-    a+=(2/sqrt(beta))*dW/dx
+    a+=(2/sqrt(p.beta))*dW/dx
 
     maxeig(SymTridiagonal(a,b))
 end
 
 #Sample program
 #t=10000 #number of trials
-#v=[StochasticAiryProcess(0.001, 10, 2) for i=1:t]
+#v=[AiryProcess(0.001, 10, 2) for i=1:t]
 #binsize=.2
 #grid=[-5:binsize:2]
 #x=hist(v, grid)

src/StochasticProcess.jl

@@ -4,7 +4,7 @@ importall Distributions
 using Catalan
 export Semicircle
 
-type Semicircle <: ContinuousUnivariateDistribution
+immutable Semicircle <: ContinuousUnivariateDistribution
     mean::Float64
     radius::Float64
     Semicircle(mu, r) = r > 0 ? new(float64(mu), float64(r)) : 

src/densities/Semicircle.jl

@@ -1,11 +1,15 @@
 #Computes the Tracy-Widom distribution by directly solving the
 #Painlevé II equation
+importall Distributions
 export TracyWidom
 using ODE
 
 include("gradient.jl")
 
-function TracyWidom(t::Real)
+immutable TracyWidom <: ContinuousUnivariateDistribution
+end
+
+function pdf(d::TracyWidom, t::Real)
   t0 = -8.0
   t>t0 ? (return 0.0) : nothing
   t>5  ? (return 0.0) : nothing
@@ -18,3 +22,12 @@ function TracyWidom(t::Real)
   f2=(F2[end]-F2[end-1])/(ts[end]-ts[end-1])# the density at t
 end
 
+
+#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))
+end
+