rmt book chapter 2 (2.3 to 2.7)

Author: vtjnash

demos/book/1/bellcurve.jl

@@ -6,20 +6,19 @@
 #Theory:      The normal distribution curve
 
 ## Experiment
-t=1000000
-dx=.2
-v=randn(t)
-grid=-4:dx:4
-count=hist(v,grid)[2]/(t*dx)
+t = 1000000
+dx = .2
+v = randn(t)
+x = -4:dx:4
+grid,count = hist(v,x)
 
 ## Theory
-x=grid
-y=exp(-x.^2/2)/sqrt(2*pi)
+y = exp(-x.^2/2)/sqrt(2*pi)
 
 ## Plot
 using Winston
 p = FramedPlot()
-h = Histogram(count, step(grid))
+h = Histogram(count/(t*dx), step(grid))
 h.x0 = first(grid)
 add(p, h)
 add(p, Curve(x, y, "color", "blue", "linewidth", 2)) 

demos/book/2/haar.jl

@@ -0,0 +1,47 @@
+#haar.jl
+#Algorithm 2.4 of Random Eigenvalues by Alan Edelman
+
+#Experiment:    Generate random orthogonal/unitary matrices
+#Plot:          Histogram eigenvalues
+#Theory:        Eigenvalues are on unit circle
+
+## Parameters
+t = 5000        # trials
+dx = 0.05       # binsize
+n = 10          # matrix size
+
+function haar_experiment(t,dx,n)
+    v = Complex{Float64}[] # eigenvalue samples
+    for i = 1:t
+        # Sample random orthogonal matrix
+        # X = QR(randn(n,n))[:Q]
+        #  If you have non-uniformly sampled eigenvalues, you may need this fix:
+        #X = X*diagm(sign(randn(n,1)))
+        #
+        # Sample random unitary matrix
+        X = QR(randn(n,n)+im*randn(n,n))[:Q]
+        #  If you have non-uniformly sampled eigenvalues, you may need this fix:
+        #X = X*diagm(sign(randn(n)+im*randn(n)))
+        append!(v, eigvals(full(X)))
+    end
+    x = ((-1+dx/2):dx:(1+dx/2))*pi
+    x, v
+end
+x, v = haar_experiment(t,dx,n)
+grid, count = hist(angle(v), x)
+
+## Theory
+h2 = t*n*dx/2*x.^0
+
+## Plot
+using Winston
+p = FramedPlot()
+h = Histogram(count, step(grid))
+h.x0 = first(grid)
+add(p, h)
+add(p, Curve(x, h2, "color", "blue", "linewidth", 2)) 
+if isinteractive()
+    Winston.display(p)
+else
+    file(p, "haar.png")
+end

demos/book/2/patiencesort.jl

@@ -0,0 +1,21 @@
+#patiencesort.jl
+#Algorithm 2.5 of Random Eigenvalues by Alan Edelman
+
+function patiencesort(p)
+    piles = similar(p, 0)
+    for pi in p
+        idx = 1
+        for pp in piles
+            if pi > pp
+                idx += 1
+            end
+        end
+        if idx > length(piles)
+            d = length(piles)
+            resize!(piles, idx)
+            piles[d+1:end] = 0
+        end
+        piles[idx] = pi
+    end
+    return length(piles)
+end

demos/book/2/spherecomponent.jl

@@ -0,0 +1,33 @@
+#spherecomponent.jl
+#Algorithm 2.3 of Random Eigenvalues by Alan Edelman
+
+#Experiment:    Random components of a vector on a sphere
+#Plot:          Histogram of a single component
+#Theory:        Beta distribution
+
+## Parameters
+t = 100000      # trials
+n = 12          # dimensions of sphere
+dx = 0.1        # bin size
+
+## Experiment
+v = randn(n, t)
+v = (v[1,:] ./ sqrt(sum(v.^2,1)))'
+grid, count = hist(v, -1:dx:1)
+
+## Theory
+x = grid
+y = (gamma(n-1)/(2^(n-2)*gamma((n-1)/2)^2))*(1-x.^2).^((n-3)/2)
+
+## Plot
+using Winston
+p = FramedPlot()
+h = Histogram(count/(t*dx), step(grid))
+h.x0 = first(grid)
+add(p, h)
+add(p, Curve(x, y, "color", "blue", "linewidth", 2))
+if isinteractive()
+    Winston.display(p)
+else
+    file(p, "spherecomponent.png")
+end

demos/book/2/tracywidomlis.jl

@@ -0,0 +1,34 @@
+#tracywidomlis.jl
+#Algorithm 2.7 of Random Eigenvalues by Alan Edelman
+
+#Experiment:    Sample random permutations
+#Plot:          Histogram of lengths of longest increasing subsequences
+#Theory:        The Tracy-Widom law
+
+## Parameters
+t = 10000       # number of trials
+n = 6^6         # length of permutations
+dx = 1/6        # bin size
+
+## Experiment
+include("patiencesort.jl")
+function tracywidomlis(t, n, dx)
+    v = zeros(t) # samples
+    for i=1:t
+        v[i] = patiencesort(randperm(n))
+    end
+    w = (v-2sqrt(n))/n^(1/6)
+    return hist(w, -5:dx:2)
+end
+grid, count = tracywidomlis(t, n, dx)
+
+## Plot
+h = Histogram(count/(t*dx), step(grid))
+h.x0 = first(grid)
+include("../1/tracywidom.jl")
+add(p, h)
+if isinteractive()
+    Winston.display(p)
+else
+    file(p, "tracywidomlis.png")
+end

demos/book/2/unitarylis.jl

@@ -0,0 +1,30 @@
+#unitarylis.jl
+#Algorithm 2.5 of Random Eigenvalues by Alan Edelman
+
+#Experiment:    Generate random orthogonal/unitary matrices
+#Theory:        Counts longest increasing subsequence statistics
+
+## Parameters
+t = 200000      # Number of trials
+n = 4           # permutation size
+k = 2           # length of longest increasing subsequence
+
+include("patiencesort.jl")
+
+function unitarylis(t, n, k)
+    v = zeros(t) # samples
+    ## Experiment
+    for i = 1:t
+        X = QR(complex(randn(k,k), randn(k,k)))[:Q]
+        X = X*diagm(sign(complex(randn(k),randn(k))))
+        v[i] = abs(trace(X)) ^ (2n)
+    end
+    z = mean(v)
+    c = 0
+    for i=1:factorial(n)
+        c = c + int(patiencesort(nthperm([1:n],i))<=k)
+    end
+    return (z, c)
+end
+
+unitarylis(t, n, k)