lessons/Numerical computing in Python.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:9e0719c23b99e944a42b64d32e2ac53818f41e631768bd1a379bf126d9a8eef8"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Numerical computing in Python \n",
      "===="
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Set up inline plotting\n",
      "%matplotlib inline\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Numpy\n",
      "----"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import numpy.random as npr\n",
      "import scipy.linalg as la\n",
      "np.set_printoptions(precision=4)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Array indexing"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create a 10 by 6 array from normal deviates and convert to ints\n",
      "n, nrows, ncols = 100, 10, 6\n",
      "xs = np.random.normal(n, 15, size=(nrows, ncols)).astype('int')\n",
      "xs"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "array([[ 95,  73,  99, 114, 114,  91],\n",
        "       [108,  90, 114,  71,  87, 108],\n",
        "       [ 77, 108, 130, 113,  81, 101],\n",
        "       [126,  84, 103, 102, 106, 109],\n",
        "       [104, 111,  89, 111, 102, 102],\n",
        "       [105,  89,  64, 108,  90,  97],\n",
        "       [112,  86, 116,  97,  77, 118],\n",
        "       [ 97,  71,  88, 105, 117,  73],\n",
        "       [112, 104,  90, 100, 110, 143],\n",
        "       [ 86, 100,  78, 110,  87,  91]])"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Use slice notation\n",
      "print(xs[0,0])\n",
      "print(xs[-1,-1])\n",
      "print(xs[3,:])\n",
      "print(xs[:,0])\n",
      "print(xs[::2,::2])\n",
      "print(xs[2:5,2:5])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "95\n",
        "91\n",
        "[126  84 103 102 106 109]\n",
        "[ 95 108  77 126 104 105 112  97 112  86]\n",
        "[[ 95  99 114]\n",
        " [ 77 130  81]\n",
        " [104  89 102]\n",
        " [112 116  77]\n",
        " [112  90 110]]\n",
        "[[130 113  81]\n",
        " [103 102 106]\n",
        " [ 89 111 102]]\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#  Indexing with list of integers\n",
      "print(xs[0, [1,2,4,5]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 73  99 114  91]\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Boolean indexing\n",
      "print(xs[xs % 2 == 0])\n",
      "xs[xs % 2 == 0] = 0 # set even entries to zero\n",
      "print(xs)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[114 114 108  90 114 108 108 130 126  84 102 106 104 102 102  64 108  90\n",
        " 112  86 116 118  88 112 104  90 100 110  86 100  78 110]\n",
        "[[ 95  73  99   0   0  91]\n",
        " [  0   0   0  71  87   0]\n",
        " [ 77   0   0 113  81 101]\n",
        " [  0   0 103   0   0 109]\n",
        " [  0 111  89 111   0   0]\n",
        " [105  89   0   0   0  97]\n",
        " [  0   0   0  97  77   0]\n",
        " [ 97  71   0 105 117  73]\n",
        " [  0   0   0   0   0 143]\n",
        " [  0   0   0   0  87  91]]\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Broadcasting, row, column and matrix operations"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# operations across rows, cols or entire matrix\n",
      "print(xs.max())\n",
      "print(xs.max(axis=0)) # max of each col\n",
      "print(xs.max(axis=1)) # max of each row"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "143\n",
        "[105 111 103 113 117 143]\n",
        "[ 99  87 113 109 111 105  97 117 143  91]\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# A funcitonal rather than object-oriented approacha also wokrs\n",
      "print(np.max(xs, axis=0))\n",
      "print(np.max(xs, axis=0))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[105 111 103 113 117 143]\n",
        "[105 111 103 113 117 143]\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# broadcasting\n",
      "xs = np.arange(12).reshape(2,6)\n",
      "print(xs, '\\n')\n",
      "print(xs * 10, '\\n')\n",
      "\n",
      "# broadcasting just works when doing column-wise operations\n",
      "col_means = xs.mean(axis=0)\n",
      "print(col_means, '\\n')\n",
      "print(xs + col_means, '\\n')\n",
      "\n",
      "# but needs a little more work for row-wise operations\n",
      "row_means = xs.mean(axis=1)[:, np.newaxis]\n",
      "print(row_means)\n",
      "print(xs + row_means)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(array([[ 0,  1,  2,  3,  4,  5],\n",
        "       [ 6,  7,  8,  9, 10, 11]]), '\\n')\n",
        "(array([[  0,  10,  20,  30,  40,  50],\n",
        "       [ 60,  70,  80,  90, 100, 110]]), '\\n')\n",
        "(array([ 3.,  4.,  5.,  6.,  7.,  8.]), '\\n')\n",
        "(array([[  3.,   5.,   7.,   9.,  11.,  13.],\n",
        "       [  9.,  11.,  13.,  15.,  17.,  19.]]), '\\n')\n",
        "[[ 2.5]\n",
        " [ 8.5]]\n",
        "[[  2.5   3.5   4.5   5.5   6.5   7.5]\n",
        " [ 14.5  15.5  16.5  17.5  18.5  19.5]]\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# convert matrix to have zero mean and unit standard deviation using col summary statistics\n",
      "print((xs - xs.mean(axis=0))/xs.std(axis=0))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[-1. -1. -1. -1. -1. -1.]\n",
        " [ 1.  1.  1.  1.  1.  1.]]\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# convert matrix to have zero mean and unit standard deviation using row summary statistics\n",
      "print((xs - xs.mean(axis=1)[:, np.newaxis])/xs.std(axis=1)[:, np.newaxis])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[-1.4639 -0.8783 -0.2928  0.2928  0.8783  1.4639]\n",
        " [-1.4639 -0.8783 -0.2928  0.2928  0.8783  1.4639]]\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Linear algebra\n",
      "\n",
      "In general, the linear algebra functions can be found in scipy.linalg. You can also get access to BLAS and LAPACK function via scipy.linagl.blas and scipy.linalg.lapack."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = np.array([[1,2],[3,4]])\n",
      "b = np.array([1,4])\n",
      "print(A)\n",
      "print(b)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[1 2]\n",
        " [3 4]]\n",
        "[1 4]\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Matrix operations\n",
      "import numpy as np\n",
      "import scipy.linalg as la\n",
      "from functools import reduce\n",
      "\n",
      "A = np.array([[1,2],[3,4]])\n",
      "print(np.dot(A, A))\n",
      "print(A)\n",
      "print(la.inv(A))\n",
      "print(A.T)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 7 10]\n",
        " [15 22]]\n",
        "[[1 2]\n",
        " [3 4]]\n",
        "[[-2.   1. ]\n",
        " [ 1.5 -0.5]]\n",
        "[[1 3]\n",
        " [2 4]]\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = la.solve(A, b) # do not use x = dot(inv(A), b) as it is inefficient and numerically unstable\n",
      "print(x)\n",
      "print(np.dot(A, x) - b)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 2.  -0.5]\n",
        "[ 0.  0.]\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Matrix decompositions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = npr.normal(100, 15, (6, 10)).astype('int')\n",
      "print(A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 95  82 106 105 101  78 100  94 125  88]\n",
        " [109  99 101 112 105  99 103  87 108 102]\n",
        " [ 98  95 111 119  98  93 110  82 101 104]\n",
        " [106 110  86 123 100  88  95  80 106 114]\n",
        " [116 110  95  75 126  97  91 102  76 114]\n",
        " [123  99  80 118 113  94  94  82 106  89]]\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "P, L, U = la.lu(A)\n",
      "print(np.dot(P, A))\n",
      "print(np.dot(L, U))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[  98.   95.  111.  119.   98.   93.  110.   82.  101.  104.]\n",
        " [ 116.  110.   95.   75.  126.   97.   91.  102.   76.  114.]\n",
        " [ 123.   99.   80.  118.  113.   94.   94.   82.  106.   89.]\n",
        " [ 109.   99.  101.  112.  105.   99.  103.   87.  108.  102.]\n",
        " [ 106.  110.   86.  123.  100.   88.   95.   80.  106.  114.]\n",
        " [  95.   82.  106.  105.  101.   78.  100.   94.  125.   88.]]\n",
        "[[ 123.   99.   80.  118.  113.   94.   94.   82.  106.   89.]\n",
        " [ 106.  110.   86.  123.  100.   88.   95.   80.  106.  114.]\n",
        " [  95.   82.  106.  105.  101.   78.  100.   94.  125.   88.]\n",
        " [ 116.  110.   95.   75.  126.   97.   91.  102.   76.  114.]\n",
        " [ 109.   99.  101.  112.  105.   99.  103.   87.  108.  102.]\n",
        " [  98.   95.  111.  119.   98.   93.  110.   82.  101.  104.]]\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Q, R = la.qr(A)\n",
      "print(A)\n",
      "print(np.dot(Q, R).astype('int'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 95  82 106 105 101  78 100  94 125  88]\n",
        " [109  99 101 112 105  99 103  87 108 102]\n",
        " [ 98  95 111 119  98  93 110  82 101 104]\n",
        " [106 110  86 123 100  88  95  80 106 114]\n",
        " [116 110  95  75 126  97  91 102  76 114]\n",
        " [123  99  80 118 113  94  94  82 106  89]]\n",
        "[[ 95  82 106 105 101  78 100  94 125  88]\n",
        " [108  99 101 112 105  99 103  87 108 102]\n",
        " [ 98  95 111 119  98  93 110  82 101 104]\n",
        " [106 110  86 123 100  88  95  80 106 114]\n",
        " [116 110  95  75 126  97  91 102  76 114]\n",
        " [123  98  80 118 113  94  94  82 106  89]]\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U, s, V = la.svd(A)\n",
      "m, n = A.shape\n",
      "S =  np.zeros((m, n))\n",
      "for i, _s in enumerate(s):\n",
      "    S[i,i] = _s\n",
      "print(reduce(np.dot, [U, S, V]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[  95.   82.  106.  105.  101.   78.  100.   94.  125.   88.]\n",
        " [ 109.   99.  101.  112.  105.   99.  103.   87.  108.  102.]\n",
        " [  98.   95.  111.  119.   98.   93.  110.   82.  101.  104.]\n",
        " [ 106.  110.   86.  123.  100.   88.   95.   80.  106.  114.]\n",
        " [ 116.  110.   95.   75.  126.   97.   91.  102.   76.  114.]\n",
        " [ 123.   99.   80.  118.  113.   94.   94.   82.  106.   89.]]\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "B = np.cov(A)\n",
      "print(B)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 181.3778   42.7778   58.5111   23.0889 -121.9778   44.7556]\n",
        " [  42.7778   48.5      52.8333   68.7778  -28.8889   80.2222]\n",
        " [  58.5111   52.8333  110.3222   74.2444  -82.1333   36.1333]\n",
        " [  23.0889   68.7778   74.2444  183.9556  -24.1778  129.7333]\n",
        " [-121.9778  -28.8889  -82.1333  -24.1778  283.0667    9.3778]\n",
        " [  44.7556   80.2222   36.1333  129.7333    9.3778  219.5111]]\n"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "u, V = la.eig(B)\n",
      "print(np.dot(B, V))\n",
      "print(np.real(np.dot(V, np.diag(u))))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ -2.0102e+02  -8.7193e+01  -8.1063e+01   2.8856e+01  -1.2287e+01\n",
        "   -1.7071e-01]\n",
        " [ -1.2770e+02   4.9743e+01  -2.0215e+00   8.8173e+00   1.0113e+01\n",
        "    2.7248e+00]\n",
        " [ -1.6736e+02  -1.2602e+01   3.7575e+01   4.3426e+01   2.1998e+01\n",
        "   -9.9552e-01]\n",
        " [ -1.9952e+02   1.3677e+02   5.2859e+01   9.3565e+00  -2.5173e+01\n",
        "   -6.4364e-02]\n",
        " [  2.4073e+02   2.0714e+02  -3.1078e+01   3.7915e+01  -5.0583e-01\n",
        "   -6.5118e-02]\n",
        " [ -1.9546e+02   1.8346e+02  -3.9716e+01  -3.5476e+01   1.2268e+01\n",
        "   -7.6681e-01]]\n",
        "[[ -2.0102e+02  -8.7193e+01  -8.1063e+01   2.8856e+01  -1.2287e+01\n",
        "   -1.7071e-01]\n",
        " [ -1.2770e+02   4.9743e+01  -2.0215e+00   8.8173e+00   1.0113e+01\n",
        "    2.7248e+00]\n",
        " [ -1.6736e+02  -1.2602e+01   3.7575e+01   4.3426e+01   2.1998e+01\n",
        "   -9.9552e-01]\n",
        " [ -1.9952e+02   1.3677e+02   5.2859e+01   9.3565e+00  -2.5173e+01\n",
        "   -6.4364e-02]\n",
        " [  2.4073e+02   2.0714e+02  -3.1078e+01   3.7915e+01  -5.0583e-01\n",
        "   -6.5118e-02]\n",
        " [ -1.9546e+02   1.8346e+02  -3.9716e+01  -3.5476e+01   1.2268e+01\n",
        "   -7.6681e-01]]\n"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C = la.cholesky(B)\n",
      "print(np.dot(C.T, C)) # dot(C.T.conj(), C)\n",
      "print(B)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 181.3778   42.7778   58.5111   23.0889 -121.9778   44.7556]\n",
        " [  42.7778   48.5      52.8333   68.7778  -28.8889   80.2222]\n",
        " [  58.5111   52.8333  110.3222   74.2444  -82.1333   36.1333]\n",
        " [  23.0889   68.7778   74.2444  183.9556  -24.1778  129.7333]\n",
        " [-121.9778  -28.8889  -82.1333  -24.1778  283.0667    9.3778]\n",
        " [  44.7556   80.2222   36.1333  129.7333    9.3778  219.5111]]\n",
        "[[ 181.3778   42.7778   58.5111   23.0889 -121.9778   44.7556]\n",
        " [  42.7778   48.5      52.8333   68.7778  -28.8889   80.2222]\n",
        " [  58.5111   52.8333  110.3222   74.2444  -82.1333   36.1333]\n",
        " [  23.0889   68.7778   74.2444  183.9556  -24.1778  129.7333]\n",
        " [-121.9778  -28.8889  -82.1333  -24.1778  283.0667    9.3778]\n",
        " [  44.7556   80.2222   36.1333  129.7333    9.3778  219.5111]]\n"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Least squares solution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Suppose we want to solve a system of noisy linear equations\n",
      "\n",
      "$$\n",
      "y_1 = b_1 x_1 + b_0 \\\\\n",
      "y_2 = b_1 x_2 + b_0 \\\\\n",
      "y_3 = b_1 x_2 + b_0 \\\\\n",
      "y_4 = b_1 x_4 + b_0 \\\\\n",
      "$$\n",
      "\n",
      "Since the system is noisy (implies full rank) and overdetermined, we cannot find an exact solution. Instead, we will look for the least squares solution. First we can rewrrite in matrix notation $Y = AB$, treating $b_0$ as the coefficient of $x_0 = 1$:\n",
      "\n",
      "$$\n",
      "\\left(\n",
      "\\begin{array}{c}\n",
      "y_1 \\\\\n",
      "y_2 \\\\\n",
      "y_3 \\\\\n",
      "y_4 \n",
      "\\end{array} \n",
      "\\right) = \\left( \\begin{array}{cc}\n",
      "1 & x_1 \\\\\n",
      "1 & x_2 \\\\\n",
      "1 & x_3 \\\\\n",
      "1 &x_4 \\end{array} \\right)\n",
      "\\left( \n",
      "\\begin{array}{c}\n",
      "b_0 \\\\\n",
      "b_1 \n",
      "\\end{array}\n",
      "\\right)\n",
      "$$\n",
      "\n",
      "The solution of this (i.e. the $B$ matrix) is found by multipling the psudoinverse of $A$ (the Vandermonde matrix) with $Y$ \n",
      "\n",
      "$$\n",
      "(A^\\text{T}A)^{-1}A^\\text{T} Y\n",
      "$$\n",
      "\n",
      "This comes from setting the first derivative of $(Y - X \\beta)^2$ to 0 and solving for $\\beta$.\n",
      "\n",
      "Note that higher order polynomials have the same structure and can be solved in the same way and are also linear equations in $\\beta$.\n",
      "\n",
      "$$\n",
      "\\left(\n",
      "\\begin{array}{c}\n",
      "y_1 \\\\\n",
      "y_2 \\\\\n",
      "y_3 \\\\\n",
      "y_4 \n",
      "\\end{array} \n",
      "\\right) = \\left( \\begin{array}{ccc}\n",
      "1 & x_1 & x_1^2 \\\\\n",
      "1 & x_2 & x_2^2 \\\\\n",
      "1 & x_3 & x_3^2 \\\\\n",
      "1 & x_4 & x_4^2 \\end{array} \\right)\n",
      "\\left( \n",
      "\\begin{array}{c}\n",
      "b_0 \\\\\n",
      "b_1 \\\\\n",
      "b_2\n",
      "\\end{array}\n",
      "\\right)\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Set up a system of 11 linear equations\n",
      "x = np.linspace(1,2,11)\n",
      "y = 6*x - 2 + npr.normal(0, 0.3, len(x))\n",
      "\n",
      "# Form the VanderMonde matrix\n",
      "A = np.vstack([x, np.ones(len(x))]).T\n",
      "\n",
      "# The linear algebra librayr has a lstsq() function \n",
      "# that will do the above calculaitons for us\n",
      "\n",
      "b, resids, rank, sv = la.lstsq(A, y)\n",
      "\n",
      "# Check against pseudoinverse and the normal equation\n",
      "print(\"lstsq solution\".ljust(30), b)\n",
      "print(\"pseudoinverse solution\".ljust(30), np.dot(la.pinv(A), y))\n",
      "print(\"normal euqation solution\".ljust(30), np.dot(np.dot(la.inv(np.dot(A.T, A)), A.T), y))\n",
      "\n",
      "# Now plot the solution\n",
      "xi = np.linspace(1,2,11)\n",
      "yi = b[0]*xi + b[1]\n",
      "\n",
      "plt.plot(x, y, 'o')\n",
      "plt.plot(xi, yi, 'r-');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "('lstsq solution                ', array([ 6.0064, -2.0662]))\n",
        "('pseudoinverse solution        ', array([ 6.0064, -2.0662]))\n",
        "('normal euqation solution      ', array([ 6.0064, -2.0662]))\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGh9JREFUeJzt3Xl8VOW9x/EPgpIq1wW5ta6lr7QqKKhY1yodqwKKdb1a\nq9W6ttclwQWLBdG0ivZWrUpc6lKtK72tV62aytLqiKLFDdndRq0booiyiElJmPvHEyBAIMnMmTmz\nfN6vV15MJiczvyPk65Pfec7zgCRJkiRJkiRJkiRJkiRJUt7cCcwFprd47lhgJtAE9IujKElSx+wP\n7MaqYb4jsD3wFIa5JBWELm18/Rmg52rPvZabUiRJmVov7gIkSdkzzCWpBLTVZslYZWVlOpVK5erl\nJalUpYBvd/Sbsh2Zd1rbF1KpFOl0umQ/Lrvssthr8Pw8N8+v9D6AykzCuK0wHwM8B+wAvA+cBhzZ\n/HhvoA54IpM3liRFp602y4/X8vwjURciScqcF0AzlEgk4i4hp0r5/Er53MDzK1dr7XlHIN3c/5Ek\ntVOnTp0gg2x2ZC5JJcAwl6QSYJhLUgkwzCWpBBjmklQCDHNJKgGGuSSVAMNckkqAYS5JJcAwl6QS\nYJhLUgkwzCWpBBjmklQCDHNJKgGGuSSVAMNckkqAYS5JJaCtML8TmAtMb/Fcd2AC8AYwHtg0N6VJ\nktqrra2J9gcWA/cAfZqf+y0wr/nPYcBmwMWtfK/bxkmKTF3dREaPHk9DQxe6dm2kunoAgwf3j7us\nyGW6bVyXNr7+DNBztecOB77f/PhuIEnrYS5Jkairm8iQIeNIpUateC6VGgFQkoGeiUx65lsQWi80\n/7lFdOVI0ppGjx6/SpADpFKjqK2dEFNFhSfbC6Dp5g9JypmGhtabCPX1nfNcSeFqq83SmrnAN4CP\ngS2BT9Z2YE1NzYrHiUSCRCKRwdtJKndduzYCkOApurGYx/khABUVTXGWFYlkMkkymcz6ddrTZO8J\nPMaqF0A/A/6H0CvfFC+ASsqhv9/zEEvOGckuixdzFrfwBIdSWTmcG24YVHI980wvgLb1DWMIFzt7\nEEbklwJ/Bf4MbAe8CxwHfNHK9xrmkrKzdCmMHg1XXcWbBx3K0M+2YsHSCioqmqiqOrjkghxyF+bZ\nMMwlZS6ZhHPOgW23hdpa+M534q4oL3I1NVGS8uujj2DoUJg0Ca6/Ho48EjrlctxZGrydX1JhWLoU\nfvc76NsXvvUtmDULjjrKIG8nR+aS4vf006GlsvXW8NxzsP32cVdUdAxzSfGZMye0VJ59Fq67zpF4\nFmyzSMq/pUtDePftC9/8ZmipHH20QZ4FR+aS8mvixNBS2XLLMCLfYYe4KyoJhrmk/JgzBy66KIT5\nddc5Eo+YbRZJudXYGKYY9ukT5ozPng3HHGOQR8yRuaTcadlSmTTJlkoOGeaSMrLOzSI+/ji0VJ5+\nOswddySec4a5pA5b22YRnZoaOfTtaTBqFJxxRpil0q1bjJWWD8NcUoe1tlnElqlB9Drxv2DvvvDM\nM7DjjjFVV568ACqpw1puFrEFH3M3J/MAJ/DH7frD+PEGeQwMc0kd1rVrI51ppIrRTKcPH7EVvZjN\n5G13tDceE9sskjrssgO/wfXJrfnw333oz0ReoxeVlcOpqhoUd2lly/XMJbXf3Lnwi1/Ak0/yyomn\nM3xKI/UNXUp6s4h8cz1zSbnT2Ag33wyXXw6nnQazZ9OvWzfGxl2XVjDMJa3bs8+GG3969Ag3AfXq\nFXdFaoVhLql1c+fCsGHwj3/AtdfCscd6cbOAZTObZQgwHZjR/FhSKWhsDHtu7rwzfP3r4caf444z\nyAtcpiPznYEzgD2ApcBY4HEgFVFdkuIwaRKcfTZsvnm4Fb9377grUjtlOjLfEZgM1ANNwNPA0VEV\nJSnP5s6FU06BH/0Ihg8PrRWDvKhkGuYzgP2B7sCGwGBgm6iKkpQnjY1w440rWyqzZ4dAt6VSdDJt\ns7wG/A8wHvgSmAIsW/2gmpqaFY8TiQSJRCLDt5MUuUmTwiyV7t1tqcQomUySTCazfp2o/vd7JfAe\n8PsWz3nTkFSIPvkkzFKZMAGuucaReIHJ9KahbGazfL35z+2Ao4AHsngtSbnWsqXSo0doqRx/vEFe\nIrKZZ/4gsDlhNsvZwMJIKpIUveeeCy2VTTeFp56CnXaKuyJFzLVZpFL2ySdw8cVhWdqrr3YkXgTi\naLNIKlRNTXDTTWEE3r17uPHnxz82yEuYt/NLpeb550NLZeONQ0tl553jrkh5YJhLpeLTT0NLZezY\nMEvFlkpZsc0iFbumprA87U47hQucs2fbUilDjsylYtaypfLkk7ZUyphhLhWjli2Vq692JC7bLFJR\naWqCW24JLZVNNgktlRNOMMjlyFwqGpMnh+Vpu3WzpaI1ODKXCt2nn8IZZ8BRR8EFF0AyaZBrDY7M\npULV1AS33QaXXQYnnhhaKptsstbD6+omMnr0eBoautC1ayPV1QMYPLh/HgtWnAxzqRAtb6lstFHY\nKKJPn3UeXlc3kSFDxpFKjVrxXCo1AsBALxO2WaRCMm8enHlmaKmcf35YZ7yNIAcYPXr8KkEOkEqN\norZ2Qq4qVYExzKVC0NQEt94aZql06xZaKj/5SbtnqTQ0tP5Ldn195yirVAGzzSLF7YUXwo0/FRVh\nw4i+fTv8El27Nrb6fEVFU7bVqUg4MpfiMm8e/OxncOSRUF0NEydmFOQA1dUDqKwcscpzlZXDqao6\nOIpKVQQcmUv51tQEd9wBI0eGG37amKXSHssvctbWjqS+vjMVFU1UVQ3y4mcZcXMKKZ9efDHMUuna\nNaw3vssucVekAuPmFFIh++wz+PnP4YgjoKoKnnnGIFekDHMpl5Ytg9tvh969wwXOWbPg5JNdS0WR\ny6Zn/kvgJ8AyYDpwKtAQRVFSSXjxxTBLZf31wx6cjsSVQ5mOzHsCZwL9gD5AZ+D4iGqSitvylsrh\nh4cwt6WiPMg0zBcCS4ENCaP7DYEPoypKKkotWypdu4ZZKj/9KaxnN1O5l2mbZT5wLfAe8BUwDvh7\nVEVJReell8IovEsXGDcOdt017opUZjIN80rgPEK7ZQHwF+BE4P6WB9XU1Kx4nEgkSCQSGb6dVKDm\nz4cRI+CRR+A3v4GTTloxEncVQ7VHMpkkmUxm/TqZXlL/EXAwcEbz5ycBewPntDjGeeYqXcuWwZ13\nhiA/7ji4/PKwmXKz1lYxrKwcwQ03DDTQtU75nmf+GiG8v9b8pgcBszJ8Lam4vPwy7LNPCPNx46C2\ndpUgB1cxVP5lGuZTgXuAl4Bpzc/dFklFUqGaPx/OOgsGDw5/PvvsWnvjrmKofMvmMvtvgZ0IUxN/\nSpjdIpWeZcvgD38Is1Q6dw6zVE45ZZ2zVFzFUPnmQlvSurzySlhLZb314IknYLfd2vVt1dUDSKVG\nrNYzH05V1aBcVaoy50JbUmvmz4dLLoGHHoKrrspovnhd3URqaye0WMXwYC9+qk2ZXgA1zBWrfE7f\na9d7LVsGd90VZqkccwxccQVstllO6pFak2mY22ZRbPK5CXG73mt5S6VTJ/jb36Bfv0hrkHLJ+4wV\nm3xO31vne33+ebh789BDw84/kyYZ5Co6jswVm3xO32vtvTqxjMS706BXLzj66LA8bffukb+3lA+G\nuWKTz+l7q7/XrkzhJs5hs4/fg3/Uwe67R/6eUj7ZZlFs8rkJ8fL32pTPqeVcxjKIx/5zM96+9z6D\nXCXB2SyKVd6m7y1bxtQLhrPtLTfyzOY7cu+OCU698HCnCqrgODVRWpspU8IFzqamsInyd78bd0XS\nWrmhs7S6L74Imycfcgicfjo8/7xBrpJlmKv0LFsGf/xjmKXS2BhmqZx+ujv+qKQ5m0Wl5dVXQ0tl\n6VJ47DFH4iobDlVUGpa3VAYOhFNPhX/+0yBXWTHMVdyWLYO77w4tlaVLQ0vljDNsqajs2GZR8Zo6\nNbRU/v1vWyoqew5fVHy++AKqq2HAgLA0rS0VyTBXEUmn4Z57QkuloSG0VM4805aKhG0WFYupU+Hc\nc6G+Hh59FPbYI+6KpIKS6ZBmB2BKi48FQHVURUkrfPEFDBkSWionnRRaKga5tIZMw/x1YLfmj92B\nJcDDURUlrdJSqa8PLZWf/SxsqCxpDVG0WQ4CUsD7EbyWBNOmhVkq9fXw17/CnnvGXZFU8KK4cnQ8\n8EAEr6Nyt2ABnHceHHzwypaKQS61S7Yj8w2AHwLDWvtiTU3NiseJRIJEIpHl26kkpdNw330wbBgc\ndhjMnAk9esRdlZQXyWSSZDKZ9etkuwTuEcBZwKBWvuYSuGrbtGlhlsqSJWF52r32irsiKVZxLYH7\nY2BMlq+hcrS8pXLQQXDCCTB5skEuZSGbMN+IcPHzoYhqUTlY3lLp1Qu+/DLMUvnv/3aWipSlbHrm\nXwI2NtV+06eHWSpLlsDDDzsSlyLkfdDKvQUL4Pzz4cADbalIOWKYK3fSabj//tBSWbQozFKxpSLl\nhGuzKDdmzAgtlUWL4KGHYO+9465IKmmOzBWthQvhggvgBz+A44+HF180yKU8MMwVjXQaHnggtFQW\nLgwtlbPOsqUi5YltFmVvxoxw48/ChfDgg7DPPnFXJJUdR+bK3MKFcOGFoaVy3HGhpWKQS7EwzNVx\nLVsqn38eRuZnn21LRYqRbRZ1zMyZYZbKggW2VKQC4shc7bNoEQwdCgccAMceCy+9ZJBLBcQw17ql\n0zBmTGipzJ+/cv64LRWpoNhm0drNnBlmqXzxBfz5z7DvvnFXJGktHJlrTYsWwUUXQSIBxxwTZqkY\n5FJBM8y1UjoNf/pTaKnMm7dy/ngXf4GTCp0/pQpmzQrBPX8+/O//wve+F3dFkjrAkXm5W95S+f73\n4eijwywVg1wqOoZ5uUqnwwi8Vy/49FNbKlKR8ye3HM2eHYJ73rzQI99vvxVfqqubyOjR42lo6ELX\nro1UVw9g8OD+MRYrqT0M83KyeDH8+tdw111w6aVhVcMWI/G6uokMGTKOVGrUiudSqREABrpU4LJp\ns2wKPAjMBmYBLlpdqNLpME+8Vy+YOze0VKqq1mipjB49fpUgB0ilRlFbOyGf1UrKQDYj8xuAvwH/\n1fw6G0VSkaLVsqUyZswqLZXVNTS0/s+hvt67PaVCl+nIfBNgf+DO5s8bgQWRVKRoLF4Mw4ZB//5w\nxBHw8svrDHKArl0bW32+oqIpFxVKilCmYf4t4FPgLuAV4HZgw6iKUhbSafjLX0JLZc4cmD4dqqvb\nNUulunoAlZUjVnmusnI4VVUH56paSRHJtM3SBegHnAu8CFwPXAxc2vKgmpqaFY8TiQSJRCLDt1O7\nvPZaaKl88klYb3z//Tv07csvctbWjqS+vjMVFU1UVQ3y4qeUQ8lkkmQymfXrdMrw+74BPE8YoQPs\nRwjzw1ock06n01mUpnZbvBiuuAL+8Ae45JKwqqHzxaWi1KlTJ8ggmzNts3wMvA9s3/z5QcDMDF9L\nmWrZUvnoo9BSGTLEIJfKUDY/9VXA/cAGQAo4NZKK1D6vvRamF86dC/ffHy50SipbmbZZ2sM2Sy7Y\nUpFKWqZtFlOgwK24vb6+Mwd+MYOhc57lawMHwLRpsOWWcZcnqUAY5gVs+e31nVMnU0sVWzKH07Ya\nwE+OP5PBBrmkFlw1sYDd+rvHOT21jGfZjyc4hH68wp8+utfb6yWtwZF5IUqn4aGHuH3SrUzgh/Rl\nGh+zciTu7fWSVmeYF5o33gizVD78kKv7HM61L927xiHeXi9pdbZZCsWXX8Lw4WHj5EGDYMoUDqg5\n09vrJbWLI/O4pdPw8MNw/vlhIaxp02CrrQBvr5fUfs4zj9Mbb4RFsD74AG68EVy7Rip7+b6dX9n4\n8ksYMSK0VAYMgClTDHJJWTHM86l5lgq9e8O774aWygUXwPrrx12ZpCJnz7yDMt7w+M03wyyV99+H\nu+92JC4pUoZ5B2S04fGSJXDllfD738Mvfxl65I7EJUXMNksHdGjD4+WzVHr3hrffhqlT4cILDXJJ\nOeHIvAPaveHxm2+GEfi//gV33QUHHJCH6iSVM0fmHdDmhsdLlsDIkbDPPnDggWE0bpBLygPDvAPW\nuuHxuQfBI4+Elspbb4UQHzrUloqkvPGmoQ6qq5tIbe2EFXdkDjumFwc8fF+YanjjjfCDH8RdoqQi\nlulNQ4Z5ppYsgd/8Bm6+GYYNC3tvbrBB3FVJKnLuNJQv6TQ8+iicdx7stRe8+ipss03cVUkqc9mE\n+bvAQqAJWArsGUVBBe2tt8II/O234Y47wkVOSSoA2VwATQMJYDdKPciXLIFLL4W99w53bk6dapBL\nKijZtlly2XOPXzoNjz0WRuN77mlLRVLByibM08DfCW2WW4HbI6moUKRSIcTfegtuvx0OOijuiiRp\nrbIJ8+8Bc4D/BCYArwHPtDygpqZmxeNEIkGiGBaX+uqrMEvlppvgF78Iqxw6S0VSjiSTSZLJZNav\nE1Wb5DJgMXBti+eKb2ri8pbKd78L114L224bd0WSyky+pyZuCHQGFgEbAQOAX2X4WvF7++0Q4m++\nCbfdZktFUtHJdDbLFoSWyqvAZOBxYHxUReXNV19BTU24uLl8/02DXFIRynRk/g6wa5SF5N3jj4eV\nDXffPWzbZktFUhErvztA33473L35+utw661w8MFxVyRJWSufVRO/+gp+9avQUtl339BSMcgllYjy\nGJk//ni4wNmvH7zyCmy3XdwVSVKkSjvM33knhPjrr8Mtt8CAAXFXJEk5UZptlvp6+PWvYY89wq4/\n06YZ5JJKWumNzOvqwmh8111tqUgqG6UT5u+8E2apzJ4dbsUfODDuiiQpb4q+zfLEwxO4+9vfZ8H2\nO3HXrIU88dubDXJJZaeow/yFmt/S+7gf8R+pzenbOJvT3nqKqqFPUVc3Me7SJCmvinMP0HffhfPO\n48MJEzl9yQOMY9AqXx44cCRjx16em/eWpBzKdKGt4hqZ19fDFVeEVQ332INTdj97jSAPh3WOoThJ\nik/xhPkTT8DOO4cZKi+/DCNGsN7XWj+0oqIpv7VJUswKP8zffReOOiosilVbGzaL+OY3AaiuHkBl\n5YhVDq+sHE5VlbfpSyovhTs1sb4errkGrr8ezj8fxoyBiopVDhk8uD8AtbUjqa/vTEVFE1VVg1Y8\nL0nlojAvgI4dC1VVoa1y3XXQs2ekhUlSocr3TkO58a9/hVH49OkwejQcckjcFUlSUSiMnnlDA4wa\nFTaK6NcvhLlBLkntFv/IfOzYcHGzd2946SVbKpKUgWzDvDPwEvAB8MMOfefylsq0aaGlcuihWZYi\nSeUr2zbLEGAW0P4rnQ0NcOWVoaWy224wY4ZBLklZyibMtwEOBe6gvVdex42DPn3ghRfgxRdh5Mg1\nphtKkjoumzbLdcBFwMZtHvnee6Gl8uqroaUyeHAWbytJWl2mI/PDgE+AKaxrVN7QAFddFWao7LIL\nzJxpkEtSDmQ6Mt8XOJzQZqkgjM7vAU5ueVDNVltBjx5w4okk+vcnYUtFklaRTCZJJpNZv04Ud4B+\nHxjKmrNZ0unHHoPDDovgLSSpPMS9BG7rs1kMcknKi8Jcm0WSylTcI/NWDRx4iVu4SVIe5PR2/vHj\nryCVCuuNuyytJOVOzhfaSqVGUVs7IddvI0llLS+rJronpyTlVl7C3D05JSm3ch7m7skpSbmX0wug\nAweOdE9OScoD55lLUgEpyHnmkqT8MMwlqQQY5pJUAgxzSSoBhrkklQDDXJJKgGEuSSXAMJekEmCY\nS1IJMMwlqQQY5pJUAjIN8wpgMvAqMAu4KrKKJEkdlmmY1wMHALsCfZsf7xdVUcUgmUzGXUJOlfL5\nlfK5gedXrrJpsyxp/nMDoDMwP/tyikep/4Mq5fMr5XMDz69cZRPm6xHaLHOBpwjtFklSDLIJ82WE\nNss2QH8gEUVBkqSOi2pzipHAV8A1LZ57C6iM6PUlqVykgG/n6816AJs2P/4aMBE4MF9vLklaVaZ7\ngG4J3E1o06wH3Av8I6qiJEmSJHXQnYTZLNPXccxo4E1gKrBbPoqKUFvndyLhvKYBkwhz7otJe/7+\nAPYAGoGjc15RdNpzbglgCjADSOa+pEi1dX49gLGEGWczgFPyU1ZktiXMkptJqL96LccVa7605/zy\nmi/7E/4Dru0f1KHA35of7wX8M5fF5EBb57cPsEnz40GU3vlBuIfgSeBx4Jh8FBWRts5tU8IP0jbN\nn/fIR1ERauv8alh5Z3YP4DMyb6vG4RuE2XIA3YDXgV6rHVPM+dKe8+tQvmS7NsszwOfr+PrhhN46\nhNv/NwW2yPI986mt83seWND8eDIrg6FYtHV+AFXAg8CnuS8nUm2d2wnA/wEfNH8+L+cVRaut85sD\nbNz8eGNCmDfmuqgIfUz4rQJgMTAb2Gq1Y4o5X9pzfh3Kl1wvtLU18H6Lzz+g+AKvvU5n5SihVGwN\nHAHc0vx5OsZaovYdoDvhV92XgJPiLSdytwM7AR8RflUfEm85WelJ+C1k8mrPl0q+9KT182upzXzJ\nx69dq89lL6VAWO4A4DTge3EXErHrgYsJf2ediO6+hEKwPtCPMKV2Q8Io6J+E/mspGE4Y+SUI93tM\nAHYBFsVYUya6EX4zHEIYwa6u2POlrfODduZLrsP8Q0Kjf7ltmp8rJX0Jo6BBtN2yKDa7A39qftwD\nOARYCjwaW0XReZ/QWvmq+WMiIexKJcz3BUY1P04B7wA7EH4LKRbrE1ph9wGPtPL1Ys+Xts4P8pwv\nPWnfBdC9Ka4LFMv1ZO3ntx3hTte981ZN9HrS9mwWgLsortkssO5z2xH4O+EC74bNx/XOT1mR6cna\nz+93wGXNj7cgtCC656GmqHQC7gGuW8cxxZwv7Tm/vObLGEJP7t+Ekc5pwM+bP5a7sbmgqYRfa4tJ\nW+d3B+HC0pTmjxdiqDEb7fn7W67Ywrw95zaUMKNlOmuf+lao2jq/HsBjhJ+76YQLvsVkP8L6T6+y\n8ufrEEonX9pzfsWeL5IkSZIkSZIkSZIkSZIkSZIkSZKy8f86OrRjB8TZaAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1090c4e90>"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# As advertised, this works for finding coeefficeints of a polynomial too\n",
      "\n",
      "x = np.linspace(0,2,11)\n",
      "y = 6*x*x + .5*x + 2 + npr.normal(0, 0.6, len(x))\n",
      "plt.plot(x, y, 'o')\n",
      "A = np.vstack([x*x, x, np.ones(len(x))]).T\n",
      "b = la.lstsq(A, y)[0]\n",
      "\n",
      "xi = np.linspace(0,2,11)\n",
      "yi = b[0]*xi*xi + b[1]*xi + b[2]\n",
      "plt.plot(xi, yi, 'r-');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG11JREFUeJzt3XucjnX+x/EXQzOhJRFy+E2NbKXTZBe/crgThkVSbWe1\nHbXVjKVW5VBTKh2pwY8Owq5SSsm4w6C5DaFNjRLlcK+2WqcOimIGM9fvj+89Y4xh7vN139f9fj4e\n12Puue7Tx+2az3zn8z2BiIiIiIiIiIiIiIiIiIiIiIiIiEjUpAAfAWuA9cAY3/mGwCJgI5AHNLAl\nOhER8Vsd39dawCqgE/A0MMx3/n7gSRviEhGRINQBPgbaAl8BTXznm/q+FxGRGFYTU2bZg2mRA+yq\ncH+NSt+LiEgMq48ps1zMkcn7p+iHIyIiZWoF8NhfADfQDtiBKa9sB5oBOys/OC0tzfJ6veGIUUQk\nkXiB1oE+qWY19zfi0EiV44EeQCEwF7jJd/4mYM4R0Xi9WJalI0zHww8/bHsMTjr0eerzjNUDSAs0\nkUP1LfNmwHRM0q8J/BNY4kvos4Bbga+Bq4J5cxERCY/qkvla4IIqzv8EdA9/OCIiEozqyiwSI1wu\nl90hOIo+z/DS52m/GhF8bctX/xERET/VqFEDgsjNapmLiDiAkrmIiAMomYuIOICSuYiIAyiZi4g4\ngJK5iIgDKJmLiDiAkrmIiAMomYuIOICSuYiIAyiZi4g4gJK5iIgDKJmLiMQAt7uAK7vdG/TzlcxF\nRGzmdhcwPHMuw/KXBf0aSuYiIjabPM7NC1tW82mVewH5R+uZi4jYaf9+VjY7G+9P7bmRf2CRBFrP\nXEQkjpSUwA03YCXV4GamYoWQkpXMRUTsUFoKt98Ou3bxy4sT+J+07JBerroNnUVEJNwsC4YMgQ0b\nIC+P3nXr8sJxyYwfP4qFC4N7SdXMRUSi7aGHIDcX8vOhQYPD7gp2D1C1zEVEoumZZ+Ctt6Cg4IhE\nHgolcxGRaHnxRZg0ySTyxo3D+tJK5iIi0fDaazB6NCxdCi1ahP3lVTMXEYm0996DQYNgyRJo2/aY\nD1XNXEQkFi1ebIYgzp9fbSIPhcaZi4hEyooVcN11MHs2tGsX0beqLpm3BPKBdcAXQJbvfDbwHVDo\nO3pFKD4RkfhUWAgDBsA//wmdO0f87aqryzT1HWuAesAnwGXAVcAeYOwxnquauYgkpi+/hG7dYOJE\nuPzygJ4aqZr5dt8B8CvwJdC87D0DfTMREcfbsgV69oSnngo4kYcikJp5KpAOrPJ9nwl8BkwBwjfy\nXUQkXm3dCt27wwMPwI03RvWt/U3m9YC3gcGYFvok4FTgfGAb8FxEohMRiRc//AA9epiRK3ffHfW3\n92doYm1gNjADmOM7t7PC/a8AuVU9MTs7u/y2y+XC5XIFE6OISGz75Rfo1Qv69zet8gB4PB48Hk/I\nIVRX964BTAd+BIZUON8M0yLHd/6PwHWVnqsOUBFxvr17ISMDzjsPxo+HGqF1JwbbAVrdEzoBBcDn\nQFlmHg5ciymxWMAWYBCwo9JzlcxFxNmKi+HSS6FpU5g6FWqGPnUnUsk8FErmIuJcBw/C1Veb22++\nCbXCM6Fe0/lFRKKltBRuvRV++82suxKmRB4K+yMQEYknlgVZWWY8+YIFkJxsd0SAkrmISGBGjIBV\nq8wKiHXq2B1NOSVzERF/jRljyipLl0L9+nZHcxglcxERf0ycCFOmmF2CGjWyO5ojKJmLiFRn+nR4\n8kmTyE85xe5oqqShiSIixzJ7NtxzD3zwAZx5ZsTfTuPMRUTCbeFCGDjQfE1Pj8pbapy5iEg4LVsG\nN9wAc+ZELZGHQtvGiYhUtno1XHEFvP46XHSR3dH4RS1zEUl4bncBOTl5FBfXos2BbeSsf5uUqVPM\nkrZxQslcRBKa213A4MEL8Xof5zS8zKArI07uRLfaDeljd3ABUJlFRBJaTk4eXu/jNOc7FtOdxxjJ\n2J3vMX78IrtDC4iSuYgktOLiWrTkGz6gG//HXbzInQAUFSXZHFlglMxFJKGdWrKTArowmTt5lr+X\nn09JKbExqsApmYtI4vriCyatf4uXG5/FOIaWn05LG05mZvx0foI6QEUkUa1eDX37kjL+BS6s34KM\n8aMoKkoiJaWEzMxe9OnTxe4IA6IZoCKSeJYtM+PIX3oJLrvM7mgOoxmgIiL+yMuD6683E4LiaBx5\ndVQzF5HEMWeOmaL/7ruOSuSgZC4iieK11+DOO2H+fOjUye5owk7JXESc76WXYNgws9Vbu3Z2RxMR\nqpmLiLONHQs5OWart9at7Y4mYpTMRcSZLAsefdR0dC5bBi1b2h1RRCmZi4jzWBb8/e9m5EpBATRp\nYndEEadkLiLOUloKd90FhYXg8UDDhnZHFBVK5iLiHAcPwl/+At9+C4sXwwkn2B1R1CiZi4gzFBfD\nNddAUZEZflinjt0RRZWGJopI/Nu7Fy69FGrWNBODEiyRQ/XJvCWQD6wDvgCyfOcbAouAjUAe0CBS\nAYqIHNPu3dCrl+nkfPNNSE62OyJbVJfMDwBDgLZAR+Bu4EzgAUwybwMs8X0vIhJdP/4Il1wCZ58N\n06ZBrcStHFeXzLcDa3y3fwW+BJoDlwLTfeenA7G17JiION/27eByQbduMHGiKbEksED+9alAOvAR\n0ATY4Tu/w/e9iEh0/Oc/0LkzXH01PPkk1Ijkat7xwd+/SeoBs4HBwJ5K91m+4wjZ2dnlt10uFy6X\nK+AARUQOs2kTdO8OQ4bA3/5mdzQh83g8eDyekF/Hn19ntYF5wHzged+5rwAXpgzTDNNJekal52lz\nChEJr7VrTWfnI4/AbbfZHU1EBLs5RXVllhrAFGA9hxI5wFzgJt/tm4A5gb6xiEhAPv7YrEH+7LOO\nTeShqC77dwIKgM85VEp5EPgXMAtoBXwNXAX8XOm5apmLSHgUFMCVV8Irr5jx5A4WbMtce4CKSGxb\nuNDsDjRzpqmVO1ykyiwiIvZ5910YONDM6kyARB4KJXMRiU0zZpjVDxcsgIsusjuamJe406VExHZu\ndwE5OXkUF9ciOfkgWVk96dOnC0yeDI89ZrZ5O+ssu8OMC6qZi4gt3O4CBg9eiNf7ePm5tLQR5Hbd\nxZkfzDdL2Kal2RihPYKtmatlLiK2yMnJOyyRg8VAby3qbZ8FX62BFi1siy0eqWYuIrYoLq7YlrR4\njnu5jDlknnuzEnkQlMxFxBbJyQcBqEkJk7mTC1nBxeRT9LsUmyOLTyqziIgtsrJ6sn3TfYzespG6\n/EYPFtEkbQyZmb3sDi0uqQNUROzxzTfs7tqNDw+cwHOn9aVWHcjM7GFGsyQwzQAVkfjx0UcwYADc\nd59Z/VBL2JbTaBYRiQ9vvAGZmfDqq9Cvn93ROIaSuYhEh2WZpWunTTOTgc491+6IHEXJXEQib98+\nuOUW2LIFVq2Cpk3tjshxNDRRRCJr+3a4+GJTF8/PVyKPECVzEYmczz6DDh3gT3+C116D44+3OyLH\nUplFRCIjN9eUViZMMBsvS0QpmYtIeFkWjB1rjnnzTMtcIk7JXETCZ/9+swb56tWwciW0amV3RAlD\nyVxEwuPHH+GKK6B+fVi+HOrVszuihKIOUBEJ3YYN0LEjtG8P77yjRG4DJXMRCc3ixdClCzz4IDz9\nNCQl2R1RQlIyF5HgTZ4MN9wAs2aZkStiG9XMRSRwJSVw771ms+Xly6F1a7sjSnhK5iISmN274Zpr\n4MABM2LlxBPtjkhQmUVEArFlC1x4IaSmwvvvK5HHECVzEfHPihUmkQ8aBBMnQu3adkckFajMIiLV\nmzEDhg6F6dOhd2+7o5EqKJmLyNGVlsJDD8Hrr5sVD9u2tTsiOQp/yiyvAjuAtRXOZQPfAYW+Qzuw\nijjN3r1mgSyPx2zzpkQe0/xJ5lM5MllbwFgg3XcsCHNcImKnrVuha1ezZO2SJdC4sd0RSTX8SebL\ngF1VnNcOrCJO9OmnZqXDAQNMjTw52e6IxA+hjGbJBD4DpgANwhOOiNjq3XchIwOefx6GDze7A0lc\nCLYDdBLwqO/2aOA54NbKD8rOzi6/7XK5cLlcQb6diESUZcFTT5khhwsWQLt2dkeUMDweDx6PJ+TX\n8ffXbiqQC5wTwH2WZVlBByYiUVJcDHfcAV98AXPnQvPmdkeU0GqYv4YC/pMo2DJLswq3B3D4SBcR\niQNudwF/vngon5/chuULPmLB8CeUyOOYP8l8JrAC+D3wLXAL8BTwOaZm3hUYEqkARST83O4Cpt/x\nCs943mHe7uvpsnM999xfgNtdYHdoEqRI9m6ozCISi0pKmHbGJfTa/BW38zLz6Fd+V0bGKBYsGG1j\ncBJsmUUzQEUSydatMHAg53//Ne34hK0cXlYpKtLGEvFKC22JJIr5880ola5debD9dUckcoCUlBIb\nApNwUMtcxOn27zdbur31Frz5JnTpwj3uAjb9ewRe7+PlD0tLG05mplbmiFdK5iJOtnkzXHstnHIK\nFBbCSScB0KdPFwDGjx9FUVESKSklZGb2Kj8v8UcdoCJONXMmZGXBww/D3XdrNmecUAeoiBi//QaZ\nmfDhh7BoEZx/vt0RSRSoA1TEST77DP7wB7MO+SefKJEnECVzESewLLOuSvfuMGIETJsG9erZHZVE\nkcosIvHup5/g1lvhm2/MPp2nn253RGIDtcxF4tny5ZCeDqeeqkSe4NQyF4lHJSXwxBOmtDJlCvTp\nY3dEYjMlc5F4s3UrXH+9uf3JJ1rpUACVWUTii9sNF1wAl1wCixcrkUs5tcxF4kFxsZmSP3u2mZbf\nubPdEUmMUTIXiXWbN8M110DLlmZKfsOGdkckMUhlFpFY9tpr8L//CzffDO+8o0QuR6WWuUgs+vVX\nuOceWLXK1MbPO8/uiCTGqWUuEmvWrDHrjtesaUarKJGLH5TMRWKFZcGECdCjh1np8NVXoW5du6OS\nOKEyi0iUud0F5OTkUVxci+Tkg2Rl9aRPx7ZmSv5338HKldC6td1hSpxRMheJIre7gMGDFx62w8+J\nXwyk24E8jr9xIMyaBccdZ2OEEq9UZhGJopycvPJEXpMSRvEo47Yu5omW3eDZZ5XIJWhqmYtEUXGx\n+ZFLZQuvcgul1KQdn9DmhJdsjkzinVrmIlFUp3Yxw3iKj/kj8+lNT/LYximkpJTYHZrEObXMRaJl\n5Ure2DyLNcdD+33/YgunAZCWNpzMzF42ByfxTslcJNJ+/tmsq/Lee/xu3Dj21G1KmwlTaVWUREpK\nCZmZvejTp4vdUUqci+R23ZZlWRF8eZEYZ1lmdMqQIdC/P4wZAw0a2B2VxLgaNWpAELlZLXORSPj3\nv+Guu8za47Nnm/VVRCLInw7QV4EdwNoK5xoCi4CNQB6g5oYIwIED8OST0L49dOtmpuMrkUsU+JPM\npwKVe2cewCTzNsAS3/ciiW3FCrNxxNKl8PHHMGwY1K5td1SSIPyty6QCucA5vu+/ArpiWuxNAQ9w\nRqXnqGYuiWHXLtPBmZsL48bBn/8MNSLZHSVOFmzNPNhx5k0wiRzf1yZBvo5I/LIsmDkT2rY1Kxyu\nWwdXXaVELrYIRweo5TuOkJ2dXX7b5XLhcrnC8HYiMcDrNR2c27apg1NC4vF48Hg8Ib9OKGUWF7Ad\naAbkozKLJIL9++G558wxbJgZdqi6uIRRtIcmzgVuAp7yfZ0T5OuIxI8PP4RBg6BVK9PBeeqpdkck\nUs6f7D8T09nZCFMffwh4D5gFtAK+Bq4Cfq70PLXMJeZUuZZ4dbMvd+2C++8Htxuefx6uvFJ1cYmY\nSLbMrz3K+e6BvpmInapaS9zrHQFQdUK3LHjjDRg6FC6/HNavh/r1oxWuSEA0nV8SRkbGSPLyHqvi\n/CgWLBh9+MmyDs7t2+Gll6BDhyhFKYku2kMTReJO2VrilRUVJR36Zv9+eOIJk7x79IDVq5XIJS5o\nbRZJGMnJB6s8X76W+PLlpoMzNdUk8dTUqMUmEiq1zCVhZGX1JC1txGHn0tKGM/QvHeCOO+Dqq+GR\nR2DePCVyiTtqmUvCKOvkHD9+FEVFSaQkH2TMuUmkD7kdrrhCHZwS19QBKolp82b4619h5051cEpM\nUQeoiD9274bsbOjYETIy1MEpjqFkLonht9/g6aehdWuzccTq1XDffZqKL46hmrk4W3GxKaOMGQMX\nXQQeD5x1lt1RiYSdkrk404EDMG0ajB4N551npuKnp9sdlUjEKJmLs5SUwOuvm7r4aaeZDZU7drQ7\nKpGIUzIXZygtNeuKP/wwNGwIU6aA1s+XBKJkLvHNskwJZdQoSEqCsWPNKBWtaigJRslc4pNlwZIl\nMHKkGakyejT0768kLglLyVziz/LlJolv22am3191ldmDUySBKZlL/Fi92pRTvvzS1MYHDoRauoRF\nQJOGJB6sXQsDBsBll8Gll8LGjXDzzUrkIhUomUvs2rABrr3WrCvepQts2mTWUznuOLsjE4k5SuYS\ne77+Gm65BTp1gnPPNYtiDRkCxx9vd2QiMUvJXGLHf/9rtmpr1w5atDAt8QcfhHr17I5MJOYpmYv9\ndu40myafc45J3Bs2wKOPQoMGdkcmEjeUzMU+u3bBiBFw5plw8CCsW2dWNmzUyO7IROKOkrlE3+7d\nZpLP6aebVnlhIeTkQLNmdkcmEreUzCV6fvkFnnnGJPGNG2HVKnj5ZWjVyu7IROKeBupK5H36KUya\nBG+/Db16wQcfQNu2dkcl4ihK5hIZe/fCG2/A5MmmlDJoEHz1FTRpYndkIo6kDZ0lvNavhxdfhBkz\nzM4+d95pVjFMSrI7MpG4EOyGzmqZS+j274d33jGt8A0b4LbbTKemauEiURNqy/xrYDdQAhwA2le4\nTy1zp9uyxeyvOXUqnH22aYX3769NkkVCYFfL3AJcwE8hvo7Ei5ISsxnE5Mnw8cdw442wdCn8/vd2\nRyaS0MJRZtFuADHO7S4gJyeP4uJaJCcfJCurJ336dAnsRbZuNVuxvfwyNG9uFryaPdvv9VLCEoOI\nHFU4WuaLMWWWF4GXQ45IwsrtLmDw4IV4vY+Xn/N6RwBUn0xLS80wwsmTzderr4a5c+H886MXg4hE\nRdmUvcbAGqBzhfsssV/PniMss8fa4UdGxsijP+mHHyzr2Wct6/TTLevccy1r0iTL2r07ujGIJChM\nIzlgobbMt/m+fg+8i+kAXVZ2Z3Z2dvkDXS4XLu2WHnXFxVX/FxcVVRoqaFmwcqVphc+dazaBmD4d\nOnYMeV9Nv2MQSUAejwePxxPy64SSzOsAScAeoC7QE3ik4gMqJnOxR3LywSrPp6SUmBt79pgx4ZMn\nw759ZkTKuHFw0knRi0EkgVVu6D7yyCNHf/AxhLI2SxNMK3wN8BEwD8gL4fUkArKyepKWNuKwc2lp\nwxnRJ9Uk7latzC73Y8eaGZpDh4Y1kR8rhszMHmF9H5FEphmgCcDtLmD8+EWU/lZKz5/XcusBLyfu\n3QN33GF29DnllKjFUFSUREpKCZmZPdT5KVKFYMeZK5k7XVEReDyQmwtvvgkdOpgWee/e2hBZJAZp\nOr8csmOHmdiTm2uGFJ5zDvTrZyb5nHqq3dGJSASoZe4ElgWff26Sd26uWR+lZ0+TwHv31s49InFE\nZZZEU1QE+fkmec+bZ9ZD6dfPHJ07w3HH2R2hiARBZZZEsH37ofJJfj6cdx707QsLF8IZZ4Q8HlxE\n4pda5rHMsmDNGtPyzs2FTZvM2uB9+5rySZiHEIqI/VRmcYp9+0yn5bx55khJMaWTvn1N+UTLy4o4\nmsos8WzbtkPJOz8f0tNNAl+0yCwtq/KJiFRDLXM7WJbZiaes89LrNeWTfv3MhscNG9odoYjYRGWW\nWLdvn5k2X5bA69Y9NPrkootUPhERQGWW2GJZZjOHNWtMC3zVKigogAsuMMk7Px/atLE7ShFxELXM\nQ1VaCps3m6Rd8bAsU/tOT4d27aBHDzjxRLujFZEYpzJLNOzfD+vWHZ60P//cDBEsS9zp6WYnnubN\ncb+/TFuliUhAVGYJtz17TJmkrFRSWGimyZ922qGkffnlJnFX0eLWVmkiEk2ObZkHtIHwjh1Hlkm2\nboWzzz68xX3OOX5vYJyRMZK8vMeqOD+KBQtGh/JPExEHi8mWeUbGSFtKC0dtFVsWfc5qeWTiLio6\nlLD794fsbNNBGcISsdoqTUSiKaLJPC/vseiXFiyLV8bmYnnvpAOraMNG0ikk3VvIHwY8B00bH0rc\nt99uvrZqFfaJOdoqTUSiKaJllrJNpkMqLVgW7N4N338PO3earxWPKs7tO1DKttIWfE9jvKRRaNI5\ndS/M470Pnw7jP/HoqvrrIC1tOC+80Es1cxE5qpgss5Q5rLRgWfDzz9Um5PJzP/wAycnQuDGcfLL5\nWnY0b246ICudv+yyx6uuV5+wNBr/XODQXyLjx4+qsFWaErmIREZEW+YvkEljvueMhqtIb/G7Q8m5\nTp3Dk3LlJF3xXKNGZrGpAKhVLCLxKiZb5l7S2HLyfzgt8x7oe4lJ0I0aRXzjBLWKRSTRRLRlnpEx\nUruwi4gEQDNARUQcINhkXjP8oYiISLQpmYuIOICSuYiIAyiZi4g4gJK5iIgDhJLMewFfAZuA+8MT\njoiIBCPYZJ4ETMAk9LOAa4EzwxWUHMnj8dgdgqPo8wwvfZ72CzaZtwc2A18DB4A3gP5hikmqoB+W\n8NLnGV76PO0XbDJvDnxb4fvvfOdERMQGwSZzTe0UEYkhwU7n7whkY2rmAA8CpcBTFR6zGUgLOjIR\nkcTkBVpH681q+d4wFTgOWIM6QEVE4lJvYAOmBf6gzbGIiIiIiIg/k4dyfPd/BqRHKa54Vd3n6QJ+\nAQp9x8ioRRZ/XgV2AGuP8Rhdm/6r7vN0oWvTXy2BfGAd8AWQdZTHRe36TMKUWVKB2lRdO/8T8L7v\ndgdgVSQDinP+fJ4uYG5Uo4pfnTE/AEdLPro2A1Pd5+lC16a/mgLn+27Xw5SsQ8qdoa7N4s/koUuB\n6b7bHwENgCYhvq9T+TsZK5KbijjJMmDXMe7XtRmY6j5P0LXpr+2YxhrAr8CXwCmVHhPQ9RlqMvdn\n8lBVj2kR4vs6lT+fpwVciPmz633McgoSHF2b4aVrMzipmL94Pqp0PqDrM9QNnf2dPFT5t7UmHVXN\nn8/lU0y9bS9mRNEcoE0kg3I4XZvho2szcPWAt4HBmBZ6ZX5fn6G2zP+L+c8r0xLz2+NYj2nhOydH\n8ufz3IP5YQGYj6mtN4x8aI6kazO8dG0GpjYwG5iB+cVXWVSvT38mD1Us4ndEnUzH4s/n2YRDv63b\nY+rrcnSp+NcBqmvTP6kc/fPUtem/GsA/gHHHeEzUr8+qJg8N8h1lJvju/wy4INIBxbnqPs+7MUOZ\n1gArMP/JUrWZwFZgP6b2eAu6NkNR3eepa9N/nTBLoKzh0FDO3uj6FBERERERERERERERERERERER\nEREREREREZEy/w8RXyOeesk2DAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10b162850>"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# It is important to understand what is going on,\n",
      "# but we don't have to work so hard to fit a polynomial\n",
      "\n",
      "b = np.random.randint(0, 10, 6)\n",
      "x = np.linspace(0, 1, 25)\n",
      "y = np.poly1d(b)(x) \n",
      "y += np.random.normal(0, 5, y.shape)\n",
      "\n",
      "p = np.poly1d(np.polyfit(x, y, len(b)-1))\n",
      "plt.plot(x, y, 'bo')\n",
      "plt.plot(x, p(x), 'r-')\n",
      "list(zip(b, p.coeffs))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 26,
       "text": [
        "[(4, -104.9127846499792),\n",
        " (8, 338.10593687644877),\n",
        " (1, -273.52108982003318),\n",
        " (6, 26.406479516940809),\n",
        " (6, 45.218269345183614),\n",
        " (9, 3.4534675105311239)]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHxJJREFUeJzt3Xl4VPW9x/E3BEiKWhBQBKs3ml5cqIg7CuK4QNDcuuN6\nq0W95T7lJtSllyUIQQSqXlASLVUoFdCKFoUCo5CIjEEtYEFkc4Ep3IIILhHFKwkkmfvHb0IWssx2\nljnzeT3PeWZyMnPO1yN8+Z3v+S0gIiIiIiIiIiIiIiIiIiIiIiIiImKbDGA1sB7YAkwO7y8AdgEf\nhLdBTgQnIiKRax9+bQOsAvoB44AHHItIRETqaR3BZ34Iv7YD0oBvwj+3siQiERGJWiTJvDWmzLIX\nWAFsDu/PBT4E/gh0tCQ6ERFJuA6YMosPOB7TMm8FPIpJ6CIi4pA2UXz2W8APnA8E6uyfCSxu+OGs\nrKxQMBiMKzgRkRQUBH4a7ZdaKrN0obaE8iNgAKb3ygl1PnMDsPGIaIJBQqGQtlCIcePGOR6DWzZd\nC10LXYtmtkmTALKiTeTQcsu8GzAbk/RbA3OB5cAcoDcQArYDQ2M5uYiIhH39NUyZEvPXW0rmG4Fz\nG9l/V8xnFBGRI02eDIMHwx/+ENPXo6mZS4x8Pp/TIbiGrkUtXYtaXroWfn8phYXFVFS0IT29kry8\ngeTk9G/+S//8J/zpT7BpU8zJ3Mq+4qFQKGTh4UVE3MXvL2X48GUEgxMP78vKymfatOzmE/qQIXDi\nifDoo7Rq1QpiyM1K5iIiCZKdPYbi4kcb2f8wS5dOaPxLmzbBFVfA1q3QoUPMyTySQUMiIhKBiorG\nK9fl5WlNf2n0aBg1Cjp0iOvcqpmLiCRIenplo/szMqoa/8I778CGDfCXv8R9brXMRUQSJC9vIFlZ\n+fX2ZWWNJjd3wJEfDoVgxAh45BFIT4/73GqZi4gkSM1DzqKihykvTyMjo4rc3EGNP/xctAj274c7\n70zIufUAVETEbpWV0KsXPPEE5OQAtV0ai4sngh6AiogkgTlz4Ljj4JprgNoujY31hImUWuYiInY6\ncAB69DAPPfv0ARp2aVTXRBER9ysqggsuOJzIoekujdHQA1AREbt8842pk69cWW93U10ao6GWuYiI\nXX73O7jhBjj99Hq7G+vSGC21zEVE7LBrF8ycaQYJNVC3S+OyZbEdXg9ARUTscN99pgfL5MnNfizW\nuVnUMhcRsdqWLfDXv5rJtCyimrmIiNXy883Q/Y4dW/5sjFRmERGx0nvvwW23waefQkZGix/XFLgi\nIm4TCsHIkTB+fESJPB5K5iIiVvH7oawM7rJ+2WQlcxERK1RVmVb5pEmQ1sziFAmi3iwi4mkxLbCc\nCC+8YB54/vzn1p8LJXMR8bDGFlgOBs1IS0sTenk5jB0Lf/4ztLKyn0mtlsosGcBqYD2wBajp7d4J\nKAE+BYoB6/rbiIjEqLCwuF4iBwgGJ1JUVGLtiSdNgvPPh759rT1PHS0l83LgcqA30Cv8vh8wEpPM\newDLwz+LiLhKTAssx2vzZpg+HQoLrTtHIyJ5APpD+LUdkAZ8A1wLzA7vnw1cn/jQRETiE/UCy/Gq\nroZf/cqs63niidacowmRJPPWmDLLXmAFsBnoGv6Z8GtXS6ITEYlDVAssJ8Kzz5rXoUOtOX4zInkA\nWo0ps3QAlmFKLXWFwtsRCgoKDr/3+Xz4fL5YYhQRiUlUCyzH67PPzEPPQABaR97rOxAIEAgE4j59\ntI9ZHwYOAPcBPmAP0A3TYj+9wWc1nF9EUseNN8LPfmZKLHGwajh/F2p7qvwIGAB8ACwC7g7vvxtY\nGO2JRUQ8Y8ECMzPi6NGOhdBS9j8L84CzdXibCzyB6Zr4CnAysAO4BdjX4LtqmYuI9337LfTsCS++\nCJddFvfhYm2Za9ZEEZF4DBsGBw/CjBkJOZwWpxARsdt775kSy+bNTkeiibZERGJy8CD8x3/AU0/B\nscc6HY2SuYhITB57DE49FQYPdjoSQDVzEZHoffKJmXdl3To4+eSEHlorDYmI2KFmyP7YsQlP5PFQ\nMhcRicasWWaK22HDnI6kHpVZREQitWcP9OoFb75pXi2gfuYiIla79Vbz0HPy5JY/GyP1MxcRsdKS\nJbB2LTz/vNORNEotcxGRluzfbybRmjULrrzS0lOpzCIiYpXf/Ab27bOlVa4yi4iIFdasgXnzXDFk\nvznqmigi0pRDh8yQ/SlToHNnp6NplsosIiJNeeAB2LoVFi2CVlamy1oqs4iIJNKrr5oZEdeutS2R\nx0MtcxGRhrZuNXOvvP46nH++rafW3CwiIonwww9w001mLU+bE3k81DIXEakRCsGQIVBZCXPnOlJe\nUc1cRCRes2bB++/D6tX4X19JYWExFRVtSE+vJC9vIDk5/Z2OsElK5iIiAOvXw8iRUFqK/+11DB++\njGBw4uFfB4P5AK5N6KqZi4js2wc33wxFRXDGGRQWFtdL5ADB4ESKikocCrBlSuYiktpq6uRXXw23\n3QZARUXjRYvy8jQ7I4uKyiwiktqmTIHdu82Q/bD09MpGP5qRUWVXVFFrqWV+ErAC2AxsAvLC+wuA\nXcAH4W2QRfGJiFhn5Ur4n/+BV16B9PTDu/PyBpKVlV/vo1lZo8nNHWB3hBFrqfvLCeFtPXA0sBa4\nHrgF2A9Mbea76pooIu61dy+cdx7MnAmDjmyP+v2lFBWVUF6eRkZGFbm5A2x5+GnXFLgLgaeBvsD3\nwJRmPqtkLiLuVFUFAwZAv35mcJCL2DECNBM4B1gV/jkX+BD4I9Ax2hOLiDhm7Fho3RrGjXM6koSJ\n9AHo0cB8YDimRT4dqPnnbAKmhX5vwy8VFBQcfu/z+fD5fLFHKiKSCH4/zJljJtBKc753SiAQIBAI\nxH2cSJrybYElwBvAU438PhNYDJzVYL/KLCLiLjt2wEUXwWuvmYm0XMiqMksrTBllC/UTebc6728A\nNkZ7YhERW1VUwODBZpSnSxN5PFrK/v2AUmADUNPMHg3cDvQO79sODAX2NviuWuYi4h7DhsGePTB/\nvqvnJ9eCziIiTZk6FZ57Dlavhg4dnI6mWZo1UUSkMc89Z+ZcKS11fSKPh5K5iHjXiy+afuSBAJx0\nktPRWEplFhHxpoUL4T//E5Yvh549nY4mYiqziIjUKCmBX/0K3ngjqRJ5PJTMRcRb3nkH7rgDFiww\nc6+kCM1nLiLesXYt3Hgj/PnPZt6VFKJkLiLesHkz5OTAjBlmEq0Uo2QuIslv2zYYOND0J7/uOqej\ncYSSuYgkt5074aqrzAyId9zhdDSOUTIXkeS1d69J5Hl5pvdKClMyF5HkVFZmauN33AEPPOB0NI7T\noCERST7795sW+aWXwhNPuHrirGhpoi0RSQ0//ADXXAOnnw7Tp3sqkYOSuYikgh9+gJtvhk6dzGpB\nrb1XKbZjDVAREefs2gX9+0OXLvD8855M5PHQ1RAR91u1yiz3dsstMHs2tNFMJA3pioiIu82ZAw89\nBLNmwb/9m9PRuJaSuYi4U1UVjBplFl9esSJlZj+MlZK5iLjPd9/B7bfDgQNmqbfOnZ2OyPVUMxcR\nd9m2Dfr0gcxMWLZMiTxCapmLiO38/lIKC4upqGhDenoleXkDycnpD2+9ZVrk48ebVYIkYupnLiK2\n8vtLGT58GcHgxMP7sk4dzcKBB/jZgpfgpZfg8ssdjNBZWjZORJJCYWFxvUTehkM89I8yjn7hNfhw\nFZx6qoPRJa+WauYnASuAzcAmIC+8vxNQAnwKFAMdrQpQRKzn95eSnT0Gn6+A7Owx+P2llp2roqK2\nDdmZryhhAN3Zza/PvleJPA4ttcwPAfcD64GjgbWYJD4k/Po4MAIYGd5EJMk0VvYIBvMBTB07wdLT\nKwHoySYWcS0vcytjeJQBRxck/FyppKWW+R5MIgf4HvgIOBG4Fpgd3j8buN6S6ERcws6Wq90alj0A\ngsGJFBWVWHK+4cOu4NHOAwjgYyyPMJrJnJL1MLm5qbfUWyJFUzPPBM4BVgNdgb3h/XvDP4t4kt0t\nV7vVLXvUVV6elviTvf8+1xQ8xIUntmLEGbewK20b2RkPk5s7yBPX0kmRJvOjgVeB4cD+Br8Lhbcj\nFBQUHH7v8/nw+XxRByjitKZbrg97IgHVlD0aysioStxJvv0W8vNh/nx44gm6/Pu/80ePTV0bq0Ag\nQCAQiPs4kSTztphEPhdYGN63FzgBU4bpBnzR2BfrJnORZGVry9UBeXkDCQbz63cVzBpNbu6g+A8e\nCpkEfv/9cPXVsGWLmb5WDmvY0B0/fnxMx2kpmbcC/ghsAZ6qs38RcDfwWPh14ZFfFfEGW1quDqq5\nuygqepjy8jQyMqoSU/bYvh2GDYN//hPmzYN+/RIQrTSlpfucfkApsIHaUsooYA3wCnAysAO4BdjX\n4LsaNCSe0Oggl6zRTJumOm+jDh2CqVPNcm4PPmi2du2cjippaKUhEQv5/aUUFZXUabkOUCJvzHvv\nwdCh8JOfwDPPqN94DJTMRcQ533wDI0fCkiXw5JMweLDn1ua0i5aNExHb+Ret4PGzruPrE05m0Rtr\nWfbkDLMakBK57ZTMRSR633/P5qF5nHXTjfTd9BVXH1zOdTv/zrDR73pqQFUyUZlFEqLJKU3FWz77\nDIqKYOZMStscx4i9s1jFxfU+kp39MEuXTnAowOSnWRPFMV4fISnA+vWmh8qSJfCLX8CaNYy9Zw6r\n9l58xEe90v8+2ajM4lF2ziVi99weYpPqanj9dbjySrOQcs+eEAzCtGlw6qme73+fbNQy9yC7W8pe\nHyGZcsrL4YUXTEu8XTvTT/zWW4/oK27pyFGJmpK5B9k9l0g8LTTV2l3kyy/h97+H6dPhvPPg6afN\nij9N9EyxbOSoxETJ3IPsbinH2kJTrd0FduyAxYth0SJYs8Z0K3zrLTjzzIi+npPTX/+vXELJ3IPs\nrmXG2kLz+myE4MI7j+pqWLcO/vpXk8A//9zUw4cNg4UL4aijnItN4qJk7kFO1DJjaaF5vdbumjuP\n8nLT2l60yLTCf/xjuPZaU1Lp0wfSvHG9U52SuQclSy3T670hHL3z+PJL8PtNAl++HM4+2yTwFSug\nRw9rzy2OUDL3qGSoZXq9N4Rtdx6HDpl5wtetM9v778NHH8FVV8H118Nzz0GXLok9p7iOkrk4Jlnu\nIGJlyZ3HgQOwYYNJ2h98YF63bIHMTDj3XLZkdOD5qp+y/uyraPU95HXOJEeJPCUomYujkuEOIlZx\n3XmEQrBvH2zcWNviXrcO/vEPOO00OPdcsw0ZAr16wVFH1anRP3P4MOodlDo0N4uIhZqcB728HHbt\nMqvw7Nxptpr3Na/V1XDWWXDOObXJu2dPSE9v9FzZ2WMoLn60kf2aKyWZaG4WSSmu6/IXCsH+/VBW\nZub2Dm85ZWXkXHGMmaBq504Y9wbcs9O0urt3h5NPhpNOMq/nnGMeUtbs69gxqqlkvd47SJqnZC5J\nJ+Fd/qqqTCKOZNu374iETVmZ2f+jH8Gxx5qtU6fa98ceC//yL3DppbWJu2tXaN301Eix/GPl9d5B\n0jwlc0k69bv8hWjPDxwI/ppFkyeQ0zHNJNh9++pvDfd9+21tgq6oMINljjmm5e2MMxpP2B07Jmyd\ny1j/sfJ67yBpnmrm4k7V1fDVV6Y8sXt3ve29+QHSy46mO7vpwlccoi3fcCwH2x/ilN7/ahJr3a0m\n2dZsHTqYrSZBt2/fbCvZbvHUvrVWafJTzVySS3W1edC3ZYvZduyon7T37DEJt3v3+lvv3hSvKcNf\n9t/spjtf0YWDmAeC2Zd640FfPLVvL/cOkuYpmYu1qqpg+/bapF2zffyxaSWfeabZevQwM/TVJO0T\nTmiy18YFJ/fkheEL2e3RcoJq3xILJXNJnB07agex1GyffgrHH1+btH0++PWvTe25Q4eYTuP1wUaq\nfUssIqnLzAJygC+As8L7CoD7gC/DP48Cljb4nmrmXhYKwbZt8PbbZistNQ8SL7zQ9IWuSd6nnQZH\nH+10tElHte/UFWvNPJIvXAp8D8yhNpmPA/YDU5v5npK5l4RCpqVdWlqbvNPS4LLLoH9/89qjR1T9\nokXkSFY+AF0JZDZ2zmhPJkmkutrMAVI3eR9zjEna2dkwaRKccoqSt4hLxFMzzwXuAv4OPAjsS0hE\n4pzycnjjDXj5ZSguhuOOM8n7xhvhqafMgBcRcaVYk/l04JHw+wnAFODehh8qKCg4/N7n8+Hz+WI8\nnVjm4EEoKTEJfPFi6N3bLN775JPQrZvT0Yl4XiAQIBAIxH2cSO+RM4HF1NbMI/mdauZuVVlpVp55\n+WWzVNgZZ5gEfvPNSuAiDrN70FA34PPw+xuAjTEeR+xSVQUrV5oE/uqrpt59661QUKDyiYgHRJLM\nXwIuA7oAOzE9WXxAbyAEbAeGWhSfxKO6Gv72N5PA5883kzvddhusXm2SuYh4huZm8aIDB2DuXFP3\nbt3aJPBbb9XajyJJQHOziJnP5Pe/h2efhYsugunTTW8UdR8UD3DdHPYuo2TuBZs2wdSpsGAB3H67\n6RN+2mlORyWSMAmfw96DVGaxScJbFaGQ6Qs+dapZJ/K//guGDoXOnRMXtIhLpNKSeCqzuFhCWxXl\n5fDiiyaJt2kDDzxgauJNzDAo4gVaEq9l7pmR38Pqr4xjBIMTKSoqifwgX3wB48dDZia89hoUFsL6\n9XD33Urk4nmaFrhlSuY2iKtV8fXXcP/9pga+ezesWAF+P1x5pR5sSsrIyxtIVlZ+vX1mWuABDkXk\nPiqz2CCmVkV5OTz9NDz2GAwebBZz6NrVoghF3M3rc9gngh6ARimWB5mN1cyzskYzbVojfxirq80g\nn9GjoVcvk8xPP92K/xQRcSE9ALVBrA8yI25VlJbCQw+ZnirPP2/6iIuIREAt8yhY1j3qk09gxAjz\nQHPSJNM7xUWrxYuIfWJtmStjRCHh3aO++AKGDYN+/aBvX1MXv+MOJXIRiZqyRhQS1j3qwAGYPNms\nkdm2rUniv/0tZGQkIEoRSUVK5lGIu3tUdTXMmWO6Ga5bB6tWmRV8NGrTNn5/KdnZY/D5CsjOHoPf\nX+p0SCIJoQegUYire9TWrXDPPXDoEMybB5dcYnG00pDm9xAv0wNQq1VVmdb35MkwdqypkadpCLIT\nUml+D0le6proRlu2wJAhcNRRsGYNnHqq0xGlNM3vIV6mmrkVDh2CiRNNP/F77oE331QidwHN7yFe\npmSeaOvXw4UXmvU2164109Kqq6EraH4P8TLVzBOlosK0xv/wB3j8cTOboSbCch2/v5SiopI6D7AH\n6OGnuEqsNXPXJPOkXhJqzRpTTsnKMku1de/udEQikqSS+gFo0nYZO3AAxo2D2bNNj5XbblNrXEQc\n4YpibkIWb4hS3INH3n0XeveGHTtgwwaz9qYSuYg4xBUtc7u7jMV1J1BZaVrjs2aZ+cZvusmSGEVE\nohFJy3wWsBfYWGdfJ6AE+BQoBjrGE4TdXcZivhP47DO44gpTI1+/XolcRFwjkmT+J2BQg30jMcm8\nB7A8/HPM7O4yFtOdwLJlcP75MHAgLF2qVX9ExFUiKbOsBDIb7LsWqFk5YTYQII6EbveSUFHdCdSU\nVZ5/Hl56CXw+S2ISEYlHrDXzrpjSC+HXuJupOTn9beu5kpc3kGAw/4hl3HJzG9yA7N5tHmy2a2dm\nOXSgNZ7UXTZFxDaJeAAaCm9HKCgoOPze5/Phc0mrNqI7geJiM/Bn2DAYNcqRybGStsumiEQsEAgQ\nCATiPk6kfekygcXAWeGfPwZ8wB6gG7ACaLjqcHKOAK2shIICU1Z54QVHyyqa5U8k9dg9aGgRcDfw\nWPh1YYzHcZe6ZZW1ax1/yOnELH8q64gkp0iS+UuYh51dgJ3AWOB3wCvAvcAO4BaL4rOPC8oqDdnd\nZVNlHRFpTCgpHDoUCuXnh0InnhgKrVjhdDT1LFnydigra3QIQoe3rKxRoSVL3rbkfAMH5tc7V82W\nnT3GkvOJyJFo4hlkS1wxAtQxe/aY+VTatnVFWaUhu7tsavEGkeSVusl83Tq4/nqzEtDYsa4oqzTG\nzi6bWrxBJHm5YqIt273yCmRnw5NPwvjxrk3kdtPiDSLJyzXzmduiutp0O5wzBxYuNLMeSj1avEHE\nWUm/OIXlvv8e7roLvvgCXnsNjj/e6YhERI4QazJPjTLLjh3Qty8ceywsX65ELiKe4/1kvnIlXHyx\nWdZt5kxIT3c6IhGRhPN2b5YZM2DMGJg710xdKyLiUd5M5pWVcP/9UFJiWuY9ejgdkYiIpbyXzMvK\n4JZbzECgVaugY1yLIImIJAVv1cw/+gguugjOPhuWLFEiF5GU4Z1k/vrrcNllkJ8PU6ZoIJCIpBRv\nlFmmTYPHHjMDgS65xOloRERsl9zJvKoKHnjAPOh87z3IzHQ6IhERRyRvMv+//4M774TvvjOJXPVx\nEUlhyVkz37sXLr8cOnSApUuVyEUk5SVfMv/oI+jTB3JyzDqd7do5HZGIiOOSq8yyYoVZTOLxx80S\nbyIiAiRTMp87Fx58EObNgyuucDoaERFXcX8yD4VgwgSYNcu0zHv2dDoiERHXcXcyP3gQhg6FjRvN\n0PwTTnA6IhERV3JvMt+3D26+Gdq3h7ffhqOOcjoiERHXcmdvlv/9X+jXD848ExYsUCIXEWlBvC3z\nHcB3QBVwCLgw3oBYuxauvRZ++1v4zW/iPpyISCqIN5mHAB9QFn8omJkOhwyB556DG25IyCFFRFJB\nImrmiVkUevp0eOQRk9AvuighhxQRSRWJaJm/iSmzPAvMiPoI1dUwapSZ8fCddyArK86QRERST7zJ\nvC/wOXAcUAJ8DKys+WVBQcHhD/p8Pnw+X/1vl5fDL38Ju3aZybI6d44zHBGR5BIIBAgEAnEfJzEl\nEmMc8D0wJfxzKBQKNf3psjK47jro1g3mzIGMjASGIiKSnFq1agUx5OZ4uia2B44Jvz8KGAhsjOib\n27ebRST69DHD85XIRUTiEk+ZpSuwoM5xXgSKW/zW+++bFnl+PgwbFsfpRUSkRiLLLA0dWWZZtAju\nuw9mzjR9yUVEpJ5Yyyz2Ded/5hmYOBH8frjgAttOKyKSCqxP5tXVMGIELF4M774Lp5xi+SlFRFKN\ntcm8vBzuugv27DFdDzt1svR0IiKpytqJtq66CtLSoLhYiVxExELWJvN+/eDFF9X1UETEYvb2ZhER\nkWY5MWhIRERcQslcRMQDlMxFRDxAyVxExAOUzEVEPEDJXETEA5TMRUQ8wNJknp09Br+/1MpTiIgI\nFs/NUlz8KMFgPgA5Of2tPJWISEqzvMwSDE6kqKjE6tOIiKQ0W2rm5eVpdpxGRCRl2bI4RUZGlWXH\n9vtLKSwspqKiDenpleTlDVRJR0RSjuXJPCtrNLm5gyw5tt9fyvDhywgGJx7epxq9iKQii3uzPMy0\naYMsS6yFhcX1EjmoRi8iqcnSlvnSpROsPDwVFY2Hrxq9iKSapB40lJ5e2eh+K2v0IiJuFE8yHwR8\nDGwFRiQmnOjk5Q0kKyu/3j5Tox/gRDgiIo6JdaWhNOAT4CrgM+B94HbgozqfsWWlIb+/lKKiEsrL\n08jIqCI3d4DrHn4GAgF8Pp/TYbiCrkUtXYtauha17F5p6EJgG7ADOATMA66L8Vhxycnpz9KlEwgE\nCli6dILrEjmYP6hi6FrU0rWopWsRv1iT+YnAzjo/7wrvExERB8SazLVSs4iIi8RaM+8DFGAeggKM\nAqqBx+p8ZhuQFXNkIiKpKQj81K6TtQmfMBNoB6wHzrDr5CIikjhXY3q0bMO0zEVERERExEmRDB4q\nDP/+Q+Acm+JyQkvX4k7MNdgAvAv0si8020U6qOwCoBK40Y6gHBLJtfABHwCbgIAtUTmjpWvRBViK\nKd1uAn5pW2T2mgXsBTY28xlb82YapsySCbSl8dr5NcDr4fcXAausDsohkVyLi4EO4feDSO1rUfO5\nt4AlwE12BWezSK5FR2Az8JPwz13sCs5mkVyLAmBy+H0X4GtsmqrbZpdiEnRTyTzqvBnv3CyRDB66\nFpgdfr8a8we3a5zndaNIrsXfgG/D71dT+5fXayIdVJYLzAe+tC0y+0VyLe4AXsWM1wD4yq7gbBbJ\ntfgc+HH4/Y8xybzxSZiS20rgm2Z+H3XejDeZRzJ4qLHPeDGJRTuQ6l5q/+X1mkj/XFwHTA//7NWx\nC5Fci38FOgErgL8Dv7AnNNtFci1mAD2B3ZjywnB7QnOdqPNmvLcvkf4FbNif3Yt/caP5b7ocuAfo\na1EsTovkWjwFjAx/thWxj3lwu0iuRVvgXOBKoD3mDm4Vpl7qJZFci9GY8osPM06lBDgb2G9dWK4V\nVd6MN5l/BpxU5+eTqL1VbOozPwnv85pIrgWYh54zMDXz5m6zklkk1+I8zG02mNro1Zhb70WWR2ev\nSK7FTkxp5UB4K8UkMK8l80iuxSVAzYozQWA7cBrmjiWV2J43Ixk8VLeQ3wfvPvSL5FqcjKkZ9rE1\nMvtFO6jsT3i3N0sk1+J04E3MA8L2mIdiZ9oXom0iuRZTgXHh910xyb6TTfHZLZPIHoDaljcbGzw0\nNLzVeDr8+w8xt5Ne1dK1mIl5oPNBeFtjd4A2iuTPRQ0vJ3OI7Fo8hOnRshHIszU6e7V0LboAizG5\nYiPm4bAXvYR5LnAQc2d2D6mbN0VERERERERERERERERERERERERERERERESkMf8PFXJlbhluER0A\nAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10b1d6fd0>"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext version_information\n",
      "\n",
      "%version_information numpy, scipy, functional"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 26
    }
   ],
   "metadata": {}
  }
 ]
}