Added code from Chapter 1 of Alan's book

Author: jiahao

demos/.gitignore

@@ -0,0 +1,2 @@
+*.png #Ignore output plots
+

demos/README.md

@@ -0,0 +1,2 @@
+This folder contains code snippets from the upcoming book "Random Eigenvalues" by Alan Edelman.
+

demos/book/1/bellcurve.jl

@@ -0,0 +1,28 @@
+#bellcurve.jl
+#Algorithm 1.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=1000000
+dx=.2
+v=randn(t)
+count=hist(v,[-4-dx/2:dx:4+dx/2])
+count/=sum(count)*dx
+
+## Theory
+x=[-4:dx:4]
+y=exp(-x.^2/2)/sqrt(2*pi)
+
+## Plot
+#XXX Currently Histogram() must start plotting from 0.
+#For now, add offset to theoretical curve
+#Filed: https://github.com/nolta/Winston.jl/issues/28
+#cjh 2013-03-16
+using Winston
+p = FramedPlot()
+add(p, Histogram(count, dx))
+add(p, Curve(x+4+dx/2, y, "color", "blue", "linewidth", 2)) 
+file(p, "bellcurve.png")

demos/book/1/gradient.jl

@@ -0,0 +1,15 @@
+#Calculate numerical derivative using:
+#forward first-order finite difference for first element
+#backward first-order finite difference for last element
+#central first-order finite difference for all intermediate elements
+function gradient(y::Vector{Float64}, x::Vector{Float64})
+    dy = zeros(length(y))
+    n = length(dy)
+    dy[1] = (y[2]-y[1])/(x[2]-x[1])
+    dy[n] = (y[n]-y[n-1])/(x[n]-x[n-1])
+    for i=2:n-1
+        dy[i] = (y[i+1]-y[i-1])/(x[i+1]-x[i-1])
+    end
+    return dy
+end
+

demos/book/1/largeeig.jl

@@ -0,0 +1,39 @@
+#largeeig.png
+#Algorithm 1.4 of Random Eigenvalues by Alan Edelman
+
+#Experiment:  Largest eigenvalue of random Hermitian matrices
+#Plot:        Histogram of the normalized largest eigenvalues
+#Theory:      Tracy-Widom as n->infinity
+
+include("tracywidom.jl")
+
+## Parameters
+n=100        # matrix size
+t=5000       # trials
+v=[]         # eigenvalue samples
+dx=.2        # binsize
+## 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
+
+count=hist(v,-5:dx:2)
+count/=t*dx
+
+## Theory
+(t, f) = tracywidom(5., -8.)
+
+## Plot
+#XXX Currently Histogram() must start plotting from 0.
+#For now, add offset to theoretical curve
+#Filed: https://github.com/nolta/Winston.jl/issues/28
+#cjh 2013-03-16
+using Winston
+p = FramedPlot()
+add(p, Histogram(count, dx))
+add(p, Curve(t+5, f, "linewidth", 2, "color", "blue")) 
+file(p, "largeeig.png")
+

demos/book/1/semicircle.jl

@@ -0,0 +1,29 @@
+#semicircle.m
+#Algorithm 1.2 of Random Eigenvalues by Alan Edelman
+
+#Experiment:  Sample random symmetric Gaussian matrices
+#Plot:        Histogram of the eigenvalues
+#Theory:      Semicircle as n->infinity
+## Parameters
+n=1000       # matrix size
+t=10         # trials
+v=[]         # eigenvalue samples
+dx=.2        # binsize
+## 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
+end
+count=hist(v,-2.4:dx:2.4)
+count/=n*t*dx
+## Theory
+x=[-2:dx:2]
+y=sqrt(4-x.^2)/(2*pi)
+
+## Plot
+using Winston
+p = FramedPlot()
+add(p, Histogram(count, dx))
+add(p, Curve(x+2.4+dx/2, y, "color", "blue", "linewidth", 2)) 
+file(p, "semicircle.png")

demos/book/1/tracywidom.jl

@@ -0,0 +1,32 @@
+#tracywidom.m
+#Algorithm 1.3 of Random Eigenvalues by Alan Edelman
+
+#Theory:      Compute and plot the Tracy-Widom distribution
+
+##Parameters
+t0=5.         # right endpoint
+tn=-8.        # left  endpoint
+
+##Theory: The differential equation solver
+#Compute Tracy-Widom distribution
+using ODE
+include("gradient.jl")
+function tracywidom(t0::Float64, tn::Float64)
+    function deq(t::Float64, y::Vector{Float64})
+        yout = [y[2]; t*y[1]+2y[1]^3; y[4]; y[1]^2]
+    end
+    
+    y0=[airy(t0); airy(1, t0); 0; airy(t0)^2] # Initial conditions
+    (t, y)=ode23(deq, [t0, tn], y0)           # Solve the ODE
+    F2=exp(-y[:,3])                           # the distribution
+    f2=gradient(F2, t)                        # the density
+    return (t, f2)
+end
+(t, f2) = tracywidom(t0, tn)
+
+## Plot
+using Winston
+p = FramedPlot()
+add(p, Curve(t, f2, "linewidth", 2))
+file(p, "tracywidom.png")
+