Kyle Cranmer added file figures/james-stein/JamesStein.v3.ipynb  almost 9 years ago

Commit id: 6211b1f6f2af4fc99c221604937063724d690238

deletions | additions      

         

{  "metadata": {  "celltoolbar": "Slideshow",  "kernelspec": {  "display_name": "Python 2",  "language": "python",  "name": "python2"  },  "language_info": {  "codemirror_mode": {  "name": "ipython",  "version": 2  },  "file_extension": ".py",  "mimetype": "text/x-python",  "name": "python",  "nbconvert_exporter": "python",  "pygments_lexer": "ipython2",  "version": "2.7.6"  },  "name": ""  },  "nbformat": 3,  "nbformat_minor": 0,  "worksheets": [  {  "cells": [  {  "cell_type": "markdown",  "metadata": {  "slideshow": {  "slide_type": "slide"  }  },  "source": [  "# A Quick Demo of the James-Stein Estimator\n",  "\n",  "Kyle Cranmer, June 25, 2015\n",  "\n",  "[![](https://i.creativecommons.org/l/by/4.0/88x31.png)]( https://creativecommons.org/licenses/by/4.0/)"  ]  },  {  "cell_type": "markdown",  "metadata": {  "slideshow": {  "slide_type": "slide"  }  },  "source": [  "## The Problem\n",  "\n",  "\n",  "Consider a standard multivariate Gaussian distribution for $\\vec x$ in $n$ dimensions centered around $\\vec\\mu$\n",  "\n",  "\\begin{equation}\n",  "f(\\vec{x}|\\vec{\\mu}) = \\prod_{i=1}^{n} \\frac{1}{\\sqrt{2\\pi}} \\exp\\left(-\\frac{(x_i- \\mu_i)^2}{2}\\right)\\;.\n",  "\\end{equation}\n",  "\n",  "\n",  "We want to compare the mean squared error\n",  "\\begin{equation}\n",  "E[||\\bar{\\vec x}-\\vec\\mu||^2])\n",  "\\end{equation}\n",  "\n",  "of the sample mean (an unbiased estimator)\n",  "\\begin{equation}\n",  "\\overline x_i = \\frac{1}{m} \\sum_{j=1}^m x_{ij}\n",  "\\end{equation}\n",  "\n",  "to the James-Stein estimator\n",  "\\begin{equation}\n",  "x_{JS} = \\left( 1 - \\frac{n-2}{||\\bar{x}||^2} \\right) \\bar{x} \n",  "\\end{equation}\n",  "\n",  "\n",  "(where $i$ is the component of the vector and $j$ is an index over the samples).\n"  ]  },  {  "cell_type": "code",  "collapsed": false,  "input": [  "%pylab inline"  ],  "language": "python",  "metadata": {  "slideshow": {  "slide_type": "slide"  }  },  "outputs": [  {  "output_type": "stream",  "stream": "stdout",  "text": [  "Populating the interactive namespace from numpy and matplotlib\n"  ]  }  ],  "prompt_number": 1  },  {  "cell_type": "code",  "collapsed": false,  "input": [  "def jamesSteinEstimator(x):\n",  " return x*(1. - (x.size-2.)/sum(x*x))"  ],  "language": "python",  "metadata": {  "slideshow": {  "slide_type": "fragment"  }  },  "outputs": [],  "prompt_number": 2  },  {  "cell_type": "code",  "collapsed": false,  "input": [  "def RMS(x,mean):\n",  " return sqrt(sum((x-mean)*(x-mean)))"  ],  "language": "python",  "metadata": {  "slideshow": {  "slide_type": "fragment"  }  },  "outputs": [],  "prompt_number": 3  },  {  "cell_type": "code",  "collapsed": false,  "input": [  "def doRun(nDim,nSamp=1000):\n",  " mean=zeros(nDim)\n",  " data = zeros(nDim*nSamp).reshape(nDim,nSamp)\n",  " for i in range(nDim):\n",  " mean[i] = 0.5 #set value of the mean's components here \n",  " data[i,:] = random.normal(mean[i],1,nSamp).T\n",  " data=data.T\n",  "\n",  " avRMS_js = 0\n",  " avRMS_mle = 0\n",  " for x in data:\n",  " avRMS_js += RMS(jamesSteinEstimator(x),mean)\n",  " avRMS_mle += RMS(x,mean)\n",  " avRMS_js /= nSamp\n",  " avRMS_mle /= nSamp\n",  " return (avRMS_js, avRMS_mle)"  ],  "language": "python",  "metadata": {  "slideshow": {  "slide_type": "fragment"  }  },  "outputs": [],  "prompt_number": 4  },  {  "cell_type": "code",  "collapsed": false,  "input": [  "nDimToTest = linspace(2,20)\n",  "av_js = zeros(nDimToTest.size)\n",  "av_mle = zeros(nDimToTest.size)\n",  "for i,nDim in enumerate(nDimToTest):\n",  " [av_js[i], av_mle[i]] = doRun(int(nDim))"  ],  "language": "python",  "metadata": {  "slideshow": {  "slide_type": "subslide"  }  },  "outputs": [],  "prompt_number": 5  },  {  "cell_type": "code",  "collapsed": false,  "input": [  "plt.plot(nDimToTest,av_js,c='b') #James Stein in blue\n",  "plt.plot(nDimToTest,av_mle,c='r') #MLE in red\n",  "plt.ylabel('Mean Squared Error')\n",  "plt.xlabel('Dimensionality')\n",  "plt.legend(('James-Stein', 'MLE'), loc='upper left')\n",  "plt.savefig('james-stein.pdf')"  ],  "language": "python",  "metadata": {  "slideshow": {  "slide_type": "fragment"  }  },  "outputs": [  {  "metadata": {},  "output_type": "display_data",  "png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEPCAYAAABCyrPIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xe8HFX5x/HPpAcSCAkllMildwRpiQS4UoQgYgHpSlNA\naYrATxEkICIKioKFpgIqCCpVAwhI6J0k1ISWhCSUAKGEIiV5fn88Z925m9l7d++d3Zmd+32/XvPa\nvbuzs8/e3T1n55TngIiIiIiIiIiIiIiIiIiIiIiIiIiISC70BSYB1yfc1w68Fe6fBJzYvLBERCSu\nXxOe42jgSWBolftvB3ZtQhwiItKJPg0+/krAzsBFQFRln2q3i4hIEzW6QjgbOA5YWOV+Az4NTAEm\nAOs2OB4REcnALsBvwvV2kvsQhgKLhevjgKcbH5aIiCRpZHPN6cBXgY+BQcASwD+Ar3XymOnAJsC8\nitufBVZrQIwiIkX2HLB61kFU2obkM4TlKFdKmwMzqjzeGhBTT43POoAE47MOoIrxWQeQYHzWASQY\nn3UACcZnHUCC8VkHUMX4rANIUFfZ2YxRRiWlwA4Nl+cDuwPfxM8i3gP2amI8IiIS06wK4fawgVcE\nJb+h3M8gIiIZavQooyKbmHUACSZmHUAVE7MOIMHErANIMDHrABJMzDqABBOzDqCKiVkHEGcwMOsY\nGiWPfQgiIrlkMMjgHuosO1tlUpiRHOs8YKkmxyLpewMYnnUQIkVgXlb+BegbwR60Tjlfs2q1nM4c\nikHvo0hKDE42uM9gMAX9bqlCKDa9jyIpMNjbYIbByPJNxaMKodj0Por0kMEYg7kGG3a8uXhUIRSb\n3keRKgx2NrjUYEur0h9g0GbwosHnFr2reFQh5N9WwNRuPlbvo0gCgzXDr/7TDJ42mGRwsJVzwGGw\nhMFj5ksNJByieFq1QpgIHJx1EJ0Yiw9NexN4HbgL2DTcdwBwZ5PiyPv7KNJ0BkMMnjA4JPzdx2BH\ng+sNXjU402ANgwkGv6ty9lDI71arVgi3AQdlHUQVS+AVwZ74B2kQsAOwQbj/AFQhiGTCIDK40uCi\npILeYNVQIcwzuMmgf/VDFU+rVwjDgH8Cc/G5E9cDK8b2mwj8CLgbmA9cByyNjyV+C3gAWDm2/9rA\nzfiv+qnAV2L37Qw8AbwNzAa+WyW2TfHx/0nWAd7Hc0zNp5x9diBwFjATeBn4HV6RgKc4nxU7xozw\n3FPwiuevVJ85mff3UaSpDI41eNDK369q+w2o1q9Q3qV4Wr1CGA58CX9zhwBXAlfH9puIrwWxCv7L\n/QngGWBbfE3qS4A/hH0Xxwve/fHUIxsBr+KVBMBLwJbh+pLAxlViGwq8BlwM7MSiE/z2Z9EzhLOB\na/AKbghecZ0e7munY4UwHbgPH/62FL6M6qEky/v7KNI0BtsavGTwiXQOVzzdrBDM0tm6rVqT0UZ0\nXPPhNuD7sb/PAv4V+3sXYFK4vidwR8Xxzgd+GK7PxNscl6ghvrWBP+IF+UfAtcCy4b4D6FghRMA7\nwKqx28YAz4fr7SxaIewT+/un+BlFkkJ+aEXqZTAqVAbbpXfI2hU8uV0UpbP12GC80J6BNwHdjv96\njx/7ldj1/+LNS/G/h4TrKwNb4M09pW0ffG0JgN3wZqMZ+JnH6HD7DXjzz3xg73DbVOBAYBSwPrAC\n8Msqr2EZfGTDw7HnvQFv2qrm5dj192OvQUQqhOahfwC/iODWLGJo5noIvVUEHAusiS8CNBc/Q3gk\n3JdUg3dWq7+AVyifrXL/Q8AX8aamI/HmqU/gS5R2ZhreNHVIlRhewwv1dfFmqXrpLEB6NfMfVZ8D\nFuBn25XbD/Az/LOyilEVQnMMwQvTt/D+hJMT9omqXK/0L+AMYD/ginDbRvgv/+fwZFb/DM81H//w\nJVkL/3BeAczBzxL2Bu4N978CrISPXvgIWAhciJ9BHIH3W6wIrAf8u5N4a3lNIoUUOny3Bg7D++pu\nBj7Ay4TKbTawc5Thj6eCNxnlwkK8EB2M/8q+B29qqXzTreJ6tfvn42cHe+EF+UvAT4AB4f798Pb7\nt/Bf+/tWiWs+3vR0P/7r5F7gUcqjkm7FO7dfptx89X/4+tb3hePfjJ/5JL2GSkmvSaSQDJYKE8We\nxPvO7gVWjWCPCL4awZci2CGCMRFsEMEqEWwV+fdSutCqo4weBnbNOogWkPf3UaQmBosZ/NbgTYPL\nDLbuYlhoE0IqnlasENbD14kelXUgLSDP76NITQxGGjxg8Bcrj9bLWiG/W61WIfwUbw88IutAWkRe\n30eRmhhsaDDT4KSMzwgqFfK71WoVgtRH76O0LPNspHPN+/XyppDfLVUIxab3UXLBPIHcaIMzDJ4y\neNLgdIPNLWEQjsGR5mmnx2QRbw0K+d1ShVBseh8lMyEf0I7mGUNfNHjcPN30pqEiOD1UDHNCh/Fn\nQ+fxr8O+bVm/hk7k7rvVF0+7cH2V+8/B8/ZMoXreHVUIxab3UTJhcIB5xtB7DI4zWKOTfdcyON7g\nboMPDG40zziQZ7n7bh2DZ+28LuG+nYEJ4foW+Pj2JKoQik3vozSdwX4Gs82z+9b72KXMf+zmXa6+\nWysBtwCfIfkM4Tw8WVvJVMo5eeJUIRSb3kdpKoOvhOahdbOOpcFyldzubOA4fLZukhXpmCFzNl6J\niIg0hMEXgHOBnSKfSSxBI3MZ7YKnPJiEp0aupnLMbrUabXzs+sSw5dUMYHk8e+jrsdsnAZ/E1z04\nBa8MT0p4/EJ8Ulv8f3EKGSa9EikC82bqC4FxkadqKZp2Oi9vM3M6XuBNx/PtvAtcWrHPeXQcu1uU\nJqPpwFN0nJi2Af76FuAprP8InFrl8QvpuO5A0eX1fZQCMdguzBcY3fXehZHL79Y2JPchxDuVR1Oc\nTuXpeCrbB2K3nQWcgBf2pQrhR1UerwpBJEUhp9Bcg62yjqXJctWHEFcK7FDKSylOwFfcehZfQOZb\nTYyn0e7DVy1bGx+NsCfw5zoen6fp7yItyaDNvLn178De0aLLwkoL6tYZgqW0hmY34p2OL4H3A7zp\nbCfgJrxiqPUM4S06roq2QzfiaBU6Q5DUGAwy2MvgZoPXDM4xTzbZG9X13Sr0AjlR9mln/4T/IlkF\n7z+pJ56NKa9XLCJdMNgQ+Dq+pOwk4PfANZEvQSs1KHSFkAMv4IX6OOCghPv1y1ikh8xTTf8MP4u+\nCNg08pF+UidVCI13MDAMX0Iz/v+Owt+DYrctwJerLN0vIlWEZHPfwJteLwXWinz1Pym4VhxltG3C\n7f3oOOx0YcV2R9hvIf7Bnh/bftHYkDOV1/dRcsrgUwb3G9xlPqRbkhXyu9VqFYLUR++j1MRgSYNz\nDV42ONC0LnxXCvndUoVQbHofexmDVQ3WqvMxWxvMMrjAYESjYiuYQn63VCEUm97HXsRgsMFUgzcM\nTjTo38X+fcN+Lxns2Kw4C6KQ3y1VCMWm97EXMTjb4HKDUQY3GDxisFGVfUca3GIw0Tw3mNSnkN8t\nVQjFpvexlzBoN195bET4OzJfpGauwakGA2L7bhf2PdU0IrK7CvndUoVQbHofewGDoQbTzTMhV963\ngsF1Bo8ZbGFwivl6BdtnEWuBFPK7pQqh2PQ+9gKhM/j3ndwfGexj8HpoJhrZzPgKqpDfrWoval64\nT1trb/Mq31jJl54WzgbjDGaYJ3zsat/+pomZaelVFYKINJjB+gYLDc4zWKwbjx9uvnZx0mRNaay6\nys5WqYWN1olVpFDM07bPBNrw0UB7RfBYnY+fF8FRjYlQOlHIslNnCCIZMFjNPIX0EuZt/F8zeNXg\n8FqadQx2M3i6O2cWkopClp2FfFEieWdwvlWs22GwhsFDBtcaLJ3wmMhgWYNtzVNMjGlexFKhkGVn\nIV+USJ4ZrGgwr0qhP8DgzNA3cJDBsQYXGdwdHjMvXD8si9jlf9SHICI9ZyHDbgTHdLLPjvjSt9OB\np4Cp4fLVSD/k8qCQZac+WCJNZLB0+JW/YtaxSI/UVXYqdayIJDka+FsEc7IORKSSzhBEmsR8zYHX\nDFbNOhbpMZ0hiEiPfBO4MfL1wEVyR2cIIk1gsFgYKrpe1rFIKnJVdg4C7gcmA08CP0nYpx14C5gU\nthMT9snVixIpKoMjDa7OOg5JTe7KztIMxX7AfcDYivvbgeu6OEbuXpRIqzFYyeCHBscbbG4VawyE\nuQUvGGyWVYySurrKzmYsOvFeuBwA9CU5s2XhxsmK5EFILzEWOBJfW+ByYCFwEbCywb3A7WHbAJga\nwYMZhSs51wf4dArHmAzMB36WcP82wOvAFGACsG7CPjpDEKmD+brFBxtMNpgWmoKWqNhnhMEXzZe0\nfMRggfn3UYoj9bJzckrHWRJvMmqvuH0o5WalccDTCY81YHxsqzyGiAQGB4YEdP8y2MlqHE1oseUr\npWW107GsTL1COAvYnXSadU4Cju1in+nA8IrbdIYgUgODPczXIdYoIYEGlJ3v4G2OH+HNPvOBt2t8\n7NLAsHB9MHAHsF3FPstRrmw2B2YkHEcVgkgXDLY3X6z+k1nHIrmRq7JzA+ARvNnpUeC4cPuhYQM4\nHHg87HMPMDrhOLl6USJ5Y7BZqAy2yjoWyZWGlJ1fAH6ONx99vhFP0AVVCNLrGPQz+HRXfQAGaxm8\nZLBrs2KTlpF62XkGcCtwEHAwcDPJE8waSRWC9DoGxxm8Y/BkGDE0KGGfFc0Xrz8wixgl91IvOx/D\n5w+U9KWO9VRTogpBehWDUSHB3Oph5bEbwlnACRYGXZgvXv+4wfFZxyu5lXrZ+SgwIvb3iHBbM6lC\nkF7F4B8GJ1fctr7BH8M6BeeEFcnOMk3slOpSLzv3BmYCFwOX4KOA9kr7SbqgCkF6DYNxBs8mNRGF\n+1cwOMPg9FrnGEivlWrZ2QfYE1gB71jeFVg+zSeokSoE6RXCDOPnzCdpivRU6mXnw2kfsBtUIUiv\nYDDe4O9ZxyGF0ZBRRscCo/DOrNLWTKoQpPAM1ggdyaOyjkUKo66ys5bOqBkJBzWau7yeoY4zKbDQ\nMXwDcEvk831E0pBq2VnqQ8iazhCk0Ax2D0NI+2cdixSK+hBEWonBUINZBltnHYsUTupNRmcArwFX\nAO/Gbk9a6KZR1GQkLStMJNsEX5Dqvwnb0cBSEeyfWZBSVHWVnd3tQwBYpdYnSYEqBGkJ5jP51wXG\nxLYV8CSP7+NzC0rb4HD5DrBdBK9kEbMUWiHLTjUZSa4ZRAY/N3jL4BmDSwwOM/ikdUz9ItJMqZWd\n8fwoX6m47/S0nqRGqhAk1wwOCp3Cy2Udi0hMamXnpCrXk/5uNFUIklshx9CrButkHYtIhbrKTuVB\nEekBgyHA34DvRvBU1vGI9IQqBJFuCpPJfgvcG8GlWccj0kgLKK+h/HHseunvZlKTkeSOwYEGTxgs\nnnUsIlUUsuws5IuSbPVkVrDBeqHfYN00YxJJWSHLzkK+KMmOwdYGHxhcZ7BP6Auo9bGLh2UtNZFM\n8q6QZWchX5Rkx+B6g+8a7GfwrzB/4AqDL1VbmCb22IvNF4wSybvUU1fkQSFn20k2DNYGbgfaIp89\njPnSsLvhqwFuDDwAfAQsDNuCcDkYz/S7WdQxlYtIHqWeuiIPVCFIagzOA16JKtYsjt2/ArARPgqv\ncusL3BXBnCaFK9ITqZWd79BxZFF8e7uGxw8C7gcmA08CP6my3znAM8AU/JdZEjUZSSoMljF4w2DZ\nrGMRaYLUy87TgG8BS4Ttm8CPanzsYuGyH3AfMLbi/p2BCeH6FmGfJKoQJBUGPzS4MOs4RJok9bLz\n0Rpv68xiwIMsOkTvPDouwDOV5FwwqhCkxwwGGbysoaLSi6SeuuJdYD+87bQvsC/enFTr8SfjaX1v\nw5uO4lYEZsX+ng2sVOOxReq1L/BItOjnUETwppyu7AP8Cvhl+PvucFstFuKdc0sCNwHtwMSKfSo7\nPKrVaONj1ycmHEekqpBm4hjgqKxjEWmg9rDl3knAsRW3nYcP8ytRk5E0hME4gymm0WrSu6Redq4F\n3Ao8Ef7eEDixhsctDQwL1wcDdwDbVewT71QejTqVpUEMbjb4WtZxiDRZ6mXnHfgIoNIaCBHlyqEz\nG+DLBk7GO6GPC7cfGraSXwPP4sNOP1XlWKoQJFFYqex0g52r/foPq5bNMRjQ7PhEMpZ62flQuIwv\nijM57SfpgioESWSwt8HUsN1ksH7CPpcYfC+L+EQylnrZeQOwOuUKYfdwWzOpQpBFGAwzeNFgtEF/\ngyMN5hqcV5p4ZrCCwTyD4VnHK5KB1MvO1fA+hPeAF/FRRm1pP0kXVCHIIgx+E9JQxG8bbnC2wWsG\nx4eF78/NKkaRjKVadvYFzgrXh+AzlbOgCkE6MNjc4CWDparcv6bBtQYfm/+oEemNUi877yP7oXqq\nEOR/DPoZPGI+YbKrfZdpRkwiOVVX2VnLxLTJwLX4QuLvxZ7kqvriEknN4cCbwF+62jGCVxsfjkgx\n1PLL/+JwWVnTHJhuKJ1S+msBwDzdyRRgywimZR2PSM4VsuxUk5EAYPA3g1OzjkOkRaTeZDQYOBjP\nEDk49gQH1ReXSM+Yz2zfGM04FmmIWrKd/gnPL7QTnlBuFLVnOxVJhXkK9V8Dh5eWvRSR5ivNSi6t\ngdAfXwmtmdRkVBAG2xncaHCE1ZDq3GCAwWdCU9FfmxGjSIGkvh7Ch+HyLTw/0TA0lE+6IeQa+hnw\nMLAZnn30AYPvG6wT22+UwSEGVwNzgTPw/FnfyiJuESn7Bj7tfxtgOj6M77Amx6AzhAIw2MXgUQs/\nREK6iW0NzjWYFfIRPWbwqsFfDPbTPAKRHilk2VnIF9WbhKykD5rnwqp2/6YGW5jPkBeRnkt9lNHJ\nVQ6uoX9Sj3HAIKpMaIz8s/VQ0n0i0hy1VAjvUq4IBgO7oDVppQ6h7+Bk4EeRL6sqIgUxELi9yc+p\nJqMWZrCTwRNW2yAGEUlPw8vO4fgKZ82kCqFFhb6Bew32zDoWkV4o9T6Ex2LX++ALj6j/QGr1WWBJ\n4O9ZByIiPdcW21bCJ6Y1m84QWlA4O7jHYO+sYxHppVI/Q3i74u+hFX/Pq+cJpVfZHl/A5sqsAxGR\ndMzAR4a8HraF4bbpwPNNikFnCC0mnB3cZbBv1rGI9GKpl50X4lkmS8YBF6T9JF1QhdBiQs6iaZpk\nJpKpusrOWhZOeBxYv4bbGqmQizy0AoPv43mGPsDzWn0Q2/6LNym+iee6il9eBlwQwZ8zCFtEXF1l\nZy19CC8CJ+Jf7AjYB5hT4/FHAZfiI5MMP7M4p2KfdnyJzlLz0z+A02o8vjSQwXrAUcD/4fNPKrfh\n+GCDJfGkh/HLGSg7qUjhjMAL8Ulh+xVeENRiJLBRuD4EX/JwnYp92oHrujiOmowyYHCBwUlZxyEi\n3dbQsnM4PZtteg2wXcVt7cD1XTxOFUKTGSxt8Ib52Z2ItKbUys6TKf+aHwjchg8xnQvs0I3jtQEz\n8TOFuG3w0UtTgAn4Up2VVCE0mcEJBr/POg4R6ZHUys4nKXdGHIIvn9kXryQerPNYQ/BMll9MuG8o\nvjwi+AimpxP2MWB8bGuv8/mlDmGVsjkGG2Ydi4jUpZ2OZWVqFcKk2PWr6LgoziRq1x+4Cfh2jftP\nZ9E+Cp0hNJHBPgb/yToOEemx1MrO+/AlM5fBm4pWjd03rcZjRPgoo7M72Wc5ymcim+OjUyqpQmiS\n2EI2u2Ydi4j0WGpl52i84J9Hx5EmnwMur/EYY/GZzZMpj1IaBxwaNoDD8XkNk4F7wvNWUoXQJAaf\nNnhWE8pECqGQZWchX1QeGVxpPvdARFpfIcvOQr6ovDFY2eB1WzSBoYi0prrKTq1gJXGHA5dEMD/r\nQEREqtEZQg3CL/xjDbYyX/+6nscOMXjNYJVGxSciTZd6cjuALfGJZaXcR4aPHmoWJbergcEf8AJ9\ncTwP0aN4R/3dwN0RvNLJY78FbB/Bl5sRq4g0RV1lZy07/hkfcjoZWBC7/cj64uoRVQhdMM8b9RSw\neuT9AIsDm+GV+ZbAGOBl4Maw3RHB++GxfcJjvxHBHVnELyINkXrZ+VTaB+wGNRl1weA0g992cn8f\ng00MfmBwp8F8gxsMjjY4zOARy/59FpF0pV52/g1YIe2D1kkVQicMFjd41WCNOh4zzGA3gwsNZhjs\n2cgYRaRZLAJbA+xAGrCm8jJ4XqMH8EVRCE+imaz5sT/eR/BMrQ+IfBGbf4RNRFqW9QM+CWyFTwYe\nC3wE3FnvkWppImivcvvEep+sB9SHUEWYUTwVOCjqxgdARPLMhuJJ6ko/wA3P/mCxbRTwAv79v8u3\naGZs/1RXTJtY68EkE58H3sA/CCJSCBYBuwO/AG7GM0V/iA8AicLWJ2xzIJrXrMjG4Omu38FPQxbi\n6+g2k/oQqggdxGr/FykMWwPsJrDHwMb29GCphBTzMN5ZOQlvnjgQOCPtJ+mCKoQEBpuHDuFazvRE\nJBMWgW0AdizYlWC/BPsm2HZgo8BCxggbDHYK2Gtg3wXrn8aT17NzLW1LDwOb4JOcSgumTKa8VnIz\nqA8hgfki9vdHnacXF5GmsxH4ypI7Ap8F/ouvC3MPnvJ/zdg2DHgWWBJvjfkORLPTCoSU+xDexZfQ\nnAL8DJ/cpMI5YwYr4x+4Q7KORaS4rB+wEp4BoC122YZP/oRyeVi6HBgeMxGvBH4M0bOdPMdQYHWg\nH0T1rkaZqloK9jY85cEA4DvAEvgEqE5eYOp0hlDBvLNpYQTHZh2LSHHYYvjwzR3CtjZe/s3AV3OM\nX8aTQFrsciHwFEQfNiPiLjSk7FwMWCvtg9ZBfQgxBksazDP4RNaxiLQ26wu2Cdj3wG4Fmw92J9gP\nwcaADcw6wh5KvezcFV85bUb4e2PgurSfpAuqEGJCRtPLso5DpLXYsqEj9ztgfwR7COxdsKfAzgH7\nPNgSWUeZstTLzkfwTo9JsdseT/tJuqAKgf+td9xm8IJ5R7+IVGUDwMaFwv9lsHlgt4OdC3YI2OjQ\nfl9kqaeu+AhPcxC3sJ4nke4JyebWBrbG2zW3xvty/hL56C8R6cD6A9sCe+CTuaYCVwKnADP9N5VU\nU0uF8ASwb9h3DXy93XsaGVRvZ17o/x4fsvYunpL6NvxD/WykMyaRwPrg/Ztb4D+avoDn9LoSGA/R\nrOxiK6bFgdOBh8L2Y2BQk2PoVQWgwbcNbjbPUSIiQJjgtTzYLmA/ArsZ7E2w58EuBzsKbOWso8yZ\nhqyYlrVeM+zUYCm8E/8zkZ+diRSc9cH7KUeEbWk85f4n8B9Fo8L1lfChnpOB+4H7gAcgmptB0K0i\ntRXTru/kYM1Of92bKoSzgKERHJp1LCLOVgOOBl7Dm2Oe9suoxpxmNgBfdXGNsK0ZLlfEK4Cl8Fxp\nr4fneB14EZiFZ/GcVd6i91J6Ub1FahXCq8Bs4HK8No7vb8DtNRx/FL728rLhMRcA5yTsdw4wDngP\nOICOI5pKz1f4CsH8S/MAsH7kM8JFMmT98YmP3wXOx7+DpQJ9dbwQfwYvJ/rifV8Dw2Xp+gi84J8V\n9o1VKMzCC/95EH3UrFfVy6RWdvbDC+lL8QL6NHzh9nqMpJzzaAjeFLJOxT47AxPC9S3w08BKvaIP\nweAKgxOzjkMkTMp6DGwCWFvC/RHYimDtYPuAfQXsC2GY57aepdM2B1srnCFINhpSdg7Ef7m/BhzR\ng+NcA2xXcdt5dEzfPBVP/hRX+ArBYIzBLPNZ4SIZsWFgvwN7EWxPL/ilhaU6D2EQ8DlgLzyn0a+A\nq7sVlj9+Y8rNTyWl08mS2Xjn0SvdfJ6WE+Yb/Bw4MfJmM5EGseXx71x/yk07pevL42eo1wPrQfRG\nVlFKNjqrEP6ENxFNAE4FHuvB8wwB/o53TL2TcH/lr5CkWm187PpEirWS227AYODPWQciRWTL4atv\n7Yl/p6fjE04/jF1+CLwP7AHR3RkFKj3XTvVlj3tkIT7EK2mrZ8W0/ngK2G9Xuf88/AykpFc1GRkM\nMHjOFm1KE+kBGw52cBir/wbYn8P4fbXn9y6pNRn16WEg4L/8fw88Cfyyyj7X4f0SfwVG42kyek1z\nEXA48FQEt2YdiLQC64sP1hiFj9VfBh+3H79cBlgN+Df+g2sCRO9nEq60lEZ3GI3F0y48SrmmOoFy\n2ubzw+WvgZ3wNA0H4gn14go57NRgOH5G1B55pSkSY5viTT0rU56gtTwwD+93m4MPD38VH/ARvz4N\novkJB5XeJbV5CHlSmAohdCCPxFdeOgx4L/JLEUL+/a/gZ83L4X15T1OenDUHog+yi09aTGHKzriW\n60Mw6GewicGRBr8xmGDwlMF7Bq8Y3GdwqfmkPen1bBTYj8FeAbsp5Obvm3VU0vJST38tNQg5iEYD\nWwKfBjYDZuKZYR8HbgSeB2ZE3jQmvZoNxCdtboGna94KH2W2NUTTsoxMeq9WOZXI9WmP+XDa04AH\ngbvxSuC+CDSOu1ewAfhIuaPxNcdnU27iKV1/BU/5sEXY1sfTN9yPf16uUpu/NID6EJop9AlMBb4W\nLTrpTgrNhgLfAL6Dt/OfiS81uxLeARy/XB54Dv+M3A88DJHOFKXR6io71WTUcxvj/8cHsg5E0mCD\n8IWJ1sEzbs4J22yIwqRKG4kvFHUIcAvwRYjiK9hNbWLAIqlRhdBzewN/1SpmrcwGAjvgM3l3Aabg\nS5RugKd5CJt9hFcOI4HLgM0hej6TkEV6sVwWtgZ9QkK6DbKORepli4HtDHZxWHz9DrAjQq6fpP0j\nsKXA1gcb0dRQRbovl2VnT+XyRRmMNR9BJLlni4PtAHYa2F1g74DdCXa0p3EWKSQNO22ivfEFhCRT\nNgRvxhkW25YMlyPxocAb4ut63A6cAtyjTl2RjnI5cidB7kYZmVemLwJjIh89Ipmw3fF8PW8Cb4XL\n+PXX8EWX7tPyi9IL5a7sTEPumowMdjQNM82QDQD7JdjzYJtkHY1ITqnJqEnUXJQZGwVciSdy20QL\nuYj0Lrk6QzAYZPCGefph6THrG9bfPQlsPNjWYSho0r6fBXsJ7HiwNFK0ixRZrsrOtOTqRRl8yeC2\nrONobTaxZ8bZAAANPklEQVQSbH+wy8BeA3sC7OdgPwN7EGw+2L/Bvg+2RWgiGg82B2ybrKMXaRG5\nKjvTkqsXZXCl+SxVqYtFoRKYFFbx+jvY10MTUOW+w8B2Df0EU8D+Czax+jwBEUmQq7IzLbl5UQZD\nDd4y0OSkutg6oUB/CGxbsDr7r2yY0kGL1C03ZWeacvOiDPYz+GfWcTSPjQnb8G4+frGQ5/9VsMNV\nqIs0lUYZNdjeeB6bgrN+wE/x1bvmAmuCfYBn9ZwWtqfx7J4zgTcgqvjw2c748qj3AxtC9FKTgheR\nbmiVCQsNmVwRUlevAmwStk3xhW7OBK6MYGHF/iPwRW5WjOCdtOPJD1sGuAL4ENgHonne/s9IPKf/\nWmFbE1/vd2WgD14xlLYVgPWAwyH6d9NfgoiAJqZ1eaDI4HsGNxu8bjDb4FqDHxrsHLYHDR42+KzF\n/pkGh5iPfy8w2wxsZmjmqaN5x4aBfTJ0BB8ZtkGNi1NEapCb5vY0pVkhrBcqgc+Z/+JN2icy2N1g\nmsGt5sthYnCbwZfSiiV/7KDQ1l/g1yjSq6hC6OJA/2ferl3Lvv3DWcFsg2sM5hnk+FevDQi/0tcG\nawtj/Yf5L3Xr5LTRBoD9DmyqjwYSkYJQhdDFge402KnOxywWKpJj04ojfdYWJnQ9Ewr2GWAvg70Z\nxvAb2Efh+nt4+ue3w3yA+WDXgi2Z9asQkVTlqkL4A764+GNV7m/Hs1JOCtuJVfZL5UUZjAhzCHL8\nK7877HNgr4AdU/1MwPqEM4FBYSjoELAl8EVfhnd+BiEiLSpXFcJW+JrDnVUI19VwnLQqhP0Mrk7j\nWPlg/cBOB5sFtmXW0YhI7uRqHsKdQFsX+zTzl+kuFGZSmY3Es60uwDN+zs04IBGRLrVR/QxhG+B1\nfFHzCcC6Vfbr8RlC6CCeV21kUb5YBDYWbBuwTcDWAlvR2/itH54NdA7YqZr5KyKdyNUZQlceAUYB\n7wHjgGvwyU5JxseuTwxbPcYCz0bwcp2PazJbHLgA2AJfkW1IxbY4MA/4KkQ3ZhWliORSe9hyq43q\nZwiVpgNJOXPSOEP4ucEPe3qcxrI1wR4Hu9g7fhP3iXRWICI1ylWnMnReISxHuQ9hczwvTpI0KoRp\nBp/q6XEax74MNhfsEI34EZGU5KpCuBxv9vgQmAUcBBwaNoDDgceBycA9wOgqx+nRizJY02CO5TKn\nh/UDOzPMG9gs62hEpFDqKjtzWEAm6lGCJoNjgLWickXUBDYQOAD4Ht7u/3RsK2UKfQe4BPgA2Bei\n15oXn4j0Akpul/Dg/xh8Pq1guni2QXje/1lgN+BrCSwfRgwdAnYW2HVg08De10ghEWmgXDUZpaXb\nL8pgmMHbBlU6adNig8GOApsN9k+wzWt4TOFqbhHJFVUIFQ/co/ErnNl+YC+GfECbNPa5RERqpgqh\n4oF/MjgszWBiR+8Pdm5IKLdpY55DRKTbVCHEHtTX4DXzyW8ps+XA7gjNQ8PSP76ISI+pQog9aEvz\nIa0ps83BXggdwn3SP76ISCpUIcQe9BODH3WyRz98QZmvg50P9khII3012BFg6y7a8WsH46uKfbE7\nMYmINJHmIcQe9DhwcAT3x25dB58gNwbYCJgNPAg8EC5fxPMebRe2AcB/wrYZ8BngixBN7farERFp\nDs1DCA9YxeAVgz7erGM7g92EryL2Y7BtqWmFMFsV7Btgl4ccQ0t0I34RkSyoySg84Ij3Gfin0PQz\nLTQH7R9mEIuI9AYtlf66ZidzcjSAD9ddnHe3Hsr8zYbx5tojmPeJ/nw4YAF9Fyyg78ely4X0+fhF\nHl3+O5zdB18u8+vAXRAVsrYUEUlDq7Qt2ZsswQL62ixGvTuXZefOY/jzb7HkEx/T760+LBwYtgER\nNrAPCwcA0TkcdepkPvVM1sGLiGSkmH0Ip3HC6lkHISLSYgrZKlLIFyUi0mB1lZ2aVCUiIoAqBBER\nCVQhiIgIoApBREQCVQgiIgKoQhARkUAVgoiIAKoQREQkaHSF8AfgFeCxTvY5B3gGmAJs3OB4REQk\nI1vhhXy1CmFnYEK4vgVwX5X98jhTuT3rABK0Zx1AFe1ZB5CgPesAErRnHUCC9qwDSNCedQBVtGcd\nQIJczVS+E3ijk/t3BS4J1+8HhgHLNTimtLRnHUCC9qwDqKI96wAStGcdQIL2rANI0J51AAnasw6g\nivasA+iprPsQVgRmxf6eDayUUSwiIr1a1hUCLJqaNY/NQyIihdeMPNltwPXABgn3nQdMBP4a/p4K\nbIN3RMc9C6zWmPBERArrOSBXSwe0UVun8miqdyqLiEiLuxx4EfgQ7ys4CDg0bCW/xs8ApgCfanaA\nIiIiIiLSQkYBtwFPAI8DR2UbTgd9gUl4/0geDAP+DjwFPIk3wWXt+/h79xhwGTAwgxiSJkcOB24G\nngb+jf/v8hDXmfj7NwW4ClgyBzGVfBdYiP/vmqlaTEfi/6vHgZ/mIKbNgQfwMuFBYLMmx1StrMzD\nZz01I4GNwvUhwDRgnezC6eAY4C/AdVkHElyCN8kB9KP5hUmlNuB5ypXAFcD+GcSRNDnyZ8Dx4fr/\nAWc0OyiS49qB8si/M2h+XNUmko4CbgSm0/wKISmmz+CFXP/w9zI5iGkisGO4Pg4vnJupWlmZh896\nw1wDbJd1EPhciVvwD2YezhCWxAvfPBmOfyiXwiuo64HtM4qljY5f3qmUJ0CODH9noY3qAy6+BPy5\neaH8TxuLxvQ3YEOyqRBg0ZiuBLbNII64NjrGdDmwR7i+N9m8d3HX4N+3uj7reZiHUKs2vFa+P+M4\nAM4GjsNPofNgFeBV4I/AI8CFwGKZRgTzgJ8DL+ADC97EK9E8WI7y0OZXyOfs+IMoj8DL0hfwCaOP\nZh1IzBrA1vioxInApplG475H+fN+Jt5cmpU2ymVlXZ/1VqkQhuDt40cD72Qcyy7AXLytsBnzOGrR\nDx+h9dtw+S7+Ac3SasC38Q/nCvh7uG+WAVVh5G8y5A/wkXmXZRzHYsAJwMmx2/Lwme+Hn3mOxn+Y\nXZltOAD8Hm+3/wTwHbyfIQtDgH/gZeX8ivvy+FmvW3/gJrxwyYPT8SG004GX8ML30kwj8lPB6bG/\nxwL/zCiWkj2Bi2J/fxX4TUaxtLFok9HIcH158tVkdABwNzCo2cEEbZRj2gD/VTk9bB8BM4BlM4wJ\n4AZ8AmvJs8CIZgbEojG9HbseAW81NRqXVFbW9VnP+xlChNe8TwK/zDiWkhPwTrZVgL2A/wBfyzQi\neBmvpNYMf2+PjzbI0lT8F9xg/H3cHn8f8+A6yh3c++PtrXmwE/6L9wvAfzOOBbzAWw7/rK+CNx19\nCj9DztI1lPsQ1gQGAK9nFw7glVKpktoWH9XTTNXKyrx+1rtlLN5OPxlvopmEf2nyYhvyM8rok/hw\nt6yGLCY5nvKw00sojwpppsrJkQfiHaO3kO1QvKRJm88AMyl/1n+bUUwfUP5fxT1P8zuVk2LqD/wJ\n/1w9TPOzjCZ9pjbF2+wnA/fS/LVdqpWVefisi4iIiIiIiIiIiIiIiIiIiIiIiIiIiGRpAT7G+nF8\nzPUxlNMpbAL8KqO47m7AMS8GdgvXLwTWDtdPaMBziYi0nHh+lmXwNMjjswml4f4IfDnh9socNSKp\nyHvqCpHOvAocAhwR/m6nnI58PD47+g48/86XgbPwrJ034AnSwM8qJgIP4Tn/S3lfJuK54+/H03iP\nDbevF26bhM8KXy3cXkq6GOHZLh8Lz1VKidwejvk3fGGXeHrkk/DFVR4Dzq/yWieGWM/A04FMCsc4\nBU9kVvJj8rWQlIhIwyT9On4DP1top2OFcAe+st2GwHuUFzC5Cs8X1B+4h3JitD3xfDDgC5ycGa6P\nw89EAM4F9gnX+1FOQleKazc8RUCEJ4GbiVcy7XgK8BXCffcAW4bHLBV7LZfiGXWh4xnCbZTXHI//\nD1bG0zeA/8B7tuJ4IjXr1/UuIi3J8DOBBXh/Qx88EyT4L/E2PDHaepTXaeiL56gpuSpcPhL2By/I\nf4AvknQVXgDHjcXTVhueBO52fDnFt/GzgNLxJ4dj3o0nQzsOTzc9PMRba7bamXhit43wiucRvIIU\nqZsqBGl1q+KF/qsJ930YLhfiqZuJ/d0P/6X+BPDpKsf+IFwuoPxduRxfmGUXfAGbQ+m4XKKx6JoB\npRz0H8RuW4BXQIPwtOCbAHPwtQfqTX19EZ5gbTmyy8MvBaA+BGllywDn4c04lWpZyGVaOMbo8Hd/\nYN0uHrMqvjbAucC1+JoBcXfiTU99wrG3xs8MqsVTKvxfxxc3+UoNcX9Exx9zV+OZLTelfBYkUjed\nIUirKXWo9gc+xtvcfxHui68IVbk6VOVKUYYXrLsD5+Dpwvvhy6MmrdtQevwewH7hsS/hnbjx+68G\nxuAdzoY3Bc3FFzxPWq3qTXxI6eP4uha1LBF7Ad5h/TC+8NBH+Locb1R5DhER6SX64JXkal3tKNIZ\nNRmJtLZ18YV1bgGeyzgWERERERERERERERERERERERERERERKb7/B+VlaO/qM+tNAAAAAElFTkSu\nQmCC\n",  "text": [  ""  ]  }  ],  "prompt_number": 6  },  {  "cell_type": "code",  "collapsed": false,  "input": [],  "language": "python",  "metadata": {},  "outputs": [],  "prompt_number": null  },  {  "cell_type": "code",  "collapsed": true,  "input": [],  "language": "python",  "metadata": {},  "outputs": [],  "prompt_number": null  }  ],  "metadata": {}  }  ]  }