{
 "metadata": {
  "name": "",
  "signature": "sha256:f8ea0ad2b3856d5d590c8799f3ea5b138e8f3202148c03b13acf15e7ebddfa26"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Linear Regression**\n",
      "\n",
      "Consider a set of $n$ data points $(x_1,y_1)$, ... , $(x_n,y_n)$. Linear regression uses the data points to predict a linear relationship between $x$ and $y$ (and for the single variable case is equivalent to least-squares fit).\n",
      "\n",
      "Assume that the relation between $x$ and $y$ is of the form $y=ax+b$.\n",
      "The values of $a$ and $b$ are determined by considering how much that assumed relation misses for the data points, i.e., $(ax_i +b)-y_i$ for the $i^{\\rm th}$ data point. To minimize the error, we define a \"cost function\":\n",
      "\n",
      "$$C(a,b) =\\sum_{i=1}^n ((ax_i +b)-y_i)^2 = ((ax_1 +b)-y_1)^2 + \\ldots + ((ax_n +b)-y_n)^2\\ .$$\n",
      "\n",
      "The values of $a$ and $b$ that minimize this function are determined by taking the derivatives and setting equal to zero:\n",
      "\n",
      "$$0={dC\\over db} = \\sum_{i=1}^n 2((ax_i +b)-y_i) = 2(a\\sum_{i=1}^n x_i + nb -\\sum_{i=1}^n y_i )\n",
      "=2n(a\\overline{x}+b-\\overline{y})$$\n",
      "$$0={dC\\over da} = \\sum_{i=1}^n 2((ax_i +b)-y_i)x_i = 2(a\\sum_{i=1}^n x_i^2 + b\\sum_{i=1}^n x_i - \\sum_{i=1}^n y_ix_i)= 2n(a\\overline{x^2}+b\\overline{x}-\\overline{xy})$$\n",
      "\n",
      "having set $\\overline{x}= {1\\over n}\\sum_{i=1}^n x_i$ (the average value of the $x_i$),\n",
      " $\\overline{y}= {1\\over n}\\sum_{i=1}^n y_i$ (the average value of the $y_i$), and similarly\n",
      "using the additional averages $\\overline{xy}={1\\over n}\\sum_{i=1}^n x_iy_i$ and $\\overline{x^2}={1\\over n}\\sum_{i=1}^n x_i^2$.\n",
      "Then  the first equation above becomes \n",
      "\n",
      "$$b = \\overline{y} - a \\overline{x}$$\n",
      "\n",
      "and substituting for $b$ in the second equation gives\n",
      "$0 = a\\overline{x^2} + (\\overline{y} - a \\overline{x})\\overline{x} - \\overline{xy}$ . Solving for $a$ gives the result\n",
      "\n",
      "$$a= { \\overline{xy} - \\overline{x}\\,\\overline{y} \\over \\overline{x^2} - \\overline{x}^2}\n",
      "= {\\overline{(x-\\overline{x})}\\overline{(y-\\overline{y})} \\over\n",
      "\\overline{(x-\\overline{x})^2}}={{\\rm Cov}[x,y]\\over {\\rm Var}[x]}$$\n",
      "\n",
      "The covariance Cov$[x,y]={1\\over n}\\sum_{i=1}^n (x_i - \\overline{x})(y_i - \\overline{y})$ measures the average of how much $x$ and $y$ are both on the same side of their means ($+$ contribution) or on opposite sides ($-$ contribution), so is overall positive when they're correlated, and negative when anti-correlated. The covariance of a random variable with itself is just the ordinary variance: Cov$[x,x]$ = Var$[x]$, measuring the average square of the deviation from the mean.\n",
      "\n",
      "The two formulae above for $a$ and $b$ are readily calculable for any set of data points $(x_1,y_1)$, ... , $(x_n,y_n)$. Let's generate $n=100$ synthetic data points, with $x$ values between 0 to 10 (recall `rand(n)` generates `n` random values between 0 and 1):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n=100\n",
      "x=10*rand(n)\n",
      "y=1.5*x + normal(0,2.5,size=n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The $y$ values in the above are given by multiplying the $x$ values by 1.5 (a linear slope), and then adding noise by pulling 100 independent numbers from a normal distribution (with mean 0 and standard deviation 2.5). The data looks like this:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "figure(figsize=(6,5))\n",
      "plot(x,y,'o');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAE4CAYAAACpJRM9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHLlJREFUeJzt3X+QXWddx/HP18Q2AVpTYunGNNrOHX+gyEDHqVXH7B2G\n7I1GC/xDYezYUeqo4G5Fh4HuJmRDykh1UDaLINrCFEGEQWXKri4bxd3F0TqggTZtGTRDx7ZkU4nS\nUKZhSP36x72bvZu9d/fc8/M557xfM3fYPXt/PDlDP/c5z/k+z2PuLgBAuX1X0Q0AACRHmANABRDm\nAFABhDkAVABhDgAVQJgDQAVECnMz22Nm/2hmD5vZSTMb6xx/oZkdN7OvmNm8me3ItrkAgF4sSp25\nmQ1JGnL3L5rZCyT9m6RXS/oVSV939983s7dKusrd35ZpiwEA60Tqmbv7srt/sfPzM5IelbRb0s2S\n7us87T61Ax4AkLOBx8zN7DpJL5f0r5KucfcznT+dkXRNai0DAEQ2UJh3hlj+StId7v7N7r95e7yG\ntQEAoABboz7RzL5b7SD/c3f/VOfwGTMbcvdlM9sl6akeryPgASAGd7eoz41azWKS7pX0iLu/p+tP\n90u6rfPzbZI+delrOw3i4a7Dhw8X3oZQHpwLzgXnYuPHoKL2zH9G0q2SHjSzE51jd0p6l6RPmNkb\nJD0m6bUDtwAAkFikMHf3f1L/Xvwr02sOACAOZoDmqNlsFt2EYHAuVnEuVnEu4os0aSjRB5h51p8B\nAFVjZvK0b4ACAMJGmANABRDmAFABhDkAVABhDgAVQJgDQAUQ5gBQAYQ5AFQAYQ4AFUCYA0AFEOYA\nUAGEOQBUAGEOABVAmANABRDmAFABhDkAVABhDgAVEHVDZwAondnZJR07Nq9vf3urLr/8gsbGRnTg\nwN6im5UJwhxAJc3OLumOOz6jU6feefHYqVMTklTJQGeYBUAlHTs2vybIJenUqXdqevp4QS3KFmEO\noJK+/e3eAw/nz2/JuSX5IMwBVNLll1/oeXzbtudybkk+CHMAlTQ2NqJGY2LNsUZjXKOj+wpqUbbM\n3bP9ADPP+jMAoJfZ2SVNTx/X+fNbtG3bcxod3Veam59mJne3yM8nzAEgPIOGOcMsAFABhDkAVABh\nDgAVQJgDQAUQ5gBQAYQ5AFQAYQ4AFUCYA0AFEOYAUAGsZw4AOch6owzCHAAylsdGGQyzAEDG8tgo\ngzAHgIzlsVEGwywAainPzZ7z2CiDMAdQO3lv9jw2NqJTpybWfF57o4z9qX0G65kDqJ1W66Dm5+/q\ncfyQ5uaOZvKZg26UMeh65vTMAdROEZs9HziwN9NdjrgBCqB2qrjZM2EOoHaquNkzY+YAain0zZ4z\n2dDZzD4o6YCkp9z9xzvHJiXdLum/O0+7093neryWMAeAAWW1ofOHJF1aQ+OS/tDdX955rAtyAEA+\nIlWzuPvnzOy6Hn+K/K0BACHKc/JQlpKWJo6a2S9L+oKk33X3b6TQJgDIRd6Th7KUpJrl/ZKul/Qy\nSaclvTuVFgFATvJYACsvsXvm7v7Uys9mdo+kT/d77uTk5MWfm82mms1m3I8FgNQUMXmon4WFBS0s\nLMR+fewwN7Nd7n668+trJD3U77ndYQ4AoQhp8tClHd0jR44M9PpIwyxm9jFJ/yzph83scTP7VUl3\nm9mDZvYlScOS3jzQJwNAwao0eYhJQwBqLa3JQ2lXxbDQFoBSCKUkcOUzjx2b1/nzW3Xs2Pya41GE\nUBVDmAPIXQjhl2Zb+lfFHMrt38NCWwByF1JJYBptCaEqhjAHkLsQwm9FGm0JoSqGMAeQuxDCb0Ua\nbQmhKoYxcwC5y2NPzDzbsjIuPj19qKsqZn+u4/+UJgIoREjriYfUlhWZrGeeBGEOAIOjzhwASiLN\nWnvCHAAy1iu0JaVaa88wCwBkqNekpEZjQldeeUYnTtyz7vmt1iHNzR1lmAVAb6FMny+zOOew36Sk\nq656fc/nx621J8yBGghp+nxUoX35xD2H/SYluV/W83jsWnt3z/TR/ggARRoZmXDJ1z1arYNFN62n\nmZlFbzTG17S10Rj3mZnFwtoU9xz2e90NN7zRh4bevObY0NBvX/w3drIzctYyAxSogZCmz0cR0tot\nK+Kew36zQ3/xF39M0tOSDkma7PzvudjtY5gFqIGQps9HEeKXT5xzuDJUtG3b17Vz5y3atWuXdu++\nQqOj+3Xs2LyWl+9d8/zlZcVeaZGeOVADIawdMogQv3wGPYcrY+zz83fp4Yc/oLNnP65nn33+xdml\naX9h0TMHaiCEtUP66XWjM6S1W1YMeg43W+M87S8swhyoiQMH9gYR3t36VYhMTbU0NdUK7stnkHO4\nWc877S8swhxAYTbqvc7NHS08vJPYrOed9tUSYQ6gMCHe6ExLlJ53mldLhDmAwoR4ozMted+nYG0W\nAIXpvW7JuKamih8fLxrrmQMolRA3hggBYQ4AFTBomDNpCAAqgBugAIIS2mqJZUGYA9hQv3DNInTL\nuFRvKAhzoOY2CuV+4fr5z5/URz7yZOqhu9kUePRHmAM1tllPuF+4vve9t+js2Y+vO/72t78pUehW\neRJR1rgBCtTYZuuG9wvXCxe29zz+yCPf1OzsUuz2VHkSUdYIc6DGNusJ9wvXrVuf7fO6H0i0gUTZ\nluoNCcMsQMmkeeNxs55wv/VFbr11WHff/Zs6f/79Xa8al7Rf589/NlZbpLCX6g0dYQ6USNJqj0u/\nCH7qp75vw8WgNgrX+++/XSdOHJK0RdJzkvZL2qtt25Jt7RbiUr1lwAxQoERarYOan7+rx/H2krEb\n6b0OyoRuvXW3Hnjg9MDT6VlXJVuDzgClZw6USJJqj343Ox94YPMvgl4YEgkLYQ6USJJqjyzK/hgS\nCQfVLECJJKn2oOyv2uiZAyWSZGgjxE2SkR5ugAI1wtrh5cF65gBQAVSzACXF0q9IgjAHAsDSr0iK\nYRZAxfeKk0wGQjUxzAIMKIReMUu/IinCHLUXwoYIedSAR736KPoqBfEQ5qi9EHrFg9SAxwnbqFcf\nIVylICZ33/Qh6YOSzkh6qOvYCyUdl/QVSfOSdvR5rQMhGxmZcMkveSz6zp2v9eHhwz4yMuEzM4uZ\nt2NmZtFbrYM+PHzYW62DPT9zZmbRG43xNW1tNMY3bV/vf6N7q3Uw1vOQvU52Rspod4/cM/+QpGlJ\nH+469jZJx939983srZ3f35b86wXI1/pe8ZK2bv0LnT37cS0uto/k0TuNss5J3CGhqFcfIVylIJ5I\nYe7unzOz6y45fLOk4c7P90laEGGOErp0ivzJk4/23N9ys8DMY6w5bthGGZOfnV3SyZOPSpqUdEHS\niKS9656HMCUZM7/G3c90fj4j6ZoU2gMUortX3GxOXuyRd9soMPMaa457o3SzMfmV9q/9EpvoPG+O\n9VtKIJUboO7uZta3mHxycvLiz81mU81mM42PBTIRJzDzqoiJu1hW99XHE088peXlb2j79l06dmy+\nb/uld2rnztdpauqN3PzMwcLCghYWFmK/PkmYnzGzIXdfNrNdkp7q98TuMAdCFycw8xprTrJq4spz\n2j3wD+jsWenkyfYVxPbt3+r5mpe85EcI8pxc2tE9cuTIQK9PEub3S7pN0t2d//1UgvcCghE1MLvH\nyNtjzetlMdacZEOIflcQO3fe0vP5jJWXR6QwN7OPqX2z83vN7HFJb5f0LkmfMLM3SHpM0muzaiSQ\nt80Cc/0Y+ZK2bv0NXbjwJxefE+Ja4f2uIIaGdmjHjo3H1A8d+rAee+wZuV+u669/vo4efR299oBE\nrWZ5fZ8/vTLFtgClsb6Hu1cXLkg7d75OL3nJj2w6/FHULMt+9wOuvfZFGh3d1/NqZHZ2Sbfffp+W\nl4ck3SNJOnFCuv3239E99zCZKBiDFKXHeYhJQ6ig4eHDPSfXDA8f3vS1cSf+pKH3Z9+54We3JxIx\nmShvymjSEIAuSdZSKXItmDg3UPsNzUhMJgoJYQ7EkGQ/zaJnWQ56A7XfF5fEDdKQEOZADElKBPNY\nITFNY2MjevDB+7S8PCFp9ctraOjNGh19TXENwxpsTgHkrNds0UZjXFNT0b4MirBazfItSZfp+utf\noHe845Zg21sFbOgMlMDs7JKmp4939er3EYxYgzAHKoANIsC2caiUOoYaG0QgDsIcuchyd5yqCWEb\nO5QPYY7MxQ3luoZa0aWLKKfvKroBqL7+oXx8w9fVNdTKVrqIMBDmyFyWu+NU0djYiBqNiTXH2hOS\n9hXUIpQBwyzIXFa741RVkglJqC9KE5G5JJNk0qjHrmNFDMqPOnMEqahJMr2/SCY0NdUKItD5okE/\nhDnQpdU6qPn5u3ocP6S5uaOZfvZmQR36Fw2KxaQhoEtRFTFRyjHLVnrJVUTYCHPkpogwKKoiJkpQ\nl6n0sq4TuMqEMEcuigqDoipiogR1mUovy3YVUUeEOXJRVBgUVeYXJaiTftGkdaUT5X3KdBVRV4Q5\nclFkGAy6s04/g4TnalC3JM1L2qrt2x/VTTcNr2mXFO+LJq0rnajvU6ariNoaZMPQOA+xoTN8ZVPg\n8m4IHGcT5sOH/9i3b//1TDZuTut8Rn2fOBtBIxmxoTOK1qsHW/bZnHGGif7lX76mZ5/9k4FeE1Va\nVzpR34dZqeEjzJGqfpftU1MtTU21ShsGccIzy6GltIY9BnmfuMNVlDTmgzBHqjbqwc7NHS3tf8Rx\nwjPLcea0rnSyvmKipDE/hDlSVdWqhzihl2VQpjXskfXwCSWN+SHMEUu/S+eqVj3ECb2sgzKtKp20\n3qeXqn65h4gwx8A2unTO40ZnUWOwcUIvy6Asg6p+uYeIMMfANhsXl7LrjTIGWy5lr2IqE1ZNxMCa\nzUktLk6uOz48PKmFhdXjWfSgi1wFEfEUtfxx2bFqIjIX5dI5qx40Y7DlU/ehprywBygGFmWPyrib\nOG+GMVigN3rmGFiUKo2setCMwQK9EeaIZbNL56x60EwrB3rjBigykWQT5yIw5Ryh4QZoiVUpUMrU\ng6bcEVVAzzwQbO5bHModEaJBe+ZUswQiq+qPoszOLqnVOqhmc1Kt1kHNzi4V3aS+KHdEFTDMEoh+\ngfLAA4+r1TpYqiGXLIctshiKotwRVUCYB6JfoDz99B7Nzx8t1RhuVivlZfUlQbkjqoAwD0SvQJHG\nJbUDJc1lQ9Ps3fZ6rzSGLXq9b1ZfEmW6WQv0Q5gHojtQHnjgcT399B61g3w1UNIYw02zd9vvva68\n8kzP50cdtuj3vtu2fb3n89M4L0w5R9nV8gZoqDfnDhzYq7m5o/rJn7xW0lF1B7mUzhhumjda+72X\n2eWbTveP877Ly9/o+XzGtoEa9szLUFOc5RhumpUb/d7riiuu1jve8YrYwxb93nfXrl3asYOxbaCX\n2oV53ttYxRmfznIMN83KjY3ea6Wtx47N6/z5rTp2bF5StC/Mfu+7e/cVGh3dx9g20Iu7Z/pof0Q4\nhocPu+TrHsPDh1P/rJmZRW80xtd8TqMx7jMzi6l/VrI23RmrTRu9V5J/e5ptBMqqk52Rs7Z2PfM8\na4pD3Mw2zV7/Ru/Vah2M/W+nugQYXOIwN7PHJJ2T9Jyk77j7jUnfM0t51hSHOrMwzcqNfu+V9N9O\ndQkwmDR65i6p6e7/k8J7ZS7PXl+dZxbW+d8OFCGtYZbIi8GEIK9eX51nFob6b6/SypRAt7R65n9v\nZs9J+oC7/1kK71kJdR77DfHfXoayVCCuxEvgmtkudz9tZldLOi5p1N0/1/V3P3z48MXnN5tNNZvN\nRJ8JxMFStwjZwsKCFhYWLv5+5MiRgZbATXU9czM7LOkZd3931zFP8zOAuJrNSS0uTq47Pjw8qYWF\n9ceBIuW605CZPU/SFnf/ppk9X9KIpCNJ3hPISog3ZRnDR1qSjplfI+lvzGzlvT7q7vOJW4XYCIf+\nQrspyxg+0pQozN39q5JellJbkFDUcKhr4Id2UzbESWUor9rNAK2yKOFQ995gSJORQp1UhnKq5RK4\nRchj2d2vfe2ZnseffPKbF3+u2l6jZRbiGD7Ki555DvLqDZ8+fbrP8eWLP9MbDEdoY/goN8I8B3mN\njQ4N7dDZsxOS1m49NzT0PRd/ozcYjtDG8FFuhHkO8uoN7959tR5+eETSIUlb1F77bL+uvXZ1CIXe\nYFhCGsNHuRHmOcirN9wO6s9sGNT0BoFqSnUGaM8PYAZozzHzRmNcU1Pph+js7JKmp493BfW+4IK6\nrqWRwCAGnQFKmOekDCGbh95fbBOammrV8nwA/RDmHfT+eiv6vLDYFRBNrmuzhKruE2P6CeG8UBoJ\nZKOSk4bSnBiTx2SfvIQwYYjSSCAbleyZp9X7C6Enm6YQesWURgLZqGSYp9X7q9pCSCH0iimNBLJR\nyTBPq/cXQk82TaH0ipkoA6SvkmGeVu8vhJ5smvLqFXdXzJw794Sky3TllS+iqgjIUGVLE9OQ52Sf\nqlh7zpYkfUbda8VQUw5EQ515yqo02SePGvO1deQHJVFTDsRBnXnK8hzfzTJs86rMWXufoVr3HICQ\nEeaByDps86rMWXufoVr3HICQVXLSUBllPaEnr8qcsbERNRoTnd9GJE2s+Xu7emZfqp8JgJ55bjYb\nQsk6bPOqzLm0YubcuTMye5OuuOJqasqBDBHmOYgyhJJ12OZZY04dOZA/wjwHUcarsw5bZl4C1Va5\nMC96iddeogyh5BG2UXvMIZ5DABurVJiHujBW1CGUEIYnQj2HADZWqWqWEJZ47WVthUdbqFUdoZ5D\nABsrvGee5iV9qAtjlWm8OtRzCGBjhYZ52pf0IS+MVZbx6pDPYRaKPt9Aatw900f7I3obGZlwydc9\nWq2DfV+zkZmZRW80xte8V6Nxp8/MLMZ6vzTMzCz6yMiEDw8f9pGRiQ3b0rv947m2P8RzmJUQzjfQ\nTyc7I2dtoT3z9iX9kqR5tS8SLkgaiX1JH9pwxqBXHiFshhHaOcxSCOcbSEuhYd5e63rtEqnShM6d\nOxP7PUOoCFkxaFiEMl4d0jnMUijnG0hDwTdAL9PaIJekd8rsTUU0JnWDhsVG49WM7aavbvcHUG2F\nhvmVV76o5/Errrg655ZkY9Cw6DcL9KabrqX2OwOhbKMHpKHQMK96z2jQsOg3Xs3YbjbqdH8A1Vdo\nmFe9ZxQnLHqNV//BH3y253MZ202uLvcHUH2FhnkdekZphMXqFczayp9z55YTtg5AVbAHaAnMzi7p\n9tvv0/LykLpvGA8N/Y7uuefVlfryA9BWuQ2dqeJou+GGN+rEifetO87myEA1VWpDZ1bwW9Wv8odx\ncwBS4KsmsoLfqqpX/gBIJugwZ4beqjItowsgf0EPs9AbXVWHyh8A8QV9A7TXmHmjMa6pKUIMQLVV\nspplevp4V290H0EOoPIqF+ZlQQklgDRVqjSxLCihBFC0xNUsZrbfzL5sZv9hZm9No1FlQwklgKIl\nCnMz2yLpvZL2S/pRSa83sxen0bAyoYQSQNGS9sxvlPSf7v6Yu39H0l9KelXyZpULJZQAipY0zHdL\nerzr9yc6x2qFCT0Aipb0Bmj1y1QiYEIPgKIlDfMnJe3p+n2P2r3zNSYnJy/+3Gw21Ww2NyzlK1OZ\n39q2ut7yllcE21YA4VpYWNDCwkL8N3D32A+1vwxOSbpO7d2ZvyjpxZc8xy81M7Pojca4S37x0WiM\n+8zM4oZ/C02Z2gqgXDrZGTmPE08aMrOfk/QeSVsk3evuv3fJ3/3Sz2i1Dmp+/q5179VqHZK79/1b\nknW7s+jtb/TvYI1xAEnkPmnI3f9O0t8N8po4pXxJyvyymtRDSSKAUBSyBO5GpXxZlPllNamHkkQA\noSgkzDcq5cuizC+rHjQliQBCUcjaLFFK+dIs88uqB01JIoBQFL5qYh5liKyLDqBsSrVqYl6rDdKD\nBlB1hfbMKe0DgN4G7ZkXuqEzpX0AkI5Cw5zSPgBIR6Fhnldp3+zsklqtg2o2J9VqHdTs7FKq7w8A\nRSv0BmjSG5NRKmHY0g1AHRRemhhX73LDCU1NtdaENDdZAZRRqW6AJhF1ij43WQHUQWnDPGpIc5MV\nQB2UNsyjhjTrpwCog0JvgCYxNjaiU6cm1k3RHx3dv+Z5zP4EUAelvQEqtW+CTk8f7wrpfYQ0gEoY\n9AZoqcMcAKqqNtUsAIBVhDkAVABhDgAVQJgDQAUQ5gBQAYQ5AFQAYQ4AFUCYA0AFlHY6/6WirG0O\nAFVViTBnAwoAdVeJ6fxZbEBBTx9AkQadzl+JnnnaG1DQ0wdQNqW+AbqyUfODD36559/jbkARdRcj\nAAhFaXvma3vPS5ImJG28tnlUbDUHoGxKG+Zre88rQx+HdNVV/6Ubb/z+RBtQsNUcgLIpJMzTuLm4\nvve8V9JevfSlk5qbm0zUvqi7GAFAKHIP87RuLmbZe2arOQBlk3tpYlplhL2+FBqNcU1NEboAyi/4\n0sS0bi7SewaAVbmHeZrDIwcO7CW8AUAF1JmPjY2o0ZhYc6x9c3Ff3k0BgMooZDr/7OySpqePdw2P\n7KOHDQBdBh0zr8TaLABQNYOGeamn8wMA2ghzAKgAwhwAKoAwB4AKIMwBoAIIcwCoAMIcACogdpib\n2aSZPWFmJzoP1ocFgIIk6Zm7pD9095d3HnNpNaqqFhYWim5CMDgXqzgXqzgX8SUdZok8Own8H7Ub\n52IV52IV5yK+pGE+amZfMrN7zWxHKi0CAAxswzA3s+Nm9lCPx82S3i/pekkvk3Ra0rtzaC8AoIdU\nFtoys+skfdrdf7zH31hlCwBiyGWnITPb5e6nO7++RtJDSRsDAIgnyU5Dd5vZy9SuavmqpF9Pp0kA\ngEFlvp45ACB7mc0ANbP9ZvZlM/sPM3trVp9TBma2x8z+0cweNrOTZjZWdJuKZGZbOhPNPl10W4pk\nZjvM7JNm9qiZPWJmNxXdpqKY2Z2d/z4eMrO/MLPLi25TXszsg2Z2xswe6jr2wk4BylfMbD5KtWAm\nYW5mWyS9V9J+ST8q6fVm9uIsPqskviPpze7+Y5JukvSmmp+POyQ9ovYQXZ1NSfpbd3+xpJdKerTg\n9hSiU0Dxa5Ju6BRRbJH0uiLblLMPqZ2V3d4m6bi7/5Ckf+j8vqGseuY3SvpPd3/M3b8j6S8lvSqj\nzwqeuy+7+xc7Pz+j9n+031dsq4phZtdK+nlJ96jGk87M7Hsk/ay7f1CS3P2Cuz9dcLOKck7tDs/z\nzGyrpOdJerLYJuXH3T8n6X8vOXyzpPs6P98n6dWbvU9WYb5b0uNdvz/ROVZ7nV7IyyX9a7EtKcwf\nSXqLpP8ruiEFu17Sf5vZh8zs383sz8zseUU3qgju/j9qz1P5L0lfk/QNd//7YltVuGvc/Uzn5zOS\nrtnsBVmFed0vn3sysxdI+qSkOzo99Foxs1+Q9JS7n1CNe+UdWyXdIOl97n6DpG8pwqV0FZlZQ9Jv\nS7pO7SvWF5jZLxXaqIB4u0pl00zNKsyflLSn6/c9avfOa8vMvlvSX0n6iLt/quj2FOSnJd1sZl+V\n9DFJrzCzDxfcpqI8IekJd/985/dPqh3udfQTkv7Z3c+6+wVJf632/1fq7IyZDUntOT2SntrsBVmF\n+Rck/aCZXWdml0m6RdL9GX1W8MzMJN0r6RF3f0/R7SmKu4+7+x53v17tG1yfdfdfLrpdRXD3ZUmP\nm9kPdQ69UtLDBTapSF+WdJOZbe/8t/JKtW+Q19n9km7r/HybpE07gEkmDfXl7hfM7LckfUbtO9P3\nunst79R3/IykWyU9aGYnOsfuZNng2g/HjUr6aKfDc0rSrxTcnkK4+5c6V2hfUPteyr9L+tNiW5Uf\nM/uYpGFJ32tmj0t6u6R3SfqEmb1B0mOSXrvp+zBpCADKj23jAKACCHMAqADCHAAqgDAHgAogzAGg\nAghzAKgAwhwAKoAwB4AK+H+EQ4NWNMrH8QAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c32dd10>"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now calculate the best linear fit slope $a$ and $y$-intercept $b$ using the formulae derived above:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xbar=mean(x)\n",
      "ybar=mean(y)\n",
      "a=(x.dot(y)/n-xbar*ybar)/(x.dot(x)/n-xbar*xbar)\n",
      "b=ybar-a*xbar\n",
      "print a,b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.45955689947 0.191167623632\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The slope is close to the generating value, and together with the above data points looks like this:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "figure(figsize=(6,5))\n",
      "plot(x,y,'o')\n",
      "xr=arange(11)\n",
      "plot(xr,a*xr+b,linewidth=3);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAE4CAYAAACpJRM9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9N/DPSQgJu4ACChY0Siu4gMpOlkchQ0Wt1p87\nFutSf1SBR0WEZFLSJmEXTIJrxRZr0bo9VkmFUDVM2BHZF5XITgLIFsAkQHKeP2YCM5ntzszd5s7n\n/XrlZXLvzNyTiXzn3HO+53yFlBJERBTd4oxuABERRY7BnIjIAhjMiYgsgMGciMgCGMyJiCyAwZyI\nyAIUBXMhxOVCiK+FEFuEEJuFEGNcx9sJIRYLIb4XQpQIIS7StrlEROSLUJJnLoToBKCTlHK9EKIl\ngLUA7gLwewA/SSmnCyFeBNBWSjlB0xYTEZEXRT1zKWWllHK96/tTALYB6AzgTgDzXA+bB2eAJyIi\nnYU8Zi6E6AagN4BVADpKKQ+6Th0E0FG1lhERkWIhBXPXEMvHAMZKKU+6n5PO8RruDUBEZIAmSh8o\nhEiAM5D/Q0r5qevwQSFEJyllpRDiUgCHfDyPAZ6IKAxSSqH0sUqzWQSAuQC2Silfdjv1GYCRru9H\nAvi08XNdDeKXlJg0aZLhbTDLF98Lvhd8LwJ/hUppz3wQgBEANgoh1rmOTQQwFcAHQojHAewCcF/I\nLSAioogpCuZSyqXw34sfol5ziIgoHFwBqqP09HSjm2AafC8u4HtxAd+L8ClaNBTRBYSQWl+DiMhq\nhBCQak+AEhGRuTGYExFZAIM5EZEFMJgTEVkAgzkRkQUwmBMRWQCDORGRBTCYExFZgOJdE4mIKDTF\nxQ4UFpagtrYJEhPPYcyYDAwfnqrJtRjMiYg0UFzswNixi1Benn/+WHl5FgBoEtA5zEJEpIHCwhKP\nQA4A5eX5KCparMn1GMyJiDRQW+t74KOmJl6T6zGYExFpIDHxnM/jSUl1mlyPwZyILKO42AGbzY70\n9BzYbHYUFzsMa8uYMRlITs7yOJacnInRo4dqcj1OgBKRJeg94RhMwzWLirJRUxOPpKQ6jB49TLO2\ncD9zIrIEm82OkpI8H8ezsXBhrgEtigz3MyeimKT3hKPZMJgTkSXoPeFoNgzmRGQJek84mg3HzInI\nMoqLHSgqWuw24TjUkMlPNYQ6Zs5gTkRkQpwAJSKKQQzmREQWwGBORGQBDOZERBbAYE5EZAEM5kRE\nFsBgTkRkAdw1kYhML9RamnrW3jQLBnMiMrVQt7Y121a4euEwCxGZWqi1NPWuvWkWDOZEZGqhbm0b\nq1vhMpgTkamFurVtrG6Fy2BORKYW6ta2sboVLndNJCLTc9/a9uTJw5CyFq1bd/GbqWKFrXC5BS4R\nWZavTJXk5CwUFNh0CdZ6pjyGGsyZmkhEUcN/pkq25sHc7CmPHDMnoqhhZKaK2VMeGcyJKGoYmali\n9pRHBnMiihpGZqp4f5A4ANixceNO2Gx2FBc7NG9DIBwzJ6Ko0TA2XVSU7ZapMkyTMevGk50DBlyG\n8vIs11CLA8AiAPk4dgwoKfEeP9d7fxhmsxARNeIva2bEiM5YubICq1f/gGPH3vd6ns2WjYULc1XJ\numFBZyKiCPmb7Fy5sgILF+bi+ut/5fN5DePnRkyWMpgTETUSbLIz2ESsEZOlDOZEZDnFxQ7YbHak\np+eENTkZLFgHm4g1IutG0QSoEOJtAMMBHJJSXuc6lgPgCQCHXQ+bKKVcqEUjiYiUUmNxz5gxGW6T\nnU7OYD3M43X8TcQGe74WFE2ACiFSAJwC8I5bMJ8E4KSUclaQ53IClIh0Y7PZUVKS5+O4c3JSqUj3\nd4n0+Zos55dSlgkhuvm6ntILERHpQa3x6uHDUyNKJYz0+aGKdMx8tBBigxBirhDiIlVaREQUAe5n\nHrrXAFwBoBeACgAvqdIiIqIIxOp+5mGvAJVSHmr4XgjxFoDP/T02Jyfn/Pfp6elIT08P97JERAHp\nuUpUTaWlpSgtLUVVbRWW710e8vMVrwB1jZl/7jYBeqmUssL1/bMA+kgpH/LxPE6AEhEFsfv4bkxd\nOhVvr38bZ+rOADlQfwJUCPEegDQAFwsh9gKYBCBdCNELgASwE8BToTefiCi2/XDkB0xZOgX/2PgP\nnKv3Pd6vBPdmISIywNbDW5Fflo/3N7+Pelnvca5f535Y9eQqlo0jIjKr9ZXrkefIwyfbPoGEZ2xM\n65qG7NRs3HLFLYiLi2PZOCIis1m9fzXyHHn4/HvvXJGM5AzYU+xI6ZoS9uszmBMRaahsdxnyyvJQ\nUl7ide6O7nfAnmpH3859I74OgzkRURgCFZ+QUuLLnV8i15ELx27PTb4EBO7pcQ+yUrLQq1Mv1drD\nYE5EFCJ/m3lJKSG6n0JeWR5W7lvp8Zw4EYcHr30QmSmZ6HFJD9XbxGBORBQir+IToh7lTW/CQ0se\nwMm1lR6PbRLXBL+7/neYMHgCrm5/tWZtYjAnIgrR+c28RB3Q8wMgNR/osAUn3R7TNL4pHu/9OMYP\nGo9uF3XTvE0M5kQWpHcx4ViTkFQL9Po7kDIZaP+Dx7lmTZrhqZuewriB49C5dWfd2sRgTmQxahRn\nIN9qz9Xi7+v/jk1pfwfOeA6niLNN8dvL78ErD89Gx5YddW8bFw0RWYxaxRnoguqz1Xjr27cwffl0\n7Kva53Eu/mwiulX2Re4dL+DB39yh2jU1KU5BRPqKZJjEiGLCgUTzkM+pM6fw2prX8NKKl3Dw9EGP\nc+2btcdzA57D032eRpukNga18AIGcyKTiXSYxEzFGaJ1yOdEzQnMWT0Hs1fOxpHqIx7nOrboiHED\nx+F/b/5ftGza0qAWeou00hARqcwr7Q1AeXk+iooWK3q+mYozRPq76O3Iz0fwp6//hK4vd4X9a7tH\nIO/cqjMKhxVi59idGDdwnKkCOcCeOZHpRDpMYqbiDGYb8vHn4KmDeGnFS3h1zas4ffa0x7luF3XD\nxMETMfKGkUhskmhQC4NjMCcyGTWGSfQuJuyPmYZ8fNlftR8zls/Am2vfRPW5ao9z3dt3R+bgTDx0\n3UNIiE8I6XWNmCdgMCcymTFjMlBenuUxPOEcJhlmYKvCo9XvEmmw3HV8F6YtnXahqo+bnpf0hD3V\njnt73Iv4uNDvIALNEwDQLMgzmBOZjJmGSSKlxe8SyaRqoKo+vTv1RnZqNn7zq98gToQ/nehvniA7\n+wlUVXXUbDKYeeZEFFXCyaPfcmgLJi+d7LOqT/8u/ZGdmo1fX/VrCKE4rduv9PQcLFmS43W8bdsH\ncOzY+4rbzTxzIrK0UCZV11WsQ35ZPj7e9rHXOfeqPmoE8Qb+5gkA35Onak0GM5gTUVRRMqm6ev9q\n5DpyseD7BV6PU6OqTyD+5glat26BY8cCtzsSDOZEMSpaV2YGmlQt212GXEcuFv/onceuZlWfQPzN\nEwDA2LHaTWxzzJwoBvmaRExOzkJBgS0qAnpxsQNFRYtRUxOPxKRzSB3ZDiW1/9atqk+43NvtDPJD\n/b7foY6ZM5gTxSArbMYlpcR/fviPIVV99MAJUCIKKlpWZvpSL+vx6fZPkefIw7rKdR7n9KrqY0YM\n5kQxyOwrM32pq6/DB1s+QH5ZPrYc3uJxTu+qPmbEYE4Ug6JplenZurP456Z/YnLZZPxw1BxVfcyI\nY+ZEMSqUyTgjNFT1mbpsKnYd3+VxrmXTlvjjzX/EcwOeM6Sqjx44AUpEUS1QVZ82iW0wpt8YjO03\nFu2btzeohZFRmhLKCVAiikqnzpzC69+8jpnLZ/qs6vNs/2fxTN9nTFHVJ1xaFutgz5yIDBWNVX3C\nFUpKKHvmRBQVjvx8BAWrClC4qhAnak94nOvcqjNeHPQinrjxCTRLaGZQC9WnZUoogzkR6erQ6UN4\naflLePWbV3HqzCmPc9FS1SdcWqaEMpgTkS60quoTTbRMCeWYORFpavfx3Zi6dKomVX0iYdRGY0pT\nQpmaSESmEKyqjz3Vjrt+dVdEVX3CFQ0bjTGYE5Ghth7eivyyfF2q+oQrGjYaYzYLERlifeV65Dny\n8Mm2TyDh2YHTqqpPuKJ5ozF/GMyJKCKr969GniMPn3//ude5hqo+VZslpo8qQW5tmSkKYUTjRmPB\nMJgTUVjKdpchrywPJeUlXufu6H4HslKy0K9LP7+rHtes2YwVKw4YUukomjYaU4pj5kSkmJQSX+78\nEnmOPCzZvcTr/D3X3AN7qt2jqo+/8elmze5HdfW/zv+s9wSk2Tca45g5EakuWFWfB659AJmDM9Gz\nQ0+v5/obn66uvsbj5/LyfBQVZesWUIcPTzVV8I4UgzmRS7QWOFaLr9//17cNDljV55HrH8HEwRMD\nVvXxNz4NeI9PR/MEpNEYzImg7W520cDr9xd1WF93D5pteRS7q3d6PLZpfFM81usxvDj4RUVVfXyN\nTzdr9hSqqx/2emw0T0AajcGcCEBhYYlHsAH0v+030vnfP+4scN18IGUyDl38PeC26j6pSRKeuukp\nvDDwhZCq+jS8f0VF2efHp/v3vwHvvrsI5eUX3lu9JyCtdifGYE4Ea+Ydh6L6DICb3gQGTwXaevbE\nWyS0wNN9ng6pqo+vQNl4MU6fPg6PAD969DBdJz+tdifGYE6E6Mk7Vrs32VDVZ/WAQiDxpOfJmja4\n8qfrsLrgU6+qPoHaoTRQGjkBacU7MQZzIkRH3rGavUmvqj7uu83+3B5Y8Sy6HT6Cwpl3+QzkgdoR\nDYHSindiioK5EOJtAMMBHJJSXuc61g7AvwB0BbALwH1SyuMatZNIU77GdfW87VdCjSAZqKrPRU3a\nocOO3uiwpy9aNK/B6Jl3+XzdYO2IhkAZLXdioVDaM/8bgCIA77gdmwBgsZRyuhDiRdfPE1RuH5Fu\nzJ53HEmQPFp9FC+vfNlvVZ/xg8bjyRufVFTVJ1g7oiFQRsOdWKgUBXMpZZkQolujw3cCSHN9Pw9A\nKRjMiTQTTpDUoqpPsHZEQ6CMhjuxUEUyZt5RStlQQvsgAGXT3EQUllCC5P6q/Zi5fCbeWPuGV1Wf\nq9tdjcyUTDx83cNeVX2UTLAGa0e0BMpw78TMmtKoygSolFIKIfxuwJKTk3P++/T0dKSnp6txWaKY\noiRIBqvqk5WShft63uezqk8oWSjB2mH2IatwaZnSWFpaitLS0rCfr3ijLdcwy+duE6DbAaRLKSuF\nEJcC+FpK+Ssfz+NGW0Qa23F0ByaXTQ6rqk9DT3P16r04fvxyABkALgQmMxVsMJqeRS303GjrMwAj\nAUxz/ffTCF6LLMyst6VWsPXwVkwum4z3Nr/nVdWnX+d+yE7Nxm1X3+a3IISvniaQ5fqv829kpiwU\no5k5U0dpauJ7cE52XiyE2AvgTwCmAvhACPE4XKmJWjWSoldOzquYPn2Ja4e8cwAyUF6+CED0rrQz\ng/WV65Fflo+Pt37sVdUntWsqslOzcesVtwat6uMrzRDIB5CNhmBupiwUo5k5U0dpNsuDfk4NUbEt\nZDHFxQ5Mn77RY89qIAvl5TYUFS2OqWCu1t1JoKo+Q68cCnuqHaldlb+uv54m4Oxpmi0LxWhmztTh\nClDSTGFhCaqrX2901NnrM8NtqV7UmDRbumcpch25Pqv63N79dthT7OjXpV/IbfPX02zb9jv07Ztt\nyiwUI5k5U4fBnDQTqNdnhttSvYS7clNKia92foVcR67Pqj4D26Zhzr2z0fvS3mG3zV9Ps6BglCkC\nlBmZNVOHwZw046/X16zZNowe/bTOrTFOqJNmUkp8seML5Dpyvar6oD4O2PwAUJaJg63n40DySfQe\nHn7bzNzTpNAwmJNm/BUlGD8+LaaChdJJs3pZj39v/zfyyvLwbcW3HudEfRzkhpFA2UTgqLOqT/lh\ndTavMmtPk0LDYE6a8d3rezho4LBaKmOwSbO6+jp8uPVD5JflY/OhzR7Pbajqs25OIlYtetnrtWNp\n7oECYzAnTYXa67Ni0QB/QxkZwwZg3vp5mLx0Mr4/8r3Hcxqq+owbOA5dWneBrcju87Vjae6BAlO8\nAjTsC3AFKIVAzxV2jel1R1B7rhbzNszD1KVTsfO4Z1WfuHMJ+MXBm5B724sYcdddHm1r/CHnnKjk\n+DZgvbs5QN8VoESqM2qFnR53BA1VfaYvn459Vfs8zsWdSUT9ivGoXzkWu6rbI+erLLRNaHf+2npO\nVEZbYLTi3VxYpJSafjkvQaRMRkaWBKTXl81mV/T8BQuWyIyMLJmWNklmZGTJBQuW6HLdQE7WnpQz\nls2QHWd0lMiBx1e7ae3kVY+lSSQd0+TaoVqwYIlMTs70aEdycqbi99EIWv7tjOSKnYpjLXvmZCqR\nrLCLpIemxR3BiZoTeGXNK5i1YpZXVZ8OLTpg3IBxGNVnFG7PmIkdNRepeu1wRUPJt8bMvF+KnhjM\nyVQiGU6IJBBdSB90ACiB85/GOVRVVYb8OxytPoqClQUoXF2I4zWelRR9VfUx034f0RgYzfT+GYnB\nnEwn3LznSALRmDEZ2LjxcVRWdoJzywGniornUFzsUNSeQ6cPYdaKWXhlzSs+q/pMGDQBj/Z61Kuq\nj5n2+4jGwGim989IDObkIdomv9xFEoiGD0/FpZe+j8pKz559ZeWsoD37AycPYMayGSFX9XG/NmCO\nVZjRGBjN9P4ZicGczov2rIBIA1Hr1h18HvfXs999fDemLZuGuevmelX16XFJD9hT7H6r+jSm5G5E\njw/aaA2MXMXKYE5uonHyy12kgUhpz37H0R2YUjYF72x8x6uqT69OvWBPsePua+72WdUnXHp+0DIw\nRicGczovGie/GoskEAXr2W87vA35ZflhV/WJRLR/0JL2GMzpvGiY/Ap3qEHJ8/z17Lvc1Ab3fnhv\nxFV9ImGFD1p30Tw3Y1YM5nSenpNf4fxjDneoIZTnuffs1+xfg7yyPHz2xmderxlOVZ9IRMMHrVLR\nPjdjWqGsMArnC1wBGlUWLFgibTa7TEubJG02uyYr/8JdZRjuSr9Qn1e2u0za/mHzWq2JHMjb598u\nV+xdEfbvHi7f79lEU6/M9MeqKzbVBq4ApUjoMfkV7vhvuEMNSp4npcTXu75GriMXpbtKvR7722t+\nC3uKXVFVHy2GEKI1y8QXqw0ZmQWDOeku3H/M4Q41BHqelBILdyxEriMXK/at8DgfJ+Jwf8/7kZWS\nhZ4dega8RgMthxCskmVipSEjM1Evd4pIoXD/MY8Zk4Hk5CyPY84x/aEhP+/K5Im4eUQS+vy1D26b\nf5tHII8X8Xi016PY9vQ2zL9nvuJADgS661is+DWsLty/IwXGnjnpLtyJ1nCHGtyfV10jcLzzZlTd\nsAH55T96PC4hLgGP9X4MLw56EVe0vSKcX41DCApYacjITBjMSXeR/GMOd6jB9uuBONJlFyaXTcZ3\nR74D3FbdJzVJwh9u/ANeGPQCurTuEvJru+MQgjJWGTIyE1YaIks7U3cG89bPw5SlU7yq+rRIaIFR\nN4/C8wOfR6eWnVS5HisCkVpCrTTEYE6WVH22GnPXzcW0ZdO8qvq0TmyNMX3HYGz/sbi4+cWqX7u4\n2IGiosVudx1DGcgpZAzmFNNOnzmN1795HTNXzETlKc+9yNs1a4dn+z+LZ/o+g4uSvItBEJkJa4BS\nTKqqrcKc1XMwe+Vs/PTzTx7nOrTogOcHPI9RN49Cq8RWBrWQSFsM5hTVAlX1uazVZRg/cDyevOlJ\nNE9oblALifTBYE5RKVBVn65tumLC4An4fa/fe1X1IbIqBnOKKgdOHsDM5TPx+jeve1X1uardVcgc\nnIkR14/wW9WHyKpiPphzK87osPv4bkxfNh1z181FbV2tx7kel/RAVkoW7ut5H5rExfz/0hSjYvr/\nfG7FaX7lR8sxZekUzNswT7eqPkTRKKZTE202O0pK8nwcz8bChbkGtIgabDu8DZOXTsb8TfO9qvr0\n7dwX2anZGH71cPznP2W8syJLYmpiCLiPhvlsqNyA/LJ8fLT1I6+qPim/SEF2ajaGXDkEQghFd1Yc\nRqNYEdPBnPtomMf5qj7feVf1GXLlEGSnZntV9Qm2LzqH0SiWxPRAI7fiNN6yPcsw7N1h6PtWX69A\nPvzq4Vjx+AosfmSxz/Jsge6siosdGDnyFW5HSzEjpnvm3IrTGGpV9fF3Z3Xy5GGMHbsIR45c4/O8\n2YbROBREaojpYA5wK049NVT1ySvLw/K9yz3ONVT1yUzJxLUdrlX0ev72RZey1nXM7vN5ZhpG41AQ\nqSWms1lIH/WyHp999xnyHHlYW7HW41y8iMcjNzyCiYMnonv77iG/tvsOhVVV+wA0xY8//owTJ7oA\nuAzAfgDm3Y6WGVXkD7NZSBE9bu3r6uvw0daPkF+Wj02HNnmcU6OqD3DhzspXDxfIAtAZQDaAeLRv\nvx0FBX80TSAHlGVUcRiGlGAwj0Fa39qfqz+H+ZvmX6jq40bNqj7ufGW2OHvk2QByXT1ycwVyIHhG\nFYdhSCldslmKix2w2exIT8+BzWZHcbFDj8uSH1oVHT5TdwZ/XftX/HLOLzHy05EegbxFQguMGzAO\nO8fuRMGvC1QN5ID/Hm6bNnths2WbamjFXbCMKhaIJqV06ZmzZ2Euai+WqjlXg7nfOqv67K3a63FO\n66o+Dfz1cPv3v9zUY8/BMqq4sI2U0iWYB1rYQZ70GB9Va7HU6TOn8cbaNzBj+QzDq/r4y2wZPXqY\n5teOVKCMKi5sI6UMGzNnz8KbXuOjkQa+qtoqvLL6FcxaOcs0VX2sumYgmj+kSF+6pCYC3tdg6pU3\nPdPUwik6fLT6KApXFaJgVYEmVX2YteEbC0THJt1TE4UQuwBUAagDcFZK2bfxY5KT2bNQQs/x0VAW\nSx06fQizV8zGK2tewckzJz3OqVXVh1kb/nFhGymhxjCLBJAupTzq7wEFBTbL3f5qwWzjoxUnKzBj\n+QxdqvoE2zSLiAJTa8w84K0AexbKmGV8dM+JPZi2dJquVX2YtUEUGbV65v8VQtQBeENK+VcVXjMm\nGT2JZ2RVH7PdlRBFGzWC+SApZYUQ4hIAi4UQ26WUZe4PyMnJOf99eno60tPTVbisNRlxF6O0qo8Q\niudiQmaWuxIio5SWlqK0tDTs56uazSKEmATglJTyJbdjUb/RllWzLDYe3Ig8R56iqj56MEPWhlX/\n1hR9dM1mEUI0BxAvpTwphGgBIAPAnyN5TbOxYpbFNwe+QZ4jD//+7t9e54ZeORT2VLvPYhBaM3pu\nxYp/a4odEfXMhRBXAPh/rh+bAPinlHJKo8dEdc/cSluULtuzDHlleVi4Y6HXudu7346slCz079Lf\ngJaZg5X+1hT9dO2ZSyl3AugVyWuYXTRkWQQaGmio6pPnyMPXu772eq7Sqj6xIBr+1kT+cAvcIMye\nZeFvaEBKifhf/uy3qs8D1z6AzMGZ6Nmhp95N9mCmMWqz/62JAmEwD8LsWRZei21EPcoT+uDhJQ+i\nam2Fx2MjreqjNrONUZv9b00UCIN5EEbnfgdzfmhA1AE9PgZS84COm1Dl9pgmIgGdKq5F552pOLCs\nJX5oUonuw40P5mZb9Wn2vzVRIAzmCqidZaHm0EJC4hng+n8AKZOBS7Z7nEtqkoSh7YZj42uXYvem\nIuxzHTdLhoYZx6iNzqghCheDuc7UGlo4U3cG72x4B5v/zzyg9oDHOXE2AXd3uQevjngZv7unCLs3\neWZoGNH79fUBxjFqIvVYMpibaVKtsUiHFgJV9Yk/l4iulX2QO/wFPHTXnQDM0fv19wE2YkRnjlET\nqcRywdyISbVQPjzCDa7hVvUxQ+83O/sdlJd3ApAD4ByADJSX52PlymzuqEmkEssFc70n1UL98Ag1\nuAar6jNuwDiM6jMKLZu29Pn8AQMuQ1nZ/6K6+vXzx/Ts/RYXO7BtWwIA96Ee5/tTUxPPMWoilZgi\nmKs5LKL3sEKoHx5K09+OVR9DwaoCn1V9OrfqjPGDxuOJG58IWNWnuNiBd9/dj+rqhwBkA4hHs2bb\nMGJEmm4BtLCwBDU1rzU6mg8gG0lJujSBKCYYHszVHhZR2vNV6wMk1A+PYOlvh08fxqwVs3xW9Ums\nboMrDwzE5Hufw139hgRtm+cHjfP1q6uBlSuzlfxqqvD3/iQl7cbo0U+E9FpmngshMprhwVztYREl\nPV81P0DCGZP2NbTQUNXnjbVv4OezP3ucS6hqi7NfvYTajSOwrT4B45ZnISGuadC2GjH52TjgVlVV\n+nxcjx6tQnqvzbbAiMhsDA/magccJQs/An2ANJxX2vuLdNVgsKo+8cuuxqb3PwLqL7xPSj/s9J78\n9BVwO3V6Dp06PY7KyrnnjyUnZ+Ivf7k/pNc22wIjIrMxPJhrEXCCTar5+wDZv/9kyL2/cFcNKq3q\nc8sHf/EI5A2UfNjpvTzdM+A6AJSgsrI1WrXahd69n0Dr1l3CzlgxQ4olkZkZHsyVDouoOVbq7wOk\noqICR478y+OYkt5fKBkZ23/ajsllzqo+ddLzA8tXVZ9IPuz0Xp5+IeA6ACyCc6ITOHkSqKrKQm7u\nLWFf2wwplkSmJqXU9AuATEubJDMysuSCBUukLwsWLJE2m12mpU2SNpvd43ELFiyRycmZEpDnv5KT\nM/2+lhK+X3Oi7NnzDx7HGr7S0iaFfa0GGyo3yPs+vE+KHCGRA4+vlLdTZMmOEllfX6+4rZH8/lrJ\nyMhytTHL5/tos9nDfu1oeh+I1OAMz8pjrS498yVLcgD4H7II1LPVYqzUX4+1sLAEW7Z4Pz6S3l+g\nqj5DrhyC7NTsgFV9omXzp+JiBw4frkRS0ijU1HT0+ZhIhkSi5X0gMoquwyzhBGGtxkr9fYCoNcYc\nqKrP8KuHw55qV1zVx+iFNcGGuS5MfL4F5xDLKz5fJ9IhEaPfByIz033MPNQgrOdYaaS9P2nBqj5K\nUgJ95bM7V3lyzxUivegezEMNwnpnZITT+5NSYlH5IuQ6cr2q+kAKdDjUE1cdGIjHbn44qgI5oGyY\ny/vuyXm8bdsHcf31v+SQCJEOdA3m7kFYaYaKmcdK62U9Pv/uc+SV5eGbA994nItDPFrs6ImTX3yI\nQ0e64xBfjnBTAAALr0lEQVSAsZujb5GLkmEu33dPqejbdzEWLszRpmFE5EGXYJ6WluMRhENdzaek\nt6znUu+6+jp8vO1j5JflY+PBjR7nEuIS8Fjvx7DlzSZY+vkcj3PRuMhFyTAXy60RGU+XYF5amuPx\ns9oZKnot9T5Xfw7vbXoPk5dOxvafvKv6/OHGP+CFQS+gS+suSJ+Z4/M1om2Ri5JAbea7J6JYYcii\nIbUzVLRe6t1Q1WfK0in48diPHudaJLTAH/v8Ec8NeA6dWnY6f9wqi1yUBmpmmhAZy5Bgrnag0yp9\nMVBVn9aJrTGm7xiM7T8WFze/2Ou5Vhp6YKAmMj9DgrnagU7tD4dwq/q449ADEelJOFeNangBIaSv\naxQXO1BUtNgt0A0NO9D5GjNPTs5EQUFowTNYVZ/nBzyPUTePQqvEVmG1k4hIKSEEpJRC8eONCuZq\ni+TD4Vj1MRSuKkTBqgIcqznmcS6xtiW67BmEbkduwLPPDGfPmoh0EbPBPByBqvp0aNoJ9Y5e+Gnx\np0BdIgAgOTkLBQU2BnQi0hyDuQIVJyswc/lMvL72da+qPpcmdka+LRfzJ36P/y6a4vVcmy0bCxfm\n6tVUv1hCjcjaQg3mhu9nrqc9J/Zg+rLpeOvbt7yq+uBQD6AsC82qN6HDdck4W7Pb52uYIU+cJdSI\nqLGYCOY/HvsRU8qcVX3O1p/1PFnRC3DYge13AzIOP8KZgZKY6Ptuwgx54iyhRkSNWTqYB6vq8/MX\nydj8yT8BeN7J1NTE44UXbjFtnjhLqBFRY5YM5hsPbkR+WT4+3PIhJDx72Cm/SEF2ajaGXDkEwz7M\nxmZ4D0klJdWZOk/cKqtLiUg9URHMlU72Barqc+sVtyI7NRtp3dLOHwu2eMmsKx+ttLo0GE70Eilj\n+mCuZLJv+d7lyHXkhlzVx8y970Citd2h4kQvkXKmT0202ewoKcnzOp5hs2PCa7ci15Frmao+5Mnf\n395Xeih78GQ1lktN9J7sWwJc9Tq+/MVXKHnHM6MjTsTh/p73IzMlE9d2uFa/RpImlE70sgdPFAXB\n/MJknwR+OQVInQN0roD7VF8c4vG7Xo9g4uCJ6N6+uxHNJA0onehlqiZRFATzp0ffig3n7sbB7j8C\nnTyr+qAuAVj3ewyKS8DfJs3x/QIUtZRO9DJVk8jEwfx8VZ/dk3Ew1bOqD84mAd8+CSwbD1R1QVxa\njiFtJG0pnehlqiaRCYN5oKo+8XUJqFs5FljxPHDqQlUf/qO1LiXpobGUqknkj2mCec25Gry97m1M\nWzYNe07s8TjXOrE1RvcdjR5VffGn+atQ7hbI+Y+WYiVVkygQw1MTT585jTfXvokZy2eg4lSFx7m2\nSW3xbP9nMbrf6PNVfdQsakFEZFZRswVuVW0VXl3zKmatmIXDPx/2OMeqPkQU60yfZx6oqs9lrS7D\n+IHj8eRNT6J5QnPN2sAFJkRkNboF88OnD2P2ytmYs3qOV1Wfrm26YsLgCXi016NIapKkaTu4wISI\nrCjiYRYhxDAALwOIB/CWlHJao/PyuYXP+azqk9w2GZkpmXjk+keQEJ8QUTuUCmWJOBGRUXQdZhFC\nxAOYA2AIgP0A1gghPpNSbnN/3KyVszyed83F1yArJQv3X3s/msTpO9Kj1wITDuUQkZ4ijaR9AeyQ\nUu4CACHE+wB+A2Cbrwff0PEG2FPt+O01v0WciIvw0uHRY4EJh3KISG+RRtTOAPa6/bzPdcxDn8v6\n4LMHPsO6p9bhf3r8j6aBvLjYAZvNjvT0HNhsdhQXOzzOjxmTgeTkLI9jzlz1oaq1wf9eIYtVuwYR\nkbtIe+aKBtxXPbEKQige+gmbkh6xHgtMuFcIEekt0mC+H8Dlbj9fDmfv3MOf//zn89+np6cjPT09\nwsv6pnT3vMZLxBt682qNb3OvECIKVWlpKUpLS8N/ASll2F9wfhiUA+gGoCmA9QCuafQYqZe0tEkS\nkF5faWmT/D5nwYIlMjk50+PxycmZcsGCJWG3w/drTozoNYkotrhip+J4HFHPXEp5TgjxDIBFcKYm\nzpWNMln0FE6PWIu9sLlXCBHpLeK8QCnlFwC+UKEtEQtn9zytxrfNWgyaiKzJNLsmqiGcHjHHt4nI\nCgzZaMtMC2p8ZcAkJ2eioIDDIkRkHNNvtGW2BTUc3yYiK9C9Z869UczHTHdKRORk+p45F9SYi9nu\nlIgoPLpvkMIJR3Ph1gNE1qB7MNdjbxRSjndKRNag+zALJxzNhXdKRNZgeEFnMhZTM4nMKWoKOpN5\nFBc7UFS02O1OaSgDOZHBGMyJiCzA9KmJoTJTDrSZ2kJE5M7UwdxMOdBmagsRUWPGFOJUyEw50GZq\nCxFRY6YO5mbKgTZTW4iIGjN1MDdTDrSZ2kJE1Jipg7mZVouaqS1ERI2ZPjXRTDnQZmoLEVkb88yJ\niCwg1GBu6mEWIiJShsGciMgCGMyJiCyAwZyIyAIYzImILIDBnIjIAhjMiYgsgMGciMgCGMyJiCyA\nwZyIyAIYzImILIDBnIjIAhjMiYgsgMGciMgCGMyJiCyAwZyIyAIYzImILIDBnIjIAhjMiYgsgMGc\niMgCGMyJiCyAwZyIyAIYzImILIDBnIjIAhjMiYgsgMGciMgCGMyJiCyAwZyIyALCDuZCiBwhxD4h\nxDrX1zA1G0ZERMpF0jOXAGZJKXu7vhaq1SirKi0tNboJpsH34gK+FxfwvQhfpMMsQpVWxAj+j3oB\n34sL+F5cwPcifJEG89FCiA1CiLlCiItUaREREYUsYDAXQiwWQmzy8XUngNcAXAGgF4AKAC/p0F4i\nIvJBSCkjfxEhugH4XEp5nY9zkV+AiCgGSSkVD2U3CfciQohLpZQVrh/vBrAp0sYQEVF4wg7mAKYJ\nIXrBmdWyE8BT6jSJiIhCpcowCxERGUuzFaBCiGFCiO1CiB+EEC9qdZ1oIIS4XAjxtRBiixBisxBi\njNFtMpIQIt610Oxzo9tiJCHERUKIj4QQ24QQW4UQ/Y1uk1GEEBNd/z42CSHmCyESjW6TXoQQbwsh\nDgohNrkda+dKQPleCFGiJFtQk2AuhIgHMAfAMAA9ADwohLhGi2tFibMAnpVS9gTQH8DTMf5+jAWw\nFc4hulhWAOA/UsprAFwPYJvB7TGEK4HiSQA3upIo4gE8YGSbdPY3OGOluwkAFkspuwP40vVzQFr1\nzPsC2CGl3CWlPAvgfQC/0ehapielrJRSrnd9fwrOf7SXGdsqYwghugC4DcBbiOFFZ0KINgBSpJRv\nA4CU8pyU8oTBzTJKFZwdnuZCiCYAmgPYb2yT9COlLANwrNHhOwHMc30/D8BdwV5Hq2DeGcBet5/3\nuY7FPFcvpDeAVca2xDCzAbwAoN7ohhjsCgCHhRB/E0J8K4T4qxCiudGNMoKU8iic61T2ADgA4LiU\n8r/GtspwHaWUB13fHwTQMdgTtArmsX777JMQoiWAjwCMdfXQY4oQ4nYAh6SU6xDDvXKXJgBuBPCq\nlPJGAKeh4FbaioQQyQD+L4BucN6xthRCPGxoo0xEOrNUgsZUrYL5fgCXu/18OZy985glhEgA8DGA\nd6WUnxrdHoMMBHCnEGIngPcA3CKEeMfgNhllH4B9Uso1rp8/gjO4x6KbASyXUh6RUp4D8Amc/6/E\nsoNCiE6Ac00PgEPBnqBVMP8GwNVCiG5CiKYA7gfwmUbXMj0hhAAwF8BWKeXLRrfHKFLKTCnl5VLK\nK+Cc4PpKSvk7o9tlBCllJYC9QojurkNDAGwxsElG2g6gvxCimevfyhA4J8hj2WcARrq+HwkgaAcw\nkkVDfkkpzwkhngGwCM6Z6blSypicqXcZBGAEgI1CiHWuYxO5bXDMD8eNBvBPV4enHMDvDW6PIaSU\nG1x3aN/AOZfyLYA3jW2VfoQQ7wFIA3CxEGIvgD8BmArgAyHE4wB2Abgv6Otw0RARUfRj2TgiIgtg\nMCcisgAGcyIiC2AwJyKyAAZzIiILYDAnIrIABnMiIgtgMCcisoD/D/vAGSOvteqjAAAAAElFTkSu\nQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c273190>"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "These equations for linear regression are commonly enough used that they're already defined as a function [linregress](http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.linregress.html) in the standard module `scipy.stats`, which returns the same values for the slope and intercept:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import linregress\n",
      "slope, intercept, r_value, p_value, std_err = linregress(x,y)\n",
      "\n",
      "print 'a =',slope,', b=',intercept,', \\xcf\\x83=', std_err\n",
      "print \"Pearson's r =\",r_value, ', p-value=',p_value"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "a = 1.45955689947 , b= 0.191167623632 , \u03c3= 0.0858274566043\n",
        "Pearson's r = 0.864232048349 , p-value= 5.35418781675e-31\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`linregress()` returns as well three other quantities:\n",
      "\n",
      "First, the `std_err` is a measure of the stdev of the slope estimation (under the\n",
      "[assumption](http://en.wikipedia.org/wiki/Simple_linear_regression#Normality_assumption) that the errors $(ax_i+b)\u2212y_i$ are normally distributed),\n",
      "$\\sigma= \\sqrt{{1\\over n-2}C(a,b) \\big/ n{\\rm Var}[x]}$ :"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sqrt((1./(n-2))*(a*x+b -y).dot(a*x+b -y) / (x-xbar).dot(x-xbar))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "0.085827456604302785"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The other two quantities returned by `scipy.stats.linregress()` are known as the\n",
      "[Pearson's r](http://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient)\n",
      "correlation, and its associated p-value. Pearson's r correlation (or $\\rho$, for Greek rho) between two variables is related to the above covariance, Cov[x,y], except normalized so that it ranges between a minimum of -1 (for completely anticorrelated) to a maximum of +1 (for completely correlated):\n",
      "$$\\rho_{x,y}={{\\rm Cov}[x,y]\\over \\sigma[x]\\sigma[y]}\n",
      "={E[(x-E[x])(y-E[y])]\\over \\sigma[x]\\sigma[y]}\\ .$$\n",
      "Dividing by $\\sigma[x]\\sigma[y]$ keeps it symmetric in $x$ and $y$, and also means that when $y=x$ (maximally correlated), it becomes ${\\rm Cov}[x,x]/\\sigma[x]\\sigma[x]={\\rm Var}[x]/{\\rm Var}[x]=1$, and similarly is equal to -1 if $y=-x$ (since from the definition, ${\\rm Cov}[x,-x]=-{\\rm Var}[x]$). This also agrees with the output of `linregress()`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "(x.dot(y)/n-xbar*ybar)/sqrt((x.dot(x)/n-xbar*xbar)* (y.dot(y)/n-ybar*ybar))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "0.86423204834868328"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Pearson's r measures something slightly different from the slope. If two variables have data points  falling on the same straight line then it will be 1 (or -1) independent of the slope, and it will be between -1 and 1 if there are deviations from the straight line fit, again independent of the slope, as illustrated in this figure (from Wikipedia):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from urllib2 import urlopen\n",
      "from IPython.display import Image\n",
      "url='http://upload.wikimedia.org/wikipedia/commons/8/86/Correlation_coefficient.gif'\n",
      "Image(urlopen(url).read(),width=425)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {
        "png": {
         "width": 425
        }
       },
       "output_type": "pyout",
       "png": "R0lGODlhwAMSAnAAACH5BAEAAFoALAAAAADAAxIChv///wAAANDQ5wAAgIyMxr293qen0+np9E1N\npvDw95qazcfH49nZ7OHh8GhotHx8vbKy2YL/gkb/Rvr/+vH/8SD/IAD/AG//b9H/0Qn/CZ7/nqz/\nrMf/x3n/eRz/HOj/6Er/Sjj/OOz/7Cr/Kmb/ZtX/1RP/E5T/lLX/tQX/BcP/w4b/hgAcAFP/UwDR\nAPb/9i7/Llj/WN//3wDsAABBAO6cANr/2mH/YeP/433/fUH/Qb7/vg7/Dk//T8z/zJD/kBf/FwDD\nADP/M7D/sABhACX/Jaf/p6L/onT/dGv/a5n/mYv/i13/XTz/PABTAAD2AAA4AAAJAABrAADjAADM\nAABPALn/uQAOAAD6AABdAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAf/gACCg4SFhoeIiYqLjI2Oj5CRkpOUlZaXmJmam5ydnp+goaKjpKWm\np6ipqqusra6vsLGys7S1tre4ubq7vL2+v8DBwsPExcbHyMnKy8zNzs/Q0dLT1NXW19jZ2tvc3d7f\n4OHi4+Tl5ufo6err7O3u7/Dx8vP09fb3+Pn6+/z9/v8AAwocSLCgwYMIEypcyLChw4cQI0qcSLGi\nxYsYM2rcyLGjx48gQ4ocSbKkyZMoU6pcyTJWgJYwY0Z7KbMmAJo2c+oEhnPnyp4+gwqVBXRoyaJG\nkyolhXTpx6ZOo0q1BHVqxqpWs2pFhHXrxK5ew1oFK9Yh2bJojZ5Nm3At27c2/93CJSh3rl2Vde/+\ny6u3r0i+fvUBDkxY4+DC9Q4jXixRMWN4jh9LVhh58rrKljMLxKzZHOfOoAWHBvh5tOl5pU93S626\n9WXXomHLLsh6drXatnOv1k0PN+/f1nwDZyZ8uHFnxY8fS668uTHmzoVBj07d1/Tqva5j325LO3dc\n3r+LdxV+H4MBAwjsLD++PSr2+AqgH2BgPcPz6ZElMIBgAIL67uUDnz0KzDeAAPYpJB96ABYjgIHo\nIZBAgPcMOA8D/RnIQIIIFTgfgsUsAOF8CFBoj4XwHEDAiAP4hOI0GEK4YSMCqNdKjYXg598BABgw\n3wIm9tbOfg6g90ABrxzA4v8AD7hoEAMrjkhjkS2uIgCVhUQ5wIwAzGdjkPG8aMuDIzowoSJLSjmJ\nj/4VwCZ9ThKUwJJNJpJAARmi9wiRRiIpyZ15VimIjl92OR+Y8ohJi45lLpImhJQQYMCMb/qpk6LN\nsPnfgnAesh+LjpAJoQM8OvKpmoJ4eCAhBiIa5joZPsBjjOgBmcijrWKiZamXFiQppfNZOgiUZdrK\nCKOjNkLsqMYKkusgz7raDqaxcMorobFEmxO10GjJJQACPDCiAbw2EuuGB+TZrCHhjluuIKLWOSiJ\n0r7D7StUKlCItqsoiZ4DQd3rjLbMSmLtsF6iaaAD6xLCqbAioldoveoIzIr/v6uyeuituOqJyJzJ\nGhLxAPrGeRDGAGvMICX57rsxIgY2iIiWzUIwn8wUp2PxKiOfKQjIHsPcsaDsjnvIwwHrA/R8KRci\nasnQLlzAu4tgDGLUQR8y6tSI5FmuuOh9mzM6O6uiqtMJc9xxIm/+eIi3Secj6s1HB1vIAlrOp4DY\nifRMyNKK4D3i3oYQ/PLYOqdDJdEAaClsJ4sfTkieL8SNT9voXU3IroccoCrTDRtyNiEjT9z55/82\nG+3IUCNOdjoaDjI336lgXKLl90SOns+Tz0f1IAdg7t8ii2dpdyPBj3i7oVmrGrrr45SNCsZMCrJA\nnq2zgi3uBNneCJ78uowe/68jD/C7IuC/zDkAKENfMTqc0lnt8SYXtL2yWCbCqKwA0EoyJQzIHwBG\nVp8DgM187kvcOVSFuurFAm7cGwjSIIE39OFKXpSo4CACZSAIJFCB5ljcATxnpOetwmtCkd4xILgJ\nBo4Ig5pI14g8+MHXnSN8vFGhMVDICRGSkEkmvETw+vOf89UwHDokxf2Mk0Rz4LAUlTviOZooCk7R\n8DhUHMcSpRiRLILCec3xYjg4hTMuPkSMZkTjN1RFOzMyRI1ShGM3RBhFN57RjsvBox5BIcca9nGP\nufljAgUJSNkQEnqHLGRrEok4RirSNI7MWSQfCZpJ1suSlMwMJl21yUxKpv+TYAKlJxcjShOVcpSE\nOWWAVInKvrCyPa9spV1iKR5ayvIttuROLm+Jll1ix5e8DAswqTPMYGqlmM5BpjGnokzlNHOZTnkm\nE6F5S2kOx5rUHAo2f7PNbNbPm6zoJji3Nc5SdEAChxBnOWWiztGYwAIkMEQ719mSeYLGBxawgBIK\nYU964sWforBCPn1AiH4C9CQG7UwHLGACEQwioQclCUQ1IwELCOGhEd3jRDMjggrAUxAbzehTRAoK\nGbxznyElKUdSahmBWsAHLFXpVWT6iRXk8wo05WJMLdMCC0whp3EEqieEYIEWCPWDO7WMDGZggRUc\n1X1JtQwU8rmBp7ouqnr/qYFWtSpPJTBUBprYKletqhus2kWsW+UnAEhggQo49BJoHStZZ2NWuMQ1\nrRgFAFHRaYm7ynWurqnrW/xag4IKQgTv7EBf/QpYQzZWEYQ1rCDwaYGqUiKyj11kZhGBWZAS4ggD\nvSxjN3sawbKlszcpBFvdOgnUkraSrzXEXdVaiIo2oY6QmG1sR2PawYpVnoYQgQcscAHR4nW3oelt\nYZAig3zuE7nVUW4qEbGB0EI3OtINDFQWaoIPXDeZ32VEVSp60fBi0byKwIpHjYpe4GTXL1iRQQb0\n2V5u1jedjHCpFe4bSP7SdhE2zYB3/Qub97rSET1lLYE1u+DUOmK9DVaN/4H1chbEWiACEYZkhNei\n3wwnd8OR8GoGSuDhzkz4LnJhQltxW+LHnHiWkyBqE1psmRfPpS4WViyNXQziSVDWsjtGjI3hwhfQ\nWoDEQV5ujynB3QEnGb5LpsRen6zdKE/iBR6NJ5UpbGUfO3fLKO7yJFCQTxWA+cZingR333rmtAwZ\nl5uo6AhY3OakgDURb2ZLZbJcZ6fYlAR3Bm6DK2PSpvZZKQvNJ6AFvWDMkNkCZj60UUSQ6KISNK+N\n9oRNLRBoSQtFBJu2gAQunWc3f6Knc/a0UlbwTovCNM2YGIEFQKDqpbA6n1TYb6Y/8YH5OrXWSlFC\nq4UA5Pum5sfAXooTPP/a1mKjlzU/YCiSk61NAGyA2RV4LlqhW5vVspnaJtsAURmq29jiZsrg5p4V\nhODazOImx+n+piDa/VjfIDve5JTtaF8rHCNfGt/sPARhL/Dtep+CrR4oOMB/InC/KrrTgE2Obem8\ncJQ05a4yUHE+mRDobQM1OcIlbsXrmQiPA4DS+bR0uWXKnObSd+T/lIQILjBfekcUOtV9KcxTwhcR\nRMDmB51Ok3eOULjuW6XXIS/RTTIYoANUOxBe+kia7tcKoIDlrZDvy6UOksPcFdvPzah3Osz1kWLC\n49fOZ7bF/ooAO7nsG2ENCsZtgh9Q3JvlSbDC4W4R3Kw7nyZYwd2hyZ7/qPN9pqVQQUUtkIEV7H2Z\n7LEwhg+PkeT4YPEZiMDjeQkfslO+ItApQU8V/XZjDkjE//58Y1ghg9HDs/TVpMVqB696yrjiA2zN\nZwsg3koUybj2q3+FCHJfVN57EkXwBv4dY0HpVo9ali+6t/LfSAuU5/P5oxSTkac9/bbc4tYW1TUl\nFTX07nsfF+AXgvgLiSl0m/8g3BK22p2NR6//dhNY/uj7DSKwI2Cb/mlkdCaXCZQVdvu3Gb+Qdm2l\nBJs3SIu1cpjwaJF2gAEhPeKWTx6gBLR3VQ8YV52wZhRYgcPwdwzleDrVgQOoCXK2gSEoJMTgA3Rn\ngn6EgvfXCXzWgv2Q/0WXl3IdAHtjAxhOVwmF9ms4GBvH4AOuR3BQJYA16AmPtn5FeCLLIAPEt2iN\ndHZNWFOMZ3xR+CrMQIUpx3E/GAyo1oUV8gwzl3ItkHqIokKyxlePkIJmqA2ENHPzZWnSokIi4Gtx\n6IFzyA2JBGopJwETaErDIH2QdXR/eA2RFAHDVojukUTRZgLcx1mKuIi3gQ3yZ1FXt0rF4G2LEISY\niBza8AOtVgEA+EvG4H6WCIGjOA2gpIBrV0vGkHyt6IevmIndIIsGGF1HSFWJKIe5OBPgMHeAx4Cq\neAxepXMll4XDKA2vRIKB54vIcAEWkHDPSA60tIP5pHnglQzkxYLZGP8c5sCNGdABDUhXyhByxTWO\n4OBLrZdyJOCDtqFGLteL7kiH6hCPD+dezJBzbJiP2IBMuKeGdyaMsMUM3MWFAgmL7jB8auiKCckM\nSteQ2SBNzSeKPOYMhmeRuggPIqCRjPFHhYaPHgkN2ySSpAQNnneSxCgPhAWFJhYNNmUCDOmSyZCS\nDoeKvCUNeoeTKEkPd6UE/zeR0NCRQLkM4mRyvJiOUCYNFkaESakMG3WBDGWSSqYMA9iSU5lHACGN\nMriSybByyxiQXTkMMQWDgBeW04UMqLVaZ4kMScWNFtCDWXkMrkVU5RWXxIBVSJhySviUeKmItsiX\nwVBXYNiPXOaWl1j/AsBomIe5EIlpAWIIY4wpkct4k5DZHQ2RhrpnlmKBRhIpCAu5mdbxEHb4mXCm\nlc5YCKxomrlgYIJ4fZDoFZ2Uf1oGm7FJEY6YT+pXFqBUgLq5mxWxib9ZEwGQnMq5nMzZnM75nNAZ\nndI5ndRZndFJAxYwA1BgndzZnd75neAZnuLJnI01nuZ5nuiZnuoZnkTAVBZABTSwnvI5n8kJQtsA\ngsNZC6VmbUU5Ft9Qkfk5C/spCLKYitX2DTcYoC7hEdKIlRGkDSWpoAv6EWC5FKzElRKqCgN6CHTJ\nlvLGDTWpmRkqChuKCNxoAuiYQuOQYCMaTibBjyZAAk4Zc+JAVOzV/6LvgRL8WJccEg5DiKM5mhJg\nqGP5Ng4uZaBA2gkl+ggiMKMWZw5eZZNJyhSw5g2gOKWhsKRb4UuviaWcoKXHhA6F6aWbAKZZAUyI\nSKaZYKb+mQ4AqaZlWqXhcKVweglsykzsIAFEWqdUwad++qeAGqiCOqiEWqiGeqiImqiKuqiMGoJ3\nOpWP2qiSOqmUWqmWeqmYmqmauqmc2qme+qmgGqqiOqqkWqqmeqqomqqquqqs2qqu+qqwGquyOqu0\nWqu2equ4mqu6uqu82qu++qvAGqzCOqzEWqzGeqzImqzKuqzM2qzO+qzQGq3SOq3UWq3Weq3Ymq3a\nuq3c2q3e+q3gGv+u4jqu5KpIkVquM4mutHCuPHd8VsWuT5pJ8Np1mzWvOfiu9eoMCGlqT2WvR9EM\no9lL+OpuzKCSbXpU/ipRBXuJ/Iqw+boMBoun/fqwrBmwoTmwBgexDCuwE0uwC7uvwImx5amvraln\nIhtx7tqxGesaDcAf/kEu63qyc5WwvlA+EdIAAgoOLUtEMKsNEWM6/ECzpCMxhuEa8WMgywMLp2Sz\n/oGz2EBGCBgTUIt4p6EjCMAjbxJEQcoNR0sv15A3WnuGMAG2Rasau1I5QAO0qyBKVou1brMIGsQK\ncQtAHNRGRqgS/hM2ZSsUwUMlD3BFo7BET5QKonS2AJC2gSNAqbD/AIqrMoqwLJAStRwhOTmSN5QL\nekPBtANAKo4yNIsARljDOK0ASoJ7uQcAAYHyCH1rJICruqiLQ5crCNRjIDDkD8q0ukzSuqUQu+wj\nP3u7E8iyMJ27NgqTNcwjui66DaDrLJKTPKgCtyzCuY7gvJFbOMZLCJoCATazMpKbDporvWqDK4/A\nu9m7vZ1CtTsRKwfwAg2gLsUrvn0zHw9QRw3gtUrLDYYTNJALOo4QvEyjLJa7uSbEu78iCG+iu/eq\nDv77L8MLv41AwJNiwPOBwBQhtIlgvteSNqBgvoBbOkSxDSMzv4NQvxEyQAd0M+v7COrLvu7bNyfM\nICmsNo5guCKo/w4r3L5vKzTE+8DXOzO+I44NYcGIkC+4Nbi68iGEwMEfrA1KPAgeXDBALAgYjDBE\nG74CDAAU57mH8AJGLCDqMMXzUsWgoMWGwMW8+xURAbKV8AJWIz7IG7qPkggvkCdOKwgHFLam0Ela\nojlKHDMsgx4KUMRn7MfvmyaIgDJR7ILoQMQq88aaQMaFgMhxBxEWWwl+AzyDOzSOfLyC3MPJiw1z\nPB91DAB3fLybyzWOwMZI3MiF7ACobL0OjDaATBfooMqZ48aFbMhWvCSI8DQr9RAR2wijIwgvsEVu\nHMfhO7SzvMTZYCCdXCWCAyGE0wiXLLuDG80GMs0yzAhN3L3lUP/NvTvIkDy8jdDN6LsQwcwIi4Nb\nWgIBiRwJ0bK896sNq6M3g8DGDcQwjDDMWGzMhYDPxbLNivAC6+PN5MDPxazBWjPOxcsIBO07v+wQ\n6TzQBlLH5WNEmPDQ42PNDJyz9Hw4bIRbL3AqSEs888HOE8yCIy08Seu4jMBDBj0O67w5KW0KZ7xB\nED3JEr2xlbA/E3I99iwKBMQ+BxTDEwrK69M+F5y6clzRTmwgGH0Irxu7Nx3Ow0PL52DG6GHRUH0K\nVe09Ed0QE50IXQshIjwKHEQ3MeuzdGNA82G3w0Ilg+fTA4Q9ACTXC73JYZwe74wa6EDXQL3MNu3J\nhEAofW17wIz/iy2kNy9kCjIEIRRMHtyQ1ugR2feMxXNrCGVNu4ddOS+Q2ZJgvo9DGvBzQYdNCqLd\nEQamxn88Pj/0AHisCWzsskV0C6D02B3ECS5k1g/01rSBDrtNu7nAQr9bE86Mxa7g2cjNmdsw2zwb\n1ZLgQx4C27KQJwlw2pCxyL7z2rHNCtaN3fCnE/7sGd+g3OCd1wOg3LZw3l7oRCe93L3A3gghxA4z\nwdMCEOOdFLuU3yOZEy8gz/bJD1Z0pucw4EYpEyJ03/7w3zksFbnE4LXyYToRRfJNDSGV4AdLDhie\nruMXEBRO4FlNzBLe4So7svIqs2RF34mB4iyuRyredxT7SC+O/7keS+IOW+My3uI3buNCNeM6zeMf\nF+PmquM9nrI7buJAnlM+XtzsR+RBfuIljrJQfuRSnuQ0teTn3ORRPrNGXuQ4PuRbnuJd/uQrm+Nh\n7uQBeOZqXn9ofuVjruRCruVUzuVT7uVlDuZzLuZ1TuZIngqsHVhtjnWqUMmAvuZ2fgpjXWOBjnSD\nztMFtugkBR2JPhlY3kVxTgqTPg1/LrZ5DumZkOnRQOicnhFOqtqXjumKXQ6gHuAXIQR7ynSnjurH\ndQ6rbkMZQVQVUOHLd+cWUetTpBHlF6+8XhG+Tt4YEQH5JKJmN+wUUezlkGfLKJOw/uUXIepSaBE5\nh6RTF+vNXv+yQWsRJfBOUkmjfY7nElFouckSlX4R667gEtFRohZw1A5IJyYCeomc3G5HJ/aT8s7s\nLi4RbNVQcZHvbmRgmwaaDDfvGgUR0d4rCv/vDqECj1mk5S7nC+EDfNijVW7mC2FhN+rw/i4Jm/6N\nCvEBw0VrD+rptzjr49Fb966iD59bjk5MC8GiBxryfWjtzqQQ1pgByk5yMZ/zqf4ddYXszHjzFd9a\nM49dCJFz0q7xdM6EI39eBnGkFkrwjqDzO18Q4c6jVx/0Ms/yLV8Qhfbx+o31JzgQ8A6HX4/zBT8Q\nLwCg0YT2QSUQPSXwDk73R8RSAf/zIN8SVuD3ZvHm+3DwIK7/EiKwAS2QWKa+5/vg9GHaEXIoAyuw\neGrnoBVM+N4w9YUgnJG/EbqlAh3AbNfnoUxO7+Sg9YNAWa8u9IWl7hxBWHfIUC1gBboe3o7fDREr\nA8PFBMbl7dsO+jvZAQjf+FbODQYrAm/4+5x/+sTuV/QYEu3ODrFksHIfCc4exBiR+CSQAdm/ENP/\nGuKgkndf6vqm+s7fEJRv+Rbw/Yh9/Nsgkj0v+Oc/9P9aEMLoA6OfcqK2At6F/oAAIDhIWGh4iJio\nuMjY6PgIGRkQSVlpeYmZqbnJ2empOfkpOppYY3pa47hiYeGTiQpLKjtLWwsQapurGwlrKiiy0WLC\namHSgvJS/xi7y9zsXIj7LD1NXW3NGH1NvcxoxWq12ZuqTV7+mG2eTiluKkFsUdGhok5ff26Pn69P\njr5/7W3hRydu/goy62fQHrsarCT8EJEwojmEEitavEiIIsZaJYZ12AiSnsaQ0xZugEgyZa2RKlu6\ndMbyZaYPw1rIvLkrJs5RMpRIWLgz6CWdQosabUT0qKIXQixIUAo1U9KohggC8NGhKSugVLsKmuo1\n7EuwXkW4qyA2rSKyXcW92EBiGCtjJ8Wpjcr2rl6JeaO2KIZyr96+ShcSM9EBHCGrgncSbgy53mOj\nHYrJiHx3ctGFTVZcxixWM+jR1ETvVMLKFemwpm++QMF1tf/X1rJrryS9gdUG211pq+zpzkJs3kp9\nEz/OybhKH6xWIIeqHCNWrayaDH8eNDr27fcwM7fwkXtR7QmBxX3X4qQgu+Kzt39/ELOIChaYwHcf\nVoTPDIcTH2J8n0rkBQjfgBaJ0NRTBN5kYD3TvdPZZwtS1eCE2FUY0VmBWSigUi9Ycd5cIKjHYW8l\nnjgUZH9VsCGKIGFYzQc+vVMBCYq5SCGOOiLV2CoWSLgjRjA+40ME1FkghGdBmrhkkxkJhlorTgr5\nkgggysVKei1OadSQXA62F0C7fVmRl7To9xcxHthIJpNt6mhmPQA592ZEcYryIDER1ukmnyfemY4M\nHvmZEKD/m1hxAZYZSEAioXg5+qdaH9BnE6T+GGqJjMHN5Z+lfXoaIKbXMIUkqPuI+kgJKxwpRAdA\nmgodrASiSo1Z8MiaD62KwIWlBVriOhuwBYpFQjEfCCuZPR9sAMI7JrCJbLDtAeinrs+sYsKr0fKj\njgwdNPFOktpuG6t47FlqbTNRokCuOrSi0EGvLSixZbuPcnddtV3lZsGY9pZj6AdHSMDfXBfc+G9o\n5rLjabq6lNBcwhNNI0ME4BKTpGoSq+VwRPny2XEtgoK3McDObJAoevSWvFfIBn1cp8uyzOcUy9zm\nAkyac0Frc8sLn+uozKMgWGrP10zGmLebIrmCxkb7/DMq/6AKLYqGT1vzGHtWxPuOBCtfDRnV/lAL\n8lHFelAv2PFhYthhPKsdNtxxF+VjCXJPQ5hhrTp9d2Ni951TUVEiDHjgbLNzbOGY/a34bUGJ2Xgz\nFFmVs3AMRz435p++9B2dmueCkF0yrHDxVpd/Dibq994kgwcWkKA66P+x0/XKQMeuMO5HMV7JpDXr\nTks/Cz3r7yBkA18c8uO51pQQys+CjgwwPw819Qze1Cxa1o+Cy9b0Tb99WryHD8mQx29ygbHkfxJA\nMFiaAP76m8u/EYy3exJBaoucT/7oS1fQtPX0gn6CGR8B18KM+FlicPsb4AEB4L2urWBL/Hvg7iz4\nogServ8TkFOEAnGXs/e1wAoGxCDeTEilXXwQEh2xgOdKscHn+Y9GHdBYCVEIExxaBEMrdMTIKtXA\n++luBfQhhte2dEMdGk6JhdKgEC1BMwUxooeNw1Z6EMjEyCQRhfZz4CZI5TxHUFFxMuAbIraYxeCl\n0SDmk1onVpQ2RDzxgGhcIynqaME6FssEcYThKZiIRzuyT5Cn2oiPzEhIbCQyM4vExxYZ2MjuRFJa\nk0xHEr9TvEpiUZOr46Q2btg5TypSlMkj5dEqIgNKmXKTq8RPK08YEd9J8ZXQoKUrbSk5ibgjjLh8\nUi9lEkg6RgSOv/RlMVsSTALyLn0/OuYgkunMW2ykgrj/Y1z+LGC3aEpTmyTZ4hyR9zcUfIOb2yRn\n/S4yRsyJDSBHMCc0nZnEdGKGmk3URwvDQ853HjOeMXzON8ukjx+a8ysDPadF5NkYhNpDaPSZZT4L\nmsKKKHQvExUJPmzFy4Hqs5je9OJ2KuoufPzFBIkr6EZ/iUZ6kgaklrTHHsfFzZP2Uqb+7GdE0xEl\nRD4Uojt83j/5Qg9+Ec6kPO2pT93YoXSEsqgEZaqdnNoJa30Hn0ylqS2t+jxdicB19oFqOb1aSLBK\npRxEc2hVxXoptGKCVmdBK1Zf+VZwkmOkfSSqWnN1V0uIqjIZgGlR47pKwOoOUzm9q2BJedjYGUqc\n/cpr/2I9+VjUAWpOef1qZVt62fJR456XjawmPau5OPkOiI7NrEVN6wgzkcqsagXtJF0bOS/ZikWm\nha38AoDb3Op2t7ztrW9/C9zgCne4xC2ucY+L3OQqd7nMTe4zvERX1Np2kdMt3JDyly13Nne73O2u\nd78L3vCK97s4qe62xjteJ7ACCuhtr3vfC1/wpsS8gqTv3TDUQdRaVr9Y4+8oc7FU/tp3jQOGW4NE\nwAqq6rfAWWQw2AxEM9IK2L83o/AZc9E8CzdVw6Xh8CEGZLWqINWtHu5wiY0piz3W66d2PXEOXbxh\nUmBrXCylpYOVeOOnRQeSImaxO2H8YhgrJ7899qhXc4qMQyT3zDgBnp2PdwrkJbrYNyIYlAdtetYo\nS/nEtIlwI54MZS3bQsk2882qHqHSMItZjUAm8/rcnMcowzl8c14z/epsPTzbmXx63rOf/wzoQAt6\n0IQutKEPjehEK3rRjG60ox8N6Qn1OdKUrrSlWTbpS18t05rutKc/DepQi3rUpC61qU+N6lRfIhAA\nOw==\n",
       "prompt_number": 7,
       "text": [
        "<IPython.core.display.Image at 0x10c3eff90>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The Pearson r is also used sufficiently that it has its own dedicated implementation as `scipy.stats.pearsonr()`, giving the same values:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import pearsonr\n",
      "pearsonr(x,y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "(0.86423204834868272, 5.3541878167471908e-31)"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The `p_value` returned by `linregress()` (or `pearsonr()`) is a measure of how significant it is that the correlation is non-zero. One way of assessing this is to randomly permit the y values and to recalculate correlation:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import random\n",
      "yr=random.sample(y,100)\n",
      "figure(figsize=(6,5))\n",
      "plot(x,yr,'o');\n",
      "print pearsonr(x,yr)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(-0.032637855742317069, 0.74717917356035102)\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAE4CAYAAACpJRM9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG09JREFUeJzt3X+MZWV9x/HPV1AWaiMVDYtKCrmpKVZbNNVgjcyNgbmk\n26D9B2s1pYY2ptVdahMDuzvrTAUTf0Tr7Jgaq2CwtlijDYG5dZm1OjM0VYMV5HfUDbSA7ELXWqFh\nidBv/7h3d2Z275177j0/nuc85/1Kbnbm7szc555z7uc8P88xdxcAoN6eF7oAAID8CHMASABhDgAJ\nIMwBIAGEOQAkgDAHgARkCnMzO9vMvmVm95rZPWa2o//8i81sv5n90MyWzOz0cosLABjEsswzN7Ot\nkra6+51m9kJJ/y7pbZLeLem/3P1jZnaVpF9x96tLLTEA4ASZaubuftDd7+x//ZSk+yW9XNKlkm7o\n/9gN6gU8AKBiY/eZm9k5kl4r6buSznT3Q/3/OiTpzMJKBgDIbKww73exfE3Sle7+5Pr/815/DdcG\nAIAATs76g2b2fPWC/O/c/ab+04fMbKu7HzSzsyQ9PuD3CHgAmIC7W9afzTqbxSRdJ+k+d//Uuv+6\nWdLl/a8vl3TT8b/bLxAPd83OzgYvQywPtgXbgm2x+WNcWWvmb5L0Lkl3mdkd/ed2SvqIpK+Y2RWS\nHpJ02dglAADklinM3f1fNbwWf1FxxQEATIIVoBVqt9uhixANtsUatsUatsXkMi0ayvUCZl72awBA\nasxMXvQAKAAgboQ5ACSAMAeABBDmAJAAwhwAEkCYA0ACCHMASABhDgAJyHzVRACT63ZXtXfvkp55\n5mSdcsqz2rFjWtu2XRi6WEgIYQ6UrNtd1ZVX3qoDBz587LkDB3ZLEoGOwtDNApRs796lDUEuSQcO\nfFgLC/sDlQgpIsyBkj3zzOAG8JEjJ1VcEqSMMAdKdsopzw58fsuW5youCVJGmJeo211VpzOjdntO\nnc6Mut3V0EVCADt2TKvV2r3huVZrl7ZvvzhQiZAiBkBLwqAXjjq6vxcW9ujIkZO0Zctz2r79Eo4D\nFIrrmZek05nR0tK1A57fo337rglQIgB1wvXMI8GgF4AqEeYlYdALQJUI85Iw6AWgSvSZl6jbXdXC\nwv51g14XM+gFIJNx+8wJcwCIEAOgANBAhDkAJIAwB4AEEOYAkADCHAASQJgDQAIIcwBIAFdNTAD3\nl0TTcMyfiDCvOS61i6bhmB+Mbpaa4/6SaBqO+cEI85rjUrtoGo75wQjzmuNSu2gajvnBCPOa41K7\naBqO+cG4amICuNQummbUMZ/CbBcugQtMKIUAwODZLq3Wbs3Pd2q1P8cNc6YmAmK6W0qGz3bZk/S+\npM8cENPdUtLU2S6EOaDmBkCKmjrbhTAHVJ8A6HZX1enMqN2eU6czo253NXSRotPU2S6N7TNnsAvr\n7dgxrQMHdh83aLZL27dfErBUG9Gvn83RbbGwsGfdbJdLTthGyWWAu5f66L1EXBYXV7zV2uWSH3u0\nWrt8cXEldNEQ0OLiinc6Mz41Neudzkx0x8P09O4Nx+zRR6czE7potVOHDOhnZ+asbWTNvKmj3djc\ntm0XRr3/6dcvTooZ0Mg+cz4UqKO69OvXQYoZEH2YlzHgw4cCddTUgb0ypJgBmbpZzOx6SdskPe7u\nr+k/NyfpTyQ90f+xne6+r8jClTXgU4fBLuB4WQf2MFqKGZBpOb+ZvVnSU5K+uC7MZyU96e6fHPG7\nnuU1Bul0ZrS0dO2A5/do375rJvqbR3E9E6DZYs+AUpbzu/ttZnbOoNfL+kKTKLNfK/bBLgDlSi0D\n8vaZbzezH5jZdWZ2eiElWifFfi0AKEOeMP+MpHMlnS/pMUmfKKRE6zDgAwDZTDzP3N0fP/q1mX1e\n0i3DfnZubu7Y1+12W+12O9NrMOADoCmWl5e1vLw88e9nvp55v8/8lnUDoGe5+2P9r98v6fXu/ocD\nfm/iAVAAaKpSBkDN7EZJU5JeYmYPS5qV1Daz8yW5pAclvWeC8gIACsCdhgAgQuPWzKNfAQoAGI0w\nB4AEEOYAkIBGXgIX6UruhgNARoQ5ksGdeNBkdLMgGcNvOLA/UImA6hDmSEaKNxwAsqKbBcngwmyI\nSdXjN4Q5kpHiDQdQTyHGb1gBiqTEfsMBNEMRN9Yp5dosQOw2NmldH/jAW5ILcaZd1keI8RvCHLXX\nhCmJTXiPKQkxfsNsFtReE6YkNuE9piTEjXWomaP2mjAlsQnvMSUhbqxDmBeAvsywmjAlsQnvMTVV\n3zCaMM+JvsziTHpSTH1KYre7qieeOKgtW/5MR4585tjzKb1H5EeY5zS8L3MPYT6GPCfFlO8Vu7Zd\nPi9pVdIebdnyH3rVq35ZH/rQ25N4jygGYZ4TfZnFyHtSrLpJW5WN2+VCSRfqyBHppS+lsoCNmM2S\nE32ZxeCkOBjbBVkR5jmFmIKUIk6Kg7FdkFWwbpZUZoCk3F9bpdQHMSfFdkFWQa7NMmiwq9Xarfn5\nDiHYYFxXZTC2SzONe22WIGFexEVoACBl44Z5kD5zBnUAoFhBwpxBHQAoVpAwZwYIABQr2M0pGNQB\ngOFqMQCK9KUy9RQIhTsNITguPgZUjxWgKBw3UgCqR5ijcEw9BapHmKNwTD0FqkeYo3BMPQWqx2wW\nFOL42StvfOPL9J3vPMbUU2BCzGZB5YbNXuHCaUB16GZBbsxeAcIjzJEbs1eA8OhmQW7MXgEGq3Il\nNGEemToug+duOMCJql4JTZhHpK7L4Ll1HnCi4WNJewjz1FW984u0bduF0ZcRqFLVY0mEeUQYSCxP\nHbuvMLkY9nfVY0mEeUQYSCxHXbuvMJlY9nfVY0lJrgCN4aw8iUEHYau1S/Pz9D/nwQ3EmyWm/Z3n\nJjyNXwEay1l5EqEGEut68suK7qtmiWl/VzmWlFyY13kQUap+ILGKk1/okwXdV5sLvX+K1tj97e6l\nPnovUZ2pqVmX/ITH1NRspeWoi+np3QO3V6czU8jfX1xc8VZr14a/3Wrt8sXFlUL+/uRl2FlpGWIV\nw/4pWir7u5+dmbM2U83czK6XtE3S4+7+mv5zL5b0j5J+VdJDki5z95+VccIZR2PPyhMqu0kaQ0uJ\nefDDxbB/itbU/Z21m+ULkhYkfXHdc1dL2u/uHzOzq/rfX11w+cbGasTxlH3yi6X/knnwg8Wyf4rW\nxP2dKczd/TYzO+e4py+VNNX/+gZJy4ogzJt6Vp5U2Sc/WkpxY/+kI88A6Jnufqj/9SFJZxZQnkI0\n8aw8qbJPfkWfLFIbrAutqpYs+618hcxmcXc3s6GTyefm5o593W631W63i3hZFKTMk1+RJ4s6TzuN\nVRUtWfZbNsvLy1peXp749zMvGup3s9yybgD0AUltdz9oZmdJ+pa7//qA3/Osr5EFZ/jmimkxCLJj\nv02mykVDN0u6XNJH+//elONvZcIZfnwpnfxSHaxLHfutGlmnJt6o3mDnS8zsYUkflPQRSV8xsyvU\nn5pYViGPSnEaVZliOvkVcVJhsK6e2G/VyDqb5R1D/uuiAssyUhln+JRqrseL5eRX1EmFaaf1xH6r\nRq2W8xd9ho+p5lqGWJq3RZ1UmHZaT+y3atQqzIs+w8dScx1l0tZDLM3bIk8qTDutJ/Zb+WoV5kWf\n4WOpuW4mT+shluZtLCcVIGW1CnOp2DN8HUImT+shluZtLCcVIGW1CPOyBinrEDJ5Ww8xNG9jOakA\nKYs+zMscpKxDyNSh9ZBFDCcVIGWV3DZuamp24hp1iqvHxmlppHAruZSnfwJlifK2cSsrc5Imq1HX\nYZByHOO2NOrQetjMsPd7++336Nvf/gkBDxRlnDtZTPKQlOsONmXfCadqqb2fUYa931NPvSypu9uk\nYHFxxaend/vU1KxPT+9mfwSmMu40VKRxa9R1GKQcR2otjVGGvd+nnz5vw/cxzu9vktQX0DVB5WE+\n7sBd3bsZjpfKgGZWw96vdOL7TfWENkoMYwp1WUCH4SoN80lr1CnNhEitpTHKoPd76qnv0dNPv/OE\nn031hLaZWGrETWsxpqiSMJ+amqt9jbooqbU0Rhn0fi+44Lf0pS/dqgMH1t5zyie0zcRSI65TizGG\nlkyMKgnz5eW5Kl6mNlJqaWQx6P2+/vWrWljYo0ceeVwHD/5Mp556lvbuXTr2800RS424Li3GWFoy\nMYp+0RDSdPSDd+WVt+rw4c/q8GHpnnua98GMpUZclxZjLC2ZGBHmCIYPZlw14jq0GGNpycSIMBd9\ncOsN2haSStk+fDDrUyOORSwtmRg1Pszpg1szaFvcddcVkl6kgwc/eey5orYPH8yeOtSIYxFDSyba\nyt84K4wmefReIl5NW5G5mcHborzts7i44q3WruNWgu5k5SE2tbi44p3OjE9NzXqnM1Pp8TL4mC1n\n9bJiXwEaG5r6awZvi/K2D10M8Ym21rlOyJZMzOM8jQ9zmvprBm+LcrcPXQzxoMtxtJgrf88LXYDQ\nduyYVqu1e8NzvT64iwOVKJxB22Lr1p9o69a/3PBcU7dP6obXOvcHKlF8Yq78Nb5mTlN/zeBt8ccD\nnmvm9kldzLXOWMQwADtMJTenKPs1AOSX4o1gytDtrmphYf+6ys3FpVRuxr05BWEOQFIad7VKCWEO\nYGJV1ToxGmEOAAkYN8wbP5sFAFJAmANAAghzAEgAYQ4ACWj8oiHErQ7XCgFiQJgjWlwrBMiOqYmI\nFisSkVWKLbhxpyZSM2+Auh7oXCsEWdCC6yHME1fnAz3mK9QhHjFfY7xKzGZJXJ0va8rliZFFlS24\nbndVnc6M2u05dToz6nZXC3+NSVEzT1yduyq4PHHaiur+q6oFF3srlzBPXN27KrgTUZqKDMaqrjEe\ne3cOYZ64mC+mj/LEPuhdZDBW1YKLvZXb+DCP/aDPi66K5om9O0AaPxhHfU6raMFF38p191IfvZeI\n0+Liirdau1zyY49Wa5cvLq6ELhowsenp3RuO6aOPTmcmdNGOGaeMsXxOB5djZ2nl6Gdn5qxtdM08\na1Mv9do70hJ7d4A0XvdfLH3VsbdyGx3mWQ76OjRZgfWi7w7QeMEY08kp5gH5Rod5loM+lloB0lZk\n668ug95Zg7EOJ6cYNDrMsxz0MdUKkKaiW3+xdweMqy4np9AaHeZZDnpqBc0SYnykjNZfzN0B40rt\n5FSW3GFuZg9J+rmk5yT9wt3fkPdvVmnUQU+toDlCjY/Q+hstpZNTWYqombuktrv/tIC/FR1qBc0R\nanyE1h+KUFQ3S+Zr7tYRtYJmCFVDpvWHIhRVM/+GmT0n6bPu/rkC/iZQuVA1ZFp/KEIRYf4md3/M\nzF4qab+ZPeDut63/gbm5uWNft9tttdvtAl4WKFbIGjKtPywvL2t5eXni3y/0tnFmNivpKXf/xLrn\nvMjXAMrU7a5qYWH/uhryxdGELCuR05B1P1Z62zgzO03SSe7+pJn9kqRpSX+V528CIcVaQ2YlchrK\n3I957zR0pqTbzOxOSd+VtOjuSzn/JoDj1PmOUVhT5n7MVTN39wclnZ+7FEBJUumaYC56Gsrcj41e\nAYrixRSeKXVNMBc9DaXux3GulzvJQxFfzxzFiuW600fV4breWVV9LW2UY5z9KK5njlBiu8Jknbsm\nBrVw5uc7zEWvuTLXFBDmKExs4VnXrolh3UPz8x3t23dNwJKtiak7rW7KmjFFmKMwsYVnXZfJh27h\njArqlMYiUkKY11SMNaPYwrOuy+RDtnCyBHXokw0GSy7MYwy5osVaM4oxPGNdBLSZkC2cLEFdxMmm\nCZ/TqiUV5rGGXNFirhlVHZ4phkLIFk6WoM57smnK57RqSYV5zCFXpNgGGkNJNRRCtnCyBHXek01d\nPqd1qyjUIsyzbtSsIVe3nXS82AYaQ9ksFI7+f133cajuoSxBnfdkU4fKSB0rCtGH+TgbNUvI1XEn\nHS+2gcai5T15P/rok7Xfx6FkDeo8J5s6VEbq0nrYYJwVRpM8lHMF6Dir+LKsrkplVeDi4op3OjM+\nNTXrnc5MMisBx1lFOmxfnnHGZRPv48XFFZ+e3u1TU7M+Pb07me0akzqsZp2amh14DE1NzVZWBsW4\nAjRPt8Y4TbIstYo6NPGyqOMsjSzGqRENa6Fs2XK6Dh8+8W+P2scptNrqIMZZT8erQ+vheJWEeZ4P\nyLgbdVTI1XEnNUkRJ++9e5d0770n/o1R+7iWTeuair0ykqUrM7axt0rCPM8HpOj+4dT7m+uuqJP3\nJPu4qFZbbB9yjG9U6yHGVlywAdCsH5Cim2R1aOI1WREn20n3cRGtthg/5JjMZq2HSVtxpZ7ox+lg\nn+QhqZIBRwau0hFqcLeIgblUBtixuUkGSMe9RLRiHABttcrt1qA2lJZQ/alFtNpSGWDH5iZpxZU9\nJlNJmB9/HeYLLniF9u5d0sc//s1CmhoMXKEoeU8kDLA3wyTdgWWf6CsJ8/UfkDJq0dSGEAsG2MOq\navB5klZc2Sf6ygdAy6hFUxuCFMcsEgbYw6m6u3XcVlzZJ/rKw7yMWjS1IcQ0bhL7HOpUxd7dWvaJ\nvvIwL6MWTW0IsX+QUb46dLeWeaKvPMzLqkVTG2q2OnyQUa6md7dWHubUolGGpn+QQXer9eaml/gC\nZl72awCD+sxbrV2an6ei0CTd7qoWFvavqyheXNv9b2Zyd8v884Q5UpH3gxzDbBjgKMIcmMDgmv1u\nzc93CHQEQZgjKnWp7XY6M1paunbA83u0b981AUqE2JV9bI8b5tHfNg71FdPc71GYDYNxxHhsPy/I\nq6IRhs/93h+oRMMxGwbjiPHYJsxRmjrVdnfsmFartXvDc71pbRcHKhFiFuOxTTcLSlOn2i7rHzCO\nGI9twnwMdRnMy6KK91K3RRysIkZWMR7bhHlGMQ54TKqq90JtF6mK8dhmauIm1tde77nnfh0+/F5J\nG3dWHaeuMQ0PiB9TEwsyqPYqHR0gWwv0GAfzRolx8CaLlLq5gKIR5kMMmnokfVjSHq0P8xgH80aJ\ncfBmlJS6uYAyMDVxiGG1V2mt9lrXqWt1nIYX47xeICbUzIcYVns944wH9OpXz0Ux4DGpGAdvRqlr\n1xBQFcJ8iGFTj+bn/zzq0MuqbtPw6tg1BFSJMB+ijrXXlMU4rxeICVMTURsp3XgAGIVL4AJAAphn\nDiA41gRUjzAHUCjWBISRe565mV1iZg+Y2Y/M7KoiCgWgvlgTEEaumrmZnSTp05IukvSopNvN7GZ3\nv7+IwgGoH9YEhOlmytvN8gZJP3b3hyTJzL4s6a2SCHOgoZq+JiBUN1PebpaXS3p43feP9J8DELlu\nd1Wdzoza7Tl1OjPqdlcL+bt1vFxEkUJ1M+WtmTPnEKihMmuPTV9wF6qbKW+YPyrp7HXfn61e7XyD\nubm5Y1+322212+2cL1seplShCYbXHvcUcrzX7XIRRZq0m2l5eVnLy8uTv7C7T/xQ72RwQNI5kl4g\n6U5J5x33M14Xi4sr3mrtcsmPPVqtXb64uBK6aEChpqZmNxznRx9TU7Ohi1Z7g3Nk59g50s/OzHmc\nq2bu7s+a2fsk3aretWGv8xrPZCm7tgLEoumDlGUK1c2Ue9GQu39d0tcLKEtwTKlCU3DhsnKF6GZi\nBeg61FbQFE0fpEwRF9paZ9AIf+8a5hzkwKSYVDAZLrSVA7UVoFhcp6U61MwBlKbTmdHS0rUDnt+j\nffuuCVCiNbG3GKiZA4hGrJMKUmwx5L5qIgAME+ukghSv7EiYAyhNrNdpibXFkAfdLABKE+ukglhb\nDHkwAAqgceowDZkbOgNABt3uqhYW9q9rMVwcTZBLhDkAJKGWUxNjn+85iRTfE4B4BQ/zFOd7pvie\nAMQt+NTEFOd7pvieAMQteJinON8zxfcEIG7BwzzF+Z4pvicAcQse5rGuEMsjxfcEIG5RTE2Mfb7n\nJFJ8TwCqwzxzAEjAuGEevJsFAJAfYQ4ACSDMASABhDkAJIAwB4AEEOYAkADCHAASQJgDQAIIcwBI\nAGEOAAkgzAEgAYQ5ACSAMAeABBDmAJAAwhwAEkCYA0ACCHMASABhDgAJIMwBIAGEOQAkgDAHgAQQ\n5gCQAMIcABJAmANAAghzAEgAYQ4ACSDMASABhDkAJGDiMDezOTN7xMzu6D8uKbJgAIDs8tTMXdIn\n3f21/ce+ogqVquXl5dBFiAbbYg3bYg3bYnJ5u1mskFI0BAfqGrbFGrbFGrbF5PKG+XYz+4GZXWdm\npxdSIgDA2DYNczPbb2Z3D3hcKukzks6VdL6kxyR9ooLyAgAGMHfP/0fMzpF0i7u/ZsD/5X8BAGgg\nd8/clX3ypC9iZme5+2P9b39f0t15CwMAmMzEYS7po2Z2vnqzWh6U9J5iigQAGFch3SwAgLBKWwFq\nZpeY2QNm9iMzu6qs16kDMzvbzL5lZvea2T1mtiN0mUIys5P6C81uCV2WkMzsdDP7qpndb2b3mdkF\nocsUipnt7H8+7jazfzCzU0KXqSpmdr2ZHTKzu9c99+L+BJQfmtlSltmCpYS5mZ0k6dOSLpH0Kknv\nMLPzynitmviFpPe7+29IukDSexu+Pa6UdJ96XXRNNi/pn939PEm/Ken+wOUJoj+B4k8lva4/ieIk\nSX8QskwV+4J6Wbne1ZL2u/srJf1L//tNlVUzf4OkH7v7Q+7+C0lflvTWkl4reu5+0N3v7H/9lHof\n2peFLVUYZvYKSb8r6fNq8KIzM3uRpDe7+/WS5O7Puvv/BC5WKD9Xr8JzmpmdLOk0SY+GLVJ13P02\nSf993NOXSrqh//UNkt426u+UFeYvl/Twuu8f6T/XeP1ayGslfTdsSYL5a0kfkPR/oQsS2LmSnjCz\nL5jZ983sc2Z2WuhCheDuP1Vvncp/SvqJpJ+5+zfCliq4M939UP/rQ5LOHPULZYV505vPA5nZCyV9\nVdKV/Rp6o5jZ70l63N3vUINr5X0nS3qdpL9x99dJ+l9laEqnyMxakv5C0jnqtVhfaGbvDFqoiHhv\nlsrITC0rzB+VdPa6789Wr3beWGb2fElfk/Qld78pdHkC+R1Jl5rZg5JulPQWM/ti4DKF8oikR9z9\n9v73X1Uv3JvotyX9m7sfdvdnJf2TesdKkx0ys61Sb02PpMdH/UJZYf49Sb9mZueY2QskvV3SzSW9\nVvTMzCRdJ+k+d/9U6PKE4u673P1sdz9XvQGub7r7H4UuVwjuflDSw2b2yv5TF0m6N2CRQnpA0gVm\ndmr/s3KRegPkTXazpMv7X18uaWQFMM+ioaHc/Vkze5+kW9Ubmb7O3Rs5Ut/3JknvknSXmd3Rf24n\nlw1ufHfcdkl/36/wHJD07sDlCcLdf9BvoX1PvbGU70v627Clqo6Z3ShpStJLzOxhSR+U9BFJXzGz\nKyQ9JOmykX+HRUMAUH/cNg4AEkCYA0ACCHMASABhDgAJIMwBIAGEOQAkgDAHgAQQ5gCQgP8HItGf\nMcNCLgsAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c29ce10>"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now the correlation is just .04, with a large probability (p-value) that the small non-zero value is a statistical fluctuation. One way to view the p-value for the original data is as the probability that one would get that large a correlation coefficient by chance if the variables were actually uncorrelated. Equivalently, that means the probability that the correlation coefficient would be as large as measured after a random permutation of the y values."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Logistic Regression**\n",
      "\n",
      "The last thing discussed in lecture was a modification of the regression scheme when trying to predict a binary outcome rather than a linear relation. This makes use of the so-called \"logistic function\":\n",
      "\n",
      "$${\\rm logit}(x) = {1\\over 1+{\\rm e}^{-x}}$$\n",
      "\n",
      "For large positive values of the argument $x$, ${\\rm e}^{-x}$ becomes very small so the function evaluates to 1, whereas for large negative values of $x$, ${\\rm e}^{-x}$ becomes very large and the function evaluates to zero. The value for $x=0$ is logit(0)=.5, and that marks the middle of the transition region between 0 and 1. \n",
      "\n",
      "$x$ can be considered to be any linear combination of feature values $f_i$ being used to make the binary prediction, $x=\\alpha_1 f_1 + \\ldots + \\alpha_n f_n + \\beta$, where the weights $\\alpha_i$ and $\\beta$ are determined by some training procedure (to provide the best fit to labelled training data). Intermediate values of the logit function can be interpreted as the probability of the label being 1 for those values of the features.\n",
      "\n",
      "For a single feature $x=\\alpha f_1+\\beta$, the logit is equal to .5 for $f_1=-\\beta/\\alpha$, and the width of the transition region is proportional to $1/\\alpha$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def logit(x): return 1/(1+exp(-x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f1=arange(-10,10,.1)\n",
      "plot(f1,logit(f1),label='logit($f_1$)')\n",
      "plot(f1,logit(10*f1-30),label='logit($10f_1-30$)')\n",
      "yticks(arange(0,1.1,.1))\n",
      "grid('on')\n",
      "legend(loc='upper left');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW5+PHPQ9jLknhRkIjiBrJUg/QGtyq4RgEpbRXR\nUhFt8fqDulyugFbF6i2VukDFerlWjKUiWkC2ilSRgLSCrF6WgKCsgVJAAgkEksDz++NM0kmYZGaS\nmZxzZp736zUvcs6cnHnmy8mTb575nu9XVBVjjDGJoZ7bARhjjIkdS+rGGJNALKkbY0wCsaRujDEJ\nxJK6McYkEEvqxhiTQMImdRGZLCL7RGRdNcf8TkS2iMiXItIttiEaY4yJVCQ99beArKqeFJHbgItU\n9WLg58DrMYrNGGNMlMImdVX9DDhUzSG3A28Hjl0OpIpI69iEZ4wxJhqxqKmnA7uCtncD58TgvMYY\nY6IUqw9KpdK2zT1gjDEuqB+Dc+QB7YK2zwnsq0BELNEbY0wNqGrljnOVYpHU5wDDgGkicgWQr6r7\nqggsBi9nAMaMGcOYMWPcDiMhJENbHjlxhO++/l2euvYphnQbQj2J/I90VfjnP2HDBti4Eb75BrZv\nh23bnH9LSqBtW2jd2nns2DGGPn3G0Lo1nHkmpKZCy5bQooXzb8uW0KgRSMRpKrlJlA0VNqmLyLvA\ndUArEdkFPAM0AFDVSar6oYjcJiJbgaPAfVFHbaK2fft2t0NIGMnQlv/11//i5gtu5oHLH6j2OFUn\nWf/977BsGaxb5yRzVejSBTp3hosugquvhvbtnccZZ1RM0IMHb+eZZ+L6dkw1wiZ1VR0YwTHDYhOO\nMSbW1uxdw1+2/IUND2047TlVJ2kvWABLl8Lnn0NKClx1FVx5JfTv7yTz1q2tZ+0XsSi/GBcMHjzY\n7RASRqK35YdbPuSOznfQsnFLAEpL4dNPYcYMmD/fSeJZWTBgAEyYAO3a1S6BJ3p7ep3UVZ1bRNRq\n6sbUvRv/eCO/yHyY1of78sc/wvTpTtnkzjuhd2/o2NF64V4mInX+QWmtRPshgEkMXvoFn5OTQ8+e\nPd0OIy4OHTnB0u3LOfDatRQehCFDnBLLBRfE7zUTuT39wPWkDt76ATfxZ7/I46+gACZOhBemLafx\nLZcwdkxLbrkF6tkUfgnP9fJL4E+LOonBeIP9n8fPsWPwu9/Byy/DTTdBWv9nadryKONuGud2aKaG\noi2/2O9tYxKAKsya5Qw5XL0aFi+Gd96B9YWLuP78690Oz9QhS+om6eXk5LgdQq188w306QNPPAFv\nvQXvvw+dOjnPrd67msz0zDqNx+/t6XeW1I3xKVV4+23o0QOuvRbWroVevf71fMGJAk7qSdIap7kX\npKlzltTDaN++PQsXLqz1ebp27cqSJUuqPWb06NFMmDChRufv0aMHGzdurNH3Jjs/jtTIz4e774bf\n/tYZcz5yJDRsWPGYvII80pun1/kH035sz0RiST0MEYnJD8X69eu59tprAecXxaefflrh+f379zNl\nyhQefPDBGp1/xIgRPP3007WO03jfhg1w+eXO7fkrVsB3vxv6uLwjebRt3rZugzOus6TuglCjP7Kz\ns+nduzeNGjUK+T2vvPIKTz75JJMmTQr5fN++fVm0aBH79oWcS81Uw0814E8+cUosv/oVvPYaNGlS\n9bF5BXmkt0ivu+AC/NSeiciSehRyc3Pp2bMnaWlpdO3alblz55Y/t3r1arp160aLFi248847GTBg\nAE899VT582VlnEGDBrFz50769u1L8+bNefHFFwGYP38+1113XcjXPXz4MO+//z79+vXjmmuuCXlM\n48aN6d69OwsWLIjhOzZe8oc/wE9+4twR+pOfhD9+T8Ee0pvXfVI37rKkHqGSkhL69u1LVlYW+/fv\n59VXX+Wee+5hy5YtFBcX079/f4YMGcKhQ4cYOHAgs2bNqlC2KSvjTJkyhXPPPZd58+ZRUFDAiBEj\nAKc807Fjx5CvvXz5cjIyMsjMzKRLly5VxtipUye+/PLL2L7xJOCHGvBvfuM8lixxPhSNRN6RPFeS\nuh/aM5FZUo/QsmXLOHr0KKNGjaJ+/fr06tWLPn36MHXqVJYtW8bJkycZPnw4KSkp9O/fn8zM6IaR\n5efn07x589P2L1++nAkTJlBaWsoHH3xQ7TmaN29Ofn5+VK9rvO+3v4XJk52E3qFD5N+XV2A19WQU\nNqmLSJaIbBKRLSIyMsTzaSLygYh8KSLLRaTqrmQNiMTmUVt79uyhXbt2Ffadd9555OXlsXfvXtLT\nK/aI2rVrF9Vdk2lpaRQUFJy2v0ePHjRp0oRHHnmE/v37V3uOI0eOkJZmw9ei5eUa8CuvwKRJsGiR\nsxBFNKymnpyqTeoikgJMBLKAzsBAEelU6bAngNWqehnwU6BmY/KqoBqbR22lp6eza9euCol6x44d\nnHPOOZx99tnk5VVcwW/nzp1VjpoJtf/SSy9l8+bNIY/Pzc2lc+fOHD58mJkzZzJ27Ngqj7vssssi\nfUvG437/e3j1VSehp9cgN1tNPTmF66lnAltVdbuqlgDTgH6VjukELAJQ1c1AexE5M+aRuqxHjx40\nbdqUcePGUVJSQk5ODvPmzeOuu+7iyiuvJCUlhYkTJ1JaWsrs2bNZsWJFledq3bo1X3/9dYV9t912\nG4sXLz7t2H379tGqVStEhJYtW9K9e3eKi4tPO+748eOsXr2am266qfZvNsl4sQY8fz48/zwsXOjM\nbx6tk6dOsq9wH2c3Pzv2wYXhxfZMJuGSejqwK2h7d2BfsC+BHwKISCZwHs7i0wmlQYMGzJ07l/nz\n53PmmWcybNgwpkyZQocOHWjQoAEzZ87kzTffJC0tjXfeeYc+ffrQsPLdIAGjR4/m+eefJy0tjZdf\nfhmAQYMG8eGHH3L8+PEKxy5fvpyrr746bHxz586lV69etGnTpvZv1rhq40a4915nlMv559fsHP88\n+k9SG6fSMCX0NWgSV7ipdyMpXPwGmCAia4B1wBrgZKgDBw8eTPv27QFITU0lIyMj8khdsm3btvKv\nO3fuXGW9sHv37qxZs6Z8u0ePHtx+++0hz3P77bdXeA6gVatW/PSnP2XSpEk8/PDDrFq1ijfeeIMz\nzjiDAQMGhI3zpZdeYvLkyZG+LdeVtWNZr87N7eD/U7fj6dq1J337wv335+D8QVaz881ZMIcWe1uU\nv69kbU8/bufk5JCdnQ1Qni+joqpVPoArgI+CtkcDI8N8zzagWYj9GkpV+/1m8eLFunfvXi0pKdHs\n7Gxt2rSp/uMf/6jx+dasWaPPPfecvvTSSxX2b9u2TceMGVPbcF3ltf/zRYsWuR2CqqqWlKhed53q\nyJG1P9fsTbO19zu9a3+iGvBKeyaKwM9Ltbk6+BGup74SuFhE2gN7gAFAhYWoRaQlUKSqxSLyM2Cx\nqhZG/+vF3zZv3sydd97J0aNHufDCC5k+fTqtW7eu8fkyMjJO+0umsLCQGTNmsGrVKtavX0/Xrl1r\nG7bBOzXg556DRo3g17+u/bncnCLAK+2ZrMIukiEitwLjgRTgTVUdKyJDAVR1kohcCWTjlGrWA/er\n6uEQ59FQr2ULJiQf+z8/3dKlcMcdzlzoZ8fgs81ffvpLGtRrwDM9n6n9yYyrYr5IhqrOV9WOqnqR\nqo4N7JukqpMCX38eeP4SVf1xqIRujJe5Pa46Px8GDYL//d/YJHRwPiht3azmfynWhtvtmezsjlJj\nXKQKDz0Et94KffvG7rwFxQU0b3j6Hcom8Xli4Wlj3ORmDfi995zFLVaujO15C04U0LyRO0ndauru\nsqRujEvy8+Gxx2DGDGjaNLbntp568rLyi0l6btWAf/lLp+Ry5ZWxP3dhcaFrPXWrqbvLeurGuGDF\nCqeHvmFDfM5fcMJ66snKeuom6dV1DfjkSXjwQRg3zlmSLh4Kigto1rBZfE4ehtXU3WVJ3Zg69vvf\nQ4sWka1eVFNuflBq3GVJ3SS9uqwB5+c7d45OnBibef5DOaWnOFZyzLWeutXU3WVJPYyytUVrq2vX\nrixZsqTaY0aPHs2ECTGdjr7O9ejRg40bN7odhme98ALcfjtUsyphrR0tPkrTBk2pJ/bjnYzsg9Iw\nytYWra3169eXf92+fXsmT57M9ddfX75v//79TJky5bR51teuXcuf/vSn8gWqAWbNmsXGjRupV68e\n6enpDBo0qNbxVWX27NkUFhby9ddf06pVKx566KFqYxgxYgRPP/0006dPj1tMsVZXNeA9e5y7RuO9\njGxBsbulF6upu8uSugtCzX2SnZ1N7969adSoUfm+l19+maVLl9KyZcvyfYcPH+a5555j1apVAFx5\n5ZXceuuttGrVqsL5XnnlFQ4cOMC5557L0KFDaxRnfn4+AwYMID8/n0aNGtGqVSt69+5NampqlTH0\n7duXBx98kH379tVqQrNE9OyzcP/9cE6cVxuwkS/Jzf4+i0Jubi49e/YkLS2Nrl27Mnfu3PLnVq9e\nTbdu3WjRogV33nknAwYM4Kmnnip/vqyMM2jQIHbu3Enfvn1p3rx5eQ98/vz5XHfddRVe77HHHqNf\nv4oLTS1ZsoTOnTuXb1922WUsWrSowjGHDx/m/fffp1+/flxzzTU1fr+pqamsWrWKxo0bIyKUlpai\nqtXG0LhxY7p3786CBQtq/Lp1rS5qwJs3w8yZMGpU3F/K1ZEvYDV1t1lSj1BJSQl9+/YlKyuL/fv3\n8+qrr3LPPfewZcsWiouL6d+/P0OGDOHQoUMMHDiQWbNmVSjblJVxpkyZwrnnnsu8efMoKChgxIgR\ngFOe6dix42mvW7lHv3v3blJTU8u3U1NT2bJlS4Vjli9fTkZGBpmZmXSpZfG27PuXLl1Kz549ad++\nfdgYOnXqxJfxrjH4zJNPwogR8RvCGMxGviS3sEldRLJEZJOIbBGRkSGebyUiH4nIWhFZLyKD4xKp\ny5YtW8bRo0cZNWoU9evXp1evXvTp04epU6eybNkyTp48yfDhw0lJSaF///5kZmZGdf78/HyaNz/9\nB7FyPT8/P5/GjRuXbzds2JDCwn9NX798+XImTJhAaWkpH3zwQZTvMrSZM2fy2muv8dJLL0UUQ/Pm\nzcnPz4/Ja9eFeNeA166Fv/8dhg+P68uUKywudLX8YjV1d1VbUxeRFGAicCOQB6wQkTmqmht02DBg\njaqOFpFWwGYR+ZOqlsYiQHk2NuO+9Jnazd+9Z88e2lVaAfi8884jLy+PvXv3kl5pufd27dpFNWd4\nWloaBQUFp+2vfI7mzZtz8ODB8u2ioqIKtesePXrQpEkTHnnkkWp76ePGjaOoqCjkc/fee2+FZbR+\n+MMfcvPNN9OtWzc+/vjjsDEcOXKEtLS0qt9skvnNb5w5XmI9v0tV3P6g1Lgr3AelmcBWVd0OICLT\ngH5AcFLfC1wa+LoFcDBWCR1qn4xjJT09nV27dqGq5b3nHTt2cMkll3D22WeTl5dX4fidO3dy0UUX\nhTxXqNE0l156KZs3b6Z79+7VHnvhhReyMmhKvwMHDnD55ZdXOCY3N5fOnTtz+PBhFi5cyObNmxk9\nenSFYx5//PEw7xj+8pe/8Otf/5q//e1vNGvWjLPOOovp06fTpUuXamPIzc3lpz/9adjze0VOTk7c\nepdbt8LChfDGG3E5fUhuf1Aaz/Y04YUrv6QDu4K2dwf2BXsD6CIie4AvgYdjF5539OjRg6ZNmzJu\n3DhKSkrIyclh3rx53HXXXVx55ZWkpKQwceJESktLmT17NitWrKjyXK1btz5t6OJtt93G4sWLTzu2\nck/92muvLR91As4HtDfccEP59r59+2jVqhUiQsuWLenevTvFzgrGUUtJSSn/4VRVdu3axaWXXsp1\n111XZQzHjx9n9erV3HTTTTV6zUQzbhz8x39AiMpa3NgMjcktXE89km7yE8BaVe0pIhcCH4vIZap6\nWi1h8ODB5X/Wp6amnrYGp5c1aNCAuXPn8tBDDzF27FjOOeccpkyZQocOHQCn7vzAAw8wevRobr31\nVvr06UPDhg1Dnmv06NEMHz6cxx9/nKeeeorHHnuMQYMG0a1bN44fP15er544cSLvv/8+u3bt4tln\nn+XRRx+lRYsWPP744zz//POcOnWKxx9/nLPOOqv83MuXL+fqq6+OyXvOysrim2++4dVXX2XHjh08\n+eST3HzzzQBVxjB37lx69epFmzZtqj23l1ZvL1vBPdbnP3AApk/vyVdf1e37KThRwMHcg+Q0yqmT\n16u8Ha/2TJbtnJwcsrOzASqUQSNW3arUwBXAR0Hbo4GRlY75ELg6aHsh8L0Q56pupeyEk5mZqdnZ\n2VF9zxNPPKHjx4+v0eutXLlShw4dqqNHj9a1a9eW79++fbuOGTOmRuesiR49euiGDRuqPSZR/88r\n+8//VP3FL+r+dR/76DH97d9+W/cvbOIi8PNSba4OfoTrqa8ELhaR9sAeYAAwsNIxm3A+SP2biLQG\nOgLfRP/rxd+WLFlChw4daNWqFe+88w7r168nKysrqnP893//d41fPyUlhXPOOYemTZty2WWXle/X\nOl7gedmyZXX6erEQjxrwt9/C5MnOyJe6VlBcQIeGHer+hQOspu6uapO6qpaKyDBgAZACvKmquSIy\nNPD8JODXwFsi8iVOjf5xVf02znF7zubNm7nzzjs5evQoF154IdOnT6/TOyozMjJOK2cVFhYyY8YM\nVq1axfr16+natWudxZPsJk1y5ng599y6f20b/ZLcpK56ciKioV4r1C3zJrEl+v95aSlccAHMmgWV\nBibViT5T+zC0+1D6dozhStbGNYGfl4jHdtsdpcbE2Ny50K6dOwkdrKee7Cypm6QX67lKJk6EYcNi\nesqouH1Hqc394i5L6sbE0MaNzuNHP3IvhoIT7k7oZdxlSd0kvViO1HjtNfj5z6GKWxTqhNvlFxv5\n4i6bT92YGDl8GN59F4LWQ3GF29MEGHdZT90kvVjVgP/4R7jpJmjbNianq5FTeoqi0iK+0/A7rsVg\nNXV3eaKnHovl4oxxkyr8z//A66+7G4etT2pcT+qJPF7Z+EMsasBffAElJfD979c+ntrwwmReVlN3\nl/06NyYGJk+G++4Dt//otJEvxpK6T1ndMnZq25bHjsGf/wxemELe7ZEvYNem2yypG1NLM2fCFVdA\neuWVBlxQWFxoPfUkZ0ndp6xuGTu1bcvJk2HIkNjEUltFJUU0bVBH6+ZVwa5Nd1lSN6YWvvkG1q2D\nvh6ZO+t46XGa1G/idhjGRWGTuohkicgmEdkiIiNDPD9CRNYEHutEpFREUuMTriljdcvYqU1bvv02\n3H03NGoUu3hqo6i0iMb1G7sag12b7qo2qYtICjARyAI6AwNFpFPwMar6oqp2U9VuOCsj5ahqfrwC\nNsYrTp2C7Gxn1ItXHC897npSN+4K11PPBLaq6nZVLQGmAf2qOf5u4N1YBWeqZnXL2KlpWy5ZAqmp\n4KWldr1QfrFr013hkno6sCtoe3dg32lEpClwCzAjNqEZ421Tp8I997gdRUVFJe6XX4y7wiX1aG73\n7AsstdJL3bC6ZezUpC1PnIAZM+Cuu2IfT214ofxi16a7wk0TkAe0C9puh9NbD+UuwpReBg8eTPv2\n7QFITU0lIyOj/E+1sgvBtiPbXhtY0dgr8STb9osv5pCeDuee6414yraPnzpOkwZNPBOPbUe/nZOT\nQ3Z2NkB5voxGtWuUikh9YDNwA7AH+AIYqKq5lY5rCXwDnKOqRVWcK+Qapcb40V13Qc+e8OCDbkdS\n0Yi/jqBNszaMuGqE26GYGInpGqWqWgoMAxYAG4H3VDVXRIaKyNCgQ38ALKgqoRuTSAoLYf58+PGP\n3Y7kdF4ovxh3hR2nrqrzVbWjql6kqmMD+yap6qSgY95W1bvjGaipqOzPNVN70bbl7NlwzTXQqlV8\n4qkNL4x+sWvTXXZHqTFRmjrVueHIi7xw85FxV7U19Zi+kNXUTQLYvx8uugjy8qCZB+fN+tH7P+Lu\nrnfzo84urnxtYiqmNXVjTEUzZsCtt3ozoYPV1I0ldd+yumXsRNOW06fDgAHxi6W2ikqKaNLAaurJ\nzJK6MRHavx9WroSsLLcjqZr11I0ldZ8qu2nB1F6kbTl7NtxyCzTx8My2Xkjqdm26y5K6MRGaPt2b\nY9ODFZUWuT6k0bjLkrpPWd0ydiJpy2+/hc8/dz4k9TIv9NTt2nSXJXVjIjBnDtx4o3dHvZTxQlI3\n7rJx6sZEoE8fZ5rdgQPdjqR6aS+kse3hbaQ2tsXHEoWNUzcmxg4fhs8+g9693Y4kPJtP3VhS9ymr\nW8ZOuLacO9eZkbFFizoJp8ZUleKTxTRKcXfBVLs23WVJ3Zgw/DDqBeDEyRM0qt8IkYj/UjcJKGxN\nXUSygPFACvAHVX0hxDE9gVeABsABVe0Z4hirqRvfOXYM2rSBHTsgLc3taKp3qOgQF/zuAg6NPOR2\nKCaGoq2pV7vykYikABOBG3FWQVohInOCF8kQkVTgNeAWVd0tIh6ckNSYmvnkE/je97yf0MFGvhhH\nuPJLJrBVVberagkwDehX6Zi7gRmquhtAVQ/EPkxTmdUtY6e6tpwzB26/ve5iqQ2vJHW7Nt0VLqmn\nA7uCtncH9gW7GDhDRBaJyEoRGRTLAI1xy6lTzoekffu6HUlk7G5SA+EXno6kCN4AuBxnHdOmwOci\nskxVt9Q2OFM1m18jdqpqyy++gDPPhAsvrNt4asorPXW7Nt0VLqnnAe2Cttvh9NaD7cL5cLQIKBKR\nJcBlwGlJffDgweWrY6emppKRkeGJ1btt27ZDbb/xBtx+u3fiCbe9/p/ry5O6F+Kx7Zpt5+TkkJ2d\nDVCeL6OiqlU+cJL+10B7oCGwFuhU6ZhLgE9wRsc0BdYBnUOcS03sLFq0yO0QEkZVbdmli+rnn9dt\nLLXxydef6PVvX+92GHZtxlggd1abq4Mf1fbUVbVURIYBCwJJ+01VzRWRoYHnJ6nqJhH5CPg/4BTw\nhqpujP7XizHe8fXXcOAAZGa6HUnkvFJ+Me6yuV+MCWH8eFi/Hv7wB7cjidyMjTOYun4qM+6c4XYo\nJoZs7hdjYsBPQxnL2OgXA5bUfavsgxVTe5Xb8tAhZ9m6G290J56a8kr5xa5Nd1lSN6aS+fOdCbya\nNnU7kuh4Jakbd1lS96myoVCm9iq3pR9LL+BMu+uF8otdm+6ypG5MkOJiWLDAWRTDb6ynbsCSum9Z\n3TJ2gttyyRLo2NGZmdFvvJLU7dp0lyV1Y4L4tfQCzugXLyR1465w0wQYj7K6ZeyUtaWqk9TnzXM3\nnpo6XnqcJg2spp7srKduTMC6dVCvHnTp4nYkNeOV8otxlyV1n7K6ZeyUtWVZ6cWvq8F5Janbteku\nS+rGBPi5ng52R6lxWFL3Katbxk7Pnj3Zswe2boXvf9/taGrOKz11uzbdZUndGJwPR7OyoEEDtyOp\nOa8kdeMuS+o+ZXXL2MnJyfF96QUCd5R6YPSLXZvuCpvURSRLRDaJyBYRGRni+Z4iclhE1gQev4xP\nqMbER1GRc9NRVpbbkdSO9dQNhBmnLiIpwETgRpyl7VaIyBxVza106GJV9Xk/x1+sbhk7J070JDMT\nUlPdjqR2vJLU7dp0V7ieeiawVVW3q2oJMA3oF+I4nw4CM8b/o17K2OgXA+GTejrOwtJldgf2BVPg\nKhH5UkQ+FJHOsQzQhGZ1y9g4eRJmzsyhb1+3I6k9r/TU7dp0V7hpAiJZf2410E5Vj4nIrcAsoEOo\nAwcPHly+OnZqaioZGRmeWL3bj9tr1671VDx+3W7YsCdnnAE7duSwY4f78dRmu/CrwvKk7oV4bLtm\n2zk5OWRnZwOU58toVLtGqYhcAYxR1azA9mjglKq+UM33bAO6q+q3lfbbGqXGc0aNgvr14fnn3Y6k\n9ho+15CC0QU0qt/I7VBMDMV6jdKVwMUi0l5EGgIDgDmVXrC1iHNjtYhk4vyi+Pb0UxnjPYlSTz95\n6iSlp0ppmNLQ7VCMy6pN6qpaCgwDFgAbgfdUNVdEhorI0MBhPwbWichaYDxwVzwDNo6yP9dMzW3Z\n4qxHWliY43YotXbi5Aka12+MeGDiGrs23RV26l1VnQ/Mr7RvUtDXrwGvxT40Y+Jr7lzo29eZmdHv\nikpsLnXjSIDLOTmVfcBiaq6s9JIIbemVudQhMdrTzyypm6R08CCsXg033OB2JLHhleGMxn2W1H3K\n6pa1M38+XH89NGmSGG3ppaXsEqE9/cySuklKiTLqpczx0uN2N6kBwoxTj+kL2Th14xEnTkDr1rB5\ns/NvIli6cymjPhnF0iFL3Q7FxFisx6kbk3AWL4bOnRMnoYONfjH/Ykndp6xuWXOzZ1NhrpdEaEsv\njX5JhPb0s7Dj1I1JJKpOPf2vf3U7ktiy0S+mjPXUfcrGAtfM6tXOiJdLLvnXvkRoSy8l9URoTz+z\npG6SyuzZ8IMfgAfupo+potIiGqd4I6kbd1lS9ymrW9bMrFnQr9IyL4nQllZTN2UsqZuksW0b7NsH\nV1zhdiSx56Xyi3GXJXWfsrpl9GbPhj59ICWl4v5EaEsvDWlMhPb0s7BJXUSyRGSTiGwRkZHVHPfv\nIlIqIj+MbYjGxMasWU49PRHZHaWmTLVJXURSgIlAFtAZGCginao47gXgI2wR6jphdcvoHDwIa9bA\njTee/lwitKWXyi+J0J5+Fq6nnglsVdXtqloCTAP6hThuODAd2B/j+IyJib/85V8TeCUiL03oZdwV\nLqmnA7uCtncH9pUTkXScRP96YJdN8FIHrG4ZnbKhjKEkQlt6afRLIrSnn4VL6pEk6PHAqMBsXYKV\nX4zHFBXBJ59A795uRxI/Xiq/GHeFmyYgD2gXtN0Op7cerDswLbA2YivgVhEpUdU5lY5j8ODBtG/f\nHoDU1FQyMjLKf6uX1eFsO7Lt8ePHW/tFuL1wIbRvn8P69aGfD64BeyHemmzv+r9dbD22FbriejyJ\n0J5ut192djZAeb6MRrVT74pIfWAzcAOwB/gCGKiquVUc/xYwV1VnhnjOpt6NoZycnPILwlTvZz9z\nZmV89NEp3EauAAAR00lEQVTQzydCW97yp1t47IrHuOWiW9wOJSHa00tiOvWuqpYCw4AFwEbgPVXN\nFZGhIjK0dqGa2rAfmsicPOlM4FX5LtJgidCWXiq/JEJ7+lnYWRpVdT4wv9K+SVUce1+M4jImJpYv\nh7POggsucDuS+PLSzUfGXXZHqU8F1y1N1ULN9VJZIrSll0a/JEJ7+pnNp24SlirMmOE8Ep2Xyi/G\nXdZT9ymrW4a3dq0zxe5ll1V/XCK0pZduPkqE9vQzS+omYU2fDj/+ceLNnR6K9dRNGUvqPmV1y+qp\nwp//7CT1cBKhLb00oVcitKefWVI3CWn9eiguhu7d3Y6kbtjoF1Om2puPYvpCdvORqUPPPANHj8KL\nL7odSfyVniql0fONKH2qFEmGWlOSienNR8b4VVk9PRmUlV4soRuwpO5bVres2saNcOQIZGZGdrzf\n29JrH5L6vT39zpK6STgzZsCPfgT1kuTqtnq6CZYkl33isbHAVfvzn52kHim/t6WX7iYF/7en31lS\nNwllwwY4dAiuvtrtSOqO18ovxl2W1H3K6pahvfsu3HVXdKUXv7ell+4mBf+3p9+FvfRFJEtENonI\nFhEZGeL5fiLypYisEZFVInJ9fEI1pnqqMHUq3H2325HULS/deGTcV+2EXiKSAkwEbsRZBWmFiMyp\ntEjGJ6o6O3D8d4EPgIviFK8JsLrl6ZYvh0aNICMjuu/ze1t6rfzi9/b0u3A99Uxgq6puV9USYBrO\nItPlVPVo0GYz4EBsQzQmMlOnwsCByTHXSzAb/WKChUvq6cCuoO3dgX0ViMgPRCQXZzGNX8QuPFMV\nq1tWVFoK773nJPVo+b0tvdZT93t7+l24pB7Rff2qOktVOwF9gSm1jsqYKH36KZx3Hlx8sduR1D2v\nDWk07gq3SEYe0C5oux1Obz0kVf1MROqLyL+p6sHKzw8ePLh8dezU1FQyMjI8sXq3H7fL9nklHre3\nX345J3AHafTfX7aCu5feTzTbRaVFfLvxW3JSvXE9+L093d7OyckhOzsboDxfRqPaCb1EpD6wGbgB\n2AN8AQwM/qBURC4EvlFVFZHLgT+r6oUhzmUTepm4KCqCtm2dMept27odTd0bv2w82w5tY8KtE9wO\nxcRBTCf0UtVSYBiwANgIvKequSIyVESGBg77EbBORNYAE4C7aha6iUbZb3YDc+Y4U+zWNKH7vS29\nVn7xe3v6Xdg1SlV1Ps4HoMH7JgV9PQ4YF/vQjInMW2/Bffe5HYV7bPSLCWbzqRtf27XLWYM0Lw+a\neKezWqdGfjyStCZpjLpmlNuhmDiw+dRNUvnjH2HAgORN6GB3lJqKLKn7lNUt4dQpmDwZhgyp3Xn8\n3pZFpUU0qt/I7TDK+b09/c6SuvGtzz5zeujf+57bkbirsLiQZg2buR2G8QhL6j4VPF49Wb31ltNL\nr+20AH5vy8LiQpo3bO52GOX83p5+Z0nd+FJBAcyaBT/5iduRuK+guIDmjbyT1I27LKn7VLLXLd9/\nH3r1grPOqv25/N6WBScKPNVT93t7+p0ldeM7qvD66/DAA25H4g0FxQVWUzflbJy68Z3ly52FMLZs\nSZ7FpavT9qW2fPGzLzinxTluh2LiwMapm4Q3cSL8v/9nCb2M1z4oNe6yHwufSta65b59MG9ebKcF\n8HNbqipHS456qvzi5/ZMBJbUja/84Q9wxx2QluZ2JN5wtOQojes3JqVeituhGI+wmrrxjdJSOP98\np6d+2WVuR+MNewv20m1SN/4x4h9uh2LixGrqJmHNnu0kdUvo/2IjX0xlESV1EckSkU0iskVERoZ4\n/h4R+VJE/k9E/iYil8Y+VBMsGeuWEyfCsGGxP6+f27KwuNBzNx75uT0TQdikLiIpwEQgC+gMDBSR\nTpUO+wa4VlUvBZ4D/jfWgZrktmIFfP019O/vdiTe4rUbj4z7IumpZwJbVXW7qpYA04B+wQeo6ueq\nejiwuRywAbNxlmzza4wdCyNGQIMGsT+3n9vSi1ME+Lk9E0EkST0d2BW0vTuwryr3Ax/WJihjguXm\nwt/+ZneQhmI9dVNZ2OXsgIiHrIhIL2AIcHWo5wcPHly+OnZqaioZGRmeWL3bj9vjx49PmvZ74QXo\n0yeHL76I3+rtZbzwfqPZLmjuJHWvxOP39vTCdk5ODtnZ2QDl+TIaYYc0isgVwBhVzQpsjwZOqeoL\nlY67FJgJZKnq1hDnsSGNMZSTk1N+QSSynTshI8Opp8drbLqf2/Klv7/E7iO7eSXrFbdDKefn9vSi\neAxpXAlcLCLtRaQhMACYU+lFz8VJ6D8JldBN7CXLD82LL8L998f3ZiM/t6UXR7/4uT0TQdjyi6qW\nisgwYAGQArypqrkiMjTw/CTgaSANeF2cFQtKVDUzfmGbZLB/P0yZAhs2uB2JdxUUF9D6O63dDsN4\nSETj1FV1vqp2VNWLVHVsYN+kQEJHVR9Q1X9T1W6BhyX0OAuuWyaqsWOd2Rjbto3v6/i5LQtOeG/0\ni5/bMxFE8kGpMXVuxw54+23rpYdTUGyjX0xFNk2ATyV63fLpp+Ghh6BNm/i/lp/b0sapm8qsp248\nZ906+OgjZxEMU73C4kKb+8VUYD11n0rkuuUTT8Do0dCiRd28np/b0os3H/m5PROB9dSNp3z2Gaxf\nD9Onux2JP3ix/GLcZfOpG884eRJ69IBHH4V77nE7Gn9o+1JbVvxsBektqpu5w/iZzadufOv116FZ\nM2cYo4mM9dRNZZbUfSrR6pZ798KzzzqJXSLuk8SGX9vylJ7iWMkxvtPgO26HUoFf2zNRWFI3nvDY\nY/Dzn0OnyjP1myodKzlm65Oa01hN3bju449h6FDnA9KmTd2Oxj9sfdLkYDV14ytHjzo3GU2caAk9\nWlZPN6FYUvepRKlbPvooXHUV3HabezH4tS29OEYd/NueicLGqRvXfPABLFwIa9a4HYk/eXHaXeO+\niHrqIpIlIptEZIuIjAzx/CUi8rmIHBeR/4x9mKYyv8+vkZcHDz4I77xTd3eOVsWvbflt0bekNk51\nO4zT+LU9E0XYnrqIpAATgRuBPGCFiMxR1dygww4Cw4EfxCVKk1BOnYJ774Vhw+CKK9yOxr/yCvJI\nb243HZmKIumpZwJbVXW7qpYA04B+wQeo6n5VXQmUxCFGE4Kf65a/+hWcOOHM8eIFfm3LvCPeTOp+\nbc9EEUlSTwd2BW3vDuwzJmrvvQfZ2c7cLik2vLpW9hTuoW3zOK8gYnwnkg9KYza4fPDgweWrY6em\nppKRkeGJ1bv9uF22zyvxRLK9aRM89VRPPvkEcnNzyM31RnxlK7i73T7Rbq9fvp57vnuPZ+Ip2/Zr\ne3plOycnh+zsbIDyfBmNsDcficgVwBhVzQpsjwZOqeoLIY59BihU1ZdCPGc3HyWx3bud+vnEifAD\n++QlJjq91onpd0yny1ld3A7FxFE8bj5aCVwsIu1FpCEwAJhT1etH+sKmdsp+s/vBgQPOOPThw72Z\n0P3UlsH2FHiz/OLX9kwUYcsvqloqIsOABUAK8Kaq5orI0MDzk0SkDbACaAGcEpGHgc6qWhjH2I0P\nfPst3HQT9O4Njz/udjSJo7C4kJKTJZ4c0mjcZXO/mLg5dAhuvBF69YLf/rbuZ19MZJsPbKbPu33Y\nMtzW/Et0NveL8YQDB+CWW+D737eEHg82Rt1UxZK6T3m5bvnVV3DllXDDDfDKK95P6F5uy6p4tZ4O\n/mzPRGJJ3cTUkiVO73zUKBg71vsJ3a+8euORcZ/V1E1MqMKkSfD00zB1qlNLN/Hzi/m/4IK0C3jk\nikfcDsXEWbQ1dZul0dTawYPwwAOwYwd89hl07Oh2RIkvryCPa869xu0wjAdZ+cWnvFK3XLgQMjLg\nggvg88/9mdC90pbR2FOwx7PlFz+2ZyKxnrqpkb174b/+y6mhv/EGZGW5HVFyyTuSR3oLbyZ14y7r\nqftU8Bwwdam4GCZMgEsvhXbtIDfX/wndrbasqcPHD3Po+CHPjn7xW3smGuupm4gUF8Nbb8Gvfw2d\nOzu180sucTuq5LRkxxJ6pPegYUpDt0MxHmQ9dZ+qq7plfr7TM+/QwVl+bto0mD8/sRK632rAi7Yv\nolf7Xm6HUSW/tWeisaRuTqPqrBs6dCicf77zAei778JHHzk3FRl3Ldq+iF7nezepG3fZOHVT7quv\nnOQ9bRoUFcGQIfCzn8HZZ7sdmSlz8NhBzp9wPgcfP0iDlAZuh2PqgI1TNxE7dsypjc+f7zwKCuCO\nO2DyZGfuc7sb1HsW71jMVe2usoRuqhS2/CIiWSKySUS2iMjIKo75XeD5L0WkW+zDNJVFW7dUhZ07\nneXkHn4Y/v3f4cwz4bnnnH/ffddZyGLCBKfEkkwJ3U814A+3fOjpejr4qz0TUbU9dRFJASYCNwJ5\nwAoRmaOquUHH3AZcpKoXi0gP4HXA1oiPs7Vr14YcOqYK+/bB5s2wYUPFR716cNVVzuOVV6B7d2jS\npO5j95qq2tJrlu1exryv5rHuP9a5HUq1/NKeiSpc+SUT2Kqq2wFEZBrQD8gNOuZ24G0AVV0uIqki\n0lpV98Uh3qRWVOQk7H37YPnyfN54w/k6Lw+2bYPt251b9Zs1g4sugq5doUsXZ7Whrl2hTZvk6oFH\nKj8/3+0QwioqKeK+2ffx6q2vcuZ3znQ7nGr5oT0TWbikng7sCtreDfSI4JhzgKRJ6qdOwcmTUFrq\nPE6ehJISOH7ceRQV/evrqvYVFsLhw87jyJGK/x4+7AwtLC6Gs86C1q2dfU2bOl936QJ9+kD79nDe\neU5SN4nh8PHDLPh6AU9++iTXnnstd3S5w+2QjMeFS+qRDlep3P8L+X1nPdKHygNgTtsO/latdCKt\nfHKt4hynR6CVTlb5vKHOGf6h5a8tcvqjXopT8kipB/VSlHr1nO16KYF9QY+U+tCgOdRPg5QUpX59\nqN8AmqVAan24sAHUT6G8pXf/73p2XLeUHWVhHww8VhGSRvxf6byviI+N8LzxOGeszrv5480sPG9h\nTM8Z8tgoznu89DgHjh3g26Jv+V7b7/Haba9x84U3R/z9btq+fbvbISS1aoc0isgVwBhVzQpsjwZO\nqeoLQcf8D5CjqtMC25uA6yqXX0TExjMaY0wNxHJI40rgYhFpD+wBBgADKx0zBxgGTAv8EsgPVU+P\nJihjjDE1U21SV9VSERkGLABSgDdVNVdEhgaen6SqH4rIbSKyFTgK3Bf3qI0xxoRUZ3eUGmOMib+4\nzv0iIneIyAYROSkil1d6bnTghqVNIuKPT4A8RETGiMhuEVkTePh8Alx3RHJznYmciGwXkf8LXJNf\nuB2Pn4jIZBHZJyLrgvadISIfi8hXIvJXEUkNd554T+i1DugPLAneKSKdcerznYEs4PciYpOLRUeB\nl1W1W+DxkdsB+U3QzXVZONfiQBHp5G5UvqdAz8A1mel2MD7zFs61GGwU8LGqdgAWBrarFddEqqqb\nVPWrEE/1A95V1ZLAjU1bcW50MtGxD59rp/zmOlUtAcpurjO1Y9dlDajqZ8ChSrvLb+4M/PuDcOdx\nq3fcFucmpTK7cW5iMtEZHphv581I/iwzpwl145xdh7WjwCcislJEfuZ2MAkg+O78fUDrcN9Q61ka\nReRjoE2Ip55Q1blRnMo+sa2kmrZ9EmeOnV8Ftp8DXgLur6PQEoVdc7F3taruFZEzgY9FZFOgB2pq\nSVU1kvt9ap3UVfWmGnxbHtAuaPucwD4TJNK2FZE/ANH8AjWOytdhOyr+BWmipKp7A//uF5EPcEpc\nltRrbp+ItFHVf4jI2cA/w31DXZZfgutsc4C7RKShiJwPXAzYJ+VRCPwHl+mP86G0iU75zXUi0hDn\nw/s5LsfkWyLSVESaB77+DnAzdl3W1hzg3sDX9wKzwn1DXBfJEJH+wO+AVsBfRGSNqt6qqhtF5H1g\nI1AKPGTLIkXtBRHJwCkhbAOGuhyP71R1c53LYflZa+ADcaYCrQ+8o6p/dTck/xCRd4HrgFYisgt4\nGvgN8L6I3A9sB+4Mex7LpcYYkzhsbLgxxiQQS+rGGJNALKkbY0wCsaRujDEJxJK6McYkEEvqxhiT\nQCypG2NMArGkbowxCeT/A5lx7TuEWHmSAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c32e3d0>"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}