Merge pull request #1 from vtjnash/jn

Code from Alan's RMT book

Author: jiahao

demos/book/1/bellcurve.jl

@@ -6,21 +6,25 @@
 #Theory:      The normal distribution curve
 
 ## Experiment
-t=1000000
-dx=.2
-v=randn(t)
-grid=[-4:dx:4]
-count=hist(v,grid)/(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, dx)
-h.x0 = grid[1]
+h = Histogram(count/(t*dx), step(grid))
+h.x0 = first(grid)
 add(p, h)
 add(p, Curve(x, y, "color", "blue", "linewidth", 2)) 
-file(p, "bellcurve.png")
+if isinteractive()
+    Winston.display(p)
+else
+    file(p, "bellcurve.png")
+end
+

demos/book/1/largeeig.jl

@@ -10,28 +10,36 @@ include("tracywidom.jl")
 ## Parameters
 n=100        # matrix size
 t=5000       # trials
-v=[]         # eigenvalue samples
 dx=.2        # binsize
-grid=[-5:dx:2]
-## Experiment
-for i=1:t
-    a=randn(n,n)+im*randn(n,n) # random nxn complex matrix
-    s=(a+a')/2                 # Hermitian matrix
-    v=[v; max(eig(s)[1])]         # Largest Eigenvalue
-end
-v=n^(1/6)*(v-2*sqrt(n))        # normalized eigenvalues
+function largeeig_experiment(n,t,dx)
+    ## Experiment
+    v=Float64[] # eigenvalue samples
+    for i=1:t
+        a=randn(n,n)+im*randn(n,n) # random nxn complex matrix
+        s=(a+a')/2                 # Hermitian matrix
+        push!(v,max(eig(s)[1]))    # Largest Eigenvalue
+    end
+    v=n^(1/6)*(v-2*sqrt(n))        # normalized eigenvalues
 
-count=hist(v,grid)/(t*dx)
+    grid=-5:dx:2
+    count=hist(v,grid)[2]/(t*dx)
+    (grid, count)
+end
+(grid, count) = largeeig_experiment(n,t,dx)
 
 ## Theory
 (t, f) = tracywidom(5., -8.)
 
 ## Plot
 using Winston
 p = FramedPlot()
-h = Histogram(count, dx)
-h.x0 = grid[1]
+h = Histogram(count, step(grid))
+h.x0 = first(grid)
 add(p, h)
-add(p, Curve(t+5, f, "linewidth", 2, "color", "blue")) 
-file(p, "largeeig.png")
+add(p, Curve(t, f, "linewidth", 2, "color", "blue"))
+if isinteractive()
+    Winston.display(p)
+else
+    file(p, "largeeig.png")
+end
 

demos/book/1/semicircle.jl

@@ -4,28 +4,39 @@
 #Experiment:  Sample random symmetric Gaussian matrices
 #Plot:        Histogram of the eigenvalues
 #Theory:      Semicircle as n->infinity
+#See also:    RandomMatrices.SemiCircle Distribution
+
 ## Parameters
 n=1000       # matrix size
 t=10         # trials
-v=[]         # eigenvalue samples
 dx=.2        # binsize
-grid=[-2.4:dx:2.4]
-## Experiment
-for i = 1:t
-    a=randn(n,n)      # n by n matrix of random Gaussians
-    s=(a+a')/sqrt(2*n)# symmetrize and normalize matrix
-    v=[v;eig(s)[1]]   # eigenvalues
+function semicircle_experiment(n,t,dx)
+    ## Experiment
+    v=Float64[] # eigenvalue samples
+    for i = 1:t
+        a=randn(n,n)      # n by n matrix of random Gaussians
+        s=(a+a')/sqrt(2*n)# symmetrize and normalize matrix
+        append!(v,eig(s)[1])# eigenvalues
+    end
+    grid=-2.4:dx:2.4
+    count=hist(v,grid)[2]/(n*t*dx)
+    (grid, count)
 end
-count=hist(grid)/(n*t*dx)
+grid,count = semicircle_experiment(n,t,dx)
+
 ## Theory
 x=[-2:dx:2]
 y=sqrt(4-x.^2)/(2*pi)
 
 ## Plot
 using Winston
 p = FramedPlot()
-h = Histogram(count, dx)
-h.x0 = grid[1]
+h = Histogram(count, step(grid))
+h.x0 = first(grid)
 add(p, h)
-add(p, Curve(x+2.4+dx/2, y, "color", "blue", "linewidth", 2)) 
-file(p, "semicircle.png")
+add(p, Curve(x, y, "color", "blue", "linewidth", 2)) 
+if isinteractive()
+    Winston.display(p)
+else
+    file(p, "semicircle.png")
+end

demos/book/1/tracywidom.jl

@@ -2,6 +2,7 @@
 #Algorithm 1.3 of Random Eigenvalues by Alan Edelman
 
 #Theory:      Compute and plot the Tracy-Widom distribution
+#See also:    RandomMatrices.TracyWidom(t) for tn=-8
 
 ##Parameters
 t0=5.         # right endpoint
@@ -28,5 +29,9 @@ end
 using Winston
 p = FramedPlot()
 add(p, Curve(t, f2, "linewidth", 2))
-file(p, "tracywidom.png")
+if isinteractive()
+    Winston.display(p)
+else
+    file(p, "tracywidom.png")
+end
 

demos/book/2/bellcurve2.jl

@@ -0,0 +1,35 @@
+#bellcurve2.jl
+#Algorithm 2.1 of Random Eigenvalues by Alan Edelman
+
+#Experiment:    Generate random samples from the normal distribution
+#Plot:          Histogram of random samples
+#Theory:        The normal distribution curve
+
+
+## Experiment
+t = 10000     # Trials
+dx = 0.25       # binsize
+v = randn(t,2)  # samples
+x = -4:dx:4   # range
+
+## Plot
+using Winston
+module Plot3D
+    using Winston
+    # Warning: the 3D interface is very new at the time
+    # of writing this and is likely to not be available
+    # and/or have changed
+    if Winston.output_surface == :gtk
+        include(Pkg.dir("Winston","src","canvas3d_gtk.jl"))
+    else
+        include(Pkg.dir("Winston","src","canvas3d.jl"))
+    end
+end
+edg1, edg2, count = hist2(v, x-dx/2)
+#cntrs = ((first(x)+dx/2), (last(x)-dx/2), length(x)-1)
+#x,y = meshgrid(cntrs, cntrs)
+cntrs = (first(x)+dx/2):dx:(last(x)-dx/2)
+x = Float64[x for x = cntrs, y = cntrs]
+y = Float64[y for x = cntrs, y = cntrs]
+Plot3D.plot3d(Plot3D.surf(x,count./(t*dx^2) *50,y))
+Plot3D.plot3d(Plot3D.surf(x,exp(-(x.^2+y.^2)/2)/(2*pi) *50,y))

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/hist2d.jl

@@ -0,0 +1,4 @@
+#hist2d.jl
+#Algorithm 2.2 of Random Eigenvalues by Alan Edelman
+
+#see Base.hist2 for a simple implementation

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,43 @@
+#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
+require("patiencesort.jl")
+function tracywidomlis(t, n, dx)
+    ## single processor: 1x speed
+    #v = [patiencesort(randperm(n)) for i = 1:t]
+
+    ## simple parallelism: Nx speed * 130%
+    v = pmap((i)->patiencesort(randperm(n)), 1:t)
+
+    ## maximum parallelism: Nx speed
+    #grouped = floor(t/(nprocs()-2))
+    #println(grouped)
+    #v = vcat(pmap((i)->[patiencesort(randperm(n)) for j = 1:grouped], 1:t/grouped)...)
+
+    w = (v-2sqrt(n))/n^(1/6)
+    return hist(w, -5:dx:2)
+end
+@time grid, count = tracywidomlis(t, n, dx)
+
+## Plot
+using Winston
+p = FramedPlot()
+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