{
 "metadata": {
  "name": "",
  "signature": "sha256:ca4f849257508b73796124cb5acc25d968b0c2b5d35f16d6741cb92dec093808"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Model validation: A simple example"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Type of model validation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "* \"R2 is not enough\"\n",
      "* Looking at the residals\n",
      "* Cross validation when you can"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "R2 and graphical methods"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy as sp\n",
      "import scipy.linalg as lin\n",
      "import matplotlib.pyplot as plt\n",
      "from __future__ import division, print_function"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# function to compute the R2 from Y and X\n",
      "# from scipy import linalg as lin\n",
      "\n",
      "def projR(to_proj, X_):\n",
      "    \"\"\" \n",
      "    projects the 2d array `to_proj` onto the residual space of `_X`\n",
      "    \n",
      "    The svd of _X is computed to have a orthonormal basis of the\n",
      "    space of the columns of _X\n",
      "    \"\"\"\n",
      "    u_, _, _ = lin.svd(X_, full_matrices=False)\n",
      "    return (to_proj - u_.dot(u_.T.dot(to_proj)))\n",
      "\n",
      "## try this function with the mean:\n",
      "# ones_n1 = np.ones((n,1))\n",
      "# y = np.random.normal(3,1,size=(n,1))\n",
      "# print(\"this should be close to 3: \", y.mean())\n",
      "# print(\"this should be close to 0: \", projR(y, ones_n1).mean())\n",
      "\n",
      "\n",
      "def R2(Y, X, rm_intercept=False):\n",
      "    \"\"\"\n",
      "    Compute the coefficient of determination R2 \n",
      "    Y: numpy array shape (n,1) : the data\n",
      "    X: numpy array shape (n,p) : the design matrix\n",
      "    \"\"\"\n",
      "    (n,p) = X.shape\n",
      "    SST = (Y**2).sum()               # total sum of square\n",
      "    C = (Y.sum())**2 / n             # sum of square of the mean\n",
      "    ones_n1 = np.ones((n,1))\n",
      "    Xdemeaned = projR(X, ones_n1)\n",
      "    # compute the svd to get an orthonormal basis \n",
      "    u_, _, _ = lin.svd(Xdemeaned, full_matrices=False)\n",
      "    \n",
      "    # compute SSReg noticing that Y^t u u^t Y = \\sum (u^t Y)^2\n",
      "    SSReg = ((u_.T.dot(Y))**2).sum()\n",
      "    if rm_intercept:\n",
      "        return SSReg / (SST - C)\n",
      "    else: \n",
      "        return SSReg / SST\n",
      "\n",
      "\n",
      "def SSR(Y, X):\n",
      "    \"compute the sum of square of the regression on X\"\n",
      "    (n_, p_) = X.shape\n",
      "    if n_ < p_:\n",
      "        raise ValueError(\"n must be >= to p\")\n",
      "    u_, _, _ = lin.svd(X, full_matrices=False)\n",
      "    return (((u_.T * Y.T).sum(axis=1))**2).sum()\n",
      "\n",
      "def Proj(X):\n",
      "    \"\"\" This function returns the projector on the \n",
      "        space of the columns of `X`\n",
      "    \"\"\"\n",
      "    (n_, p_) = X.shape\n",
      "    if n_ < p_:\n",
      "        raise ValueError(\"n must be >= to p\")\n",
      "    u_, _, _ = lin.svd(X, full_matrices=False)\n",
      "    return u_.dot(u_.T)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# define a true model\n",
      "\n",
      "def quadratic_model(x, param=[1,-8,16]):\n",
      "    return param[0]*x**2 + param[1]*x + param[2]\n",
      "\n",
      "model = quadratic_model"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 16 # number of points\n",
      "noise_level = 1./8\n",
      "x = np.arange(n)\n",
      "true_y = model(x)\n",
      "noiseSD = (true_y.max() - true_y.min())*noise_level\n",
      "y = true_y + np.random.normal(0, noiseSD, size=(n,)) \n",
      "plt.plot(x,y,x,true_y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 65,
       "text": [
        "[<matplotlib.lines.Line2D at 0x69ed8d0>,\n",
        " <matplotlib.lines.Line2D at 0x69edb10>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD9CAYAAABZVQdHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVOX3B/DPIJCaKYsyKINBCsK44I5LJqm4i4iKYiWh\n5q+w3EVMS20R0KzMotWFckVRwEISNcxERERxAQUVlHVcAJdQh4H7++P5jiINMMudDc779bov4c6d\n+xwpz1yee+55BBzHgRBCSMNkou8ACCGEaA8leUIIacAoyRNCSANGSZ4QQhowSvKEENKAUZInhJAG\nTKkkP2PGjM1CoVDStWvXCzVfW79+/SITE5OqkpISK/m+kJCQZU5OTtkuLi6XDx06NJzPgAkhhChP\nqSQfEBCwJT4+fmTN/Xl5efYJCQmeL7/88g35voyMDPHu3bunZGRkiOPj40cGBgaGV1VV0W8MhBCi\nB0ol30GDBh23tLQsrbl/4cKFX65duzao+r6YmJjxfn5+O83MzCocHBxyO3bseDUlJaUvXwETQghR\nntpX2DExMeNFIlF+t27dzlffX1hY2E4kEuXLvxeJRPkFBQV2mgRJCCFEPabqvKm8vLz5mjVrPkxI\nSPCU7+M4TlDb8QKB4D+9ExTtI4QQUr+68m1Nal3JX7t2rUNubq6Dm5tbuqOjY05+fr6oV69eZyQS\nidDOzq4gLy/PXn5sfn6+yM7OrqCWQA1+W7lypd5joDgpToqTYpRvqlIryXft2vWCRCIR5uTkOObk\n5DiKRKL8tLS0nkKhUOLl5RW7a9euqVKp1DwnJ8cxOzvbqW/fvinqjEMIIUQzSiV5Pz+/nQMGDEjK\nyspytre3z9uyZUtA9derT72IxeIMX1/fSLFYnDFq1KiD4eHhgTQ1Qwgh+iFQ5/Kfl4EFAk5fY6si\nMTERHh4e+g6jXhQnvyhOfhlDnMYQIwAIBAJwKszJU5InhBAjomqSp4eUCCGNXmUlsG8f+7OhoSRP\nCGn0zp0DJk4ERo0C7tzRdzT8oiRPCGn0kpMBf3+gRw+gd2/gzBl9R8QfSvKEkEYvORl49VUgLAxY\nvx4YORLYtEnfUfGDbrwSQho9JycgOhro3Jl9f/kyMGECMGgQsHEj8MIL+o2vOrrxSgghKrh9G7h1\nC3B1fbbPxQVISQFKSliiv3lTf/FpipI8IaRRO3UK6NsXMKmRDV96CdizB/D1BdzdgSNH9BOfpijJ\nE0IateRkoF8/xa8JBMDixcCOHcCbb7I5e2ObZaYkTwhp1OpK8nKvv86mb/btY6WW9+/rJjY+UJIn\nhDRalZUsebu713+svT3w99+AjQ2b3snI0H58fKAkTwhptDIzAaEQaN1aueNfeAH44Qdg6VJg8GA2\nZ2/oKMkTQhqtkyeB/v1Vf19AAPDnn0BQEJuzl8n4j40vlOQJIY2WMvPxtenZE0hNBS5eBIYNAyQS\nfmPjCyV5QkijpUmSBwBra+CPP1gtfe/e7HyGhp54JYQ0SmVlgEgElJYCZmaany82Fpg1C1i9Gnj3\nXVZ+qQ30xCshhCjh9Gk25cJHggcALy/gxAkgPJzN2T96xM95NUVJnhDSKGk6VaOIkxM7r1QKDBgA\nXL/O7/nVQUmeENIoJSerV1lTnxdfBLZvB95+G5g2Tf9PyCqV5GfMmLFZKBRKunbtekG+b8mSJetc\nXV0z3dzc0n18fPbdu3evlfy1kJCQZU5OTtkuLi6XDx06NFwbgRNCiLo4jiV5ZR6CUodAAMybBxw7\npr25eWUpleQDAgK2xMfHj6y+b/jw4YcuXbrUOT093c3Z2TkrJCRkGQBkZGSId+/ePSUjI0McHx8/\nMjAwMLyqqop+YyCEGIzsbKBFC6BdO+2N8Vj2GPuydkLfBSZKJd9BgwYdt7S0LK2+z9PTM8HExKQK\nANzd3U/l5+eLACAmJma8n5/fTjMzswoHB4fcjh07Xk1JSenLf+iEEKIebczH1xR8OBhRmVHaHUQJ\nvFxhb968ecbo0aPjAKCwsLCdSCTKl78mEonyCwoK7PgYhxBC+KDtJB9/NR5RmVH4adxPEOh5vsZU\n0xN8/vnny83NzaXTpk3bUdsxAoFA4e8rq1atevq1h4cHPDw8NA2HEELqlZwMvPWWds59699bmBEz\nA9t9tsOqmZXG50tMTERiYqLa79coyW/duvXtuLi40UeOHBkq32dnZ1eQl5dnL/8+Pz9fZGdnV6Do\n/dWTPCGE6MK//wJXrrAaeb5xHIcZMTMw3W06Xnd8nZdz1rwAXr16tUrvV3u6Jj4+fuS6deuWxMTE\njG/atOlj+X4vL6/YXbt2TZVKpeY5OTmO2dnZTn379k1RdxxCCOFTairQtat21m0NPx2O4ofF+OT1\nT/g/uZqUupL38/PbeezYscF37txpbW9vn7d69eqVISEhy6RSqbmnp2cCAPTv3/9keHh4oFgszvD1\n9Y0Ui8UZpqamsvDw8MDapmsIIUTXtDUff+nWJaw6tgonZpyAeRNz/gdQE/WuIYQ0KhMmAFOmAFOn\n8nfOx7LH6PtzX8xzn4eZPWfyd2IFVO1dQ0meENJocByrjT95EnBw4O+88+PnI/9+PvZM3qP1ahpV\nk7zG1TWEEGIsbt5kif7ll/k7p7xcMv3ddL2XSypCSZ4Q0mjI5+P5ysV8l0tqA7UbIIQ0Guou96eI\nNsoltYGSPCGk0eCzssYQyyUVoRuvhJBG4ckTwNISuH2btQPWxKVbl+AR4YETM07A2dqZnwCVRCtD\nEUKIAmfPAp06aZ7gH8sewy/KD6FDQ3We4NVBSZ4Q0ijwNVUTfDgYztbOmNFjhuYn0wGqriGENArJ\nycCoUZqdw9DLJRWhK3lCSKOg6XJ/8nLJX71/NdhySUUoyRNCGryiIuDBA7bQtjo4jsPM2JkGXy6p\nCE3XEEIaPPl6rurOsISfDkfRgyJE+ep/pSdVUZInhDR4mtx0NdTuksqi6RpCSIOnbpI3tnJJRehh\nKEJIgyaTARYWQH4++1MVuuwuqSzqQkkIIdVcuMC6Tqqa4I2xXFIRSvKEkAbt5EnVp2qMobuksmhO\nnhDSoKk6H2/M5ZKKUJInhDRoqiZ5ebmkoXeXVJZSSX7GjBmbhUKhpGvXrhfk+0pKSqw8PT0TnJ2d\ns4YPH36orKzs6YxXSEjIMicnp2wXF5fLhw4dGq6NwAkhpD537wLFxYBYrNzx8nLJHRN3GGW5pCJK\nJfmAgIAt8fHxI6vvCw0NDfb09EzIyspyHjp06JHQ0NBgAMjIyBDv3r17SkZGhjg+Pn5kYGBgeFVV\nFf3GQEgjVVGhv7FPnQL69AGaNKn/2IZQLqmIUsl30KBBxy0tLUur74uNjfXy9/ePAAB/f/+I6Oho\nbwCIiYkZ7+fnt9PMzKzCwcEht2PHjldTUlL68h86IcTQVVayypZLl/QzvipTNcbWXVJZalfXSCQS\noVAolACAUCiUSCQSIQAUFha269evX7L8OJFIlF9QUGCn6ByrVq16+rWHhwc8PDzUDYcQYoCys1nf\nmG++AX78UffjnzwJzJtX/3HRl6OxL3Mfzr17zuDKJRMTE5GYmKj2+3kpoRQIBJxAIKj1yabaXque\n5AkhDU96OjBwIBAZCYSEAFY6rEasrARSUljPmrpcuXMF7xx4B7/7/W6Q5ZI1L4BXr16t0vvVnisX\nCoWS4uJiWwAoKipqa2NjcwsA7OzsCvLy8uzlx+Xn54vs7OwK1B2HEGK8zp0Dhg8Hxo8Hfv5Zt2Nf\nvgy0acO22jx48gATdk/AmiFr4C6q59PASKmd5L28vGIjIiL8ASAiIsLf29s7Wr5/165dU6VSqXlO\nTo5jdna2U9++fVP4CpgQYjzOnQO6d2dTJt9+q9ubsPXNx3Mch4CYAAxsPxDv9HpHd4HpmFLTNX5+\nfjuPHTs2+M6dO63t7e3zPvnkk4+Dg4NDfX19Izdt2jTTwcEhNzIy0hcAxGJxhq+vb6RYLM4wNTWV\nhYeHB9Y1lUMIabjkSb59e+CVV4D9+wFfX92MXV+SX5e0Dnn387DdZ7tuAtITalBGCNEKiQRwdWW1\n6gIBsG8f8MUXQFKSbsbv2hXYsgXo3fu/rx2+fhhv7X8LKbNSYN/K/r8HGDBVG5RR/TohRCvS09lV\nvLxYZfx4VmmTooPJ2/v3gZwcwM3tv6/dKLuBN/e9iR0+O4wuwauDkjwhRCvkUzVyTZoA778PbNig\n/bFTUoAePQAzs+f3P6p4BJ9IHwQNDGoQfWmUQUmeEKIV587990p65kzg4EGgsFC7Yyuaj+c4DoFx\ngXCycsKCfgu0G4ABoSRPCNGKmlfyAOvpPm0aEB6u3bEVJfkfUn9AamEqNnltMrgHnrSJbrwSQnhX\nXg5YWwP37gHmNfp8ZWUBr74K3LgBNGvG/9gcx2rj09MBu/89a5+UlwTvXd44MeMEnKyd+B9Uh+jG\nKyFE7y5eBFxc/pvgAcDZmTUN27FDO2Nfu8Y+POQJvvhhMXz3+GLz+M1Gn+DVQUmeEMI7eWVNbebN\nYzdgtfHLfHIy0L8/+7qisgKT90zGrJ6zMNZ5LP+DGQFK8oQQ3im66VqdpydbYPuvv/gfu/pyf4sT\nFqPlCy3x8eCP+R/ISFCSJ4TwTtFN1+oEgmdX83yT33Tddn4b/sj6A9smbIOJoPGmOrrxSgjhVVUV\n0KoVcPMmYGlZ+3Hl5azXfHIy0KEDP2OXlwOtWwNHM85h3B5PHJ1+FF2FXfk5uYGgG6+EEL26do1V\n1tSV4AGgeXNWN79xI39jnzkDuPQowbQYH2wctbHBJXh1UJInhPCqvqma6ubMAX79lbUh4EPSyUrc\n9ZiGCa4TMLXLVH5OauQoyRNCeFVfZU119vbsJuyWLfyMvTV3FV5s9QRhw8L4OWEDQEmeEMKr+ipr\napo/ny0PWFmp2bj7M6OR1TwCEWN3w9SEl0XvGgRK8oQQXqkyXQOwShhra+CPP9Qf88qdK5gV8w5a\n/bkHvV1t1D9RA0RJnhDCm9u3gYcPAQcH5d8jELCr+a+/Vm9M+RJ+kyzXYNAr7mhEbWmUQkmeEMKb\n9HQ2VaNqop00CbhyBTh/XrX3VV/C78Ur79S5ElRjRUmeEMIbVadq5MzNgcBA1R+OWpe0Djfv3cTG\nURvrXe6vsaIkTwjhjSqVNTXNns2WCLx9W7njD18/jK+Sv0KUbxQElU2Rns4an5HnaZzkQ0JClnXu\n3PlS165dL0ybNm3HkydPXigpKbHy9PRMcHZ2zho+fPihsrIyCz6CJYQYNlUra6pr0waYOBH46af6\nj625hF96OuDkBLRood7YDZlGST43N9fh559/fictLa3nhQsXulZWVjbZtWvX1NDQ0GBPT8+ErKws\n56FDhx4JDQ0N5itgQohhevwYuHoVEIvVP8e8eWxBEam09mPKK8r/s4Rf9aZk5HkaJfmWLVveNzMz\nqygvL28uk8lMy8vLm7dr164wNjbWy9/fPwIA/P39I6Kjo735CZcQYqguXWJX002bqn+Orl1ZH/q9\nexW/LquSYereqehi0+W5JfxoPr52Gj0xYGVlVbJo0aL17du3v9msWbNHI0aM+NPT0zNBIpEIhUKh\nBACEQqFEIpEIFb1/1apVT7/28PCAh4eHJuEQQvRI3ZuuNc2bB3z2GeDn93yVDsdxmHtwLh7JHmHv\nuL3PLeGXnAysXKn52IYoMTERiYmJ6p+A4zi1t6tXr3ZwdXXNuHPnjnVFRYWpt7f3/t9+++1NCwuL\n0urHWVpaltR8LxuaENJQvP8+x61fr/l5ZDKOe+UVjktKen5/yPEQrtv33bh7j+89t7+oiOMsLDiu\nslLzsY3B/3Kn0nlao+ma1NTU3gMGDEiytra+a2pqKvPx8dl38uTJ/ra2tsXFxcW2AFBUVNTWxsbm\nlibjEEIMnyaVNdU1aQLMnfv8w1Hbzm/D96nfI25aHFq+0PK540+dAtzdAROqFVRIox+Li4vL5eTk\n5H6PHj1qxnGc4PDhw8PEYnHGuHHjDkRERPgDQEREhL+3t3c0P+ESQgxRVdWzB6H4EBAAJCQAeXnA\nketHsOjQIsRNi4NdS7v/HEvz8XXTeNGQtWvXBkVERPibmJhU9ezZM+2XX36Z9eDBg5d8fX0jb968\n2d7BwSE3MjLS18LCouy5gWnREEIajOvXgcGDWVLmy/z5wINmF3DAaigiJ0fCw8FD4XEeHsCyZcCI\nEfyNbchUXTSEVoYihGhs3z7WLvjAAf7O+c/5fAz+dQA2+62Ffy/FveFlMrY4SX2rUDUktDIUIUTn\n+Kqskbv3+B7eOzEKLvfm4nFq7Yt/XLwIiESNJ8Grg5I8IURjfN10BQBppRQTdk/A6w6v45upi7Bh\nA1DbL/00H18/SvKEEI1p0s6guiquCjNiZsCiqQW+GvEVhgwRwNSU3YRVhJJ8/SjJE0I0UlIClJYC\nr7yi+bmWH12O66XXsd1nO5qYNIFAwB6Oqq07JSX5+lGSJ4RoJD0d6NZN8zr18NPhiMqIQqxfLJqZ\nNXu6f9o04PRp1m++upISoLAQ6NJFs3EbOkryhBCN8HHTNeZyDD77+zPEvxmP1s1bP/das2asDfHG\njc+/59QpoHdv9vAUqR0leUKIRjRN8sn5yZh1YBZipsbgFUvFcz6BgcD27UBZtadtaKpGOZTkCSEa\n0aSyJvtuNibsnoCt47eij13tK360aweMHg1s2vRsHyV55dDDUIQQtUmlQKtWbH68WbP6j6/u1r+3\nMGDTAAQNDMLsXrPrPT4lBfD1ZT3rTUwAKys2Ty9U2OO24aKHoQghOpORwapqVE3w/0r/xbid4+DX\n1U+pBA8AffuyK/rYWJbcrawaX4JXByV5Qoja1JmPl1XJ4BflB5fWLvjE4xOV3jt/PutOmZwM9O+v\n2riNFSV5QojaVE3y3P8W/ngse4yfx/383MIfypgwAcjJAX78kebjlUVJnhCiNlWTfNiJMCTlJWGv\n716YNzFXeTwzM+D991n5JCV55dCNV0KIWjju2c1PG5v6j992fhtWHF2BpJlJaPdSO7XHLSlhV/QJ\nCYC56p8TRo9aDRNCdOLGDTYvXlhY/7FHrh/BtH3TcHT6UXS26az94BowVZO8Rgt5E0IaL2Wnas5L\nzsMvyg+RkyMpwesBzckTQtSiTJK/WnIVY3aMwTejvql1ZSeiXZTkCSFqqS/JXy25iiERQ/DRax9h\napfaF/4g2kVJnhCilrqSvDzBr3hthdIPOxHt0DjJl5WVWUyaNGmvq6trplgszjh16pR7SUmJlaen\nZ4Kzs3PW8OHDD5WVlVnwESwhxDCUlQG3bwMdOvz3NUrwhkXjJD9v3rwNo0ePjsvMzHQ9f/58NxcX\nl8uhoaHBnp6eCVlZWc5Dhw49EhoaGsxHsIQQw3D+PNC163/b/FKCNzwalVDeu3evVY8ePc5ev379\nuf6gLi4ul48dOzZYKBRKiouLbT08PBIvX77s8tzAVEJJiNH65hsgMxP4/vtn+yjB64ZOSyhzcnIc\n27RpczsgIGBLenq6W69evc58/fXX8yUSiVAoFEoAQCgUSiQSicI2QqtWrXr6tYeHBzw8PDQJhxCi\nI+fOAe7uz76nBK89iYmJSExMVPv9Gl3Jp6am9u7fv//JpKSkAX369Dk9f/78r1966aUH33777ful\npaWW8uOsrKxKSkpKrJ4bmK7kCTFaPXuyq3h3d0rwuqbTVsMikShfJBLl9+nT5zQATJo0aW9aWlpP\nW1vb4uLiYlsAKCoqamtjY3NLk3EIIYajogK4fJnNyVOCN3waJXlbW9tie3v7vKysLGcAOHz48LDO\nnTtfGjdu3IGIiAh/AIiIiPD39vaO5iNYQoj+Xb4MtG8PFD6mBG8MNO5dk56e7jZr1qxfpFKpeYcO\nHa5t2bIloLKysomvr2/kzZs32zs4OORGRkb6WlhYlFV/H03XEGKcfvsN2HXoKi70pASvD9SgjBCi\nVTMWX8W+l4Zg7VhK8PpAy/8R0kh88gnrBKlLV0uuYofZEPg7UoI3FpTkCTFCx48DK1cCX36puzHl\nN1mbnFiBZcMpwRsLSvKEGBmOA4KDgc8+Y/Pj9+9rf0x5gg/ssgIvZc2Gra32xyT8oCRPiJH5/XeW\n2IODAU9PYOtW7Y5XvUyyy5PZcHPT7niEX5TkCTEilZXAsmXAmjWsb8zcucDGjUBVlXbGq1kHr+qa\nrkT/KMkTYkS2bwdatQLGjmXfDxgAtGwJxMfzP5aiB50oyRsfSvKEGIknT4CPPwZCQwHB/wroBAJ2\nNf/NN/yOVduTrJTkjQ8leUKMxI8/Al26AIMGPb9/yhSWfC9f5mec2hL8/ftAURHg7MzPOEQ3KMkT\nYgQePGDz8GvW/Pe1pk2B2bPZ3Lym6upFc+EC+5Cp2UOeGDZK8oQYgS+/BIYNA7p1U/z6u+8CO3aw\nFZvUlV6cjtcjXq+1VcG5c6DKGiNESZ4QA3f7Nptz/+ST2o9p1w4YNQrYskW9MeKvxsPzN098OfzL\nWp9kpfl440RJnhAD9/nnwLRpwCuv1H3c3LnAt9+yMktV/HTmJwTEBCB6ajQmd55c63GU5I0TNSgj\nxIDl5gK9egEZGYBQ4fpqz3AcW8Tjo4+AcePqP3cVV4UPj3yIfZn7EPdGHDpadaz1WJmMlWpKJMBL\nL6n2dyD8ogZlhDQgq1YBc+bUn+ABVk45bx6wYUP9xz6WPYZflB9O5J3AyZkn60zwAHDlCiASUYI3\nRpTkCTFQFy8CBw8Cixcr/57Jk4FLl9hWmzvldzD016EwEZgg4a0EWDe3rve86el009VYUZInxEAt\nX87607Rsqfx7zM1ZpU1t5ZTZd7PRf1N/DH55MLb7bEdT06ZKnZfm440XJXlCDFBSEkus772n+nv/\n7/+A3buB0tLn95+4eQKDtgxC0IAgrBm6BiYC5f/5U5I3XpTkCTEw8lbCq1axB51UZWvLetts2vRs\n3+6LuzFh9wREeEfgnV7vqBwPJXnjRUmeEANz8CBw9y4wfbr655CXU1ZUcAj7JwxLEpYg4a0EjOg4\nQuVzFRWxRN+unfrxEP3hJclXVlY26dGjx9lx48YdAICSkhIrT0/PBGdn56zhw4cfKisrs+BjHEIa\nuqoq1kr48881ax/Qpw/Q1k6GsT+8i12XduHkzJNws1Xvzqn8Kl6gdNEeMSS8JPkNGzbME4vFGQKB\ngAOA0NDQYE9Pz4SsrCznoUOHHgkNDQ3mYxxCGrqdO4HmzYHx4zU7z/0n91HuPRap2Xn4++2/YdfS\nTu1zUWWNcdM4yefn54vi4uJGz5o16xd5gX5sbKyXv79/BAD4+/tHREdHe2s6DiENnVTKHmSq3kpY\nHfn38zFoyyD0c3HEC/tikXNFs+J2mo83bqaanmDBggVfrVu3bsn9+/efFnpJJBKhUCiUAIBQKJRI\nJBKFj3KsWrXq6dceHh7w8PDQNBxCjNbPPwOdOgGDB6t/jnPF5zBu5zjM7TsXiwcshn26ABs3snOr\nfc5z7MOH6EdiYiISExPVPwHHcWpvBw4cGBsYGPgdx3H466+/PMaOHXuA4zhYWFiUVj/O0tKypOZ7\n2dCEEI7juAcPOM7WluPS0tQ/R1xWHNdmbRsu8mLk0323bnGchQXH3b6tflzNmnGcVKp+XIRf/8ud\nSudpja7kk5KSBsTGxnrFxcWNfvz4cdP79++3fOutt34TCoWS4uJiW1tb2+KioqK2NjY2tzQZh5CG\n7uuvgddfB3r0UO/9P535CSsTVyJ6ajQG2A94ur9NG8DbG/jlF1aWqaoLFwCxGDAzUy8uon+8NSg7\nduzY4C+++GLxgQMHxgUFBa21tra+u3Tp0rDQ0NDgsrIyi5o3X6lBGSHMnTuAiwuQnAx0rLuFzH9U\ncVVYdmQZ9mfur7XJWFoau5GbkwOYqnhZ9/33wJkz7EOCGAa9NiiTV9cEBweHJiQkeDo7O2cdPXp0\nSHBwcCif4xDSkISEAL6+qif4x7LHmLp3Kk7crLvJWM+egIMDEB2temxUWWP8qNUwIXp08yaborl4\nEWjbVvn33Si7galRU+Fg4YAt47fU24Nmzx7Wz+bvv1WLr18/YN26/64rS/SHWg0TYkRWr2YNxVRJ\n8FEZUejzcx/4uPgo3WRswgTWm/7sWeXHqaxkc/K1LTlIjIPGJZSEEPVkZgIHDgDZ2cod/6jiERYd\nWoT4q/H4fdrv6GvXV+mxTE2BwEC2jKCySwRmZ7M+OK1aKT0MMUB0JU+InixfDgQFKZdEM29nwv0X\nd5Q8KsHZ/zurUoKXmzWLzcvfUrLWjR6CahgoyROiB6dOAadPs1Wf6sJxHDalbcJrW1/DPPd52Dlx\nJ1o1Ve/SunVrYOJE5R+MOneObro2BHTjlRAd4zhgyBDgzTeBmTNrP+7e43t49493cfHWReyetBvi\nNmKNxz5/Hhg9mpVT1lf7PmoU62fv5aXxsIRHdOOVEAN36BBr3+vvX/sxpwtOo+dPPWHR1AIps1J4\nSfAAu4nq5ARERdV/LE3XNAyU5AnRoeqthBU9mFTFVeGLpC8wZscYrB22Ft+P+R7NzJrxGsPcuewG\nbF2Ki4EnTwB7e16HJnpASZ4QHdqzhyV3H5//vnbr31sYs2MM9mXuQ8o7KZgonqiVGMaNAwoL2T2B\n2qSnUw/5hoKSPCE6UlEBrFihuJXwketH0OPHHuhh2wPH3j4GBwsHrcVhaspu+Na22DdAUzUNCdXJ\nE6IjmzYBjo7spqucrEqGlYkrsfXcVkR4R2DYK8N0EsvMmUCHDmxaxtb2v6+fOweMHKmTUIiW0ZU8\nITpQXg58+inrUyN3o+wGBm8djDOFZ5A2O01nCR4ArKyAKVOAH39U/Lp8uoYYP0ryhOjAunXAq68C\nvXqx7+WtCbw7eSPujTgIWyhcV0erPvgA+OEHtiJVdeXlrMTS1VXnIREtoOkaQrToyRNg8WIgLg44\nfPj51gQH/A7AXeSut9g6d2bbnj3AG28823/xImt9bG6ut9AIj+hKnhAtuXEDeO01ID+f9WR/3IK1\nJrj76C7O/t9ZvSZ4OUXllHTTtWGhJE+IFsTFAX37sj7xe6OqsPvqj3ht62uY6z4XuybuUrs1Ad/G\njGGLlpye/Ps6AAAZyElEQVQ69WwfJfmGhaZrCOGRTAasXAn8+it7qrTZK2cwcPMcCAQCJPonorNN\nZ32H+JwmTYD33wc2bAB27GD70tPZhxNpGKh3DSE8KS4G/PxY4vxucwm+vrAc+zP3I2RoCPy7+8NE\nYJi/OJeVsdLOS5eetRa+eROwtNR3ZEQR6l1DiB4cO8YqZwa9VoUpYZswOFKMJoImyJyTiYAeAQab\n4AHAwgKYNo1V2ly7BlhbU4JvSAx2ukYqBR48AB4+ZFv1r2t+L//aygr47DPVFysmRF1VVaw88quv\ngBXfncG20jkwOWeCg28cRI+2PfQdntI++ADw8ACcnWk+vqHRaLomLy/Pfvr06b/eunXLRiAQcLNn\nz/5p7ty535SUlFhNmTJl940bN152cHDIjYyM9LWwsCh7bmCBgBs/nqs1YXMc8NJLQIsWz/6Ub9W/\nr/51ZCTrsvfVVxr/XAipV0kJ6yRZfK8ELnOWIyHP8Kdm6jJyJHD1KmuBvGqVvqMhtVF1ukajJF9c\nXGxbXFxs271793MPHz5s0atXrzPR0dHeW7ZsCWjduvWdoKCgtWFhYUtLS0stQ0NDg2sEyu3fzylM\n1i1aAC+8oHo8paWAuzvr8hcQoPZfi5B6paYCkyZXwcl3My60WYFJ4kn49PVPYdnMeOc54uJYtc3+\n/YC3t76jIbVRNcmD4zjetvHjx0cnJCQM69Sp0+Xi4mIhx3EoKiqy7dSp0+Wax7Kh+ZeZyXFt2nDc\niRNaOT1p5KqqOC48nOMsxKmcU5g71/+X/lxaYZq+w+JFZSXH+fhwXFGRviMhdflf7lQ6L/M2e52b\nm+tw9uzZHu7u7qckEolQKBRKAEAoFEokEonCZ7ZXVfud0MPDAx4eHhrH4eICbN0KTJrEan+pHzbh\ny8OHgP+7JUhsshym0/djmafxTs0oYmKi3GIiRLcSExORmJio9vt5KaF8+PBhi8GDBx/76KOPPvX2\n9o62tLQsLS0tffp7q5WVVUlJSYnVcwNruYRy7Vo2R3/8ONCM3zUXSCN08VIVPIM2416vFZjeexJC\nPI17aoYYL52XUFZUVJhNnDgx6q233vrN29s7GmBX78XFxbYAUFRU1NbGxkbJ9eH5s2QJu6qfOZPd\nxCVEXZ9uOoMe3w5Ai0GbceK9g/jB61tK8MRoaJTkOY4TzJw5c5NYLM6YP3/+1/L9Xl5esREREf4A\nEBER4S9P/rokELBV6bOzgbAwXY9OGoLC0hJ0DnoPq7PH4OPR/4crS/8xqrJIQgANp2v++eefV197\n7bW/u3Xrdl4gEHAAEBISsqxv374pvr6+kTdv3mxfVwmlNqdr5AoKWA+RH38Exo7V+nCkAajiqhB2\naDM+TlyBl/+dhCMrPsXLNnTlTgyDTksoNSEQCLjViasxp88cWDe31upYycmAlxeQmAiI+Vn0njRA\nsioZojKi8GFcKG5eb4YFzt8hbEEPWueUGBSjamtw494NOG10woI/FyDvXp7WxunXjz2VOH48e4CF\nkOoeVTzCD6k/oNO3nbAybiNKoj5B4vR/sHYhJXhi/PSa5Dd5bcL5986jiaAJ3H5ww9vRbyPjdoZW\nxvL3Z6vUT5nCOgUSUva4DCHHQ+C4wRF/ZP+Bt1pE4N5X/yDxp3EYOKBhlEUSovf/k0UtRfhi+Be4\nNvcanKycMCRiCMbvGo+kvCTex1q7ltUCL1nC+6mJESl8UIighCB0+KYDMu9k4vD0w5hadQA/rXgV\nhw4Bbm76jpAQ/hhcq+FHFY+w5dwWfJH0BUQtRVg6cClGO42GgKffm6n1QeOVdTcL65LWISojCm+5\nvYWF/RbiZYuXEREBfPghkJBA92yI4TOqG691jS2rkmHPpT0IOxGGKq4KQQODMKXzFJg1MdN47MxM\nYPBgICYG6N9f49MRA3e64DTCToTh7xt/I7BPIN7v+z5aN28NANi8Gfj4Y7b+qouLngMlRAkNJsnL\ncRyHP6/9ibATYcgpzcGi/osws+dMNDdrrtH4f/wBzJ7NWh+IRBqdihggjuOQcD0BYSfCkH03G4v6\nL8KsnrPwovmLT4/56Sfg00+BI0dYi11CjEGDS/LVnco/hbATYTiRdwJz+szRuPwyLIytVE+tDxoO\neRlk2IkwSCulCBoYBL8ufv/5DfD774HQUJbgO3bUU7CEqKFBJ3m5y3cuY13SOuzP3A//7v5Y2G8h\n7Fup3omM41jvbI4Dtm8HlcsZsceyx9h6bivWJa1D2xZtsXTgUoxxHqOwedjGjcD69cDRo8Arr+gh\nWEI00CiSvFzB/QJ8lfwVtpzbgnHO47B4wGJ0semi0jkePQJeew2YOBEIDq7/eGJYSh+V4ofUH/BN\nyjfo3a43lg5cilfbv1rr8V99xZL80aOAg4Pu4iSEL40qycuVPipF+OlwfHf6O7R5sQ18xb6Y3Hky\nnK2Vm2jNz2cVN9T6wDiUPS5DzOUYRGZE4viN4/B28UbQwKB6P+C/+IKtY3r0KNC+vY6CJYRnjTLJ\ny1VWVeJE3gnsydiDvRl7YfOijdIJ/+RJ9kQstT4wTDUT+9BXhmKyeDLGOo9Fyxda1vv+0FBg0ybg\nr7/oRjsxbo06yVenTsLfuhX4/HMgJYVWq+c4tkh1kyb6i0HTxC732WfAtm3sCr5dOy0GTIgOUJJX\nQJWEv3AhcPEiW+/SlLd1s4yLVAq88Qbw5AkQG6vbsflK7HKrVwO7d7MqmrZttRAwITpGSb4e9SV8\nmQwYPRro3JndpGtsnjwBfH3Z15cvs5/B6NHaHbN6Yv/n5j8Y4jhEo8QOsN9EVq4E9u1jCV6ocAFK\nQowPJXkV1Jbwh4smY9oIZ3z4YeNqffD4MasyatoU2LmTPea/aBFw/jxgbs7vWLUl9nHO4/DSCy9p\ndG6OA5YvB37/nSX4Nm14CpoQA0BJXk01E36rJjbI+9MXPy30ht8w1wazWHNtHj0CJkwAWrVi89dm\nZixZjh4NDB8OLFig+RiShxLEX43XSmKX4zhg6VL2AZWQALRuzctpCTEYlOR5IE/46+L24ODVODRv\nfRvdbbujd7ve6NW2F3q16wVna+cGk/jLy1llkY0NEBHx/L2IzEz2HEFGhmpXxJKHEpwpOoPUwlSc\nKTqDM4VnUF5RjsEOg3lP7HIcx37zOHaMJXgrq/rfQ4ixoSTPsw8+ACT3SzB7ZRrOFJ55mrjulN9p\nEIn/4UPWZ799e9asS1E1zYIF7Er/hx8Un6O2hN6rXS/2s/nfz8fRwpG3bqI1cRwwbx5bBezPP6k6\nijRcBpPk4+PjR86fP//rysrKJrNmzfpl6dKlzy2nbSxJ/t9/ge7d2cpS3t7P9pc8KkFaEUv8qUWp\nOFN4xugS/4MHwJgxgJMTa9ZVW7lkaSng6grExwNtO+o/oddUVQW8/z5w9iyLsVUrnQxLiF4YRJKv\nrKxs0qlTpyuHDx8eZmdnV9CnT5/TO3fu9HN1dc2sFqhRJHkAOHECmDyZ3YCsa45XnvirJ0BDTfz3\n7wOjRgFdurBmXSYKwql+hR518gwu3zuDFy30m9BrqqoC3nuPlb0ePAi0VK8YhxCjYRBJ/uTJk/1X\nr169Mj4+fiQAhIaGBgNAcHBwaLVAjSbJA8DixUBeHqu5VsXd8rvsir/ozNPEn3c/DxZNLWDVzAqW\nTS1h2czy+T8V7fvfny3MW6iUUJ/InqD0cSlKH5U+/bOgpBRrviyFzcsl6D+kFGVPnn9d/ucLpi+g\nV9te6N2uN3oIe2Hl7F74dKEjJk82nE5uCxYAqansuYaX+J3iJ8QgqZrktfK4T0FBgZ29vf3TlblF\nIlH+qVOn3LUxlq58+inQsycQGfmsjlwZ1s2t4dnBE54dPJ/uk1ZKFSbVp0n4QQEu3Lqg8HVppVTh\nB4RZEzOF55JVyZ77oGjRxBJpJywh6mCJEZ6WsGzWHt2buin8UHnR7MXnPlBs1rCS0rFjDaM18969\nbOGXtDRK8ITURitJXiAQKHWJvmrVqqdfe3h4wMPDQxvh8KJZM1Z54uXFVpXS5OEa8ybmELYQQthC\n9ZPU9gEhrZTWm6jv3gWGDQMChgLr1qjeWvn114HevVmb3hUrVA6dV9evA4GBbPEXCwv9xkKINiUm\nJiIxMVHt92tluiY5ObnfqlWrVsmna0JCQpaZmJhUVb/5amzTNXLLlwOXLgH79xtX//nbt4GhQ9mN\n1jVqJHi5nByW6NPT9dfo68kT4NVXWeuF+fP1EwMh+qLqdI1W7v717t07NTs72yk3N9dBKpWa7969\ne4qXl5eOu6Box8cfs6vIbdv0HYnyJBJ2Fe7trVmCBwBHR3ajU5+995cuBezsWMkkIaRuWpmuMTU1\nlX377bfvjxgx4s/KysomM2fO3FS9ssaYvfAC8Ouv7CnQIUNYsjFkRUUsTj8/9gHFh+Bgtuh1UhIw\nYAA/51TW/v1AdDQrlzSm36QI0Rd6GEpNn3zCetDHxRlusikoYAne3x/48EN+z71tG7BhA1sIXVH5\npTbk5gJ9+wIHDrBFXghpjAxiuqYxWLYMuHWLLURhiG7eZDeIZ83iP8EDwLRprP3Br7/yf25FpFJg\nyhQ2VUMJnhDl0ZW8Bi5eZHPdqanAyy/rO5pncnPZFfwHH/DTWKw2KSmsqdnly9ovYVy4ELh6lZVM\nGupvToToAl3J61CXLqwh1owZ7MlLQ3D9OuDhwZK7NhM8wKZOhg1jN3O1KTYWiIpiK3dRgidENXQl\nryGZjJXzvfUWMGeOfmPJzmZlkh9+CLz7rm7GLCwEunVjc/MdOvB//hs32IfJ/v26v8lLiCEyiLYG\nSg3cQJI8AFy5AgwcyDogduyovxiGDgVWrWLz8LoUEsKmbvbv5/e8FRWszbGPD7BkCb/nJsRYUZLX\nk6++YkvNJSbqfvHrM2fYk7iffw68/bZuxwbYilJiMetkOWwYf+ddsoT1sT9wQHcVPIQYOpqT15N5\n89h88YYNuhtTJgM++4x1k9ywQT8JHmDLBa5fz54+lcn4Oefvv7NmcBERlOAJ0QRdyfPo2jVW3vfP\nP+xhIW3KygKmT2dVLVu26K/FgBzHsat4Hx/N703k5bHWCVFR7H4HIeQZupLXow4d2ENS/v78XdHW\nxHFAeDi7Cfnmm2wVJH0neID9FvP118Dq1awRmroqKoCpU1llECV4QjRHV/I8q6piLQ+GDeO/v0tB\nASvXLC0FfvsN6NSJ3/PzYc4cNr2ycaN67w8OZs3P/viDpmkIUYRuvBqAmzeBXr2Av/5itfR82LWL\nzfvPmcNKJE210nVIc3fvsqUCjx5V/e9+8CAwezbrD6/KouGENCaU5A3Epk3Ad9+x+nEzM/XPU1LC\n+qanp7Or9969+YtRWzZuZA8wHTqk/MNL+fns7xYZycomCSGK0Zy8gZgxA2jbVrOnQePj2YNGbduy\nq1tjSPAAexCroICVPipDJmNdMj/4gBI8IXyjK3ktKiwEundnybpnT+Xf9++/bE3ZuDhWOTNkiPZi\n1JY//2RTS5cusfbMdVm+HDh9mv2caB6ekLrRlbwBadcO+PJLVm3z5Ily7zl5kn0wlJezKRpjTPAA\nMGIEe0CqvucG/vyT1cJv20YJnhBtoCt5LeM4Vjvu6lr31I1UysoP5XP5EyfqLkZtyc4G+vdn3Tpt\nbf/7emEhu0G9cydrqkYIqR/deDVAEgng5sba5CrqhX7pEmtwZmcH/Pyz4oRorIKCgDt3gM2bn98v\nk7Ey06FDgY8+0k9shBgjmq4xQEIhqzjx9wcePXq2v6qKtQPw8GAVNLGxDSvBA8CKFaw08vTp5/ev\nXs3KQLWxoAkh5Bm6ktehqVPZ1fr69Wxhj7ffBior2Zz0K6/oOzrt2bwZ+OUX4MQJVlJ5+DD7wEtL\nYx+AhBDl6exKfsmSJetcXV0z3dzc0n18fPbdu3evlfy1kJCQZU5OTtkuLi6XDx06NFzdMQxBYmIi\nb+f67js2/7xsGdCnDzBmDOtayUeC5zNOvr39NrvnsHMnEBWViOnTWc2/ISd4Q/55Vkdx8scYYlSH\n2kl++PDhhy5dutQ5PT3dzdnZOSskJGQZAGRkZIh37949JSMjQxwfHz8yMDAwvKqqyminhfj8D29t\nzebcjx8HjhxhrXT5aktsyP+DmpiwKpulS4FFixIxe7bhVw0Z8s+zOoqTP8YQozrUTr6enp4JJiYm\nVQDg7u5+Kj8/XwQAMTEx4/38/HaamZlVODg45Hbs2PFqSkpKX74CNnZjxrAuld266TsS3Ro4kD3o\nJBDQjVZCdImXK+zNmzfPGD16dBwAFBYWthOJRPny10QiUX5BQYEdH+MQ47ZlC+ucqetFVQhp1DiO\nq3UbNmxYQpcuXS7U3GJjY8fJj/nss8+W+/j4RMm/f//99zdu27btDfn3M2fO/CUqKsqn5rkBcLTR\nRhtttKm+1ZW3a2519jJMSEjwrOv1rVu3vh0XFzf6yJEjQ+X77OzsCvLy8uzl3+fn54vs7OwKar5X\nlbvDhBBC1KP2dE18fPzIdevWLYmJiRnftGnTx/L9Xl5esbt27ZoqlUrNc3JyHLOzs5369u2bwk+4\nhBBCVKF2V/IPPvhgo1QqNff09EwAgP79+58MDw8PFIvFGb6+vpFisTjD1NRUFh4eHigQCDj+QiaE\nEKI0VeZ2+NoOHjw4slOnTpc7duyYHRoaulQfMdS33bx5097Dw+MvsVh8qXPnzhc3bNgwV98x1bXJ\nZLIm3bt3Pzt27NgD+o6ltq20tNRi4sSJe11cXDJdXV0zTp482U/fMdXc1qxZs0wsFl/q0qXLBT8/\nvx2PHz9+Qd8xcRyHgICAzTY2NpIuXbpckO+7e/eu1bBhwxKcnJyyPD09D5WWlloYYpyLFy9e5+Li\nktmtW7f0CRMm7CsrK2tliHHKty+++GKRQCCounv3rpWhxvnNN9984OLiktm5c+eLQUFBYXWdQ+dB\ny2SyJh06dLiak5PjIJVKzdzc3M5lZGS46vuHWXMrKiqyPXv2bHeO4/DgwYMWzs7OVwwxTvm2fv36\nhdOmTds+bty4WH3HUts2ffr0iE2bNs3gOA4VFRWmhvCPvfqWk5Pj4OjoeF2e2H19fXdv3brVX99x\ncRyHv//+e1BaWlqP6v/YlyxZsjYsLCyI4ziEhoYuXbp0aaghxnno0CHPyspKE47jsHTp0lBDjZPj\n2MXdiBEj4h0cHHIMIckrivPo0aOvDxs2LEEqlZpxHIdbt261qescOn9IKSUlpW/Hjh2vOjg45JqZ\nmVVMnTp1V0xMzHhdx1EfW1vb4u7du58DgBYtWjx0dXXNLCwsbKfvuBTJz88XxcXFjZ41a9YvnIHe\n0L53716r48ePD5oxY8ZmADA1NZW1atXqnr7jqq5ly5b3zczMKsrLy5vLZDLT8vLy5oqKBvRh0KBB\nxy0tLUur74uNjfXy9/ePAAB/f/+I6Ohob/1E94yiOGt7pkafFMUJAAsXLvxy7dq1QfqISRFFcX7/\n/ffvLVu2LMTMzKwCANq0aXO7rnPoPMkXFBTY2dvb58m/N4Y6+tzcXIezZ8/2cHd3P6XvWBRZsGDB\nV+vWrVsi/4dkiHJychzbtGlzOyAgYEvPnj3T3nnnnZ/Ly8ub6zuu6qysrEoWLVq0vn379jfbtWtX\naGFhUTZs2LDD+o6rNhKJRCgUCiUAIBQKJRKJxIAbRTDVn6kxNDExMeNFIlF+t27dzus7lrpkZ2c7\n/f3336/169cv2cPDIzE1NbXONeN0nuSN7Sbsw4cPW0yaNGnvhg0b5rVo0eKhvuOp6ffffx9rY2Nz\nq0ePHmcN9SoeAGQymWlaWlrPwMDA8LS0tJ4vvvjiv6GhocH6jqu6a9eudfj666/n5+bmOhQWFrZ7\n+PBhi+3bt7+h77iUIRAIOEP/t/X5558vNzc3l06bNm2HvmOpqby8vPmaNWs+XL169Ur5PkP99yST\nyUxLS0stk5OT+61bt26Jr69vZF3H6zzJ16yjz8vLs6/+hKwhqaioMJs4cWLUm2++uc3b2zta3/Eo\nkpSUNCA2NtbL0dExx8/Pb+fRo0eHTJ8+/Vd9x1WTSCTKF4lE+X369DkNAJMmTdqblpamwqKI2pea\nmtp7wIABSdbW1ndNTU1lPj4++5KSkgboO67aCIVCSXFxsS0AFBUVtbWxsbml75hqI3+mxlA/NK9d\nu9YhNzfXwc3NLd3R0TEnPz9f1KtXrzO3bt2y0XdsNYlEonwfH599ANCnT5/TJiYmVXfv3rWu7Xid\nJ/nevXunZmdnO+Xm5jpIpVLz3bt3T/Hy8orVdRz14ThOMHPmzE1isThj/vz5X+s7ntqsWbPmw7y8\nPPucnBzHXbt2TR0yZMjRX3/9dbq+46rJ1ta22N7ePi8rK8sZAA4fPjysc+fOl/QdV3UuLi6Xk5OT\n+z169KgZx3GCw4cPDxOLxRn6jqs2Xl5esREREf4AEBER4W+oFyK1PVNjSLp27XpBIpEIc3JyHHNy\nchxFIlF+WlpaT0P84PT29o4+evToEADIyspylkql5tbW1ndrfYM+7hjHxcWNcnZ2vtKhQ4era9as\nWabvO9iKtuPHj78qEAiq3NzcznXv3v1s9+7dzx48eHCkvuOqa0tMTBxsyNU1586dc+vdu/dpQyql\nq7mFhYUFyUsop0+fHiGvYND3NnXq1J1t27YtNDMzk4pEorzNmzcH3L1712ro0KGHDamEsmacmzZt\nmtGxY8fs9u3b35D/O3rvvffCDSVOc3PzJ/KfZ/XXHR0drxtCdY2iOKVSqdmbb775W5cuXS707Nnz\nzF9//eVR1zn0tmgIIYQQ7TPaPu+EEELqR0meEEIaMEryhBDSgFGSJ4SQBoySPCGENGCU5AkhpAH7\nf7bO84WoPbF3AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x6730150>"
       ]
      }
     ],
     "prompt_number": 65
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, we would like to know if the data are linear or quadratic with x.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# fit a linear model y = XB + E\n",
      "X = np.hstack((np.ones((n,1)),x.reshape((n,1))))\n",
      "print(X)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[  1.   0.]\n",
        " [  1.   1.]\n",
        " [  1.   2.]\n",
        " [  1.   3.]\n",
        " [  1.   4.]\n",
        " [  1.   5.]\n",
        " [  1.   6.]\n",
        " [  1.   7.]\n",
        " [  1.   8.]\n",
        " [  1.   9.]\n",
        " [  1.  10.]\n",
        " [  1.  11.]\n",
        " [  1.  12.]\n",
        " [  1.  13.]\n",
        " [  1.  14.]\n",
        " [  1.  15.]]\n"
       ]
      }
     ],
     "prompt_number": 66
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# fit the data\n",
      "beta_h = lin.pinv(X).dot(y)\n",
      "resid = y - X.dot(beta_h)\n",
      "plt.plot(resid)\n",
      "print(R2(y, X))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.789569189143\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD9CAYAAABZVQdHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XdcVFf6P/DPKJDEuBZUQBgMRCQwYK9pBgXU+I1I0KCY\ngi3Zr8YkGmNLdjea70YwpmiK7ibRhDTLJq4Q14aFaGJBRWyIoMJKEYwFu9Lu74/zm4gIOOXeuTN3\nPu/Xa16rw8y9z0Z45vCc55yjkyQJRESkTY3UDoCIiJTDJE9EpGFM8kREGsYkT0SkYUzyREQaxiRP\nRKRhVif5srKyFsOHD/8xODj4qMFgyNq9e3fv8+fPu0dGRqYGBgbmDBgwYGNZWVkLOYIlIiLzWJ3k\nX3vttYWDBw9ee/To0eCDBw92CgoKyk5MTJwZGRmZmpOTExgeHr45MTFxphzBEhGReXTWLIa6ePFi\n865du+4/efLkgzWfDwoKyv7ll1+e8PT0LC0pKfEKCwtLy87ODrI6WiIiMotVI/m8vDz/Nm3a/D5m\nzJivunXrlvHiiy9+cfXq1ftLS0s9PT09SwHA09OztLS01FOecImIyBwu1ry5srLSJSMjo9unn346\nqWfPnnsmT568oHZpRqfTSTqd7o5fF+p6joiI7k6SJJ2pr7VqJK/X6wv1en1hz5499wDA8OHDf8zI\nyOjm5eVVUlJS4gUAp0+fbuvh4XGmnkDt/vH222+rHgPjZJyMkzEaH+ayKsl7eXmV+Pr6FuTk5AQC\nwKZNmyJCQkKODBky5OekpKR4AEhKSoqPjo5ebc19iIjIMlaVawDgk08+eeXZZ5/9vry83K19+/Yn\nvvrqqzFVVVWNY2NjVy5ZsmScn59f/sqVK2PlCJaIiMxjdZLv3LnzgT179vSs/fymTZsirL22PQgL\nC1M7BJMwTnkxTnk5QpyOEKMlrGqhtOrGOp2k1r2JiByVTqeDZKuJVyIism9M8kREGsYkT0SkYUzy\nREQaxiRPRKRhTPJERBrGJE9EpGFM8kREGsYkT0SkYUzyREQaxiRPRKRhTPJERBrGJE9EpGFM8kTk\n9G7eBObNA6qr1Y5EfkzyROT0MjOBmTOBH35QOxL5MckTkdM7cADo0gWYNQu4elXtaOTFJE9ETi8z\nExg9GujbV5RttIQnQxGR03v0UWDuXODBB8WIPiMDeOABtaOqm7knQzHJE5FTq64GmjcHTp0CWrYE\nZs8GsrOB5cvVjqxuDnX835Urat6diAg4cQJo1UokeACYPh3YsQP49Vd145KLqkk+NBRYt07NCIjI\n2R04AHTufOvvTZoAiYnA5MnaaKlUNcl/8QUwaRIQFweUlqoZCRE5q8xMUYevKS4OcHMDkpLUiUlO\nqib5yEjg0CExwdGxI7B0KcAyPRHZkrF9siadDli4EHjrLeDyZXXikovdTLxmZgIvvgg0bQr8859A\nYKAqYRGRk/H1BX75RXTW1BYfD3h7AwkJto+rPg418VpTly7Arl1AdDTwyCPA3/8OlJerHRURadm5\nc8ClS4CfX91fT0gQZeWTJ20alqzsJskDQOPGwGuviR7VXbuArl3FLDcRkRIOHAA6dQIa1ZMJvb2B\nKVNEx42jsqskb9SuHfDzz6JfdfhwYOJE4OJFtaMiIq2pa9K1ttdfB/btEyUdRyRLkq+qqmrctWvX\n/UOGDPkZAM6fP+8eGRmZGhgYmDNgwICNZWVlLcy9pk4HPPMMcOQIUFkJhIQA//63HNESEQl1TbrW\ndt99wHvviSpDVZVt4pKTLEl+4cKFrxkMhiydTicBQGJi4szIyMjUnJycwPDw8M2JiYkzLb12y5bA\n55+L3eFmzQKefhooLJQjaiJydpmZt/fI12f4cKBZM9EB6GisTvKFhYX6tWvXDh4/fvyXxhnflJSU\nqPj4+CQAiI+PT1q9enW0tffp2/fWooWuXYHPPnPMT1Uisg/l5UBurqgS3I1OByxYAPztb45XOrY6\nyU+ZMuWj+fPnT2vUqNEfa8NKS0s9PT09SwHA09OztLS01NPa+wDAPfeIOv0vvwDLlgGPPSb67ImI\nzJWVBfj7i3KMKbp1AwYPFp1/jsTFmjevWbPmKQ8PjzNdu3bdn5aWFlbXa3Q6nWQs49Q2e/bsP/4c\nFhaGsLA6L3EHgwHYtk20NvXvD7z0EvCXv5j+j0VEZMqka23vviu2Y3npJaBDB2Xiqi0tLQ1paWkW\nv9+qxVBvvvnm3G+//fZ5FxeXyhs3btx76dKlZjExMav27NnTMy0tLczLy6vk9OnTbfv167c1Ozs7\n6LYby7QLZXGxmBDJzBRJ38TPCSJyclOmAG3bmt8eOW+eaO1OTlYmrrux6WKouXPnvllQUOCbl5fn\nv3z58pH9+/ff8u233z4fFRWVkpSUFA8ASUlJ8dHR0autuU9DvL2Bf/0L+OADMSl77pxSdyIiLbFk\nJA+IjcsOHwY2bZI/JiXI2idvLMvMnDkzMTU1NTIwMDBny5Yt/WfOnJko533qEhUlHkuWKH0nInJ0\nknTn7pOmuuce4P33RbKvrJQ/NrnZzd41cti3D4iJEftDu1g120BEWnbqFNC7N3D6tGXvlyQgPPzW\nYk1bcti9a+TQvTug1wMpKWpHQkT2zNJSjZGxpXLOHODCBfniUoKmkjwAvPoq8PHHakdBRPbM0lJN\nTZ06iXnAd96RJyalaC7Jx8QAx4+Lf0QiorpYO5I3+r//A777TpwJa680l+RdXYEJE4BPPlE7EiKy\nV6bsWWOKNm3EditTp1p/LaVoauLV6PffxaEjublA69aK3IKIHNSlS6I//tIlsb25tcrLxQKpjz8G\nBg2y/np349QTr0Zt2ojDR778Uu1IiMjeHDokkrIcCR4QZ8F+8IFYXFVRIc815aTJJA8Ar7wiNjFz\nhD5WIrIdU3eeNMdTT4ljBBcvlve6ctBsku/WTRzptVqxtbZE5IjkmnStSacDPvpIbF5mb6vuNZvk\nAbZTEtGd5GifrEtICBAbC7z9tvzXtoYmJ16NKirECewpKWIPeiJybpWVQPPmQEkJ8Kc/yX/9c+eA\n4GBg61bT9qm3BCdea3B1FUuO2U5JRIDouPP2VibBA0CrVmLb8ylTxNYH9kDTSR4AXnxRnA37++9q\nR0LknKZNA65fVzsKQYlJ19omTBBHlK5Zo+x9TKX5JN+6tVgF+8UXakdC5HwuXhQ7Nv76q9qRCEpM\nutbm6gp8+KFYIFVeruy9TKH5JA+IdspFi+yzh5VIy44dE/9rL3uvKzXpWtugQeLkqE8/Vf5ed+MU\nSb5LF6B9e1G2ISLbOXZMtDLbS5K3xUje6MMPgdRU9WvzTpHkAbZTEqkhOxt49lmxaeDZs+rGUlIi\nfpvX621zv4ceAtatEz30anKaJD90qDgoYN8+tSMhch7Z2UDHjkDfvqKtUE3GUo3aSdfWnCbJu7gA\nL7/MdkoiW8rOBoKCgIgI9Us2cu086Wic6pC88eOBgADgzBnAw0PtaIi0rbISOHlSTEC6uKhfLs3M\nBAYOVDcGNTjNSB4QCxWGDwc+/1ztSIi0Lz8f8PQEmjQBDAbg2jWR9NViy0lXe+JUSR4Q7ZSLF7Od\nkkhpxlINIOrgapZsrl8H8vLElgPOxumSfKdO4kCRn35SOxIibauZ5AF1k/zhw6Lbxc1NnfuryemS\nPMB2SiJbOHbs9iQfHg5s2QJUV9s+FmeddAWcNMkPGQIUFwN79qgdCZF2ZWeL0bORXi9ObcvMtH0s\nttizxl45ZZJnOyWR8mqXawAgMlKdko2zTroCTprkAWDcOODnn8UqOCKS17lzYnMuL6/bn1ejLl9d\nDRw8yJG803F3F6e4sJ2SSH7Genzt1aVPPAHs3AncuGG7WPLygBYtxM+8M3LaJA+Idsp//MM+tgMl\n0pLa9Xij5s3FNgc7dtguFmeedAWcPMmHhoq+2R9/VDsSIm2pqx5vZOuSjTNPugJWJvmCggLffv36\nbQ0JCTkSGhp6+OOPP34VAM6fP+8eGRmZGhgYmDNgwICNZWVlLeQJV35spySSnz0leY7kreDq6lrx\n0UcfTTly5EjIrl27+nz22WcvHz16NDgxMXFmZGRkak5OTmB4ePjmxMTEmXIFLLenngJKS4Hdu9WO\nhEg7avfI19Snj/gQOH/eNrFwJG8FLy+vki5dumQCQNOmTa8EBwcfLSoq8klJSYmKj49PAoD4+Pik\n1atXR8sRrBIaNwYmTWI7JZFcysuB//5XHNRTFzc34LHHbLP18PnzwIULwIMPKn8veyXbLpT5+fl+\n+/fv79q7d+/dpaWlnp6enqUA4OnpWVpaWupZ13tmz579x5/DwsIQFhYmVzhmGTtWfBOcPg20batK\nCESaceIE4OsL3HNP/a8xlmyGDVM2loMHxURvIweefUxLS0NaWprF79dJMpxNdeXKlaZPPPHEL3/9\n61//Lzo6enXLli0vXLhwoaXx6+7u7ufPnz9/WwOTTqeT5Li3XCZMENsPz5mjdiREju3f/waWLhXr\nUOpz8KBI8Lm5ysayYIG4x2efKXsfW9LpdJAkyeSjT6z+fKuoqHAdNmzYT88///y30dHRqwExei8p\nKfECgNOnT7f18PA4Y+19lPbKK8A//wncvKl2JESOraF6vFHHjsDly2I7YiU5+6QrYGWSlyRJN27c\nuCUGgyFr8uTJC4zPR0VFpSQlJcUDQFJSUrwx+dszg0F8461cqXYkRI6tvh75mnQ6sWHZ5s3KxuLs\nk66AleWaX3/99bG+fftu69Sp00GdTicBQEJCwqxevXqlx8bGrjx16lQ7Pz+//JUrV8a2aNGi7LYb\n21m5BhC/Xs6ZIzYuc7ZzIInk0qcP8P77YnK1IV99BWzcCCxbpkwc5eVipevZs+LgEq0wt1wjS03e\nEvaY5KuqxF7z330HPPyw2tEQOR5JAlq2BI4fB1q3bvi1p04BPXqI/aOUmBg9eBAYMQI4elT+a6vJ\n5jV5LTG2U3JxFJFlzpwRu7zeLcEDQLt24gPh0CFlYnHmnSdrYpKvZcwYYP16oKhI7UiIHE9DK13r\nouTqV066CkzytbRoAYwaJTYuIyLzmDLpWlNEBJCaqkwsnHQVmOTrMGmS2ILYltuhEmmBuSP5sDDg\nt9/kb12WJI7kjZjk6xAcLL45VqxQOxIix2JKj3xNLVuK9uWdO+WNo6hIzLHVPrTEGTHJ1+PVV4GF\nC8WIgIhMY+5IHlCmLs9J11uY5Ovx5JPApUu2PdyAyJFdvw4UFwP+/ua9T4lzXw8cYD3eiEm+Ho0a\nia0O2E5JZJrcXLHRn4uZ2x4+/DCQlQWUld39tabiSP4WJvkGjB4NbNggDiUmooaZW483uuce4JFH\nACs2WrwDJ11vYZJvQPPmwODByi27JtISS+rxRnLW5S9fBgoLxep1YpK/q9Gjga+/VjsKIvtnbo98\nTXIm+UOHgJAQ88tGWsUkfxfh4WJvDaWWXhNphTUj+U6dRFm0oMD6ODjpejsm+bto3Bh44QUgKUnt\nSIjslyQBOTmWj+QbNRIDKjlG85x0vR2TvAni48XOlJWVakdCZJ+KioCmTcW2IJaSq2TDSdfbMcmb\n4KGHRGvYhg1qR0Jkn6ypxxsZk7w1CxCrqoDDh0X5hwQmeRPFx3MClqg+1tTjjfz8gD/9SSRpSx0/\nDnh6As2aWReLljDJm2jECHGKDXvmie5kaY98bdaufuXOk3dikjdRixaiZ375crUjIbI/cozkAevr\n8px0vROTvBnYM09UNzlq8gDQrx/w66/ifFZLsH3yTkzyZoiIEBswHTmidiRE9uPKFVHGbNfO+mu5\nu4uVqrt3W/Z+juTvxCRvBvbME90pJwfo0EH8fMjB0pLNmTNiJ0w5Pmy0hEneTPHxwLffsmeeyEiu\neryRpUneWKrR6eSLRQuY5M0UFAQ88IDotCEi+erxRo8+Chw8CFy8aN77WKqpG5O8BTgBS3SL3CP5\ne+8F+vQBfvnFvPdx0rVuTPIWGDFCrH49f17tSIjUJ1ePfE2WlGw4kq8bk7wFWrYUxwOyZ56cXVWV\nOBFK7r3bzU3yN24AJ0+KQ8HpdkzyFho9ml02RKdOAa1bi83J5NS1q+iWKSoy7fVHjgABAeKUKbod\nk7yFIiPF6TNZWWpHQqQeuSddjRo1Avr3BzZvNu31LNXUT7Ekv379+kFBQUHZHTp0yJ03b94Mpe6j\nlsaNgeef52ie1HPokPqtvErU443MKdlw0rV+iiT5qqqqxpMmTfp0/fr1g7KysgzLli2LO3r0aLAS\n91ITe+ZJLUVFQO/ewOrV6sYhd2dNTeZsPcyRfP0USfLp6em9AgICjvv5+eW7urpWjBw5cnlycvJQ\nJe6lpuBgwNcXSE1VOxJyNn/7G9CqlenlDKUomeQffFC0Ux492vDrJEn01XMkXzdFjrotKiry8fX1\n/eO0Rr1eX7h79+7etV83e/bsP/4cFhaGsLAwJcJRlLFn/skn1Y6EnMXBg8CaNcC//gWMH69uLErV\n5I2Mo/mGumby88XEb+vWysWhprS0NKSlpVn8fkWSvE6nM+lsl5pJ3lGNHAnMmgVcuCBaK4mUNn06\n8Je/AI89BpSViQ4XNfZrKSsDrl4FfHyUu0dEhCiJvvpq/a/R+nF/tQfAc+bMMev9ipRrfHx8igoK\nCnyNfy8oKPDV6/WFStxLbS1bAgMHAitWqB0JOYONG4ETJ4A//9n8DhS5HTsmRvFK7hXTvz+wbRtQ\nUVH/a3hQSMMUSfI9evTYm5ub2yE/P9+vvLzcbcWKFSOioqJSlLiXPeA2B2QLVVXAtGnAvHmAm5t4\nTq7Dry2hZD3eqHVroH17ID29/tdw0rVhiiR5FxeXyk8//XTSwIEDNxgMhqwRI0asCA4Ovsv0ieOK\njBS/Mt9tgojIGt98I85AffrpW8+Fh4uRvDWHX1tK6Xq80d0+yNg+2TCdpMZ3B0TdXq17K2H6dPHr\nc2Ki2pGQFl29KhLqTz+J1sma2rcHkpOB0FDbxhQTA8TFAc88o+x9UlOBd94Btm+/82tlZYBeL3as\nlGs/e3un0+kgSZLJRTKueJWJsWe+qkrtSEiLPvxQbMFbO8EDt0bztmaLcg0gJpgzM4HLl+/82sGD\nQMeOzpPgLcEkL5OQENFloFZ9lLSrpARYsABISKj762rU5SsqxIZgHToof6/77gN69RITsLWxHn93\nTPIy4gQsKWH2bPG99eCDdX+9X7+7d6DILS8P8PYWi5Vsob4PMib5u2OSl9HIkcC6daJOSCSHrCxg\n1Srgrbfqf02bNuIDYM8e28Wl5J41dakvyXPS9e6Y5GXk7g4MGMCeeZLPjBnAzJnie6shtq7L26oe\nb9Stm9ivp6Tk1nMVFaKjrWNH28XhiJjkZcaSDcllyxaxT/rLL9/9tbauy9s6yTduLMpSNf8/Hjsm\n9o66/37bxeGImORlNmCA2Evj2DG1IyFHVl0NvPGGaMk15SCMxx8H9u0TrZa2YKse+Zpqf5CxHm8a\nJnmZubgAzz3HfebJOt9/L1a1mtqDfv/9QPfudfeSK8HWNXngzq2Htb5njVyY5BUQHy9WJzpyz3xV\nFVBernYUzun6dTHR+v775u0LY6u6/Nmz4vvDw0P5e9UUECDKNsbfkrlnjWmY5BUQGgq0bav+Xt+W\n+u03ccbmuHFqR+KcFi4EevYUi4DMYaskb6zHK7kxWV10OrGFiHE0z3KNaZjkFeKIE7Bnz4rEPmIE\nMGkS8PPPYlRJtvP772IEb8n2GL16iR0qz56VP66a1KjHGxlLNqdPi7+3batOHI6ESV4hI0cCa9c6\nRs98dTXw5Zdi1W6zZqI3+6WXgB49RN8/2c477wCjRlm2ktTVVUzAbt0qf1w1qVGPN+rfH0hLA/bu\nFaN4W/824YiY5BXSqpX41XLlSrUjadiBA6IssGQJsGED8NFHItEDYtLvX/9SNz5ncuwYsGyZONrP\nUrZopbR1+2RNHh6An5/4fmU93jRM8gqKj7ffLpvLl4HXXxcfRKNHizp87frm00+LkTxLNrYxc6bY\nzdSaY+xsUZdXM8kD4oMsJYX1eFMxySto4EBRI83JUTuSWyQJ+PFHcWbmhQtisc1LL4ltkmvz8BBt\neevX2z5OZ7N9O7B/f8PH3JkiNFR8gOfnyxLWHW7eBAoK6t9HxxYiIsT/MsmbhkleQa6u9tUzf+IE\nMHiw2PDqhx+Ar74S+5405Jln7L/k5Oiqq4GpU4F337V+wy+dTtnR/IkTwAMP3DqZSg2PPy62OVBr\n8tfRMMkrzB565m/eFBN6vXuLiav9+8UPiiliYliyUdrKlSLRx8XJcz0l6/Jql2oAsfBr3z4xiKK7\nY5JXWMeOgKen2IdEDampIobMTCAjQ5wRas4PB0s2yrp5E5g1S7RN1lUys0R4uPh+q66W53o12UOS\nJ/MwyduAGj3zxcWijfOll8SpQqtWAe3aWXYtdtko59NPxYdwWJh813zgAdEhdfiwfNc0YpJ3PEzy\nNhAXB/znP+IcSqVVVooVk507i2XgR44ATz1l3TVjYkTPP0s28jp3Tix6mjdP/msrVZc/doy1cEfD\nJG8DrVqJHzqlR8O7d4tVj8nJolvj738HmjSx/roeHmKiiyUbef397+K3pOBg+a8dHi5/XV6S1F3t\nSpZhkrcRpUo21dVitP6//yv62qdOFSM4uX+lZslGXidOiIPfZ89W5vr9+gG//irvkYAlJWLb41at\n5LsmKY9J3kYGDQKOHwdyc627TlmZWJk6e7bow3d3B4YOFSP2rCzg2WeVWerNko28Zs0CpkxRbifH\n1q2B9u3Fb3dyYT3eMTHJ24irq0jA5vTMV1eLxL1kCTB+vNhbxtcXSEgQXRkvvywWWh0/LiZXW7RQ\nLn5PT1Gy2bBBuXs4i507xWPKFGXvExEhb12e9XjHpJOMO/Db+sY6naTWvdVy8KCYBM3Pr7tdrqxM\njLx27gR27RJ/dncHHn741qNTJ3EwiRoWLxYlgO+/V+f+WiBJwKOPAn/+s1hDoaQNG0TdX66DRCZP\nFoOMqVPluR5ZRqfTQZIkk39fZ5K3se7dgffeEzXT7Oxbo7qdO4FTp8TXjQm9Tx/bH8zQkNJSMZI7\nfRq47z61o3FMP/0kEu/eveIADCVdvSp+AyspAZo2tf56gwYBr7wC/M//WH8tspy5SV6lMaHzGj1a\nPK5cuX2UPnGi6Je251V8np7iMJENG4DoaLWjcTzl5cCMGcA//qF8ggfEytAePYBt28R2FtZSc4th\nshyTvI2NGycmxLp3F0nT0cTGii4bJnnzLV4MBAbe2mDLFox1eWuT/LVr4jcCPz9ZwiIbsnjiddq0\nafODg4OPdu7c+UBMTMyqixcvNjd+LSEhYVaHDh1yg4KCsjdu3DhAnlC1oUkT8QPniAkeuNVlc+OG\n2pE4lrIysQHZe+/Z9r5yLYrKzRWDE1v8BkLysjjJDxgwYOORI0dCDhw40DkwMDAnISFhFgBkZWUZ\nVqxYMSIrK8uwfv36QRMnTlxUXV3NLh6N8PQUW7yyy8Y0kgQcOiQmWocOFVsB21LPnmKi/8wZ667D\n9knHZXHyjYyMTG3UqFE1APTu3Xt3YWGhHgCSk5OHxsXFLXN1da3w8/PLDwgIOJ6ent5LroBJfVwY\n1bDqaiA9XdTfAwNFR5W3t2XntlrLxQXo29f6IwFZj3dcstTkly5dOjYuLm4ZABQXF3v36dNnl/Fr\ner2+sKioyKeu982usdwvLCwMYXLu0kSKiYkB3nxTlGys3f9cKyorRXvpqlXAv/8tJj2HDQOWLxfr\nC9Q8i9S49fCIEZZfIzsbePJJ+WIi06WlpSEtLc3i9zeY5CMjI1NLSkq8aj8/d+7cN4cMGfIzALz7\n7rtvubm5lY8aNeqH+q6j0+nq7JWcrdSablKUl9etLpuhQ9WORj03b4p696pV4jg6X1/xAbhxozL7\n0VgqPBxYsMC6a2RnK794i+pWewA8Z84cs97fYJJPTU2NbOjrX3/99ei1a9cO3rx5c7jxOR8fn6KC\nggJf498LCwv1Pj4+RWZFRXbPWLJxtiR/9arYqG3VKjEBHRIiRux/+Yv9dp4YDGI7ipMnLTu2r7pa\nrKzmalcHJUmSRY9169YNMhgMR37//ffWNZ8/cuSIoXPnzpk3b950O3nypP+DDz54orq6Wlf7/eLW\n5KhOn5akFi0k6fp1tSNR3vnzkvTNN5IUHS1JzZpJUmSkJC1eLEnFxWpHZrpnn5Wkzz+37L2nTkmS\nt7e88ZDl/n/uNDlXW1yTf+WVVz4pLy93i4yMTAWAhx9+eOeiRYsmGgyGrNjY2JUGgyHLxcWlctGi\nRRPrK9eQ4/LyEnvWb9wIREWpHY38SkvFls2rVonVyP36iVLMkiViEZujCQ8Xv4G8+KL57+X2wo6N\n2xqQxT77TCTA775TOxL5nDwp9mjZvl1MNMbEiOX8cmwLoCbjlhmlpeYfM/jJJ8DRo8CiRcrERuYx\nd1sD9q+TxYYNEydeaWFh1M2bYrFSr17AY4+J/Xl++AEYPtzxEzwgjn5s2VJskmcu9sg7NiZ5sljN\nko0j27pV/P9ITwf27QOmT9dma6ilWw+zR96xMcmTVRx5YVRpKfD882LDuHnzRA3+gQfUjko5lh4J\nyJq8Y2OSJ6sMGwasWeNYJZvqauCf/xS7frZtK45PdIZW0H79gN9+E7thmuryZeDCBbEGgBwTkzxZ\nxctLHGTiKCWbzEzgkUfE+aqbN4sNw7RQczeFu7vYZmHXrru/1ujYMfEecydryX7wn46s5gglm8uX\nxYrNgQNFG+G2bWIk72zMrcuzHu/4mOTJasaSzc2bakdyJ0kCfvxRrPq8eFGUZsaNc96Rqbl1edbj\nHZ+TfquTnNq2tc+SzcmT4qi62bNFO+TSpUDr1mpHpa5HHwUOHAAuXTLt9WyfdHxM8iQLeyrZ1Ox5\nf+IJICMDePxxtaOyD02aiP8u27aZ9nomecfHJE+yGDYM+Pln9Us2xp73XbvEYdkzZgBuburGZG9M\nPS2qqgo4cUJMvJLjYpInWbRtKyYy1SrZGHve4+PF4RwpKfa7K6TaTJ18/e9/AQ8PMfonx8UkT7Ix\nHvJtSzXm0RYkAAAN6UlEQVR73r28gKwscci4mod02Lvu3cVeNqWlDb+Ok67aIMvJUESAKNn87W+i\nZHPPPcrfr7JS7C1z5ozoGOnUSfl7aoGLCxAWBmzZAsTF1f861uO1gSN5kk3btuKg6tRU5e9VXQ2M\nHy9Wb/7yCxO8uUxppWSPvDYwyZOsbNFlI0nAtGlAbq7ogXd1VfZ+WmQ897Wh3b45ktcGJnmS1bBh\nYtJTyS6befPEBO+aNZwUtFRQEFBRIbpn6sOavDYwyZOsvL2VLdl88QXw+efiEPGWLZW5hzPQ6Rru\nsrlwQZwL27atbeMi+THJk+yU6rL56Sfg7bfFKN7bW/7rO5uG6vLGejy7lBwfkzzJTomFUZs3AxMm\nAGvXAgEB8l3XmYWHi8Vj1dV3fo31eO1gkifZeXsDISGWHVBRl717Ravfjz8CXbrIc00C9Hqxl09m\n5p1fYz1eO5jkSRHPPAOsXGn9dbKzgSFDgC+/BPr2tf56dLv6tjjgSF47mORJEXKUbAoKxP7viYlA\nVJR8sdEtxlbK2tgjrx1M8qQIHx/rSjZnzwIDBgCvvir2oyFlhIUBO3bc/mFcUQHk53PuQyuY5Ekx\nli6MunwZGDxY7EEzdar8cdEtLVsCwcHAzp23njt5UnxI22JrClIekzwpxrgwypyDo2/eBGJixHbB\nc+cqFxvdUrsuz3q8tjDJk2J8fMSxe6YujKqqEtsFN2sG/OMf7NG2ldp1edbjtYVJnhRlaslGkoCX\nXwbOnQO+/x5o3Fj52Eh45BHg0CFxBi7AkbzWMMmTooYPN61k89e/in741auBe++1TWwk3Hcf0KeP\n2M0TYI+81lid5D/44IOpjRo1qj5//ry78bmEhIRZHTp0yA0KCsreuHHjAGvvQY7LWLJpqMtmwQIx\n2l+3DvjTn2wXG91irMtLEkfyWmNVki8oKPBNTU2NfOCBB/5rfC4rK8uwYsWKEVlZWYb169cPmjhx\n4qLq6mr+xuDEGirZfPst8OGHom7fpo1t46JbjJuVnT0r5kJat1Y7IpKLVcn39ddf//C9996bXvO5\n5OTkoXFxcctcXV0r/Pz88gMCAo6np6f3si5McmT1ddn85z9iX/j164F27dSJjYRu3YCiIrGXDTcm\n0xaLk3xycvJQvV5f2KlTp4M1ny8uLvbW6/WFxr/r9frCoqIiH2uCJMem14vEUbNks307MHo0kJws\nyjmkrsaNgX79gEWLWI/XmgbPeI2MjEwtKSnxqv38u++++1ZCQsKsmvV2SZLq/ezX6XR1nj8ze/bs\nP/4cFhaGsLAwE0ImR2TcfnjwYODAATG6//57oHdvtSMjo/BwYNIk8W9E9iMtLQ1paWkWv18nNXT+\nVz0OHz4cGh4evrlJkybXAKCwsFDv4+NTtHv37t5fffXVGACYOXNmIgAMGjRo/Zw5c97u3bv37ttu\nrNNJltybHFNhoVjg9NtvIpl89JFI/GQ/srPF6tfkZO4VZM90Ol2Dg+o7Xi9HovX398/bt29fd3d3\n9/NZWVmGUaNG/ZCent6rqKjIJyIiYtPx48cDao/mmeSdz6OPAocPiw3HJkxQOxqqTZJE6WzdOsDP\nT+1oqD7mJvkGyzVm3PSPbG0wGLJiY2NXGgyGLBcXl8pFixZNrK9cQ85l6lQxomeCt086HXD0qNpR\nkNxkGclbdGOO5ImIzGbuSJ7960REGsYkT0SkYUzyREQaxiRPRKRhTPJERBrGJE9EpGFM8kREGsYk\nT0SkYUzyREQaxiRPRKRhTPJERBrGJE9EpGFM8kREGsYkT0SkYUzyREQaxiRPRKRhTPJERBrGJE9E\npGFM8kREGsYkT0SkYUzyREQaxiRPRKRhTPJERBrGJE9EpGFM8kREGsYkT0SkYUzyREQaxiRPRKRh\nViX5Tz755JXg4OCjoaGhh2fMmDHP+HxCQsKsDh065AYFBWVv3LhxgPVhqictLU3tEEzCOOXFOOXl\nCHE6QoyWsDjJb926tV9KSkrUwYMHOx0+fDj0jTfeeB8AsrKyDCtWrBiRlZVlWL9+/aCJEycuqq6u\ndtjfGBzlH55xyotxyssR4nSEGC1hcfJdvHjxhFmzZiW4urpWAECbNm1+B4Dk5OShcXFxy1xdXSv8\n/PzyAwICjqenp/eSK2AiIjKdxUk+Nze3w7Zt2/r26dNnV1hYWNrevXt7AEBxcbG3Xq8vNL5Or9cX\nFhUV+cgRLBERmUmSpHofERERqaGhoYdqP5KTk6NCQ0MPvfrqqwslSUJ6enpPf3//k5IkYdKkSZ98\n9913zxqvMW7cuC9/+umnmNrXBiDxwQcffPBh/qOhvF374YIGpKamRtb3tcWLF0+IiYlZBQA9e/bc\n06hRo+qzZ8+29vHxKSooKPA1vq6wsFDv4+NTVPv9kiTpGro3ERFZz+JyTXR09OotW7b0B4CcnJzA\n8vJyt9atW5+NiopKWb58+cjy8nK3vLw8/9zc3A69evVKly9kIiIyVYMj+YaMHTt26dixY5d27Njx\nkJubW/k333zzAgAYDIas2NjYlQaDIcvFxaVy0aJFE3U6nSRfyEREZDJzajtyPdatWzfooYceyg4I\nCMhNTEycoUYMd3ucOnXKNywsbKvBYDgSEhJyeOHCha+qHVNDj8rKysZdunTZ/9RTT/2sdiz1PS5c\nuNBi2LBhPwYFBR0NDg7O2rlzZx+1Y6r9mDt37iyDwXAkNDT0UFxc3A83bty4R+2YJEnCmDFjlnp4\neJSGhoYeMj537tw594iIiNQOHTrkREZGbrxw4UILe4zzjTfemB8UFHS0U6dOB55++ulVZWVlze0x\nTuPj/fffn6rT6arPnTvnbq9xfvzxx68EBQUdDQkJOTx9+vR5DV3D5kFXVlY2bt++/fG8vDy/8vJy\n186dO2dmZWUFq/0fs/bj9OnTXvv37+8iSRIuX77cNDAw8Jg9xml8fPDBB6+PGjXq+yFDhqSoHUt9\njxdeeCFpyZIlYyVJQkVFhYs9/LDXfOTl5fn5+/ufNCb22NjYFV9//XW82nFJkoRt27Y9npGR0bXm\nD/u0adPemzdv3nRJkpCYmDhjxowZifYY58aNGyOrqqoaSZKEGTNmJNprnJIkBncDBw5c7+fnl2cP\nSb6uOLds2dIvIiIitby83FWSJJw5c6ZNQ9ew+SKl9PT0XgEBAcf9/PzyXV1dK0aOHLk8OTl5qK3j\nuBsvL6+SLl26ZAJA06ZNrwQHBx8tLi72VjuuuhQWFurXrl07ePz48V9KdjqhffHixebbt29/fOzY\nsUsBwMXFpbJ58+YX1Y6rpmbNml1ydXWtuHbtWpPKykqXa9euNamraUANjz/++PaWLVteqPlcSkpK\nVHx8fBIAxMfHJ61evTpanehuqSvOyMjI1EaNGlUDQO/evXcXFhbq1YnulrriBIDXX3/9w/fee2+6\nGjHVpa4461ujVB+bJ/mioiIfX1/fAuPfHaGPPj8/32///v1de/fuvVvtWOoyZcqUj+bPnz/N+INk\nj/Ly8vzbtGnz+5gxY77q1q1bxosvvvjFtWvXmqgdV03u7u7np06d+kG7du1OeXt7F7do0aIsIiJi\nk9px1ae0tNTT09OzFAA8PT1LS0tLPdWO6W6WLl06dvDgwWvVjqMuycnJQ/V6fWGnTp0Oqh1LQ+pb\no1Qfmyd5R5uEvXLlStPhw4f/uHDhwteaNm16Re14aluzZs1THh4eZ7p27brfXkfxAFBZWemSkZHR\nbeLEiYsyMjK63X///VcTExNnqh1XTSdOnGi/YMGCyfn5+X7FxcXeV65cafr9998/q3ZcptDpdJK9\n/2y9++67b7m5uZWPGjXqB7Vjqe3atWtN5s6d++acOXPeNj5nrz9PlZWVLhcuXGi5a9euPvPnz58W\nGxu7sqHX2zzJ1+6jLygo8K25QtaeVFRUuA4bNuyn55577rvo6OjVasdTlx07djySkpIS5e/vnxcX\nF7dsy5Yt/V944YVv1I6rNr1eX6jX6wt79uy5BwCGDx/+Y0ZGRje146pp7969PR555JEdrVq1Oufi\n4lIZExOzaseOHY+oHVd9PD09S0tKSrwA4PTp0209PDzOqB1Tfb7++uvRa9euHWyvH5onTpxon5+f\n79e5c+cD/v7+eYWFhfru3bvvO3PmjIfasdWm1+sLa69ROnfuXKv6Xm/zJN+jR4+9ubm5HfLz8/3K\ny8vdVqxYMSIqKirF1nHcjSRJunHjxi0xGAxZkydPXqB2PPWZO3fumwUFBb55eXn+y5cvH9m/f/8t\nxnZWe+Ll5VXi6+tbkJOTEwgAmzZtiggJCTmidlw1BQUFZe/atavP9evX75MkSbdp06YIg8GQpXZc\n9YmKikpJSkqKB4CkpKR4ex2IrF+/ftD8+fOnJScnD7333ntvqB1PXTp27HiotLTUMy8vzz8vL89f\nr9cXZmRkdLPHD8661ii1atXqXL1vUGPGeO3atU8GBgYea9++/fG5c+fOUnsGu67H9u3bH9PpdNWd\nO3fO7NKly/4uXbrsX7du3SC142rokZaW9oQ9d9dkZmZ27tGjxx57aqWr/Zg3b950YwvlCy+8kGTs\nYFD7MXLkyGVt27YtdnV1Ldfr9QVLly4dc+7cOffw8PBN9tRCWTvOJUuWjA0ICMht167df40/RxMm\nTFhkL3G6ubndNP73rPl1f3//k/bQXVNXnOXl5a7PPffct6GhoYe6deu2b+vWrWENXUMnSXZdxiMi\nIis47D7vRER0d0zyREQaxiRPRKRhTPJERBrGJE9EpGFM8kREGvb/AP22Iwb+X4IqAAAAAElFTkSu\nQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x6542e10>"
       ]
      }
     ],
     "prompt_number": 67
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X2 = np.hstack((X, (x*x).reshape((n,1))))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 69,
       "text": [
        "array([[   1.,    0.,    0.],\n",
        "       [   1.,    1.,    1.],\n",
        "       [   1.,    2.,    4.],\n",
        "       [   1.,    3.,    9.],\n",
        "       [   1.,    4.,   16.],\n",
        "       [   1.,    5.,   25.],\n",
        "       [   1.,    6.,   36.],\n",
        "       [   1.,    7.,   49.],\n",
        "       [   1.,    8.,   64.],\n",
        "       [   1.,    9.,   81.],\n",
        "       [   1.,   10.,  100.],\n",
        "       [   1.,   11.,  121.],\n",
        "       [   1.,   12.,  144.],\n",
        "       [   1.,   13.,  169.],\n",
        "       [   1.,   14.,  196.],\n",
        "       [   1.,   15.,  225.]])"
       ]
      }
     ],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# fit the data\n",
      "beta_h = lin.pinv(X2).dot(y)\n",
      "resid = y - X2.dot(beta_h)\n",
      "plt.plot(resid)\n",
      "print(R2(y, X2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.895987535281\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD9CAYAAABZVQdHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXtYlGX6x7/DyUOUeGJAB4UEQtAEMyWslhUhTSXU1sRK\n1sPullubHSy13dLdDTHXFCutLU06qbVtQKkkHuiwpmZq7YoKurAyHEZFMA/ocHh/f9y/NxEHmMN7\nHO7Pdc0FzMz7PLfIfOee73M/92MQBAEMwzCMe+KhdgAMwzCMfLDIMwzDuDEs8gzDMG4MizzDMIwb\nwyLPMAzjxrDIMwzDuDGSiHxjY6NnTEzMwQkTJnwGAGfPnu2RmJiYHx4eXpSUlLSttrbWT4p5GIZh\nGMeQROQzMzOfiIyMLDQYDAIAZGRkzE9MTMwvKioKT0hI2JGRkTFfinkYhmEYx3BZ5M1ms2nLli33\nzp49+21BEAwAkJubm5yWlpYFAGlpaVnZ2dkprs7DMAzDOI7LIv/kk0+uWLZs2TwPD48m8T6LxWI0\nGo0WADAajRaLxWJ0dR6GYRjGcbxcufjzzz8f7+/vfyomJuZgQUFBvK3nGAwGQbRxWt7vytwMwzAd\nFdE1sQeXMvndu3fH5ebmJoeEhJSkpqZu2Llz56iHH374PaPRaKmqqgoAgMrKykB/f/9TrQSq+duL\nL76oegwcJ8fJcXKM4s1RXBL59PT0hWVlZUElJSUhGzdunDpq1Kid77333sPJycm5WVlZaQCQlZWV\nlpKSku3KPAzDMIxzSFonL1ow8+fPz8jPz08MDw8v2rlz56j58+dnSDkPwzAMYx8uefLN+cUvfvHl\nL37xiy8BoEePHme3b98+Wqqx1SQ+Pl7tEOyC45QWjlNa9BCnHmJ0BoMzHo8kExsMglpzMwzD6BWD\nwQBBqYVXhmEYRtuwyDMMw7gxLPIMwzBuDIs8wzCMG8MizzAM48awyDMMw7gxLPIMwzBuDIs8wzCM\nG8MizzAM48awyDMMw7gxLPIMwzBuDIs8wzCMG8MizzAM48awyDMMw7gxLPIMwzBuDIs8wzCMG8Mi\nzzAM48awyDMMw7gxLPIMw8jKhAnAxYtqR9Fx4TNeGYaRjbo6oGtX4McfgcGD1Y7GPeAzXhmG0QwV\nFfT15El14+jIsMgzDCMb5eX0lUVePVwS+cuXL3ceMWLE3ujo6EORkZGFCxYsWAIAZ8+e7ZGYmJgf\nHh5elJSUtK22ttZPmnAZhtETosj/73/qxtGRcUnkO3fufHnXrl2/PHToUPSPP/54665du375zTff\n3JmRkTE/MTExv6ioKDwhIWFHRkbGfKkCZhhGP5SXA/7+nMmrict2TdeuXS8BgNVq9WlsbPTs3r17\nTW5ubnJaWloWAKSlpWVlZ2enuDoPwzD6w2wG7riDM3k18XJ1gKamJo+hQ4ceOHHixIBHH310TVRU\n1GGLxWI0Go0WADAajRaLxWK0de2iRYt+/j4+Ph7x8fGuhsMwjIYoLwfi4oBXX1U7Ev1SUFCAgoIC\np6+XrITy3Llz3e65554vlixZsmDSpEn/rKmp6S4+1qNHj7Nnz57tcc3EXELJMG5PXByQng4kJVGt\nvLe32hHpH9VKKLt163Zu3Lhxm7///vvbjEajpaqqKgAAKisrA/39/U9JNQ/DMPqhvBzo3x8wGq8u\nwjLK4pLInzlzppdYOVNXV9clPz8/MSYm5mBycnJuVlZWGgBkZWWlpaSkZEsRLMMw+qGpCaiqAvr0\nIaFnX14dXPLkKysrA9PS0rKampo8mpqaPB5++OH3EhISdsTExBycMmXKR2vXrp0VHBxc+tFHH02R\nKmCGYfTB6dPATTcBnToB/fpxhY1acFsDhmFk4cABYOZM4NAhYMECwNcXeP55taPSP9zWgGEYTVBe\nDvTtS9/368d2jVqwyDMMIwvNRb5/f7Zr1IJFnmEYWeBMXhuwyDMMIwstRf7kSYCX4ZSHRZ5hGFlo\nLvI33QT4+ADV1erG1BFhkWcYRhaaizzAZZRqwSLPMIwstBR53hClDizyDMNIzqVLwOXLQI9mHas4\nk1cHFnmGYSSnvJzaGRiabdnhMkp1YJFnGEZyWlo1AJdRqgWLPMMwklNeDphM197Hmbw6sMgzDCM5\nnMlrBxZ5hmEkx5bIBwQAtbVAXZ06MXVUWOQZhpEcWyLv4UEWTlmZOjF1VFjkGYaRHFsiD3AZpRqw\nyDMMIzmtiTxviFIeFnmGYSSlqQmwWIDAwOsf40xeeVjkGYaRlFOnAD8/akjWEi6jVB4WeYZhJKU1\nqwbgMko1YJFnGEZS2hJ5rWbydXXue/4sizzDMJLSlsgHBQFmM/n2WuLIESA9HTh/Xu1IpIdFnmEY\nSTGbWxf5Ll2Abt1oYVZLFBfT18OH1Y1DDljkGYaRlLYyeUCbvrwo8v/5j7pxyIFLIl9WVhb0y1/+\ncldUVNThQYMG/WfVqlV/AICzZ8/2SExMzA8PDy9KSkraVltb6ydNuAzDaB17RF5rvnxxMRAVBfz7\n32pHIj0uiby3t3f9ihUrnjx8+HDUnj17Yl9//fXfHzlyZGBGRsb8xMTE/KKiovCEhIQdGRkZ86UK\nmGEY4rPPtOkhtyfyWtwQVVQETJzImfx1BAQEVEVHRx8CAF9f3wsDBw48Ul5e3jc3Nzc5LS0tCwDS\n0tKysrOzU6QIlmEY4swZ4Fe/Av71L7UjuR69ZvKTJrlnJu8l1UClpaXBBw8ejBkxYsRei8ViNBqN\nFgAwGo0Wi8VitHXNokWLfv4+Pj4e8fHxUoXDMG7N3/8OXLkC/Pe/akdyLRcuAPX1tBmqNfr3B3bt\nUi6m9qipAaxWIDoaaGykRWGjTcVSh4KCAhQUFDh9vSQif+HCBd/Jkyd/kpmZ+cSNN954zQdIg8Eg\nGAwGwdZ1zUWeYRj7qK8HVq8GUlOBkhK1o7kWMYtvfuxfS7S28FpcDISFUcyDBpFloyWRb5kAL168\n2KHrXa6uqa+v9548efInDz/88HspKSnZAGXvVVVVAQBQWVkZ6O/vf8rVeRiGIT75hETpvvu0K/Jt\nobUNUUVF9PsEgMGD3c+ycUnkBUEwzJo1a21kZGTh3LlzV4r3Jycn52ZlZaUBQFZWVpoo/gzDuM6q\nVcAf/gDcfLP27Bp7RL5nT7KatLJoXFwMhIfT9yzyLfjXv/418v33339o165dv4yJiTkYExNzMC8v\nb8z8+fMz8vPzE8PDw4t27tw5av78+RlSBcwwHZnvvgMqK4HkZCAkRJ+ZvMGgrcVX0a4Brto17oRL\nnvydd975TVNTk803iu3bt492ZWyGYa4nMxN47DHA05My4oYGWjjs3l3tyIjyciA0tP3nib58VJT8\nMbVHURF9MgJI5A8fprYLHm6yVdRN/hkM4/5UVgKbNwMzZ9LPBgNZNlrK5u3J5AHt+PKCcG0m360b\nvXlq6XfqKizyDKMT3niDKmqaZ+1as2zsFXmtVNicPg14eZGwi7ibZcMizzA64MoV4M03r9oKInoW\neS1k8s2zeBF3W3xlkWcYHbBxI23WiYi49n4tVdg0NtKpULaO/WuJVlobNC+fFBk8mDN5hmEURBBo\nwbVlFg9oK5O3WIAePQBv7/afq+VMftAgzuQZhlGQf/0LuHgRGDPm+se0JPL2WjUAYDIBVVVUHaQm\nzWvkRSIi6NPRlSvqxCQ1LPIMo3EyM4HHH7dd0hcSApSWauOkJUdE3tsb8Pena9TEVibfqRP9Xo8d\nUycmqWGRZxgNc/IksHMnkJZm+/GuXakZWGWlsnHZorycMnR7UbuMsmX5ZHPcafGVRZ5hNMzrrwPT\npwM33tj6c7Ri2TiSyQPql1FWVNDv9aabrn+MRZ5hGNm5dAlYt46smrbQSoWNoyKvdibfWhYPuFet\nPIs8w2iU998H4uJIxNuCM3nnsFU+KcKZPMMwsiIIV7tNtoeeRV6rmXxICFBdDZw7p2xMcsAizzAa\nZMcO6k0zalT7z9WzXaNmJm+rfFLEwwOIjKRmZXqHRZ5hNIiYxbd1wpKIFjL58+epjNPWImZriJm8\nYPPcOPlpK5MH3MeyYZFnGI1x4gTw7bfAgw/a93yTidoJqLl5x55j/1rSrRs1B6upkS+u1mhspE8/\nbbVFdpf2BizyDKMxXn0VmDWLauDtwcuLhF5Nf9tRq0ZErcXXsjLqPNnW79hd2huwyDOMhvjpJ+Dd\nd4Hf/96x60JC1PXlnRV5tcoo2/LjRUS7Ri07SSpY5BlGQ2RlAaNHA0FBjl2nti+vt0y+PT8eoLYL\nHh7UY0fPsMgzjEZoaiKr5oknHL9W7ROizGZ9ZfJt1ciLGAzusfjKIs8wGmHrVtpmHxfn+LV6tWvU\nzOTbs2sAFnlNUV8P7NqldhQM4zyZmZTFO1KhIqJXu0ZNT769TB5wj/YGbiPyr75KG0e++UbtSBjG\ncQoLKWN84AHnrlfbrtFTJl9fT3O21y4C4ExeM1RUAOnpwF//Cjz6KP0nMq7x/vv0O2WU4dVXgd/9\njnqZO0OvXlQnr8Y2/IYG4MwZICDA8WsDAqhO/vJl6eNqjdJSoE8f+37XUVHAkSNUV69XXBb5mTNn\nrjMajZbBgwf//H539uzZHomJifnh4eFFSUlJ22pra/1cnactnnmGXiALF9J/XmamnLO5P42NwAsv\nAJs2qR1Jx6Cmhs5wfeQR58cwGNSzbKqq6E3Gy8vxaz096ROA2Sx9XK1hrx8P0BqJv7822kY4i8si\nP2PGjHfy8vKuOZgsIyNjfmJiYn5RUVF4QkLCjoyMjPmuztMau3bR8WgLF9If+uuvAxkZtNmBcY5P\nP6WNIidOUN02Iy9vvw2MH+9cJtwctUTeWatGRGnLxl4/XkTvlo3LIn/XXXd93b1792s2Jufm5ian\npaVlAUBaWlpWdnZ2iqvz2KK+njaNrFwJ3HAD3RcaCjz2GDB3rhwzdgyWLweeew647TZgzx61o3Fv\nGhooMbGn22R7qNWozFWRV3rx1Z7yyebofeerEx+w2sdisRiNRqMFAIxGo8VisRhtPW/RokU/fx8f\nH4/4+HiH5snMpD+QlBZvIfPn07vvli3Avfc6FntHZ/duwGIBJk4Evv+efk5KUjsq9yU3lyzG2293\nfayQEMpSlUaPmfz48fY/f/Bg4J//lC+e9igoKEBBQYHT18si8s0xGAyCwWCwuTG4ucg7itlMtsy3\n315fcta5M2VHjzxC5U/29gBhgFdeoU9Bnp7AyJG8viE3YtmkFISEANu2STOWI0iRye/eLV087eGM\nXbN4sXzxtEfLBHixg8HIUl1jNBotVVVVAQBQWVkZ6O/vf0rqOZ55hippWvvPSkoChg0DliyRemb3\n5b//BQoKgJkz6ec77gD27iVLgZGeQ4fodz5pkjTj6dWuUTKTv3KFDj0PDrb/mvBwik/JCiApkUXk\nk5OTc7OystIAICsrKy0lJSVbyvF37CDxWbCg7eetWAGsWQMcOybl7O7LypXA7NmAry/93LMnvXj1\nvhlEq6xaBcyZA3h7SzNecDCVByrdUEtPnvyJE/Sm4kglkI8PrfUdOSJfXLIiCIJLt6lTp24IDAys\n8Pb2tppMprJ169bNqK6u7pGQkLA9LCysKDExcVtNTY1fy+toase5ckUQIiIEITvbvuevWCEIo0YJ\nQlOTU9N1GM6eFYTu3QXBbL72/lmzBOG119SJyZ05dUoQ/PwE4fRpacft3VsQKiqkHbM9wsMFobDQ\n+esvXhSETp0EobFRuphaIztbEMaNc/y61FRByMqSPh5n+H/ttFujXfbkN2zYkGrr/u3bt492dWxb\nrFxJH0uTk+17/mOPAevXUx1yqs1IGQB4801ajGqZkcXFAdu3O976lmmbN98EJk+m+nIpES2bwEBp\nx20NQXA9k+/alerRT51yvYy0PRz140X03N5AVztey8qAl1+mj7n29vfw8iLL5pln3ONQXjmwWmnH\n5dNPX//YyJHKLop1BKxW+puUasG1OUrXyv/0E70WHTn2zxZKWTaOlk+K6LlWXlci//TTlFEOGODY\ndXfcAYwbB/zxj/LEpXc2bgQGDgSGDLn+sfBw4MIFytYYafjkE+CWW0g4pEZpkXc1ixdRavHV2Uye\nRV4B8vOB/fupBt4ZliwBPv6Yar+ZqwgCbX566inbjxsMZNlwNi8dmZnSbH6yhdIVNlKJvFKZvCMt\nDZrTrx99alHjPFpX0YXIX7lC3npmJtCli3Nj9OxJQv/oo/puNiQ1O3bQzuExY1p/zsiR1DqCcZ29\ne8l7njBBnvE5k2+dixeB6mrHT90C6ISoqCh9+vK6EPkVK+jd19UXRloalUP9/e/SxOUOiFm8Rxt/\nCXFxLPJSkZlJCYunpzzjqyHyJpPr4yiRyR8/Tp902vpbbwu9tjeQfcerq5w8Cfztb8C+fa6P5eFB\nC16jRtEGFKPNZgsdh8OHgYMHqSFZWwwbRv3OL1682iOIcZyKCjr9afVq+eYICqKukFYrJTRyU15O\nGa6rKJHJO2vViAwezJm8LDz1FPD44/Y1+LeHwYOBX/8amDdPmvH0zCuv0Gaczp3bfl6XLsCttwLf\nfadMXO7KmjXAtGmAn4yNt729qReOUpuL9OTJO7voKqLXxVdNi/wXX1Cm+eyz0o774ou0fd+Fnj+6\np6qKmi49+qh9z2fLxjUuXyabUK4F1+YoadlIJfK9egF1dVTJJReuirxYK6/0jmJX0azIX7lCGfyq\nVc4vtraGry9tqpozhz7WdkRef52Omuvd277nc7288wgC8N57wNChVDopN0pW2Egl8gYDWTZyZvPO\n1siL9O5Np0nprZxYs5788uVUuz1unDzjT5wIrF1LloWzZZl65dIl4I03HDsPNy6O+to0NTm/cOXu\nNDRQBn306NXbkSP01WBof+1DKpTK5OvrqVpFqrUt0ZePjJRmvJa46skDVy0bKRablUKTIv+//5HI\n798v3xwGA+3yHD4cmDrVsa50eicrizaIOZJVBgQA3buTYMn1ItQLFy5Q07vmIn70KDW/CggAIiIo\nQRkxgiq6IiIoC7R3l7arhIQAOTnyz1NVRUfjSVUpJGcmf+4cFQ642u5BtGzGjpUmLiXQpMg/+ST1\nNA8JkXeem2+meZ54QpkXhRZoaqKS1LfecvxasV6+I4i8IFBLWltZeXU1ZYQREXT71a/oa1iYNs4u\nUMqukcqqEZFz8bW4mDpJuvpGO3iw/tbyNCfyW7cCP/4IfPihMvPNm0eVI7m59jc90zOff059Ru6+\n2/FrxcXX3/xG+ri0QmMjLUZv2kT+q5iVR0RQ9hYRQWKkZctKKbtGapHv1492tsuBq4uuIoMHkwOg\nJzQl8pcvU/XBq6+2X9YnFZ06Ud3yrFlAQoL714EvX049gJzJaEaOpDUMd6WpCfjd76gn+4kT0neI\nVAp/f6pU+ekn1xuHtYXeMnlX/XiAPsUeO0brL470pFcTTeUjf/sbeV5K+10JCZSl/vWvys6rNPv3\nU4Z3//3OXR8VRVvyT0l+zpf6CALtySgsBLKz9SvwAL2BK5HNm83SZ/JybYiSKpO/4Qbah3D8uOtj\nKYVmRL60lMoaV6xQZ/7ly4G336YXubuyfDmtPzh7EpGHBxAb656llC++CHz5JR3+Lp6MpWeUEHmp\nM3mTidZB5Dhu0tXyyeborb2BZkR+7lxacFWryiUwEHjhBaqd19tmB3s4eZIOeZ4927Vx3LFeftky\n6lC6bZu8u1GVRI8i7+NDVUgVFdKNKSKVXQPor72BJkR+82bqo/LMM+rGMWcO+Zjvv69uHHKQmUnt\nHLp1c20cd+tI+cYb1G5g+3b7N4bpASUqbKQWeUAeX766mtZbpLLg9NbeQHWRb77Y2qmTurF4etKL\n/tln9dk3ujXOnaMjEKU4iWj4cODQIdqRrHfefx946SUSeKnFSm3kzuSlOPbPFnL48qJVI9U+Bb0d\nBai6yL/8MhAd3XY/cyUZPhxISQGef17tSKTj7beBpCR6AbmKry+VEer98JVPP6Xy2S++kK75nZaQ\nW+Rra2ltR+r1Czk2REm16CoSFkaLzpcuSTemnKgq8iUl1JtGrcXW1khPJxGQor2x2tTXk1Vj6/xW\nZ9F7s7IvvqBSyc2b3Xdjlyjycq0vyZHFA/LYNVL68QC9uYWH66dIQ1WRnzuXytakyDClpHt3YOlS\n9zhF6h//oBf8sGHSjannxdevvwYeeojexIcOVTsa+bjxRir3s1jkGV8ukZfDrpE6kwf0ZdmoKvJH\njkibYUrJww/TC2XNGrUjcR7x/Fapf8fi4qveqpD27wcmTwY2bKB/g7sjp2Wjp0xeyvJJET0tvsom\n8nl5eWMiIiKOhoWFFS9duvQ5W8/RwmJraxgMwGuvAX/5C1keeuSrr4Dz54Hx46UdNyiI/t/0tCHk\n8GH6Pbz1FjB6tNrRKIOcFTZyZ/JSJRCCIE8m3+FFvrGx0fOxxx57LS8vb0xhYWHkhg0bUo8cOTKw\n5fPuuUeO2aVj0CDq1Lh5s9qROMfy5bT3QI4+K3qybI4fp4XnV14B7rtP7WiUQ4+ZvJ8f/b3W1koz\nnsVCLVK6d5dmPJEOb9fs27dveGho6PHg4OBSb2/v+qlTp27MycnR5ctr5kxg3Tq1o3CcY8eAPXuA\n6dPlGV8v9fJlZUBiIu1onTZN7WiURY8iD0jry8uRxQP0afbiRarB1zqytNgpLy/vGxQUVCb+bDKZ\nzHv37h3R8nmLFi36+fv4+HjEx8fLEY5L3H8/ZcOVla73olaSFSuARx6Rr/VtXBzw5pvyjC0VFgtZ\nM48/Dvz2t2pHozwhIcDGjfKMLafIi758dLTrY8nhxwNk54rtDeSWrYKCAhS40N9YFpE3GAx2OWrN\nRV6r+PrSYt1770l/1qxcnD5NrXKPHpVvjiFDKNuqqZH+o7AU1NSQRZOaShVcHRE9evKA9Jm8lOWT\nzRHbG8gt8i0T4MWLFzt0vSx2Td++fcvLysqCxJ/LysqCTCaTWY65lGDGDLJs9FJN8sYbwKRJ0h3L\nZgsvL+D224Fvv5VvDmc5f546mY4eTTZNR6VfP/oEKnXhgNVKb6L+/tKOKyLlhii57BpAP4uvsoj8\nsGHD9hcXF4eVlpYGW61Wn02bNj2QnJycK8dcShAXRwKvRUFryeXLdEi3EtmrFhdf6+ro8Jdbb6XW\n1UoduadFvL3pOMKysvaf6wiVlZRASHXsX0ukLKOUy64B9NONUhaR9/Lyanjttdceu+eee76IjIws\nfOCBBzYNHDjwiBxzKYHBoJ8F2A8+AGJiqPe73Ght8dVqpaP4+vSh/Q0dWeBF5LBs5LRqAOnsmqYm\nOvwlNNT1sWwh2jVa/4Qv29kmY8eO3Tp27Nitco2vNNOn0xb4zEztnh4lCFQmmJmpzHyxsbTBqL7e\n+R71UtHYSBvYPD2pGZtcWabekKPCRm6RlyqTLy+nkswbb3R9LFv06EFjnzxJMWsV1RuU6YXAQODO\nO6lNgFbJyyOvPCFBmfn8/Kj//w8/KDNfazQ10bmz1dW04Kz2G46W0KPIBwbS/6WrnU7l9ONF9FAv\nzyLvAFq3bFw5v9VZ1LZsBIFKXI8dA3JylDsbWC/IZdeYTNKO2RxPT7LczC6Wasjpx4voYfGVRd4B\nxo2jssTiYrUjuZ5Dh6gX0NSpys4bF6fu4uvKldR0bPNm7dpoaqLHTB6QxpdXIpNnkXczfHyog+H6\n9WpHcj2vvEKbfnx8lJ1XzWZlly9Tt9APPnCfY/ukRq8iL4UvL2eNvAjbNW7IjBlAVpa2WhCXlwOf\nf0490pXm5ptp4VXqzoH28OGH1C544HVdkRiRgADgwgW6SQVn8leJjCRbSMtNDFnkHWTQIPILt21T\nO5KrvPwykJamzs5Tg0GdenmxkujJJ5WdV28YDLQ4LlU2Lwh00LbWM/mGBvo3DxggXUy26NKF+tgU\nFck7jyuwyDuBlhZgzWZqufCczWbOyqDG4uv27SRgHaVtsCtIadnU1FCbabl6Iom4uuv15EnasNWl\ni3QxtYZYL69VWOSdYOpUID8fOHNG7UjoIOrZs+ljuVqocRygmMXzhqf2kbLCRgmrBnDdrlHCqhHR\n+uIri7wT+PnRARQffqhuHKWlwEcfqd84behQelGdP6/MfIWFwMGDHa91sLNImckrKfJlZc4v6CtR\nPimi9fYGLPJOMnMmsHatulua//IXYM4coFcv9WIA6ON7TAywd68y861cSefvck28fehR5G+4gW6n\nTjl3vdKZPNs1bkh8PPDTT5RRqkFxMW3+0UobXaUsm9OngY8/JpFn7ENKu8ZsVkbkAdcWX5UonxQJ\nDaWmbVJWMEkJi7yTeHhcbUGsBosXA088oZ1e7kpV2LzxBvX3l6vNrTsSEkLWnhSfOpXK5AHXfHkl\nM3lPTyAigmxELcIi7wJpacCGDbQpR0kKC6mE84knlJ23Le64g44blHP/wJUrwOrVXDbpKDfdRJba\n6dOuj6WkyDubyVut5OeHhEgfU2toefGVRd4F+venRcfsbGXnXbSIetTcdJOy87ZF795U4XP4sHxz\nbNhAfeKVaKPsbkhl2eghky8pod46Su7+ZpF3Y5Sumf/hB+Crr4DHHlNuTnuRs15eEOjcWq2sQegN\nqRZf9ZDJK+nHi2i5vQGLvIukpADff6/ctv4XX6SNT1psxiXn4uvOnbSLMSlJnvHdHSlE/soVKjbo\n3VuamNrD2Q1RSvrxIpzJuzFdutDmqKws+efav59ujzwi/1zOIOfi64oVvPnJFaSwayoqyJLzUEg1\nnLVrlKyRF+nTh9YCnC35lBMWeQmYORN45x06vEJOXngBWLhQma3aznDLLcC5c1ROJiVHjwLffQc8\n+KC043YkpMjklbRqAKqguniRbo6ghl1jMGi3Xp5FXgKGDqVjwL78Ur45du+mRc1Zs+Sbw1U8PKjK\nRmrLJjOTPr1o9c1ND+hR5A0Gav7lqGWjhl0DaNeyYZGXACUO+v7Tn+jWqZN8c0iB1JZNdTWwcSPt\n7GWcp39/EumGBufHUFrkAccXX+vqAIuFrB6l0Wp7AxZ5iXjwQeCzz8iukJqCAvIm09KkH1tqpF58\nffNNYOKylXGQAAAVJElEQVRE6ijIOI+PD/0Oy8qcH0MNkXfUlz9xglore3nJFlKrsF3j5vTqRW1v\nN26UdlxBoAz+hRf0cUD17bfTH3pdnetjWa3Aa6/x5iepcNWy0UMmr4YfLzJoEFmqcq/NOQqLvITI\nYdmILY31sujYtSv9sX/3netjbdpEG58GD3Z9LMb1Chs9ZPJq+fEAdaf183P9RCupcVrkP/74419F\nRUUd9vT0bDxw4MDQ5o8tWbJkQVhYWHFERMTRbdu2dZjK5qQkauAk1a5PMYtftIj6Y+gFKSwbPvlJ\nejpCJq9G+WRztLj46rTIDx48+N+ffvrpxLvvvvur5vcXFhZGbtq06YHCwsLIvLy8MXPmzFnd1NTU\nIT4xeHmRb/7OO9KMt3kz2R6/+pU04ymFFIuvX35JPYHGjJEmJsY1kVfq2L+W6CmTB9xM5CMiIo6G\nh4dfd7JhTk7OfampqRu8vb3rg4ODS0NDQ4/v27dvuGth6ocZM+g4PlcP9m1qoix+8WLlNp9IRVwc\nibwr3uQrrwBz5+rv365lXLFrqqvJilO6jNVkojcXexvfqenJA9psbyD5GnRFRUWf2NjYPeLPJpPJ\nXF5ebvP9f9GiRT9/Hx8fj/j4eKnDUZywMNoUtHkztTxwlk8/JYFzZQy16NOHmqcdOwYMHOj49UVF\n1NFS6kXsjo4rmbwaVg1AJcO9etEGO5Op7edeuADU1qoTp8jgwcDSpdKOWVBQgIKCAqevb1PkExMT\n86uqqq47PTQ9PX3hhAkTPrN3EoPBYLOTdXORdyfEBVhnBbqxkXrULF2q3238omXjjMhnZgK//a38\nh0V3NAICqPfMxYuO9z5SS+SBq5ZNeyJfXAwMGKDup7+ICCrjtFql64LZMgFevHixQ9e3KfL5+fmJ\njgbUt2/f8rKysiDxZ7PZbOrbt2+5o+PomfvvpwXDykogMNDx6z/6iHbQ3nuv9LEphbj46ugO3bNn\nqaWwnC2LOyoeHrSQWVJCtoIjqCny4uLryJFtP09tqwagIymDg+lTrFaqwiR5zxME4ed8Mzk5OXfj\nxo1TrVarT0lJSUhxcXHY8OHD90kxj17w9QUmTSJv3lEaGqia5i9/0W8WDzi/+Pr3vwMTJjj35si0\nz803O2fZaCGTbw+1F11FtLb46rTIf/rppxODgoLK9uzZEztu3LjNY8eO3QoAkZGRhVOmTPkoMjKy\ncOzYsVtXr149pzW7xp0RLRtHj1z74AP6WJ2QIE9cSjFoEH2SOXPG/mt485P8OOvLl5e3b5fIhb1l\nlGqXT4porb2B0yI/ceLET8vKyoLq6uq6VFVVBWzdunWs+NjChQvTjx8/Hnr06NGIe+655wtpQtUX\ncXEk8N9+a/819fVUTaP3LB6guv4RIxzL5j/+mD5uR0fLF1dHJyTEuQobzuTtR2vtDbhATSbEpmWO\n1My/8w4tHN19t3xxKYkjlo24+YlPfpIXPdo19mbyWvDkATeya5j2mT4d+Mc/7OuHfeUK8Ne/Uhbv\nLjiy8/Xrr6kETs+LzXrAFbtG7Uy+LeuzpoY2z2mhkV1ICB2a/tNPakdCsMjLSGAgcOedJPTt8dZb\nlAHExsofl1LExgIHD9IbWHusWMGbn5RAFHlH1oouX6Y34F695IurLfz86GtbHV5Fq0YLNqenJxAZ\nqZ0KMX5JyYw9Tcvq6oAlS4A//1mZmJTixhvphXfwYNvPO34c+OYb+uTDyIufH3UzdWRBvKKCEha1\nBNRgaP+8V6348SJasmxY5GVm3DjgyBH6I2yNNWtokfK225SLSynssWxWrQJ+8xttHk7ujjhq2ahp\n1Yi0t/iqFT9eREvtDVjkZcbHB3joIWD9etuPX7gAvPwyVdW4I+0tvtbWAu+/D/z+98rF1NHRo8i3\nt/jKmXzrsMgrwMyZQFaW7SZLr74KxMdrZ3ec1IwcSZl8ax7wW2/Rpx21RaQj4WijMrNZ/f+f9jJ5\nrdTIi4i18o7uk5EDFnkFGDSImnZt23bt/efOUdmgm7bwAUAvTk9P26JSX09WDW9+UhZ3y+QFQXt2\nTcD/d/yyWNSNA2CRVwxbNfMrVwJjx1JTI3fFYGjdsvnkE9oXMHTo9Y8x8qFHkW8rkz9zhqqyevZU\nNqa2MBi0Y9mwyCvE1KmUyYtVDWfPklXz4ovqxqUEthZf+eQn9XDUrtGCyLeVyWvNqhHRSnsDFnmF\n8PMDxo8HPvyQfl6+HJg4kTJZd0f05ZuzezdtYBk/Xp2YOjL9+5PPbu9BHFoQ+cBASpCs1usf09qi\nq4hW2huwyCvIzJnA2rW0G+6NN4A//lHtiJQhOprsgdraq/e98grwxBP6OrvWXejUCejdm4S+PZqa\nqNFcnz7yx9UWXl4k9LZi1pofL8J2TQckPp62Oj/0ENk3/furHZEyeHsDw4bRaU8AWQVffgn8+teq\nhtWhsdeyOXOGNrV17ix/TO3R2oYorWbyQ4YA6elqR8EirygeHiRsX34JLFyodjTK0nzxddUqOkzE\n11fdmDoy9i6+asGqEenf3/biq1Y9+a5dgUSHj12SHsnPeGXa5rHHKKvVygtHKUaOpHWIc+eAd98F\nfvxR7Yg6NnoUeVuZvCBQWwwtirxW4ExeYXr2pM0/HY3YWGDfPlqLGDNGvQMoGMJeu0ZrIt8yk6+s\npHYY3bqpE5MeYJFnFKFHD3qR/vnP3DNeC+gxk7dVRqlVP15LsF3DKEZcHIn9sGFqR8I4IvJaaX9t\nK5PXqh+vJVjkGcVYuJAOKmfUp08f2qdw6RItELaGljJ50ZMXhKttj7VaPqkl2K5hFCMkhLMureDh\nQfZHaWnbz9OSyPv60hvS6dNX72O7pn1Y5Bmmg2KPZaMlkQeu9+XZrmkfFnmG6aC0V2FTV0d2jpYa\nfzUvo2xqovhDQ9WNSeuwyDNMB6W9TL68nLx7LZybKtJ8Q1RZGb0B8YlibeO0yM+bN2/ZwIEDjwwZ\nMuSHSZMm/fPcuXM/V6ouWbJkQVhYWHFERMTRbdu2JUkTKsMwUmKPyGvJqgGuzeTZj7cPp0U+KSlp\n2+HDh6N++OGHIeHh4UVLlixZAACFhYWRmzZteqCwsDAyLy9vzJw5c1Y3NTXxJwaG0Rjt2TVaFXkx\nk2c/3j6cFt/ExMR8Dw+PJgAYMWLEXrPZbAKAnJyc+1JTUzd4e3vXBwcHl4aGhh7ft2/fcKkCZhhG\nGsRMvrUj6rQo8s0XXjmTtw9J6uTXrVs3MzU1dQMAVFRU9ImNjd0jPmYymczl5eU2/1QWNTv3Lj4+\nHvHx8VKEwzCMHXTvTqWUZ8/aXlwtL6fMWUs0z+SLi6mzq7tTUFCAgoICp69vU+QTExPzq6qqAlre\nn56evnDChAmfAcBLL730vI+Pj3XatGkftjaOwWCwmSsscufDTRlGB4iWTWsiHxenfExt4e8PnD9P\nVT8dJZNvmQAvXrzYoevbFPn8/Pw2G2WuX7/+11u2bLl3x44dCeJ9ffv2LS8rKwsSfzabzaa+ffuW\nOxQVwzCKIFo2t99+/WNatGs8PICgIHpj+t//6E2KaRunPfm8vLwxy5Ytm5eTk3Nf586dL4v3Jycn\n527cuHGq1Wr1KSkpCSkuLg4bPnz4PmnCZRhGStqqsNGiyAPky3/1FZ0UpYXDTLSO0578448//qrV\navVJTEzMB4A77rjj29WrV8+JjIwsnDJlykeRkZGFXl5eDatXr57Tml3DMIy6hITY7u3f1ARUVal/\n7J8t+vUDduzoGFaNFBiE1pbW5Z7YYBDUmpthGGLrVmDFCmDbtmvvt1jojNJTp9SJqy0WLwZWrgSm\nTQNef13taJTHYDBAEAS7t6hx/TrDdGBas2u0atUAlMnX1nImby8s8gzTgQkOprrzxsZr79eyyPfv\nT1+5xbB9sMgzTAemc2egVy8S9eaYzdoVebF2nzN5+2CRZ5gOji3LRsuZfFAQZfPBwWpHog9Y5Bmm\ngxMScn0PGy2LfKdOdNiJt7fakegDFnmG6eDcfLO+MnnGMVjkGaaDoze7hnEMFnmG6eDoza5hHINF\nnmE6OC3tmosXAauVulQy+odFnmE6OH36ULvhujr6WYvH/jHOwyLPMB0cT08qSywtpZ/ZqnEvWOQZ\nhrnGsmGRdy9Y5BmGuabChkXevWCRZxiGRd6NYZFnGObnYwABFnl3g0WeYRjO5N0YFnmGYX7eECUI\nLPLuBos8wzDo0YO+VlfTqVCBgerGw0gHizzDMDAYKJvfu5d2uvr4qB0RIxUs8gzDACCR//prtmrc\nDRZ5hmEAUIXNN9+wyLsbLPIMwwCgTP6771jk3Q2nRf5Pf/rTX4YMGfJDdHT0oYSEhB1lZWVB4mNL\nlixZEBYWVhwREXF027ZtSdKEqg4FBQVqh2AXHKe0dMQ4Q0Ko+6QcIq+H36ceYnQGp0X+2WefffmH\nH34YcujQoeiUlJTsxYsXvwgAhYWFkZs2bXqgsLAwMi8vb8ycOXNWNzU16fYTg17+4zlOaemIcd58\nM31lkXcvnBbfG2+88bz4/YULF3x79ep1BgBycnLuS01N3eDt7V0fHBxcGhoaenzfvn3DpQiWYRj5\nEA/GNplUDYORGC9XLn7++edfeu+99x7u0qVLnSjkFRUVfWJjY/eIzzGZTOby8nJ2+RhG43TpQll8\nUFD7z2X0g0EQhFYfTExMzK+qqgpoeX96evrCCRMmfCb+nJGRMf/YsWO3vPPOOzMef/zxV2NjY/c8\n+OCDHwDA7Nmz37733nu3TJo06Z/XTGwwtD4xwzAM0yqCINh9pEubmXx+fn6iPYNMmzbtw3vvvXcL\nAPTt27e8+SKs2Ww29e3bt9yVIBmGYRjncNqTLy4uDhO/z8nJuS8mJuYgACQnJ+du3LhxqtVq9Skp\nKQkpLi4OGz58+D4pgmUYhmEcw2lPfsGCBUuOHTt2i6enZ+OAAQNOrFmz5lEAiIyMLJwyZcpHkZGR\nhV5eXg2rV6+ew9YMwzCMSgiCoPht69atY2655ZajoaGhxRkZGc+pEUN7t5MnTwbFx8fvioyMPBwV\nFfWfzMzMP6gdU1u3hoYGz+jo6IPjx4//TO1YWrvV1NT4TZ48+R8RERFHBg4cWPjtt9/Gqh1Ty1t6\nevqCyMjIw4MGDfp3amrqh5cvX+6kdkyCIGDGjBnr/P39LYMGDfq3eF91dXWP0aNH54eFhRUlJiZu\nq6mp8dNinM8888yyiIiII7feeusPEydO/GdtbW03LcYp3v72t789bTAYmqqrq3toNc5Vq1Y9HhER\ncSQqKuo/zz777NK2xlA86IaGBs8BAwYcLykpCbZard5Dhgw5VFhYOFDtX2bLW2VlZcDBgwejBUHA\n+fPnfcPDw49pMU7xtnz58qemTZv2wYQJE3LVjqW12/Tp07PWrl07UxAE1NfXe2nhxd78VlJSEhwS\nEvJfUdinTJmyaf369WlqxyUIAr766qu7Dhw4ENP8xT5v3ryXly5d+qwgCMjIyHjuueeey9BinNu2\nbUtsbGz0EAQBzz33XIZW4xQESu7uueeevODg4BItiLytOHfu3PnL0aNH51utVm9BEHDq1KnebY2h\n+Calffv2DQ8NDT0eHBxc6u3tXT916tSNOTk59ykdR3sEBARURUdHHwIAX1/fCwMHDjxSUVHRR+24\nbGE2m01btmy5d/bs2W8LGl3QPnfuXLevv/76rpkzZ64DAC8vr4Zu3bqdUzuu5tx0000/eXt711+6\ndKlrQ0OD16VLl7raKhpQg7vuuuvr7t271zS/Lzc3NzktLS0LANLS0rKys7NT1InuKrbiTExMzPfw\n8GgCgBEjRuw1m82qV+LbihMAnnrqqVdefvnlZ9WIyRa24lyzZs2jCxYsWOLt7V0PAL179z7d1hiK\ni3x5eXnfoKCgMvFnPdTRl5aWBh88eDBmxIgRe9WOxRZPPvnkimXLls0TX0hapKSkJKR3796nZ8yY\n8c7QoUMP/OY3v3nr0qVLXdWOqzk9evQ4+/TTTy/v16/fyT59+lT4+fnVjh49ervacbWGxWIxGo1G\nCwAYjUaLxWIxqh1Te6xbt26mWImnNXJycu4zmUzmW2+99Ue1Y2mL4uLisK+++uru2NjYPfHx8QX7\n9+8f1tbzFRd5vS3CXrhwwff+++//R2Zm5hO+vr4X1I6nJZ9//vl4f3//UzExMQe1msUDQENDg9eB\nAweGzpkzZ/WBAweG3nDDDRczMjLmqx1Xc06cODFg5cqVc0tLS4MrKir6XLhwwfeDDz54UO247MFg\nMAhaf2299NJLz/v4+FinTZv2odqxtOTSpUtd09PTF4rtWQDtlnk3NDR41dTUdN+zZ0/ssmXL5k2Z\nMuWjtp6vuMi3rKMvKysLMplMZqXjsIf6+nrvyZMnf/LQQw+9n5KSkq12PLbYvXt3XG5ubnJISEhJ\namrqhp07d46aPn36u2rH1RKTyWQ2mUzm22+//TsAuP/++/9x4MCBoWrH1Zz9+/cPi4uL292zZ89q\nLy+vhkmTJv1z9+7dcWrH1RpGo9EiblasrKwM9Pf3P6V2TK2xfv36X2/ZsuVerb5pnjhxYkBpaWnw\nkCFDfggJCSkxm82m22677ftTp075qx1bS0wmk1ncXHr77bd/5+Hh0VRdXd2ztecrLvLDhg3bX1xc\nHFZaWhpstVp9Nm3a9EBycnKu0nG0hyAIhlmzZq2NjIwsnDt37kq142mN9PT0hWVlZUElJSUhGzdu\nnDpq1Kid77777nS142pJQEBAVVBQUFlRUVE4AGzfvn10VFTUYbXjak5ERMTRPXv2xNbV1XURBMGw\nffv20ZGRkYVqx9UaycnJuVlZWWkAkJWVlabVRCQvL2/MsmXL5uXk5NzXuXPny2rHY4vBgwf/22Kx\nGEtKSkJKSkpCTCaT+cCBA0O1+MaZkpKSvXPnzlEAUFRUFG61Wn169uxZ3eoFaqwYb9myZWx4ePix\nAQMGHE9PT1+g9gq2rdvXX399p8FgaBoyZMih6Ojog9HR0Qe3bt06Ru242roVFBT8QsvVNYcOHRoy\nbNiw77RUStfytnTp0mfFEsrp06dniRUMat+mTp26ITAwsMLb29tqMpnK1q1bN6O6urpHQkLCdi2V\nULaMc+3atTNDQ0OL+/Xr9z/xdfToo4+u1kqcPj4+V8TfZ/PHQ0JC/quF6hpbcVqtVu+HHnrovUGD\nBv176NCh3+/atSu+rTHa7F3DMAzD6Bvd9nlnGIZh2odFnmEYxo1hkWcYhnFjWOQZhmHcGBZ5hmEY\nN4ZFnmEYxo35P6FUAlB8DZAtAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x6c5db10>"
       ]
      }
     ],
     "prompt_number": 70
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Cross-validation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "1. pick a cross validation scheme (here, leave one out)\n",
      "2. make test and train datasets\n",
      "3. fit the model on the train data\n",
      "4. measure how well is this model doing on the left out data (square prediction error)\n",
      "5. repeat across (test,train) and cumulate results"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def make_test_train(y, lmodels, scheme='leave_one_out'):\n",
      "    \"\"\"\n",
      "    Parameters\n",
      "    ----------\n",
      "    y : numpy array\n",
      "        the data\n",
      "    lmodels : list of numpy arrays\n",
      "        each element is a design matrix\n",
      "    Returns\n",
      "    -------\n",
      "    \n",
      "    \"\"\"    \n",
      "    if scheme=='leave_one_out':\n",
      "        for i, y_test in enumerate(y):\n",
      "            y_train = np.hstack((y[:i],y[i+1:]))\n",
      "            m_loo = []\n",
      "            for mX in lmodels:\n",
      "                m_loo.append(np.vstack((mX[:i,:],mX[i+1:,:])))\n",
      "            yield i, y_test, y_train, m_loo"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 71
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = len(y)\n",
      "x = np.arange(n)\n",
      "X1 = np.hstack((np.ones((n,1)),x.reshape((n,1))))\n",
      "X2 = np.hstack((X1, (x*x).reshape((n,1))))\n",
      "models_list = [X1, X2]\n",
      "\n",
      "for k in range(2,8):\n",
      "    models_list.append(np.hstack((X2, np.random.normal(size=(n,k)))))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 72
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y_split = make_test_train(y, models_list)\n",
      "\n",
      "predicted_error = np.zeros(len(models_list),dtype='float')\n",
      "\n",
      "for i, y_test, y_train, models_train in y_split:\n",
      "    for midx, md in enumerate(models_train):\n",
      "        # fit the model on train\n",
      "        bh = lin.pinv(md).dot(y_train)\n",
      "        # compute the difference between y_test and predicted\n",
      "        predicted_error[midx] += (y_test - models_list[midx][i,:].dot(bh))**2\n",
      "\n",
      "print(predicted_error/n)\n",
      "plt.plot(predicted_error/n)  "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 1216.81644631   419.52053563   520.91121718   551.83232086   851.68596693\n",
        "   800.49557814   843.35431674  1096.09708688]\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 73,
       "text": [
        "[<matplotlib.lines.Line2D at 0x6c531d0>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD9CAYAAABdoNd6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHX6B/DPCGPm/cqAM9QYDOIoIV6QLhiloNlKWIZi\nJYmubeyW3dN1t2R3FaxtV22j9ldoZCXSDbCURWunrBRUNM0JGQ2UGWA0EcUrCOf3x7cpcgFhGDgz\nzOf9ep0XeGbOOQ+lD2ee8/0+X4UkSSAioq6tm9wBEBFRx2OyJyJyA0z2RERugMmeiMgNMNkTEbkB\nJnsiIjfQYrJPSEhYq1KprEFBQQds+/785z//NTg4+NvRo0fvmzRp0mdlZWW+tteSk5OX6HQ6U2Bg\nYFFeXl6Ubf+ePXvGBgUFHdDpdKZFixat7pgfhYiImiVJUrPbl19+GV5YWBgyatSoA7Z9Z86c6WP7\nfs2aNY/Onz//TUmScPDgQX1wcPC+2tpaZUlJidbPz+9wQ0ODQpIkjB8/viA/Pz9UkiTceeedm7ds\n2TK1pety48aNGzfHbi3e2YeHh28fMGDAqcb7+vTpU2P7/uzZs70HDx78IwBkZ2ffHRcXt0GpVNZp\ntdpSf3//w/n5+RMqKip8ampq+oSGhhYAwNy5c9/OysqK6YhfXERE1DRPew5aunTp8vXr1z947bXX\nXigoKAgFgPLy8qFhYWE7be/RaDRmi8WiViqVdRqNxmzbr1arLRaLRd3+0ImIqLXsSvbLly9funz5\n8qUpKSmLH3/88VXr1q2b54hgFAoFezcQEdlBkiRFS6+3azTOnDlz3tu1a9d4QNyxN35YazabNRqN\nxqxWqy1ms1nTeL9arba0ELDLbi+88ILsMbhj7Ixf/o3xy7u1RpuTvclk0tm+z87OvjskJGQvAERH\nR+dkZGTMrq2t7V5SUjLMZDLpQkNDC7y9vSv79u17Jj8/f4IkSYr169c/GBMTk9XW6xIRkf1aLOPE\nxcVt+OKLL2778ccfB/v6+pYlJSW9sHnz5mmHDh0a7uHhUe/n53fktddeewQA9Hq9MTY2NlOv1xs9\nPT0vp6amJtrKMqmpqYkPPfTQWxcuXLh22rRpm6dOnZrbGT8cEREJitZ+BOgMCoVCcqZ42spgMCAi\nIkLuMOziyrEDjF9ujF9eCoUC0lVq9kz2REQurjXJnu0SiIjcAJM9EZEbYLInInIDTPZERG6AyZ6I\nyA0w2RMRuQEmeyIiN8BkT0TkBpjsiYjcAJM9EZEbYLInInIDTPZERG6AyZ6IyA0w2RMRuQEmeyIi\nN9Bisk9ISFirUqmsQUFBB2z7nnnmmZdGjBjxfXBw8Lf33HPPR6dPn+5ney05OXmJTqczBQYGFuXl\n5UXZ9u/Zs2dsUFDQAZ1OZ1q0aNHqjvlRiIioOS0m+3nz5q3Lzc2d2nhfVFRU3sGDB0d+++23wQEB\nAcXJyclLAMBoNOo3btw4y2g06nNzc6cmJiam2prpP/LII6+lpaXNN5lMOpPJpLvynI1x7RIiIsdr\nMdmHh4dvHzBgwKnG+yIjI7d269atAQAmTJiQbzabNYBYfDwuLm6DUqms02q1pf7+/ofz8/MnVFRU\n+NTU1PQJDQ0tAIC5c+e+nZWVFdPcNUtL2/0zERHRFVpccPxq1q5dmxAXF7cBAMrLy4eGhYXttL2m\n0WjMFotFrVQq6zQajdm2X61WWywWi7q5cz755DIEB4vvIyIiXHpdSCKijmAwGGAwGNp0jN3Jfvny\n5Uu7d+9eO2fOnPfsPUdTBg5chmXLHHlGIqKu5cob4aSkpKseY1eyf+uttx7avHnztM8++2ySbZ9a\nrbaUlZX52v5sNps1Go3GrFarLbZSj22/Wq22NHfuL7+0JyIiImpJm4de5ubmTn3ppZeeyc7OvrtH\njx4Xbfujo6NzMjIyZtfW1nYvKSkZZjKZdKGhoQXe3t6Vffv2PZOfnz9BkiTF+vXrH4yJiclq7vwn\nTwLl5fb+OERE1JQWk31cXNyGm2+++ZtDhw4N9/X1LVu7dm3Co48++srZs2d7R0ZGbg0JCdmbmJiY\nCgB6vd4YGxubqdfrjXfeeeeW1NTURIVCIQFAampq4oIFC97U6XQmf3//w1OnTs1t7prh4cD27Y79\nIYmI3J1CcqKxjgqFQvr73yX88APw6qtyR0NE5BoUCgVsQ92b43QzaCdOZN2eiMjRnO7Ovq5OwsCB\nQEkJMGiQ3BERETk/l7yz9/QEbroJ+OoruSMhIuo6nC7ZAyzlEBE5GpM9EZEbcLqavSRJuHgRGDwY\nqKgA+vSROyoiIufmkjV7AOjRAxg7FtixQ+5IiIi6BqdM9gBLOUREjuS0yT48nMmeiMhRnLJmDwA1\nNYCPD/Djj6KsQ0RETXPZmj0gHszq9UBBgdyREBG5PqdN9gDr9kREjsJkT0TkBpy2Zg8AVVWAVit6\n3CuV8sVFROTMXLpmDwADB4pkv3ev3JEQEbk2p072AEs5RESO0GKyT0hIWKtSqaxBQUEHbPvef//9\n+0aOHHnQw8OjvrCwcEzj9ycnJy/R6XSmwMDAory8vCjb/j179owNCgo6oNPpTIsWLVrdlgCZ7ImI\n2q/FZD9v3rx1ubm5UxvvCwoKOvDxxx/PmDhx4q9SsNFo1G/cuHGW0WjU5+bmTk1MTEy11ZAeeeSR\n19LS0uabTCadyWTSXXnOloSHi3bHDQ1t+bGIiKixFpN9eHj49gEDBpxqvC8wMLAoICCg+Mr3Zmdn\n3x0XF7dBqVTWabXaUn9//8P5+fkTKioqfGpqavqEhoYWAMDcuXPfzsrKimltgD4+oinawYOtPYKI\niK7k6agTlZeXDw0LC9tp+7NGozFbLBa1Uqms02g0Ztt+tVptsVgs6ubOs2zZsp+/j4iIQERExM+l\nnKAgR0VLROS6DAYDDAYDAODIkdYd47Bk7yiNk73NxInAp58Cv/9958dDRORsbDfCW7cCb7wBAElX\nPcZho3HUarWlrKzM1/Zns9ms0Wg0ZrVabTGbzZrG+9VqtaUt57Y1RXOiKQFERLIym4G5c4F3323d\n+9uV7BsP4o+Ojs7JyMiYXVtb272kpGSYyWTShYaGFnh7e1f27dv3TH5+/gRJkhTr169/MCYmJqst\n19Fqxdq0hw+3J1oioq6hrg6YNQt47DEgIqKVB0mS1Ow2e/bsDT4+PuVKpbJWo9GUpaWlJXz88ccx\nGo2mrEePHhdUKlXl1KlTt9jev3z58j/6+fkdHj58eFFubu4U2/7du3ePHTVq1AE/P7/Djz766Jrm\nrifCadqcOZL05pvNvkxE5DaeekqSpk2TpPp68eefcmeL+dyp2yU09u9/A998A6Snd3JQRERO5OOP\ngSeeAPbsAQYNEvtcvl1CY5xcRUTu7sgR4OGHgczMXxJ9a7lMsg8MBM6eBY4dkzsSIqLOd+ECMHMm\n8PzzQGho2493mWSvUIi7++3b5Y6EiKjzLVoEBATYPwTdZZI9wFIOEbmnt98Wue/NN8WNrz2Y7ImI\nnNh33wFPPQV88IFYrtVeLpXsb7wRqKgAjh+XOxIioo5XUyPq9H//OzBqVPvO5VLJ3sMDuOUW0QWT\niKgrkyRg4ULRQSA+vv3nc6lkD7CUQ0TuITUVKCoC1qxxzPmY7ImInMyuXUBSEvD++8C11zrmnC4z\ng9amtlasTWuxAP36dVJgRESdpKoKGDsWePll4J57WndMl5pBa9O9u5hQ8PXXckdCRORYDQ2ik+WM\nGa1P9K3lcskeYCmHiLqmF18Ud/YrVzr+3Ez2REROwGAAVq0CNm4ElErHn9/lavYAcP48MGQIcOIE\n0LNnJwRGRNSBKitFnX7dOiAqqu3Hd8maPSASfHAwsHPn1d9LROTMLl8G4uKABQvsS/St5ZLJHmAp\nh4i6huefFyvxPf98x16nxWSfkJCwVqVSWYOCgg7Y9lVVVQ2MjIzcGhAQUBwVFZVXXV3d3/ZacnLy\nEp1OZwoMDCzKy8v7+XfUnj17xgYFBR3Q6XSmRYsWrXZE4Ez2ROTqPv0UWL9erCPr4dGx12ox2c+b\nN29dbm7u1Mb7UlJSFkdGRm4tLi4OmDRp0mcpKSmLAcBoNOo3btw4y2g06nNzc6cmJiam2mpIjzzy\nyGtpaWnzTSaTzmQy6a48pz1uuQUoKBDj7omIXM3Ro0BCArBhA+Dl1fHXazHZh4eHbx8wYMCpxvty\ncnKi4+Pj0wEgPj4+PSsrKwYAsrOz746Li9ugVCrrtFptqb+//+H8/PwJFRUVPjU1NX1CQ0MLAGDu\n3Llv245pj379RG/nPXvaeyYios516RJw333As88Ct97aOdf0bOsBVqtVpVKprACgUqmsVqtVBQDl\n5eVDw8LCfn5kqtFozBaLRa1UKus0Go3Ztl+tVlssFou6ufMvW7bs5+8jIiIQ0cLS6bZSzk03tfWn\nICKSz9NPA2o18OST9h1vMBhgMBjadEybk31jCoVCUigUDh272TjZX83EiUBaGvDcc46MgIio42zc\nCGzeLKoS9i5EcuWNcFJS0lWPafNoHJVKZa2srPQGgIqKCh8vL6/jgLhjLysr87W9z2w2azQajVmt\nVlvMZrOm8X61Wm1p63WbEh4u2ibU1zvibEREHevQIeAPfxANzvr3v/r7HanNyT46OjonPT09HgDS\n09PjY2Jismz7MzIyZtfW1nYvKSkZZjKZdKGhoQXe3t6Vffv2PZOfnz9BkiTF+vXrH7Qd015DhgA+\nPsD+/Y44GxFRxzl/XixEsnw5MGaMDAFIktTsNnv27A0+Pj7lSqWyVqPRlK1du3beyZMnB06aNGmb\nTqcrjoyMzDt16lR/2/uXL1/+Rz8/v8PDhw8vys3NnWLbv3v37rGjRo064Ofnd/jRRx9d09z1RDht\ns3ChJK1a1ebDiIg6TUODJM2dK0kPPCC+d7SfcmeL+dwl2yU09u67wEcfAR9+2EFBERG105tvir43\n+flAr16OP39r2iW4fLIvKxMfiY4ft/9hBxFRR9m3D4iMFCMHR4zomGt02d44jfn6Ar17i+W7iIic\nyenTok6/Zk3HJfrWcvlkD7B1AhE5H0kC5s0DpkwRjc7kxmRPRNQBVq0SZeZ//EPuSASXr9kDgMkE\n3H67+A/Luj0Rye2bb4CYGPFAdtiwjr+eW9TsAcDfX0ysKi2VOxIicncnTgCzZ4vZ/Z2R6FurSyR7\nhUKUcrZvlzsSInJn9fXAAw8Ac+YA06fLHc2vdYlkD7BuT0Ty+9vfgIsXxVdnw2RPROQAW7cC//43\nkJEhVp5yNl0m2Y8cCZw8CVRUyB0JEbkbsxmYO1fM6PfxkTuapnWZZN+tm1i9inV7IupMdXXArFmi\nm+Xtt8sdTfO6TLIHWMohos63ZIlYOW/JErkjaZkTVpbsN3Ei8PbbckdBRO7i44+BDz4QC5F0c/Jb\n5y4xqcqmrg4YNEiMtx840HFxERFd6cgRsSTqpk3AhAnyxuI2k6pslEogLAz46iu5IyGiruzCBdHg\n7M9/lj/Rt5bdyX716tWLgoKCDowaNeq71atXLwKAqqqqgZGRkVsDAgKKo6Ki8qqrq39eeCs5OXmJ\nTqczBQYGFuXl5UU5IvimsG5PRB1t0SJApxMPZV2FXcn+u+++G/Xmm28u2LVr1/hvv/02+JNPPvnN\nkSNH/FJSUhZHRkZuLS4uDpg0adJnKSkpiwHAaDTqN27cOMtoNOpzc3OnJiYmpjY0NHTIpwomeyLq\nSOvXAwaDWJDElXpx2ZVwi4qKAidMmJDfo0ePix4eHvW33XbbFx9++OG9OTk50fHx8ekAEB8fn56V\nlRUDANnZ2XfHxcVtUCqVdVqtttTf3/9wQUFBqCN/EJvQUMBoBM6e7YizE5E7++474Mknxcp4ffvK\nHU3b2JXsR40a9d327dvDq6qqBp4/f77n5s2bp5nNZo3ValWpVCorAKhUKqvValUBQHl5+VCNRmO2\nHa/RaMwWi0XtmB/h13r0ECtX7djREWcnIndVUyPq9C+9BAQFyR1N29k19DIwMLDoueeeWxkVFZXX\nq1evc6NHj97n4eFR3/g9CoVCUigUzQ6tae61ZcuW/fx9REQEIiIi2hyfrZQTGdnmQ4mI/ockAQsX\niombDz0kdzSAwWCAwWBo0zF2j7NPSEhYm5CQsBYAli5dulyj0ZhVKpW1srLS29vbu7KiosLHy8vr\nOACo1WpLWVmZr+1Ys9msUavVlqbO2zjZ22viRGD58nafhogIAPDaa6I8vHOn3JEIV94IJyUlXfUY\nux+SHj9+3AsAjh07dt1HH310z5w5c96Ljo7OSU9PjweA9PT0+JiYmCwAiI6OzsnIyJhdW1vbvaSk\nZJjJZNKFhoYW2Hvtq7npJjHJ4eLFjroCEbmLXbuAF14Qk6euvVbuaOxn9539zJkzPzh58uQgpVJZ\nl5qamtivX7/TixcvTomNjc1MS0ubr9VqSzMzM2MBQK/XG2NjYzP1er3R09PzcmpqamJLJZ726tMH\n0OvF/6Tw8I66ChF1dVVVQGws8PrrYqilK+tSM2gbe+opMYt26VKHnI6I3ExDA3D33WIlvH/+U+5o\nWuZ2M2gb43h7ImqPF18EfvwRWLlS7kgco8ve2Z88KdZ/rKpyzoUEiMh5ffGFaFu8axfg63v198vN\nre/sBw0Crr8e2LtX7kiIyJVUVoo1ZN96yzUSfWt12WQPsJRDRG1z+TIQFwfMnw9MnSp3NI7FZE9E\n9JMXXgA8PMTXrqbL1uwBoLwcGDVKPGRx9oUFiEhen34KPPwwUFgIeHnJHU3buHXNHgCGDhW1+4MH\n5Y6EqHlffQXMmwesWweUlckdjXs6ehRISAAyMlwv0bdWl072gCjlcBFyclZnzwIPPACo1UBuLhAS\nAgwfDvz+92LJu+pquSPs+i5dEhOnnn4auPVWuaPpOF26jAMA6enAli3iNzaRs3n8ceDUKfH3FBAT\neb79Fti2TWzffCNmg0+eLLabbwauuUbemLuChgbg0CHRHfeDD4Du3cUvV1fqT99Ya8o4XT7Zl5SI\nTnUWi+v+j6SuaedOYMYM0SN90KCm33PxokhItuRvNIqEb0v+wcF8HtUaZ84A+fniv+WOHeL7fv1E\nH62bbhKdLPv0kTtK+zHZQ7Qmve464L//FdOeiZxBba1Yd+FPfwJmz279cdXVYpUkW/I/eRK4445f\nkv+wYR0WsstofNe+c6f4WlIi/nuHhf2S4L295Y7UcZjsf3L//cCkSeIBDJEz+MtfxOzMnJz2feIs\nKwM+++yX5N+r1y+J//bbgcGDHRezs7raXftNN4lPQEql3JF2HCb7n7z+uvgN/9ZbDj81UZsZjcBt\nt4nZ3RqN484rSWLkmS3xf/ml+DRrS/633gr07Om468nBHe/aW4PJ/idGI/Cb3wA//ODwUxO1SUOD\nSLoPPAAkJnbstWprgYKCX5L/vn1ijWZb8h87Vkwgcma8a28dJvufSJIYO1tY2LV6XZDr+de/gI0b\nRaOtzn6weuaMuNu3JX+LRZR6bMlfp5N3EAPv2u3HZN/IPfeIxYLnzOmQ0xNd1bFjInF99RUQGCh3\nNEBFBfD55yLxb90qfvnYEv+kSYBK1bHX512743Rosk9OTl7yzjvvPNCtW7eGoKCgA+vWrZt37ty5\nXrNmzdp49OjR620rVfXv37/a9v61a9cmeHh41K9Zs+axqKiovCYC7rBkv2oVUFQk6vdEnU2SRCnx\nppvECBxnI0lAcfEvd/0Gg/gUbEv+EycCvXvbf37etXesDkv2paWl2jvuuOPz77//fsQ111xzadas\nWRunTZu2+eDBgyMHDx7847PPPvviypUrnzt16tSAlJSUxUajUT9nzpz3du3aNd5isagnT568rbi4\nOKBbt24NVwTcYcm+sFDUSY3GDjk9UYveew9ISQF27xYTeJzd5ctiHWdb8t+1SyTmSZNE8g8NbfmO\nm3ftnas1yd6uZT369u17RqlU1p0/f76nh4dH/fnz53sOHTq0PDk5eckXX3xxGwDEx8enR0REGFJS\nUhZnZ2ffHRcXt0GpVNZptdpSf3//wwUFBaFhYWGdtlZ7cLCoUZ44AQwZ0llXJRKN+J58UgyzdIVE\nD4gFfyZMENvSpcC5c6L8tG0b8Ic/iMEOEyf+cuffrVvzd+0LF4q+P7xrl5ddyX7gwIFVTz311MvX\nXXfdsWuvvfbClClT/hMZGbnVarWqVCqVFQBUKpXVarWqAKC8vHxo48Su0WjMFotF7ZgfoXU8PMRM\n2q++ErMWiTrLE0+IZ0WhoXJHYr9evYApU8QGiJsmW71/9WpRBrLdsS9cyLt2Z2RXsj9y5IjfqlWr\nHi8tLdX269fv9H333ff+O++880Dj9ygUCkmhUDRbk2nutWXLlv38fUREBCIiIuwJsUm2/vZM9tRZ\ncnPFDcZ338kdiWMNGSKW7Zs1S+5I3JPBYIDBYGjTMXYl+927d4+7+eabvxk0aNBJALjnnns+2rFj\nx03e3t6VlZWV3t7e3pUVFRU+Xl5exwFArVZbysrKfh70aDabNWq12tLUuRsne0ebOBF49NEOOz3R\nr9TUAL/7HfB//yfujIkc5cob4aSkpKseY9dI38DAwKKdO3eGXbhw4VpJkhTbtm2brNfrjdOnT9+U\nnp4eDwDp6enxMTExWQAQHR2dk5GRMbu2trZ7SUnJMJPJpAsNDS2w59rtMW6cGHFw+nRnX5nc0Z/+\nBEREAFFRckdCZOedfXBw8Ldz5859e9y4cbu7devWMGbMmMKFCxf+X01NTZ/Y2NjMtLS0+bahlwCg\n1+uNsbGxmXq93ujp6Xk5NTU1saUST0fp3h0YP160jb3zzs6+OrmTHTuAzMyuV74h1+U2k6psli0T\nixUkJ3foZciNXbokRqI8/zxr2tQ53H5ZwqaEh3MRcupYKSmAn59Y/YjIWbjdnf25c6JPzokTrt8B\nkJzPwYOiTu/ojpZELeGdfRN69QJuvFHM6CNypPp6YMEC0aueiZ6cjdsle+CX8fZEjpSaKiYSPfyw\n3JEQ/S8meyIHOHoUSEoC3niDa8KSc3K7mj0g1vH09RXrd7pKrxJyXpIE3HWXaMexdKnc0ZA7Ys2+\nGf37i+Xa9uyROxLqCt57TzTZe/ZZuSMhap5bJntAlHK2b5c7CnJ1J04ATz0FvPkmG3+Rc3PrZM+6\nPbXXE08A998vZmYTOTO3rNkDwPHjwPDhote4sy+6TM5pyxbR233/fjY6I3mxZt8CLy+xmMKBA3JH\nQq7I1tHy3/9moifX4LbJHmAph+y3dClwxx1ilSYiV8Bkz2RPbbRjB/DBB8DLL8sdCVHruXWytzVF\nc6LHFuTkLl0C5s8XS/ENHCh3NESt59bJ/rrrRDO0Q4fkjoRcRXIyoNMBM2fKHQlR29i1eElXYivl\nBAbKHQk5u4MHgVdfBfbtAxQtjnsgcj523dkfOnRoeEhIyF7b1q9fv9Nr1qx5rKqqamBkZOTWgICA\n4qioqLzq6ur+tmOSk5OX6HQ6U2BgYFFeXp7TLNTGuj21Rn29KN/87W+AWi13NERt1+5x9g0NDd3U\narWloKAg9JVXXnl08ODBPz777LMvrly58rlTp04NSElJWWw0GvVz5sx5b9euXeMtFot68uTJ24qL\niwO6devW8KtgOnGcvU1xMTBpEnDsGO/WqHlr1gAffgj8979sdEbOp1PG2W/btm2yv7//YV9f37Kc\nnJzo+Pj4dACIj49Pz8rKigGA7Ozsu+Pi4jYolco6rVZb6u/vf7igoCC0vdd2BJ0OqKsTXQuJmnL0\nqOhRz46W5MraXbPPyMiYHRcXtwEArFarSqVSWQFApVJZrVarCgDKy8uHhoWF7bQdo9FozBaLpckP\nw8uWLfv5+4iICERERLQ3xBYpFL/0ydFqO/RS5IIkSfSnf+opICBA7miIBIPBAIPB0KZj2pXsa2tr\nu2/atGn6ypUrn7vyNYVCISkUimZrMs291jjZdxZb3f7BBzv90uTk3n0XqKwEnn5a7kiIfnHljXBS\nUtJVj2nXh9ItW7bcOXbs2D1Dhgw5AYi7+crKSm8AqKio8PHy8joOAGq12lJWVuZrO85sNmvUarWl\nPdd2JD6kpaawoyV1Je1K9hs2bIizlXAAIDo6Oic9PT0eANLT0+NjYmKybPszMjJm19bWdi8pKRlm\nMpl0oaGhBe0L3XFGjRL/sCsr5Y6EnMmiReLT3rhxckdC1H52j8Y5d+5cr+uvv/5oSUnJsD59+tQA\nQFVV1cDY2NjMY8eOXafVakszMzNj+/fvXw0AK1as+OPatWsTPD09L69evXrRlClT/vM/wcgwGscm\nOlr8w77vPlkuT07m00+Bxx4TjfJ69pQ7GqKWtWY0jtu2OL7S3/8uRl288ooslycnUlMjPu2tXSuG\n5RI5O7Y4bgPW7cnmj38USZ6JnroS3tn/pK5ONLY6epQNrtzZN9+IvjcHDwIDBsgdDVHr8M6+DZRK\nICwM+PpruSMhuVy6BCxYIGbLMtFTV8Nk3whLOe5txQoxceree+WOhMjx3L7rZWMTJwLPPit3FCSH\n774DUlPZ0ZK6LtbsG7lwARg8GLBagd69ZQuDOll9PXDLLUBCArBwodzRELUda/ZtdO21wJgxYtk5\nch//+hdwzTWiXk/UVTHZX8HWFI3cQ2kp8Ne/sqMldX38630FPqR1H7aOlk8/zY6W1PWxZn+FmhrA\nxwc4eVJ8tKeua/164B//AAoK2OiMXBtr9nbo0wcYMQLYtUvuSKgjHT8u7ujZ0ZLcBZN9E1jK6foW\nLQLi44GxY+WOhKhzMNk3gcm+a/vkE/HJTYZ1cohkw5p9E06eBG64QXz15LSzLuXMGdHR8q23gDvu\nkDsaIsdgzd5OgwYBvr5iNiV1LUuWAFFRTPTkfnjf2gxbKYerFHUdX38NZGWJ1ghE7sbuO/vq6ur+\nM2fO/GDEiBHf6/V6Y35+/oSqqqqBkZGRWwMCAoqjoqLyqqur+9ven5ycvESn05kCAwOL8vLyohwT\nfsdh3b5ruXiRHS3Jvdmd7BctWrR62rRpm7///vsR+/fvvzEwMLAoJSVlcWRk5Nbi4uKASZMmfZaS\nkrIYAIx4rHPbAAAQoUlEQVRGo37jxo2zjEajPjc3d2piYmJqQ0ODU5eQwsPFTNqGBrkjIUdYvlwM\nqWVHS3JbkiS1eauuru43bNiwH67cP3z48KLKykqVJEmoqKjwHj58eJEkSVixYsWSlJSU52zvmzJl\nSu6OHTvCrjxehOM8/Pwk6cABuaOg9tq/X5IGD5Yki0XuSIg6xk+5s8W8bVfNvqSkZNiQIUNOzJs3\nb923334bPHbs2D2rVq163Gq1qlQqlRUAVCqV1Wq1qgCgvLx8aFhY2E7b8RqNxmyxWNRNnXtZo/Fw\nERERiIiIsCdEh7D1yRk1SrYQqJ3q60X5ZsUKYOhQuaMhcgyDwQCDwdCmY+xK9pcvX/YsLCwc869/\n/esP48eP3/X444+vspVsbBQKhaRQKJodR9nca8ucaPDzxInAf/4DPPKI3JGQvV55RXQznT9f7kiI\nHOfKG+GkpKSrHmNX3Vyj0Zg1Go15/PjxuwBg5syZHxQWFo7x9vaurKys9AaAiooKHy8vr+MAoFar\nLWVlZb62481ms0atVlvsuXZnsj2kdYKh/2SHkhLgb39jR0siwM5k7+3tXenr61tWXFwcAADbtm2b\nPHLkyIPTp0/flJ6eHg8A6enp8TExMVkAEB0dnZORkTG7tra2e0lJyTCTyaQLDQ0tcNyP0TGGDROr\nFv3wg9yRUFvZOlo+8wyg08kdDZH87B5n/8orrzx6//33v1tbW9vdz8/vyLp16+bV19d7xMbGZqal\npc3XarWlmZmZsQCg1+uNsbGxmXq93ujp6Xk5NTU1saUSj7NQKH65u/fzkzsaaov164ETJ4Ann5Q7\nEiLnwHYJV/H660B+PrBundyRUGsdPw4EBQFbtoiVx4i6OrZLcABOrnI9jz0GPPQQEz1RY2yXcBUj\nRojmWWYzoNHIHQ1dzaZNwO7dwNq1ckdC5Fx4Z38VCgVw661cl9YVnDkD/P73YvRNz55yR0PkXJjs\nW4GlHNeweDEwZQpw++1yR0LkfFjGaYWJE8XydeS8vvoKyM5mR0ui5vDOvhWCg0XN/sQJuSOhK508\nKRL9ggVitiw7WhI1jXf2reDpCdx8s0gqM2bIHY37kSTxy/b77/93u3RJPESfPRu45x65IyVyXkz2\nrWRrisZk33EuXxazlY3GXyf0oiKgVy+R1EeMAEaOBO67T3zv4yMeohNRyzipqpW+/hpYtEgM66P2\nuXABOHTof+/SjxwRyduW1BtvLM8QNa81k6qY7Fvp0iWxNm15OdC3r9zRuIZTp5ouvVRUiPYTVyb0\ngAAOmSSyR2uSPcs4rXTNNcD48cA33wBTp8odjfOQJPELsKmkfv48EBj4SzJfuFB8veEG8RyEiDoP\n/8m1gW28vTsm+/p60TLYlshtdfWiIqBHj1/foc+YIb6q1aynEzkLlnHa4LPPgBdeEKNyuqqLF4Hi\n4v+9Sz98GFCpmq6nDxwod9RE7o01ewc7d04kvBMnxOpHXcX27cA//wns3w9YLKKPf+NkrtcDw4ez\nnk7krFizd7BevUTr3Px8QMalcR1m717gj38UpZilS8U6rX5+gFIpd2RE5GicQdtG4eGu3yfn0CEg\nNha46y5g+nTx5wULxMNUJnqirsnuZK/VaktvvPHG/SEhIXttSwxWVVUNjIyM3BoQEFAcFRWVV11d\n3d/2/uTk5CU6nc4UGBhYlJeXF+WI4OXgyk3Rjh0TSf3WW0Wvd5MJSEwEuneXOzIi6mh2J3uFQiEZ\nDIaIvXv3hhQUFIQCQEpKyuLIyMitxcXFAZMmTfosJSVlMQAYjUb9xo0bZxmNRn1ubu7UxMTE1IaG\nBpf8VHHLLaKMU1srdyStd+IE8MQTQEiIeOZQXCw6RPbqJXdkRNRZ2pVwr3wgkJOTEx0fH58OAPHx\n8elZWVkxAJCdnX13XFzcBqVSWafVakv9/f0P235BuJoBA0Rdu7BQ7kiu7vRp4PnnRXmmvh44eBBY\nvpyzUYnckd0PaBUKhTR58uRtHh4e9Q8//PC/f/vb375htVpVKpXKCgAqlcpqtVpVAFBeXj40LCxs\np+1YjUZjtlgs6qbOu2zZsp+/j4iIQIQTPgm19ckJC5M7kqZduAC8+irw4ovAtGnAnj2AVit3VETk\nKAaDAQaDoU3H2J3sv/7661t8fHwqTpw4MSQyMnJrYGBgUePXFQqFpFAomh1H2dxrjZO9s5o4EUhP\nB555Ru5Ifq2uTizH99e/AhMmAAaDGDZJRF3LlTfCSUlJVz3G7jKOj49PBQAMGTLkxIwZMz4uKCgI\nValU1srKSm8AqKio8PHy8joOAGq12lJWVuZrO9ZsNmvUarXF3mvLLTxcTKyqr5c7EqGhAXjvPTEm\n/sMPgY8/Fl+Z6InIxq5kf/78+Z41NTV9AODcuXO98vLyooKCgg5ER0fnpKenxwNAenp6fExMTBYA\nREdH52RkZMyura3tXlJSMsxkMulsI3hckUolNrlXRZIkscD26NFi4Y433gDy8kQPHyKixuwq41it\nVtWMGTM+BoDLly973n///e9GRUXljRs3bndsbGxmWlrafK1WW5qZmRkLAHq93hgbG5up1+uNnp6e\nl1NTUxNbKvG4AtsQzOBgea5vMIgJUTU14qHr9OnsQ0NEzWO7BDu9845Y8/T99zv3urt3i9muhw8D\nf/mLWKHJw6NzYyAi59KadgkuOdbdGdju7Dvrd9P33wMzZwJ33y26Sn7/PXD//Uz0RNQ6TPZ2uu46\n0QytuLhjr3P0KDBvHnDbbUBoqJj1+rvfcdYrEbUNk307dGTrBKsVeOwx0dZAoxFJ/tln2XmSiOzD\nZN8OHdEUrbpa1OT1eqBbN1Gu+etfgX79HHsdInIvTPbt4Mg7+/PngZQUQKcDKitF++FVqwAvL8ec\nn4jcG5N9OwQEiJWdjh61/xy1tUBqKuDvL/rtbN8OpKWJZwJERI7CZN8OCoX9d/f19cD69aJJ2aZN\nYsvMFH8mInI0Jvt2sjVFay1JArKyxGSs118H3noL2LIFGDu2w0IkIuKyhO01caLoMNkan38uZr1e\nuACsXCk6UnLWKxF1Bs6gbaf6emDwYLGOq0rV9HsKCkSSP3pUzHqdNUuMtCEicgTOoO0EHh5imb+m\nSjkHD4rZrvfcI9Z8NRqBuDgmeiLqfEw7DnDlQ9qSEmDuXOCOO8QvApMJWLiQi3kTkXyY7B3Aluwr\nKoA//AEYNw644QaR5J96SrRVICKSE2v2DlBXBwwcKPrVPPSQWMx7yBC5oyIid9Gamj2TvYN8+ilw\n442Ar+/V30tE5Egd/oC2vr7eIyQkZO/06dM3AUBVVdXAyMjIrQEBAcVRUVF51dXV/W3vTU5OXqLT\n6UyBgYFFeXl5Ue25rjO66y7gyBGD3GHYra2LFzsbxi8vxu/82pXsV69evUiv1xttq06lpKQsjoyM\n3FpcXBwwadKkz1JSUhYDgNFo1G/cuHGW0WjU5+bmTk1MTExtaGjocs8LXPkvjCvHDjB+uTF+52d3\nwjWbzZrNmzdPW7BgwZu2jw85OTnR8fHx6QAQHx+fnpWVFQMA2dnZd8fFxW1QKpV1Wq221N/f/3BB\nQUGoY34EIiK6GruT/RNPPPHPl1566Zlu3bo12PZZrVaVSqWyAoBKpbJarVYVAJSXlw/VaDRm2/s0\nGo3ZYrGo2xM4ERG1nl3tEj755JPfeHl5HQ8JCdlrMBgimnqPQqGQWlpUvLnXFC7ePyApKUnuEOzm\nyrEDjF9ujN+52ZXsv/nmm5tzcnKiN2/ePO3ixYs9zpw50/fBBx9cr1KprJWVld7e3t6VFRUVPl5e\nXscBQK1WW8rKyn4ep2I2mzVqtdpy5Xmv9jSZiIjsJElSuzaDwXDbb37zm02SJOGZZ555MSUl5TlJ\nkpCcnLz4ueeeS5EkCQcPHtQHBwfvu3TpUvcffvhh2A033HCkoaFB0d5rc+PGjRu31m0O6XppK8ks\nXrw4JTY2NjMtLW2+VqstzczMjAUAvV5vjI2NzdTr9UZPT8/LqampiS2VeIiIyMHk/m0jSRK2bNky\ndfjw4UX+/v4m2ycDV9nmzZu31svLyzpq1KgDcsdiz3bs2DHfiIiI/+r1+oMjR478bvXq1Y/JHVNb\ntgsXLvQIDQ3NDw4O3jdixAjj4sWLk+WOqa3b5cuXPUaPHr3X9gnZ1bbrr7++NCgoaP/o0aP3jh8/\nvkDueNqynTp1qv+99977QWBg4PcjRoww7tixI0zumFq7FRUVDR89evRe29a3b9/TLf37lT3gy5cv\ne/j5+R0uKSnR1tbWKoODg/cZjcYRcsfV2u3LL78MLywsDHHVZF9RUeG9d+/e0ZIkoaampndAQMAh\nV/rvL0kSzp0711OSJNTV1XlOmDBh5/bt22+VO6a2bC+//PKTc+bMeXf69Ok5csdiz6bVaktOnjw5\nUO447Nnmzp2bnpaWliBJ4u9PdXV1P7ljsmerr6/v5u3tXXHs2DHf5t4j+8SmgoKCUH9//8NarbZU\nqVTWzZ49OyM7O/tuueNqrfDw8O0DBgw4JXcc9vL29q4cPXr0PgDo3bv32REjRnxfXl4+VO642qJn\nz57nAaC2trZ7fX29x8CBA6vkjqm1mpqv4opcMfbTp0/32759e3hCQsJaAPD09Lzcr1+/03LHZY9t\n27ZN9vPzO+Lr61vW3HtkT/YWi0XdOECOwZdPaWmpdu/evSETJkzIlzuWtmhoaOg2evTofSqVynr7\n7bf/V6/XG+WOqbWamq/iahQKhTR58uRt48aN2/3GG2/8Vu54WqukpGTYkCFDTsybN2/dmDFjCn/7\n29++cf78+Z5yx2WPjIyM2XPmzHmvpffInuz5oNY5nD17tvfMmTM/WL169aLevXuflTuetujWrVvD\nvn37RpvNZs2XX345sbm5H86m8XwVV7wztvn6669v2bt3b8iWLVvufPXVV3+/ffv2cLljao3Lly97\nFhYWjklMTEwtLCwc06tXr3O2Fi+upLa2tvumTZum33fffe+39D7Zk/2VY/DLysp8G8+2pY5XV1en\nvPfeez984IEH3omJicmSOx579evX7/Rdd9316e7du8fJHUtr2OarDBs2rCQuLm7D559/fsfcuXPf\nljuutvLx8akAgCFDhpyYMWPGx67SCkWj0Zg1Go15/PjxuwBg5syZHxQWFo6RO6622rJly51jx47d\nM2TIkBMtvU/2ZD9u3LjdJpNJV1paqq2tre2+cePGWdHR0Tlyx+UuJElSzJ8/P02v1xsff/zxVXLH\n01Y//vjjYFt31QsXLly7devWyJCQkL1yx9UaK1as+GNZWZlvSUnJsIyMjNl33HHH52+//fZcueNq\ni/Pnz/esqanpAwDnzp3rlZeXFxUUFHRA7rhaw9vbu9LX17esuLg4ABB175EjRx6UO6622rBhQ1xc\nXNyGq75R7qfIkiRh8+bNdwYEBBzy8/M7vGLFiiVyx9OWbfbs2Rt8fHzKu3fvfkmj0ZStXbt2ntwx\ntWXbvn37rQqFoiE4OHifbQjXli1bpsodV2u3/fv3B4WEhBQGBwfvCwoK2v/iiy8+I3dM9mwGg+E2\nVxyN88MPPwwLDg7eFxwcvG/kyJHfudq/33379gWPGzdu14033vjtjBkzPnK10Thnz57tNWjQoB/P\nnDnT52rvdarFS4iIqGPIXsYhIqKOx2RPROQGmOyJiNwAkz0RkRtgsicicgNM9kREbuD/ARQeabXQ\n0NzRAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x4628110>"
       ]
      }
     ],
     "prompt_number": 73
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "References\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Cross validation:\n",
      "http://en.wikipedia.org/wiki/Cross-validation_%28statistics%29\n",
      "\n",
      "Coefficient of determination:\n",
      "http://en.wikipedia.org/wiki/Coefficient_of_determination\n",
      "\n",
      "Testing\n",
      "* http://en.wikipedia.org/wiki/Jarque%E2%80%93Bera_test#Jarque.E2.80.93Bera_test_in_regression_analysis\n",
      "* http://en.wikipedia.org/wiki/Lack-of-fit_sum_of_squares\n",
      "\n",
      "Model checking / specification in fMRI:\n",
      "* https://web.njit.edu/~loh/Papers/Sinica.Loh.etal.2008.pdf\n",
      "* http://www.sciencedirect.com/science/article/pii/S1053811902910943\n",
      "* http://www.sciencedirect.com/science/article/pii/S1053811908012056\n",
      "* http://www.sciencedirect.com/science/article/pii/S1053811903001496\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}