Tracy_Widom.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'RandomMatrix'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-2f762be5acda>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0msys\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'../../'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0mRandomMatrix\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mRM\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     13\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'matplotlib'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'inline'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'RandomMatrix'"
     ]
    }
   ],
   "source": [
    "from sympy import *\n",
    "import numpy as np\n",
    "from sympy import symbols, Function, integrate, diff\n",
    "from scipy.special import airy\n",
    "from scipy.integrate import solve_ivp, quad\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import random\n",
    "\n",
    "import sys\n",
    "sys.path.append('../../')\n",
    "import RandomMatrix as RM\n",
    "\n",
    "%matplotlib inline\n",
    "sns.set()\n",
    "\n",
    "np.random.seed(140)\n",
    "random.seed(140)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defining the Tracy Widom distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Tracy-Widom distribution was defined by Craig A. Tracy and Harold Widom in their paper [*On Orthogonal and Symplectic Matrix Ensembles*](https://arxiv.org/abs/solv-int/9509007). It is an important distribution for studying the distribution of the fluctuations of the maximum eigenvalue of a Gaussian Ensemble. Each Gaussian ensemble has a corresponding Tracy Widom distribution. We will denote the Tracy Widom distribution for Gaussian Ensemble with Dyson index $\\beta$ as $TW_\\beta$. Let us study the distribution before its application.\n",
    "\n",
    "Suppose $q$ is the solution to the Painlevé differential equation\n",
    "\n",
    "$$q'' = sq + 2q^3 $$\n",
    "\n",
    "with the boundry condition $q \\sim Ai(s)$ as $s \\rightarrow \\infty$. Define $I(t)$ and $J(t)$ to be the following\n",
    "\n",
    "$$I(t) = - \\int_{t}^{\\infty} (x-t) q^2(x) dx$$\n",
    "\n",
    "$$J(t) = \\int_{t}^{\\infty} q(x) dx$$\n",
    "\n",
    "If we define $F_\\beta$ to be the CDF of $TW_\\beta$ then \n",
    "\n",
    "$$F_2(t) = e^{I(t)}$$\n",
    "\n",
    "$$F_1(t) = \\sqrt{F_2(t) e^{-J(t)}} $$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Density of Tracy Widom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To solve for the density of the Tracy Widom we would need to solve for a closed form solution to the Painlevé differential equation. We will instead approximate the curve using the Scipy's ODE solver. In order to define an initial condition, let $t_0$ be a large value so that the condition $q(t_0) = \\text{Ai}(t_0)$ approximates our asymptotic boundry condition well. Using the Fundamental Theorem of Calculus, we can differentiate $I$ and $J$ to make them a solution to a differential equation our ODE solver can solve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\int\\limits_{t}^{\\infty} \\left(- t + x\\right) q^{2}{\\left(x \\right)}\\, dx$"
      ],
      "text/plain": [
       "Integral((-t + x)*q(x)**2, (x, t, oo))"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Definition of I(t)\n",
    "t, x = symbols(\"t x\")\n",
    "q = Function(\"q\")\n",
    "I = integrate((x-t)*q(x)**2, (x, t, oo))\n",
    "I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\int\\limits_{t}^{\\infty} - q^{2}{\\left(x \\right)}\\, dx$"
      ],
      "text/plain": [
       "Integral(-q(x)**2, (x, t, oo))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# I'(t)\n",
    "I_prime = diff(integrate((x-t)*q(x)**2, (x, t, oo)), t)\n",
    "I_prime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle q^{2}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "q(t)**2"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# I''(t)\n",
    "I_prime2 = diff(integrate((x-t)*q(x)**2, (x, t, oo)), t, 2)\n",
    "I_prime2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus we can solve for $I(t)$ by solving for the differential equation\n",
    "\n",
    "$$ \\frac{d}{dt} \\begin{pmatrix}\n",
    " \\text{I}(t)\\\\\n",
    "\\text{I}'(t)\n",
    "\\end{pmatrix} =\n",
    "\\begin{pmatrix}\n",
    " \\text{I}'(t)\\\\\n",
    "q^2\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "such that\n",
    "\n",
    "$$I(t_0) = - \\int_{t_0}^{\\infty} (x-t_0) q^2(x) dx $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\int\\limits_{t}^{\\infty} q{\\left(x \\right)}\\, dx$"
      ],
      "text/plain": [
       "Integral(q(x), (x, t, oo))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J = integrate(q(x), (x, t, oo))\n",
    "J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - q{\\left(t \\right)}$"
      ],
      "text/plain": [
       "-q(t)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# J'(t)\n",
    "diff(J)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus we can solve for $J$ by solving\n",
    "\n",
    "$$ \\frac{dJ}{dt} = -q(t)$$\n",
    "\n",
    "such that \n",
    "\n",
    "$$J(t_0) = \\int_{t_0}^{\\infty} q(x) dx$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to improve accuracy, we will use the solver to simultaneously find the values of $q(t)$, $I(t)$ and $J(t)$ in one system of differential equations. This strategy was suggested by Alan Edelman and Per-Olof Persson in their paper [*Numerical Methods for Eigenvalue Distributions of Random Matrices*](https://arxiv.org/abs/math-ph/0501068). Note that we must now slightly change our initial conditions for $I(t)$ and $J(t)$ to avoid using $q$ (since we won't have an approximation for it yet). Luckly, we decided to choose $t_0$ large enough so that $q(t) \\approx \\text{Ai}(t)$ for $t \\geq t_0$, so we may replace $q$ in our intial conditions with $\\text{Ai}$. The system of differential equations becomes \n",
    "\n",
    "$$\\frac{dy}{dt} = \\begin{pmatrix}\n",
    "q' \\\\\n",
    "tq+2q^3 \\\\\n",
    "I' \\\\\n",
    "q^2 \\\\\n",
    "-q\n",
    "\\end{pmatrix}$$\n",
    "\n",
    "where\n",
    "\n",
    "$$y = \\begin{pmatrix}\n",
    "q \\\\\n",
    "q' \\\\\n",
    "I \\\\\n",
    "I' \\\\\n",
    "J\n",
    "\\end{pmatrix}$$\n",
    "\n",
    "such that\n",
    "\n",
    "$$y(t_0) = \\begin{pmatrix}\n",
    "\\text{Ai}(t_0) \\\\\n",
    "\\text{Ai}'(t_0) \\\\\n",
    "\\int_{t_0}^{\\infty} (x-t_0) \\text{Ai}^2(x) dx \\\\\n",
    "\\text{Ai}(t_0)^2 \\\\\n",
    "\\int_{t_0}^{\\infty} \\text{Ai}(x) dx\n",
    "\\end{pmatrix} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we define the function we will pass into the ODE solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t, y):\n",
    "    d0 = y[1]\n",
    "    d1 = t*y[0]+2*y[0]**3\n",
    "    d2 = y[3]\n",
    "    d3 = y[0]**2\n",
    "    d4 = -y[0]\n",
    "    return np.array([d0, d1, d2, d3, d4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0 = 4\n",
    "tf = -6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the initial condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0_0 = airy(t0)[0]\n",
    "y0_1 = airy(t0)[1]\n",
    "y0_2 = quad(lambda x: (x-t0)*airy(x)[0]**2, t0, np.inf)[0]\n",
    "y0_3 = airy(t0)[0]\n",
    "y0_4 = quad(lambda x: airy(x)[0], t0, np.inf)[0]\n",
    "y0 = np.array([y0_0, y0_1, y0_2, y0_3, y0_4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve the ODE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = solve_ivp(f, (t0, tf) ,y0, max_step=.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the density of $TW_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "F2 = np.exp(-sol.y[2])\n",
    "f2 = -sol.y[3] * F2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the density of $TW_1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "F1 = np.sqrt(F2 * np.exp(-sol.y[4]))\n",
    "f1 = 1/(2*F1)*(f2+sol.y[0]*F2)*np.exp(-sol.y[4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the densities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEXCAYAAACtTzM+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXgV1f348ffMXbNvZCMhIQGSACHsi4AsCoJIxB3rQm0rXazar7+2T/laa621td9qa9VWrdrWulVbilK0ICiLyr4atoRsBBKy78vN3eb8/ghEkBCy3C035/U8eeDemTnnM/dCPnPOmTlHEUIIJEmSpEFP9XYAkiRJkm+QCUGSJEkCZEKQJEmSzpIJQZIkSQJkQpAkSZLOkglBkiRJAmRCkDzgiSeeYNmyZSxbtozMzEwWLVrU+bq9vd0jMbz88svcd999na+FEFx55ZUsX778gv2uu+46tm3bxrPPPsv7779/UTl1dXWkp6e7JUan08n3vvc9Fi1axJtvvtnncsrKyrj33nvJzs5myZIlfPDBB706/tlnn+Xxxx+/5PacnBxuv/12li1bRnZ2NmvXru1zrJKPEZLkQfPnzxc5OTker/f48eNi6tSpwul0CiGEOHz4sFi+fLmYMWOGqKmpEUIIUVlZKcaPHy8sFssly6mtrRVpaWluibGsrExkZmYKh8PRr3K+/e1vi7/97W9CCCEqKirElClTREtLy2WPKy8vFw888IDIysoSv/jFL7rcR9M0MXfuXLF9+/bOY2bMmCGKi4v7FbPkG/TeTkiSlJmZydVXX01ubi5PP/00eXl5vPvuu9jtdhobG1m5ciV33HEHAH/+859577330Ov1JCcn85vf/IYHH3yQxYsXd17tv/jii9TX1/Pwww931pGRkYHZbCYvL4/Ro0ezZcsW5s2bR2FhIZs3b+bWW29l165dTJ8+HbPZzKpVqxg1ahTf+ta32LhxI8888wwBAQFkZmZeEPuf/vQnPvzwQ3Q6HSkpKfzsZz8jOjqau+++m7Fjx7Jr1y5qa2tZsWIFtbW17NmzB4vFwh/+8IcLWhotLS3ce++9OBwObrrpJp5//nmqqqr47W9/i8ViwWAw8D//8z/MmTOHNWvWsHr1aiwWC8HBwbzxxhsXxFRUVMT3v/99AGpqaggNDcVoNFJQUMAPf/jDiz7/FStWcPPNN7N69WomT55MamoqTU1NXX5XNpuN73//+8ycOROAuLg4IiIiqKioYPjw4b385iWf4+2MJA0uXbUQ0tLSxHvvvSeEEKKlpUXcdtttoq6uTgghxMGDB8WECROEEEJ8/PHH4pprrhENDQ1CCCF+/etfixdeeEFs2rRJ3HzzzUIIIZxOp5g/f74oLCy8qO6HH36488r5xhtvFMePHxcffPCB+O53vyuEEGLVqlXirbfeEkII8ZOf/ES8+uqrorq6WkyePFnk5+cLIYR46aWXOlsIq1evFsuXLxetra1CCCGee+458c1vflMIIcRdd90l7r//fiGEEIcOHRJpaWnik08+EUII8atf/Uo88sgjF8V3+vTpznOtq6sTV1xxhTh06JAQQogTJ06IadOmiVOnTol///vfYurUqaK5ubnbz/jOO+8UU6dOFRs2bOhyv0t57rnnLtlC+Kp33nlHzJ07t9tWlTRwyBaC5BOmTJkCQFBQEC+99BLbtm3j5MmT5Obm0tbWBsDOnTtZvHgxYWFhAPzv//4v0NH3/sQTT5Cbm0tlZSWJiYmkpqZeVMecOXNYu3Yt1157LQ0NDWRkZDB06FAee+wxnE4nu3fv5oEHHrjgmP3795OWlsbIkSMBWL58Ob///e8B+PTTT7npppsIDAwEOq60X3rpJWw2GwALFy4EYNiwYQBceeWVACQlJbFnz55uP4+cnBySkpIYP348AKNGjWLSpEns2bMHRVFIT08nODi42zLefPNNysrK+OY3v4nZbCYhIaHbFkJvvfzyy7z++uu8+uqrmM3mXh8v+R6ZECSfcO6XakVFBcuXL+e2225j8uTJLF68mC1btgCg0+lQFKXzmKamJpqamkhMTOT2229n9erVVFVVcfvtt3dZx8yZM3n88cfZsmULc+bMASA0NJT09HQ2btxIUFAQQ4cOveAYRVEQ5033pdd/+V9GfGUaME3TcDgcna+NRuMF2w0GQ48/D03TLnpPCIHD4cBgMHR+XpdisVgASEhIIDs7m23btvHoo4+6ZADYZrOxatUqCgoKeOedd0hMTOx3mZJvkHcZST7lyJEjREZGct9993HllVd2JgOn08nMmTPZtGkTLS0tADz//PO89tprANx66618/PHHHD16tPPK/KtCQkJISUnh7bffZv78+Z3vz5s3j5deeom5c+dedMyUKVMoKCggNzcXgDVr1nRumz17NmvWrOlswbzxxhtMnTr1okTQF+PHj6e4uJicnBwA8vPz2bt3L9OmTevR8WvWrEHTNGw2G3v27HHpnVEPPvggLS0tMhn4IdlCkHzKrFmzWL16NYsXLyYgIICsrCwiIyMpKSlh7ty5FBQU8LWvfQ2AkSNH8stf/hKAqKgoMjMzGTFiRLdX4nPmzOGFF15gxowZne/NmzePp556ip/97GcX7R8ZGcnTTz/Nj370IwwGA1OnTu3cdsstt1BeXs6tt96KpmkkJyfz9NNPu+RziIyM5Nlnn+WXv/wl7e3tKIrCk08+SUpKCgcPHrzs8bGxsdxxxx00NDQwbdq0PnUJnVNZWcm3v/1tXn75ZUpLS9myZQvDhw/v/B4AfvSjH3V2iUkDlyK+2u6VpAGorq6OW265hbfeeov4+Hhvh+NVV111Fc8++yzjxo1zWZk//vGPefjhh4mIiHBZmZLvkV1G0oD3z3/+kyVLlrBixYpBnwzcwWKxMHv2bJkMBgHZQpAkSZIA2UKQJEmSzpIJQZIkSQJkQpAkSZLOkglBkiRJAgb4cwj19a1oWu/HxKOigqmtbXFDRL5LnvPgIM95cOjrOauqQkRE0CW3D+iEoGmiTwnh3LGDjTznwUGe8+DgjnOWXUaSJEkSIBOCJEmSdNaA7jKSJEnqLSEE9fXV2GztwMDsaqqqUrucEfccnU5PcHA4AQGXHi/oikwIkiQNKi0tjSiKQmxsIooyMDtJ9HoVh6PrhCCEwG630dBQDdCrpDAwPw1JkqQ+slhaCAkJH7DJ4HIURcFoNBEeHk1LS0OvjvXPT0SSJOkSNM2JTuf/nSMGgxGn03H5Hc/j/5+KNKi1WOz85/Nith+pIGFIEBnJEWSmRDIqMeyC1dekwWUwfPd9OUeZECS/pGmC97YW8I+NebTbHExOi6au2cqHO0/ywY6TLJ6exK3zRgyKXwyS7/vgg/f517/eQdM0xo+fyA9/uOqS/zbb2tr4v//7JSdPFiOEYMWKb7JgwSKXxCETguSX1nxaxH93lTAuNYrb5o8gIbpjQfq2dgertxawYfcpzEYd189K8XKk0mBXXV3FX//6Cm+88S6BgUE8+uj/snv3TmbMmNnl/m+++RqxsXH84hdPUl9fxz333MGkSVOIjIzqdywyIUh+59jJOtbvKuGa6cncPn/EBdsCzXruWpSOzaHx/mfFmAw6Fk1L8lKkkgRnzpQxdGgCQUHBWCwWWlqaCQ+P4J577rho38cf/zUTJ04iJaXjQiYiIpLQ0FBqa2tlQpCkr2pqs/HKB8eIiwpk5bJMmpssF+2jKgrfWJKBze7k3c0FxEYGMmHkEC9EK3nb9sPlfJ5T7payZ2fFM2vc5VfwczqdAGzd+gm//vXjLFy4iLS0dF577e0u909KGt552+knn2zEbreTkpLqkphlQpD8hhCCv354nFaLnYduHY/ZpKf5EvvqVJVvXz+Wsr/u4d1P8slMiUSvkzfdSd4zb97VzJ49lz/84Sn++MdnOHBg/0X7PP74r0lKGg7A5s0f89xzv+Ppp59Hr3fNr3KZECS/setoJTmFtdyxYBRJsSGX3V+vU1l+1Uj+8K8cthwsY+GUYR6IUvIls8b17Cre3azWdgD0ej0337ycRx9dxRtv/POS+//zn+/w5pt/5/e//xMjRox0WRzykkjyC0II/rurhIToIK6enNjj48alRjF2eAT/+byYFovdjRFK0qUVFhZQVFQAwIED+0hJGXHJfT/9dCvvvPMWL774F5cmA5AJQfITh4tqKatp5drpSb26lVRRFJZfNYo2q4MPdpx0X4CS1I0RI0byxhuv8c1v3smWLR9z330/uOS+f/nLn7FarfzkJ/+Pe+65g3vuuYPc3GMuiUN2GUl+4b+7ThEZamLa6NheH5sYE8yVWUP5ZH8p8yclEBsR6IYIJenSTCYzP//5Ez3a9+9//0e3cxn1h2whSANeYVkjJ043cM3UpD4PDN94ZQqKorBp72kXRydJA4dMCNKAt2H3KYLMeuaM7/vgYFiwiakZ0ew4UkG7rXfzv0hSf0yaNIU//vFlb4cByIQgDXAVdW0cOFHN/EkJmI396wGdPymRdpuTXUcrXRSdJA0sMiFIA9pnX5xBVRWuntz/W0ZHDA1lWEwwWw6WIcTAXDhFkvpDJgRpwNKEYM/xSsamRBIWZOx3eYqiMH9SAqerWig80+SCCCVpYJEJQRqwCssaqW2yMr0PdxZdyowxsZiNOrYcKHVZmZI0UMiEIA1Ye45VYdCrTBjlunmIzEY9MzPj2JtbRXObzWXlStJAIBOCNCA5NY29uZWMHxFFgMm1j9PMn5iAwynYeaTCpeVK0qV88MH7fP3rt3P33bfx9NNP9mgMq7CwgLvuus2lcbg1Iaxbt44lS5awcOFC3nrrrUvut3XrVq666ip3hiL5mdySBpra7Ewf47ruonMSooNJig1mb26Vy8uWpK86tx7CCy+8yuuvv0tjYyO7d+/s9pj//vcDfvjDB2hvv3g23/5w25PKlZWVPPPMM6xZswaj0cjtt9/O9OnTGTnywrk3ampq+L//+z93hSH5qd3HKzEbdWSN6P8c8F2ZmhHDv7cVUdNoYUhYgFvqkLzPfmI79rxP3VK2IX0OhrRZl92vt+shREYO4dNPt/LYY7/iiSd+7tKY3dZC2LFjBzNmzCA8PJzAwEAWLVrEhg0bLtrvkUce4f7773dXGJIfsjs09udVMyktGoNe55Y6pmTEALAvt9ot5UvSOeevh7Bs2WKGDk3oXA/hqz9JScMJDg7mN795mtjYOJfH4rYWQlVVFdHR0Z2vY2JiyMnJuWCf119/nTFjxjB+/Hh3hSH5oSPFtVisDrd0F50TGxFIcmwI+/KqWDxdrqjmrwxps3p0Fe8JvV0PwR3clhC6GhQ5fxbKEydOsHHjRl577TUqKvo2eBcVFdzn+KKjLz9fvr/xl3PO3VzQMVXFlMvPXdSfc543ZRh///AYmk5HbOTAmfDOX77n3ujNOVdVqej1vnM/jU6nYrNZ0etV9Hojy5d/jZ/+9Ce8/fa/Lnsc0O25qKraq8/GbQkhNjaWffv2db6uqqoiJiam8/WGDRuorq7m5ptvxm63U1VVxR133MHbb3e9bFxXamtb0LTeP1EaHR1CdfWl1tLyT/5yzkII9h2vZHRyBPV1rd3u299zHj0sDICNO4oHTCvBX77n3ujtOWua5paZQvvK6dQoKMjnxIkTpKaOZM+ePQwfntptjHq9itPZsb27/TRNu+CzUVWl2wtpt6XJmTNnsnPnTurq6rBYLGzcuJE5c+Z0bn/wwQf56KOPWLt2LS+//DIxMTG9SgbS4FRa3Up9s5Vxqe4ZTD5fTHgAw+NC2Jsr5zaS3Ks36yG4k1tbCA899BArVqzAbrdzyy23kJWVxcqVK3nwwQcZN26cu6qW/NjholoAMj2QEACmjo7hX1sKqW6wEB0u7zaS3KM36yGcEx8/lNWr17k0DrcukJOdnU12dvYF773yyisX7ZeYmMjmzZvdGYrkJ44U1TIsJpiIEJNH6pua3pEQ9uVVce30ZI/UKUne4jsjK5J0GRarg/zSRjJTIz1W55DwABKjg8kpqPVYndLgItdDkKQ+OHayHqcmyPJQd9E5WSOiyC9tpK3d7tF6JfcZDNOb9+UcZUKQBozDRbUEmHSMSAjzaL1ZI6LQhODoyXqP1iu5h6rqcDr9f1U8u92GTte7UQGZEKQBQQjB4aJaxiRH9nnd5L4akRBKkFlPTkGNR+uV3CMgIJjm5gaE8J1bT11JCIHNZqWhoZrg4PBeHevWQWVJcpUzNWdvN53t2e4iAJ2qkpkaxeGiWjQhUM97wFIaeIKDw6ivr6ayshQYmF1HqqqiaZdOaDqdnpCQCAICgnpVrkwI0oBwuKgOgMwUzw0ony9rRBS7j1VSUtFMSnyoV2KQXENRFCIjYy6/ow9z1wOIsstIGhByT9UTFxlIZKjZK/WPS41CAb6Q3UaSH5MJQfJ5mibIL20gPal3/aGuFBxgIDUhlJxCefup5L9kQpB83qmqZixWp1cTAkDWiCGcrGimsVUurSn5J5kQJJ+Xd6oBgPRhEV6NY/zZxXgOy1aC5KdkQpB8Xt6pBmIjAjw2XcWlDIsJJizYyNGTdV6NQ5LcRSYEyadpmuDEae+OH5yjKApjkiM5drIObRA86SoNPjIhSD6ttLqFNquD9CTvdhedM2Z4BM1tdkqrWrwdiiS5nEwIkk/L7Rw/8H4LAWDM8I7nII7JaSwkPyQTguTT8k7VEx1u9trzB18VEWIiPiqQYyVyHEHyPzIhSD5LE+fGD3yju+icMcMjOXG6AbsPLcMoSa4gE4Lks8qqW2ltd/hMd9E5Y4ZHYLNrFJY1ejsUSXIpmRAkn5V3qqOf3hfuMDpfRlIEqqLIbiPJ78iEIPms/NJGokJNDAnzrbWMA0x6UoeGyoFlye/IhCD5rIKyRkYm+lbr4JwxwyMoLm+iVa6iJvkRmRAkn1TX1E59s5WRHl4drafGDI9ECMgtafB2KJLkMjIhSD6p4OyA7YgE31x7IHVoKCaDTo4jSH5FJgTJJxWUNmI0qCRGB3s7lC7pdSqjEsM6J96TJH8gE4LkkwrKGkmND/X4+sm9kZEcwZmaVprkdNiSn/Dd/23SoGW1OTlV2cIIHx0/OOfc7bC5p+TdRpJ/kAlB8jknK5rQhPDZAeVzkmNDMBl1sttI8hsyIUg+58sBZd9OCHqdSlpiuGwhSH5DJgTJ5xSUNhIfFUhwgMHboVxWRlI45bVtNLZYvR2KJPWbTAiSTxFCUHimyedbB+dkJHdMvJcru40kP6D3dgCSdL7KegstFrvXxw+05mq0ulK0xgq0phoUczBqeHzHT+QwFLXjWiopNhizUUfeqXqmj4n1asyS1F8yIUg+paC0Y/zAGwlBCIGz7Bi2nPU4S498ucEYAPZ2OLtsphISjTFzIYb0K9EZA0gbFi5bCJJfkAlB8ikFZY0EmvTERQV6tF5n7Wnat72KVlOCEhCKccpN6BPHoobGopiDEQ4bWlMVWu0p7Me3Yt35NtZ972GafD0Zw0aTU1hLfbOViBCTR+OWJFeSCUHyKUVnmkgZGoqqKB6pTwiBPe9TrNvfRDEGYprzDQyjZqLoLhzQVvRGdJGJ6CITMYyaibOqCOuBtVh3vcvU6AzWK+PIO13PjDFxHolbktxBJgTJZ7TbHJTVtDBx1HCP1Cecdto/fQ1H/nZ0CWMwz/8OamDPuqp0MakELPof7LnbYMfbrAo/yYE8HYxZ7OaoJcl95F1Gks8oqWhGiI6J49xNaE7aP3kJR/52jJOWEXDtj3qcDM5RFAXj6HkE3vRz7LogZlb/E8epQ26KWJLcz60JYd26dSxZsoSFCxfy1ltvXbR906ZNZGdnc91117Fq1SpsNjknzGBWdKYJgBQ3JwQhNNq3vorj5H5MM+/ENOXGzruG+kIXkUDe6Hs544igbePz2E8ecGG0kuQ5bksIlZWVPPPMM7z99tusXbuWd999l4KCgs7tbW1tPP744/ztb3/jww8/xGq18t5777krHGkAKDrTRHS4mdBAo9vqEEJg/fx1HAU7MU69GWPmQpeUOyIlgReaF9IeNJT2TX/CXrzfJeVKkie5LSHs2LGDGTNmEB4eTmBgIIsWLWLDhg2d2wMDA9m8eTNDhgyhra2N2tpaQkN9c+57yTOKyptIHere203tRzZiP74V44TrME3Mdlm5yXHBaPoAPom4FTV6OO2b/4yz+qTLypckT3BbQqiqqiI6OrrzdUxMDJWVlRfsYzAY2LZtG/Pnz6e+vp7Zs2e7KxzJx9U3W6lvtpIa776LAmdVEdbd/0SfPBHj1FtcWrZOVRmZGMaxUgsBi36AEhCCZeOzaG3y+QRp4HDbXUbi7EM851O6uJVw7ty57N69m9///vc89thj/O53v+txHVFRfV88JTo6pM/HDlS+fM755c0ATBob59I4z5XltLRQ9u6L6EMiSbj5f9AFuH7hnUkZsbyx/jjmyCGE3/4wZ/7+Uxyb/0T83Y+j6t3XDfZVvvw9u4s8Z9dwW0KIjY1l3759na+rqqqIiYnpfN3Q0MCRI0c6WwXZ2dk89NBDvaqjtrYFTbs48VxOdHQI1dXNvT5uIPP1cz6UW4lOVQg1qi6L89w5CyFo3/gcjuY6Aq//KXUtAlpc/1kkRgUAsPNQGZPTozHNW0n7pucpe/8FzHO/5fL6uuLr37M7yHPuOVVVur2QdluX0cyZM9m5cyd1dXVYLBY2btzInDlzOrcLIfjxj3/MmTNnAFi/fj2TJk1yVziSjys600RSbDAGvc7lZdvzPsVRchDTtNvQxaS6vPxzUuJDMepV8k53TIdtSJmMcWI29rzPsBftdVu9kuQqbksIsbGxPPTQQ6xYsYIbbriBpUuXkpWVxcqVKzl8+DARERH88pe/5Dvf+Q7XX389J0+e5Mc//rG7wpF8mKYJTlY0kxrv+gFlzdKEdfc/0cWnYxh3jcvLP59epzIiIYwT581rZJy8DDU6lfbPXkNrqXNr/ZLUX259Ujk7O5vs7Avv5HjllVc6/75gwQIWLFjgzhCkAaCsphWr3emWB9KsO/8B9nZMV369yzEsV0sfFs7az4tpa7cTaDagqHoCrvo2rf/+Oe3bXiVgyY9QFPk8qOSb5L9MyeuKznTMcOrqhNBW9EXH8wYTlqILH+rSsi8lPSkcAZw4O2srgBoWh2nmHTjLjmE/vNEjcUhSX8iEIHldcXkTQWY9MREBLitTOGzUbHgZJSwO44TrXFbu5aQODUWvUy7oNgIwpM9BnzwR6941aE1VHotHknpDJgTJ64rONJMSH+rSLh3b4Y9w1Fdgnr0CxYO3fBr0OlLjQzsHls9RFAXT7BWgqrR//nqXt2VLkrfJhCB5ldXmpKymhRQXPpCmWZqwHfqQwLRp6BPGuKzcnkpLiqCkogWL1XHB+2pQBKapN+MsPYKjcLfH45Kky5EJQfKqksqOGU5dOaGd7cB/wGEjcv6dLiuzN9KTwtGEoLCs8aJthjFXo0anYN35NqK9xQvRSdKlyYQgeVXnDKcuaiFojZXYj23BkDEX45BEl5TZWyOHhqFTFfJOXzxthaKqmK+8B9HegnXPv7wQnSRdmkwIkledrGgiKtRMWJBr+vmte1eDTo9x8jKXlNcXJqOO4XEh5F1inWXdkGQMmQux536Ks6bEw9FJ0qXJhCB5VdGZJlLiXTMni7OqCEfRXoxZi1EDw11SZl+lJYVTXN6E1e7scrtp0vUo5uCOriM5wCz5CJkQJK9parNR09jusvED6/73UUzBGLO8v4xl+rBwnJqgqItxBADFFIRxyk04y/NwFO/rch9J8jSZECSvOVneMX7giimvndXFOE/nYMhajGJ03fMMfTUqMRxFoctxhHMMGXNRI4dh3f0uwiFXC5S8TyYEyWuKzjShKJAc1/8uI9uB/4ApCOPYq10QWf8FmPQkxV56HAE6BphNM+9ANNdgO/yRB6OTpK7JhCB5zcmKZoYOCcJs7N+UWs6aEhwlBzFmXuMTrYNz0oeFU3imCbtDu+Q++qGj0Q+fjO3Qh2iWJg9GJ0kX61FCeOCBB9ixY4e7Y5EGESFEx4ByXP+7i2wH14EhAGOmb02UmD4sHIdTo7i8+1/0pmm3gMPWcR6S5EU9SgjXXHMNL7zwAosWLeIvf/kLDQ1yWUCpf2oa22mx2Ps9oOysK8VRvA/juIUopiAXRecao4Z13OmUd6q+2/3U8HgMGXOwH9ss5zmSvKpHCSE7O5s333yTF154gdraWm699VZ+/OMfk5OT4+74JD9V7KIBZdsX60Fvwpjp3rUO+iI4wEBidFC3A8vnGCctA0WHdd8aD0QmSV3r8RiCpmmUlJRw8uRJHA4HUVFRPPbYYzz11FPujE/yU0VnmtDrVBKi+35Vr7XU4SjYhSFjDorZ9Wsku0J6UgQFZY04nJceR4COeY6M467BUbBLPqwmeU2PEsIzzzzD3LlzefXVV1myZAkbN25k1apVvPnmm6xevdrdMUp+qLi8ieTYYPS6vt/XYD/6MaD5ZOvgnIykcGx2jZPll1//1jhhCZiC5JQWktf06PaOuro6XnnlFTIyMi54PzAwkN/97nduCUzyX05No6SimTkT+r5ojbBZsB3bgj5lKmpotAujc620s+MIuafqGZnY/RKhijEQ08SlWHe9i6M8D318uidClKROPbo8czqdFyWDBx54AIDZs2e7PirJr5VVt2JzaP0aP7Dnfgp2i088ldydkEBjxzjCZQaWzzGMuQolMBzbvjVySgvJ47ptIfz85z+nsrKS/fv3U1f35QLhDoeDoqIitwcn+aeicwPKfbzDSGgObEc2ootPRxeT6srQ3CI9KYLPcs7gcGqX7SJT9CaME5Zi3fEmzrJj6BPHeihKSbpMQrjlllvIz88nLy+PRYsWdb6v0+mYOHGi24OT/FPxmY4lM6PD+/YQmaNoH6KlFuOsu1wcmXtkJIXzyf5STpY3X7bbCMAwei62L/6Ldd8adAljXLqSnCR1p9uEMG7cOMaNG8esWbOIjY31VEySnysubyJlaN+XzLQd2YQSFosuabyLI3OPc+MIeacvP44AoOgMGCddj/Wz13CezkE/QM5TGvi6TQg/+MEPePbZZ7n33nu73L5unXyyUuqddpuDsppWJoeAPy8AACAASURBVKX1bSDYWVWEVlWIaeadKMrAmHklJNBIQnQQuacauO6Knh1jSJ+N7dCHHa2EYVmylSB5RLcJYeXKlQD87Gc/80gwkv8rqehYMrOv4we2ox+DwYwhbWDdzJAxLILPDvdsHAFAUfWYJi+jfeurOE8dQp8su2gl9+v2X2ZmZiYA06ZNIz4+nmnTptHQ0MCePXsYPXq0RwKU/Mu5AeXhfbjDSGtrxFG4B0PaLJ+axK4n0s89j1Bx+ecRztGPvAIlNAbr/vflHUeSR/Sozf3oo4/yyiuvUFhYyOOPP05ZWRk//elP3R2b5IeKzzQxJMxMaGDvl8y0524FzYFxrG9NYtcTaUk9m9fofIqqwzQxG62mBOepL9wVmiR16lFCOHLkCI899hibNm3ixhtv5Mknn6SsrMzdsUl+qLi8qU/dRUJzYD+2BV1iJmp4vBsic6/Qc+MIJT1PCAD6UVeghETLVoLkET1KCEIIVFVl+/btzJgxAwCLxeLWwCT/09hipbbJ2qcH0hzF+xFtDRgzF7ohMs/ISIogv/Ty8xqdT1H1Z1sJJ3Gelq0Eyb16lBCSkpJYuXIlpaWlTJs2jR/+8Iekp8vH6qXeOTd+0Jcpr+1HP0EJiUY3bJyrw/KYMckR2BwahZdYZ/lS9Gkzz7YS1spWguRWPUoITz75JEuXLuWNN97AYDAwZcoUfv3rX7s7NsnPFJc3oSoKybG9WzLTWVeKs+IExjHzB8ytpl1JT+pYZ/l4L7uNFFWPceJStOpinKcPuyk6SephQggMDGTKlCk0NTVx9OhRsrKy5NQVUq8VnWkiMSYIo0HXq+Psx7aATo8+/Uo3ReYZgWYDw+NCep0QAAyjZqEER2E9IFsJkvv0aLbTp556ijfffJOoqKjO9xRF4ZNPPnFbYJJ/0bSOJTOvGBvXq+OEvR17/nb0qdNQzb1rWfiijOQINu45jdXmxGTseWJUdHqME67D+vnrco4jyW16lBDWr1/Pxo0b5fQVUp+V17bSbnP2+g4je/5OsLdjHHOVmyLzrDHJkazfdYoTpQ2MS426/AHnMaRfie3gOmwH/yMTguQWPeoyio+Pl8lA6pfCMx0DyiMSLj+XzzlCCOzHN6NGJaHGjHBXaB41MjEMnar0qdtI0Rkwjl+CszwPx5lcN0QnDXY9SghXXHEFv/3tb9m/fz9Hjx7t/JGkniosayTIrCc2oudPGGtVhWi1pzvWCPCTuXxMBh0jEsL6lBAADBlzUQJCsR1Y6+LIJKmHXUZr1nQs/L1hw4bO93oyhrBu3TpefPFF7HY799xzD3feeecF2z/++GOef/55hBAkJiby5JNPEhbW8ytIaeAoOtPEiISwXv1itx3b3DFv0cgZbozM88YkR7D282JaLHaCAwy9OlbRGzGOvxbrrndxVuSjixvlpiilwahHCWHz5s29LriyspJnnnmGNWvWYDQauf3225k+fTojR44EoKWlhccee4x///vfxMbG8uyzz/L888/zyCOP9Louybe1tTs4U9PK1NExPT5GtLfgKNqDIX0OisHsxug8LyM5gvc/LybvVAOT03s/66th9FXYDv0X68H/EHjtD90QoTRY9ajLqLW1lccff5yvf/3rNDQ08Oijj9La2trtMTt27GDGjBmEh4cTGBjIokWLLmhh2O12Hnvssc6xifT0dMrLy/txKpKvKq5oQgAjhva89WfP3w5OB4bR890XmJekDg3FZNBxvKTu8jt3QTGYMIxbhPP0YZxV8vZvyXV6lBCeeOIJQkJCqK2txWQy0dLSwqOPPtrtMVVVVURHf3n1ExMTQ2VlZefriIgIFizomKSsvb2dl19+ufO15F8KyxpRgJQeTlnRMZi8FTVmBLqoYe4Nzgv0OpVRw8I4drJv4wgAxrFXgykI20G5JonkOj3qMjp+/DhPPvkk27ZtIyAggKeffpqlS5d2e0xXD8901X/c3NzMfffdR0ZGBjfeeGMPw+4QFRXcq/3PFx098O9p7y1vnXNpTRuJsSEkD4vo0f6WU0dpaSgneun3CelnzL76PU/PHMpf/nMEodMRExnYhxJCMExbSv1n7xKq1WKKHd65xVfP2Z3kObtGjxKCql7YkHA6nRe991WxsbHs27ev83VVVRUxMRf2IVdVVfGtb32LGTNm8PDDD/c05k61tS1oWu+f2oyODqG6uufz0vsDb52zEILjxbVMTIvucf2Wnf8FYwCWmCza+xGzL3/Pw2OCAPjswGnmjB/apzJE6hzYtZbKT/5BwML7Ad8+Z3eR59xzqqp0eyHdoy6jqVOn8tRTT9He3s5nn33G/fffz7Rp07o9ZubMmezcuZO6ujosFgsbN25kzpw5ndudTiff/e53ufbaa/npT3/qN7cVSheqqrfQ2u5gZA+fP9Dam3EU7euYqkFvcnN03jM0KpCIEBNHivs2jgCgmIIwjl2Ao3g/zjo5Hb3Ufz1KCD/60Y8IDAwkJCSEP/zhD2RkZLBq1apuj4mNjeWhhx5ixYoV3HDDDSxdupSsrCxWrlzJ4cOH2bx5M8eOHeOjjz5i2bJlLFu2TC6644cKz3TM7NnTJ5QdeZ+D5p+DyedTFIWxwyM5frKuT63ccwxZi0BvlGMJkkt022V09913X3DlbjZ33P538OBB7r33Xl5//fVuC8/OziY7O/uC91555RUAxo0bR26ufNrS3xWWNWE26hgaFXTZfYUQ2HK3oosdhS4ywQPRedfYlEg+P1zOyYrmPq8xrZpDMI69GlvOerTJy2AQ9qVLrtNtC+Guu+7izjvvJCYmhsDAQFasWME3vvENIiIiSE5O9lSM0gBWWNZISnwoqnr5LkFneS6isRLD6HnuD8wHjBkegQIcLa7tVzmGrMWgGrAe/MA1gUmDVrcthEWLFgHwl7/8hXfeeadzIHnevHksX77c/dFJA1pbu4PT1S1kzxzeo/3tx7eCKQh96lS3xuUrQgKNJMWFcLS4juxZKX0uRw0IxTBmPvYjm7DXlQN9v/tOGtx6NIZQX1+P1WrtfN3a2kpjY+9WfZIGn6IzjQgBo4aFX3ZfzdKEo3g/hlEzUfRGD0TnGzJTIik804TF6uhXOcbx14Kqo2HHGhdFJg1GPbrtdOnSpdx2220sXLgQIQQbNmzgtttuc3ds0gCXX9qIotCjNZQdJ7afHUye5/7AfEhmSiQf7iwht6SeiWm9n8biHDUwHMPoeTQf3kzQ6GtRQ/teljR49aiF8IMf/IAf/OAHNDU10dzczKpVq7j33nvdHZs0wOWXNpAUE0KAqfvrjs7B5Lg0dBH+P5h8vhEJYZgMOo6c7Pvtp+cYxy9BUVRsh+QdR1Lf9KiFALBgwQI5tYTUYw6nRtGZph49dNU5mDzxeg9E5lv0OpWMpHCOFNUihOjX8zhqUAQhExfStP8jjBOyZStB6rWBu2K55NNOV7Vgc2iMTLz8A2n2Y1sG1WDyV2WNHEJ1QzsVdW39Liv8ihtAVWQrQeoTmRAkt8g/3QDAqMTuB5Q1SxOOk/vPPpk8eAaTz5d1dinNLwr6d/spgD40CkPGPOx529GaqvtdnjS4yIQguUV+aSNDwsxEhHQ//YQ97zPQnINuMPl8UWFmEqODyCmscUl5xgnXyVaC1CcyIUguJ4Qgv6zxsq0DITTsx7eii09HF9G3Cd78xfiRQ8gvbaSt3d7vstSgiPNaCVUuiE4aLGRCkFyuqsFCU6uNUZcZP3CWHUM0V/v9vEU9kTUiCqcm+jXZ3fk6Wgkq1gP/cUl50uAgE4LkcvmnOx5avFxCsB/bgmIOQZ8y2RNh+bQRQ8MIMuvJKez/OAKcbSWMuQpH/na0BrkSodQzMiFILldQ1kCQWU/8kEtPaKe11uMoOYg+bTaKrncLzfsjVVUYlxrF4aLafs1+ej7jhOtAZ8S6f61LypP8n0wIksvlnW5kZEIYajf31NvzPgWhYRzEg8lflTUyiuY2O8XlTS4pTw0IxZi5EEfhbpx1pS4pU/JvMiFILlXfbKWyro30pEsvlyk0J/bj29AljEUNi/VgdL4tMyUKRYEvXNRtBGDMWgwGM7Z977msTMl/yYQguVTuqY6F40cnXzohOE4dQrTWYRhzlafCGhCCAwyMTAgjp8A1t58CKOZgjFmLcJzcj7O62GXlSv5JJgTJpXJL6gk06RkWc+kpmO1HN6MERaJPnuDByAaGiaOiOVXVQnWDxWVlGsctQjEFY92z2mVlSv5JJgTJpXJP1ZOeFH7JBXG0hgqcZUcxjJ6Houo8HJ3vm5TeMf/QgROue8pYMQZgnJiNs+wojrJjLitX8j8yIUguU9vYTnVDOxndjB/Yjm8BRYchY44HIxs4YsIDSIoJZr8LEwKAYcx8lKBIrHv+hRCuuYtJ8j8yIUguc278IOMS4wfCYcWe9xn6lMmogZdfNGewmpQeTWFpIw0t1svv3EOK3ohpyo1o1cU4ive5rFzJv8iEILlMbkk9wQEGEqK7fv7AUbAbbG0Yxl7t4cgGlslp0QjgoItbCfpRs1AjhmLb+2+E5nRp2ZJ/kAlBcgkhxJfjB108fyCEwHb0E9SIBHRxaV6IcOAYOiSI2MhAl3cbKaqKcerNaI0VHZMKStJXyIQguUR1Yzu1TdZLjh84K/PRakswjF3Qr0VgBgNFUZiSHk1uSQMtlv5Pdnc+ffIkdHFp2PatQdhcdyeT5B9kQpBcIrek+/ED+5FNYArCMGqmJ8MasCalRaMJwRcufCYBOpKNacZyhKUJW856l5YtDXwyIUgukXuqntAgI0OjAi/aprXU4ijejyF9Doqh+/URpA7D40KIDDWxP8/1i9zoYkagHzEd2xcb0FrrXV6+NHDJhCD1mxCC4yX1ZCSFd9kdZD+2GRAY5WByjymKwuS0GI4U17pkjYSvMk27BYSGde8al5ctDVwyIUj9drqqhcYWG5kpURdtEw4btuNb0SdPQg0Z4oXoBq4ZY2NxOAX73NBKUEOiMWQuxHHic5w1JS4vXxqYZEKQ+u1wUcdkbJmpkRdtsxfsBGsrhsyFng5rwBseF0JMRAC7j1W6pXzTxKUo5mCsO96SD6tJgEwIkgscLqojKSaY8OALxweEENgPb0KNHIYuPt1L0Q1ciqIwY0wsuSX11De77iG1zvJNQRin3oyz4gSOwt0uL18aeGRCkPqlrd1BYVkjmakXdxc5S4+g1ZdizFosbzXtoxlj4xDgtlaCIX0O6pBkrLv/ibC7PulIA4tMCFK/HC+pw6kJxnXRXWTLWY8SGI5+xHQvROYf4iIDGR4X4raEoKgqppl3IVrrsB36wC11SAOHTAhSvxwuqiPApGNEwoXrJztrSnCWHcOQeQ2KTu+l6PzDjLFxlFQ2U17b6pby9XGj0I+cgS1nPVpTlVvqkAYGmRCkPhNCcLioljHJkeh1F/5TsuVsAIMZ45h53gnOj0wbHYOiwM6j7mklAJimLwdVT/v2N+UA8yAmE4LUZ2dqWqlvtl50d5HWUoujcDeGjLkoxosfVJN6JzzYxOjkCHYdrUBz0y9rNSgC05QbcZ7OkbOhDmIyIUh9drioDoBxXxlQth3ZBIBx3DUej8lfzcqMp6axnbwS9z1ZbBi7ADUqqeM2VDnP0aDk1oSwbt06lixZwsKFC3nrrbcuud9PfvIT1qyRT0wONIeLakkYEkRkqLnzPdHegv34VvQjpqEGX3znkdQ3k9OjCTLr2fbFGbfVoag6zFfeg2hrxLpP/n8cjNyWECorK3nmmWd4++23Wbt2Le+++y4FBQUX7fPd736XDRs2uCsMyU1a2+2cON1A1oguWgf2dowTsr0UmX8yGnRcMTaOAyeqaW6zua0eXUwqhjHzsR/9GGd1sdvqkXyT2xLCjh07mDFjBuHh4QQGBrJo0aKLfvGvW7eOq6++mmuvvdZdYUhu8kVBDU5NdK4BDCBsbdiObEI/fDK6yAQvRuef5kwYisMp2HGkwq31mKbejBIQRvu2vyCcDrfWJfkWtyWEqqoqoqO//GURExNDZeWFd0nce++93Hrrre4KQXKj/XnVRISYSIkP7XzPdnQz2NowTpKtA3dIjA5mREIo2w6dceudQIopCPOV96DVlWI7uM5t9Ui+x203iHf1D9bVT6tGRQX3+djo6BAXRjIwuOqcLVYHR4vruGZGMrExHQlBs7Vz6uhGAkZMJG70OJfU4wr+9j0vnZ3Ks+8eorrFztgung4HF51z9JVUlR2g5dAHRE+agyl2eP/LdCN/+557wh3n7LaEEBsby759X96+VlVVRUxMjEvrqK1tQdN6f6UUHR1CdXWzS2Pxda485725VdgcGmOTwjvLtOV8hNbWhDJ2ic98tv74PWckhBFg0rF2awExIcaLtrv0nCfdhlJ4iPL3nyPwhp+hqL75gKE/fs+X09dzVlWl2wtpt3UZzZw5k507d1JXV4fFYmHjxo3MmTPHXdVJHrQ/r4qQQAOjEsMBEA4rti/+i27oaHRxo7wcnX8zGXXMGBvH3twqmtw4uAygmIMxzbobraYE28EP3VqX5BvclhBiY2N56KGHWLFiBTfccANLly4lKyuLlStXcvjwYXdVK7mZ3eHki8JaJo6KRlU7ugBtRz5BWBoxTrnRy9ENDgsmJ+Jwamw5UOb2ugypU9GPvALbgbU4Kwsuf4A0oLm1DZidnU129oUDjK+88spF+/3mN79xZxiSCx0trsdqczLl7N1FwtaG7YsP0Q3LQh+X5uXoBof4qCCyRkSx+UAp105PwmjQubU+8+y7aa04gWXznwm6+XEUY4Bb65O8Rz6pLPXK/hNVBJr0ZCRHAGfnLLK2Ypp6s5cjG1wWTR1Gc5udXW6aBfV8ijEQ8/xvI1pqaN/xttvrk7xHJgSpx+wOJ4fyaxg/cgh6nYpmacKW8xH61GnohiR7O7xBJSM5gqSYYD7ac8pt8xudTx+fjnH8dThOfIa9cI/b65O8QyYEqccOFdTS2u5gZmYcALZDH4LTjmnKTV6ObPBRFIVF05Ior23jyNklTN3NOOUG1JgRtH/6V7QG9z4cJ3mHTAhSj32eU05ESMfMm1pTFfajn2BIm40aHuft0AalqaNjiAgx8dGe0x6pT1H1BCy4D1Qdlo//iHDIFdb8jUwIUo/UN1s5UlzLrHFxqKqCddc7oNNjnCpbB96i16ksmJLI8ZJ6CsoaPVKnGhxFwFXfQasro/3zNz1Sp+Q5MiFIPbLzaAVCdEzD7Cg7huPkAYwTlqIGhns7tEFt/sQEQgINvPdpkcfq1A/LwjhxKY4Tn2E7vtVj9UruJxOCdFlCCLYfLmdUYhgx4SasO99GCYmW6x34ALNRz3UzkjleUk+uG9dK+Crj5BvRJWZi3f4GjvI8j9UruZdMCNJlFZ1pory2jVnj4rHnbkOrK8U0/TYU/cVTJ0ieN29iAuHBRt77rMhjy18qqkrA1d9DCYmmfdMf0ZqqPVKv5F4yIUiXtf1wOUaDypThAdj2rkEXn44+ZYq3w5LOMhp0LJ05nPzSRg7mee4Xs2IKInDR/yA0J5aPnpWrrPkBmRCkblmsDnYfr2RyWgzKgX8h7BZMs+52+cy1Uv9cmTWUqFAzb2447rFWAoAaHkfAgvvQGs5g+fhPcv2EAU4mBKlbn+eUY7E6WZzUjCN/B8bxS9BFJno7LOkrDHqV62cPJ/90A7s98PTy+fSJmZiu/DrO0iMdi+oIzaP1S64jE4J0SU5NY9O+02QkBBJxbDVqWBzGiXLxG181a1w8o4aF8+6WAixWz16pGzPmYpxyE46CnVh3/9OjdUuuIxOCdEkHTtRQ09jO8qjjiOZqTHO+IQeSfZiqKHz3piwaW2ys237S4/UbJ2ZjGHM19pwNWA/J6bIHIpkQpC4JIdiw+xSTwxuILPsUQ8Zc9PHp3g5Luoy0pAiuzIpn077TnKlp9WjdiqJgmnkn+hHTse35F7ac9R6tX+o/mRCkLuWXNlJRUcNy8zaUkGhMM273dkhSD908bwQmg463Pz7h0QFm6Lgd1Tz/2+hTp2Ld9S62nI88Wr/UPzIhSF36aHcJXwvZjcneQsBV35Vz4A8goYFGbpyTyrGT9Xx+uNzj9SuqDvNV30GfMgXrrn/IlsIAIhOCdJFTlc0YTu1hvL4Y45Qb0MWkejskqZfmT0ogIymctz/Op6rB888HKKoe89Xf7WwpWPf8y+OtFan3ZEKQLvLxx7u5JWg3xKZhHH+dt8OR+kBVFL513RhUBV794Bia5vlfxoqqx3zV9zCMnoft0IdYP/0rQnN6PA6p52RCkC5w7MRp5jWtRTGYCVrwPRRV/hMZqKLCzNy1MJ2C0kbW7y7xSgyKqmKa/XWMk5Zhz/sMy8bn5BPNPkz+b5c6OR0OnFv/TITaSvDiB1GDIrwdktRPM8bGMiUjhvc/K6ag1DNTZH+VoiiYptyIafYKnKcP07b2l2hNVV6JReqeTAhSp5L1f2M4pVSl3YRpaJq3w5FcQFEUvr44nahQM3987zB1Te1ei8U45ioClvwIra2R1vd+gaPsmNdikbomE4IEQOv+D4ku385BXRZp8+S4gT8JMht44JYsbHYnz//7MFa79/rx9QljCLrx56iBYVj++xTW/e8jNDnVha+QCUHCdnwr2v5/cdCazNBF35AT1/mhhCFBfOf6sZyqbOavH3p2AryvUkNjCLzhUfQjZ2Lb/z6WD3+L1uq5tRykS5MJYZCzF+yi/bO/c8yWQFXmHYxMlOMG/mr8yCHcMm8Ee3OreHtTvleTgmIwEzB/JeZ59+KsLqJ19SPY83fIW1O9TO/tACTvsefvoH3rq5RosXxkvpb/vXKUt0OS3Gzx9CQaW21s3HsanU5h+VUjvdoiNKTNRhczAsu2v9C+5WX0RXsxXfl1uTSrl8gWwiBly/mI9i0vU6FP4M9N87knezwGvfzn4O8UpSMJLJicyMa9p1m9tdDrV+VqeDyB2Q9jmvE1HKVHaH13FbacDQhNrq3gabKFMMgIoWHb+29shz6kJmwMTxVPZNncUSTFhng7NMlDFEXhawtG4dQE63eforXdzl3XpKPXee+CQFFVjFmL0CdPoH3HW1h3vYM991NMV3wNXWKmHNfyEJkQBhFhbcWy5WWcp76gLm46TxwbxdTRcVw7I9nboUkepigKd16TRqBZz4c7S6huaOe+GzMJMhu8GpcaFkvA4odwnjpE+463saz/Hbr4DEzTbkEXO9KrsQ0Gso9gkLCWF9K65uc4S4/QMPomfpWXTtqwyLPTG8irr8FIVRRunjuCb103mhOnG/jV6/sprW7xdlgoioI+eSJBt/0a08w70RrO0Lb2CdrW/x5HeZ7Xu7j8mSIG8KdbW9vSpzlaoqNDqK5udkNEvkdoDmxfrMd2YC2KOZSa8St4elMT4SEm/veuSV6/InSnwfQ9n9PXc847Vc+L7x+hzerklrmpLJg6zGcuFIS9HduRTdgPb0S0N6PGjsSYtRh98kQUVSe/515QVYWoqOBLbpcJwY85q4tp3/ZXtLrTBI2+gkPh1/DShhKiQk388PYJDAnz7ymtB8v3fL7+nHNTq43X1udyqKCG0ckRrFicTmxEoIsj7DvhsGLP+wxbzkeI5mqUwHAMGXOJm7WEeqvJ2+F5lEwIXZAJoWtacw3W/e/jyN+OEhCGadbd7GpN4G8fHGNkYhgP3DSOkED/XwrT37/nrvT3nIUQfJZTzj8+ycfh0Fg4ZRhLZw4n0Ow7w41Cc+I8lYPt+Gacp48AoBuagWHUTPQpk1GMvpPE3EUmhC7IhHAhraUWW84G7Me2gAKGMVfTMnIRr28+xZHiOqZkxLBy6WgMep23Q/UIf/2eu+Oqc25ssbLm0yI+zyknKMDA4ulJzJswlEAf62LUmqowlO2j8YttiKZKUPXoEkajT56EPnmC307QKBNCF2RC6Liic1acwH5kE46TBwCBIf1KdBOu59P8dlZvKwQBX79uDNPSh/hMv7An+NP33FOuPueSimZWby3g6Ml6zEYdcycMZd7EBJ/qSoqODqGqqgmtugh70V4cJw8gzs6mqkYmoksYiz5xLLrYUX6z8p9MCF0YrAlBCIFWV4qjcDf2wt2I5mowBWHMmIuSPpftxQ7W7yqhprGdsSmRfH1ROqNHxQzoc+6Lgf4994W7zrmkopkNe06x93gVmhCMTAxjZmYck9KiCfVy9+NXz1kIgVZfhuPUFzjLjuIsPwGaAxQFNXIYuthR6GJSUIekoIbHD8g1PwZkQli3bh0vvvgidrude+65hzvvvPOC7cePH+eRRx6hpaWFKVOm8Itf/AK9vud9lYMpIWit9Tgr8nGWHcFRehTRUguKii5hDPrUaZwOzGBvfiO7j1XS2GojJT6UpVckM2HUEBRFGZDn3F/ynF2vvtnKzqMVbD9cTnltGwqQOjSUrBFRjB4eSXJsiMefeL/cOQuHFWdFAc7K/I7/Q1WFYD87DbjeiBqRgBqRiC4yETU8DjUsDiVkCIrqu12rAy4hVFZW8rWvfY01a9ZgNBq5/fbb+f3vf8/IkV8+XLJ06VKeeOIJJkyYwMMPP0xmZiZ33HFHj+vwx4QgHDa0pmq0hjK0+jNotadxVhchzs0GaQxAN3QMlsg08pVUjlU6OH6yntqmdnSqwrjUKBZOSSQjOeKCpzt9+ZzdRZ6z+wghOFXZwhcFNXxRWENxeUedep1KSnwIKfGhJEYHMywmmLioQEwG9/1y7e05C01Da6xAqy7GWVOCVl+KVleKsDR9uZOqQwmOQg0ZghoyBCUoCjUoAiU4EiUgDCUwDMUcjKJ4p3XhroTgtlsHduzYwYwZMwgP75ikatGiRWzYsIH7778fgLKyMtrb25kwYQIAN910E88991yvEoIvE5oGThvCYQOHDWFvB3s7wmZBWFs7ftpbEJZGRFsjWls9ornmgn+UAgV7QBRNAcOoCZrKKWc0RxtDOXPYgs2uAacJMusZlRjOstkp7Fe0tgAACSJJREFUTEob4nODfpJ/UhSF5LgQkuNCuH52Ck2tNvJLG8kvbaCgrJEtB8uwO75c5yAixERsRABRYWYiQsxEhpgICzYSEmgkNNBAcIABs0nvkTEuRVXRRQxFFzEUQ9qszvc1SxNaYyWisaIjYTRVozXX4Cg5dGGy+LKgjqRgDkYxh3Tc3WQKQjEFoRgDOn4MAWA0o+jNYDCh6E2gN6LojR1/6gygN4Ci84npOdyWEKqqqoiOju58HRMTQ05OziW3R0dHU1lZ6a5wOjXV1XL8779EtVvo+JV7vnOvz3+/4+8KAgUBQqAiUNBQ+PLvqtBQ0dDhPPvTs0U/WoWZFmGmUQugxhFDnZZKnRZMhTOcKmco9rNfkU5VCA82ER9lYm5SJEOHBDIyIYz4IUGDaqBY8k2hQUYmp0czOb3j/7SmCSrr2zhd1UJFXRtV9RYq69s4drKehhYrXfVLKIDZpCfQpMNk1GM26jAZdBj0Kka9ikGvw6BX0OlUDDoVnaqg0ymEhgTQbrGhUxUURUFVQFEV1HN/VxSUc3+erUih4z2gY9u5//FKCAohEDQKzruQVjUnelsTemsjBnszensLelszOntrx09rK7rGRnQOCzqHBVWz9erzEygIVYdQ9AhVj1B0Z1/rQFE7Xitqx99VHYZFd0N4aq+/p8txW0Loqifq/Ax4ue090V3T51L0tFN69kOF7lICcF56EIpy3ja187WGriM1KCoaKpqi70gLig6nosepGnAqBhw6E07F2PGnPhCnPgCnIRC9wYBBp2Iw6DAbdSQYdIw06wkOMBAcaCQ00EhUmJmwYBOq2r9f/NHRg28CO3nO3hMbG0pWxsXvO50a9c1W6pvbaWyx0dhipbnNTovFRmubnTarA4vVQbvVQbvNSavVQV2zE5tDw+FwYndoOJwaDqfA4dS6TC7upQPCzv50TUXDpDgwKzZMigMTdkyKHaPixKA4MOLAoDg7fnCiV5zo0dArHReTekVDjxMVgU45d7EpUBQ7KnbKj1ey8PrxLj8ztyWE2NhY9v3/9u4mJKo1AOP432val4R0r0q1KPrSRUjQImtRCEnZlE7YQluMFZQJZcmlnCwKKtNKcNG0iAiEMIihZiT6oBBalG6chUQi1SL6sBxL6lKaH3fOXTgIca/czDnzMs7zW82Zg2eed6GP5z3vOdPePrYdDAZJT0//Yf/Hjx/Htnt7e3/Y/zN+7RrCDNb/WRdTc8sjg8N8Ghye1DE0nx4fYmnMqTOmkTpjGvwxuSWsc39PIRj8i79DFqGQNbrKyBo9S7EIrzoK/52wrNFti/A/feHXozsZez3epdXxymdSnTTBRktISCA7KyO2riGsXbuWixcv0tfXx8yZM3nw4AGnT58e279gwQKmT59OIBBg1apV+P1+1q1bZ1ccEZmiEn9LYFrib8TJ/ZbAxGdTfpZtl8gzMjKorKzE5XLhdDrZsmUL2dnZ7Nmzh6dPnwJQX19PbW0t+fn5DAwM4HK57IojIiL/QzemxQmNOT5ozPHBrmWnsXeLnoiI2EKFICIigApBRETCVAgiIgLYuOw0GiZzo9Zkb/KKRRpzfNCY48OvjPn/fiamVxmJiEjkaMpIREQAFYKIiISpEEREBFAhiIhImApBREQAFYKIiISpEEREBFAhiIhImApBRESAOC6EYDDI3r17cTqdFBcX8/btW9ORoqazs5MVK1aYjhEVgUCAoqIiCgsLKS0t5d27d6Yj2eb27dts3ryZvLw8mpqaTMexncfjweFw4HA4OH/+vOk4UXXu3DncbnfEjxu3hXDkyBFyc3Px+/0UFhZSX19vOlJUDAwMcOrUKYaHJ/cdzbHi8OHD1NTU0NzczNatWzlz5ozpSLbo6emhoaGB69ev09zczI0bN3j58qXpWLZpbW3l8ePH+Hw+/H4/z5494+HDh6ZjRUVbWxs+n8+WY8dlIfT19dHV1UVxcTEARUVFHDp0yHCq6Kirq2Pnzp2mY0TF0NAQBw8eJCsrC4DMzEzev39vOJU9WltbycnJITU1lVmzZrFx40bu379vOpZt0tLScLvdJCcnk5SUxJIlS+ju7jYdy3afP3+moaGBffv22XL8uCyEN2/eMH/+fM6ePUtBQQEVFRUkJSWZjmW7lpYWvn//zqZNm0xHiYrk5GQKCwsBCIVCeDweNmzYYDiVPYLBIGlpaWPb6enp9PT0GExkr2XLlrFy5UoAXr16xd27d1m/fr3hVPY7ceIElZWVzJkzx5bjx/Tjr3/GvXv3qK2t/eG9hQsX0tnZyYEDBzh27Bherxe32821a9cMpYys/xrz4sWL+fr1K42NjWZC2Wy8MTc2NjI0NITb7WZkZISysjJDCe31Xw8tTkiY+o+EfvHiBWVlZVRVVbFo0SLTcWzl9XqZN28ea9as4datW7Z8Rlw+/vr169ds27aNQCAAjM6r5+Tk0NHRYTiZfbxeL5cvX2b27NkAdHV1kZWVRVNTEykp43/pdqz79u0b5eXlpKamUl9fT3JysulItvD5fLS3t1NTUwPApUuXsCyL/fv3G05mn0AgQEVFBdXV1TgcDtNxbLdr1y56e3tJTEzky5cv9Pf343Q6qa6ujtyHWHEqPz/fevTokWVZlnXnzh2rpKTEcKLoWr58uekIUVFeXm4dP37cCoVCpqPY6sOHD1Zubq716dMnq7+/3yooKLA6OjpMx7JNd3e3tXr1aqu1tdV0FCNu3rxpVVVVRfy4U37KaDwej4eTJ09y4cIFUlJSqKurMx1JIqyzs5OWlhaWLl2K0+kERufWr1y5YjhZ5GVkZFBZWYnL5WJ4eJjt27eTnZ1tOpZtrl69yuDg4A+/t8XFxZSUlBhMFfvicspIRET+LS5XGYmIyL+pEEREBFAhiIhImApBREQAFYKIiISpEEQibPfu3fT19ZmOITJhKgSRCHvy5InpCCK/RIUgEkFHjx4FoLS0dMo+WVWmLt2YJhJhmZmZtLW1MXfuXNNRRCZEZwgiIgKoEEREJEyFIBJhiYmJjIyMmI4hMmEqBJEIy8vLY8eOHTx//tx0FJEJ0UVlEREBdIYgIiJhKgQREQFUCCIiEqZCEBERQIUgIiJhKgQREQFUCCIiEqZCEBERAP4BGjTrffqoPsgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(sol.t, f2, label=\"β=2\")\n",
    "plt.plot(sol.t, f1, label=\"β=1\")\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"density\")\n",
    "plt.title(\"Tracy Widom for β=1,2\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tracy and Widom derived the distribution in their paper and proved its relationship to the fluctuations of the largest eigenvalue of Gaussian ensembles. Let the CDF for the $N \\times N$ Gaussian ensemble with Dyson index $\\beta$ be denoted as $F_{N, \\beta}$. Also define $F_\\beta$ and $f_\\beta$ to be the CDF and density of the Tracy Widom distribution for Dyson index $\\beta$ respectively. Tracy and Widom proved the following theorem \n",
    "\n",
    "$$ F_\\beta(s) = \\lim_{N \\rightarrow \\infty} F_{N, \\beta} \\left(2 \\sigma \\sqrt{N} + \\frac{\\sigma s}{N^{1/6}} \\right)$$\n",
    "\n",
    "We can take the derivative to find the density of the distribution\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{aligned}\n",
    "f_\\beta(s) &= \\lim_{N \\rightarrow \\infty} f_{N,\\beta} \\left(2 \\sigma \\sqrt{N} + \\frac{\\sigma s}{N^{1/6}} \\right) \\frac{\\sigma}{N^{1/6}} \\\\\\\\\n",
    "&= \\lim_{N \\rightarrow \\infty} f_{N, \\beta} \\left(\\frac{s+2\\sqrt{N}N^{1/6}}{\\frac{N^{1/6}}{\\sigma}} \\right) \\frac{1}{\\frac{N^{1/6}}{\\sigma}}\n",
    "\\end{aligned}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "\n",
    "Recall the formula for the density of linear transformations $f_{aX+b}(s) = f_{X}(\\frac{s-b}{a}) \\frac{1}{|a|}$. This implies that if $\\lambda_{\\max}$ is the greatest eigenvalue and $\\sigma$ is the standard deviation of it's off diagonals, then \n",
    "\n",
    "$$\\lim_{N \\rightarrow \\infty} \\frac{N^{1/6}}{\\sigma} \\lambda_{\\max} - 2 \\sqrt{N} N^{1/6} \\sim \\text{Tracy Widom}_{\\beta}$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GOE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall from the first chapter the off diagonal terms in a GOE have standard deviation $\\frac{1}{\\sqrt{2}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "trials = 10000\n",
    "GOE_samples = []\n",
    "N = 400\n",
    "σ = 1/(np.sqrt(2))\n",
    "for _ in range(trials):\n",
    "    A = RM.Generate_GOE(N)\n",
    "    λs = np.linalg.eigvals(A)\n",
    "    λ_max = max(λs)\n",
    "    GOE_samples.append(N**(1/6) * (λ_max/σ - 2*np.sqrt(N)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true,
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEXCAYAAACzhgONAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deVxU9R7/8deswAiI4IC4r4mJuGRKuGVpqIma6c2yvN6bVvZLi5sWt80ly65Zdq+Vt7Rbv59LqakQVzOXVpNMXHHBQnMDBQRR9tnO7w+uk4QiKMNhhs/z8eChc86cOe8vMHzmfL/nfI9GURQFIYQQ4hq0agcQQghRt0mhEEIIUSkpFEIIISolhUIIIUSlpFAIIYSolBQKIYQQldKrHUCI65k7dy67du0C4NixYzRr1gxvb28AVq1a5fy/qw0dOpRnn32WQYMGAbB9+3YeffRRZs+ezbhx4wA4cOAAU6ZMYfv27YwaNYply5bh7+9f7nU++ugjfv31V9544w2X5ExPT2fRokXs3bsXvV6PzWZjyJAhPPXUUxgMBgBsNhtLliwhMTERjUYDQK9evXj66acJCAgA4JFHHiE9PR0/P79yrz9lyhSGDBnikuyibpJCIeq8l156yfn/u+66iwULFtClS5daz9G/f3927tzpLBTffPMNAwcO5Ouvv3YWip9++on+/fuj0WhISEio9YyZmZk88MADTJs2jXnz5qHRaCgsLCQuLo5//OMfzu/l3/72NwwGAytXriQgIACr1conn3zCAw88wNq1a/H19QXgueeek6IgpOtJuL/w8HCefvppoqOjSUlJ4fPPP2fs2LGMGjWKgQMHsnLlSudzP/jgA4YMGcLw4cP5P//n/5Cfn89f/vIXVq1a5XzO4sWLef311yvsp3///vz888/Ox9988w3PP/88e/fupaioCICkpCQGDBgAQMeOHcnNzcVqtTJz5kwGDx7MuHHj2LNnj/M1zp07xxNPPEFMTAzDhw9n6dKlAJw5c4a7776bF198kXvvvZfo6Gi2bdvGY489xqBBg3jmmWdwOBwVMn744Yfcc889/OlPf3IeKTRo0ICXX36ZVq1aAbB371727dvHP/7xD+fRg8FgYPLkybRt25bPPvvsxn4QwmNJoRBuz2q1MnDgQL766ivatm3LmjVr+PDDD4mPj2fhwoW8+eabAGzbto1169axatUq/vvf/9K8eXOWL1/O+PHjWbNmDQAOh4M1a9Y4jxCudPvtt3Pq1Cny8vI4evQoDRs2pE2bNkRERPDjjz9isVg4cOAAffv2LbfdypUrOXHiBBs2bOA///kPZ8+eda6bPn06vXv3JjExkU8//ZQvvviCDRs2AGXF4q677mLDhg1ERkby2muv8fbbb7NhwwaSk5PZt29fhYzJyckV9g8QHBzMI488ApQVim7duqHXV+xQiIqKYvfu3c7H8+fPZ+TIkeW+Lly4cN2fifAs0vUkPELPnj2Bsk/P//73v/nuu+84ceIEqamp5T7tDxkyhIYNGwLw97//HQC73c7cuXNJTU0lMzOT5s2b07Zt2wr7MBqN9O7dm+TkZNLS0rjzzjsBGDhwINu3b8ff35/OnTs7u20uS0pKYvjw4RiNRoxGIzExMRw9epSioiL27NnDf/7zHwD8/PwYPXo033//PV27dsVgMHDXXXcB0LJlS7p37+587eDgYC5evFgho6IoziMJgKVLl5KYmAjA+fPnnUXIZrNd9ftosVjKbS9dTwLkiEJ4CJPJBJR15YwaNYr09HRuu+02nnnmGedzdDpduT+Cly5d4syZM+h0OsaNG8fnn3/O2rVrnUcTL774ovNT9KeffgpAv3792LVrF99++y0DBw4EYMCAAezZs4effvrJWTwqo9PpgLKjlz9OteZwOJx/xA0GQ7m8lweiK9O9e/dy3WOTJk0iISGBhIQEzp8/j8PhoEePHhw4cIDi4uIK2+/cuZPu3btfdz+ifpFCITzKwYMHCQwM5Mknn6Rfv3588803QNlRQ1RUFFu2bKGgoACARYsW8cknnwAwduxYtm7dyqFDhxg8eDAAr732mvOP7IMPPgiUjVP8+OOPZGRkOAfUW7RoAcDWrVud4xNX6tevH/Hx8ZSWllJaWsrGjRsB8PX1pWvXrqxYsQKA/Px84uPjiYqKuuH2T5kyhS+//JL4+Hjsdruz7Zf3qdVq6datG7169SIuLs55VGK321m8eDEnTpy4arebqN+k60l4lD59+vD5558zZMgQfHx8iIiIIDAwkJMnTzJgwADS0tKcf/Tbt2/Pq6++CkBQUBDh4eG0a9eu0k/uLVq0wGaz0adPn3Kf9vv168fmzZtp165dhW3GjRvHqVOnGD58OAEBAc5BZYAFCxYwZ84c1q1bh8ViISYmhtGjR5Oenn5D7W/SpAmrVq3i3Xff5aOPPkKn01FcXEznzp1ZvXq1c/D6zTff5OOPP+bhhx8GyrqcevfuzWeffVbudNj58+ezePHicvsYPHgwTz311A3lE+5JI9OMCwG5ubmMGTOGFStWEBoaqnYcIeoU6XoS9d7q1asZNmwYEyZMkCIhxFXIEYUQQohKyRGFEEKISkmhEEIIUSkpFEIIISolhUIIIUSlPPI6igsXCnE4qj9GHxTkS05OgQsS1R2e3kZpn/vz9DbWxfZptRoaNWpwzfUuLRSJiYksXrwYq9XKxIkTGT9+fLn1W7Zs4V//+hcOh4MuXbowZ84cjEYj8fHxLFiwgKCgIADuvPNOYmNjq7xfh0O5oUJxeVtP5+ltlPa5P09vo7u1z2WFIjMzk4ULF7Ju3TqMRiPjxo2jd+/etG/fHoCioiLmzJnD+vXrady4MbGxsaxfv54HHniAlJQU4uLiGD58uKviCSGEqCKXjVHs2LGDyMhIAgICMJlMREdHs2nTJud6k8nE119/TePGjSkqKiInJ8d5J7CUlBTi4+MZMWIE06dPv+osmUIIIWqHywpFVlYWZrPZ+Tg4OJjMzMxyzzEYDHz33XcMHDiQCxcuOOfRN5vNTJ06lYSEBEJDQ5kzZ46rYgohhLgOl12Z/e9//5vi4mLn2MKaNWtISUm55h/9t99+m/T0dN56661yyy9evMigQYOc90wWQghRu1w2RhESEkJycrLzcVZWFsHBwc7HeXl5HDx40HkUERMTQ2xsLPn5+axdu5aJEycCZTdiudqduCqTk1NwQ4NFZrMf2dn51d7OnXh6G6V97s/T21gX26fVaggK8r32elftOCoqiqSkJHJzcykuLmbz5s3079/fuV5RFGbMmEFGRgYAX375JT169MBkMrF06VL2798PwPLly533BxBCCFH7XHpEERsby4QJE7BarYwZM4aIiAgmT57MtGnT6NKlC6+++iqPP/44Go2G9u3bM3v2bHQ6He+88w6zZs2ipKSE1q1bM3/+fFfFFEIIcR0eOXusdD1dm6e38Y/t8/P3wdvr989DJaU28i9VvAWou/D0nx94fhvrYvuu1/XkkVdmC3GZt5eemGcTnI8T3xpJ3XqLClH3yVxPQgghKiWFQgghRKWkUAghhKiUFAohhBCVkkIhhBCiUnLWkxBX8LTTaYWoCVIohLiCnE4rREXS9SSEEKJSUiiEEEJUSgqFEEKISkmhEEIIUSkZzBb1nuJwYE8/iD0zjbMXTvKc/ylO24M4bjVju5QDGNWOKISqpFCIes2Rd47i75biyEwDjQajuRWXHD50MZwi0iuN04t3Y+x5P4bwwWi0cgAu6icpFKJeUhQF66GtlO5cA3oD3gMeRd/2doKbmpnybAIaFJro8phzWwZFP32K9bdd+Nw5GW3DELWjC1Hr5COSqJcsexMp3bECXbNONBgzF0PHfmgM3s71ChrO2hsR8qe/4z3wMRx5ZylKnIfjUpaKqYVQhxQKUe9Yj/6AJXkd+g5R+EQ/g7ZBo2s+V6PRYOgQhSnmBRS7laINb+IovFCLaYVQnxQKUa8Upe2h5PuP0TXrjHf/v6LRaKq0nS6wGaahz6KU5FO88U0cJXK9tqg/pFCIeqOhppDM9W+jDWyOz+Cn0Oj0+Pn7YDb7Ob/+yGK1O9c16dyVoPuew3Epi5Jvl+KBdxEW4qpkMFt4lMt/2CtSGNtgJzjsZUXC6ANcfW6nKxkNugrrC3qNpTTpU6xHv8cYNsAl7RCiLpFCITzK1f6wA3QznKSL8QyNBkzA6h98U/swhA/GdnIfpUmfom/aCe1Nvp4QdZ1Lu54SExMZNmwYgwcPZsWKFRXWb9myhZiYGO69917i4uKwWCwAZGRkMH78eIYMGcKUKVMoLCx0ZUzh4UyaUsY0+JlTtiB8ug+ttKupKjQaLd53TgI0ZV1QDkfNBhaijnFZocjMzGThwoWsXLmShIQEVq1aRVpamnN9UVERc+bM4eOPP2bDhg2Ulpayfv16AGbPns1DDz3Epk2bCA8P5/3333dVTFEPxPjswaQp5bPCO/DyMhLzbILz60ZpfYPw7vMw9nO/YE39tubCClEHuaxQ7Nixg8jISAICAjCZTERHR7Np0ybnepPJxNdff03jxo0pKioiJycHf39/rFYru3btIjo6GoDRo0eX206I6rDmZhDplcYPpWGk2wNr9LX1HaLQNbkFy+54FIvc3Eh4LpcViqysLMxms/NxcHAwmZmZ5Z5jMBj47rvvGDhwIBcuXKBv375cuHABX19f9Pqy4ROz2VxhOyGq6sL2z7GhZWtxeI2/tkajwStyHErxJSz7NtT46wtRV7hsMPtqpw5e7Zz1AQMGsHPnTt5++21mzZrFc889V6XtKhMU5Fut51/pRvut3Ul9aCNAsPYiBQd/4IeSMPIVnxp73XLfP3NXsn7tR+HBzTTpF4Pev3GN7adK+/dQnt5Gd2ufywpFSEgIycnJzsdZWVkEB/9+dkheXh4HDx6kb9++AMTExBAbG0tgYCAFBQXY7XZ0Oh3Z2dnltquKnJwCHI7qn+NuNvuRne3ZF1J5ehuvfANG+xxAozfwdUnnGt3HH79/SsRIlNQkMjb9X3wGPlaj+/ojT//5gee3sS62T6vVVPoB22VdT1FRUSQlJZGbm0txcTGbN2+mf//+zvWKojBjxgwyMjIA+PLLL+nRowcGg4GePXuyceNGAOLj48ttJ0RVhGjz6GH8Df+eQymowaOJq9H6NcbYJRrbrzuw555x6b6EUIPLCkVISAixsbFMmDCBUaNGMXz4cCIiIpg8eTIpKSk0atSIV199lccff5wRI0Zw4sQJZsyYAcDMmTNZvXo1w4YNIzk5mWeeecZVMYWHGuRzECt6AiJHXv/JNcAYMRT0XjJWITySSy+4i4mJISYmptyyJUuWOP8/aNAgBg0aVGG7Zs2asWzZMldGEx6sgaaEHsYT7CjtQCeTf63sU+Pti+HWgVhTvsLR8z65CE94FJnrSXicO7x+Ra9xsL2kY63u1xgxBDQ6LPs31up+hXA1KRTCoygOO329jnLU2oRMR0Ct7ltrCsDQsR/Wo9tlKnLhUaRQCI9S9MsuGumK+L4kTJX9G7sOBcWBJeUrVfYvhCvIpIDCo1xM/pJcewMOWZvXyv78/H3w9vr9bVTi70Npu95YD3+DV48RaIymWskhhCtJoRAew56bTsnJg2wv7YHiooPlq01j/sfZao1dorGlJWH95UeM4YNdkkOI2iSFQngM69HvQavnp9L2LtvHtaYxv5LO3BptcDssh7Zh6Hw3Go308Ar3Jr/BwiMoDju2tCRMHW6jUPFWOw7G8EEoF89hP3NI7ShC3DQpFMIj2M8cRCm+hF+XO9WOAoC+ze1ofPyxHNqqdhQhbpp0PQm3UmHwuNRG/qVirL/8iMbLF1P77sBZ1fJdOYaRe9s95G1fi+NSllyAJ9yaFArhVq52j+tLpYXYTu7BEDYAjc6gYrryYxgNNTAnSIvl0Da873hQ1VxC3AzpehJuz3p8F9htGDr0UTtKORcVEw1u6YXtlx9R7Fa14whxw6RQCLdn+3UH2oBQtOY2akepwK/b3SilBdhO7lU7ihA3TAqFcGvWC+ewn/sF/S19qn2Dq9rg0yYCTYNArKnfqx1FiBsmhUK4tcIjSQAY2vVWOcnVabQ6DB37YT9zCEf+ebXjCHFDpFAIt5Z/eAdeoe0Jadu2zt5e0tCx7C6O1l+2q5xEiBsjhUK4rSBtPtbM46w+5k/MswnlzoaqS7R+ZnTNbsV69AcUxaF2HCGqTQqFcFtdjScB2GdppXKS6zN07IdSkIM9/bDaUYSoNrmOQritbsaTeIW2Ize3bnY5XUnfugd4NcB6dDv65uHO5de6gFCIukQKhXBLgdoCWulzaNBpGNTh6ZSuvFJb07kPBQe+RbEUozH6AFe/gDBflaRCXJsUCuGWLnc7NQiLBHaqG6YSV16p3VqvJdbfgu3Ebgy39FU5mRBVJ2MUwi11M57ktC0QQ6MmakepshM2M/qAEKy/7lA7ihDV4tJCkZiYyLBhwxg8eDArVqyosH7r1q2MHDmSESNG8OSTT3Lx4kUA4uPj6du3LyNHjmTkyJEsXLjQlTGFm2moKaS1/rxbDGKXp8E3vD/29CM4CnLVDiNElbms6ykzM5OFCxeybt06jEYj48aNo3fv3rRvX3ZTmYKCAmbNmsXatWsJCQnhn//8J4sWLeKll14iJSWFuLg4hg8f7qp4wo2FG88AkGJpoXKS6vPr0p+87Wuwpv2EV7dhascRokpcdkSxY8cOIiMjCQgIwGQyER0dzaZNm5zrrVYrs2bNIiQkBICOHTty9mzZ9NApKSnEx8czYsQIpk+f7jzSEAIg3HCabLsfmY6GakepNkNgU7TB7bClSfeTcB8uKxRZWVmYzWbn4+DgYDIzM52PGzVqxKBBgwAoKSnhww8/dD42m81MnTqVhIQEQkNDmTNnjqtiCjfjKC3mFsM5DlpaAHVvbqeqMHSIwpF7BnvOKbWjCFElLut6UhSlwrKrTdqWn5/Pk08+SVhYGPfddx8A7733nnP9pEmTnAWkqoKCfKuZ9nd1dRqImuTObSxITUKvcZBidb9up8ua9Lqbk0krMWTsgbDOFdZf7+fjzj+/qvL0Nrpb+1xWKEJCQkhOTnY+zsrKIji4/F2+srKyePTRR4mMjOSFF14AygrH2rVrmThxIlBWcPT66sXMySnA4ahYqK7HbPYjO9uzz2J39zYqvyRT6DDym818/SfXUbmFoGvWmUspPxA4cHyF9ZX9fNz951cVnt7Gutg+rVZT6Qdsl3U9RUVFkZSURG5uLsXFxWzevJn+/fs719vtdp544gmGDh3Kiy++6DzaMJlMLF26lP379wOwfPlyBg8e7KqYwo0oDjtFackcsjbH4eZndhvaR6IU5FB65he1owhxXS49ooiNjWXChAlYrVbGjBlDREQEkydPZtq0aZw7d47Dhw9jt9v56quvAAgPD+e1117jnXfeYdasWZSUlNC6dWvmz5/vqpjCjdgz03AUF3DQcpvaUW6avlV30BkoOLwdcN+jI1E/uPTK7JiYGGJiYsotW7JkCQBdunQhNTX1qtv17NmT9evXuzKacEO2k3tBpyfV2lTtKDdNY/RB37IrhUd2oCXG7Y+QhGeT307hNmwn9+HTKpxSDGpHuWGX534ym/0I6jEQe+FF2uvPqR1LiErJXE/CLTguZqJcPIep972QbFU7zg27cu4nAzbeauLDbV4n+MXm/kdJwnPJEYVwC7bTBwAwteuucpKaY0VPg469iTCcRIdd7ThCXJMUCuEWbKcPoGnYBENgqNpRapRv576YtFY6GTLUjiLENUmhEHWeYivFnnEEfYsItaPUOJ/WXSh0GOn2v2nThaiLpFCIOs3P3wffwpNgtxHUpbfacWqcRqfngKUlXYyn0Uv3k6ijpFCIOs3bS8/q/7uWUkXP+PeOqx3HJfZaWuOtsRJmSFc7ihBXJYVC1GmKotDZkM4v1ibY0KkdxyV+tTWhwOFFd+l+EnWUFApRp1lz0gnSFXDE2kztKC7jQMsBS0vCjadxWEvVjiNEBVIoRJ1WdGwPgEcXCoC9llZ4a2wUH9undhQhKpBCIeq04mN7OWdvSK7jxqeOdwdp/+t+Kjjyo9pRhKhACoWosxRbKSWnjnDEA+Z2uh4HWvZbWlH0624Um3Q/ibpFCoWos+wZR1HsVo+YBLAq9llaoVhLsJ1OUTuKEOVIoRB1lu1MChq9kWPWELWj1Io0Wwhakz+248nXf7IQtUgKhaiz7GcO4t2yE9Z6MnelAy0NbumF7dQ+FJtF7ThCOEmhEHWSoyAHR95ZfNp2UztKrWrQKQqsJdjOHFQ7ihBOUihEnXS5n95UzwqFT6vO4NUA2/FdakcRwkkKhaiT7GcOomnQCEPjFmpHqVVWhwa/sEjsp/bRuJE3fv4+akcSQgqFqHsUhx1b+iH0zcPRaDRqx6lVRoOO+dtBsRTz95c+wNurfozPiLpNCoWocxxZx8FSjK55F7WjqOIXa6hMPS7qlCoViqlTp7Jjxw5XZxEC4H8DuRr0zW5VO4oqHGhJsbYg3HAaxe6+t30VnqNKheKee+7h/fffJzo6mo8++oi8vDxX5xL1mC39EFpzGzTenj1tR2X2W1ph0lop/u2A2lGEqFqhiImJYfny5bz//vvk5OQwduxYZsyYwYEDlf8SJyYmMmzYMAYPHsyKFSsqrN+6dSsjR45kxIgRPPnkk1y8eBGAjIwMxo8fz5AhQ5gyZQqFhYU30DThjpTSQhxZx9E376x2FFUdtYZS5DBQcOQntaMIUfUxCofDwcmTJzlx4gQ2m42goCBmzZrFm2++edXnZ2ZmsnDhQlauXElCQgKrVq0iLS3Nub6goIBZs2bx4Ycf8sUXX9CxY0cWLVoEwOzZs3nooYfYtGkT4eHhvP/++zfZTOEuDHnHQXEQFH47ZrOf2nFUY0fHQWsLin75GcVhUzuOqOeqVCgWLlzIgAEDWLp0KcOGDWPz5s3ExcWxfPlyPv/886tus2PHDiIjIwkICMBkMhEdHc2mTZuc661WK7NmzSIkpGx6ho4dO3L27FmsViu7du0iOjoagNGjR5fbTng2++mDlCp6HvznUWKeTVA7jqr2W1rhKCnAnpGqdhRRz1Xp3Lvc3FyWLFlCWFhYueUmk4m33nrrqttkZWVhNpudj4ODg8t1VTVq1IhBgwYBUFJSwocffsgjjzzChQsX8PX1Ra8vi2Y2m8nMzKxWo4KCbrxvuz58iq3LbTz1237SrCHYPfRudtWRam2KxuiDPmMf5u53OJfX5Z9fTfH0Nrpb+6pUKOx2e4UiMXXqVBYtWkTfvn2vuo2iKBWWXe2c+Pz8fJ588knCwsK47777rloUqnsufU5OAQ5Hxf1fj9nsR3Z2frW3cyd1uY2OS9nYLpwj1Xq72lHqBBs6GnToSUHqTpSeD6LR6ur0z6+meHob62L7tFpNpR+wKy0UM2fOJDMzk927d5Obm+tcbrPZOH688hvdh4SEkJz8+yyYWVlZBAcHl3tOVlYWjz76KJGRkbzwwgsABAYGUlBQgN1uR6fTkZ2dXWE74Zls6YeAsoFcUaZB2B0UHPoB+9mj9fZ0YaG+SgvFmDFj+PXXXzl69KhzzABAp9PRvXv3Sl84KiqKRYsWkZubi4+PD5s3b+bVV191rrfb7TzxxBMMHTqUJ5980rncYDDQs2dPNm7cSExMDPHx8fTv3/9G2yfciP3MQXR+QWTmNlQ7Sp3h064b6L2wHf9ZCoVQTaWFokuXLnTp0oU+ffo4B52rKiQkhNjYWCZMmIDVamXMmDFEREQwefJkpk2bxrlz5zh8+DB2u52vvvoKgPDwcF577TVmzpxJXFwcixcvJjQ0lLfffvvGWyjcguJwYEs/jG9YJJysX9N2VEZr8ELfsiu2E3tQ+kxQO46opyotFE8//TT//Oc/mTRp0lXXJyYmVvriMTExxMTElFu2ZMkSoKwIpaZe/WyOZs2asWzZskpfW3gWx/nfwFKEqW1X+DZH7Th1ir7t7diO/4z93FEI6aV2HFEPVVooJk+eDMDLL79cK2FE/WU7UzY+4dO6C/CtqlnqGn3LCNAby6Ye7yqFQtS+Sq+jCA8PB6BXr16EhobSq1cv8vLy+Pnnn+nUqVOtBBT1gz39ENqgVugayPjEH2n0/+t++i0ZxWFXO46oh6p0wd0rr7zCkiVLOHbsGHPmzCE9PZ0XX3zR1dlEPaFYS7BnptX7aTsqo297O0rxJUpOH1E7iqiHqlQoDh48yKxZs9iyZQv33Xcf8+bNIz093dXZRD1hP5sKDju6ZlIorkXfoivojBQeSVI7iqiHqlQoFEVBq9Xy448/EhkZCUBxcbFLg4n6w3bmEOgM6Jp0UDtKnaUxeKFvGUFh6k8oDofacUQ9U6VC0bJlSyZPnsyZM2fo1asXzz77LB07dnR1NlFP2NMPoWtyCxq9Ue0odZq+7e3YC/OwZ/6qdhRRz1RpCo958+axZcsWbrvtNucFcaNGjXJ1NlEPOAov4LiQgdctV58KRvxO37IrGr0R27Gf0YfKBzVRe6p0RGEymejZsyeXLl3i0KFDREREXHcKDyGqwv6/aTtkfOL6NAZvTO17/O/sJ+l+ErWnSkcUb775JsuXLycoKMi5TKPRsG3bNpcFE/WD7cwhNN5+aINaqB3FLTS4tQ+FqT9hP3cUfVM5RV3UjioVii+//JLNmzdXexoPISqjKI6y8YlmndFoqnwPrXrN1K6H8+I7KRSitlTp3RkaGipFQtQ4R+4ZlOJLcv1ENWiN3uhbdsN2fJdcfCdqTZWOKO644w7mz5/P3Xffjbe3t3N5587yBhc3zv6/aTuCI3qj93evG7moSd+uV9ncTzL1uKglVSoU69atAyh3S1IZoxA3ws/fB2+vsl+7s1mpGBo3577Z253rE98aqVY0t6FvEQEG77Kzn6RQiFpQpULx9ddfuzqHqCe8vfTEPJuAARvzGqUQ1HsY/KJ2Kvei0RvRt+pWdvZT34fRaKv0NhbihlVpjKKwsJA5c+bw5z//mby8PF555RUKCwtdnU14sLb6LAwaB6Y2XdWO4pb0bXuhlBZgT5e5n4TrValQzJ07Fz8/P3JycvDy8qKgoIBXXnnF1dmEB4C19ScAACAASURBVOtoyMCmaPFuKV0n1WGx2jGb/QjtHoXWywSnkq+/kRA3qUqF4siRI8TGxqLX6/Hx8WHBggUcOSKfZMSNCzOc5bgtGK3R+/pPrscuF4bLX0aDjphnExjx/EaSLoVSkrYLxWZRO6bwcFUqFFpt+afZ7fYKy4SoKn9NEc30F0i1NlU7Sp13uTBc/rrS7tI2KKVF2E6nqJRO1BdV+mt/++238+abb1JSUsIPP/zAU089Re/evV2dTXioWwxnAUi1hqqcxL39amuC1uSP7dhOtaMID1elQjF9+nRMJhN+fn688847hIWF8dxzz7k6m/BQYYaz5Du8ybAHqh3FrTnQ4tspCtvJfSjWErXjCA923UKxZcsWHnnkEZYuXcqZM2fw8/OjR48eeHl51UY+4WEUxUFHQwa/WJugoFE7jtvz7dwX7BZsJ/eqHUV4sEpPwI6Pj+f9999n2rRphIWFodFoSElJ4bXXXqO0tJR77rmn0hdPTExk8eLFWK1WJk6cyPjx46/6vOeff57evXszevRo534XLFjgnITwzjvvJDY29kbaJ+oYS+ZJ/LUlHLE2UzuKR/Bq3hFNg0CsaT9haH+H2nGEh6q0UCxbtoxPPvmEpk1/H3Rs164dXbt25YUXXqi0UGRmZrJw4ULWrVuH0Whk3Lhx9O7dm/bt25d7zsyZM0lKSio35pGSkkJcXBzDhw+/mbaJOqj4eNkn36MyPlEjNBot+na9sKZsQSkpQOPtq3Yk4YEq7XqyWq3lisRlbdq0obS0tNIX3rFjB5GRkQQEBGAymYiOji43BQiUHXHcfffdDB06tNzylJQU4uPjGTFiBNOnT+fixYtVbY+o44qO7+OMrRGXFJPaUTyGof0doNixHt+ldhThoSotFDqd7prrFEWp9IWzsrIwm83Ox8HBwWRmZpZ7zqRJkxg7dmyFbc1mM1OnTiUhIYHQ0FDmzJlT6b6Ee1AsxZScTpXTYmuQxWonJKwzhsbN0Zz8GT9/H7UjCQ/ksklirlZINJqqDV6+9957zv9PmjSJQYMGVWvfQUE3fvhtNnv+LKa12UaL1Y7RUPaBo/CXVAocdhmfqEGXr7MY5G0mxrQXXXEOZrP73xLA09+H7ta+SgvF0aNH6dGjR4XliqJgsVR+NWhISAjJyb9PL5CVlUVwcPB1A+Xn57N27VomTpzo3JdeX716lpNTgMNR+RHP1ZjNfmRn51d7O3dS2200m/2cF4qNNf1EPz9vfrOZr7OVqK7dljbEmPZScPAHbB2j1Y5zUzz9fVgX26fVair9gF3pX+AtW7bc8I6joqJYtGgRubm5+Pj4sHnzZl599dXrbmcymVi6dCndu3ena9euLF++nMGDB99wDlFXKIQZMvBpHY4989pdmuLGXHD4kmYNodPB7/G65Z4qH70LURWVFopmzW68iyAkJITY2FgmTJiA1WplzJgxREREMHnyZKZNm0aXLl2uup1Op+Odd95h1qxZlJSU0Lp1a+bPn3/DOUTdYNbm01hXgE/bbrBT5iZyhWRLG9rn/IQh5yS6xq3VjiM8iEsnso+JiSEmJqbcsiVLllR43htvvFHucc+ePVm/fr0ro4laFmbIAMDUrjsgU064wj5LK8b5J2P9NUkKhahRMrOfqBWdDOlk2/0wNGqidhSPVax4YWp/G7a0JLmftqhRUiiEyxmw0cFwjiNyWqzL+XUZgFJ8CfuZg2pHER5ECoVwuXaGTIwaO4eszdWO4vFM7Xug8fbD+sv26z9ZiCqSQiFcrrPhDBZFxzGr+5/fX9dpdAb07SOxndiLUiq3KxY1QwqFcClFUbjVkM6v1iZYXXvuhKDsAsfg3tHgsOGduU+u1BY1QgqFcClrbgaNdQXS7VRLjAYdY948QLqtEcn/XY+3lxRncfOkUAiXKkrbAyDTdtSyn0vb0Vp/Hsv5M2pHER5ACoVwqeK03Zy1NSTXIdNf16ZkSxvsiob8A9+oHUV4ACkUwmUUSzHFp45wWLqdal2B4sNhazMKDnyL4rCpHUe4OSkUwmVsGYfBYeOwdDupIqm0A/bCPGyn9qsdRbg5KRSiRvn5+2A2+2E2+2HIPIzGy8Rx2/VnDRY174i1GTrfQKxHvlM7inBzckqEqFHeXnpink1Ag8KrAT/SJPw2HGfl84gaHGjx63oXeTvW4ijIQesbpHYk4abkHSxcoqXuPH7aEhp0uF3tKPWaX7e7QQHr0R/UjiLcmBQK4RJdjKexKxp82nVXO0q9ZggIRtfsVqxHf0BxONSOI9yUFArhEuGG0xyzhaDzkdNi1WboNAClIAd7ukwUKG6MFApR4xprLxGqv8hBSwu1o9R7Fqud0J790Zr80Rz7Qab0EDdECoWoceHGsquBU+T6CdUZDTpGPLeRr3JaUnA0GX3JBbUjCTckhULUuHDDaTJsAeQ6/NSOIv7nx9JbALi0Z7PKSYQ7kkIhapS9KJ92+ixSrNLtVJdccPhyyNqMS/u2otitascRbkYKhahRRWm70WoUDlmk26mu+aEkDEfRJWzHd6kdRbgZKRSiRhWm/sQFu4lT9sZqRxF/8IstFENgKJbDX6sdRbgZlxaKxMREhg0bxuDBg1mxYsU1n/f888+zbt065+OMjAzGjx/PkCFDmDJlCoWFcqcud6BYiik+vo8D1pYoaNSOI/5AQYNfj2gcmWnYz59UO45wIy4rFJmZmSxcuJCVK1eSkJDAqlWrSEtLq/CcJ554gk2bNpVbPnv2bB566CE2bdpEeHg477//vqtiihpkO7UfxW5ln6WV2lHENfh1vQv0XlgOyqC2qDqXFYodO3YQGRlJQEAAJpOJ6OjoCgUhMTGRu+++m6FDhzqXWa1Wdu3aRXR0NACjR4+usJ2oO66cBFCXsR9dgwB+s5nVjiWuQefdAEPHvtjSfsJRlKd2HOEmXFYosrKyMJt//4MRHBxMZmZmuedMmjSJsWPHllt24cIFfH190evL5is0m80VthN1x+VJAO9/9nNyj+yiQcfeKDL0VacZw+8BhwProW1qRxFuwmWzxyqKUmGZRnP9fusb3e5KQUE3Pm2E2ez55/67oo2dDBl4aWw06HQHbD1e468vaobFaiekfXvO3dKTktRvaXznn/A2mdSOVYGnvw/drX0uKxQhISEkJyc7H2dlZREcfP37EgQGBlJQUIDdbken05GdnV2l7a6Uk1OAw1Gx4FyP2exHdnZ+tbdzJzXdxsu/8F2NJylweOHd8lZACkVdZTToiHk2gXb6QKb552NJ3U5+izvUjlWOp78P62L7tFpNpR+wXdZHEBUVRVJSErm5uRQXF7N582b69+9/3e0MBgM9e/Zk48aNAMTHx1dpO6EePXY6G89wwNICjVandhxRBcdswZyyBXHx50QURWaVFZVzWaEICQkhNjaWCRMmMGrUKIYPH05ERASTJ08mJSWl0m1nzpzJ6tWrGTZsGMnJyTzzzDOuiilqQCdDOt4aG/utcraT+9DwbUknrDkZ2E7uUzuMqONceoe7mJgYYmJiyi1bsmRJhee98cYb5R43a9aMZcuWuTKaqEG3ef1GvsObX6yhakcR1bDX0pq/NvwFy97/om/VvdpjgaL+kNNTxE1xlBYRbjjNXksrHPLr5FYcaAm4YxSO7OPYM46oHUfUYfLOFjel8OhODBoHyaVt1Y4iboBv14FofBpi2fdftaOIOkwKhbgpBYd+4Lzdl5Myt5Nb0uqNGCOGYE8/jD3rmNpxRB0lhULcMEdRHsW/pbDb0gZkbie3Zeh0J3g1wLJXjirE1UmhEDfMdnwXKA52l7ZRO4q4CRqjD8bwwdhO7pXJAsVVSaEQN8yaloQxpA2ZjgC1o4ibZAwfDF4NKE1ed/0ni3pHCoW4Ifa8DBxZx/EN76d2FFEDNF4NysYqTu2XsQpRgRQKcUOsqT+ARodv+AC1o4gaYgwfjMbbj9JdclQhypNCIarFz9+HxoE+OI7twNThNvS+0u3kzixWu3Oa+OCmZvwiR2FPP4Tt7FG1o4k6RAqFqBZvLz0vvrgYe+FF3tl947P0irrh8iSBl78Cew1FYwrAsmvtVWdyFvWTFApRbZFeaVx0+HDE2lTtKKKGaQ1eGHuMxH7uF2wn9qgdR9QRUihEtdjyc7nVkM7Ppe1kyg4PZQjrj7ZRU0p3rkax29SOI+oAeaeLask/8C1ajcLO0vZqRxEuotHq8Or9AMqlTKyHv1Y7jqgDXDp7rHB/fv4+eHuV/ZooDjun920lzRpCtsNf5WTClXQtItA160zpngQMt/RB49VA7UhCRVIoRKUu3xMboLPhDI/5ZfJDqdxIytNpNBq8Ih+gaO1MSncn4B31kNqRhIqk60lUWX/vI+j8Ajlgaal2FFELdEEtMYQNwHpoK/ac02rHESqSQiGqJESbR5jhLP49omUQux7x6jUGjVcDSrf/P7llaj0m73hRJf29U7EqWvy7D1Y7inChKy/AM5v98A8249X7T9gzf8X2y49qxxMqkTEKcV0+Ggu3ex1nj6UNHRs0VDuOcKHLF+BdlvjWSPS39EGX+j2lO1eX3TLVWy60rG/kiEJcV2+vNLw0Nr4vCVM7iqhlFqud4OCGNIl5AqW0EPvuNWpHEiqQIwpRKcVu5U6vw6RZgzljD1I7jqhlVx5hDPPpTPSh7/Bp0R19y24qJxO1yaVHFImJiQwbNozBgwezYsWKCuuPHDnC/fffT3R0NC+++CI2W9lVoPHx8fTt25eRI0cycuRIFi5c6MqYohL5B76jka6ILSVd1I4iVPZVcQQGc0tKvv8EpbRQ7TiiFrmsUGRmZrJw4UJWrlxJQkICq1atIi0trdxzZsyYwcsvv8xXX32FoiisXr0agJSUFOLi4khISCAhIYHY2FhXxRSVUBx28nas45QtiFSZ16nes6MjOOYplOJLlOxYqXYcUYtcVih27NhBZGQkAQEBmEwmoqOj2bRpk3N9eno6JSUldOtWdgg7evRo5/qUlBTi4+MZMWIE06dP5+LFi66KKSphO7YTW14mm4u7IPfEFgBeoe0wdrsX268/Yv1tt9pxRC1xWaHIysrCbDY7HwcHB5OZmXnN9Waz2bnebDYzdepUEhISCA0NZc6cOa6KKa5BURxY9v0Xg7klB60t1I4j6hBjj5FozW0o+e4jHPnn1Y4jaoHLBrOvNpe9RqOp0vr33nvPuWzSpEkMGjSoWvsOCrrx0/fMZr8b3tZdVKWNBalJFFzIoPGoZ1CO5tZCKuEugps0wjp2OmeWTsf2/RKaPjIHja5m/5R4+vvQ3drnskIREhJCcnKy83FWVhbBwcHl1p8///unkezsbIKDg8nPz2ft2rVMnDgRKCsoen31YubkFOBwVP+mK2azH9nZ+dXezp1UpY2Kw0bh1uVoA0Jp0CkK+G/thBN1nsVqx2jQYWjUBPOwJ8iKX0j6pmV49RpTY/vw9PdhXWyfVqup9AO2y7qeoqKiSEpKIjc3l+LiYjZv3kz//r9PJtesWTO8vLzYvbusnzM+Pp7+/ftjMplYunQp+/fvB2D58uUMHixXA9cWP38fvM/sRLl4DvPgP6PR6tSOJOqQK++I9+B/cvDrNgjLvv9ilZsceTSXFYqQkBBiY2OZMGECo0aNYvjw4URERDB58mRSUlIAWLBgAfPmzWPo0KEUFxczYcIEdDod77zzDrNmzWLo0KEcOnSIGTNmuCqm+AMjVn7buIw0awjjPkhXO46o44Lu+WvZeMU3H2LPld8XT+XSC+5iYmKIiYkpt2zJkiXO/4eFhfH5559X2K5nz56sX7/eldHENeT9FI+ftoQP8+9CznQS16M1eOEzeCpF62dRvPmfNBj1ikzx4YFkCg/h5Ci8wMWdiewpbc0pe2O14wg3ofUNxOeeaSgFuVi+/TeNA72dkwr6+fuoHU/UAJnCQziV7lgBisJ/i7urHUW4icuzzWLuTr5jCtmJi/h8dhzLCvuhoCHxrZHUrWFbcSOkUAgArCf2YPstmUZ3jidnnQxgi6r542yzy+5/mNu+WU6B4s26ottVTCZqknQ9CRRLMaU/LkPbqDkBkSPUjiPcWMM7RvFtSScGeKdyj3eK2nFEDZFCISjdtRalMA/v/hNr/MIpUb9oNBrii3ryc2lb7jXt48L2iierCPcjfxXqOVv6YayHtmG49S50Ie3VjiM8gIKGlYVRAPT67lOMBcV43TZK5VTiZsgRRT3mKL5EydcfoG0YglfvsWrHER5EQcvKwih8I+7Csjuekp9WyT233ZgcUdRTiuKg5JsPUSxFNB3/Cl4h5utvJEQ1KGhpGP0YWoORS7u/xGi7hN89UygosqkdTVSTFIp6yrL/S+xnDtJ4yGOMmb8fKJsyJfGtkeoGEx7Fy2hgzBYzd3n3YOThH7EXXEB/55NyUZ6bka6nesh2ci+WXZ+jb3s7fj3uUTuO8Hgavi4J55OCfpSk/0LhupnYs0+oHUpUgxSKeqbg1FFKvv43xpC2NL//mXJTvwvhSnstbWg64TVQFIq+mIsl9bur3m5A1D1SKOoRx6Vsctb9g/MlBp470oORf/9K7UiinvFu2h7T/bPRNelI6fcfU7L1PRwlcu12XSeFop5w5GdTtGE+it3Ov/PvJl+ROXhE7bNY7YS0aErLCTMJHPgwtpN7KVrzIrYTe9WOJiohhaIecOSdpeiLeSiWIpqMe4ksR0O1I4l66vKUHyNm/JdH1mpp9tf5aHwaUrz5nxRv/pfcWrWOkkLh4bxLMinZ8AZa7DR7ZA7ezTqoHUkIJ6+Q1pjum4mx1xhsZw5SuPoFLvywBsVaonY0cQUpFB7Gz9/HOcWzd+Zezq98mfOXrMxJv/N/p8EKUbdodHq8ug2nwZ/m4d22Oxe+/4zi1c9jPPEDvg0MascTyHUUHsfbS8+oZ9cx0rSbAd6peLe8lbf2d5YxCVHnaX2DCPnTc0x9finDrXvosPkjdEnr0XcejCHsTjRG+R1WixQKD1Ny5igz/DcQqs/jm5JO/PWhmeTv26B2LCGuynk/iyucsJl5N/8eOurP8rdW5yj5aRWle77AcEs/DLfeiS6gqUpp6y8pFB7CUZKPZXc8+Ye+xlvjwwf5d3HY2pxHZTZYUYf98X4Wv88MoOGorSmNH3gMR/ZvXNz5BYWHv8Z6cDO60I7oO0RhaNMTjVcDdYLXM/JXxM34+fvg7fX7j63oYh65P8ZjObgFrKX43z6U575qSCnStyvcn9GgI+atQ0A7fDVNWTzaiwt7tlD6/ceUbv9/6JqHo2/dA33LrmhNAWrH9VhSKNyMt5eemGcTaKbLoY/Xr/T1O4ViLaFBpyga9fsTRnMLSr9KuP4LCeFmChQfTLcPp+Edo7CcO07Boe0UHNlB6ff7KQWMTdpiaNEZR+Nb0DXpgMbgrXZkjyGFwg34+fvgZdRhPX+aC9u/Yrr/V7TQ52JRdDTo1J+Xf/Al48dG8OMeEt9qoXZcIVymfFdVAF8s+DdPPP9/6Ww4Tdjps7TP+hIciaDRoA1sjnfzMEzN2+MV0rbsQ5RdQ/6lYlXb4I5cWigSExNZvHgxVquViRMnMn78+HLrjxw5wksvvURBQQE9e/Zk9uzZ6PV6MjIymDFjBjk5ObRp04YFCxbQoEH96otU7FYcF9KxZ5/g0oXjnN6fTCNdEQA2GrOu8HZ+trRldcyfyPhWjiBE/aTRaMiwNyLD3ogtJREkvBBN5sE92DN/xZ55jOLUHynavwUAu6LB29wMxTcEbaOmaP2D0fgHo/UzozEFoNHK1QLX4rJCkZmZycKFC1m3bh1Go5Fx48bRu3dv2rf//S5qM2bMYO7cuXTr1o0XXniB1atX89BDDzF79mweeugh7r33Xt577z3ef/99ZsyY4aqoqlBspSglBSjFl1AK83AUXcBxKQvlUjaOi+dw5J0DxQ6ArUEAJ2xmtpQ0YfrfJ/H0rB9UTi9E3aQ1eqNv0QV9iy4ANG7cgEefW0EzXS7N9BcY0cgbe046lpN74cobKWl0aHwboTU1QmNqWFY4fPz/9+WHxsv3f18mNEYT6I31akJNlxWKHTt2EBkZSUBA2QBTdHQ0mzZt4qmnngIgPT2dkpISunXrBsDo0aP517/+xdixY9m1axfvvfeec/nDDz9crUKh1Vb/B6g4bOQf/B7bhYuXl4BzYksFFAXnAsVR9lhRQHGU3bnL4Shb7rCjOOzgsIHdhmK3gt0CNmvZ1aa2UhRLMditFTJodAb0/o3RteiIV/e7MDRujtHcAn3DYOa8tgWMoPcLJLhR+fPJ5bE8lsdlrna6raZhCBmEkAE88ad7eHTuZrQoNNQWMnNcByjMxV5wAVt+Ltb8C9gK8lDOnUOxlXJNGh0YvNDovdAYjKDzAr0BdAY0emPZvzoD6PSg1ZX9X6sDrY6S0lvQ+jS/9mur4Hp/MzWKi+b5/eCDDygqKiI2NhaANWvWcODAAV599VUA9u7dy/z58/n0008BOHnyJI899hjLli1jzJgxfP/99wDYbDa6devGwYMHXRFTCCHEdbisU+5q9efKQ7Vrrb/edkIIIWqXywpFSEgI58//PhNkVlYWwcHB11yfnZ1NcHAwgYGBFBQUYLfbyy0XQgihDpcViqioKJKSksjNzaW4uJjNmzfTv39/5/pmzZrh5eXF7t27AYiPj6d///4YDAZ69uzJxo0byy0XQgihDpeNUUDZ6bEffPABVquVMWPGMHnyZCZPnsy0adPo0qULqampvPTSSxQWFnLrrbcyb948jEYj6enpxMXFkZOTQ2hoKG+//TYNG8o9FIQQQg0uLRRCCCHcn1xhIoQQolJSKIQQQlRKCoUQQohKSaEQQghRKSkUf5CVlcVjjz3GqFGjGDduHGfOnFE7kkscPnyY8PBwtWO4xO7du7n//vsZOXIkf/7zn0lPT1c7Uo1ITExk2LBhDB48mBUrVqgdp8a9++673Hvvvdx7773Mnz9f7Tgu849//IO4uDi1Y1SLFIo/eO655xg4cCDx8fGMHDmSBQsWqB2pxhUXFzNnzhys1orzTXmCGTNm8Nprr5GQkEBMTAxz585VO9JNuzzJ5sqVK0lISGDVqlWkpaWpHavG7Nixg+3bt7N+/Xri4+M5dOgQW7ZsUTtWjUtKSmL9+vVqx6g2KRRXyM3NJTU1lXHjxgFw//3388wzz6icqua98cYbTJw4Ue0YLmGxWHj66acJCwsDoGPHjpw9e1blVDfvykk2TSaTc5JNT2E2m4mLi8NoNGIwGGjXrh0ZGRlqx6pReXl5LFy4kCeeeELtKNUmheIKp0+fpmnTprz++uuMGDGCadOmYTB41i1Ft23bRklJCUOGDFE7iksYjUZGjiy777LD4eDdd99l0KBBKqe6eVlZWZjNZufj4OBgMjMzVUxUszp06OCcSfrEiRNs3LiRAQMGqJyqZr3yyivExsbi7++vdpRqq7d3uPvyyy+ZN29euWWtWrXi8OHDTJ06lRdffJE1a9YQFxfHsmXLVEp5467WvrZt21JQUMAnn3yiTqgadq02fvLJJ1gsFuLi4rDZbDz++OMqJaw59WWyzF9//ZXHH3+c559/ntatW6sdp8asWbOG0NBQ7rjjDtatW6d2nGqTK7OvcOrUKe677z7n/FPFxcVERkayf/9+lZPVjDVr1vDBBx847xaYmppKWFgYK1aswNfXV+V0NaewsJApU6YQEBDAggULMBqNake6aevXryc5OZnXXnsNgPfeew9FUZz3d/EEu3fvZtq0abzwwgvce++9asepUX/5y1/Izs5Gp9Nx8eJFioqKGDVqFC+88ILa0apGEeUMHTpU+fbbbxVFUZQNGzYoDz74oMqJXOeWW25RO4JLTJkyRXnppZcUh8OhdpQac+7cOWXgwIFKTk6OUlRUpIwYMULZv3+/2rFqTEZGhtK7d29lx44dakdxubVr1yrPP/+82jGqpd52PV3Lu+++y8yZM3nzzTfx9fXljTfeUDuSqIbDhw+zbds22rdvz6hRo4Cy/vwlS5aonOzmhISEEBsby4QJE5yTbEZERKgdq8Z89NFHlJaWlnu/jRs3jgcffFDFVOIy6XoSQghRKTnrSQghRKWkUAghhKiUFAohhBCVkkIhhBCiUlIohBBCVEpOjxXiKvbt28dbb71FXl4eiqLQpEkTnn/+eTp06ODS/Xbs2JGkpCQCAwNduh8hqkMKhRB/YLFYePzxx/nPf/5D586dAUhISGDy5Mls27YNnU6nckIhapcUCiH+oLi4mPz8fIqKipzLRowYga+vL3a7nXnz5rF//34KCwtRFIW5c+dy2223ERcXh5eXFykpKZw/f56hQ4cSGBjIN998Q3Z2NnPnzuWOO+4gLi4OjUbDsWPHyM3NpU+fPrz00ksVJqBcs2YNn376KQ6Hg4CAAF5++WXatWtHcnIyb7zxBg6HA4DHH3+c6OjoWv0eifpFxiiE+IOGDRsyY8YMJk2axN13382MGTNYu3YtUVFRHDp0iKysLFatWsXGjRu57777yl31feTIEVatWsXatWv55JNPMJlMfPbZZ0yYMKHc81JTU/n444/ZuHEjx44dY9WqVeUy/Pzzz8THx7NixQri4+OZNGkSU6dOBWDRokX85S9/Yd26dbz++uv89NNPtfONEfWWHFEIcRV/+ctfGDt2LLt27WLXrl0sWbKEJUuW8Pnnn/PMM8/w2Wefcfr0aXbu3OmcZBFg4MCBGAwGzGYzJpOJfv36AdCyZUvy8vKcz7vvvvuc240cOZJt27bx8MMPO9d/++23nDx50nlvFICLFy+Sl5fH0KFDmTNnDl9//TVRUVH87W9/c/W3Q9RzckQhxB/s3r2bpUuX4uvry8CBA3nuuefYsGEDWq2WrVu3Oqctv/vuuyvMRfTHmWr1+qt/FrtynENRFLTa8m9Fh8PByJEjSUhIICEhgfXr17N27VoaNmzIuHHj+OKLL+jTpw/bt29nxIgRY7PyPAAAAUxJREFU5Ofn10TThbgqKRRC/EFgYCCLFy8mOTnZuSw7O5vi4mI2bNjAwIEDeeihh+jSpQtbt27FbrdXex9ffvklFouF0tJS1q9fz8CBA8ut79OnDxs2bCArKwuATz/9lD//+c9A2WR5R44cYfTo0bz66qtcunSJixcv3kSLhaicdD0J8Qdt2rThvffeY+HChZw7dw4vLy/8/PyYM2cOzZo1Y/r06cTExKDT6ejZsyebN292DixXlbe3Nw899BCXLl0iOjqa+++/v9z6fv36MXnyZP7617+i0Wjw9fXl3XffRaPRMH36dF5//XXeeecdtFotTz31FM2bN6/Jb4EQ5cjssULUsri4ODp06MCjjz6qdhQhqkS6noQQQlRKjiiEEEJUSo4ohBBCVEoKhRBCiEpJoRBCCFEpKRRCCCEqJYVCCCFEpaRQCCGEqNT/B3Vp3F+ZKP2LAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.ylabel(\"Density\")\n",
    "plt.xlabel(\"Samples\")\n",
    "plt.title(\"Tracy-Widom GOE\")\n",
    "plt.hist(GOE_samples, bins=70, density=True)\n",
    "plt.plot(sol.t, f1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GUE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall from the first chapter that the off diagonal terms of a GUE have standard deviation 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "trials = 30000\n",
    "GUE_samples = []\n",
    "N = 100\n",
    "σ = 1\n",
    "for _ in range(trials):\n",
    "    A = RM.Generate_GUE(N)\n",
    "    λs, V = np.linalg.eigh(A)\n",
    "    λ_max = max(λs)\n",
    "    GUE_samples.append(N**(1/6)*(λ_max/σ - 2*np.sqrt(N)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true,
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEXCAYAAACtTzM+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXRU9f3/8eedNRmyQZiEIIvskSWIoCAoNAoiYADXUvot1la+2h6l8lOUqrUUcamitNW6Ua39KigiS4wosolVE4Gg7GvYEyArgSyTWe69vz8CkQiEADO5s7wf53DIbPe+bjLJa+72uYqu6zpCCCEinsnoAEIIIYKDFIIQQghACkEIIcRJUghCCCEAKQQhhBAnSSEIIYQApBBEEJkxYwZjxoxhzJgx9OzZk+HDh9fdrqmpabIcI0aMYMWKFXW3v/nmG7p168aHH35Yd9+mTZsYNGgQuq4zZswYTpw4ccZ03n77baZOnRqwnAUFBUydOpXhw4czatQohg8fzqxZs/B6vQDk5+fTp0+fBnMtXLiQvn371n2fT/179NFHA5ZbBC+L0QGEOOXJJ5+s+/qGG25g5syZ9OrVq8lzDB48mDVr1jB06FAAvvzyS9LT01m1ahXjxo0D4LvvvmPw4MEoikJmZmaTZywsLOTnP/85kyZN4rnnnkNRFKqqqpg6dSp//etf630vz6dfv368+eabAUwrQoUUgggZPXv25MYbb2THjh3MnDmTnTt3Mm/ePLxeL8ePH2fixImMHz8egDfffJNFixZhsVho3749zz//PJMmTeLmm2/m5z//OQCvv/46x44d4/HHH683n8GDB/PCCy/U3f7yyy95++23ueuuu6iursbhcJCTk1M3nW7dupGTk0NsbCwzZswgOzubxMREEhMTiY2NBeDo0aNMmzaNgoICdF1n7Nix3HvvveTn53P33XczYMAANmzYgM/n49FHH2XevHns3buXnj178vLLL2My1V+Zf+utt7jpppu466676u5r1qwZf/rTn/jiiy/8/80XEUE2GYmQ4fV6SU9P54svvqBjx47Mnz+ft956i8WLFzNr1ixefPFFAFauXMnChQuZN28en376KW3atOH999/nl7/8JfPnzwdA0zTmz59f94n/dFdffTUHDx6kvLycnTt3Eh8fT4cOHUhLS+Pbb7/F4/GwadMmrrvuunqvmzt3Lvv372fJkiW88847HDlypO6xRx55hP79+5OVlcUHH3zAJ598wpIlS4DaTTs33HADS5YsYcCAATzzzDO8/PLLLFmyhNzcXDZs2HBGxtzc3DPmD5CUlMSvfvWrC/q+5ubmnrHJaMGCBRc0DREeZA1BhJR+/foBtZ+G33jjDb766iv279/Pjh07qK6uBiAnJ4ebb76Z+Ph4AP74xz8CoKoqM2bMYMeOHRQWFtKmTRs6dux4xjxsNhv9+/cnNzeXvLw8fvaznwGQnp7ON998Q1xcHD169CAmJqbe63Jycrjllluw2WzYbDYyMjLYuXMn1dXVfP/997zzzjsAxMbGctttt/Hf//6X3r17Y7VaueGGGwBo164dffr0qZt2UlISx48fPyOjrusoilJ3+1//+hdZWVkAlJSUsGTJkjPWKk7RNK3eY7LJSJwiawgipDgcDqB2E8zYsWMpKCigb9++PPTQQ3XPMZvN9f5Ynjhxgvz8fMxmM+PGjePjjz9mwYIFdWsHTzzxRN0n4w8++ACA66+/nnXr1rF69WrS09MBGDJkCN9//z3fffddXUk0xGw2A7V/gH86ZJimafh8PgCsVmu9vFar9bzT7tOnD2vXrq27fe+995KZmUlmZiYlJSVomkZ8fDw1NTW43e56ry0tLSUhIeG88xCRRwpBhKQtW7bQokULfv/733P99dfz5ZdfArVrAQMHDmT58uVUVlYC8Morr/Duu+8CcOedd7JixQq2bt3KsGHDAHjmmWfq/pj+4he/AGr3I3z77bccPny4bsd227ZtAVixYgVDhgw5I9P111/P4sWLcbvduN1uPvvsMwBiYmLo3bs3c+bMAaCiooLFixczcODAi17+3/3ud3z++ecsXrwYVVXrlv3UPE0mE82aNaNv37785z//qXtdYWEhS5cuPWt+IWSTkQhJgwYN4uOPP+bmm28mOjqatLQ0WrRowYEDBxgyZAh5eXl1f9w7d+7M008/DUBiYiI9e/akU6dODX4Sb9u2LT6fj0GDBtX79H799dezbNkyOnXqdMZrxo0bx8GDB7nllltISEigffv2dY/NnDmT6dOns3DhQjweDxkZGdx2220UFBRc1PK3atWKefPm8eqrr/L2229jNptxuVz06NGDjz76qG4NYObMmTz77LOMGjUKk8mE2Wxm0qRJ9O/fv25ap/YhnM5sNrNw4cKLyiZClyLDX4tIUlZWxh133MGcOXNISUkxOo4QQUU2GYmI8dFHHzFy5EgmTJggZSDEWcgaghBCCEDWEIQQQpwkhSCEEAKQQhBCCHGSFIIQQgggxM9DOHasCk278H3iiYkxlJZWBiBR8JJljgyyzJHhYpfZZFJo3rzZOR8P6ULQNP2iCuHUayONLHNkkGWODIFYZtlkJIQQApBCEEIIcZIUghBCCEAKQQghxElSCEIIIQApBCGEECdJIYiIoKs+oyMIEfRC+jwEIc5HqyjBve5jfHnfYW6eQrMOaUR36I2j29W4PRoVJ1xGRxQiaEghiLCkq15KV71H1dpPAQXrFelYvScoWr+SqO+/ILumM7+c/jwVRgcVIohIIYiw5M6ei3f7l1i6DMR+9e2YYhJxOmO57+FFjIjeyE3Rmyld/g76lXfWu0SmEJFMCkGEvNi4aKLsP76VK3fkULH9S+IHjCbxxrvrPVfDxBLXlVhRSV/3GTYP2K6RUhACpBBEGIiyW8h4OBOA5qZKHo37lOaXdea3n8WifpZZ97ysl05dSF5hsasvowe1oeL7zzC1aIO1y0ADkgsRXOQoIxE2TGhMaPY1JkUn6dbJqJgbeLZCy5snYnJ2wL32Y3Sfu8lyChGspBBE2LjWvpuO1mI+quqPtXmr8z5fUUzYB4xDryrDs2lpEyQUIrhJIYiwYELjhqit7PO1ZL2nQ6NfZ0nphqVDPzwblqBVHQtgQiGCnxSCCAu9bQdoaa5kpasncGE7iO397wJNw71uQWDCCREipBBEyNN1naFRWziqxrPF2/aCX2+KS8Lacxi+Xd+ilhwIQEIhQoMUggh5rn0baWM5xipXd/QLXDs4xd7nFrDa8WxZ5ud0QoQOKQQR8spzFlOuRZPr6XjR01DszbB2GYhvzxr0msi6Pq8Qp0ghiJCmFu+nZv9mVtd0P89hpudn7X4DqD68O7/2UzohQosUgghp3p3/RbHYyK7pcsnTMrdog7lVVzzbVqHrmh/SCRFapBBEyNI1H76963B0vRo3Nr9M09r9BvSKYtT8LX6ZnhChRApBhCy1YBt6TQUx3a/z2zQtHfqhRMfh2brKb9MUIlTIWEYiZHnz1oDNgaNTH+DIBb/e41VxOmPr3Vfj9uFJHYLnh0/RKooxxTr9lFaI4CeFIEKS7vPg278ea8erUSzWi5qGzWquGxTvlKyXxmC94md4fvgU765vsfcd64+4QoSEgG4yysrKYuTIkQwbNow5c+ac83mrV6/mhhtuCGQUEWZ8BzeCtwZLpwF+n7YpJhFzSld8e9f6fdpCBLOAFUJhYSGzZs1i7ty5ZGZmMm/ePPLy8s54XklJCX/9618DFUOEKd+eNSjRcZhbXxGQ6Vs6Xo127DBqWUFApi9EMApYIWRnZzNgwAASEhJwOBwMHz6cpUvPHFHyySef5IEHHghUDBFmYuOiSYwzox7cSGyP60hKjg/IfCwd+oGiyFqCiCgBK4SioiKczh93yCUlJVFYWFjvOf/3f/9H9+7d6d27d6BiiDATZbfw9LQ30FUvf1mln7EPwF9MjgTMKan49q5F1/WAzEOIYBOwncpn+yU6/TKFu3btYtmyZbz77rscPXr0ouaRmBhz0fl+enRJJAiXZe5hzee4Fs1+X0u/T/v0I49OpF1PydK3aOYrpVnrxg+pbbRw+TlfCFlm/whYISQnJ5Obm1t3u6ioiKSkpLrbS5cupbi4mNtvvx2v10tRURHjx49n7ty5jZ5HaWklmnbhn96czliKiysu+HWhLFyWuWWig1TrETZ523Khw1w3xulHHsUoLp5OUHDvyqHa6v/yCYRw+TlfCFnmxjOZlAY/SAdsk9HAgQPJycmhrKwMl8vFsmXLGDx4cN3jkyZN4osvviAzM5O33nqLpKSkCyoDEZncBbtwmDxs87QJ+Lwq9Wh2+1pRtT1bNhuJiBCwQkhOTmby5MlMmDCBsWPHcsstt5CWlsbEiRPZvHlzoGYrwlx13veousIuX0qTzO8HT3u8ZUfQSg82yfyEMFJAT0zLyMggIyOj3n2zZ88+43lt2rRh1SoZKkCcX/WeH9jnc+LS/TN20fls8rRjXMxafPtyMbds3yTzFMIoMpaRCBladTmewn1s917WZPOs0qOIatMN38FNTTZPIYwihSBChnqodlNjUxYCgKNzX7TSA2hVx5p0vkI0NSkEETJ8hzZjjmlOgdq8Sefr6Nz35PxlLUGENykEERJ0TcWXv4Xojn0IxOGmDc47oTWWuJaYj27B6YzF6YwlNi66STMI0RSkEERIUIv2gKcaR+ermnzedpuF1UUtOLbzB8Y+vJCMhzOJsstAwSL8SCGIkKDmbwFFIbpDmiHz3+ptg13x0dlSeP4nCxGipBBESFAP78DU8nLMUc0Mmf9ubys8upketnxD5i9EU5BCEEFP97lRi/ZiTkk1LIMXC7u9rehuLQDkrGURnqQQRNBTC/eA5sPS2rhCgNrNRk5zBUmmE4bmECJQpBBE0FOP7ARFwdyqq6E5tp08/6G7TS6aI8KTFIIIeuqR2v0His3YQz2PaTEcVePpZjlsaA4hAkUKQQQ13edBLdyDOaWb0VEA2OlNobO1EN3nNTqKEH4nhSCCmloUHPsPTtnpTcGmqNQU7DQ6ihB+J4Uggpp6eEdQ7D84Jc+bjKoruPZuNDqKEH4nhSCCmnpkB6bE9ig2h9FRAHBj44CvJa79ck0PEX6kEETQ0n0e1KK9xHZKqxtDKBjs9KXgPrIH3V1ldBQh/EoKQQQttWgvqF5eXFFJxsOZddc6NtpOb2vQNXyHtxsdRQi/kkIQQUs9sgNQ2OtLMjpKPQd8LVFs0aj5W42OIoRfSSGIoKUe3Y0tqX2TXS6zsTRMRLfvga9gm9FRhPArKQQRlHRNQy3aQ1Sb4Dj/4KeiO6ShnyhEqyg2OooQfiOFIIKSdqwAvDXYg7YQegPgk81GIoxIIYigpBbuBiCqbXCckPZTelwrzDEtsJTuliuoibAhhSCCknp0N0p0PJb44NqhfIrdZmFNaTxHtqwn4+HFcgU1ERakEERQUgvzMCd3RlGa9vrJF2K3rxVxphqSTceNjiKEX0ghiKCjVZejVxRjbtXZ6CgN2u1tBUBnq1xWU4QHKQQRdNTCPQCYk7sYnKRhpVoMx1QHXSxHjY4ihF9IIYigoxbuBrMFU8v2Rkc5D4Xdvla1w2HrcllNEfqkEETQUQvzMLfsgGK2Gh3lvHZ7WxFrqsFbfMjoKEJcMikEEVR0nweteD+m5ODef3BKnq92P4LrwBaDkwhx6aQQRFDRSg6A5gv6HcqnlGkxlKoxUggiLEghiKCiFuYBYE4KjUIA2O1LpubgVnRdMzqKEJdECkEEFbUwDyXWickRb3SURtvtbYXmqkQryzc6ihCXRApBBBW1eC/m5E5Gx7ggp/YjqHJ9BBHipBBE0NAqy9CrjmFOCq1CKNeaYWneqvb6z0KEMCkEETTUopMnpCV1NDjJhYtu1x3f0V2yH0GEtIAWQlZWFiNHjmTYsGHMmTPnjMeXL19ORkYGo0aNYurUqXg8nkDGEUFOLdoLJgumxHZGR7lgUe16gLsKrazA6ChCXLSAFUJhYSGzZs1i7ty5ZGZmMm/ePPLy8uoer66uZvr06fz73/9myZIluN1uFi1aFKg4IgRoxXsxtWwXEiek/VRU++7Aqct+ChGaAlYI2dnZDBgwgISEBBwOB8OHD2fp0qV1jzscDlatWkXLli2prq6mtLSUuLi4QMURQU7XVNTifSG3/+AUa3wSSmxL2Y8gQlrACqGoqAin01l3OykpicLC+qNCWq1WvvrqK9LT0zl27BjXXXddoOKIIKeV5YPPE5L7D04xp3RDPbpLxjUSIStgV/U42y/F2ca2HzJkCGvWrOHll19m2rRpvPTSS42eR2JizEXnczpjL/q1oSqYl/nEocNUA87UNKzNgzdnQ5p3vZLiXd+SoJRjcxq3HySYf86BIsvsHwErhOTkZHJzc+tuFxUVkZT049WvysvL2bJlS91aQUZGBpMnT76geZSWVqJpF/5pzOmMpbi44oJfF8qCfZlde7ZhcsRx24zvgB8/OGS9NMa4UBeoOqZ2dNbird9j69HckAzB/nMOBFnmxjOZlAY/SAdsk9HAgQPJycmhrKwMl8vFsmXLGDx4cN3juq4zZcoUDh8+DMDnn3/OVVddFag4IshpRXuIat2F08sg1CixTpRmLVCP7DQ6ihAXJaBrCJMnT2bChAl4vV7uuOMO0tLSmDhxIpMmTaJXr148/fTT3HfffSiKQufOnfnLX/4SqDgiiOmearTyI9jTBp//yUFMUZTa/QgFW9F1Pagv/ynE2QT0yuAZGRlkZGTUu2/27Nl1Xw8dOpShQ4cGMoIIAWrRPkDHflkX4IDRcS6Kx6vidMZyomtvSvJySDBXoMUkU3HCZXQ0IRpNzlQWhjt1hrK9dXBfMrMhNquZjIczeWhu7eU0Z740lyh7QD9vCeF3UgjCcGrRXkzxrTBHNTM6yiUr1mI5rkXT2VJ4/icLEWSkEIShdF2vPUM5RE9IO5PCHm+yXGdZhCQpBGEovaoM3XUCc1IHo6P4TZ4vmQRTNb5yWUsQoUUKQRhKLdoLgNkZumco/1SeNxkA14GtBicR4sJIIQhD/TjCaVujo/hNoRZPhRZFzcFtRkcR4oJIIQhDacV7MSWG5gin56awx5dEzUFZQxChRQpBGEbXNNTi/WG1/+CUPG8yvuPFaBUlRkcRotGkEIRhtPLD4HOH1f6DU/acus6yDGMhQogUgjCMdmqHcggPeX0uR9QETNExcsEcEVIaVQgPPvgg2dnZgc4iIkhsXDTWinxMdgdJnTuH3fDFOgpRba/AJ2sIIoQ0qhBuuukmXnvtNYYPH87bb79NeXl5oHOJMBdlt7BtXS7bK+MZ/UgWGQ9nGh3J76La9UA/UYRWdczoKEI0SqMKISMjg/fff5/XXnuN0tJS7rzzTqZMmcKmTZsCnU+EKc3rprX5GAd8LY2OEjDR7XoAcp1lEToavQ9B0zQOHDjA/v378fl8JCYmMm3aNF588cVA5hNhylO4D7Oic8CXaHSUgLEltwdbtFxnWYSMRg3HOGvWLBYuXEjbtm0ZP348f//737FarVRXV5Oens6UKVMCnVOEGffhPAAOhvEagmIyY27VVfYjiJDRqEIoKytj9uzZpKam1rvf4XBc0DWQhTilpmAXx1QHJ3SH0VECytI6FffBjWhVxzA1M+aymkI0VqM2GamqekYZPPjggwB110QW4kK4D+dxQA3ftYNTzClXAHI+gggNDa4h/PnPf6awsJD169dTVlZWd7/P52Pv3r0BDyfCk+Y6ga+8kAO+8L6Gtserkpzanf12B9Zje3A6h1Hj9slV1ETQarAQ7rjjDnbv3s3OnTsZPnx43f1ms5k+ffoEPJwIT1rxPoCwPsIIfryK2sSYFjhz1/CblSlkvTSGCqODCXEODRZCr1696NWrF4MGDSI5ObmpMokwpxbtBcXEoTA+wuh0eb5ketryiVOqjY4iRIMaLIQ//OEP/P3vf+fee+896+NZWVkBCSXCm1q8F5uzDZ7ScBrh9Nx2e2vHNepiPWpwEiEa1mAhTJw4EYA//elPTRJGhD9d11GL9uK44lqIkMPzC9TmVGtWuc6yCHoNHmXUs2dPAK655hpSUlK45pprKC8vZ+3atVxxxRVNElCEF/1EEbirsLfubHSUJqNjYo8vmc6yhiCCXKMOO33qqaeYPXs2e/bsYfr06RQUFPDEE08EOpsIQ2rRHgDsrbsYnKRp5flakWSuwHei1OgoQpxTowphy5YtTJs2jeXLl3Prrbfy3HPPUVBQEOhsIgypxfvAYsPmDJ9LZjbGqessy2U1RTBrVCHouo7JZOLbb79lwIABALhcciy1uHBq0R7Mzg4oJrPRUZrUqf0IrgNbjI4ixDk1qhDatWvHxIkTyc/P55prruHhhx+mW7dugc4mwoyu+tBKDmJyht8lM8/n1H4EKQQRzBo1ltFzzz3H8uXL6du3L1arlX79+jF27NhAZxNhRis9CJoPc1Ino6MYYpc3hV7H1qFVlGCKDe+T8kRoatQagsPhoF+/fpw4cYKtW7eSlpYmQ1eIC6YWh+8lMxtj96nrLBfIfgQRnBq1hvDiiy/y/vvvk5j445mliqKwcuXKgAUT4Uct2osSHY/SrIXRUQxxRE3A3CwB3+FtWFMHGx1HiDM0qhA+//xzli1bJsNXiEuiFu3BnNQRRVGMjmIQhajLe1K9dzO6rkfw90EEq0ZtMkpJSZEyEJdEr6lEP16IKTlyTkg7G1vbnuiu4yQox3E6Y3E6Y4mNizY6lhBAI9cQrr32Wl544QVuvPFGoqKi6u7v0aNHwIKJ8BEbF41WvotKILFrL6KdsUZHMkxMpzSOAX+bNYev3bVn+8sIqCJYNKoQFi5cCMDSpUvr7pN9CKKxouwWXv9PFkOjFP7nH9vxkEfWS2OMjmUIa0IyJWoMXa1H6wpBiGDRqEJYtWpVoHOIMHe5pZjDanM8RMYIpw3Z5W3FlbYDKGjojdtqK0STaNS7saqqiunTp3P33XdTXl7OU089RVVV1Xlfl5WVxciRIxk2bBhz5sw54/EVK1YwZswYRo8eze9//3uOHz9+4Usggp6uqbS3lLLf5zQ6SlDY5UvBYfLS1lx2/icL0YQaVQgzZswgNjaW0tJS7HY7lZWVPPXUUw2+prCwkFmzZjF37lwyMzOZN28eeXl5dY9XVlYybdo03nrrLT755BO6devGK6+8cmlLI4KStySfKMXL/jC/Qlpj5Z28PkJX6xGDkwhRX6MKYfv27UyePBmLxUJ0dDQzZ85k+/btDb4mOzubAQMGkJCQgMPhYPjw4fX2QXi9XqZNm1Z39FK3bt04ckR+QcJRTcEuAFlDOKlCj+awL4GuMhy2CDKN2odgMtXvDVVVz7jvp4qKinA6f/wDkJSUxKZNm+puN2/enKFDhwJQU1PDW2+9xa9+9atGBwdITIy5oOefzhmBR7oYtcxF3+2iUrNTokXe9/xcdvpSuM6+Eys+wL8/G3lvR4ZALHOjCuHqq6/mxRdfpKamhq+//pr333+f/v37N/gaXdfPuO9sJ+JUVFTw+9//ntTUVG699dZGxq5VWlqJpp05n/NxOmMpLo6sA/2MXGb34V0c8LUE5ESsU3Z4W5MetZ1O1tqrqPnrZyPv7chwsctsMikNfpBu1CajRx55BIfDQWxsLH/7299ITU3l0UcfbfA1ycnJlJSU1N0uKioiKSmp3nOKiooYP348qampPPPMM42JIkKM7q7CW5LPPtlcVM8ebzJe3USq9bDRUYSoc95CWL58Ob/61a/417/+RX5+PrGxsVx11VXY7fYGXzdw4EBycnIoKyvD5XKxbNkyBg/+cfwWVVW5//77GTFiBE888YScxh+m1KLaAe0OSCHU48XCHl8yqRbZbyaCR4ObjBYvXsxrr73GpEmTSE1NRVEUNm/ezDPPPIPb7eamm24652uTk5OZPHkyEyZMwOv1cscdd5CWlsbEiROZNGkSR48eZdu2baiqyhdffAHUXsNZ1hTCi1qYBygnNxmJ0+3wtmasY/3Jy2rajI4jRMOF8N577/Huu+/SunXruvs6depE7969efzxxxssBICMjAwyMjLq3Td79mwAevXqxY4dOy42twgRatEebEltcZfJCWk/tcPbGliPa99GaH210XGEaHiTkdfrrVcGp3To0AG32x2wUCI86JqGWpiHvU2q0VGC0hE1geNaNNV7NxgdRQjgPIVgNp/7urdnO4pIiNNpZYfAW0NUWxmz5+wUdnpTcO3biK5pRocRQgZSEYGjHt0NQFRbWUM4lx3e1miuSrSS/UZHEaLhfQg7d+7kqquuOuN+XdfxeDwBCyXCg3p0F0qzFlji5Aijc9npTQHAl78lYi8tKoJHg4WwfPnypsohwoyu66hHd2FO6SaHFDegUo/G1qojav4WuGq00XFEhGuwEC677LKmyiHCjF5Zgl5djrlVF6OjBD1Hpz6UZy9Cr6lEibr44ViEuFSyD0EExKn9B+ZWXQ1OEvwcXfqBruHL32J0FBHhpBCE38XGRWMt349id5DcVXYon489pRNKVCy+gxuNjiIinBSC8Lsou4Vd69aytSKB0VM+NTpO0FNMZsxt0/Ad2oSuqUbHERFMCkH4neqqJMVSzl5f0vmfLPB4VVr0GgDuKuI8R3A6Y4mNizY6lohAUgjC79z5OwGkEBrJZjVzz9sFqLrCf978kIyHM4myN2pkeiH8SgpB+F1N/nZ8uomDMqBdo9XoNvb4kulhzTc6iohgUgjC71wHt5OvtsDbuOsviZO2etvQ2lJOC1Ol0VFEhJJCEH6le924D+eR5002OkrI2eppA0B3WUsQBpFCEH6lFu4GzcduXyujo4ScYi2OIjVWNhsJw0ghCL9SD28Hk5m9XtmhfDG2etrS1XoUzV1tdBQRgaQQhF/5CrZjb90FD3JBnIux0dsOi6JRnbfe6CgiAkkhCL/RPS60kv1Et+9pdJSQtd/n5LgWTdWO74yOIiKQFILwG/XoTtA1oi+XQrhYOgqbPW2p3vMDuk+uSiialhSC8BtfwXYwW7BfJgPaXYqNnvboXje+QzLYnWhaUgjCb9TDOzAndcZktRsdJaTl+ZIxRcVgObIRpzNWhrIQTXhgKTYAABbRSURBVEbOHBJ+oddUopUexNZ3rNFRQp6GCUfXayja8DX/+9+2qJjJemkMFUYHE2FP1hCEX/iO7AR0zK1luGt/aJbaH4fJSxfrUaOjiAgihSD8Qj28Hcw2zEmdjI4SFqI7pFGjW+ltPWB0FBFBpBCEX6iHt2Nu1QXFLFsh/cFksbHVcxlptkOY0IyOIyKEFIK4ZFplKdqxAixt5HBTf/re04EYk5tU62Gjo4gIIYUgLpnv0GYAzG3TDE4SXrZ7W1Ol2ehn22t0FBEhpBDEJVMPbUZp1gJT89ZGRwkrKmY2eC6npy0fzeMyOo6IAFII4pLoqg9fwVYsbdNQFMXoOGEn19MBu+Kjaudao6OICCCFIC6JveIgeGto0fOaupOohP/s8yVRpjajcsvXRkcREUAOCRGXxHdoE6qucPdbB3FTe8x81ktjDE4VPnQU1ns60GLfRpq5TmCKjjM6kghjsoYgLolrz/fs9SXhxmZ0lLC13tMBdA3fnjVGRxFhTgpBXDStsgxP0UG2eS8zOkpYO6I2x5bUHu/uHKOjiDAnhSAumi+/9nDT7VIIARfT62doxXtRywqMjiLCmBSCuGjqoc2YY1twRE0wOkrYi037GZgseHesNjqKCGMBLYSsrCxGjhzJsGHDmDNnzjmf99hjj7Fw4cJARhF+pqtefPlbcHS6CpDDTQPN7IjD0qEv3t3Z6D6P0XFEmApYIRQWFjJr1izmzp1LZmYm8+bNIy8v74zn3H///SxdujRQMUSAqPlbwVtDs279jY4SMaxX/AzcVfj25RodRYSpgBVCdnY2AwYMICEhAYfDwfDhw8/4w5+VlcWNN97IiBEjAhVDBIh3Xy7Yoonu0MvoKBHDnJKKEpeMd/tqo6OIMBWw8xCKiopwOp11t5OSkti0aVO959x7770ArF+//qLmkZgYc9H5IvEEKn8ts676OHBwAzFdr0YxW/0yTdEwj1clKSmO8n43UbbqPeKV4xCfgs1qPuO58t6ODIFY5oAVgq7rZ9zn76ENSksr0bQz53M+TmcsxcWRdf0pfy6zL38LWk0lvta9/TI9cX42q5mMhzOJUWr4S4KJ9/7+Br99esYZP1N5b0eGi11mk0lp8IN0wDYZJScnU1JSUne7qKiIpKSkQM1ONCHfvlyw2LG0kc1FTa1Sj2aTpy3X2Paged1GxxFhJmCFMHDgQHJycigrK8PlcrFs2TIGDx4cqNmJJhAbF03LRAfawR9o1uUqklISjY4Ukb52p9LM5KFy02qjo4gwE9A1hMmTJzNhwgTGjh3LLbfcQlpaGhMnTmTz5s2Bmq0IoCi7hcl/fBu16jj/XGch4+FMoyNFpL2+JPb7WnJ8bRa6LldTE/4T0MHtMjIyyMjIqHff7Nmzz3je888/H8gYwo962w7g1U1s87QxOkoEU/iypjuXl/0X84ENWC+/yuhAIkzImcqi0XRNJc16kB3e1riRo4uMtMnTDku8E+8mOYdH+I8Ugmg01/4tNDdXk+vpaHSUiKdholnfUahHdxHrOVp3LQqPVzU6mghhUgii0So2raJKs7HF09boKAJoftVQqjUrn775OhkPZ5LxcOZZz0sQorGkEESj6O4qqneuZb2nAz7kj04wMNmjyXZ35UrbQZJMx42OI8KAFIJoFO+eNeg+D2vcnY2OIk7zZU13vJi5OXqj0VFEGJBCEI3i3fkNVmc78tUWRkcRp6nUo/mq5gr62PaTYj5mdBwR4qQQxHmpxwrQivcS2zsdGeo6+Kyq6U6NbmVk9Aajo4gQJ4Ugzsu782tQzMT2HGJ0FHEWLt3OlzU9SLMdwn047/wvEOIcpBBEg3TVh293NpZ2aZibxRsdR5zDVzWpVGp2yr76wOgoIoRJIYgG+fZ8h+46gbX7DUZHEQ1wY2NFTU9cezfgOyRDw4iLI4Ugzio2LpqWLWNQty7DltSOVn2uNTqSOI//1qRibZFCTfb76KrX6DgiBEkhiLOKslt49I+v4y0+yDt72zL6kU+MjiTOQ8VM4k2/RT9eiGfzF0bHESFICkGcU3rUNsq1aL73XG50FNFIjk59sFx+FZ7vP0GrLDU6jggxUgjirNyF+0m1HuG/NVegypnJIcPjVUkZNbH24OAfPsbpjCU2LtroWCJESCGIszq+Jgu3biHb3cXoKOIC2Kxmbns6h09PdKdqew6PT32FKHtAR7kXYUQKQZxBO15I5davyXF3waXbjY4jLsLKmh4c8rVgXLMcfJXlRscRIUIKQZzBvXY+itnKCldPo6OIi6Ri5r3K67ApPkqWvIau60ZHEiFACkHUox7djW9fLgnXjqFCl23PoaxQSyCr+iqq89bj3b7a6DgiBEghiDq6rlOzZh6KI4H4/qONjiP84Gt3KtEdeuP+7gPUsnyj44ggJ4Ug6vj25aIV5mHrdysmW5TRcYQf6Cg4Mx5AsTlwffE3tJoKoyOJICaFIIiNiyaxuR3f+gVYne1oPWiE0ZGEH2lR8aTcNRWqj6OufoOYZnI9bHF2UgiCKLuFf097Gt+xo/xtT1dGT/nU6EjCj2xWM3e+vI3/HO9PzcGtVH75b9nJLM5KCkHgOrCFIVHb+LqmG7t8rY2OIwJkvacjy109qdiwAs/6xUbHEUFICiHC6R4XxZ+8QokWyyfVVxkdRwTYElcfYnvfiOf7TNzrM42OI4KMnMIY4Wqy5+KrKOP9yuF4kG3L4U5HoeWo+3G53HjWLwKTCXufDKNjiSAhhRDBPFtX4tv1NQmD7uBAlsPoOKKJKIqJqMG/oUZT8axbgO6uwn7NXSgm2WAQ6eQdEKF8+3/Anf0+lvZ9aD74LqPjiCbk8aokJcfT9s7JxPUbgXfTUnxfvoburTE6mjCYrCFEIHtVPpWrXsee0pmUn09BMcloppHEZjWT8fCp/QdOrrNfw517czGVP0P00AcwxScbmk8YR9YQIoxavJ+yRS9Q5LbzyNYrGTN1qdGRhMG+cafS6uePo1WUUrXgKTzbV8thqRFKCiGC+A78QHXWs5isdt6ouJFKGatInOTo1Idmd8zAnNwZ99fv1p7VLBfYiThSCBFA13WOr/sM17J/YGp+Ga1//RwlWpzRsUQQ8XhVkju0p93d00i86bdoh7dT9dFU3OsWyL6FCCL7EMKcw+SiYtU7lO5ah6PL1SSNfUjGKRJnqL9fAZqbbuHFQUVU/ZCFd+fX2HqPwJo6BMUq751wJoUQpnRNw7vzv1St/YiaGjefufqyek039DVfkPXSGKPjiSB3TIuh+ahxxPcbSdmX71OT8wHeHz7B0XsYesfrMMW2NDqiCAAphDASGxeN3aJTsfkrjudk4i07TFT7HvxlYxfZRCQumM1qJmPWDqAf7c3tucGzlSvXZMKaxZhTUrF2HYSlfR+UqBijowo/CWghZGVl8frrr+P1evn1r3/NL3/5y3qPb9++nSeffJLKykr69evHX/7yFywW6agLpes6WskBqjZvYNu3y0gwVXPI14IVrsFM/+WDlPyQZXREEeIOqE7+Xfkz5v3xGtw7vqVy82pqvnobFIWoNqk4Ol9FVPte2Ft1wO2DihMuoyOLixCwv76FhYXMmjWLhQsXYrPZGDduHP3796dz5851z5kyZQozZszgyiuv5PHHH+ejjz5i/PjxgYoUNnRdR68oRj26u/YKZwVb0SuKwWTmsK8VH7ivZYe3NaCgKHLcgPAfR8sUfr7ICgylnbmUp0c62PzVKtoemgOARzcTd3kqanxbzIntMLVoiyk+GcViMza4aJSAFUJ2djYDBgwgISEBgOHDh7N06VIeeOABAAoKCqipqeHKK68E4LbbbuMf//hHRBeCrmvg86J7XejuanR3FXpNBXrVMfSqY2iVpWjlR9DKj4DPXfsiWzTmVl2x9LmFVv2G8IcnVxm7ECJCKBxUW9LiZ2OYmeUgTqmmo6WIDtZihrlr8G1diVf11j3XHJsIsU6UmJaYYpqjNGuByRGPEh2HEhVbu9nJGiUfYAwWsEIoKirC6XTW3U5KSmLTpk3nfNzpdFJYWBioOHU01wmOfvRP3JUnLn1iPz15R9fROe0+XaP2pn7y65P/axq6roKmomgqaF50nxfd6z73vBQT5tgW2BMvw9qhO6aEy4i5/AqszrbySyQMd0J3sMF7ORu8l/Ob345hzMOLcJpOcJnlGEmmE9zePQ61/Ci+w1vxVB4DznbimwK2aBRbdO3RTFY7isUOFlvtGobZhmK2gMkCZkvtGfYmM8diHbhdPlDMKCYFMIHJBIoCiql2ukrtP+XU13Dy/9O/5idfn3ZfvZjnuP9cz2+MC3qpghZ/7cXPq6Ep6wE6JfGNN97A5XIxefJkAObPn8/mzZuZPn06AN9//z0vvvgiH3zwAQAHDhzgvvvuY+nSwJ45q1Ydp/iz19Fqqi/shef8gdV/QymnfV33hlNObrpRat+kirn2jayYzOiKGbO19g2vWO21h4Ra7FgcsZiimmGOjsUc2wLVFovd/uNopG6vit1af8iJi73Pn9MK9PQjYVqhlPViXqerPtSqclzlpZg9lahVx1FrKlGrK9E91WgeF5rbhe6tQfO40X0edK+79kOTevKfz4eu+UBVOXu5hLcWQ39NQn//j1IbsDWE5ORkcnNz624XFRWRlJRU7/GSkpK628XFxfUeb4zS0ko07ULfDCZa3TmV4uLguLasQu3b+dRSqCf/957+pBqgpubkFxfH6YwNmmVuKrLMwcwO0a0hGoiv/4gCXMjoWi0THRQXnQBdBU2rWxvXT62Vo9fef+q37PQ1dzhtTV8/rVt+/LtSf63/XCkuoZQu9KWKQnyXrhf1czaZFBITz31UWMAKYeDAgbzyyiuUlZURHR3NsmXLePrpp+sev+yyy7Db7axfv56+ffuyePFiBg8eHKg4QogwpZjMtZuSfvLn7BI24AQ95ZybrS5NwDY+JycnM3nyZCZMmMDYsWO55ZZbSEtLY+LEiWzevBmAmTNn8txzzzFixAhcLhcTJkwIVBwhhBDnEbB9CE3h4jYZhdJqtf/IMkcGWebIcLHLfL5NRnJ4ihBCCEAKQQghxElSCEIIIQApBCGEECeF9EhyJtPFH3p1Ka8NVbLMkUGWOTJczDKf7zUhfZSREEII/5FNRkIIIQApBCGEECdJIQghhACkEIQQQpwkhSCEEAKQQhBCCHGSFIIQQghACkEIIcRJUghCCCGACC6EoqIi/vd//5exY8cybtw48vPzjY7UZLZt20bPnj2NjtEk1q9fz+23386YMWO4++67KSgoMDpSwGRlZTFy5EiGDRvGnDlzjI4TcK+++iqjRo1i1KhRvPDCC0bHaVJ//etfmTp1qt+nG7GF8Oijj5Kens7ixYsZM2YMM2fONDpSk3C5XEyfPh2v13v+J4eBKVOm8Mwzz5CZmUlGRgYzZswwOlJAFBYWMmvWLObOnUtmZibz5s0jLy/P6FgBk52dzTfffMOiRYtYvHgxW7duZfny5UbHahI5OTksWrQoINOOyEIoKytjx44djBs3DoDbb7+dhx56yOBUTeP555/n17/+tdExmoTH4+EPf/gDqampAHTr1o0jR44YnCowsrOzGTBgAAkJCTgcDoYPH87SpUuNjhUwTqeTqVOnYrPZsFqtdOrUicOHDxsdK+DKy8uZNWsW999/f0CmH5GFcOjQIVq3bs2zzz7L6NGjmTRpElar1ehYAbdy5Upqamq4+eabjY7SJGw2G2PGjAFA0zReffVVhg4danCqwCgqKsLpdNbdTkpKorCw0MBEgdWlSxeuvPJKAPbv389nn33GkCFDDE4VeE899RSTJ08mLi4uINMP6eGvG+Pzzz/nueeeq3df+/bt2bZtGw8++CBPPPEE8+fPZ+rUqbz33nsGpfSvsy1zx44dqays5N133zUmVICda5nfffddPB4PU6dOxefzcd999xmUMLDONmixooT/kNC7d+/mvvvu47HHHuPyyy83Ok5AzZ8/n5SUFK699loWLlwYkHlE5PDXBw8e5NZbb2X9+vVA7Xb1AQMGsHHjRoOTBc78+fN58803adasGQA7duwgNTWVOXPmEBNz7otuh7qqqip+97vfkZCQwMyZM7HZbEZHCohFixaRm5vLM888A8A///lPdF3ngQceMDhZ4Kxfv55Jkybx+OOPM2rUKKPjBNw999xDcXExZrOZ48ePU11dzdixY3n88cf9NxM9Qo0YMUJfvXq1ruu6vmTJEv0Xv/iFwYmaVteuXY2O0CR+97vf6U8++aSuaZrRUQLq6NGjenp6ul5aWqpXV1fro0eP1jdu3Gh0rIA5fPiw3r9/fz07O9voKIZYsGCB/thjj/l9umG/yehcXn31Vf785z/z4osvEhMTw/PPP290JOFn27ZtY+XKlXTu3JmxY8cCtdvWZ8+ebXAy/0tOTmby5MlMmDABr9fLHXfcQVpamtGxAubtt9/G7XbX+70dN24cv/jFLwxMFfoicpOREEKIM0XkUUZCCCHOJIUghBACkEIQQghxkhSCEEIIQApBCCHESRF72KkQABs2bOCll16ivLwcXddp1aoVjz32GF26dAnofLt160ZOTg4tWrQI6HyEuBBSCCJieTwe7rvvPt555x169OgBQGZmJhMnTmTlypWYzWaDEwrRtKQQRMRyuVxUVFRQXV1dd9/o0aOJiYlBVVWee+45Nm7cSFVVFbquM2PGDPr27cvUqVOx2+1s3ryZkpISRowYQYsWLfjyyy8pLi5mxowZXHvttUydOhVFUdizZw9lZWUMGjSIJ5988oyBFOfPn88HH3yApmkkJCTwpz/9iU6dOpGbm8vzzz+PpmkA3HfffQwfPrxJv0cissg+BBGx4uPjmTJlCvfeey833ngjU6ZMYcGCBQwcOJCtW7dSVFTEvHnz+Oyzz7j11lvrneG8fft25s2bx4IFC3j33XdxOBx8+OGHTJgwod7zduzYwb///W8+++wz9uzZw7x58+plWLt2LYsXL2bOnDksXryYe++9lwcffBCAV155hXvuuYeFCxfy7LPP8t133zXNN0ZELFlDEBHtnnvu4c4772TdunWsW7eO2bNnM3v2bD7++GMeeughPvzwQw4dOsSaNWvqBgYESE9Px2q14nQ6cTgcXH/99QC0a9eO8vLyuufdeuutda8bM2YMK1eu5H/+53/qHl+9ejUHDhyouzYHwPHjxykvL2fEiBFMnz6dVatWMXDgQP7f//t/gf52iAgnawgiYq1fv55//etfxMTEkJ6ezqOPPsqSJUswmUysWLGibqjsG2+88Ywxcn46aqrFcvbPVqfvh9B1HZOp/q+cpmmMGTOGzMxMMjMzWbRoEQsWLCA+Pp5x48bxySefMGjQIL755htGjx5NRUWFPxZdiLOSQhARq0WLFrz++uvk5ubW3VdcXIzL5WLJkiWkp6czfvx4evXqxYoVK1BV9YLn8fnnn+PxeHC73SxatIj09PR6jw8aNIglS5ZQVFQEwAcffMDdd98N1A7Wtn37dm677TaefvppTpw4wfHjxy9hiYVomGwyEhGrQ4cO/POf/2TWrFkcPXoUu91ObGws06dP57LLLuORRx4hIyMDs9lMv379WLZsWd0O3saKiopi/PjxnDhxguHDh3P77bfXe/z6669n4sSJ/OY3v0FRFGJiYnj11VdRFIVHHnmEZ599lr/97W+YTCYeeOAB2rRp489vgRD1yGinQgTI1KlT6dKlC7/97W+NjiJEo8gmIyGEEICsIQghhDhJ1hCEEEIAUghCCCFOkkIQQggBSCEIIYQ4SQpBCCEEIIUghBDipP8PdNlDOjKZCVkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.ylabel(\"Density\")\n",
    "plt.xlabel(\"Samples\")\n",
    "plt.title(\"Tracy-Widom GUE\")\n",
    "plt.hist(GUE_samples, bins=50, density=True)\n",
    "plt.plot(sol.t, f2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Longest Increasing Subsequence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Tracy Widom distribution suprisingly arises when studying sequences, as explained in [On the Distribution of the Length of the Longest Increasing Subsequence of Random Permutations](https://arxiv.org/abs/math/9810105) by Jinho Baik, Percy Deift, and Kurt Johansson. Let $S_N$ be the group of permutations of $1,2,3,...,N$. Define $\\ell_N(\\pi)$ to be the longest monotone increasing subsequence $\\pi \\in S_n$. Then\n",
    "\n",
    "$$ \\frac{\\ell_N - 2\\sqrt{N}}{N^{1/6}} \\longrightarrow TW_2$$ \n",
    "\n",
    "as $N \\longrightarrow \\infty$. \n",
    "\n",
    "Below we sample from $S_N$, measure their lengths and plot the normalized distribution where $N=5000$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [],
   "source": [
    "def longest_increasing_subsequence(X):\n",
    "    # https://en.wikipedia.org/wiki/Longest_increasing_subsequence#Efficient_algorithms\n",
    "    \"\"\"Returns the Longest Increasing Subsequence in the Given List/Array\"\"\"\n",
    "    N = len(X)\n",
    "    P = [0] * N\n",
    "    M = [0] * (N+1)\n",
    "    L = 0\n",
    "    for i in range(N):\n",
    "        lo = 1\n",
    "        hi = L\n",
    "        while lo <= hi:\n",
    "            mid = (lo+hi)//2\n",
    "            if (X[M[mid]] < X[i]):\n",
    "                lo = mid+1\n",
    "            else:\n",
    "                hi = mid-1\n",
    "        newL = lo\n",
    "        P[i] = M[newL-1]\n",
    "        M[newL] = i\n",
    "        if (newL > L):\n",
    "            L = newL\n",
    "    return L\n",
    "\n",
    "def sample_random_perm_lens(N, sample_size):\n",
    "    lengths = []\n",
    "    numbers = np.arange(1, N+1)\n",
    "    for _ in range(sample_size):\n",
    "        random.shuffle(numbers)\n",
    "        lengths.append(longest_increasing_subsequence(numbers))\n",
    "    return lengths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 5000\n",
    "sample_size = 10000\n",
    "perm_lens = sample_random_perm_lens(N, sample_size)\n",
    "χ = (perm_lens - 2*np.sqrt(N))/N**(1/6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD7CAYAAAB+B7/XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3Qb9YEv8O/M6GX5HUeS7Tjvh0NITIC0MVkaSknWBeIkLN0lQE+gXULoXhbqy3LJoZyWhbbQB5tzgBZKbu9l25vsttmGBi/FhPIq4ECJIU/HsZ13HFtybCfxQ495/O4fxgZjx69IGmn0/ZzDSeQZSd+JpS+jn2Z+IwkhBIiIyDJkswMQEVF0sdiJiCyGxU5EZDEsdiIii2GxExFZDIudiMhiWOxERBZjMzsAAHR0dMMwxn44fV5eBtraumKQKHFxm1MDtzk1jHebZVlCbm76BZcnRLEbhhhXsffdN9Vwm1MDtzk1xGKbORRDRGQxLHYiIothsRMRWQyLnYjIYljsREQWw2InIrIYFjsRkcUkxHHsRJ+XmZUGl3PgS9Pjyez/eyisofN8MN6xiJIGi50SjstpQ/kD2y+4vPKpVeiMYx6iZMNiJ0sZam+/D/f0KVWw2MlShtvb554+pQp+eUpEZDHcY6ekE1H1AV+mEtFALHZKOg67MuxwC1Gq41AMEZHFsNiJiCyGQzGUNJyIYLbdD7WjBYAAIJkdiSghsdgp4dmg42rXISx37UOGHMbJX76FR7PdqNcK8HLPFegSaWZHJEooHIqhhFaktOGR7Jdwk3sXTuoT8FzndcgrW4ejmgeXO47h3swdyJR40hHR57HYKWG5pTC+nfE2AODZ88vxfOdy1KmTkL3o6/j37mvwq87rMEHpxr1ZO5Al9ZgbliiBsNgpIUkQ+Gb6e8iWg/g/XV9Fg1YwaJ1GLR+/6rwOuXJvuadLIROSEiUeFjslpGWu/bjU0YSXehbhhD7xgusd1nx4vvM6TJQ7UZ72cRwTEiUuFjslnNDJOtyQths14Wl4L1w84vpHNB/eCV2Cq1yNCDU1xCEhUWLjUTGUcNrf3oLzwoX/7L4Koz2ksSpYgiudR+F8bRMkLIEY4n7DTUXAmR/JSljslFC05kMInTiAN4NfQgT2Ud8vDAe291yJtc3vYbHTiw/CswetM9JUBJz5kayCQzGUUCIfvwwlPRvVQxTzSGoi0+GafAnK0z6GWwrHIB1RcmCxU8LQ/Y3Qmw4ge/FKqOP6MCkhr+wuZMhhXO08FPV8RMmCxU4JI/zxy5CcGci6smzcj+H0TUNtpBBXuw5BgR7FdETJg8VOCUE/cwz6yb2wL/hbyI6LmyLgnfAlyJaDuNxxPErpiJILi50Sglr7JmBzwjF/2UU/1iG1EC16Nq5xHUTvZGFEqYXFTqYTWgTq4Y9gm7EIksN98Y8HCX8JzcUUWxum21qjkJAoubDYyXTa8U8ANQj77L+J2mN+FJ6BHsOBa5wHo/aYRMmCxU6mUxuqIaVPgFI4N2qPGYEdO8OzUeI4gVy5K2qPS5QMWOxkKqPnHPST+2CffRUkKbovx3fDxZAg8GXH4ag+LlGiG9U7qbKyEjfccAOWL1+OzZs3X3C9t99+G1/72teiFo6sKzMrDR5PJlz+3YAw4P3ycng8mRc85X88OowMHNZ8uNJ5FPwSlVLJiGeB+P1+bNy4Edu2bYPD4cCaNWuwePFizJo1a8B6Z86cwU9+8pOYBSVrcTltKH9gO/4l678hkIf7n/wYQO/sjJVPrYra89REpmNN+geYrLRF7TGJEt2Ie+zV1dUoLS1FTk4O3G43ysrKUFVVNWi9Rx55BPfee29MQpI1FSgdmGxrx0fhGTF7jt2RqdCEjEXOozF7DqJEM+IeeyAQgMfj6b/t9Xqxd+/eAev85je/wbx583DZZZeNK0ReXsa47gcgqh/dk4VVtvlyxzHoQsLHkWkxe46gcKJWnYTLHccgjOHPRE20f9dEyxMP3OboGLHYhRg8NilJn02JWl9fjx07duDFF19ES0vLuEK0tXXBMMY+BurxZKK1NbXm5LPKNns8mVhgP4kjmjfmF6PeFZ6BEsdJBI/tH3a9RPp3tcrveSy4zaMny9KwO8QjDsX4fD6cOXOm/3YgEIDX6+2/XVVVhdbWVtx88824++67EQgEcNttt405KKUW9VwAhbaz2K9Ojvlz1aqTEDTs6DrwbsyfiygRjFjsS5Yswc6dO9He3o5gMIgdO3Zg6dKl/cvvu+8+vPbaa9i+fTteeOEFeL1ebNmyJaahKfn11O8CAByIFMX8uVTYsFedgu66D2CHFvPnIzLbqPbYKyoqsHbtWqxevRorVqxASUkJ1q1bh3379sUjI1lQT+Mu+PUstBpZcXm+XeHpEJEg5tmb4vJ8RGYa1aTX5eXlKC8vH/CzTZs2DVqvqKgIb775ZnSSkWWJSBDB4wdwIDInbs/ZoOVDdmehJHwCe9SpcXteIjPwzFOKO+3UfkDXsF+N/TBMHwEZ7llXYp69CTKMuD0vkRlY7BR32ondkF0ZOKp5R145itLnfAluOYJZNn9cn5co3ljsFFfCMKCf2Av3zMthxPnllzb9MkSEgvmOk3F9XqJ4Y7FTXBmtRyBCnXDPXhT355YdLhxSC7DAfhKcO4asjMVOcaUd3w1IMtJmXm7K8++PTMYEpRuTlA5Tnp8oHljsFFdaUy1k7wwornRTnv+AWgRD4NO9diJrYrFT3IhwN4wzR2GbNM+0DJ0iDcc0D8fZydJY7BQ3WvMhQAgoheYVOwDsUydjsq2dV1Yiy2KxU9zoTbWA4oDim2lqjv2R3vlp5ttPmZqDKFZY7BQ3+ulaKAVzICl2U3MEjGy06pm4hNMLkEWx2CkujJ6zMDpOmz4M0+egWohZdj9sGH6OdqJkxGKnuNCbagHA1C9OP++gOglOScMMnoVKFsRip7jQmmoBZzrkvClmRwEANKo+aELGJfbTZkchirpRze5IdDGEENCbamErmAtJTox9iQjsOKz5cIn9NLYHgYiqX/ASZaGwhs7zwTgnJBo/FjvFnDjvh+huh7LwRrOjDHBQLcRqdw1y5G447ArKH9g+5HqVT61Cal2wjZJdYuw+kaVpCTa+3uegOgkAMJfDMWQxLHaKOf10HSR3DqTsfLOjDNCiZ6NDd/OwR7IcFjvFlBACeks9lIK5kCTJ7DhfIOGgWohiWzOEzmuhknWw2CmmRGcrRM9ZKAXxuwzeWNSpk5AmqwifbjA7ClHUsNgppvTmQwAAJT8xi71eK4AuJPQc/sTsKERRw2KnmNJb6nuPX88tNDvKkILCgRPaRASP7jU7ClHUsNgpprTmetjy50CSEveldkgrQLj5MNKkiNlRiKIicd9tlPSMnrMQ5/0JOwzTp14tAISBWbYWs6MQRQWLnWJGb6kHgIT94rTPMW0iJLsLxfZms6MQRQWLnWJGb64HbA7IE6eaHWVYOhS4pszDHBY7WQSLnWJGbzkExTcLkpz4M1ekTS+BTznPqyqRJbDYKSZEuBtG26mEH1/v455eAgCYw3F2sgAWO8WE7m8EIKAUFJsdZVTsnik4b7g4HEOWwGKnmNBb6gFZgeKdYXaUUZEkCfVqAYrtzZAgzI5DdFFY7BQTeksD5LypkGxOs6OM2iG1AJlyCAXKWbOjEF0UFjtFndA16K1HoeTPNjvKmNRrBQCAYk7jS0mOxU5RZ7SdAHQVim+W2VHG5KyRDr+ehdn8ApWS3KiOQ6usrMRzzz0HVVVx55134vbbbx+w/PXXX8fTTz8NwzCwYMECPPbYY3A4HDEJTIlP9/fOlJg9az7S84a+3FyialDzsch5BDIMGNzvoSQ1YrH7/X5s3LgR27Ztg8PhwJo1a7B48WLMmtW7N9bT04PHHnsML730EiZOnIiKigq89NJLuOWWW2IenhKT7m+ElDkR6XneYS83l4jq1QJc7arHFNsZHNO8ZschGpcRd0mqq6tRWlqKnJwcuN1ulJWVoaqqqn+52+3Gm2++iYkTJ6KnpwdtbW3IysqKaWhKXL0X1mhIumGYPg2aDwCPZ6fkNmKxBwIBeDye/tterxd+v3/AOna7He+88w6uvfZadHR04Oqrr45+UkoKoru998IaSVrsPcKFU1ouZttZ7JS8RhyKEWLwMb1DXeLsmmuuwYcffoh/+7d/w6OPPoqnnnpq1CHy8jJGve4XeTzJNYYbDYm8zV2BPegGMHHuZWZHGbcGtQBXu+pghwb107eIGf/mifx7jhVuc3SMWOw+nw+7du3qvx0IBOD1fjb2ePbsWezfv79/L728vBwVFRVjCtHW1gXDGPtJIR5PJlpbO8d8v2SW6NscatgH2Jw4J01Aso5Q12v5uFaqxXRba/8hkPH+N0/033MscJtHT5alYXeIRxyKWbJkCXbu3In29nYEg0Hs2LEDS5cu7V8uhMCDDz6I06d7j/199dVXccUVV4w5KFmD7m+E4p0BSVbMjjJuh1UfdCFxOIaS1ojF7vP5UFFRgbVr12L16tVYsWIFSkpKsG7dOuzbtw+5ubl4/PHHsX79eqxcuRLHjh3Dgw8+GI/slGCEGobRdiJpx9f7hGHHCW0i542hpDWq49jLy8tRXl4+4GebNm3q//uyZcuwbNmy6CajpKO3HgGEAcWXXGecDqVey8dy1364eLk8SkI8A4OipndGR0DxzTQ5ycWrVwsgSwIzbf6RVyZKMCx2ihrd3wg5pxCSM93sKBftmOZBRCiYw3F2SkIsdooKIQQM/+GkH1/vo0HBUc3DeWMoKbHYKSrEOT9EuAuyBYZh+jSoBZhk64Defc7sKERjwmKnqNADfePr1thjB4B6NR8AEDxxwOQkRGPDYqeo0P2NgCMNck6B2VGi5qSeh5CwI3h0n9lRiMaExU5RofsPQ/HOhCRZ5yVlQEaj6kPoOIudkot13oVkGhEJwug4BcVrnfH1Pg1qPtT2ZhhdbWZHIRq1UZ2gRDQcvfUoIISlxtf71Gu94+zuzmPInD5t0PJQWEPn+WB8QxGNgMVOF63/xCTvDJOTRF+zngs5LROv/Ncr2Pzvg4+OqXxqFVJr2ipKBhyKoYumBw5DzrXGiUlfJCAhber8T49nH/sMpERmYLHTRek/MclrvWGYPmnTFiBX6YFH5r45JQcWO10UK56Y9EWuaQsAgNP4UtJgsdNF6T8xycJ77PYJBThruDHHxml8KTmw2Omi9J+YlGudE5O+SJIk1Kv5mG1vgcRxdkoCLHa6KFY8MWkoh9QCZMhhFCrtZkchGpG1340UU/0nJlnw+PUvqld7P5FwGl9KBix2Gjcrn5j0ReeFGy16Noo5zk5JgMVO46b7GwBIljwxaSiH1ALMtPuhQDc7CtGwWOw0brr/MOTcSZAcbrOjxEW9mg+HpGOardXsKETDYrHTuAhhQPc3WuL6pqPVqOXDEBKK7RyOocTGYqdxMc62AJGelBhf7xMSDpzQ8zCHl8ujBMdip3Ex/Na7YtJoHFILMMV2Bk5EzI5CdEEsdhoX3d8IONMhZeebHSWu6tUCKJLALLvf7ChEF8Rip3HRA42fnpgkmR0lro5qHkSEwnF2SmgsdhozEe6G0XE65YZhAECHgsOqj8VOCY3FTmOmBw4DSL3x9T51WiHylXPIkbvNjkI0JBY7jZnuPwxIqXNi0hcd+nR6gWLbaZOTEA2NxU5jpvsbIU8ogmR3mR3FFM16Ds4ZaZjL4RhKUCx2GhNh6NADh6H4ZpsdxUQS6tRCFNtPQxicXoASD4udxsRoPwWoISj5c8yOYqo6tRDpcgThlqNmRyEaxGZ2AEoueks9AEDJT+U99s/G2bsaP4Fn6d8PuU4orKHzfDCesYgAsNhpjPSWBkjpEyBn5AEAMrPS4HKm3suoW7hwUpuA2cf3ovwBx5DrVD61Crz8NZlhVO/IyspKPPfcc1BVFXfeeSduv/32Acv//Oc/45lnnoEQAkVFRXjiiSeQnZ0dk8BkHiEEdH/DgGEYl9OG8ge2D7l+5VOr4hXNFHVqISafOggnLkMYdrPjEPUbcYzd7/dj48aN2LJlC7Zv347f/e53aGxs7F/e1dWFRx99FC+88AJefvllFBcX45lnnolpaDKH6GqD6O5I8S9OP1OnFgKGjlm8qhIlmBGLvbq6GqWlpcjJyYHb7UZZWRmqqqr6l6uqikcffRQ+nw8AUFxcjOZmHgZmRb0X1uD4ep+jmgeS3Ym5dh7PTollxGIPBALweDz9t71eL/z+zyZAys3NxbJlywAAoVAIL7zwQv9tsha9pQGwuyBPmGx2lISgQ0Ha1PmYZ28CIMyOQ9RvxDF2IQa/YIea+KmzsxP/9E//hLlz5+Kmm24aU4i8vIwxrf95Hk/muO+brMza5lNnDiNtcjG8Pn5/0sc96wpMbKyBR+5Eq5E1aPnF/K742k4NsdjmEYvd5/Nh165d/bcDgQC8Xu+AdQKBAP7xH/8RpaWlePjhh8ccoq2tC4Yx9j0ejycTra2pddyBWdsswt2IBE7AMfmKAc+fim/Ez0ubeQUAYJ79FN4Jzxu0fLy/K762U8N4t1mWpWF3iEccilmyZAl27tyJ9vZ2BINB7NixA0uXLu1frus67rnnHlx//fX43ve+l3LTuKaK3om/BMfXv8Ce40Wzlo1LHU1mRyHqN6o99oqKCqxduxaqquIb3/gGSkpKsG7dOtx3331oaWlBbW0tdF3Ha6+9BgCYP38+fvSjH8U8PMWP3tIASHLKTvw1nFp1Eq5x1cEBFREe9kgJYFTHsZeXl6O8vHzAzzZt2gQAWLBgAerq6qKfjBKK3tIAOW9Kyk78NZxatQjXpdVijr0Z+9UpZsch4lwxNDKhq70TfxUUmx0lIR3RvAgJ+6dHxxCZj8VOI9IDRwBdZbFfgAEZdWoBD3ukhMFipxHpzYcASLCl+IyOw6mNTEKu0oNC5azZUYhY7DQyvflQ74U1XOM/38DqDqqTAPQe9khkNhY7DUsYWu/EXxyGGdZ54cZJbQIuZbFTAmCx07CM1mOAFmGxj8J+tQjTbK3IkDgHO5mLxU7D0k73HsrKYh/ZvsgUyBK4106mY7HTsPTmOsi5kyCnDZ4HhQZq0nPRpqdjgeOk2VEoxbHY6YKEoUP3N3JvfdQk7Fcno9jeDAdUs8NQCmOx0wUZZ473Xri6YK7ZUZLGvshkOCQdc+28JgGZh8VOF6Q3942v8/j10Tqs+dBtODCfwzFkIhY7XZDWfAhydj5kd47ZUZKGARkH1CLMt5+EMHSz41CKYrHTkIShQW8+BKWQwzBjtT8yGelyBKGTB82OQimKxU5D0gNHe8fXJ11qdpSkc1AthCpkdNd/ZHYUSlEsdhqSfmo/AAm2SYOvCkTDi8COQ2oheuo+GPLSkkSxxmKnIelNtZA90yA5082OkpR2R6ZCO38GRuCw2VEoBbHYaRARCUIPHIaNwzDjtk+dDCg2qIf/anYUSkGjuoISpRa9uQ4QBpQiFvt4hYQDrukLoR7fhYkr10GSBu5DhcIaOs9zThmKDRY7DaKdOgDYHFB8swAAmVlpcDn5UhmrrEuvRqBxF7674X/jiOYbsKzyqVUY+7XpiUaH71YaRG86AKWgGJLSe2Fml9OG8ge2D7lu5VOr4hktqbhnL0JEKLjccWxQsRPFEsfYaQCjqx3G2WYeDRMFsjMNteokLHQchwTD7DiUQljsNIDedAAAoEyab3ISa/gkPA1ZcggzbQGzo1AKYbHTAFpTLaS0LMgTisyOYgkH1CKEhQ1XOI6ZHYVSCIud+glh9I6vT5oHSZLMjmMJKmzYHynCQscxKODcMRQfLHbqZwSOQATPwzblMrOjWMpHkZlIlyO8shLFDYud+mkn9gCSDNvkErOjWMohtQDnjDR82cmzUCk+WOzUTzu+G0r+bE4jEGUGZHwUnoF59iZk8kLXFAcsdgIAGF1tMNpPwjZlodlRLOmv4ZlQJIErHUfNjkIpgMVOAHr31gHANpXFHgt+IwfHtIlY7GwEwBkfKbZY7AQA0E7shpTlg5Sdb3YUy/preCYKbWcxSWk3OwpZHIudINQQ9KaDsE1dyMMcY+jjyDRoQsZifolKMcZiJ2hNBwBD42GOMRYUTuyNTMYixxEYatjsOGRhLHaCfnw34EiDUjDH7CiWVx2eg3Q5gu7a982OQhY2qmKvrKzEDTfcgOXLl2Pz5s0XXO+hhx7Ctm3bohaOYk8YBrQTe2ArWgBJ5mSfsdag5aNFz8b5miqzo5CFjVjsfr8fGzduxJYtW7B9+3b87ne/Q2Nj46B17rnnHlRV8cWabPSWQ71nm85YZHaUFCHh3VAxws2HoQeOmB2GLGrEYq+urkZpaSlycnLgdrtRVlY2qMArKytx3XXX4frrr49ZUIoNrfFDwObk+HocfRSeAcnhQuTAG2ZHIYsa8bN3IBCAx+Ppv+31erF3794B69x1110AgJqamnGFyMvLGNf9AMDjyRz3fZNVtLZZ6BqOH69B+pxF8BZMjMpj0sjCcCBzwVfRufsNTFhxFxR31pDr8bWdGmKxzSMWuxCDT6aI9iFxbW1dMIyxn7Th8WSitTW1LjAWzW3WTu6DEeyEXnTlsI+Zim+2WEu7bDnO11Shc/cbyFly04BlobAGl9PG13YKGO82y7I07A7xiMXu8/mwa9eu/tuBQABer3fMQSjxqIf/CtjTYJu8wOwoKSe9YBoaVB/y/vwSHv9D73wyfXi5QbpYI46xL1myBDt37kR7ezuCwSB27NiBpUuXxiMbxZDQVWjHdsE27Yr+a5tSfL0VmocJSjcvwkFRN2Kx+3w+VFRUYO3atVi9ejVWrFiBkpISrFu3Dvv27YtHRooB/dR+IBKEfeZis6OkrFq1CKe1HCxz7YPE+WMoikZ14HJ5eTnKy8sH/GzTpk2D1nvyySejk4piTj38V8CZDqWIF602i4CEP4fmY23Ge7jUfhL71SlmRyKL4JmnKUioIWjHP4F9+pU8Kclkn0Sm4YyegeVp+8FZHylaWOwpSD38IaCGYC/mdyVmMyDjzdClmGY7g9m2FrPjkEVwdy0FqQffhpxbBNk7s/9nmVlpcDn5cjDDh+FZKEvbi+Vp+9DQWWB2HLIAvpNTjH7mOIzWo3AuuX3A+Qgupw3lD2wf8j48/C62NCh4KzQPq901mMW9dooCDsWkGLXuHUCxwz57idlR6HPeCxWjQ3djpbtmyJMCicaCxZ5ChBqG2lAN24wv8YLVCUaFDX8KLsRUWxu6D1abHYeSHIs9hWh9X5pe8lWzo9AQPorMwGktB+1vb4HQVbPjUBJjsaeQSN3bkHMKofhmmx2FhiAgozJ4BbSOFpz/+HWz41ASY7GnCL2lAUbgCOyXfJXXNU1gteokuKZeio73tkKEu82OQ0mKxZ4iwrv/G5IzA/a5PHY9sUnIW/YtGMEuhD/8vdlhKEmx2FOA3nYS+ok9sC9YDsnuMjsOjcCZPx3ZX14Bte4daM2HzI5DSYjFngIiu18B7C44Ll1mdhQahYiqI3fpLbBle6Dt/C0m5rrg8WTC48lEZlaa2fEoCfAEJYszzvmhHfkQ9gVf5yGOScJhV1D+wH9jnn0B1me+iWd/8CR2hEoA9J4sllqXoqDx4B67xUX2vArIChwlZWZHoTGqVYvwSXgqytL2olDpMDsOJREWu4UZ51uh1r8H+5yvQHbnmB2HxuG/er6MbuHEHRl/gQM8tp1Gh8VuYeEP/rN3b/2KlWZHoXHqEmn4f11Xwyufw83pH5kdh5IEi92itKZaaMdq4Lh8BeT0XLPj0EWo1wrwemgBSp2N6Nr/rtlxKAmw2C1IGDrC1ZshZXrgWMCxdSuoCl6GI6oHra/+Cnr7KbPjUIJjsVuQWvsWjI4mOEvXQLI5zI5DUWBAxr93L4XscCFYtRFGz1mzI1ECY7FbjNFzFuGal6BMmgfbtCvMjkNRdNZIR/4tD0OEuhCs2gihhsyORAmKx7FbiBAGQm9tAjQVzr/55qA5YXiVpOQn5U1F/s0PoOX3T0L/ywvI//uHICm9v9NQWEPn+aDJCSkR8F1uIZE9VdCbDsD5lTuh5BQOWs6rJCU/h13Bzc+dwt84v4R/OPwhtj/2L3ixayl0KDx5ifpxKMYi9MARRD76A2zTF8E+9xqz41CMvR8uxn91fwkljpO4K+Mt2KGZHYkSCPfYLUCEuhB883lI7my4ln6L0/KmiHfDl0CDgn9wf4C7M9+EEfq62ZEoQbDYk5xQw+ip2gjR1Y60FQ9xPpgUszM8B6pQcFt6NU793w3I//uH4JhYNOS6HINPHSz2JCZ0FcHXn4HRegSu5ffCls8rI6WiXZGZ6DDS8d30D9D4/L/gN11X44A6edB6HINPHRxjT1LC0BF66wXop/bDtfTbsE+70uxIZKLDWj6Kvv1TBPQs3J35Fm52f8i5ZVIY99iTkAh1IfjGL6E31cJZegvsxV/pX8ZDGlOXLduDp8+Xodz9Mb7irMOl9ib8Z/dVqNcKzI5GccYGSDKRM6fQ/ccfQ3S1wfXVu2Cfc/WA5TykMbWpsGFbz5exOzIVt6ZX439kvY5PIlPxSs/lZkejOGKxJwkhDKh1f0HTX38PyDa4VzwEhWPqdAFHNB9+eq4c16Xtx9dctSjJPoHWV3tgzFkOOctjdjyKMRZ7EtDbTiD03m9g+BvhmnIplKu/DTkjz+xYlOBU2FAVXIj3Q8X427S9WLr7DeCT12GbegXs85dDKSjmobEWxWJPYHrgCCJ7q6Ad/QiSMwOur65DwZIynDnTZXY0SiKdIg1/6FmMbz74P9Hy7suIHHwb2rEaSFle2Gcuhm1WKZTcSWbHpCgaVbFXVlbiueeeg6qquPPOO3H77bcPWH7w4EE88sgj6OrqwqJFi/Cv//qvsNn4/4zxMILnoR3/BFr9+9Bb6gF7GuwLyuBcuAKSK4N7WDRuRlouim78Noy/vR3dte+j68C7CO5+BZFPKqHk5EMuvBS2yfOh5M/h+RBJbsT29fv92LhxI7Zt2waHw4E1a9Zg8eLFmDVrVv86Dz74IH74wx9i4cKFePjhh/H73/8et912W0yDW4UId0NvPQrd3wj99MHeMhcCUpYXztJbYZ+7FJJj4P1W6ekAAAhmSURBVJXpeeQLjUfvRbI//8X6QmRKxVjoOI5vztAQbngPau0bAAD7xCI4C+fAmT8dyClEyOmBnJZlTnAasxHbobq6GqWlpcjJ6b1mZllZGaqqqnDvvfcCAJqamhAKhbBw4UIAwN/93d/h6aefTvliF4YBaCGIcDdEuAci1AXRcxZGz1mIrnYY51pgnGuB6GrrvYMkQZ4wGY7LV8I2/UrIEyZfcO+cR75QtHSKNLwbnov/desqrH5gG6bbWjHdFsC002cwLfA+MuQ3P1vZmQ450wM5ywMpfQLk9FxI7hxI7mxIrsze/5zp/bNNknlG/A0EAgF4PJ99i+71erF3794LLvd4PPD7/VGOOZjR3YGWt3+JcHf3OB9BfO6vYojF4rN1Pl0uhNH7dyEAYQBCBwwDwtABQwN0DUJXAS0M6MNMyuRIg5ydDyV/DuScAijemVC8MwbsmXOvnOJNh4JGLR+NWv6nPxHIkoJ4bv1cGO2noHW0QO1ogdpxEtqJPRBaZOgHUhy9r2VHGiS7E5LdBdgckBQHoNgBxd5b/ooNkBVIcu+fHZluhHtUSLIMSH3/SQP/BABJggSp92ef3gY+//fP/dl7Y+ico1lnrMbyMJIMI2dJdJ73iw8txFCt9pnnn38ewWAQFRUVAICtW7di3759eOyxxwAAH3/8MX72s5/hP/7jPwAAx48fx/r161FVVRWTwH20821ofeUXMCIXc7GBz34LQ+4df+GFI0kSIMsAel9okqJAkpXeF+enL1ZJsUFyuCDbnL1/utKhuNIhp2XClpkLJSMX8heGVoYSUXU47MqQy8KqDmeKLku0PMmyLNqPK4SAEepG5NwZINQJvec89J7zMELdMMLdMILdMNQQRCQEIxKCUMMQegSGGoHQVMDQIDQVQtc+3THSL5jbyibecA+yLl8e9ccdcZfQ5/Nh165d/bcDgQC8Xu+A5WfOnOm/3draOmD5aLS1dcEwhv3/yxAcKLj1+2htTazZL/r2842hFuoAzmnARczY4fFk4vzZnnHfPxl5PJkJ93uOtaTZZiUPSM8DhviutW++kgv/r+YzQgh48txoDZz79BOxDgjR+ynZMACIzz4t4/N/AgM/WYsBH8Y/f0MM+JQ+9DoXZawPI0vInDVnXL9nWZaQl5dxweUjFvuSJUvwzDPPoL29HWlpadixYwcef/zx/uWTJk2C0+lETU0NrrzySvzxj3/E0qVLxxyUiFKXJEm9n3a/cI1eqx8DFquj3EacBMzn86GiogJr167F6tWrsWLFCpSUlGDdunXYt28fAODnP/85nnjiCVx//fUIBoNYu3ZtTMISEdHIRhxjj4fxDcUk0cfVKOI2pwZuc2oY7zaPNBTDaXuJiCyGxU5EZDEsdiIii2GxExFZTEKc2ijL4z/k52Lum6y4zamB25waxrPNI90nIY6KISKi6OFQDBGRxbDYiYgshsVORGQxLHYiIothsRMRWQyLnYjIYljsREQWw2InIrIYFjsRkcUkfbEHAgHcfffdWL16NdasWYNTp06ZHSluamtrMX/+fLNjxEVNTQ1uvvlmrFq1CnfccQeamprMjhQTlZWVuOGGG7B8+XJs3rzZ7Dhx8eyzz+LGG2/EjTfeiJ/+9Kdmx4mrn/zkJ9iwYUP0H1gkuTvuuENs2bJFCCHEli1bxP33329yovjo6ekRt9xyi5gzZ47ZUeLi2muvFQcPHhRCCLF161Zxzz33mJwo+lpaWsS1114rOjo6RHd3tygvLxcNDQ1mx4qp999/X9xyyy0iHA6LSCQi1q5dK3bs2GF2rLiorq4WixcvFg899FDUHzup99jb29tRV1eHNWvWAABuvvlmfPe73zU5VXw8+eSTuPPOO82OEReRSAT3338/5s6dCwAoLi5Gc3Ozyamir7q6GqWlpcjJyYHb7UZZWRmqqqrMjhVTHo8HGzZsgMPhgN1ux8yZM3H69GmzY8Xc2bNnsXHjRtxzzz0xefykLvaTJ0+isLAQP/7xj7Fy5Urcd999sNvtZseKuTfeeAOhUAhf//rXzY4SFw6HA6tWrQIAGIaBZ599FsuWLTM5VfQFAgF4PJ7+216vF36/38REsTd79mwsXLgQAHDs2DH86U9/wjXXXGNyqtj7/ve/j4qKCmRlZcXk8RNi2t7RePXVV/HEE08M+NnUqVNRW1uLf/7nf8b3vvc9bN26FRs2bMBvf/tbk1JG11DbPGPGDHR1deHFF180J1SMXWibX3zxRUQiEWzYsAGapmH9+vUmJYwdMcREq7G6in2iaWhowPr16/HQQw9h2rRpZseJqa1bt6KgoABXXXUVtm3bFpPnSOppe0+cOIGbbroJNTU1AIBgMIjS0lLs2bPH5GSxs3XrVvzqV79Ceno6AKCurg5z587F5s2bkZFx4YvbJrvu7m585zvfQU5ODn7+85/D4XCYHSnqXnrpJezatQs/+tGPAAC/+MUvIITAvffea3Ky2KqpqcF9992Hhx9+GDfeeKPZcWLuW9/6FlpbW6EoCs6dO4eenh6sXr0aDz/8cPSeJOqj9nF2/fXXi7ffflsIIcQrr7wibr31VpMTxVeqfHn6ne98RzzyyCPCMAyzo8RM35enbW1toqenR6xcuVLs2bPH7Fgxdfr0abF48WJRXV1tdhRT/OEPf4jJl6dJMxRzIc8++yx+8IMf4Gc/+xkyMjLw5JNPmh2Joqy2thZvvPEGZs2ahdWrVwPoHX/etGmTycmiy+fzoaKiAmvXroWqqvjGN76BkpISs2PF1K9//WuEw+EB79s1a9bg1ltvNTFV8kvqoRgiIhosqY+KISKiwVjsREQWw2InIrIYFjsRkcWw2ImILIbFTkRkMSx2IiKLYbETEVnM/wdQJcy7klsvjAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(χ, bins=28, density=True)\n",
    "plt.plot(sol.t, f2, label=\"β=2\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We expect the bias to be corrected as we increase $N$, but that may be very computationally expensive."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}