week4_tutorial.ipynb

{
 "metadata": {
  "name": "week4_tutorial"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Psych216A Tutorial 4: Model Fitting"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This tutorial will give an introduction to linear and non-linear fitting\n",
      "proceedures. It covers the basic algebra of linear regression and\n",
      "compares the solutions from ordinary least squares regression to the\n",
      "solutions obtained from a non-linear fitting proceedure.\n",
      "\n",
      "Written by Jason D. Yeatman jyeatman@stanford.edu\n",
      "\n",
      "Translated to Python by Michael Waskom"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].\n",
        "For more information, type 'help(pylab)'.\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "import utils\n",
      "import seaborn\n",
      "seaborn.set()\n",
      "colors = seaborn.color_palette()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Linear Models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As was emphasized in class, linear models can be fit using simple\n",
      "algebra.  This tutorial will briefly review the math underlying the\n",
      "linear model.  The optimal weights to minimize the squared difference\n",
      "between the model prediction and the measured data can be calculated\n",
      "with he following formula:\n",
      "\n",
      "$$w = (X^TX)^{-1}X^Ty$$\n",
      "\n",
      "Where w is the weights or model parameters to be estimated, X is the\n",
      "regressor matrix and y is the dependent variable.  The weights obtained\n",
      "from this formula provide the best prediction of y as a weighted sum of\n",
      "the regressors in X:\n",
      "\n",
      "$$y = Xw$$\n",
      "\n",
      "Understanding the ordinary least squares formula requires understanding\n",
      "matrix inversion.  The inverse of matrix, A, will reverse the effects of\n",
      "multiplying by A."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# For example if we have some data X\n",
      "X1 = randn(100, 2)\n",
      "plot(X1[:, 0], X1[:, 1], \"o\");"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAFRCAYAAAA1uqfwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9wVNX9//HXZpNgcHCVIlD8UO1MqcRADfaDUD+VAIUI\nAkUTNTA2MgJCpwMpf4iO4/ht63xrR50OBRln+Ej90SIiEuTTyphGDZ9l7IxU6o8y8Udx5mNrndAv\nPz6k/gDJJvf7ByZkk93s3d2755579/mY6UxdNtmzZzf3dc8573tuxHEcRwAAwJgSvxsAAECxIXwB\nADCM8AUAwDDCFwAAwwhfAAAMI3wBADAsr/A9ffq0pk2bpurqak2fPl0bNmzwql0AAIRWJN/rfD//\n/HMNHz5cX3zxhb797W9rz549+sY3vuFV+wAACJ28p52HDx8uSfr000+VSCQ0bNiwvBsFAECY5R2+\nPT09uvLKKzVmzBitWbNG48eP96JdAACEVt7Tzr0+/PBDXX/99Xr66ac1ZcqUvsd7ehwlEt1evERo\nlZZGJYl+coG+cod+co++cod+cqe0NKqSkkjm53n1gpdddpmuv/56HThwICl8E4ludXae8uplQikW\nq5Ak+skF+sod+sk9+sod+smdWKxC5eWZozWvaedjx47p5MmTkqTjx4+rtbVVixcvzudXAgAQenmN\nfDs6OrRs2TJ1d3dr7NixuvPOO/XVr37Vq7YBABBKeYXv5MmT9cYbb3jVFgAAigI7XAEAYBjhCwCA\nYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABjm2f18\ngUza4nHtaH5B3U5E0YijJfULNbumxu9mAYBxhC+MaIvHtWnrTo2sqpckJSRt2rpTkghgAEWHaWcY\nsaP5hb7g7TWyql7P7t7rU4sAwD+EL4zodiIpH0/0GG4IAFiA8IUR0YiT8vFSvoEAihBrvhYJc0HS\nkvqFSWu+knSifZeaVjb42CoA8Afha4nWl/eFuiCp9z08u3uvEj1nR7xNKxtC8d4AIFsRx3FSzwd6\n5MyZhDo7TxXyJQIvFqvQrSt+rNOjawf9W1lHi7ZsfMiHVtkpFquQJL5TGdBP7tFX7tBP7sRiFSov\nzzyuZcXNEt1pCo8oSAKA8CF8LRFN80lQkAQA4cOh3RLLlt6gE+3NSY+daN+lhroFPrUIAFAoFFxZ\nonbOLH32+RcUJAFAESB8LTK7poawBYAiwLQzAACGEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG\n+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACG\n5RW+H330kWbNmqWqqirNnDlT27dv96pdAACEVmk+P1xWVqYNGzaourpax44d09VXX61FixZpxIgR\nXrUPAIDQyWvkO3bsWFVXV0uSRo0apaqqKh08eNCThgEAEFYRx3EcL37RBx98oNraWh06dEjnn39+\n3+M9PY4SiW4vXiK0SkujkkQ/uUBfuUM/uUdfuUM/uVNaGlVJSSTj8zwpuPrkk0/U0NCgDRs2JAUv\nAAAYLK81X0nq6upSfX29GhsbtXjx4kH/nkh0q7PzVL4vE2qxWIUk0U8u0Ffu0E/u0Vfu0E/uxGIV\nKi/PHK15jXwdx9GKFSs0adIkrVu3Lp9fBQBA0cgrfP/4xz9q27Ztamtr05QpUzRlyhS1tLR41TYA\nAEIpr2nn7373u+rp6fGqLQAAFAV2uAIAwDDCFwAAwwhfAAAMI3wBADCM8AUAwDDCFwAAwwhfAAAM\nI3wBADCM8AUAwDDCFwAAwwhfAAAMI3wBADCM8AUAwLC87moEAG3xuHY0v6BuJ6JoxNGS+oWaXVPj\nd7MAqxG+AHLWFo9r09adGllVL0lKSNq0dackEcDAEAhfFC3TI7YwjhB3NL/QF7y9RlbV69ndewP/\n3oBCInxRlEyP2MI6Qux2IikfT/QYbggQMBRcoSgNNWILw+uZEo04KR8v5cgCDIk/ERQl0yO2sI4Q\nl9Qv1In25qTHTrTvUkPdAp9aBAQD084oStGIo0SKxws1YjP9eqb0Tpk/u3uvEj1n30/TyoZAT6UD\nJhC+KEpL6hcmrcFKZ0dsTSsbQvF6Js2uqSFsgSwRvihKpkdsjBAB9BdxHCd1xYRHzpxJqLPzVCFf\nIvBisQpJop9coK/c8bOfgnZJFd8pd+gnd2KxCpWXZx7XMvIF4JmwXlIFeC3g5R4AbBLWS6oArxG+\nADwT1kuqAK8x7QzAM2G9pKpQgrY+Du8QvkBItb68T7/+bbPRA3uYL6nyGuvjxY3wBQwxOcppfXmf\nHty8TRdW1kkyd2APwiVVAz+HFY31qp0zy3g7uClFcSN8AQNMj3KeemZPX/D2MnVgt3nTjVSfw4Ob\nt0mSpk2dbrQtrI8XN1ZiAANMVwF3pzmAF/uBPdXncGFlnX6zY4/xtnBTiuLGxwwYYHqUE03zl13s\nB/a0n0O34YaIm1IUuyL/UwTMMD3KWbb0Bp18d3fSYxzYh/gcooYborPT800rb1FZR4siH7eorKPF\nuvVxFA5rvoABqaqA3/vvLbpl4bUFeb3eAqLHt+22tvDJD6k+h5PvNGv12kZf2mPz+jgKi72dLcCe\nqe4Fua9+tXmznm5u0bALxsnp6daoS6sV6XxPTStv8fwAHOR+KrS2eDypGnv5D+pUO2cWfZUB3yl3\n3O7tTPhagC+1e0Hqq4GXtBw9+v900ZXLBj2vrKNFWzY+5OlrB6mf/EZfpdf/OzysrETLlt5gvCo8\naNyGL2u+QAH0XtKSGDdfziXzlBg3X//839M6/o/2Qc8t9gpk2Gngd/j06Fo9uHmb2uJxv5sWCoQv\nUACpLmmZOHO1jv3trUHPLfYKZNgp3WVZ3CTDG/zZAwWQ9pKWUyeS/psKZNiKTUAKi2pnoADS3WBg\n/Fe/orKOFiqQYT1uklFYhC9QAOlvMNBI2CIQ0l2WtWbFLT62KjwIX6AAgnCDAWAoA7/D55WXaPXa\nRqqdPZLXpUbLly/X3r17NXr0aB06dCjlc7jUKDMudXCPvnKHfnKPvnKHfnLHyKVGt99+u1paWvL5\nFQAAFJ28wvfaa6/VRRdd5FVbAAAoCgVf8y0tjfZNVyC10i93daefMqOv3KGfUmt9eZ+eemaPunvO\n3vlp2dIbdP28OZLoq0z4TrlT6vIuHRRcASgKrS/v04Obt+nCyjpJUpekBzdvUzQa0XVzv+dv41B0\nCh6+iUQ3C/QZDCxkGLgn8JL6hVTJfomiD3fop8F+/dvmvuDtdWFlnZ54+nl9b9ZM+ioDvlPuuC24\nYuRrmd79VHuvrUtI2rR1pyQRwEAe0u7Y1G24IYDyLLhaunSprrnmGv31r3/V+PHj9cQTT3jVrqKV\naj/VkVX17KcK5CkaSX1VpcslOsBTeY18n3nmGa/agS+xnyq80FtY9EVXD0sXX0q369g9Tbf52CoU\nK6adLcN+qv4Kw3p7WzyuzY8/17e+ydLFWel2HaudM8vnlqEYEb6WSb8ncIOPrSoOYVlv39H8wqDC\not6liyC9j0KYXVNT9H0AOxC+lmFPYP8Mtd4epP4v1NJFGGYFAFsQvhbi7NwfYVlvL8TSRVhmBQBb\nsJIIfCltNWzA/kqW1C/UyXd3Jz12on2XGuoW5Pw7qcIHvBWwwwpQOEvqF+pEe3PSY/mGlh9m19To\n7jU/UMXRVkU+blFZR0veSxdhmRUAbMG0M/ClMK23186Zpdo5swbtRpTrui1V+IC3CF+gnzCvt+ez\nbksVPuCtiOM4qRe6PHLmTIK9QDNgz1T3bO8rWyqCU/XTqqb1SoybP+i5ZR0t2rLxoYy/sy0eT5oV\nqJzwNb3z17/5/l7zZft3yhb0kzvs7QxkKd/gtL0iON912/6zAra/V8B2rNgAOhcmiXHz5VwyT4lx\n87Vp6061xeOuf4ftFcFeVnPb/l4B2xG+gLwJE9srgr2s5nb7Xtvica1qWq8Va+/Sqqb1WZ3MAGHG\ntDMgb4KzkBXBXqwle1nN7ea9MjUNpEf4AvImOAtVEZxtiLXF49r1Xy+qu0dyerqTgtqram437zUs\n23UChUD4AvImOAt1nXA2ITYwqKXCjDbdvFfbp+EBPxG+gLwLzkJcJ5xNiJkcbWZ6r2zMAaRH+AJf\nsnWDjWxCzKbRJhtzAOkRvgg9Wza+yFU2IWbTaDNM23UCXiN8EWphqLjNJsRsG23aOpsA+I3tJS3A\ntm3uZdtX+W6pGERt8biaf/eiEt2SnB411C0gAIfA35879JM7bC8JyK41UFNm19Toxu/Pk8SBErAV\ndYcINS+3VAQAr3AIQqh5uaUiAHiFaWeEGhW3AGxE+CL0qLgFYBumnQEAMIzwBQDAMKadYZWg70YF\nAG4QvrBGGHajgh04iYPtCF9YIyz3f+XA7y9O4hAEhC+GZDJIwrAbFQd+/4XlJA7hRsEV0uoNksS4\n+XIumafEuPnatHWn2uLxgrxeGHajGurADzPCcBKH8GPki7RMjyBsuyNPLjjw+8/kbRVznRliaQKE\nL9IyHSRh2I3Kpvvp2sZU4Jg6ict1icG2pQlOBPxB+CItP4Ik6LtRhWH0XggmA8fUSVyuM0M2rUnb\ndiJQTAhfpEWQZC8Mo/dCMB04Jk7icp0ZsmlpwqYTgWJD+CItm4Kkd2osUhJVtES6afF8aw8OQR+9\nF4JNgeOVXGeGbFqaCOPnEhSEL4ZkQ5AMnBrrUnimxvqvt5049k9Fo6WKXfSV0K292RQ4Xsl1Zsim\nGaUwfi5BQRfDemG9fKf/pVzHnPH65ydRXTDpViOXdZkWxvsqz66pUdPKW1TW0aLIxy0q62hxNTOU\n688VQhg/l6Bg5AvrhXVqrP9JxdG/vaWJ3/1B0r+Hbe2t1Dmtw//9qLq7z2jsxRepafXywL+3XGeG\nbJhR6m2HZMfSUrEhfGG9sE6N9T+pKCmJpnxO0E8wpH7LBpNu1QVfPjZwtJXp57O9FIbLZ9yz5USg\n2BC+sJ7J6zZNHrD7n1T09HSnfI5XJxh+hlE+FbW5XArD5TMIgrz+tPfv36/KykpNmDBBjzzyiFdt\nApL0XyMrO9KqiqOtnk+Nmd5KU0peb7v40mq99+q2pH/3au3Nj/fWXz7LBrms94e1RgDhktfI98c/\n/rG2bNmiSy+9VNddd52WLl2qUaNGedU2WMTvabzeqbFYrEKS1Nl5Kuffleq9ZDM686ov+q+3jYpI\n0Qu69a/27YpdONLTtTe/r+XMZ9kgl+AOSo2A339T8FfO4dvZ2SlJmjFjhiSptrZWBw4c0IIFVMmF\nTZim8dK9l6i6FBs3+PkDD9he94XNm0F4JZ9lg1yCOwg1AmH6m0Jucg7f119/XRMnTuz77yuuuEKv\nvfbaoPAtLY32jVaQWmnp2WIbW/tp13+9mHLk1Py7F3Xj9+cZbUu+fZXuvRxu26hY1eDnn1dekvRa\nNvXFUPr307CyEp1O8ZyB761Qbvz+PJ0/fJh+s2OPEt1SaVS6p+k21c6ZlfFnVzTW68HN23RhZV3f\nYyffadbdaxvTtj3bn/Hj7y8o36P+bD9O2aK3nzI+r8DtQAh0pxkhJVLXCFkt3Xu5eNRInXx396AD\n9uq1ja5+3ua+WLb0hpRhNPC9FVLtnFmuwjbVz0lKCu7VaxuH/F25/IxpQfwewVs5h+/UqVO1fv36\nvv9ub2/XvHmDz9gSie681ueKgRfrmIXkpKnEldNjvM359lW693LxxReroW5B0vWOa1bcomlTpye9\nVjZ94eeaXv9+mjZ1utYs/yLje7PVtKnTNW3q9KTHMrU7m5/x4+/Ppr8pt2w/TtkiFqtQeXnmaM05\nfGOxmKSzFc9f+9rX9NJLL+knP/lJrr8OFrNpO7x8DfVe3Ky/uu0L29b0uJbTLmH6mzIhjMVpeU07\n/+pXv9Lq1avV1dWlpqYmKp1DKky74OT7Xtz+vN8VxrBbmP6mCs22E1mvRBzHcQr5AmfOJJimyIDp\nHPeC0lcr1t4l55LByzCRj1v060ceKvjrB6WfbEBfueNXP61qWq/EuPmDHi/raNGWjYX/W8qW22ln\ni4rvgfCIRlKf09p0uQsQBH5fKlcoVDvDSkFf42FND/BGEK7bzgXhC+sMtcbjxTWQJoKdNT24FfQT\nzUIL64ks4QvrDFWslG/4mizeoMIYmYS1mMhLYT2RJXxhnUKu8VCFDJvwfXQnjCeyhC+s49UaT6rp\nvLAWbyCY+D4WL8IX1sm0xtP68j79+rfNQ66RpZvOK3VO64JLBr/m+++9p7Z4PNBn1/1PNoaVlWjZ\n0hsG7fIEu4S1mAiZEb7wxVBFJkOt8bS+vC9pn+J0a2TppvNO/uUpnWhvTvq39179rS6+fG6g19oG\nnmyclvTg5m1as/wLq94PxUXJwlpMhMwIXxjnpsgk3RrPU8/sSbpBgJR6jSzddN5FXxmjW29aoP/z\nf38pDR8rp6dbF192lb7yb1XSv1X1/Z6ghUSqk40LK+usWjukuGiwsBYTITPCF8aDJp8ik7R3gxnw\n+FDTebNravT0rr0pd6BK9AQzJLxYO8zne+DmZykuSi2MxUTIjJWFItcbNIlx8+VcMk+JcfO1aetO\ntcXjBXvNfIIimuYbO3CNbEn9Qp1ob0567ET7LjXUnb3f9FA7UA0VErbKd0etfL4Hbn+W4iLgHMK3\nyPkRNPkExbKlN+jku7uTHusfqr1m19SoaeUtKutoUeTjFpV1tCRN5w0VzkEMiVTv5+Q7zYP6JZ18\nvgduf5YtN4FzmHYucn4ETT5FJr03RH982+6Ma2RDTecNtda2o/kF3ytQs50CHvh+zisv0eq1ja6r\nnfP5Hrj9WYqLgHMI3yLnx6UO+RaZ1M6Z5cklNOnC2e+QyHXNuf/7yfYONPl8D9z+LMVFwDncUtAC\nft7SbOCBXkq+ubxtTPVVWzyeFBINdQuM9YcXt1DLtp/y+R6Y/g55XSDILQXdoZ/ccXtLQUa+RS7X\n0UjQLsXJlp8VqH4sBeQzKjU5og1iJTqQCuFrGT9CLdug4QBYWPlMAfd+fyIlUUVLpJsWz3f9meRz\nwmHqZIXLlRAW1BlaxI/LfnIRxEtxgiTTZVLp9P/+dI2t1enRtVZ+f/IRxEp0IBVGvhYJyll9qgPg\n8X+06+j772nF2rtCOQ1tUq7TuEH5/uSDvZARFoSvRYJyVj/wAHj8H+06+uGbmvi9dXLENLQXcpnG\nDcr3Jx9+V6IDXiF8LZLNWb2fBU8DD4BH//aWJn73B0nPCduIKwiKYVTI5UoIC8LXIm7P6v0ueBp4\nAIx2daZ8XphGXLbqfxLW+b8ndOJ//lOX/ceqvn8v9KgwCAWCgI0IX4u4Pav3e21v4AF3zKhYyueF\nacRlo4EnYRdcIn1ycJtO/uUpXTz6qyqNFnZU6PdJoJ/CfqkdCo/wtYybs3o/1/ZSHXBP/M9/6pOD\n23TJv5+bemYdrvBSnYRd8u8/UFlHi3Y8sVFSYTdE8Psk0C/FfNIB7xC+AVSItT23Z/KpDriX/ccq\nnfzLUyrraGEdziC/C6z8fn2/FOtJB7xF+AaQ1xWf2ZzJD3WTerdbH8IbfhdY+f36finWkw54K+R/\nJuGU6XZ52cpm0wxuC2ePXDfjCMvr+4W/AXiBkW9AeVnxmc2ZPNdZ2sPvy278fn2/8DcALxC+yGr6\nsFgPuLZItTbv53R/MV72w98AvED4Iusz+WI84NqAKlt78DeAfBG+4Ew+IKiyhc249jk7hC8kcSYf\nBFTZwlbMymSP8AUColgv7emP0ZWdmJXJXhH92QLBVqyX9vQKyv2uixGzMtlj5AujGLnkrtjX5hld\n2YtZmewRvjDGhnWhoId/Ma/NM7qyF9c+Z4/whTF+j1xsCH/kjtGVvYp9ViYXhC+M8Xvk4nf4Iz+M\nruxWzLMyuSB8YYzfIxe/wx/5YXSFMCF8YYzfIxe/wx/5Y3SFsCB8YYzfIxe/w992QS9Gy1Xry/v0\n1DN79EVXT1G9b/iL8IVRfo5c/A5/m9lcjFbIk4K2eFybH39OF1bWSbLrfSPcIo7jpL45ZQbPPfec\nfvrTn+q9997T66+/rquuuirl886cSaiz81RejQy7WKxCkugnF+grd7Ltp1VN65UYN3/Q42UdLb7e\nNWngSYEknWhvVtPKWzwJR1vft43423MnFqtQeXnmcW3Oq12TJ0/W888/rxkzZuT6KwBYwtZitKEq\n1L1g6/tG+OU87Txx4kQv2wHAR7YWoxU6HG193wi/gq/5lpZG+6YrkFppaVSS6CcX6Ct3su2nFY31\nenDztr61T0k6+U6z7l7b6GtfDysr0ekUj59XXuJJu1Y01uuhzU8rVnlj32M2vG8b8bfnTm8/ZXze\nUP84d+5cHTlyZNDjDzzwgBYtWpRby+CL3orO7h4pWiItW3qDaufM8rtZsETvd+E3O/Yo0S2VRqXV\naxt9/44sW3pDypOC1WsbPfn9tXNmKRqN6Mnte9SVcKx53wi/nAuues2aNUu//OUvKbjKQ6ELGQpd\ntGISRR/umOynQl+i1BaPJ1WoN9Qt8PT3851yh35yx23BlSfTznnmNwosSNsqFuu1pkFl4hIlNtZA\nGOVcVvD8889r/Pjxeu2117RgwQLNnz+4XB92CEpFJ/drDZ5CVyMDYZVz+N5444366KOPdOrUKR05\nckQvvviil+2Ch6KR1DMTtlV0ciAPnqCc2AG2YYerIhCUbRW9PJAzfW1GLpfq8NkAhG9RCMq2il5d\nc2nzVolhk+2JHZ8NcFbe1c6ZUO2cGVWEZ6Wuyt6VdKLgpq/YMtB8tbPbamQbPxv+/tyhn9wxWu0M\neMGrETrrkGZlU43MZwOcRfjCKl5cVsKWgfbiswHO4iuP0FlSv1An2puTHjvRvksNdQt8ahF68dkA\nZzHyRegEpcCsGPHZAGdRcGUBChnco6/coZ/co6/coZ/cKfj9fAEAQG4IXwAADCN8AQAwjIIrH/Vu\nsxcpiSpaIt20eD6FJxZiO0QAXiN8fTJwN6cusc2ejdgOEUAhMO3sE+7gEwx8TgAKgZGvT9hmLxjc\nfk5MTQPIBuHrE7bZCwY3nxNT0wCyxaHeJ2yzFwxuPiempt1ri8e1qmm9Vqy9S6ua1qstHve7SYAv\nGPn6pP82e4qUqDTKNns2crMdYhCWEGyYFmeGADiH8PVR7x182LbNbpnutGT7EoItoTfUDAHhi2Jj\nyeEBCC7blxBsmRYPwgyBbZimDy9GvoDym5a1/U49toRephkCG6bGbWLLjAUKg/BF0fPiIJdpatpP\nbiu2Cx18S+oXJvWzdHaGoGllA0GTAtP04Ub4ouhHHGE/yA0VepK5EdZQMwSrmtaH+jPIhS0zFigM\nwrfIMeII/0Eu07S4yZOPdDMEYf8McmF7IR/yQ/gWubCP+twohoPcUNPiNgRfMXwG2co0Y4FgK+Kv\nNiQ7Drx+s71audCiESfl4yaDr9g/g1Rm19SoaeUtKutoUeTjFpV1tFhVyIf8MPItcrmMOMK2Rmx7\ntXKh2TDCKvbPIB2bC/mQn4jjOKlPez1y5kyCzSMy8HOTjYFrvtK5A2+qP/rUz29W08pbjBwk2JDE\nnWz7qS0eTwq+hroFRXPQ5zvlDv3kTixWofLyzONawtcCfn+psznwrmpar8S4+YMeL+to0ZaNDxW6\nqb73VVDQT+7RV+7QT+64DV+mnZHV1BZrxACQPwqukBUbinMAIOg4ZCIrVKUCQP6YdkZWqEoFgPwR\nvsgalz8AQH6YdgYAwDDCFwAAwwhfAAAMI3wBADCM8AUAwDDCFwAAwwhfAAAM4zpfGBG22xACQD5y\nHvmuX79elZWVuuqqq7Ru3TqdOsWdLpBa720IE+Pmy7lknhLj5mvT1p1qi8f9bhoA+CLn8K2trVV7\ne7sOHjyozz77TNu3b/eyXQiRHc0vJN3/V5JGVtXr2d17fWoRAPgr5/CdO3euSkpKVFJSouuuu05x\nRjFIg9sQAkAyT9Z8H3vsMa1cudKLX4UQikYcJVI8zm0IU2N9HAi/IcN37ty5OnLkyKDHH3jgAS1a\ntEiSdP/992vEiBG6+eabU79AaVSxWIUHTQ2v0tKoJIW2n1Y01uvBzdt0YWVd32Mn32nW3Wsbs37P\nYe+r1pf3afPjz/X1VULS5sef0/nDh6l2zizXvyfs/eQl+sod+smd3n7KJOI4Tuq7o7vw5JNP6rHH\nHtMrr7yi8847L+VzenocJRLdub5EUej9sMLcT60v79NvduxRolsqjUq3LbkhqzDpFfa+unXFj3V6\ndO2gxyuOtmrb1o2uf0/Y+8lL9JU79JM7paVRlZSkXmpLel6uL9DS0qKHH35Y+/fvTxu80tkPqrOT\nSuih9J5Jhrmfpk2drmlTpyc9lsv7DXtffdGVeiH89JmerN5z2PvJS/SVO/STO7FYhcrLM0drzqtu\na9eu1aeffqo5c+ZoypQp+tGPfpTrrwLwpWgk9UQU6+NAuOQ88j18+LCX7QAgaUn9Qm3aujPp0qwT\n7bvUtLLBx1YB8Bo7XAEW6a1qfnb3XiV6zo54m1Y2UO0MhAzhC1hmdk0NYQuEHCtJAAAYRvgCAGAY\n4QsAgGGELwAAhhG+AAAYRvgCAGAYlxoBCATu9oQwIXwBWK8tHk/a+SshadPWnZJEACOQmHYGYL0d\nzS8kbbkpSSOr6vXs7r0+tQjID+ELwHrdTupbtCVS3wQKsB7hC8B63O0JYcNXF4D1ltQv1In25qTH\nTrTvUkPdAp9aBOSHgisA1uNuTwgbwhdAIHC3J4QJ084AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG\n+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACG\nEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCA\nYTmH73333acrr7xS1dXVamxs1PHjx71sFwAAoZVz+N511116++239dZbb2nChAnauHGjl+0CACC0\ncg7fESPpBX3rAAAEu0lEQVRGSJISiYQ+++wznXfeeZ41CgCAMIs4juPk+sP33nuvtmzZossvv1z7\n9u1TeXn5oOf09DhKJLrzamTYlZZGJYl+coG+cod+co++cod+cqe0NKqSkkjG5w0ZvnPnztWRI0cG\nPf7AAw9o0aJFkqTPP/9c9957ryRpw4YNubYXAICikdfIt9ehQ4d0xx136LXXXvOiTQAAhFrOa76H\nDx+WdHbN95lnnlFdXZ1njQIAIMxyDt977rlHkydP1jXXXKNEIqE77rjDy3YBABBaOYfvrl27dOjQ\nIf3pT3/SQw89pIsuuijl87ge2L3169ersrJSV111ldatW6dTp0753SQrPffcc6qqqlI0GtUbb7zh\nd3Oss3//flVWVmrChAl65JFH/G6OtZYvX64xY8Zo8uTJfjfFah999JFmzZqlqqoqzZw5U9u3b/e7\nSVY6ffq0pk2bpurqak2fPj1zDZRTYP/617/6/v/PfvYz57777iv0SwZWa2ur093d7XR3dzsrV650\ntm7d6neTrPTuu+8677//vjNz5kznz3/+s9/NsU51dbUTj8edDz/80Ln88sudo0eP+t0kK+3fv995\n4403nEmTJvndFKt1dHQ4b775puM4jnP06FHn61//etJxHed89tlnjuM4zunTp52qqirn8OHDaZ9b\n8O0luR7Yvblz56qkpEQlJSW67rrrFI/H/W6SlSZOnKhvfvObfjfDSp2dnZKkGTNm6NJLL1Vtba0O\nHDjgc6vsdO2116adscM5Y8eOVXV1tSRp1KhRqqqq0sGDB31ulZ2GDx8uSfr000+VSCQ0bNiwtM81\nsrfzvffeq7Fjx+rVV1/VnXfeaeIlA++xxx7ru5wLcOv111/XxIkT+/77iiuu4CoEeOaDDz5Qe3u7\nrr76ar+bYqWenh5deeWVGjNmjNasWaPx48enfa4n4Tt37lxNnjx50P9+//vfS5J+/vOf6+9//7uu\nvvpq3X333V68ZGBl6itJuv/++zVixAjdfPPNPrbUX276CYA5n3zyiRoaGrRhwwadf/75fjfHSiUl\nJXr77bf1wQcf6NFHH9Wbb76Z9rmlXrzgSy+9lPE5w4cP1/Lly4u+KjpTXz355JP6wx/+oFdeecVQ\ni+zk5juFwaZOnar169f3/Xd7e7vmzZvnY4sQBl1dXaqvr1djY6MWL17sd3Osd9lll+n666/XgQMH\nNGXKlJTPKfi0M9cDu9fS0qKHH35Yv/vd71gbd8nJf4+YUInFYpLOVjx/+OGHeumllzRt2jSfW4Ug\ncxxHK1as0KRJk7Ru3Tq/m2OtY8eO6eTJk5Kk48ePq7W1dcgTFU92uBrKTTfdpPfff18VFRWaOXOm\n7rnnHooc0pgwYYLOnDmjkSNHSpK+853v6NFHH/W5VfZ5/vnn1dTUpGPHjikWi2nKlCl68cUX/W6W\nNeLxuH74wx+qq6tLTU1Nampq8rtJVlq6dKni8biOHz+u0aNH6/7779ftt9/ud7Os8+qrr2rGjBn6\n1re+pUjk7J7Fv/jFL5hRGeDQoUNatmyZuru7NXbsWN1666267bbb0j6/4OELAACSGal2BgAA5xC+\nAAAYRvgCAGAY4QsAgGGELwAAhhG+AAAY9v8B4oY6ZvGrkZUAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# And we multiple that data by a matrix A\n",
      "A = array([[1, .5], [.5, 1]])\n",
      "X2 = dot(X1, A)\n",
      "plot(X2[:, 0], X2[:, 1], \"o\");"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAFRCAYAAAA1uqfwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+Q1NW55/HPdM+AgOwgIkzwGk1ViIwDCuYiJDcywMII\nKiKQCJRBSkCxvDKh9oquZdxU2ERLXYuIlndZiGUUkYuglNFyAgppgnUlEvxBjeCPPzDqHRIYAhEE\noWe++weZkR76x+nu8/39flVZJU3TfeaZ78zT55znPN8Kx3EcAQAAzyT8HgAAAHFD8gUAwGMkXwAA\nPEbyBQDAYyRfAAA8RvIFAMBjZSXf48ePa+TIkRo2bJhGjRqlpUuX2hoXAACRVVHuOd8vv/xSPXv2\n1FdffaXvfve72rBhg7797W/bGh8AAJFT9rJzz549JUlHjhxROp1W9+7dyx4UAABRVnbybW9v12WX\nXaYBAwbojjvu0AUXXGBjXAAARFbZy84d9u7dq6uvvlrPPvushg8f3vl4e7ujdLrNxltEVmVlUpKI\nkwFiZYY4mSNWZoiTmcrKpBKJisLPs/WGF110ka6++mpt3749I/mm0206fPiYrbeJpOrqHpJEnAwQ\nKzPEyRyxMkOczFRX91C3boVTa1nLzgcOHNChQ4ckSa2trdq4caOmTJlSzksCABB5Zc18W1paNGfO\nHLW1tammpkZ33nmnvvGNb9gaGwAAkVRW8h06dKh27txpaywAAMQCHa4AAPAYyRcAAI+RfAEA8BjJ\nFwAAj5F8AQDwGMkXAACPkXwBAPAYyRcAAI+RfAEA8BjJFwAAj5F8AQDwGMkXAACPkXwBAPAYyRcA\nAI+RfAEA8BjJFwAAj5F8AQDwGMkXAACPkXwBAPAYyRcAAI+RfAEA8BjJFwAAj5F8AQDwGMkXAACP\nkXwBAPAYyRcAAI+RfAEA8BjJFwAAj5F8AQDwGMkXAACPkXwBAPAYyRcAAI+RfAEA8BjJFwAAj5F8\nAQDwGMkXAACPVfo9AAAAyrU5ldKa9S+rzalQssLRzOnXalx9vd/DyonkCwAItc2plJatXKu+ddMl\nSWlJy1aulaTAJmCWnQEAobZm/cudibdD37rp+o8XXvFpRIWRfAEAodbmVGR9PN3u8UCKwLIzAKAs\nfu+3JiscpbM8Xhng6WWAhwYACLqO/db0wElyzp+o9MBJWrZyrTanUp6NYeb0a3WweX3GYweb12nG\ntGs8G0Oxykq+n376qcaOHau6ujqNGTNGq1evtjUuAEAIBGG/dVx9vRrn36CqliZVfN6kqpYmNc6f\nEdhiK6nMZeeqqiotXbpUw4YN04EDB3TFFVdo8uTJ6t27t63xAbDI7+VBRE9Q9lvH1deH6louK/nW\n1NSopqZGktSvXz/V1dVpx44dGjt2rJXBAbAnjMcxEHxh3G8NggrHcRwbL/Txxx+roaFBu3btUq9e\nvTofb293lE632XiLyKqsTEoScTJArMxki9ON836i4/0bznhuj/0btWrlo56NLWi4pszkitPG17bo\nwcdXqU/ttM7HDr2/XncvnK2G8fGbiFVWJpVIZF8NyHiejTf74osvNGPGDC1dujQj8QIIjrYcy4Dk\nHJSjI8E+vWaD0m1SZVJaENPEW4yyk+/Jkyc1ffp0zZ49W1OmTDnj79PpNh0+fKzct4m06uoekkSc\nDBArM9ni5LTnyLJOe6zjyTVlJl+cRo4YpZEjRmU8Ftd4Vlf3ULduhVNrWavyjuNo3rx5GjJkiBYt\nWlTOSwFwWRiPYwBRVdbM94033tCqVat06aWXavjw4ZKkBx54QBMnTrQyOAD2dBRV/ccLryjdfqog\nJujHMYCoslZwlcuJE+nYLj+YYtnLHLEyQ5zMESszxMmMJ8vOAACgeCRfAAA8RvIFAMBj3NUIAIpA\ni07YQPIFAENBaNFJ8o8Glp0BwJDfd/AJwu37YAfJFwAM+X0HH7+TP+wh+QKAoWRF9rYIXt3Bx+/k\nD3tIvgBgyO8WnX4nf9jDtwwADI2rr1fj/BtU1dKkis+bVNXS5GmLTr+TP+yh2hkAijCuvt636mL6\nc0cHyRcAQsTP5A97SL4AkAfnauEGki8A5BCEphqIJgquACAHztXCLSRfAMiBc7VwC8vOiD329JBL\nssJROsvjnKtFubiEEGv0ykU+nKuFW5j5Itby7ekx+wXnauEWki9ijT09FMK5WriB5ItYY08PfqLe\nIL74FYNYY08PfqHeIN6Y+SLW2NODX6g3iDeSL2KPPT34gXqDeCP5ArCKfUwz1BvEG99mANawj2mO\neoN4Y+YLwBr2Mc1RbxBvJF8A1tjcx4zD8jX1BvFF8gVgja19TG7lh6hjzxeANbb2MbmVH6KOmS8A\na2ztY3IMB1FH8gVglY19TI7hIOq4lAEEDsdwEHXMfIGIC2PVMMdwEHUkXyDCwlw1zDEcRBnLzkCE\nUTUMBBMzXyDCqBpGOU7fsuheldCcWddr5IhRfg8rEki+QITZrhoO4/4xStN1y+K4pAcfX6U75n7l\nyfc86tcay85AhNmsGuamCfGSbcuiT+00T7Ys4nCtkXyBCBtXX6/G+TeoqqVJFZ83qaqlqeSqYfaP\n48XPLYs4XGssOwMRZ6tqmP3jePGz0UkcrjVmvgCMJCucrI/TdSqasm1ZHHp/vSeNTuJwrZU18507\nd65eeeUV9e/fX7t27bI1JgABcXrRy8EDf9WX/7VK5//zjzv//mDzOjXOn+HjCOGWro1OzuqW0IKF\nsz2pdp45/dqMYi8petdaheM42T9iGPjDH/6gs88+WzfddFPO5HviRFqHDx8reYBxUF3dQ5KIkwFi\nZcZGnLpWu0rS3jf+n/r2OVvVffqqMiHNmHZN6CtQuabMeB2nzalURoezsFxr1dU91K1b4XltWTPf\nK6+8Unv37i3nJQAEVLail4v+5VZVtTRp+aMP+TQqxEXUO5y5XnBVWZns/MSE7Cork5JEnAwQKzM2\n4lSRSOb6i0jFn2vKDHEy0xGngs9zeRwAQmDja1v0m+c2qK1dSiakObOuVzIhnczyXMPfLa6OrWH8\nWG8HAVjmevJNp9vYSymAPSdzxMpMMXHqurd7UtIDy57RmJG1+v329VmLXrzc98s2tqNf2uuyFMRr\nKojdnYIYpyAy3fONUOE2gFLkamiw+6M/W2vQYXtsUWq20FUcujuhzJnvrFmzlEql1NraqgsuuEBL\nlizRzTffbGtsADyQr6GB30UvcWi20FW+Dxx+z35hT1nJ97nnnrM1DgA+caOTka1l02SFo7981qz9\nn7yjRCKp9vY2nXfhMNV4vO/spTh+4IgjCq6AmLPd0KDrPm1a0rKVayWp6AR8yXcu1LsvbVXduNs6\nH2ve/H819rrRgdwXtcHPto7wDskXiLmunYwqEyprb9fmsun7H36SkXglqW7cbdq2fbV+v323lQTv\nplI+IMShuxNIvgBkd2/X5rJprtf6y/6/adCYf814LGj7oqWuANj+MIRgIvkCsMrmsmmu12pvy/Zo\nsPZFy1kB8LvQDe5jFwGAVR13w2n9rFl73nhWH/7nGu363SOqHfTNkl/rdAeb16l/vz5Znx+kfVEK\np5APM18AVo2rr9d7u3Zp7ct/0OAxCzof//329bp0aKqoGV2uJVhJgd8XpXAK+ZB8AVj3/oefZCRe\nqfQ92XxLsDb2RTvaV351st1q1XShwqmoVmvDDMkXgHVeLLna2BfdnErp8SefV5/aaZLsVk3nK5yy\neRwL4UTyBWBdWJZc16x/uTPxdrBZNZ3rAwJdrBCwHwUAUZCrUGrGtGt8GlF2fhVFUYwFZr4ArAvL\nWVW/ZuhhWRmAe0i+AFwRhrOqM6dfm7HnK3lTNU0XK5B8AcTWuPp69erZXU+v2aDjJ9o9m6GHZWUA\n7qlwHMdx8w1OnEhz8+UCuEm1OWJlxo04RfVoDNeUGeJkprq6h7p1KzyvZeYLoCCOxgB2sb0PoKB8\nR2MAFI+ZLxBzJsvJHI0pT1SX7FE6ki8QIF7/kjZdTg7z0Ri/Ex9L9sgmBD86QDx0/JJOD5wk5/yJ\nSg+cpGUr12pzKuXae5ouJ4elaUZXfsS0K5bskQ3JFwgIP35Jmy4nj6uvV+P8G1TV0qSKz5tU1dIU\niqMxQUh8LNkjG5adgYAo5pe0raXUYpaT3Wya4dbScBASX5iX7OEevv1AQCQrsh+57/pL2uZSahCW\nk91cGjaNqZuCEGMEDzNfICBMWw7avCNOEDotFfP1FDtDDkIbxyDEGMFD8gUCwvSXtO2lVL97MJt+\nPaVUDQcl8fkdYwQPyRcIEJNf0lHbQzT9ekqd8ZP4EEQh/XEF4itqe4imX08QiqcAW5j5AiETlKVU\nW0y/nqjN+BFvJF/AQ7aO1IR5KTVXDAp9PUEongobv7t7ITeSL+AR2gyWF4OozfjdxvUWbNzPNwC4\nT6a5MMfq1sbFSg+cdMbjVS1NWv7oQ1bfK6hx8jIGpoIaq3LZjnVU42Sb6f182S0BPELBEDHwErEO\nNpIv4JEgdFvyGzHwDrEONr4NgEeidkSoFMTAO8Q62Ci4AjxCwRAx8BKxDjYKrgKAQgZzxMoMcTJH\nrMwQJzMUXAEAEFAkXwAAPMaeL1AkugYBKBfJFygCXYMA2EDyBYpg80b2fmDWDgQDyRcoQpi7BjFr\nB4KD5AsUweZt7byehYZ91g5ESVnVzlu3blVtba0GDRqkxx57zNaYgMCy1TWoYxaaHjhJzvkTlR44\nSctWrtXmVMrmcDOEedZeqs2plG5tXKx5C+/SrY2LXY0vUIyyZr4/+clPtHz5cl144YW66qqrNGvW\nLPXr18/W2IDAsdU1yI9ZaNxuRs8yO4Ks5OR7+PBhSdLo0aMlSQ0NDdq+fbuuuYa+oYg2GzeyzzcL\n7ViO/stf/6oDrX9TTU2Nzj2nuuxl6bjdjJ5ldgRZycn3rbfe0uDBgzv/fMkll+jNN988I/lWViY7\n25Ihu8rKpCQRJwNRiVX3qoSOZ3n8yBeH9PiTz6ut98Xaf6RVg//7IkmnZm2PP/m8evXsrobxYwu+\nfrY4Tb1uonr17K6n12xQuk2qTEr3NN5k9Hr5bHxti37z3Aa1tUvJhDRn1vVlv6YNFYlkrr/IiEtU\nrim3ESczHXEq+DyXx4EQC+ov1SiYM+t6Pfj4KvWpndb52KH31yvRflJ9hszUnjee1eAf/Djj3/Sp\nnaan12wo63vQMH6s1e/hxte2ZHwdJyU9+PiqzvfK9W+8uK6SiVPj6crwdyPgqpKT74gRI7R48eLO\nPzc3N2vixIlnPC+dbqMRdwFBbFjedb/spKQHlj2jo19+5euSXRBjVYqRI0bpjrlfZewd3zHvBj27\n7hU5khI5Zm3HT7Qbfe1exenXz6zP+AAhnfqQ8OSqFzRyxKgznu/ldfXDKZNyLrOfHpeoXFNuI05m\nTG+sUHLyra6ulnSq4vmb3/ymNm3apJ/97GelvhwChv0y92XbO16z/mWlJbW3t2X9N0Erjiq2gtrL\n64pb6iHIylp2/tWvfqUFCxbo5MmTamxspNI5QuJ4LCUIOoqizrtwmPZsW5Wx9BzE4qhiK6i9vq5s\nFMcBbigr+dbX12v37t22xoIAiduxlKA4fbams9v04eZHNWDAAPXr28fTWZtpA5BiK6i5roBTKLhC\nVnE7lhIkfs/WijkfW+zSLtcVcEqF4ziOm29w4kSaDfoCglrIsDmVyvilOmPaNb4v4QU1VkFTTpxu\nbVys9MBJZzxe1dKk5Y8+VPbYgnZdcU2ZIU5mXC+4QvT5PQODP9zel+W6Asrs7QwgepIV2RfD2JcF\n7GHmi1Dj/rS5bXxti379zPqiY8O+LOA+ki9Ci8b5uXXtPFVMbDgfC7iPgqsAoJDB3OmxcrswKJd8\ns203ZuKlvObt/+N/6nj/hjMedzs2YcTPnxniZIaCK0SeW4VBhZJrrtm2/vH/Nmfipc7u23LEgCYp\nQDCQfBFabjRsKJTs8rVHdBzHeuvEUtsx5rypAEVTQCDwo4jQmjn9Wh1sXp/x2MHmdZoxrfR7SudL\ndlL+2bYbM/FSX3POrOt1aPcLGY+VGxsA9jDzRWi5URhUKNnlm207jv2ZeKmz+45b9D256gWKpoAA\nIvki9E7VDFbIRu1goWRX6BiO7SM65Rz7aRg/Nutt/QD4j+SL0Cq1GClfQVWhZGcy27Y5E7c5u+dM\nNBAcHDUKAEr4zZV71Khrwpakg83r1Tj/hoyK5iD1Hi5F12vK5OuOK37+zBAnMxw1QuSVUoxkUj0c\nxd7DbtzEvpyZNLNwxB3JF6FVSjGS1zdzDwrbX3c53cXoTAZw1AghVspRoyDfNGBzKqVbGxdr3sK7\ndGvjYm1Opay9tu2vu9CRLLf+LRAVzHxhhR/LiKUUI3l90wDTuLg9G7T9dZczk47r6gNwOpIvyubn\nMmKx+7Ne3jSgmLi4sSd7OttfdzndxdzoTAaEDckXZXM7cdjmVUFVMXHxYjZo8+suZybNLQsBki8s\nYBkxu2LiErbZYDkzaW5ZCJB8YUHYEodXiolLGGeD5cyko3icCygGyRdlC2LiCMI50mLiwmwQiBc6\nXAVAFDrHeNUVyiRWQermdHpc/tb6FyUSlao+51zXPxCUe00F4cOLV6Lw8+cF4mTGtMMVyTcAuKjN\nmcQqV9vJvzevVt9zzvEloXj9gaCcaypIH168wM+fGeJkxjT5xnxXDlHU+rfDWR//5PP9Sg+cJOf8\niUoPnKRlK9dabWSRT5gaS4RprEBYkXwROfv27cv6eLde/TL+7GVCCVNFeJjGCoQVBVcInY2vbdGv\nn1mfc/m437nnaM+2VRr8gx93Pvbea/+u8wePPuO1vEooYaoID9NYgbDixwmhsvG1LXrw8VV5l48H\n9O+v8y4arg/eeFYf/ucaffDGs2pva9O5/1R3xut5lVBK6UPtlzCNFQgrZr4w4lb1a7Gv+5vnNqhP\n7bSMx7p2jeo44nPxv9zY+Zy9byzX5ztW6fx//no27OVxqDAdJQrTWIGwIvmiILd6N2d73V/8n3/X\nEyueznkcpy3HMvHpy8fZksdP77z9jMe8TihhaiwRlrHG6UgUooXki4Lc6t3c9XVbP2vWcaen/tuQ\nG+Uoe5JPJqSTWV6r6/JxruTh1i9mkoD3uC8wwow9XxTkVvVr19fd/8k7GUVS0pkVyXNmXa9Du1/I\neI7f+5EdScCvY0xxxZEohBkzXxTkVvVr19dNJJJZn3d6km8YP1aS9OSqFwKzHxm2uzpFBUeiEGYk\nXxTkVu/mrq/b3t6W9Xldk3zD+LEaOWJUWe9tE0nAHxyJQpiRfFGQW9WvXV+3pnfa14rkUpEE/BHE\nG3oApujtHABR6plabuFRoRs0BDFW2XshrzP6gOJWoVYQ4+QGGzf0iEusykWczHBjhRCJykXtRUP+\noMaqlCRQTLyKTdJBjVMQESszxMmMafJl2RnWxLnwqJRzsabxMjlS0zU5z5s9vbM4DUDwkHxhDYVH\nxTGNV6EknS05P/j4KkkKVGEagK9REgJrkhXZdzAoPMrONF6FknS25NyndpqeXrOh7DECcAe/FmEN\nDfmLYxqvQkk6Z3LOfnILQACUnHyff/551dXVKZlMaufOnTbHhJAaV1+vxvk3qKqlSRWfN6mqpcn3\nBhhBZhqvQkk6Z3LO3rMEQACUXO28Z88eJRIJLViwQI888oguv/zyrM+j2rkwqgjNxTVW+aqps1VN\nH3p/ve5eOJs9XwNxvaaKRZzMuF7tPHjw4FL/KYAi5aumztYE5e6Fs9Uwfiy/KIGAcr3aubIy2fmJ\nCdlV/mN9kDgVRqyym3rdRE29bmLnn4mTOWJlhjiZqTTc78mbfCdMmKB9+/ad8fj999+vyZMnlzYy\nAABiLm/y3bRpU9lvkE63sfRVAHsp5oiVGeJkjliZIU5mTPd8rRw1crlDJQAAkVLynu+LL76oxsZG\nHThwQNdcc42GDx+uV1991ebYUCK3mvUDAOwoOflOnTpVU6dOtTkWWGDSB7jjeSRoAPAHHa4iJl8f\n4A4dCTo9cJKc8ycqPXCSlq1cq82plNfDBYBY4sYKEWPSrP/0BN36WbP2f/KOEomk7vvFI5LEDBgA\nXEbyjZhkhaN0lsdPb9bfkaBbP2vW/r1va/APftz5d9mWqOE9tgWAaGPZOWJMmvV39ALe/8k7GYlX\nOnOJGt5jWwCIPma+EZOt1WDXZv0zp1+rZSvXKpHI3onF7fvvMqvLL9e+/f/6xSN6dt0rxAyIAJJv\nBOXrA9zx95I693i7cvP+u6bV2HGWa99ePWtOzYRFzICwY9k5psbV1+t///TfPL//rkk1dtzlukWg\n0/71DXqJGRBuzHwjrNDyrskStW0m1dhx17EtcPqHlD3bntF5F2XetpOYAeFF8o0o0+XdQkvUtplU\nY8dd1w9FH+zZo/MunqBz/6ku43nEDAgvfnwjKqjLuybV2DiVgJc/+pB+/dhDWvLTf1PF4T0Zf0/M\ngHBj5htRQV3e9WOpO+yIGRA9JN+ICvLyrtdL3VFAzIBoIfmGSDHnY7MV7RxsXqfG+TO8Gi4AIAeS\nb0gUez6WpUoACC6Sb0jkK6DKlVDzLVXSZQoA/EPytahQQisn4dksoKLLVHZ8IAHgFZKvJYUSWr6/\nn3rdxIKvb7OAqpRZdNTxgQSAlwJQ+xoNhc7Vlnvu1ub52KAeQ/JTru/P48uf1K2NizVv4V26tXEx\ndxYCYAUzX0sKJbRyE57NAqogH0PyS7bvT+tnzdp/8Lj6XDpJErNhAPaQfC0plNBsJDxbZz05hnSm\nbN+f/Z+8o8FjFmQ8FvfleQB2xHiuY1ehZeEgtVUcV1+vxvk3qKqlSRWfN6mqpSn2x5CyfX/ajh3M\n+tw4L88DsIOZryWFloWDdu6WjkmZsn1/Lqjpm/W5cV6eB2BHheM42W8easmJE2kdPnzMzbcIverq\nHpIUijj5fRzHy1h1rYCWvl6eD/oHlzBdU34jVmaIk5nq6h7q1q3wvJaZL4zF7ThO0FYrAEQHyRfG\n4ng+mOV5AG4g+cZIuUvGnA8GADtIvjFhY8m4mONS+RK93/vGAOA3km9M2FgyNj0fnC/R6x//H5d9\nYwDIhuQbEzaWjE0LkPIlesdxYrdvDABdkXxjwlZLSZMCpPyJnn1jAKBdQEx42WErWZH96HhlIv/f\nAUBc8CsvJrxsKZkv0QepzSYA+IVl5xjx6syqyd4wjSsAxBntJQOAtm3miJUZ4mSOWJkhTmZM20uy\n7AwAgMdIvgAAeIzkCwCAxyi4soi2iQAAEyRfS6J8uz0+VACAXSw7W5KvpWKYdXyoSA+cJOf8iUoP\nnKRlK9dqcyrl99AAILSY+VoS1dvtxfEevjawWgAgH5KvJbZ6JwdNVD9UuCnKWxAA7Cg5NSxevFi1\ntbW6/PLLtWjRIh07Fu+D11Ftm0gv5uJFdQsCgD0l/wptaGhQc3OzduzYoaNHj2r16tU2xxU6XvZO\n9lJUP1S4idUCAIWUvOw8YcKEzv+/6qqr9NJLL2nevHlWBhVWXvVO9pLpPXzxtahuQQCwx8qe74oV\nKzR//nwbL4UAiuKHCjfNnH5txp6vdGq1oHH+DB9HBSBI8t5YYcKECdq3b98Zj99///2aPHmyJGnJ\nkiV67733tG7duqyv0d7uKJ1uszTcaKqsTEoScTIQllhtfG2Lnl6zQek2qTIp3TTzejWMH+vZ+4cl\nTkFArMwQJzOVlUklEtm3nk5X1l2NnnrqKa1YsUKvv/66zjrrrKzPIfkWxkVtjliZIU7miJUZ4mTG\nNPmWvOzc1NSkhx9+WFu3bs2ZeKVT3yhuQZUft+oyR6zMECdzxMoMcTLj+i0FFy5cqCNHjmj8+PEa\nPny4br/99lJfCgCAWCl55vvRRx/ZHAcijo5PAPA1OlzBdXR8AoBMnDyE6+j4BACZSL5wHR2fACAT\nyReuoz80AGTi1x9cR39oAMhEwRVcR39oAMhE8oUn6A8NAF9j2RkAAI+RfAEA8BjJFwAAj7HnWyba\nJgIAikXyzcEkqdI2EQBQitAlXy9mmqZJNV/bRJIvACCXUO35diTF9MBJcs6fqPTASVq2cq02p1JW\n38e0FzFtEwEApQhV8vWqQb9pUqVtIgCgFKFKE17NNE2TKm0TAQClCFXy9WqmaZpUx9XXq3H+Dapq\naVLF502qammibSIAoKBQFVzNnH5tRiGUdCopNs6fYfV9iulFTNtEAECxQpV8vWzQT1IFALglVMlX\nIikCAMIvVHu+AABEQehmvkFFm0kAgCmSrwW0mQQAFINlZwu8av4BAIgGkq8FtJkEABSD5GsBbSYB\nAMUgPVhAm0kAQDEouLLAy+YfAIDwI/laQvMPAIAplp0BAPAYyRcAAI+RfAEA8BjJFwAAj5F8AQDw\nGMkXAACPkXwBAPAYyRcAAI+RfAEA8BjJFwAAj5F8AQDwGMkXAACPkXwBAPAYyRcAAI+VnHzvu+8+\nXXbZZRo2bJhmz56t1tZWm+MCACCySk6+d911l95991298847GjRokB599FGb4wIAILJKTr69e/eW\nJKXTaR09elRnnXWWtUEBABBlFY7jOKX+43vvvVfLly/XxRdfrC1btqhbt25nPKe93VE63VbWIKOu\nsjIpScTJALEyQ5zMESszxMlMZWVSiURFweflTb4TJkzQvn37znj8/vvv1+TJkyVJX375pe69915J\n0tKlS0sdLwAAsVHWzLfDrl27dMstt+jNN9+0MSYAACKt5D3fjz76SNKpPd/nnntO06ZNszYoAACi\nrOTke88992jo0KH6/ve/r3Q6rVtuucXmuAAAiKySk++6deu0a9cu/fGPf9RDDz2kc845J+vzOA9s\nbvHixaqtrdXll1+uRYsW6dixY34PKZCef/551dXVKZlMaufOnX4PJ3C2bt2q2tpaDRo0SI899pjf\nwwmsuXPnasCAARo6dKjfQwm0Tz/9VGPHjlVdXZ3GjBmj1atX+z2kQDp+/LhGjhypYcOGadSoUYVr\noByX/f3vf+/8/5///OfOfffd5/ZbhtbGjRudtrY2p62tzZk/f76zcuVKv4cUSLt373Y++OADZ8yY\nMc6f/vRevIDoAAADI0lEQVQnv4cTOMOGDXNSqZSzd+9e5+KLL3b279/v95ACaevWrc7OnTudIUOG\n+D2UQGtpaXHefvttx3EcZ//+/c63vvWtjN/r+NrRo0cdx3Gc48ePO3V1dc5HH32U87mut5fkPLC5\nCRMmKJFIKJFI6KqrrlIqlfJ7SIE0ePBgfec73/F7GIF0+PBhSdLo0aN14YUXqqGhQdu3b/d5VMF0\n5ZVX5lyxw9dqamo0bNgwSVK/fv1UV1enHTt2+DyqYOrZs6ck6ciRI0qn0+revXvO53rS2/nee+9V\nTU2Ntm3bpjvvvNOLtwy9FStWdB7nAky99dZbGjx4cOefL7nkEk4hwJqPP/5Yzc3NuuKKK/weSiC1\nt7frsssu04ABA3THHXfoggsuyPlcK8l3woQJGjp06Bn//fa3v5Uk/fKXv9Sf//xnXXHFFbr77rtt\nvGVoFYqVJC1ZskS9e/fWj370Ix9H6i+TOAHwzhdffKEZM2Zo6dKl6tWrl9/DCaREIqF3331XH3/8\nsZ544gm9/fbbOZ9baeMNN23aVPA5PXv21Ny5c2NfFV0oVk899ZR+97vf6fXXX/doRMFkck3hTCNG\njNDixYs7/9zc3KyJEyf6OCJEwcmTJzV9+nTNnj1bU6ZM8Xs4gXfRRRfp6quv1vbt2zV8+PCsz3F9\n2ZnzwOaampr08MMP66WXXmJv3JBTfo+YSKmurpZ0quJ579692rRpk0aOHOnzqBBmjuNo3rx5GjJk\niBYtWuT3cALrwIEDOnTokCSptbVVGzduzPtBxUqHq3x++MMf6oMPPlCPHj00ZswY3XPPPRQ55DBo\n0CCdOHFCffv2lSR973vf0xNPPOHzqILnxRdfVGNjow4cOKDq6moNHz5cr776qt/DCoxUKqXbbrtN\nJ0+eVGNjoxobG/0eUiDNmjVLqVRKra2t6t+/v5YsWaKbb77Z72EFzrZt2zR69Ghdeumlqqg41bP4\ngQceYEWli127dmnOnDlqa2tTTU2NbrzxRt100005n+968gUAAJk8qXYGAABfI/kCAOAxki8AAB4j\n+QIA4DGSLwAAHiP5AgDgsf8PPrwE33ZWMz8AAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Multiplying by the inverse of matrix A will undue the\n",
      "# linear transformation of X that was obtained by multiplying by A\n",
      "Ainv = inv(A)\n",
      "X1b = dot(X2, Ainv)\n",
      "plot(X1b[:, 0], X1b[:, 1], \"o\");"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAFRCAYAAAA1uqfwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9wVNX9//HXZpNgcHCVIlD8UO1MqcRADfaDUD+VAIUI\nAkUTNTA2MgJCpwMpf4iO4/ht63xrR50OBRln+Ej90SIiEuTTyphGDZ9l7IxU6o8y8Udx5mNrndAv\nPz6k/gDJJvf7ByZkk93s3d2755579/mY6UxdNtmzZzf3dc8573tuxHEcRwAAwJgSvxsAAECxIXwB\nADCM8AUAwDDCFwAAwwhfAAAMI3wBADAsr/A9ffq0pk2bpurqak2fPl0bNmzwql0AAIRWJN/rfD//\n/HMNHz5cX3zxhb797W9rz549+sY3vuFV+wAACJ28p52HDx8uSfr000+VSCQ0bNiwvBsFAECY5R2+\nPT09uvLKKzVmzBitWbNG48eP96JdAACEVt7Tzr0+/PBDXX/99Xr66ac1ZcqUvsd7ehwlEt1evERo\nlZZGJYl+coG+cod+co++cod+cqe0NKqSkkjm53n1gpdddpmuv/56HThwICl8E4ludXae8uplQikW\nq5Ak+skF+sod+sk9+sod+smdWKxC5eWZozWvaedjx47p5MmTkqTjx4+rtbVVixcvzudXAgAQenmN\nfDs6OrRs2TJ1d3dr7NixuvPOO/XVr37Vq7YBABBKeYXv5MmT9cYbb3jVFgAAigI7XAEAYBjhCwCA\nYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABjm2f18\ngUza4nHtaH5B3U5E0YijJfULNbumxu9mAYBxhC+MaIvHtWnrTo2sqpckJSRt2rpTkghgAEWHaWcY\nsaP5hb7g7TWyql7P7t7rU4sAwD+EL4zodiIpH0/0GG4IAFiA8IUR0YiT8vFSvoEAihBrvhYJc0HS\nkvqFSWu+knSifZeaVjb42CoA8Afha4nWl/eFuiCp9z08u3uvEj1nR7xNKxtC8d4AIFsRx3FSzwd6\n5MyZhDo7TxXyJQIvFqvQrSt+rNOjawf9W1lHi7ZsfMiHVtkpFquQJL5TGdBP7tFX7tBP7sRiFSov\nzzyuZcXNEt1pCo8oSAKA8CF8LRFN80lQkAQA4cOh3RLLlt6gE+3NSY+daN+lhroFPrUIAFAoFFxZ\nonbOLH32+RcUJAFAESB8LTK7poawBYAiwLQzAACGEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG\n+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACG\n5RW+H330kWbNmqWqqirNnDlT27dv96pdAACEVmk+P1xWVqYNGzaourpax44d09VXX61FixZpxIgR\nXrUPAIDQyWvkO3bsWFVXV0uSRo0apaqqKh08eNCThgEAEFYRx3EcL37RBx98oNraWh06dEjnn39+\n3+M9PY4SiW4vXiK0SkujkkQ/uUBfuUM/uUdfuUM/uVNaGlVJSSTj8zwpuPrkk0/U0NCgDRs2JAUv\nAAAYLK81X0nq6upSfX29GhsbtXjx4kH/nkh0q7PzVL4vE2qxWIUk0U8u0Ffu0E/u0Vfu0E/uxGIV\nKi/PHK15jXwdx9GKFSs0adIkrVu3Lp9fBQBA0cgrfP/4xz9q27Ztamtr05QpUzRlyhS1tLR41TYA\nAEIpr2nn7373u+rp6fGqLQAAFAV2uAIAwDDCFwAAwwhfAAAMI3wBADCM8AUAwDDCFwAAwwhfAAAM\nI3wBADCM8AUAwDDCFwAAwwhfAAAMI3wBADCM8AUAwLC87moEAG3xuHY0v6BuJ6JoxNGS+oWaXVPj\nd7MAqxG+AHLWFo9r09adGllVL0lKSNq0dackEcDAEAhfFC3TI7YwjhB3NL/QF7y9RlbV69ndewP/\n3oBCInxRlEyP2MI6Qux2IikfT/QYbggQMBRcoSgNNWILw+uZEo04KR8v5cgCDIk/ERQl0yO2sI4Q\nl9Qv1In25qTHTrTvUkPdAp9aBAQD084oStGIo0SKxws1YjP9eqb0Tpk/u3uvEj1n30/TyoZAT6UD\nJhC+KEpL6hcmrcFKZ0dsTSsbQvF6Js2uqSFsgSwRvihKpkdsjBAB9BdxHCd1xYRHzpxJqLPzVCFf\nIvBisQpJop9coK/c8bOfgnZJFd8pd+gnd2KxCpWXZx7XMvIF4JmwXlIFeC3g5R4AbBLWS6oArxG+\nADwT1kuqAK8x7QzAM2G9pKpQgrY+Du8QvkBItb68T7/+bbPRA3uYL6nyGuvjxY3wBQwxOcppfXmf\nHty8TRdW1kkyd2APwiVVAz+HFY31qp0zy3g7uClFcSN8AQNMj3KeemZPX/D2MnVgt3nTjVSfw4Ob\nt0mSpk2dbrQtrI8XN1ZiAANMVwF3pzmAF/uBPdXncGFlnX6zY4/xtnBTiuLGxwwYYHqUE03zl13s\nB/a0n0O34YaIm1IUuyL/UwTMMD3KWbb0Bp18d3fSYxzYh/gcooYborPT800rb1FZR4siH7eorKPF\nuvVxFA5rvoABqaqA3/vvLbpl4bUFeb3eAqLHt+22tvDJD6k+h5PvNGv12kZf2mPz+jgKi72dLcCe\nqe4Fua9+tXmznm5u0bALxsnp6daoS6sV6XxPTStv8fwAHOR+KrS2eDypGnv5D+pUO2cWfZUB3yl3\n3O7tTPhagC+1e0Hqq4GXtBw9+v900ZXLBj2vrKNFWzY+5OlrB6mf/EZfpdf/OzysrETLlt5gvCo8\naNyGL2u+QAH0XtKSGDdfziXzlBg3X//839M6/o/2Qc8t9gpk2Gngd/j06Fo9uHmb2uJxv5sWCoQv\nUACpLmmZOHO1jv3trUHPLfYKZNgp3WVZ3CTDG/zZAwWQ9pKWUyeS/psKZNiKTUAKi2pnoADS3WBg\n/Fe/orKOFiqQYT1uklFYhC9QAOlvMNBI2CIQ0l2WtWbFLT62KjwIX6AAgnCDAWAoA7/D55WXaPXa\nRqqdPZLXpUbLly/X3r17NXr0aB06dCjlc7jUKDMudXCPvnKHfnKPvnKHfnLHyKVGt99+u1paWvL5\nFQAAFJ28wvfaa6/VRRdd5FVbAAAoCgVf8y0tjfZNVyC10i93daefMqOv3KGfUmt9eZ+eemaPunvO\n3vlp2dIbdP28OZLoq0z4TrlT6vIuHRRcASgKrS/v04Obt+nCyjpJUpekBzdvUzQa0XVzv+dv41B0\nCh6+iUQ3C/QZDCxkGLgn8JL6hVTJfomiD3fop8F+/dvmvuDtdWFlnZ54+nl9b9ZM+ioDvlPuuC24\nYuRrmd79VHuvrUtI2rR1pyQRwEAe0u7Y1G24IYDyLLhaunSprrnmGv31r3/V+PHj9cQTT3jVrqKV\naj/VkVX17KcK5CkaSX1VpcslOsBTeY18n3nmGa/agS+xnyq80FtY9EVXD0sXX0q369g9Tbf52CoU\nK6adLcN+qv4Kw3p7WzyuzY8/17e+ydLFWel2HaudM8vnlqEYEb6WSb8ncIOPrSoOYVlv39H8wqDC\not6liyC9j0KYXVNT9H0AOxC+lmFPYP8Mtd4epP4v1NJFGGYFAFsQvhbi7NwfYVlvL8TSRVhmBQBb\nsJIIfCltNWzA/kqW1C/UyXd3Jz12on2XGuoW5Pw7qcIHvBWwwwpQOEvqF+pEe3PSY/mGlh9m19To\n7jU/UMXRVkU+blFZR0veSxdhmRUAbMG0M/ClMK23186Zpdo5swbtRpTrui1V+IC3CF+gnzCvt+ez\nbksVPuCtiOM4qRe6PHLmTIK9QDNgz1T3bO8rWyqCU/XTqqb1SoybP+i5ZR0t2rLxoYy/sy0eT5oV\nqJzwNb3z17/5/l7zZft3yhb0kzvs7QxkKd/gtL0iON912/6zAra/V8B2rNgAOhcmiXHz5VwyT4lx\n87Vp6061xeOuf4ftFcFeVnPb/l4B2xG+gLwJE9srgr2s5nb7Xtvica1qWq8Va+/Sqqb1WZ3MAGHG\ntDMgb4KzkBXBXqwle1nN7ea9MjUNpEf4AvImOAtVEZxtiLXF49r1Xy+qu0dyerqTgtqram437zUs\n23UChUD4AvImOAt1nXA2ITYwqKXCjDbdvFfbp+EBPxG+gLwLzkJcJ5xNiJkcbWZ6r2zMAaRH+AJf\nsnWDjWxCzKbRJhtzAOkRvgg9Wza+yFU2IWbTaDNM23UCXiN8EWphqLjNJsRsG23aOpsA+I3tJS3A\ntm3uZdtX+W6pGERt8biaf/eiEt2SnB411C0gAIfA35879JM7bC8JyK41UFNm19Toxu/Pk8SBErAV\ndYcINS+3VAQAr3AIQqh5uaUiAHiFaWeEGhW3AGxE+CL0qLgFYBumnQEAMIzwBQDAMKadYZWg70YF\nAG4QvrBGGHajgh04iYPtCF9YIyz3f+XA7y9O4hAEhC+GZDJIwrAbFQd+/4XlJA7hRsEV0uoNksS4\n+XIumafEuPnatHWn2uLxgrxeGHajGurADzPCcBKH8GPki7RMjyBsuyNPLjjw+8/kbRVznRliaQKE\nL9IyHSRh2I3Kpvvp2sZU4Jg6ict1icG2pQlOBPxB+CItP4Ik6LtRhWH0XggmA8fUSVyuM0M2rUnb\ndiJQTAhfpEWQZC8Mo/dCMB04Jk7icp0ZsmlpwqYTgWJD+CItm4Kkd2osUhJVtES6afF8aw8OQR+9\nF4JNgeOVXGeGbFqaCOPnEhSEL4ZkQ5AMnBrrUnimxvqvt5049k9Fo6WKXfSV0K292RQ4Xsl1Zsim\nGaUwfi5BQRfDemG9fKf/pVzHnPH65ydRXTDpViOXdZkWxvsqz66pUdPKW1TW0aLIxy0q62hxNTOU\n688VQhg/l6Bg5AvrhXVqrP9JxdG/vaWJ3/1B0r+Hbe2t1Dmtw//9qLq7z2jsxRepafXywL+3XGeG\nbJhR6m2HZMfSUrEhfGG9sE6N9T+pKCmJpnxO0E8wpH7LBpNu1QVfPjZwtJXp57O9FIbLZ9yz5USg\n2BC+sJ7J6zZNHrD7n1T09HSnfI5XJxh+hlE+FbW5XArD5TMIgrz+tPfv36/KykpNmDBBjzzyiFdt\nApL0XyMrO9KqiqOtnk+Nmd5KU0peb7v40mq99+q2pH/3au3Nj/fWXz7LBrms94e1RgDhktfI98c/\n/rG2bNmiSy+9VNddd52WLl2qUaNGedU2WMTvabzeqbFYrEKS1Nl5Kuffleq9ZDM686ov+q+3jYpI\n0Qu69a/27YpdONLTtTe/r+XMZ9kgl+AOSo2A339T8FfO4dvZ2SlJmjFjhiSptrZWBw4c0IIFVMmF\nTZim8dK9l6i6FBs3+PkDD9he94XNm0F4JZ9lg1yCOwg1AmH6m0Jucg7f119/XRMnTuz77yuuuEKv\nvfbaoPAtLY32jVaQWmnp2WIbW/tp13+9mHLk1Py7F3Xj9+cZbUu+fZXuvRxu26hY1eDnn1dekvRa\nNvXFUPr307CyEp1O8ZyB761Qbvz+PJ0/fJh+s2OPEt1SaVS6p+k21c6ZlfFnVzTW68HN23RhZV3f\nYyffadbdaxvTtj3bn/Hj7y8o36P+bD9O2aK3nzI+r8DtQAh0pxkhJVLXCFkt3Xu5eNRInXx396AD\n9uq1ja5+3ua+WLb0hpRhNPC9FVLtnFmuwjbVz0lKCu7VaxuH/F25/IxpQfwewVs5h+/UqVO1fv36\nvv9ub2/XvHmDz9gSie681ueKgRfrmIXkpKnEldNjvM359lW693LxxReroW5B0vWOa1bcomlTpye9\nVjZ94eeaXv9+mjZ1utYs/yLje7PVtKnTNW3q9KTHMrU7m5/x4+/Ppr8pt2w/TtkiFqtQeXnmaM05\nfGOxmKSzFc9f+9rX9NJLL+knP/lJrr8OFrNpO7x8DfVe3Ky/uu0L29b0uJbTLmH6mzIhjMVpeU07\n/+pXv9Lq1avV1dWlpqYmKp1DKky74OT7Xtz+vN8VxrBbmP6mCs22E1mvRBzHcQr5AmfOJJimyIDp\nHPeC0lcr1t4l55LByzCRj1v060ceKvjrB6WfbEBfueNXP61qWq/EuPmDHi/raNGWjYX/W8qW22ln\ni4rvgfCIRlKf09p0uQsQBH5fKlcoVDvDSkFf42FND/BGEK7bzgXhC+sMtcbjxTWQJoKdNT24FfQT\nzUIL64ks4QvrDFWslG/4mizeoMIYmYS1mMhLYT2RJXxhnUKu8VCFDJvwfXQnjCeyhC+s49UaT6rp\nvLAWbyCY+D4WL8IX1sm0xtP68j79+rfNQ66RpZvOK3VO64JLBr/m+++9p7Z4PNBn1/1PNoaVlWjZ\n0hsG7fIEu4S1mAiZEb7wxVBFJkOt8bS+vC9pn+J0a2TppvNO/uUpnWhvTvq39179rS6+fG6g19oG\nnmyclvTg5m1as/wLq94PxUXJwlpMhMwIXxjnpsgk3RrPU8/sSbpBgJR6jSzddN5FXxmjW29aoP/z\nf38pDR8rp6dbF192lb7yb1XSv1X1/Z6ghUSqk40LK+usWjukuGiwsBYTITPCF8aDJp8ik7R3gxnw\n+FDTebNravT0rr0pd6BK9AQzJLxYO8zne+DmZykuSi2MxUTIjJWFItcbNIlx8+VcMk+JcfO1aetO\ntcXjBXvNfIIimuYbO3CNbEn9Qp1ob0567ET7LjXUnb3f9FA7UA0VErbKd0etfL4Hbn+W4iLgHMK3\nyPkRNPkExbKlN+jku7uTHusfqr1m19SoaeUtKutoUeTjFpV1tCRN5w0VzkEMiVTv5+Q7zYP6JZ18\nvgduf5YtN4FzmHYucn4ETT5FJr03RH982+6Ma2RDTecNtda2o/kF3ytQs50CHvh+zisv0eq1ja6r\nnfP5Hrj9WYqLgHMI3yLnx6UO+RaZ1M6Z5cklNOnC2e+QyHXNuf/7yfYONPl8D9z+LMVFwDncUtAC\nft7SbOCBXkq+ubxtTPVVWzyeFBINdQuM9YcXt1DLtp/y+R6Y/g55XSDILQXdoZ/ccXtLQUa+RS7X\n0UjQLsXJlp8VqH4sBeQzKjU5og1iJTqQCuFrGT9CLdug4QBYWPlMAfd+fyIlUUVLpJsWz3f9meRz\nwmHqZIXLlRAW1BlaxI/LfnIRxEtxgiTTZVLp9P/+dI2t1enRtVZ+f/IRxEp0IBVGvhYJyll9qgPg\n8X+06+j772nF2rtCOQ1tUq7TuEH5/uSDvZARFoSvRYJyVj/wAHj8H+06+uGbmvi9dXLENLQXcpnG\nDcr3Jx9+V6IDXiF8LZLNWb2fBU8DD4BH//aWJn73B0nPCduIKwiKYVTI5UoIC8LXIm7P6v0ueBp4\nAIx2daZ8XphGXLbqfxLW+b8ndOJ//lOX/ceqvn8v9KgwCAWCgI0IX4u4Pav3e21v4AF3zKhYyueF\nacRlo4EnYRdcIn1ycJtO/uUpXTz6qyqNFnZU6PdJoJ/CfqkdCo/wtYybs3o/1/ZSHXBP/M9/6pOD\n23TJv5+bemYdrvBSnYRd8u8/UFlHi3Y8sVFSYTdE8Psk0C/FfNIB7xC+AVSItT23Z/KpDriX/ccq\nnfzLUyrraGEdziC/C6z8fn2/FOtJB7xF+AaQ1xWf2ZzJD3WTerdbH8IbfhdY+f36finWkw54K+R/\nJuGU6XZ52cpm0wxuC2ePXDfjCMvr+4W/AXiBkW9AeVnxmc2ZPNdZ2sPvy278fn2/8DcALxC+yGr6\nsFgPuLZItTbv53R/MV72w98AvED4Iusz+WI84NqAKlt78DeAfBG+4Ew+IKiyhc249jk7hC8kcSYf\nBFTZwlbMymSP8AUColgv7emP0ZWdmJXJXhH92QLBVqyX9vQKyv2uixGzMtlj5AujGLnkrtjX5hld\n2YtZmewRvjDGhnWhoId/Ma/NM7qyF9c+Z4/whTF+j1xsCH/kjtGVvYp9ViYXhC+M8Xvk4nf4Iz+M\nruxWzLMyuSB8YYzfIxe/wx/5YXSFMCF8YYzfIxe/wx/5Y3SFsCB8YYzfIxe/w992QS9Gy1Xry/v0\n1DN79EVXT1G9b/iL8IVRfo5c/A5/m9lcjFbIk4K2eFybH39OF1bWSbLrfSPcIo7jpL45ZQbPPfec\nfvrTn+q9997T66+/rquuuirl886cSaiz81RejQy7WKxCkugnF+grd7Ltp1VN65UYN3/Q42UdLb7e\nNWngSYEknWhvVtPKWzwJR1vft43423MnFqtQeXnmcW3Oq12TJ0/W888/rxkzZuT6KwBYwtZitKEq\n1L1g6/tG+OU87Txx4kQv2wHAR7YWoxU6HG193wi/gq/5lpZG+6YrkFppaVSS6CcX6Ct3su2nFY31\nenDztr61T0k6+U6z7l7b6GtfDysr0ekUj59XXuJJu1Y01uuhzU8rVnlj32M2vG8b8bfnTm8/ZXze\nUP84d+5cHTlyZNDjDzzwgBYtWpRby+CL3orO7h4pWiItW3qDaufM8rtZsETvd+E3O/Yo0S2VRqXV\naxt9/44sW3pDypOC1WsbPfn9tXNmKRqN6Mnte9SVcKx53wi/nAuues2aNUu//OUvKbjKQ6ELGQpd\ntGISRR/umOynQl+i1BaPJ1WoN9Qt8PT3851yh35yx23BlSfTznnmNwosSNsqFuu1pkFl4hIlNtZA\nGOVcVvD8889r/Pjxeu2117RgwQLNnz+4XB92CEpFJ/drDZ5CVyMDYZVz+N5444366KOPdOrUKR05\nckQvvviil+2Ch6KR1DMTtlV0ciAPnqCc2AG2YYerIhCUbRW9PJAzfW1GLpfq8NkAhG9RCMq2il5d\nc2nzVolhk+2JHZ8NcFbe1c6ZUO2cGVWEZ6Wuyt6VdKLgpq/YMtB8tbPbamQbPxv+/tyhn9wxWu0M\neMGrETrrkGZlU43MZwOcRfjCKl5cVsKWgfbiswHO4iuP0FlSv1An2puTHjvRvksNdQt8ahF68dkA\nZzHyRegEpcCsGPHZAGdRcGUBChnco6/coZ/co6/coZ/cKfj9fAEAQG4IXwAADCN8AQAwjIIrH/Vu\nsxcpiSpaIt20eD6FJxZiO0QAXiN8fTJwN6cusc2ejdgOEUAhMO3sE+7gEwx8TgAKgZGvT9hmLxjc\nfk5MTQPIBuHrE7bZCwY3nxNT0wCyxaHeJ2yzFwxuPiempt1ri8e1qmm9Vqy9S6ua1qstHve7SYAv\nGPn6pP82e4qUqDTKNns2crMdYhCWEGyYFmeGADiH8PVR7x182LbNbpnutGT7EoItoTfUDAHhi2Jj\nyeEBCC7blxBsmRYPwgyBbZimDy9GvoDym5a1/U49toRephkCG6bGbWLLjAUKg/BF0fPiIJdpatpP\nbiu2Cx18S+oXJvWzdHaGoGllA0GTAtP04Ub4ouhHHGE/yA0VepK5EdZQMwSrmtaH+jPIhS0zFigM\nwrfIMeII/0Eu07S4yZOPdDMEYf8McmF7IR/yQ/gWubCP+twohoPcUNPiNgRfMXwG2co0Y4FgK+Kv\nNiQ7Drx+s71audCiESfl4yaDr9g/g1Rm19SoaeUtKutoUeTjFpV1tFhVyIf8MPItcrmMOMK2Rmx7\ntXKh2TDCKvbPIB2bC/mQn4jjOKlPez1y5kyCzSMy8HOTjYFrvtK5A2+qP/rUz29W08pbjBwk2JDE\nnWz7qS0eTwq+hroFRXPQ5zvlDv3kTixWofLyzONawtcCfn+psznwrmpar8S4+YMeL+to0ZaNDxW6\nqb73VVDQT+7RV+7QT+64DV+mnZHV1BZrxACQPwqukBUbinMAIOg4ZCIrVKUCQP6YdkZWqEoFgPwR\nvsgalz8AQH6YdgYAwDDCFwAAwwhfAAAMI3wBADCM8AUAwDDCFwAAwwhfAAAM4zpfGBG22xACQD5y\nHvmuX79elZWVuuqqq7Ru3TqdOsWdLpBa720IE+Pmy7lknhLj5mvT1p1qi8f9bhoA+CLn8K2trVV7\ne7sOHjyozz77TNu3b/eyXQiRHc0vJN3/V5JGVtXr2d17fWoRAPgr5/CdO3euSkpKVFJSouuuu05x\nRjFIg9sQAkAyT9Z8H3vsMa1cudKLX4UQikYcJVI8zm0IU2N9HAi/IcN37ty5OnLkyKDHH3jgAS1a\ntEiSdP/992vEiBG6+eabU79AaVSxWIUHTQ2v0tKoJIW2n1Y01uvBzdt0YWVd32Mn32nW3Wsbs37P\nYe+r1pf3afPjz/X1VULS5sef0/nDh6l2zizXvyfs/eQl+sod+smd3n7KJOI4Tuq7o7vw5JNP6rHH\nHtMrr7yi8847L+VzenocJRLdub5EUej9sMLcT60v79NvduxRolsqjUq3LbkhqzDpFfa+unXFj3V6\ndO2gxyuOtmrb1o2uf0/Y+8lL9JU79JM7paVRlZSkXmpLel6uL9DS0qKHH35Y+/fvTxu80tkPqrOT\nSuih9J5Jhrmfpk2drmlTpyc9lsv7DXtffdGVeiH89JmerN5z2PvJS/SVO/STO7FYhcrLM0drzqtu\na9eu1aeffqo5c+ZoypQp+tGPfpTrrwLwpWgk9UQU6+NAuOQ88j18+LCX7QAgaUn9Qm3aujPp0qwT\n7bvUtLLBx1YB8Bo7XAEW6a1qfnb3XiV6zo54m1Y2UO0MhAzhC1hmdk0NYQuEHCtJAAAYRvgCAGAY\n4QsAgGGELwAAhhG+AAAYRvgCAGAYlxoBCATu9oQwIXwBWK8tHk/a+SshadPWnZJEACOQmHYGYL0d\nzS8kbbkpSSOr6vXs7r0+tQjID+ELwHrdTupbtCVS3wQKsB7hC8B63O0JYcNXF4D1ltQv1In25qTH\nTrTvUkPdAp9aBOSHgisA1uNuTwgbwhdAIHC3J4QJ084AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG\n+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACG\nEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCAYYQvAACGEb4AABhG+AIAYBjhCwCA\nYTmH73333acrr7xS1dXVamxs1PHjx71sFwAAoZVz+N511116++239dZbb2nChAnauHGjl+0CACC0\ncg7fESPpBX3rAAAEu0lEQVRGSJISiYQ+++wznXfeeZ41CgCAMIs4juPk+sP33nuvtmzZossvv1z7\n9u1TeXn5oOf09DhKJLrzamTYlZZGJYl+coG+cod+co++cod+cqe0NKqSkkjG5w0ZvnPnztWRI0cG\nPf7AAw9o0aJFkqTPP/9c9957ryRpw4YNubYXAICikdfIt9ehQ4d0xx136LXXXvOiTQAAhFrOa76H\nDx+WdHbN95lnnlFdXZ1njQIAIMxyDt977rlHkydP1jXXXKNEIqE77rjDy3YBABBaOYfvrl27dOjQ\nIf3pT3/SQw89pIsuuijl87ge2L3169ersrJSV111ldatW6dTp0753SQrPffcc6qqqlI0GtUbb7zh\nd3Oss3//flVWVmrChAl65JFH/G6OtZYvX64xY8Zo8uTJfjfFah999JFmzZqlqqoqzZw5U9u3b/e7\nSVY6ffq0pk2bpurqak2fPj1zDZRTYP/617/6/v/PfvYz57777iv0SwZWa2ur093d7XR3dzsrV650\ntm7d6neTrPTuu+8677//vjNz5kznz3/+s9/NsU51dbUTj8edDz/80Ln88sudo0eP+t0kK+3fv995\n4403nEmTJvndFKt1dHQ4b775puM4jnP06FHn61//etJxHed89tlnjuM4zunTp52qqirn8OHDaZ9b\n8O0luR7Yvblz56qkpEQlJSW67rrrFI/H/W6SlSZOnKhvfvObfjfDSp2dnZKkGTNm6NJLL1Vtba0O\nHDjgc6vsdO2116adscM5Y8eOVXV1tSRp1KhRqqqq0sGDB31ulZ2GDx8uSfr000+VSCQ0bNiwtM81\nsrfzvffeq7Fjx+rVV1/VnXfeaeIlA++xxx7ru5wLcOv111/XxIkT+/77iiuu4CoEeOaDDz5Qe3u7\nrr76ar+bYqWenh5deeWVGjNmjNasWaPx48enfa4n4Tt37lxNnjx50P9+//vfS5J+/vOf6+9//7uu\nvvpq3X333V68ZGBl6itJuv/++zVixAjdfPPNPrbUX276CYA5n3zyiRoaGrRhwwadf/75fjfHSiUl\nJXr77bf1wQcf6NFHH9Wbb76Z9rmlXrzgSy+9lPE5w4cP1/Lly4u+KjpTXz355JP6wx/+oFdeecVQ\ni+zk5juFwaZOnar169f3/Xd7e7vmzZvnY4sQBl1dXaqvr1djY6MWL17sd3Osd9lll+n666/XgQMH\nNGXKlJTPKfi0M9cDu9fS0qKHH35Yv/vd71gbd8nJf4+YUInFYpLOVjx/+OGHeumllzRt2jSfW4Ug\ncxxHK1as0KRJk7Ru3Tq/m2OtY8eO6eTJk5Kk48ePq7W1dcgTFU92uBrKTTfdpPfff18VFRWaOXOm\n7rnnHooc0pgwYYLOnDmjkSNHSpK+853v6NFHH/W5VfZ5/vnn1dTUpGPHjikWi2nKlCl68cUX/W6W\nNeLxuH74wx+qq6tLTU1Nampq8rtJVlq6dKni8biOHz+u0aNH6/7779ftt9/ud7Os8+qrr2rGjBn6\n1re+pUjk7J7Fv/jFL5hRGeDQoUNatmyZuru7NXbsWN1666267bbb0j6/4OELAACSGal2BgAA5xC+\nAAAYRvgCAGAY4QsAgGGELwAAhhG+AAAY9v8B4oY6ZvGrkZUAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# We can check this assertion formally in Python with the caveat that there\n",
      "# is some very small numerical rounding errors.\n",
      "assert all(X1 - X1b < 1e-5), \"Matrix 1 does not equal matrix 2\"\n",
      "\n",
      "# We can also get built-in assertion with consideration of floating point precision\n",
      "from numpy.testing import assert_array_almost_equal\n",
      "assert_array_almost_equal(X1, X1b)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Another way to say this is that the matrix obtained by multiplying A with\n",
      "# its inverse is a matrix that will not change the values of X. This matrix\n",
      "# is known as the identity matrix and contains 1s on the diagonal and 0s in\n",
      "# every other entry.\n",
      "I = dot(A, inv(A))\n",
      "assert all(I == eye(I.shape[0]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The formula $w = (X^TX)^{-1}X^Ty$ is used for fitting an ordinary least\n",
      "fsquares regression in any statistics package. We can show that this\n",
      "formula obtains the same results as the function ``lstsq`` (which\n",
      "lives in the ``numpy.linalg`` namespace normally."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Fabricate a dataset with 4 variables\n",
      "X = randn(100, 4)\n",
      "# Fabricate a variable y that depends on these four variables and some noise\n",
      "noise = randn(100) * 10\n",
      "y = 1 * X[:, 0] + 7 * X[:, 1] + 3 * X[:, 2] + 5 * X[:, 3] + 7 + noise\n",
      "\n",
      "# To capture the constant term, we will add a column of 1s to our regressor matrix\n",
      "X = column_stack((X, ones(100)))\n",
      "\n",
      "# Fit the model based on the ordinary last squares (OLS) formula\n",
      "w_ols = dot(dot(inv(dot(X.T, X)), X.T), y) \n",
      "\n",
      "# Note how the lack of builtin matrix multiplication for arrays is a bit annoying\n",
      "# We can also use matrix objects\n",
      "X_mtx = asmatrix(X)\n",
      "y_mtx = asmatrix(y).T\n",
      "w_ols2 = inv(X_mtx.T * X_mtx) * X_mtx.T * y_mtx\n",
      "\n",
      "# We see these are the same (once we account for the 2D nature of the matrix\n",
      "assert_array_almost_equal(w_ols.reshape(w_ols2.shape), w_ols2)\n",
      "\n",
      "# Use numpy's lstsq function to fit the model\n",
      "# lstsq returns a 4-tuple\n",
      "w_lstsq, residue, rnk, s = lstsq(X, y)\n",
      "\n",
      "# You can also do\n",
      "w_lstsq = lstsq(X, y)[0]\n",
      "\n",
      "# If you are following the MATLAB tutorial, note that the argument order for lstsq is backward relative to regress\n",
      "\n",
      "# These are also identical\n",
      "assert_array_almost_equal(w_ols, w_lstsq)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can see how well the model describes the data be computing the mean\n",
      "squared difference between our model's prediction and our data."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# First we compute our models prediction\n",
      "modelprediction = dot(X, w_ols)\n",
      "\n",
      "# Visualize this prediction by plotting it against our data\n",
      "# The imperfections arise from the noise we added to the data\n",
      "plot(modelprediction, y, \"o\", color=colors[1])\n",
      "xlabel(\"model prediction\")\n",
      "ylabel(\"original data\");"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAFgCAYAAABXHWtRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtclHXeP/7XzCB4XNRSsQhL1wAFdVgB28REw5Q0PEIn\nzMCzpmy/6LFuh03vR3aX9VXLILZiC+2guRl531oC4SmLYYssucGV1MUTaZgIpiLD9fuDZWKcGeaa\nYa65DvN6Ph49kmsO13veHN7X53jpBEEQQERERJqilzsAIiIi8jwWeCIiIg1igSciItIgFngiIiIN\nYoEnIiLSIBZ4IiIiDZK1wJvNZhiNRkyZMgUAUF9fj6SkJISEhGDq1KloaGiQMzwiIiLVkrXAr1+/\nHkOGDIFOpwMAZGdnIyQkBEeOHEFwcDDeeOMNOcMjIiJSLdkK/MmTJ7Fjxw7MnTsXrXvtmEwmpKen\nIyAgAGlpaSgpKZErPCIiIlWTrcD/6U9/wpo1a6DX/xZCaWkpwsLCAABhYWEwmUxyhUdERKRqfnKc\n9H/+53/Qt29fGI1G7N6923Jc7K65zc0CmprMEkWnHX5+BgBgrpxgnsRjrsRhnsRhnsTz8zNAr9e5\n9hqJYmnXgQMH8Omnn2LHjh24cuUKLl68iNTUVERHR6OiogJGoxEVFRWIjo62+/qmJjPq6i57OWr1\nCQzsAgDMlRPMk3jMlTjMkzjMk3iBgV3g7+9ayZali3716tU4ceIEjh07hg8//BDjxo3Dxo0bERsb\ni9zcXFy+fBm5ubkYNWqUHOERERGpniLWwbfOol+0aBGqq6sRGhqKU6dOYeHChTJHRkREpE46Nd4u\ntrGxiV06IrD7SxzmSTzmShzmSRzmSTzVdNETERGRtFjgiYiINIgFnoiISINY4ImIiDSIBZ6IiEiD\nWOCJiIg0iAWeiIhIg1jgiYiINIgFnoiISINY4ImIiDSIBZ6IiEiDWOCJiIg0iAWeiIhIg1jgiYiI\nNIgFnoiISINY4ImIiDSIBZ6IiEiDWOCJiIg0iAWeiIhIg1jgiYiINIgFnoiISINY4ImIiDTIT+4A\niIiIfEHxviJs3rEZZl0zDIIeKYkpiI8bL9n5WOCJiIgkVryvCBu2ZiEoaRAM/zm2YWsWAEhW5NlF\nT0REJLHNOzYjKGmQ1bGgpEHYsnOLZOdkgSciIpKYWdds93gTzJKdkwWeiIhIYgbBfrn1s3TYex4L\nPBERkcRSElNQk/+j1bGa/CokT0qW7JycZEdERCSx1ol0W3ZuQRPM8IMBS2cu0eYs+itXruCuu+7C\n1atX0blzZ6SkpOBPf/oT6uvr8fDDD6OsrAxRUVHYtGkTunfvLleYREREHhEfN17Sgn492broO3fu\njOLiYnz33XfYs2cP3n77bRw5cgTZ2dkICQnBkSNHEBwcjDfeeEOuEImIiFRL1jH4rl27AgAaGhrQ\n1NSEgIAAmEwmpKenIyAgAGlpaSgpKZEzRCIiIlXSCYIgyHXy5uZmGI1GlJeXY926dVi6dCkGDBiA\nw4cPo3Pnzvj1118RHh6Of//739e9TkBTk3RLC7TCz69ldiZz1T7mSTzmShzmSRzmSTw/PwP0ep1r\nr5EoFlH0ej0OHjyI48ePIzExEXfeeSdkvN4gIiLSDEXMor/11luRmJiIkpISREdHo6KiAkajERUV\nFYiOjrZ5flOTGXV1l2WIVF0CA7sAAHPlBPMkHnMlDvMkDvMkXmBgF/j7u1ayZRuD//nnn3HhwgUA\nQG1tLXbt2oWkpCTExsYiNzcXly9fRm5uLkaNGiVXiERERKolW4E/c+YMxo0bh+HDh+PBBx/EE088\ngf79+2PRokWorq5GaGgoTp06hYULF8oVIhERkWrJOsnOXY2NTezSEYHdX+IwT+IxV+IwT+IwT+Kp\nqoueiIiIpKOISXZEROS64n1F2LxjM8y6ZhgEPVISU7y6U5qUtPzZvIUFnohIhYr3FWHD1iwEJQ2y\n3I9sw9YsAFB9IdTyZ/MmdtETEanQ5h2bEZQ0yOpYUNIgbNm5RaaIPEfLn82bWOCJiFTIrGu2e7wJ\n6t8VTsufzZtY4ImIVMgg2P/z7Wfp1FYvLX82b2KBJyJSoZTEFNTk/2h1rCa/CsmTkmWKyHOU/NmK\n9xVh4Yr5mPeXuVi4Yj6K9xXJHZJDnGRHRKRCrZPNtuzcgiaY4QcDls5coolJaEr9bGqb/MeNbjSM\nm0iIwzyJx1yJwzyJo7Y8LVwxH4aEnjbHmwvrkL06R9Jzu7PRDVvwREQqxbXi3mXWNdudBaDUyX8s\n8EREKqS27mItUNvkP06yIyJSIa4V9z4lT/6zhy14IiIV8kZ3MYcArCl18p8jLPBERCokdXcxhwDs\ni48br5rPzy56IiIVkrq7mEMA6scWPBGRCkndXay2GeNkiwWeiEilpOwuVtuMcbLFLnoiIrKhthnj\nZIsteCIisqG2GeNkiwWeiIjsUtOMcbLFAk9Eqsf12kS2WOCJSNW4XpvIPk6yIyJV43ptIvtY4IlI\n1cy6ZrvHuV6bfB0LPBGpGtdrE9nHMXgiUrWUxBTLGHyrmvwqLJ25RMaofAcnOCoXCzwRqRrXa8uH\nExyVTScIgiB3EK5qbGxCXd1lucNQvMDALgDAXDnBPInHXInjK3lauGI+DAk9bY43F9Yhe3WO09f7\nSp48ITCwC/z9XWuTswVPROQidku34A1plE22SXYnTpxAfHw8hg4dirFjx+L9998HANTX1yMpKQkh\nISGYOnUqGhoa5AqRiMhGa7e0IaEn/O/uDUNCT2zYmoXifUVyh+Z1nOCobLIV+E6dOmHt2rUoLy/H\n1q1b8fTTT6O+vh7Z2dkICQnBkSNHEBwcjDfeeEOuEImIbHDd/W94Qxplk62LPigoCEFBQQCAG2+8\nEUOHDkVpaSlMJhOefvppBAQEIC0tDS+88IJcIRIR2dBit7S7Qw6enuDIoQ/PUsQYfFVVFcrLyxET\nE4NHH30UYWFhAICwsDCYTCaZoyMi+o3WuqU7OhPeUzek4Yx8z5O9wNfX1yMlJQVr165F9+7dIWZS\nv5+fwTL7khzz82v5NWGu2sc8icdcAWmzZmPNxnXoO2Wg5djZ7UeRmZphyYua8vSPXVvtDjl8XLAV\nUydPlvTcbfMkZxxq0Jorl14jQRyiXbt2DTNmzEBqaiqSkpIAANHR0aioqIDRaERFRQWio6PlDJGI\nyEpCfAIAYOMn71m6pTNTMyzH1casM9udjOXtIQelxKElshV4QRCQnp6OiIgIZGRkWI7HxsYiNzcX\nL730EnJzczFq1Cib1zY1mbluUgSuMRWHeRKPuWoREzUaMVGjrY61zYma8iRcc/BAk07y+NvmSc44\n1MCddfCyzaL/8ssvsWnTJnzxxRcwGo0wGo347LPPsGjRIlRXVyM0NBSnTp3CwoUL5QqRiMgrivcV\nYeGK+Zj3l7lYuGK+V5fcKWUmvFLi0BLZWvCjR49Gc7P9u0Dl5+d7ORoiInnIPbms9Ryvv7sBP104\nC0MnA/oG9pH8vI7i4JbDniP7JDsi8h1cBmWrvXX1ns5Ne/kXuukx7ME7Lc+VYwa7p2bkUwsWeCLy\nCrlbqkrlrXX17eXfmxcZ5D0s8ETkFSwi9rWuqz9XcRo/fX8SOoMegrkZfcy9PHqe9vKvxc17SMZJ\ndkTkW8w6+3NufL2IpCSm4PDfv0HNwROISInB0JkjEZESg18NVzw62a69/Gtt8x5qwRY8kcJpZdxa\n7UVEqu9DfNx4ZL+XhT7J1q3r21IiPdq70V7+kxOTLd33rWryq7B05hKPnJvkwRY8kYJp6c5lal4G\nJfX3IbC3/e54T/ZutJf/+LjxWDpzMZoL69BYeB7NhXWcwa4BbMETKZhU49Zy9AqoeRmU1PMHvNG7\n4Sz/nMGuPSzwRAomxeQnOWezq7WISD0JLSUxxStd5GrNP7mHBZ5IwaRo2XE2u+ukbmGruXeDlIsF\nnkjBpGjZcUmU67zRwmbrmjyNBZ5IwaRo2al9Nrsc2MImNdIJYm7ArjCNjU28u5AIarqjlZx8LU9t\nx+BbtbZGnRUsX8uVu7yVJ7UvoeTPk3ju3E2OLXgijXL0x5+tUW3g1r/kDAs8kQY5++PP8V7142RJ\ncoYb3RBpUHt//EkbuPUvOcMWPJEGcaa89nljsqTax/h9HVvwRBrEmfLaJ/XWv1raJtlXscATaZCa\n930ncaTeP57DPOrHLnoiDeJMee9r250dYOiE2dMeQkzUaI+/9/UrIqT6nnKYR/1Y4Ik0Sokz5aUe\n05VrzPj6VQsCgDUb12HxpasdPr9cy+E4zKN+LPBE5BVSFyo514Xb687uO2WgR5asubscrqMXO+5u\nz8uJecrBAk9EXtGRddtiioac68Kl7M525709cbHjzjAPN99RFhZ4IvIKd4ug2KIh55ixlN3Z7ry3\n2IsdZxdOrg7zcPMdZeEseiLyCneLoNjZ3HKOGdtbtXB2+1GPrFpwZ0WEmE1wpFgGx813lIUteCIF\n09J4prtjumJb5t64pasj13dnd/brhMzUDI/Monenq1zMxY4UrW1OzFMWFngihdLaeKa7S/fEFg25\nlwa27c729F3SXO0qF3OxI8WQhpwXWWSLBZ5IobQ4nunO0j1XioYSlwbKQczFjhStbbkvssgaCzyR\nQnGjkRYsGu5xdrEjVWubF1nKwQJPpFAcz/yNXEXD03MglDSnghdO2scCT6RQHM90nycKqafnQChx\nTgVb29om2zK5tLQ09OvXD5GRkZZj9fX1SEpKQkhICKZOnYqGhga5wiOSndQ3E3FV8b4ipC6bgweX\np2LhivmKvauYp5Z/efpmK7x5C3mbbAX+0UcfxWeffWZ1LDs7GyEhIThy5AiCg4PxxhtvyBQdkTLE\nx41H9uocvLn6LWSvzpG1uG/YmgUhvgf0Y3sq+tahniqknl7TzTXi5G2yFfi4uDj06tXL6pjJZEJ6\nejoCAgKQlpaGkpISmaIjorbU1Pr0VCH19BwIzqkgb1PUGHxpaSnCwsIAAGFhYTCZTHaf5+dnsKwz\nJcf8/Fr+cDBX7WOenNN1cvCAn6C4vAUYOkGwc7yzXyeXYk2bNRtrNq5D3ykDLcfObj+KzNQMp+9j\n72eqI++nVfzdE681Vy69RoI43CYI9n4tiUhuBsFgt2i2bX0WFBcgb9t7MOvMMAgGzJ72EBLiEzwe\ni7PzzJ72kMNC6orW99z4yXuWWeaZqRlufyZPvx+RM4oq8NHR0aioqIDRaERFRQWio6PtPq+pyeyx\nHaK0zNO7aWkV8+TcjAkzHc7or6u7bDVDXI+W+6G/+O5aXPLA/dDbEnOemKjRWHzpqtXyr8XTFyEm\narTL3+OYqNE2282KeQ9HP1Puvp9W8XdPvMDALvD3d61kK6rAx8bGIjc3Fy+99BJyc3MxatQouUMi\nIvy2jOvjgq0tY9lNOqsZ/d7adU/sebj8i0jGAv/AAw9gz549qK2txS233IJVq1Zh0aJFePjhhxEa\nGoqoqCi8+OKLcoVHRNeJjxuPqZMnA7BtcXlr1z3u7kcknmwF/oMPPrB7PD8/38uRkC8r3leEf+za\nCrPODOEaVH23Njn9crYW5zZXQWfQQzA3o9+wYPQJv8njM8Q5E51IPKcF/syZM3jrrbeQn5+P2tpa\nAIBOp8PRo0clD45ISteP5wLy7yymRsX7ivCr4QoiUmIsx374sATnvz6Dp+b/xaPn4u5+ROI5LfDP\nP/88goODcenSJWzYsAE5OTkYM2aMN2IjkpQW79Ymh807NuO2lEirY5H3x+LcR0c9nkfun64OStpz\n35c5LfBffvklysrKkJeXh3vvvRdxcXEYP348nnjiCW/ERyQZjud6hqM8BvbqKcn52ptAx8IiPyXu\nue+rnBb4Tp1adriIjIxEWVkZBg8ejPr6eskDI5Iax3M9Qyl5FFNYeAEgPfaMKYfTrWrHjBmD2tpa\nPPzww0hMTMTAgQMxZ84cL4RGJK2UxBTU5P9odawmvwrJk5JlikidlJJHZ9vpeuomNFIq3leEhSvm\nY95f5ir6hj7t4Z77yuG0Bf/yyy8DAKZMmYKqqipcuHABN998s+SBEUnN2dpuEkcp4+LOhlyU3rLU\nSte2Unp0SESBT05OxpYtLVfA3bp1Q7du3ayOEalZe2u7STwlbCzjrLAocc5F2yGDw/+qQFD8bVaP\nK+kCRCyudFAOpwX+yJEjNseqqqokCYaIyF3hA8Lwj9xtiEy7w3KsbWHxVstS7N4K17fYI+++Ez98\n2HIHzT7hN1mep7aubaX06FA7Bf7NN9/E3/72N/zrX/+y2hO+pqYGiYmJXgmOiEiM4n1F2FvxJYLu\nGIBDW0zQ6fW4cqYBD9xzv6WweKNl6creCvaGDCLvj8WhLSarAq/Grm0l9OhQOwV+woQJGDx4MJYs\nWYKXX37Zcqe3gQMHIiQkxGsBEqkRZ2t7V9ti2bY4VhYetvzbWcvSE98zV8b5HQ0Z6PS/9TSwa5s6\nwmGBHzBgAAYMGIDy8nJvxkOkelqZLKUmYsfXHbUsPfU9c2Wc39GQQfPZRjQWnmfXNnWY0zH4hoYG\nbNu2zWar2i+++ELy4IjUSOmztbWoo+PrnvqeuRKHoyGDlctX8ueEPMJpgV+1ahUaGhpQVlaGxx9/\nHO+++y4m/2fWMRHZUuJsba1LSUzB839bjWtdzZYb3nT61SB6L/zW79m5itP46fuTlvfobQ50OQ6x\n4/ycjEZSc1rgi4qKUFpaioiICCxZsgSzZs3Cfffdh2effdYb8RGpjq+uA/b0vANX38+/qz9C2+yJ\nf2zzD6LPZRD0OFdxGjUHT7Tso/+fQn+y9iSSF87AoocWi/osru6twMloJCWnBV6v10Ov1+P222/H\nsWPHEBISggsXLngjNiJV8sV1wJ6ed+Dq+9m74c1tKZGiu9hTElOw4v+twB8eG2tV6Fu58lm4twIp\nhdMCHxoaitraWkyfPh2jR49GYGAgJk6c6I3YiFRJzq5XuWbve3regavv19Fhkfi48bjlw1sAAD99\nf9KquDs7N5FSOS3wmzZtAgDMnj0b8fHxOH36NGJjY528isi3ydH12tFWdEcuDjw978DV9/PEsMgN\nv7sBAKAz2H8vzqEgtXFY4H/99VebYzfccANuuOEG/Prrr+jataukgRGRY/aKcUda0R29OPD0vIP2\n3s/eZ/fEsEjrewhm+zdL0focCtIeh3eT6969u8P/evTo4c0YiagNR3dFq71Ya/f5Ylqezu7E5oy9\nO8r9kPsVwkJCRb1ezPvV5FchLCTU7mcHgKUzF6O5sA6NhefRXFjn8rBIfNx4LJ25GH3MvfD937+y\n+Sw1Z2tUe4c38k0OW/DNzS1XsWvXrsXJkyexePFiCIKAnJwc3k2OSEaOivEPWQfQF4Nsni+m5emJ\nMezvyw/ig9c+ROf+3SE0NyPojgHYW/Elhu0b7vJwhaN5DO1diGSvzunwsEjr0ErxviJs2bkFP9f9\njBOnTyDk7tstO+Rx0yJSC6dj8Lm5uSgrK4OfX8tTX3jhBRiNRmRkZEgeHBHZclSMb+h5A2ryf3Sr\nm9oTXewV/67EHx4ba30wHG5PTrM3j+H9nR94ZY+B1nMvXDEffWdw0yJSJ4dd9K369+8Pk8lk+fqf\n//wn+vfvL2lQROSYo2Ic1DfI7W5qR13iyZOSRcdl1tkfu/Zk8fX2HgPe+ExEUnHagn/mmWcwa9Ys\n9OvXDzqdDjU1Ndi8ebM3YiMiO9qbUObu7H1PLO3zRvH19h4DvrppEWmD0wIfFxeHkydPoqSkBDqd\nDjExMdDpdN6IjYjskGqdfUeX9nmj+Hp7jwFf3LSItEMntN4HVkUaG5u4Q5QIgYFdAHA3LWeYJ/Gc\n5ap1clpr8U2elKz6sWp3PhN/psRhnsQLDOwCf3+nbXIrLPAaxl8ecZgn8Zgrx9quzw8wdMLsaQ8h\nJmq03GEpGn+exHOnwDudZEdERO27fm8CIb4H1mxcxzXzJCsWeCKiDrK3Pr/vlIGiNwoikoLD9v7r\nr78OnU4Hez34Op0OixcvljQwIpKeXDen0RpP78VP5AkOC3xpaalss+X37t2LBQsWoKmpCcuWLcNj\njz0mSxxEWubpW7z6Mi6nIyVS5CQ7o9GI9evXY8CAAbjnnnuwf/9+3HjjjZbHOclOHE5gEcdX87Rw\nxXwYEnraHG8urEP26hy7r/HVXDnT9mKp1dntR7F4+iJeLLWDP0/iuTPJTtSzr169isLCQtTW/nYz\ni9mzZ7sWnUh1dXUAgDFjxgAAJkyYgJKSEtx7772SnI9ISkruAnenW7mguAB5297DVfM1xX0eOV2/\nPr+zXydkpmZwFj3JymmB//DDD/Hss8/i3LlzGDx4ML755hskJiZKVuBLS0sRFhZm+XrIkCH4+uuv\nrQq8n5/BcuVHjvn5tfz5Zq7aJ1WeCooLkPVxNvpOGWgppFkfZ6NbtwAkxCd49FzuCDB0gr3uu85+\nnezmoqC4AC9vWoc+kwfC/z/HlPR55DZ18mRMnTwZwG8/U01NHINvD/9GideaK1c4nUX/9ttv4+uv\nv8Ytt9wCk8mEwsJC/O53v3MrQCJfkrftPfSdMtDqWN8pA7Hxk/dkisja7GkP4ez2o1bHzm4/itSp\nD9l9ft6299BncsvnOVdxGoc2m3Du6nlkvvBnFBQXtHuuguICpC6bgweXpyJ12RynzyeijnPagq+t\nrUXv3r3RtWtX1NfXIz4+HkuXLpUsoOjoaGRmZlq+Li8vx8SJE62e09Rk5piNCBzfEkeqPF01X7O0\ndNu60nRNEd+TmKjRWHzpqtUubYunL0JM1Gir+FqHGSqOVwBnDejcswsu//IrIu+PtTznxXfX4tKl\nq3a769uOT+sBCE6erwX83ROHeRJPkjH47t27o7GxEXfddReWL1+OQYMGSXo/+MDAQAAtM+lDQkJQ\nUFCAv/71r5Kdj0gqaphZ7Wz/+bbFOQJ3AAC+eWsvQu78vdXz2ruFanv3cNdqgSdSAqcF/vXXX0dj\nYyNWrFiB7OxsnDhxAhs2bJA0qHXr1mHBggW4du0ali1bZjWDnkgttHCjkuz3shCUbF2c/zB3DA5t\nMaFP+E1Wxx1NzpN6jbiSJzJeT02xkvo5LfCRkZGWf69YsULSYFrdddddqKio8Mq5iKTirTufSVU0\nivcV4eT50+iDQTaP6fS2vROOeiak7MlQ01p+NcVK2uC0wJ8/fx7btm3Dnj17cPlyyziJTqfDli3c\ngpHImY7egtUZsUXDnYuAzTs2o1OvALuPXTnTYPV1ez0TUvZkyNX9724+OVRB3uS0wD/++ONoaGjA\nhAkT4O/fMmWI94MnUgYxRcPdlqNZ14x+w4Lxw4clVhPqfvj7V3jgnvtRWXhYVM+ElD0ZcmwR25F8\ncjtb8ianBd5kMuHQoUPQ2+mSIyJ5iSka7rYcDYLeMs5+aIsJOr0eQnMz+nftg+ULH3cpTql6MuSY\nyNiRfNqjpEmXpC1Oq7bRaERlZaU3YiEiF4kpGmZds93nOGs5piSmoCb/R/QJvwkRyTEYOnMk+na5\nEf/fAteKu5RaY2yrJr8KyZOSJTtnR/PZltSxkm9z2oL/5ZdfYDQaERcXh169egHgGDyRUogZ33a3\n5Wivaz0zNQMJ8QmKWbfsrYmMbXkyn1LHSr7N6c1m3nnnHdsX6XR45JFHpIrJKd5sRhxuIiGO2vNU\nvK/IqmgkT0q2mWDn6CLA1eKi9lx5gph8Mk/iME/iubPRjSLvJucMC7w4/OURxxfy5OwiQCyl58pb\n68yd5VPpeVIK5kk8j+5kt27dOmRkZFhtG9tKp9PhpZdecj1CIpKF1Mv1lMAb68yvv4B4MPEBUe/N\nDW5IDg4LfJcuLVdW3bp1g06ngyAIVv8nUhv+kdU2R7Pbs9/P9tjGP+5cQHCDG5KLwwK/YMECAMBz\nzz3nrViIJNPeH9nWW3ySujlaMniy9hSK9xV1uJi6uzyOG9yQXJx26GdmZtq02H//+9/jvvvuQ1BQ\nkGSBEXlSe39kWeC1wdHs9k69AzxSTN3dqIYb3JBcnK6Dr6urQ15eHmpqalBTU4ONGzeioKAAMTEx\nePfdd70RI1GHubt2mdQjJTEF3//9K6tjP3xYgn6RwR75Pru7PI4b3JBcnBb4Y8eO4cCBA8jLy0Ne\nXh4OHDiA2tpafPXVV8jLy/NGjEQdxj+y2hcfNx79Am7EoS0mlG/9Jw5tMSFo+C3oE36TR77P7m5U\nww1uSC5Ou+hPnz6N/v37W74OCgpCTU0Nbr75Zpw7d07S4Ig8RQu3biXnljyyVLLvs7sb1XCDG5KL\n0wKfnJyM2bNn45FHHoEgCNi0aRNmzZqFq1evWm4+Q6R0/COrHe2thpD6++zuckNfWKZIyuN0o5uG\nhgZs3LgRn376KQRBQFJSElJTU9G5c2f88ssv6NOnj7diteBGN+JwEwlxmCfx5M6V/V3kfsTSmYsV\nVUDlzpNaME/icSc7ssJfHnGYJ/HkztXCFfNhSOhpc7y5sA7Zq3NkiMg+ufOkFsyTeB7dyW7Lli1I\nTk7G66+/brVMrnWjm8WLF7sfKRGRG7jkjEg8hwW+vLwcAFBaWsqd64hIEbgagkg8hwV+5cqVaG5u\nxqxZs3Dvvfd6MyYiIru4GoJIvHY79PV6PZ5++mkWeFIl7j2vPVwNQSSe0xH7UaNGYdu2bZg2bZo3\n4iHyCN7gQ7tcWXLGizzyZU4L/FdffYWcnBwEBQXh5ptvBtByu1iTySR5cETu4g0+iBd55OucFvh1\n69bZHOOkO1I6zrYmXuSRr3Na4MeOHQsAMJvNEAQBfn6urcMjkoOvzrZWWpe0nPHwIo98ndObzZw8\neRJpaWkICgpC//79kZ6ejlOnTnkjNiK3+eINPlq7pA0JPeF/d28YEnpiw9YsFO8r8sl4fPUij6iV\n0wK/bt06DBo0CIcPH0ZlZSUGDRqEtWvXeiM2IrfFx43H0pmL0VxYh8bC82gurNP8bOv2uqR9MR5f\nvMgjastpf3tBQQHKysqg17dcC/z5z39GVFSU5IERdZSv3eBDaV3ScsfDJXXk65wW+KFDh2L//v0Y\nM2YMAOB8amvwAAAXxElEQVTAgQMYMmSI5IERkWuU1iWthHh87SKPqC2nXfTLly/HnDlzYDQaYTQa\n8cgjjyAjI8MbsRGRC5TWJa20eIh8jdMWfGxsLI4ePWrZk37kyJEdPulHH32E5557DpWVlSgtLbXq\n8n/11Vfx2muvoVOnTvjb3/6G0aNHd/h8RFJTwux1b3ZJi/m87sSjhDwSaYXoNW/R0dEeO2lkZCS2\nbduGBQsWWB0/e/YssrKyUFRUhGPHjmHZsmX49ttvPXZeIikoaUMVb3RJFxQXiP68ru46p5Q8EmmB\nLIvaw8LC7B4vKSnBxIkTERISgpCQEAiCgPr6evTo0cPLERKJ52sbquRte6/d2fHutsB9LY9EUlPU\nrjUmkwnh4eGWr0NDQ2EymTB+vPUvt5+fAYGBXbwdnur4+bW0g5ir9nU0T7pOjt5Y0Fzu/fwMaNY3\nw95eludqzyLr42z0nTLQ0gLP+jgb3boFICE+wel7aymP/N0Th3kSrzVXLr1GgjgAAAkJCaipqbE5\nvnr1akyZMsXuawRBsDnGbXFJ6QyCAbY/udrdUMUgGNBs5/i58+cQkfJHq2N9pwzExk/esyrwBcUF\nyNv2Hsw6MwyCAbOnPYSE+ASfyyOR1CQr8AUFBS6/JjY2FoWFhZavKysr7Y79NzWZUVd3uUPx+YLW\nq2Lmqn0dzdOMCTMd3qNca7kPDOyC1KkP4sV319p83r59+9l9zZWma5Y8tB1n1wMQALz47lpcunRV\nU3nk7544zJN4gYFd4O/vWsmWvYu+bas9JiYGmZmZqK6uxtGjR6HX6zn+TrJwZTa3r22okhCfgEuX\nrtp83s07Ntt9ftsWeHvj7NmrcwD4Th6JpCZLgd+2bRuWLVuGn3/+Gffeey+MRiN27tyJfv36YdGi\nRRg3bhz8/f2Rk5MjR3ikQp5cXuXObG5f21DF0ed11AJv5Wx3O1/LI5GUZCnw06ZNw7Rp0+w+tnz5\ncixfvtzLEZGaeXp5FWdzu0dMT4YSdrcj8hWyd9ETdZSnC7Lce6irmbMWeEpiitNWPhF5Bgs8qZ6n\nCzJbmdLxtfkKRHJigSfV83RBZitTWhxnJ/IOpzebIVI6T9/UxBfvJU9E2sMWPKmeFN2+bGUSkdqx\nwJMmsCATEVljFz0REZEGscATERFpEAs8ERGRBrHAExERaRAn2RGRhSf39CciebHAE2mQO4Xa03v6\nE5G82EVPpDGthdqQ0BP+d/eGIaEnNmzNQvG+onZf196e/kSkPizwRBrjbqE265rtHudNdojUiQWe\nSGPcLdS8yQ6RtnAMnqiDlDYxzd1CzZvsEGkLCzxRByhxYpq7hZq3ciXSFp0gCILcQbiqsbEJdXWX\n5Q5D8QIDuwAAc+VER/K0cMV8GBJ62hxvLqxD9uqcDsfmruJ9RVaFOnlSskcKNX+mxGGexGGexAsM\n7AJ/f9fa5GzBE3WAWddst+Nb7olpvPkOEXGSHVEHcGIaESkVCzxRB6QkpqAm/0erYzX5VUielCxT\nRERELdhFT9QBnJhGRErFAk/UQRzvJiIlYoH3EqWtlSYiIm1jgfcCJa6VJiIibeMkOy/gTTyIiMjb\nWOC9gDfxICIib2MXvRdwrTSJwXkaRORJbMF7AddKkzPu3sOdiMgRWQp8ZmYmwsPDERUVhYyMDFy+\n/Ns+xK+++ioGDx6MIUOGYP/+/XKE53HxceOxdOZiNBfWobHwPJoL6zS3Vrp4XxEWrpiPeX+Zi4Ur\n5rMwuYjzNIjI02Tpop8wYQJefPFFAMCCBQvw/vvvIz09HWfPnkVWVhaKiopw7NgxLFu2DN9++60c\nIXqcltdKc5VAxyl1T3siUi9ZWvAJCQnQ6/XQ6/W45557sGfPHgBASUkJJk6ciJCQENx1110QBAH1\n9fVyhEguYOuz4zhPg4g8TfZJdm+++Sbmzp0LADCZTAgPD7c8FhoaCpPJhPHjrVuBfn4Gy20GyTE/\nv5biIHWudJ0cBSCo4vvkrTy1J23WbKzZuA59pwy0HDu7/SgyUzMUlUMl5EoNmCdxmCfxWnPl0msk\niANASyu9pqbG5vjq1asxZcoUAMCqVavQo0cPzJo1CwBg79b0Op1OqhDJQwyCAbbfObY+XZEQnwAA\n2PjJe5Y97TNTMyzHiYhcJVmBLygoaPfxd955B59//jmKin6bjBUbG4vCwkLL15WVlYiOjrZ5bVOT\nGXV1l22Ok7XWq2KpczVjwkzLGHyrmvwqLJ25pEPn9tayMW/lyZmYqNGIiRptdUzKmNzJr1JypXTM\nkzjMk3iBgV3g7+9ayZali/6zzz7DmjVrsHfvXnTu3NlyPCYmBpmZmaiursbRo0eh1+vRo0cPOUIk\nF0hxRzVO3JMW80ukfTrBXr+4xAYPHozGxkb07t0bAHDHHXcgK6vlj8v69evx2muvwd/fHzk5OYiL\ni7N5fWNjE6/4RFDz1fHCFfNhSOhpc7y5sA7Zq3M8ei4158ld7ubXF3PlDuZJHOZJPNW04I8cOeLw\nseXLl2P58uVejIaUSI3LxtS0E50a80tErpF9Fj2RPWpbNqa2Lm+15ZeIXMetakmR1La9r9r2AlBb\nfonIdWzBkyJJMXFPSmrr8lZbfonIdSzwpFhq2t5XjV3easovEbmOXfREHsAubyJSGrbgiTyAXd5E\npDQs8EQewi5vIlISdtETERFpEAs8ERGRBrHAExERaRALPBERkQaxwBMREWkQCzwREZEGscATERFp\nEAs8ERGRBrHAExERaRALPBERkQaxwBMREWkQCzwREZEGscATERFpEAs8ERGRBrHAExERaRALPBER\nkQaxwBMREWkQCzwREZEGscATERFpEAs8ERGRBrHAExERaZCf3AGoUfG+ImzesRlmXTMMgh4piSmI\njxsvd1hEREQWsrTgn3nmGQwfPhwjRoxAamoqamtrLY+9+uqrGDx4MIYMGYL9+/fLEV67ivcVYcPW\nLBgSesL/7t4wJPTEhq1ZKN5XJHdoREREFrIU+CeffBIHDx7Ed999h8GDB2P9+vUAgLNnzyIrKwtF\nRUXIzs7GsmXL5AivXZt3bEZQ0iCrY0FJg7Bl5xaZIiIiIrIlSxd9jx49AABNTU24dOkSAgMDAQAl\nJSWYOHEiQkJCEBISAkEQUF9fb3m+Eph1zTDYOd4Es9djISIickS2MfinnnoKOTk5CA0Nxe7duwEA\nJpMJ4eHhlueEhobCZDJh/Hjr8W0/PwMCA7t4M1yLAEMnCHaOd/brJFtMjvj5tVyKKC0upWGexGOu\nxGGexGGexGvNlSsk66JPSEhAZGSkzX/bt28HADz//POorq5GTEwMnnzySQCAINiWTp1OJ1WIbpk9\n7SGc3X7U6tjZ7UeROvUhmSIiIiKyJVkLvqCgwOlzunbtirS0NMybNw8AEBsbi8LCQsvjlZWViI6O\ntnldU5MZdXWXPResC2KiRmPxpavYsnMLmmCGHwxYPH0RYqJGyxaTI61XxUqLS2mYJ/GYK3GYJ3GY\nJ/ECA7vA39+1ki1LF/2RI0cwePBgNDU14YMPPsD06dMBADExMcjMzER1dTWOHj0KvV6vqPH3VvFx\n47ksjoiIFE2WAr9ixQocPnwYXbp0wdixYy0t+H79+mHRokUYN24c/P39kZOTI0d4REREqqcT7A18\nK1xjYxO7dERg95c4zJN4zJU4zJM4zJN47nTRc6taIiIiDWKBJyIi0iDuRU+KwP39iYg8iwWeZNe6\nv39Q0iDLLoEbtmYBAIs8EZGb2EVPsuP+/kREnscCT7Iz65rtHuf+/kRE7mOBJ9kZBPs/hn52b+tD\nRERisMCT7FISU1CT/6PVsZr8KiRPSpYpIiIi9eMkO5Jd60S6tvv7L525hBPsiIg6gAWeFIH7+xMR\neRa76ImIiDSIBZ6IiEiDWOCJiIg0iAWeiIhIg1jgiYiINIgFnoiISINY4ImIiDSIBZ6IiEiDWOCJ\niIg0iAWeiIhIg1jgiYiINIgFnoiISINY4ImIiDSIBZ6IiEiDWOCJiIg0iAWeiIhIg1jgiYiINIgF\nnoiISINY4ImIiDSIBZ6IiEiDZC3wr7zyCvR6Pc6fP2859uqrr2Lw4MEYMmQI9u/fL2N0RERE6uUn\n14lPnDiBgoICDBgwwHLs7NmzyMrKQlFREY4dO4Zly5bh22+/lStEIiIi1ZKtBf/444/jpZdesjpW\nUlKCiRMnIiQkBHfddRcEQUB9fb1MERIREamXLC34/Px8BAcHY9iwYVbHTSYTwsPDLV+HhobCZDJh\n/PjxVs/z8zMgMLCLV2JVMz8/AwAwV04wT+IxV+IwT+IwT+K15sql10gQBwAgISEBNTU1Nseff/55\nvPDCC9i1a5flmCAIVv9vS6fT2RzT63Xw95dtdEF1mCtxmCfxmCtxmCdxmCdpSJbVgoICu8cPHTqE\nY8eOYfjw4QCAkydP4g9/+ANKSkoQGxuLwsJCy3MrKysRHR0tVYhERESa5fXLpoiICPz000+Wr2+7\n7TZ888036N27N2JiYpCZmYnq6mocPXoUer0ePXr08HaIREREqid7v0jbLvh+/fph0aJFGDduHPz9\n/ZGTkyNjZEREROol+0Y3R48eRe/evS1fL1++HFVVVfi///s/xMXFWY4/88wzGD58OEaMGIHU1FTU\n1tZaHuPaeWuZmZkIDw9HVFQUMjIycPnyZctjzNVvPvroIwwdOhQGg8FmOSbzZGvv3r0IDw/H4MGD\n8dprr8kdjmKkpaWhX79+iIyMtByrr69HUlISQkJCMHXqVDQ0NMgYoTKcOHEC8fHxGDp0KMaOHYv3\n338fAHNlz5UrVxAbG4sRI0Zg1KhRWLt2LQA3ciWoxMWLFy3/XrlypfDMM88IgiAIP/30kxAaGir8\n+9//Fnbv3i0YjUa5QlSMXbt2CWazWTCbzcLcuXOFt956SxAE5up6FRUVwuHDh4WxY8cK33zzjeU4\n82TfiBEjhD179gjHjx8XQkNDhXPnzskdkiLs3btX+Pbbb4WIiAjLsRdffFFYunSpcOXKFWHJkiXC\nmjVrZIxQGc6cOSOUlZUJgiAI586dE2677Tbh4sWLzJUDly5dEgRBEK5cuSIMHTpU+Ne//uVyrmRv\nwYvVOhbf1NSES5cuoXPnzgC4dt6ehIQE6PV66PV63HPPPdizZw8A5up6YWFhuP32222OM0+26urq\nAABjxozBgAEDMGHCBJSUlMgclTLExcWhV69eVsdMJhPS09MREBCAtLQ05gpAUFAQRowYAQC48cYb\nMXToUJSWljJXDnTt2hUA0NDQgKamJgQEBLicK9UUeAB46qmnEBQUhP379yMzMxOA47Xz1OLNN9/E\nlClTADBXYjFPtkpLSxEWFmb5esiQIfj6669ljEjZ2uYrLCzM539+rldVVYXy8nLExMQwVw40Nzdj\n+PDh6NevH5YuXYqQkBCXc6WoAp+QkIDIyEib/7Zv3w6gZQ19dXU1YmJi8OSTTwIQv3Zea5zlCgBW\nrVqFHj16YNasWQB8M1di8nQ9X8wTeZa9nyFqUV9fj5SUFKxduxbdu3dnrhzQ6/U4ePAgqqqqkJWV\nhbKyMpdzJfss+rYcrZ1vq2vXrkhLS8O8efMAwGfXzjvL1TvvvIPPP/8cRUVFlmO+mCsxP1PX88U8\nORMdHW3pNQOA8vJyTJw4UcaIlC06OhoVFRUwGo2oqKjw+Z+fVteuXcOMGTOQmpqKpKQkAMyVM7fe\neisSExNRUlLicq4U1YJvz5EjRwC0jMF/8MEHmD59OgAgJiYGn3/+Oaqrq7F7926unQfw2WefYc2a\nNfj0008tcxUA5qo9ba+MmSdbgYGBAFpm0h8/fhwFBQWIjY2VOSrlio2NRW5uLi5fvozc3FyMGjVK\n7pBkJwgC0tPTERERgYyMDMtx5srWzz//jAsXLgAAamtrsWvXLiQlJbmeK8mmAHrYjBkzhIiICCE6\nOlrIzMwUzp8/b3ls3bp1wqBBg4Tw8HBh7969MkapDL///e+FkJAQYcSIEcKIESOERYsWWR5jrn7z\n8ccfC8HBwULnzp2Ffv36CRMnTrQ8xjzZ2r17txAWFiYMGjRIWL9+vdzhKMb9998v9O/fX/D39xeC\ng4OF3Nxc4eLFi8J9990n3HLLLUJSUpJQX18vd5iy27dvn6DT6YThw4db/jbt3LmTubLj+++/F4xG\nozBs2DBhwoQJwrvvvisIguByrnSCwAEQIiIirVFNFz0RERGJxwJPRESkQSzwREREGsQCT0REpEEs\n8EQa8txzz1mtV3dkzpw5eP31170Q0W/Gjh2L//3f/wUAzJs3D19++WW7z9+zZ4/VPganT5/GuHHj\nJI2RSEsUtdENEXWM2B33dDpdh3fnM5vNMBgMop/f9pxvvvmm0+cXFxfj0qVLSEhIAADcdNNN+OKL\nL9wLlsgHsQVPJBO9Xo+1a9ciKioK4eHhMJlMWLlyJYYPH47Jkyfj2LFjlud++OGHSEhIwKhRo/DK\nK69YbgF86dIlrFy5EkOGDEF8fDxOnTpleU1TUxPefvtt3Hvvvbjjjjvw7LPP4urVq5bHHa2Q1ev1\nePnllxEZGYkJEyZY3S639bGRI0di1apVuHTpEtavX4+7774bcXFxePnlly3PPXbsGObOnYuwsDCk\npaXh2rVrlsfatuYbGhrwyiuv4M4778SIESPw2GOP4dChQ8jJyUFeXh6MRiNeeuklHD9+HDfeeKPl\nPUpKSvDggw9i+PDheOKJJ3Dy5EkAwO7du2E0GrFs2TIMHToUM2bMsMolka9ggSeS0YULF/Dtt98i\nIyMDCQkJCA4OxsGDBxEREYG3334bAHD8+HE89dRTePvtt1FQUIB9+/bho48+AgDk5eXh+++/xz//\n+U/k5OTgk08+sbSSP/30U5SXl+PTTz/F/v37cebMGfzjH/8QFVd1dTW+++47PPfcc7j//vutHvvp\np5/w9ddfY+XKlcjJyYEgCCgsLERRUREKCwtx4MABAMCf//xn3H777aisrMSUKVMsxwHr1nxOTg7K\nysqwc+dOfPfdd1i5ciUiIiKwcOFCPPLIIygrK7Pce6Jtr8MDDzyAhQsXoqysDL169cKzzz5reezg\nwYOYMmUKysvLMWDAAOTm5rr0fSHSAhZ4Ihm13lNh7NixuHbtGmbPng0AGDdunOU2vx999BGSkpIQ\nEhKCHj164NFHH8WWLVsAAPn5+ZgzZw66du2K22+/Hffcc4/lvTdt2oQdO3Zg5MiRGDlyJL788kvR\ne/PPmzcPBoMBf/zjH9G3b1+ru1bNmzcPfn5+lnPk5ubCaDQiNjYWx48fx2effYZr166hoKAACxYs\nAABMmzYNN998s91zbdy4ERkZGfjd734HAOjdu7flMUe9DKWlpQgMDMSYMWOg1+uxaNEibN++HWaz\nGQBw2223Wbr2J06caMklkS/hGDyRjFqLWUBAALp164ZOnToBAPz9/XHlyhUALa3WtoVOEARLS/b6\ncfS2zzObzXjmmWfw0EMPeTTm4OBgq3NkZ2fjzjvvtHpOY2OjTTztjfl3dEPN1te3/r/tRYK/v79l\nSIPIl7AFT6RwM2fOxPbt23HixAk0NDQgLy8PycnJAICkpCTk5eXh119/RVVVFXbt2mV5XXp6OnJy\ncixj0/X19aisrBR1ztzcXDQ1NeGrr77C2bNnERMTY/d56enpWL9+PWprawG03CTj+PHj8Pf3x4QJ\nE/DWW28BaOlpaI3jeqmpqVi/fj3q6uoAwPJeAwcOxJkzZ+y+Jjo6GhcvXsT+/fvR3NyMnJwc3Hff\nfZaeBSJigSeSzfUt2rZftx2jvvXWW/H8888jLS0Nd999N+68807MnDkTQEtxjIiIwMiRIzF//nxM\nmzbN8h6TJ0/GnDlzkJ6ejmHDhiEuLs6qwLfXog4ODobRaMRf//pXbN682eFr5s6di7i4OEyfPh3D\nhg3DpEmTcPr0aQDAf//3f6OyshJhYWH45JNPbFr5rRYsWACj0YgJEyZgxIgR+K//+i8AwKRJk/Dz\nzz9bJtldP/P/gw8+QFZWFqKionD+/HmsWrXKJnf2vibyFbzZDBFZ0ev1aGhoQNeuXeUOhYg6gC14\nIrLC1i6RNrAFT0REpEFswRMREWkQCzwREZEGscATERFpEAs8ERGRBrHAExERaRALPBERkQb9/wVW\nfGQTYQxCAAAAAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Then we subract the model prediction from the measured data, square that difference and take the mean\n",
      "mse_ols = mean((y - modelprediction) ** 2)\n",
      "print \"Our mean squared error is %.4f using ordinary least squares regression\" % mse_ols"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Our mean squared error is 78.7821 using ordinary least squares regression\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now to convince you that this formula did in fact find the weights that\n",
      "best predict y as a weighted sum of the variables in X we can use a\n",
      "non-linear fitting proceedure to search for the parameters that minimize\n",
      "the squared difference between the model prediction and the data.  The\n",
      "fitting function we will use is lsqnonlin.  The following code will walk\n",
      "you through specifying a model to be fit with lsqnonlin.\n",
      "\n",
      "First we need to define our cost function.  lsqnonlin will search for the\n",
      "parameters that minimize the squared value of this function for each\n",
      "datapoint.  In this case the cost function we will use is y-X*w meaning\n",
      "that we want the algorithm to search for the weights, w, that minimize the\n",
      "difference between our data, y, and our regressors, X, multiplied by our\n",
      "weights.  Notice that we are defining this function in terms of the\n",
      "variable containing our regressors, X and the variable containing our\n",
      "dependent measure, y."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Various solvers live in the optimize package\n",
      "from scipy import optimize\n",
      "\n",
      "# Define an anonymous function that returns the mean squared error between our data\n",
      "# and the result of multiplying the predictor matrix by a weight vector\n",
      "costfunc = lambda w: mean(square(y - dot(X, w)))\n",
      "\n",
      "# Start the search with all weights at 0\n",
      "seed = zeros(5)\n",
      "\n",
      "# Fit the model by searching over a space of cost function values and finding the lowest\n",
      "# Set the limit on function evaluations to infinity\n",
      "optsol = optimize.fmin(costfunc, seed, maxfun=inf)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Optimization terminated successfully.\n",
        "         Current function value: 78.782052\n",
        "         Iterations: 633\n",
        "         Function evaluations: 1019\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that in the report we see the final function value is identical to the mean squared error from the OLS fit\n",
      "\n",
      "Plot the results to see the correspondence in the resulting weight vectors"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "axes(aspect=\"equal\")\n",
      "plot(w_ols, optsol, \"D\", color=colors[2], markersize=7)\n",
      "xlabel(\"OLS solution\")\n",
      "ylabel(\"Optimize solution\");"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAFeCAYAAACsM1VYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9UVGXiP/D3hWFAcxnDNLQohLSgBAYQWON83IMebf1V\npETowR+VUVmZyFl/rFqZa2WimbsR6X7LUnBJTwZoa1HCWipIYXYUc6VM2ZAQc1RiHGfmfv/wMDHO\nDDNO3suDvF/ncE7c586dt0/45vrMvTOSLMsyiIhIWD6dHYCIiDrGoiYiEhyLmohIcCxqIiLBsaiJ\niATHoiYiEpxiRf3dd99Br9fbvnQ6Hd544w27fcrLy6HT6Wz7LFu2TKk4RERdlkapA995552oqakB\nAFitVtxyyy1ITU112G/48OEoLi5WKgYRUZenytJHWVkZwsPDERIS4jDG+22IiDqmSlFv3rwZkydP\ndtguSRL27NmDmJgYZGdno66uTo04RERdiqT0LeQmkwm33HILDh8+jL59+9qNnT9/Hr6+vvDz88OG\nDRuwbds2lJaWOhzDapVhNluUjOmWRuMLAMzBHC6JkoU5xM3h4yN59VjFi/qjjz5CXl4e/v3vf3e4\nnyzLCA4OxokTJ+Dv7283ZjKZYTC0KhnTLZ2uBwAwB3O4JEoW5hA3h1br3cuCii99FBYWIiMjw+lY\nY2OjbY26pKQEUVFRDiVNRNTdKXbVBwC0tLSgrKwM69ats23Lz88HAGRlZWHLli3Iy8uDRqNBVFQU\ncnNzlYxDRNQlKb70cS1w6YM5RM8BiJOFOcTNIezSBxER/T4saiIiwbGoiYgEx6ImIhIci5qISHAs\naiIiwbGoiYgEx6ImIhIci5qISHAsaiIiwbGoiYgEx6ImIhIci5qISHAsaiIiwbGoiYgAGAxnsfCZ\nJ2AwnO3sKA5Y1ETU7RkMZ7Fk5gyMOt2MJTNnCFfWLGoi6tbaSjozoAcG3HADMgN6CFfWLGoi6rba\nl3SgVgsACNRqhStrFjURdUvOSrqNaGXNoiaibunVRfOR6qtxKOk2gVotUn01eHXRfJWTOWJRE1G3\nNG/ZK/jQYsY5k8np+DmTCR9azJi37BWVkzliURNRt6TT9cbSde/gfWOrQ1mfM5nwvrEVS9e9A52u\ndycl/A2Lmoi6LWdlLVpJAyxqIurm2pf1Ty0twpU0wKImIrKV9Sc39RGupAEFi/q7776DXq+3fel0\nOrzxxhsO+y1YsABhYWGIi4vDkSNHlIpDRNQhna43lq99S7iSBgCNUge+8847UVNTAwCwWq245ZZb\nkJqaardPVVUVdu/ejerqauzcuRM5OTkoLS1VKhIRUZekytJHWVkZwsPDERISYre9srISkyZNQlBQ\nEDIyMlBbW6tGHCKiLkWxM+r2Nm/ejMmTJztsr6qqQmZmpu37vn37oq6uDuHh4Xb7aTS+0Ol6KJ6z\nIxqNLwAwB3O4JEoW5hA7hzcUP6M2mUwoKSlBWlqaw5gsy5Bl2W6bJElKRyIi6lIUP6P++OOPERcX\nh759+zqMJSYm4vDhwxg9ejQAoKmpCWFhYQ77mc0WGAytSkftUNtvY+ZgDldEycIc4ubQar2rXMXP\nqAsLC5GRkeF0LDExEVu3bkVzczMKCgoQERGhdBwioi5H0TPqlpYWlJWVYd26dbZt+fn5AICsrCwk\nJCQgOTkZ8fHxCAoKwsaNG5WMQ0TUJUnylYvEAjKZzEL8swUQ459PzCFeDkCcLMwhbg5hlz6IiOj3\nYVETEQmORU1EJDgWNRGR4FjURESCY1ETEQmORU1EJDgWNRGR4FjURESCY1ETEQmORU1EJDgWNRGR\n4FjURESCY1ETEQmORU1EJDgWNRGR4FjURESCY1ETEQmORU1EJDgWNRGR4FjURESCY1ETEQmORU1E\nJDgWNRGR4FjURESCU7SoW1paMG3aNAwePBiRkZHYt2+f3Xh5eTl0Oh30ej30ej2WLVumZBwioi5J\no+TBn3/+edx2223Iz8+HRqNBS0uLwz7Dhw9HcXGxkjGIiLo0RYu6rKwMe/fuRUBAAABAp9M57CPL\nspIRiIi6PElWqCnr6+sxcuRIJCUloba2Fg8++CBmz55tK20AqKiowIMPPoiQkBCkpKRg1qxZCA8P\ndziW1SrDbLYoEdNjGo0vADAHc7gkShbmEDeHj4/k1WMVW6M2Go04evQoJk6ciPLychw6dAhFRUV2\n+8TGxuLkyZPYv38/IiMjMXv2bKXiEBF1WYqdUQNAREQEamtrAQAff/wx3nvvPRQWFjrdV5ZlBAcH\n48SJE/D397cbM5nMMBhalYrpEZ2uBwAwB3O4JEoW5hA3h1br3Wqzold9DBo0CJWVlbBardi+fTtG\njhxpN97Y2Ghboy4pKUFUVJRDSRMRdXeKvpi4cuVKTJ06FUajESNHjkR6ejry8/MBAFlZWdiyZQvy\n8vKg0WgQFRWF3NxcJeMQEXVJii59XCtc+mAO0XMA4mRhDnFzCLn0QUREvx+LmohIcCxqIiLBsaiJ\niATHoiYiEhyLmohIcCxqIiLBsaiJiATHoiYiEhyLmohIcCxqIiLBsaiJiATHoiYiEhyLmohIcCxq\nIiLBsaiJiATHoiYiEhyLmohIcCxqIiLBsaiJiATHoiYiEhyLmohIcCxqIiLBsaiJiATHoiYiEpyi\nRd3S0oJp06Zh8ODBiIyMxL59+xz2WbBgAcLCwhAXF4cjR44oGYeIqEtStKiff/553HbbbTh48CAO\nHjyIiIgIu/Gqqirs3r0b1dXVyMnJQU5OjpJxiIi6JEWLuqysDAsXLkRAQAA0Gg10Op3deGVlJSZN\nmoSgoCBkZGSgtrZWyThERF2SRqkD19fXw2g04sknn0RtbS0efPBBzJ49GwEBAbZ9qqqqkJmZafu+\nb9++qKurQ3h4uH1IjS90uh5KRfWIRuMLAMzBHC6JkoU5xM7hDcXOqI1GI44ePYqJEyeivLwchw4d\nQlFRkd0+sixDlmW7bZIkKRWJiKhLkuQrm/IaioiIsC1nfPzxx3jvvfdQWFhoG1+7di3MZjPmzJkD\nAAgPD0ddXZ3DcUwmMwyGVqVieqTttzFzMIcromRhDnFzaLXeLWIoukY9aNAgVFZWwmq1Yvv27Rg5\ncqTdeGJiIrZu3Yrm5mYUFBQ4vNhIREQKrlEDwMqVKzF16lQYjUaMHDkS6enpyM/PBwBkZWUhISEB\nycnJiI+PR1BQEDZu3KhkHCKiLsnt0kdNTQ1WrVqFiooKtLZe/qeDJEn4+eefVQkIcOmDOcTPAYiT\nhTnEzeHt0ofbR82bNw9jx47FokWL4O/v79WTEBGR99wW9alTpzB79mw1shARkRNuX0x86KGHsH79\nely6dEmNPEREdAW3Z9Svv/46zpw5g2effRY33HADAPXXqImIujO3RV1dXa1GDiIicsFtUYeGhgIA\nmpqaAFy+zZuIiNTjdo36yJEjSElJwT333IN77rkHI0aM4NuREhGpyG1Rr169GtOnT0dDQwMaGhow\nY8YMrF69Wo1sREQED4p67969mDp1Knx8fODj44MpU6Zg7969amQjIiJ4UNRJSUm2W7tlWUZBQQGS\nkpIUD0ZERJe5LernnnsO69evx80334zg4GCsX7/e9m53RESkPLdXfURGRqK8vBynTp0CAAQHByse\nioiIfuOyqC9evAh/f3/8+uuvAIDAwEAAsH3fs2dPFeIREZHLok5KSkJNTQ169erlMCZJEiwWi6LB\niIjoMpdFXVNTAwCwWq2qhSEiIkcevZjoyTYiIlKG26KuqKjwaBsRESnD5dLHBx98gKKiIhw/fhxp\naWm27Q0NDRg4cKAq4YiIqIOiHjx4MMaOHYvKykqMGzcObZ/YFR4ejmHDhqkWkIiou3NZ1NHR0YiO\njsaECRMQFBSkZiYiImrH7Q0vWVlZdt9LkgQAKCoqUiYRERHZcVvUY8eOhSRJkGUZZ86cQWFhIVJS\nUtTIRkRE8KCop0+fbvf9jBkzMGHCBKXyEBHRFdxennclPz8/NDc3K5GFiIiccHtG3f7SvIsXL2L/\n/v145plnFA1FRES/8WiNuk2PHj3w1ltvYcCAAR4dPDQ0FIGBgfD19YWfnx+qqqrsxsvLy3H//fcj\nLCwMADBx4kQsWrToavITEV33rnqN+mpIkoTy8vIOL+8bPnw4iouLvX4OIqLrncuibr/kcSVJkjy+\nPK/tRhlvx4mIujuXRd3+srwrtV1L7Y4kSUhJScHAgQPxyCOPOFwtIkkS9uzZg5iYGKSkpGDWrFkI\nDw93DKnxhU7Xw6PnVIpG4wsAzMEcLomShTnEzuENSVbwlLahoQH9+/dHbW0txo8fjy+++MLuE2LO\nnz9vW7/esGEDtm3bhtLSUofjWK0yzObOff/rtklmDuZwRZQszCFuDh8fz05yr+S2qH/55ResXbvW\nttSRnp6Op59+GjfeeONVPVF2djYiIiIwc+ZMp+OyLCM4OBgnTpyAv7+/3ZjJZIbB0HpVz3ettf02\nZg7mcEWULMwhbg6t1u3Lgk65vY56zZo1aGhowKZNm7Bx40Y0NDRgzZo1bg/866+/4vz58wCApqYm\n7Ny5E/fdd5/dPo2NjballZKSEkRFRTmUNJFIDIazWPjMEzAYznZ2FOpG3NZ7SUkJ9u7dC61WCwBY\nvXo1hg0bhhdeeKHDxzU2NiI1NRUA0KdPH8ydOxchISHIz88HcPk9RLZs2YK8vDxoNBpERUUhNzf3\nd/5xiJRjMJzFkpkzkOqrwZKZM7B03TvQ6Xp3dizqBtwW9a233opjx44hMjISAPD999/jlltucXvg\ngQMH4sCBAw7b27/J06xZszBr1qyryUvUKdpKOjOgBwK1WmSaTCxrUo3bon7qqacwYsQI/PGPf4Qs\ny9i3bx82bNigRjYiIVxZ0gAulzXAsiZVuC3q0aNHo66uDtu3b4ckSSgoKECPHp1/KRSRGpyVdBvH\nsubfC1KGR2/K1LNnT6SlpWHs2LH43//+p3QmImG8umg+Un01DiXdJlCrRaqvBq8umq9yMupO3BZ1\ncnIyDAYDTCYT7rrrLowcORIrVqxQIxtRp5u37BV8aDHjnMnkdPycyYQPLWbMW/aKysmoO3Fb1Bcu\nXIBOp8O2bdswYcIEHDt2jO/NQd2GTtcbS9e9g/eNrQ5lfc5kwvvGVq5Rk+LcFrVGo4HVakVxcTEe\neOABaDQaXLhwQY1sREJwVtYsaVKT26IeN24cwsLC8O2332LEiBH45ZdfoNF4d3cNUVfVvqx/amlh\nSZOq3Bb1Cy+8gPLycnz11VcAAIvFgnXr1ikejEg0bWX9yU19WNKkKo9OjUNDQ23/fdNNN+Gmm25S\nKg+R0HS63li+9q3OjkHdzFV/ZiIREamLRU1EJDgWNRGR4NwW9Y8//ojHH38cd955JwDg66+/dvvO\neUREdO24LeqXX34ZycnJtvf3iI6O9vjzEomI6PdzW9Rff/01pk6davucRF9fX15HTUSkIrdFfccd\nd6ClpcX2fU1NDYYMGaJoKCIi+o3bU+Pp06fjz3/+M5qamvDYY49hx44dfD9qIiIVuS3qUaNGISYm\nBps3b4bFYsErr7zCG16IiFTk0Wcmjh49Gs8++6xtW0VFBYYPH65oMCIiusztGvWkSZOQkpKCM2fO\n2LY999xzioYiIqLfuC3qyMhITJs2DcnJyTh27JgamYiIqB2P7kycOXMmVq9ejVGjRmHv3r1KZyLy\nmsFwFgufeQIGw9nOjkJ0zXh8C/no0aPx4YcfIjMzE8ePH1cwEpF32j6IdtTpZiyZOYNlTdcNt0X9\nt7/9zfbf0dHRqKio4Bo1Caf9p4UPuOEGZAb0YFnTdUOSZVnu7BDumExmGAytnZpBp7t8Cz1ziJfj\n7NmzeC79YWQG9LD7tHC1Py5LpDlhDjFzaLXe3dXt8ow6JSUFwOUPCujbt6/dV79+/Tw6eGhoKKKi\noqDX65GQkOB0nwULFiAsLAxxcXE4cuSIF38E6s7Onj2L7MkZDiUNAIFaLc+s6brgst43btwIAKiu\nrvb64JIkoby8HEFBQU7Hq6qqsHv3blRXV2Pnzp3IyclBaWmp189H3c8Lc7NxP3wdSrpNoFaL1EuX\n8Oqi+fxkFuqyXJ5RDxgwAMDls2JnX57qaGWlsrISkyZNQlBQEDIyMlBbW+t5ciIAL+Suwkew2D4d\n/ErnTCZ8aDFj3rJXVE5GdO24XTD5+uuvsWrVKlRUVMBoNAK4fKb8888/uz24JElISUnBwIED8cgj\nj2DChAl241VVVcjMzLR937dvX9TV1SE8PNw+pMbXts7UWTQaXwBgDuFy9MIbRUV49qGHMAVwWKPe\ndMmI1/+1Gb17K79GLc6cMIfIObx6rLsd5syZg4yMDCxevBj+/v5XdfAvv/wS/fv3R21tLcaPH4+E\nhAQEBwfbxmVZdjjjbns7VSJP9e7dG6sKCpE9OcNW1m0lvaqgUJWSJlKS26s+4uPjf9c6dZvs7GxE\nRERg5syZtm1r166F2WzGnDlzAADh4eGoq6tzeCyv+mAOT3K0XaKX6qvBhxazald7OMvSmZhD3BzX\n/KqPNo8//jhefPFF1NfX49dff7V9ufPrr7/i/PnzAICmpibs3LkT9913n90+iYmJ2Lp1K5qbm1FQ\nUICIiAiv/hBEAKDT9cbSde/gk5v6qF7SREpyW+/9+/dHdnY2XnzxRds2SZJgsVg6fFxjYyNSU1MB\nAH369MHcuXMREhKC/Px8AEBWVhYSEhKQnJyM+Ph4BAUF2a40IfKWTtebV3fQdcft0kdYWBjy8vLw\npz/96arXqK8VLn0wh+g5AHGyMIe4Obxd+nD7qH79+mHEiBH8nEQiok7itn0feOABZGVlYcaMGXY3\nrkRGRioajIiILnNb1G+99RYkScLnn39ut/2HH35QLBQREf3GbVHzLU2JiDqXy6K+ePEi/P39XV6K\n17NnT8VCERHRb1wWdVJSEmpqatCrVy+HMU8uzyMiomvDZVHX1NQAAKxWq2phiIjIkds7E519mgs/\n4YWISD1ui7qiosKjbUREpAyXSx8ffPABioqKcPz4caSlpdm2NzQ0YODAgaqEIyKiDop68ODBGDt2\nLCorKzFu3Djb25GGh4dj2LBhqgUkIuruXBZ1dHQ0oqOjMX78ePTp00fNTERE1I7bNWpZlrFgwQJE\nREQgIiICCxcuxOnTp9XIRkRE8KColy9fjgsXLuBf//oXNm/ejAsXLmD58uVqZCMiInhwC/mOHTvw\n7bffws/PDwCQm5uLIUOGYNWqVYqHIyIiD86o77nnHruP4vrqq69w9913KxqKiIh+4/aMuqWlBcnJ\nyYiOjoYkSThw4ABGjx6NtLQ0SJKEoqIiNXISEXVbbos6PT0d6enpdtskSYIsy/zEcCIiFbgt6unT\np6sQg4iIXOlwjfr777/H8uXLERUVhSFDhmD58uX4/vvv1cpGRETooKh37dqFe++9F42NjcjNzUVu\nbi5OnTqFYcOG4bPPPlMzIxFRt+Zy6eOvf/0rVq5ciSlTpti2jRo1ComJiVi8eDFGjBihSkAiou7O\n5Rn1qVOnkJGR4bA9PT0dDQ0NioYiIqLfuCxqrVaLQ4cOOWw/fPgwtFqtoqGIiOg3Los6PT0dOTk5\n+Oabb2zbDhw4gL/85S8Ol+sREZFyXBb1okWLkJiYiLS0NISEhCAkJARpaWkYOnQoFi9e7PETWCwW\n6PV6jB8/3mGsvLwcOp0Oer0eer0ey5Yt8+5PQUR0HXP5YqKfnx+WLl2KpUuX2m4hj4+Pv+onWLNm\nDSIjI3H+/Hmn48OHD0dxcfFVH5eIqLtw+14fwOWC9qak6+vrsWPHDjz22GO2Dx64kqvtRER0mUdF\n7a05c+bgtddeg4+P86eRJAl79uxBTEwMsrOzUVdXp2QcIqIuye0t5N4qLS1Fv379oNfrUV5e7nSf\n2NhYnDx5En5+ftiwYQNmz56N0tJSx5AaX+h0PZSK6hGNxhcAmIM5XBIlC3OIncMbkqzQ2sPChQvx\n/vvvQ6PRwGg04ty5c5g4cSLee+89p/vLsozg4GCcOHEC/v7+dmNWqwyz2aJETI+1TTJzMIcromRh\nDnFz+Ph490Z2ihV1exUVFVi5ciVKSkrstjc2NqJfv36QJAnFxcVYu3YtPv30U4fHm0xmGAytSsfs\nUNtvY+ZgDldEycIc4ubQar1bxFBs6eNKbW+Jmp+fDwDIysrCli1bkJeXB41Gg6ioKOTm5qoVh4io\ny1DljPr34hk1c4ieAxAnC3OIm8PbM2pFr/ogIqLfj0VNRCQ4FjURkeBY1EREgmNRExEJjkVNRCQ4\nFjURkeBY1EREgmNRExEJjkVNRCQ4FjURkeBY1EREgmNRExEJjkVNRCQ4FjURkeBY1EREgmNRExEJ\njkVNRCQ4FjURkeBY1EREgmNRExEJjkVNRCQ4FjURkeBY1EREgmNRExEJTvGitlgs0Ov1GD9+vNPx\nBQsWICwsDHFxcThy5IjScYiIuhzFi3rNmjWIjIyEJEkOY1VVVdi9ezeqq6uRk5ODnJwcpeMQEXU5\nihZ1fX09duzYgcceewyyLDuMV1ZWYtKkSQgKCkJGRgZqa2uVjENE1CUpWtRz5szBa6+9Bh8f509T\nVVWFyMhI2/d9+/ZFXV2dkpGIiLocjVIHLi0tRb9+/aDX61FeXu50H1mWHc60nS2RaDS+0Ol6KBHT\nYxqNLwAwB3O4JEoW5hA7hzcUO6Pes2cPiouLMXDgQGRkZODzzz/H1KlT7fZJTEzE4cOHbd83NTUh\nLCxMqUhERF2SJDtbPL7GKioqsHLlSpSUlNhtr6qqQnZ2Nj766CPs3LkTBQUFKC0tdXi8yWSGwdCq\ndMwOtf02Zg7mcEWULMwhbg6t1rtFDMWWPq7UtqSRn58PAMjKykJCQgKSk5MRHx+PoKAgbNy4Ua04\nRERdhipn1L8Xz6iZQ/QcgDhZmEPcHN6eUfPORCIiwbGoiYgEx6ImIhIci5qISHAsaiIiwbGoiYgE\nx6ImIhIci5qISHAsaiIiwbGoiYgEx6ImIhIci5qISHAsaiIiwbGoiYgEx6ImIhIci5qISHAsaiIi\nwbGoiYgEx6ImIhIci5qISHAsaiIiwbGoiYgEx6ImIhIci5qISHAsaiIiwSlW1EajEYmJiYiJiUFS\nUhJWr17tsE95eTl0Oh30ej30ej2WLVumVBwioi5Lo9SBAwICsGvXLvTs2RMXL15EXFwcxo8fjzvu\nuMNuv+HDh6O4uFipGEREXZ6iSx89e/YEAFy4cAFmsxn+/v4O+8iyrGQEIqIuT9GitlqtiI6Oxs03\n34ynn34aISEhduOSJGHPnj2IiYlBdnY26urqlIxDRNQlSbIKp7THjx/HmDFjsGnTJuj1etv28+fP\nw9fXF35+ftiwYQO2bduG0tJSh8dbrTLMZovSMTuk0fgCAHMwh0uiZGEOcXP4+EhePVaVqz5CQ0Mx\nZswYVFZW2m3/wx/+gJ49e8LPzw+PPvoo9u/fj4sXL6oRiYioy1DsxcTTp09Do9Ggd+/eaG5uxief\nfIK5c+fa7dPY2Ih+/fpBkiSUlJQgKirK6Tq22WyBwdCqVFSP6HQ9AIA5mMMlUbIwh7g5tFrvKlex\nom5oaMC0adNgsVgQHByMnJwc9O/fH/n5+QCArKwsbNmyBXl5edBoNIiKikJubq5ScYiIuixV1qh/\nL5PJLMRvQ0CM38rMIV4OQJwszCFuDm/PqHlnIhGR4FjURESCY1ETEQmORU1EJDgWNRGR4FjURESC\nY1ETEQmORU1EJDgWNRGR4FjURESCY1ETEQmORU1EJDgWNRGR4FjURESCY1ETEQnuuilqg+EsFj7z\nBAyGs50dhYjomrouitpgOIslM2dg1OlmLJk5g2VNRNeVLl/UbSWdGdADA264AZkBPVjWRHRd6dJF\n3b6kA7VaAECgVsuyJqLrSpctamcl3YZlTUTXky5b1K8umo9UX41DSbcJ1GqR6qvBq4vmq5yMiOja\n6rJFPW/ZK/jQYsY5k8np+DmTCR9azJi37BWVkxERXVtdtqh1ut5Yuu4dvG9sdSjrcyYT3je2Yum6\nd6DT9e6khERE10aXLWrAeVmzpInoetMlivrsWdc3s7Qv659aWljSRHTd6RJFnT05o8ObWdrK+pOb\n+rCkiei6o1hRG41GJCYmIiYmBklJSVi9erXT/RYsWICwsDDExcXhyJEjTveZ4hfg9mYWna43lq99\niyVNRNcdxYo6ICAAu3btwoEDB1BRUYF//vOfOHbsmN0+VVVV2L17N6qrq5GTk4OcnBynx+LNLETU\nnSm69NGzZ08AwIULF2A2m+Hv7283XllZiUmTJiEoKAgZGRmora11e0yWNRF1NxolD261WqHX63Ho\n0CG8/vrrCAkJsRuvqqpCZmam7fu+ffuirq4O4eHhHR43UKtF6qVLWPXiX/H6P/+fItmvpNH4AgB0\nuh6qPB9zdK0cgDhZmEPsHF499hrmcODj44NvvvkGx48fx5gxY3DvvfdCr9fbxmVZhizLdo+RJMnt\ncc+ZTCiRrHhj7RvQahX9IzhQ+/lcYQ57ouQAxMnCHPZEyeENVZKHhoZizJgxqKystCvqxMREHD58\nGKNHjwYANDU1ISwszOHx93601WHbn5WLS0QkFMXWqE+fPo2zZy+vITc3N+OTTz7B/fffb7dPYmIi\ntm7diubmZhQUFCAiIkKpOEREXZZiZ9QNDQ2YNm0aLBYLgoODkZOTg/79+yM/Px8AkJWVhYSEBCQn\nJyM+Ph5BQUHYuHGjUnGIiLouWQCtra1yQkKCHB0dLScmJsqrVq1yut/8+fPlgQMHyrGxsXJtbW2n\n5Ni1a5ccGBgox8TEyDExMfJLL710zXO0MZvNckxMjDxu3Din40rPhyc51JyP22+/XR4yZIgcExMj\nDx061Ok+asyJuxxqzcmFCxfkqVOnyoMGDZIjIiLkvXv3Ouyjxny4y6HGfBw5csR2/JiYGDkwMFBe\ns2aNw35Kz4cnObyZDyGKWpZluaWlRZZlWTYajfLdd98t//e//7Ubr6yslO+99165ublZLigokMeO\nHdspOXZz/wCZAAAHzklEQVTt2iWPHz9ekee+Um5urjx58mSnz6fWfLjLoeZ8hIaGys3NzS7H1ZoT\ndznUmpO5c+fKixYtkltbW+VLly7JZ8+etRtXaz7c5VDzZ0SWZdliscjBwcHyiRMn7Lar+Xemoxze\nzIcwt5Arcc21EjkAOFypooT6+nrs2LEDjz32mNPnU2s+3OUA1JkPT55LrTlxl8OT8WuhrKwMCxcu\nREBAADQaDXQ6nd24WvPhLgeg7s9IWVkZwsPDHS4HVvPno6McwNXPhzBFbbVaER0djZtvvhlPP/20\n02uuIyMjbd+3XXOtdg5JkrBnzx7ExMQgOztbkQwAMGfOHLz22mvw8XH+v0it+XCXQ635aHuulJQU\nPPDAAyguLnYYV2tO3OVQY07q6+thNBrx5JNPIjExEa+++iqMRqPdPmrMhyc51PwZAYDNmzdj8uTJ\nDtvV+vlwl8Ob+RCmqNuuuT527BjefPNN1NTU2I3LXl5zfa1zxMbG4uTJk9i/fz8iIyMxe/bsa56h\ntLQU/fr1g16v7/AsVun58CSHGvPR5ssvv8Q333yDl19+GdnZ2Th16pTduFo/I+5yqDEnRqMRR48e\nxcSJE1FeXo5Dhw6hqKjIbh815sOTHGr+jJhMJpSUlCAtLc1hTK2fD3c5vJoPL5dfFDV37lw5Ly/P\nbtsbb7xh9+JeWFhYp+Roz2q1yv369ZONRuM1fd4FCxbIt956qxwaGioHBwfLPXv2lDMzM+32UWM+\nPMnRnlLz4cycOXPkt99+225bZ/yMOMvRnpJzctddd9n+e8eOHfLDDz9sN67WfLjL0Z7SPyPbtm2T\nR48e7XRMzZ+PjnK05+l8CFHUTU1N8i+//CLLsiyfPn1aHjJkiPzTTz/Z7dP2QsDp06flTZs2KfJC\ngCc5Tp06JVutVlmWZfmjjz6SR44cec1ztFdeXu70ags15sOTHGrNR0tLi3zu3DlZlmX5559/liMj\nI12+WKTknHiSQ605GT9+vLxv3z7ZYrHIs2bNktevX283rtbPiLscav6dSU9Pl999912nY2r+neko\nhzfzIcQ9laJcc+1Jji1btiAvLw8ajQZRUVHIzc295jmu1PbPs86+Bt1ZDrXmo7GxEampqQCAPn36\nYO7cuQgJCVF9TjzJodacrFy5ElOnToXRaMTIkSORnp7eKT8j7nKoNR8tLS0oKyvDunXrbNs6Yz7c\n5fBmPiRZVvHlWCIiumrCvJhIRETOsaiJiATHoiYiEhyLmohIcCxq6nSXLl3C888/j7vuugthYWEY\nNmwYdu3aZRsvLy/H0KFDnT7273//O1JSUhAdHY27777b5YcoX62OnrO9H3/80e7VfQAYO3Ysfvjh\nh2uSgwhQ6YMDiDry0ksvobKyEh988AHuvvtuFBYWYvLkyfjyyy+dfpBEm2+//RZvvfUW9u3bh169\nesFkMil+e/KVfvjhB7z99tuYOXOmbdv27dtVzUDXP55RU6fbvHkzVqxYgSFDhsDHxwdTpkzBQw89\nhM2bN3f4uAMHDmDAgAG44YYbAABardblh0/k5eUhKSkJMTExiI2Nxblz5wBcfqOeyZMnIzo6Gjk5\nOaivr3d47JVn1+2/nzVrFg4fPgy9Xo+HHnoIwOVPNDp8+DAA4OTJk5g7dy5iYmIwZcoUVFVV2Y7j\n4+ODNWvWIC4uDsOGDcNnn33m6ZRRN8Oipk5VXV2N1tZWREdH220fN26c26JOTU2FwWDArbfeikcf\nfRQVFRVO92ttbcXixYvx6aef4sCBA9i9e7et3DMyMvDEE0+gpqYGN954I5YsWXJV+d98801ERkai\npqbG9h4X7d8/YvHixejTpw9qamrw+OOPIyMjw+7xZ86cQXV1NebNm4cXXnjhqp6bug8WNQlJdvIG\nOlfq1asXKisrUVJSgtDQUEyfPh2PPvqow349evRAXFwcpkyZgnfeeQcA4Ovri/3790On0+H//u//\n4OPjgyeffBIlJSWwWCxXldMVk8mE7du346mnnoIkSRg+fDgCAwPx1Vdf2faZOXMmJEnCfffdh337\n9uHSpUsePzd1Hyxq6lTx8fEICAjAgQMH7LaXlpY6nH26EhsbaztjLigowMWLFx322blzJxYuXIja\n2loMHjzY6RJHW+leWb4BAQF2xzxz5oxHua48rjNBQUEAAH9/f1gsFhY1OcWipk738MMPY/78+Th4\n8CDMZjM2bdqEoqIiPPzwwx0+7uuvv8aPP/4I4HIZ/uc//0FCQoLDhz1cuHABDQ0NSEpKwooVK3D7\n7bfj6NGjGDp0KM6dO4cvvvgCVqsV+fn5mDBhAjQa+9fYY2NjUV9fj4aGBrS0tODdd9+1jYWFhTm8\nzWkbrVaLcePGIT8/35bv/PnziIuL82KWqDtjUVOnW7JkCRISEpCWloZBgwbhH//4BzZt2mS74kOS\nJBw8eBAhISG2rxkzZuD06dNIT09HZGQk4uPjUVFRgddff93h+AaDAampqYiOjsYf//hHDBs2DCkp\nKQCAwsJCvPnmm4iNjcWZM2ewdOlS23O2rTVrtVqsWLECo0aNwogRIzB06FDb2G233YYxY8YgNjbW\n9mJiey+99BKampqg1+uRn5+PwsJC29iV74Ws1HsjU9fHN2UiIhIcz6iJiATHoiYiEhyLmohIcCxq\nIiLBsaiJiATHoiYiEtz/B1RAsuTNlvBQAAAAAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To verify the correspondence in mean squared error values, compute the MSE for the nonlinear solution"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mse_optim = mean(square(y - dot(X, optsol)))\n",
      "print \"OLS MSE: %.4f\" % mse_ols\n",
      "print \"Optimized MSE: %.4f \" % mse_optim"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "OLS MSE: 78.7821\n",
        "Optimized MSE: 78.7821 \n"
       ]
      }
     ],
     "prompt_number": 13
    }
   ],
   "metadata": {}
  }
 ]
}