(* Mathematica example for random matrix exercise *)

(* Generating M members of the GOE ensemble of NxN matrices, and 
   calculating the difference between the eigenvalues 
   in the middle of the range *)

(* Load packages for Gaussian random numbers and histograms *)

<< Statistics`NormalDistribution`
<< Graphics`Graphics`


(* GENERATING RANDOM GOE MATRICES *)

(* Generating M members of the GOE ensemble of NxN matrices, and 
   calculating the difference between the eigenvalues 
   in the middle of the range *)

(* Increase M, N as appropriate *)

M = 20;
Size = 4;

(* Generate M random SizexSize arrays of Gaussian (normal) random
   numbers; symmetrize them; find their sorted eigenvalues,
   keep a list of the difference between the two middle ones *)

diffs = Table[0, {i, 1, M}];
M11 = Table[0, {i, 1, M}];
For[i = 1, i <= M, i++,
  {Mat = RandomArray[NormalDistribution[], {Size, Size}];
    MatSym = Mat + Transpose[Mat];
    eigs = Sort[Eigenvalues[MatSym]];
    diffs[[i]] = eigs[[Size/2 + 1]] - eigs[[Size/2]];
    M11[[i]] = MatSym[[1]][[1]];
    }]

(* You can also do this with functional programming: it's considered
   more elegant in Mathematica world *)

Symmetrize[m_] := m + Transpose[m];
EigSorted[m_] := Sort[Eigenvalues[m]]
EigDiffs[eigs_] := eigs[[Size/2 + 1]] - eigs[[Size/2]]
diffs = Map[EigDiffs[EigSorted[Symmetrize[#]]] &, 
      RandomArray[NormalDistribution[], {M, Size, Size}]];
M11 = Map[Symmetrize[#][[1]][[1]] &, 
      RandomArray[NormalDistribution[], {M, Size, Size}]];

(* Divide out by mean value of the splittings *)
diffAve = Mean[diffs]

(* Histogram generates a histogram of diffs. HistogramScale normalizes
   the integral to one. *)

Histogram[diffs/diffAve, HistogramScale->1]


(* PLOTTING THEORY OVER HISTOGRAM *)

(* Demonstrated with histogram for diagonal element
   First plot histogram for M11; store as "histGauss" for later
   plotting with theory *)

histGauss = Histogram[M11, HistogramScale->1]


(* Now, plot expected Gaussian fit *)

sigma = 2;
theoryGauss = Plot[(1/Sqrt[2 Pi sigma^2]) Exp[-x^2/(2 sigma^2)], {x, -5, 5}]

(* and show the comparison *)

Show[histGauss, theoryGauss]


(* GENERATING RANDOM +-1 MATRICES *)

(* Symmetric matrix with integer values +-1 with 50/50 probability *)

(* Start with random matrix of +-1, no symmetry *)

MatPM = Table[Random[Integer, {0, 1}]*2 - 1, {i, 1, Size}, {j, 1, Size}];

(* Symmetrize it (not very elegant) *)

MatPMSym = Table[If[i >= j, MatPM[[i]][[j]], MatPM[[j]][[i]]], 
			{i, 1, Size}, {j, 1, Size}];

(* Take a look at it *)
MatrixForm[MatPMSym]