{ "metadata": { "name": "", "signature": "sha256:2c31d031d4a34a4da1a6cb54e920e9421e2f23406dd76356da14eabee3cac521" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Statistical Computing Languages

\n", "\n", "

Python

\n", "\n", "

Statistical Computing Languages, RSS

\n", "
\n", "\n", "

\n", "

[Neil D. Lawrence](http://staffwww.dcs.shef.ac.uk/people/N.Lawrence), [@lawrennd](http://twitter.com/lawrennd)

\n", "

University of Sheffield

\n", "

21st November 2014

\n", "\n", "

\n", "
" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# rerun Fernando's script to load in css\n", "%run talktools" ], "language": "python", "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "The history saving thread hit an unexpected error (DatabaseError('database disk image is malformed',)).History will not be written to the database.\n" ] }, { "html": [ "\n" ], "metadata": {}, "output_type": "display_data", "text": [ "" ] } ], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
My Python History (Personal)
\n", "\n", "1. Background as an Undergraduate Mechanical Engineer (1991-1994)\n", "2. Learnt BBC Basic (1983-1990), Lotus 1-2-3 (1990-1991), Quick Basic & Assembler (1991-1993), MATLAB (1994)\n", "3. Worked as Field Engineer on Oil Rigs, tried learning Borland C++ (1994-1996)\n", "4. PhD Course: Learnt C and MATLAB (1997-2000).\n", "5. Started as a lecturer in Computer Science at Sheffield 2001\n", " 1. Felt obliged to learn Java\n", " 2. Read Bruce Eckel's book \"Thinking in Java\"\n", " 3. All scripts in Book were read using Python.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
Why Python?
\n", "\n", "
\n", "[From Guido van Rossum's blog](http://python-history.blogspot.co.uk/2009/01/introduction-and-overview.html)\n", "
\n", "

\n", "Python is currently one of the most popular dynamic programming languages, along with Perl, Tcl, PHP, and newcomer Ruby. Although it is often viewed as a \"scripting\" language, it is really a general purpose programming language along the lines of Lisp or Smalltalk (as are the others, by the way). Today, Python is used for everything from throw-away scripts to large scalable web servers that provide uninterrupted service 24x7. It is used for GUI and database programming, client- and server-side web programming, and application testing. It is used by scientists writing applications for the world's fastest supercomputers and by children first learning to program." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
Ease of Use
\n", "\n", "Python has the style of a scripting language (cf `perl`, `bash`)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Open my web page and look for links (modified by code on Guido's blog)\n", "import re\n", "import urllib\n", "\n", "regex = re.compile(r'href=\"([^\"]+\\.html)\"')\n", "\n", "def matcher(url, max=10):\n", " \"Print the first several URL references in a given url.\"\n", " data = urllib.urlopen(url).read()\n", " hits = regex.findall(data)\n", " return hits\n", "\n", "hits = matcher(\"http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/index.html\")\n", "for hit in hits:\n", " print hit" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "/people/N.Lawrence/inaugural.html\n", "/people/N.Lawrence/lectureNotes.html\n", "/people/N.Lawrence/research.html\n", "http://ml.dcs.shef.ac.uk/sitran/talks.html\n", "/people/N.Lawrence/software.html\n", "/people/N.Lawrence/miscellany.html\n", "/people/N.Lawrence/notes.html\n", "./thematic08.html\n", "kappenball.html\n", "inaugural.html\n", "http://ml.dcs.shef.ac.uk/gpss/past_meetings.html\n", "./inaugural.html\n", "http://maths.dept.shef.ac.uk/maths/staff_info_452.html\n", "http://staffwww.dcs.sheffield.ac.uk/people/N.Fusi/index.html\n", "http://staffwww.dcs.sheffield.ac.uk/people/A.Damianou/index.html\n", "./3PhaseData.html\n", "./inaugural.html\n", "./com4509_6509_2012.html\n", "../A.Cowling/campus_only/teach/com3310.html\n", "http://ttic.uchicago.edu/~rurtasun/tutorials/GP_tutorial.html\n", "./mlss2012.html\n", "http://www.kyb.mpg.de/nc/employee/details/bs.html\n", "./matweave.html\n", "./thematic.html\n", "./licsb.html\n", "./com4509_6509_2011.html\n", "campus_only/lectureNotes.html\n", "./interspeech_tutorial.html\n", "http://www.sheffield.ac.uk/dcs/postgrad/resdegrees/newphdhome.html\n", "./icml_tutorial.html\n", "./thematic08.html\n", "./thematic.html\n", "http://www.cs.ucl.ac.uk/people/D.Barber.html\n", "http://www.well.ox.ac.uk/~mtitsias/publications.html\n", "./chi.html\n" ] } ], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
Well Developed Object Orientated Language
\n", "\n", "* Python has the advantage that it was designed in the mid-late 1990s\n", "* Not in the mid-late 1970s or 1980s (cf `S`, `MATLAB`)\n", "* As a result it has a coherent and functional object system.\n", " * It also has some attributes of *functional programming*.\n", "* It is easily extended and allows for operator overloading.\n", "* It provides an object system without imposing one (cf `java`)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import pandas as pd\n", "import numpy as np\n", "\n", "class my_data(pd.DataFrame):\n", " def pronounce(self):\n", " print \"Extending existing objects is easy, this encourages code reuse and good programming practice.\"\n", " \n", "a = my_data(np.zeros((3, 3)))\n", "a.pronounce()\n", "a.describe()\n" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Extending existing objects is easy, this encourages code reuse and good programming practice.\n" ] }, { "html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
012
count 3 3 3
mean 0 0 0
std 0 0 0
min 0 0 0
25% 0 0 0
50% 0 0 0
75% 0 0 0
max 0 0 0
\n", "
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 3, "text": [ " 0 1 2\n", "count 3 3 3\n", "mean 0 0 0\n", "std 0 0 0\n", "min 0 0 0\n", "25% 0 0 0\n", "50% 0 0 0\n", "75% 0 0 0\n", "max 0 0 0" ] } ], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
Development as a Numerical Platform
\n", "\n", "* Around turn of millenium, several researchers were looking for alternatives for `MATLAB`.\n", "* Python seemed appropriate because of\n", " * interactive shell, \n", " * friendly scripting nature, \n", " * well supported modern programming language\n", "* Some, like me, realised lack of *plotting*, *array types*, *numerical computing*, *fully functional shell* was prohibitive\n", "* Others, like Fernando Perez, John Hunter, Travis Oliphant set about adding this funcitonality.\n", "* Targeted more at *scientific computing* rather than *statistical computing*.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
Recent Additions
\n", "\n", "* More recently Wes McKinney added `pandas`: data anlaysis in python.\n", " * This gives data frame capabilities like in `R`. Proving very popular for data analysis.\n", "* The `scipy` stack catpures a range of scientific computing facilities (FFT, linear algebra etc)\n", "* `scikit-learn` encodes models from a machine learning perspective.\n", "* `statsmodels` encodes models from a statistics perspective.\n", "* The Jupyter notebook and the notebook viewer provide an excellent platform for sharing analysis and code. \n", " * For example my group makes its [work available via GitHub and the nbviewer](http://nbviewer.ipython.org/github/SheffieldML/notebook/blob/master/lab_classes/machine_learning/index.ipynb).\n", " * Jupyter is also multilingual and has plug ins for `MATLAB`, `R`, `Julia` ... \n", "* As a research group we moved from an entire `MATLAB` code base to Python.\n", " * See blog post about it [here](http://inverseprobability.com/2013/11/25/gpy-moving-from-matlab-to-python/)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
Software Distributions
\n", "\n", "* Code is all BSD licensed (cf `R` and `octave` which are GPL licensed)\n", "* Companies deliver optimised python distributions for scientific computing ([Enthought's Canopy](https://www.enthought.com/products/canopy/) and [Continuum Analytics's Anaconda](http://continuum.io/downloads))\n", "* I use Canopy, but we have Anaconda installed (free) on the University managed desktop in Sheffield (free). Commercial firms pay a charge for the full version." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
Python in Action: Regression
\n", "\n", "* Like most modern languages, python is perhaps more about learning the libraries than the syntax.\n", "* Here we use `pandas` to load in data." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import pandas as pd # first we import the pandas library\n", "df = pd.read_csv('../R/data/regression.csv') # now we read in the data\n", "\n", "# gender is stored as male/female. We convert it to 0/1\n", "df['Gender'] = df['Sex'].map( {'Female': 0, 'Male': 1} ).astype(int)\n", "df.describe()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
OIAgeGender
count 101.00000 101.000000 101.000000
mean 6.18495 48.940594 0.198020
std 3.17320 14.413759 0.400495
min 0.98000 -10.000000 0.000000
25% 3.75000 39.000000 0.000000
50% 5.65000 52.000000 0.000000
75% 7.90000 60.000000 0.000000
max 16.50000 73.000000 1.000000
\n", "
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 4, "text": [ " OI Age Gender\n", "count 101.00000 101.000000 101.000000\n", "mean 6.18495 48.940594 0.198020\n", "std 3.17320 14.413759 0.400495\n", "min 0.98000 -10.000000 0.000000\n", "25% 3.75000 39.000000 0.000000\n", "50% 5.65000 52.000000 0.000000\n", "75% 7.90000 60.000000 0.000000\n", "max 16.50000 73.000000 1.000000" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
`matplotlib`
\n", "\n", "* Plotting is done using the `matplotlib` library." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# import a plotting library: matplotlib\n", "import matplotlib.pyplot as plt \n", "# ask for plots to appear in the notebook\n", "%matplotlib inline\n", "fig, ax = plt.subplots(figsize=(10, 7))\n", "ax.scatter(df.Age, df.OI)\n", "ax.set_xlabel('Age')\n", "ax.set_ylabel('OI')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 5, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAmAAAAG2CAYAAADY5Dp/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdcVfXjx/HX5TIvoOIOB85Cc+Q2V2RqWrnSrLShZstM\n+2Zmy0JtmQ1X2jTzV6mVmZkNU6Pck9yjHKmEiqgo3Mu44/cHWGYOEO49917ez8fDR3A4401UvDvn\ncz4fEBERERERERERERERERERERERERERERERERERERERERERcavpwBFgy1nb6gLfAr8BC4A6BuQS\nERER8VttgUb8u4DNBvrkfXwnMMvToURERESMFuDGcy8DTpyzLQ0ok3fdMuf5uoiIiIgUUjX+fQes\nBLCT3CK2A4g0IJOIiIiIoQI9fL3pwGTgXeAR4EP+eST5t5o1a7r27Nnj4WgiIiIil2UPUKsgB5jc\nFOSMauQOtq+f9/lhoDpgAyKAP4CK5znO5XK53BxN3CU+Pp74+HijY8hl0M/Ot+nn59v08/NdJpMJ\nCtip3DkG7Hx+Brrlfdwd+MnD1xcRERExnDsL2CxgJXAVcBAYALwI9AA2ATcBL7nx+iIiIiJeyZ1j\nwO4s4HbxE3FxcUZHkMukn51v08/Pt+nnV7y4ewzY5dIYMBEREfEJvjAGTERERKTYUwETERER8TAV\nMBEREREPUwETERER8TAVMBEREREPUwETERER8TAVMBEREREPUwETERER8TAVMBEREREPUwETERER\n8TAVMBEREREPUwETERER8TAVMBEREREPUwETERER8TAVMBEREREPUwETERER8TAVMBEREREPUwET\nERER8TAVMBEREREPUwETERER8TAVMBEREREPUwETERER8TAVMBEREREPUwETERER8TAVMBERkUv4\n5ptvKF++OiEhEdxwQzdSU1ONjiQ+zmR0gAtwuVwuozOIiIiwdetWWrRoj9U6F6hPUNBzXHvtPn75\nZaHR0cRLmEwmKGCnCnRPFBEREf+QkJCA03kr0BaAnJzxrFhREpfLdeYXr0iB6RGkiIjIRZQpUwaz\neSfgzNuyg/Dw0ipfUigqYCIiIhfRq1cvYmNdhId3Ijj4MSyWm5k69S2jY4mP89b6rjFgIiLiNbKy\nspgzZw4pKSm0a9eOZs2aGR1JvMjljAFTARMREREphMspYHoEKSIiIuJhegtSRESkmMvJyWHt2rVk\nZ2fTvHlzwsPDjY7k91TAREREirGMjAzatevC7t0nCQgIJzIyldWrl1K5cmWjo/k1dz6CnA4cAbac\ns30AsAPYBoxz4/VFRETkEl599XW2b7+C9PTfOHVqFYcP38ngwSOMjuX33FnAPgI6n7OtHvAA0A24\nGnjdjdcXERGRS9i+fS+ZmTdyphI4HDeya9ceY0MVA+4sYMuAE+ds6wJ8CPye93mKG68vIiIil9Cq\nVSMslk8BK+AgJGQ6LVo0MjqW3/P0W5CdyL0Lth74AKjr4euLiIjIWYYNG0KXLtGEhFQiNDSa+vX/\nYPLk14yO5fc8PQg/FChN7oJaHYApQPvz7RgfH//3x3FxccTFxbk/nYiISDETGBjIl1/+H0eOHCE7\nO5vKlStrmaVLSEhIICEhoVDncPff4WrAAqB+3ufjgQTgzBLyfwE1gMxzjtNErCIiIuITfGEi1lXk\njgMzAS2APfy3fImIiIj4NXcWsFnASuBK4CC500/MJ/ex53bgKeBxN15fRERExCt560NePYIUERER\nn+ALjyBFREREij0VMBEREREPUwETERER8TAVMBEREREPUwETERER8TAVMBEREREPUwETERER8TAV\nMBEREREPUwETERER8TAVMBEREREPUwETERER8TAVMBEREREPUwETERER8TAVMBEREREPUwETERER\n8TAVMBEREREPUwETERER8TAVMBEREREPUwETERER8TAVMBEREREPUwETERER8TAVMBEREREPUwET\nERER8TAVMBEREREPCzQ6gIiIiMiFJCUl8c4775GebuW223rSqlUroyMVCZPRAS7A5XK5jM4gIiIi\nBkpKSqJBgxakpd2Kw1EBi2Uys2e/T9euXY2O9i8mkwkK2Kn0CFJERES80tSp75KW1guHYxLwLFbr\ndEaMGGt0rCKhAiYiIiJe6dSpDByOK87acgUZGRmG5SlKKmAiIiLilW67rTsWy0RgEbAJi2Uo/fr1\nMjpWkdAYMBEREfFa8+bN48knX8Rms9K3b29eeSUes9lsdKx/uZwxYCpgIiIiIoWgQfgiIiIiPkAF\nTERERMTDVMBEREREPMydBWw6cATYcp6vDQecQGk3Xl9ERETEK7mzgH0EdD7P9ipAR+BPN15bRERE\nxGu5s4AtA06cZ/ubwJNuvK6IiEiRstvtfPXVV7z77rts27bN6DjiBzy9GHd34BCw2cPXFRERuSx2\nu5327buSmHgCp/NqYBSffvoePXr0MDqa+DBPFjAL8Ay5jx/PuOCcGfHx8X9/HBcXR1xcnLtyiYiI\nXNC8efPYuPEUGRkrADOwioEDe6uAFWMJCQkkJCQU6hzunoi1GrAAqJ/3ZzFgzftaZSAJaA4cPec4\nTcQqIiJeYfLkyTz55HYyM6flbckkICASuz37zAScUsx5+0SsW4AKQPW8P4eAxvy3fImIiHiNNm3a\nYDLNA34D7AQGvkCzZtepfEmhuLOAzQJWAlcCB4EB53xdt7hERMTrNWrUiPfff4uIiA4EBITRqNE6\nvv76E6NjiY/z1vquR5AiIuJVXC4XDoeDwEBPv78m3k6LcYuIiIh4mLePARMRERERVMBEREREPE4F\nTERE5BJWrlxJ/fqtqFixNvfe+xAZGRlGRxIfpzFgIiIiF7Fnzx4aNmxJRsYUoAGhofHceKOZr7/+\nzOho4iUuZwyYXuUQERG5iB9//BGHoztwOwCZmR+wcGE5XC6X5gKTy6ZHkCIiIhcRHh6O2Zx81pZk\ngoMtKl9SKCpgIiIiF9GrVy8qVDhASMhdwGtYLF146aXRRscSH+et9V1jwERExGukpaUxdeo0kpNT\n6NTpem655RajI4kX0USsIiIiIh6miVhFREREfIAKmIiIiIiHqYCJiIiIeJgKmIiIiIiHqYCJiIiI\neJgKmIiIiIiHqYCJiIiIeJgKmIiIiIiHqYCJiIiIeFig0QFEREQkfxYsWMCCBYsoVy6Kxx57lHLl\nyhkdSS6TliISERHxAVOmTGPkyPFYrY8SFLSLsmUXs23bOqKiooyOVuxpLUgRERE/FRVViZMnfwDq\nAxAWdjuvv34dgwcPNjaYaC1IERERf5WVZQX+eeRot5fHZrMZF0gKRQVMRETEB9x22+2EhQ0EEoHZ\nBAXN5uabbzY6llwmPYIUERHxAVlZWQwf/iwLFvxImTKlmTTpJdq0aWN0LEFjwEREREQ8TmPARERE\nRHyACpiIiIiIh6mAiYiIiHiYCpiIiIiIh6mAiYiIiHiYCpiIiIiIh6mAiYiIiHiYCpiIiIiIh7m7\ngE0HjgBbzto2HtgBbAQmAGFuziAiIiLiVdxdwD4COp+zbRFwNdAUCAf6ujmDiIiIiFdxdwFbBpw4\nZ9tPgDPvz4/AdW7OICIiIuJVjB4Ddj+wwOAMIiIiIh4VaOC1nwdOA1+c74vx8fF/fxwXF0dcXJxH\nQomIiIhcTEJCAgkJCYU6R4FW7r5M1ci9y1X/rG39yb37dQOQeZ5jXC6Xy+3BRERERArLZDJBATuV\nEXfAOgMjgHacv3yJiIiI+DV33wGbRe4g+7LkTkfxAvA0EAwcz9tnFTD4nON0B0xERER8wuXcAfPE\nI8jLoQImIiIiPuFyCpjRb0GKiIiIFDtGvgUpIiIiXiwxMZE1a9YQHR3NLbfcQkCA7tsUFT2CFBER\nkf+YMWMmjzwyEpfrFszmjVx3XU2++Wa2Sth5aAyYiIiIFJrT6SQ8vBSZmauBukA2ERFNmTv3dTp1\n6mR0PK+jMWAiIiJSaDabjZycbKBO3pZg4GqOHDliYCr/ogImIiIi/xIeHs6VV9YnIOAlwA6swOlc\nTMuWLY2O5jdUwEREROQ/fvhhLvXqfU9AQCilSt3GnDkfUbt2baNj+Q2NARMREZELcjgcmM1mo2N4\nNQ3CFxEREfEwDcIXERER8QEqYCIiIiIepgImIiIi4mEqYCIiIiIepgImIiIi4mEqYCIiIiIepgIm\nIiIi4mEqYCIiIiIepgImIiIi4mEqYCIiIiIepgImIiIi4mEqYCIiIn7I5XKRmprKqVOnjI4i56EC\nJiIi4mfS09O57rqbiI6uQdmy0fTv/xBOp9PoWHIWFTARERE/M3ToSNauLUt2dio5Ocl88cU2pkyZ\nZnQsOYsKmIiIiJ9ZsWIdWVmDgUAgEqu1P8uWrTM6lpxFBUxERMTP1KgRg9n8c95nLkJCEqhVq6qh\nmeTfTEYHuACXy+UyOoOIiIhP2r9/Py1axJGZWROX6zRVqsDq1UuIjIw0OppfMplMUMBOpQImIiLi\nh9LS0li+fDnBwcG0a9eOkJAQoyP5LRUwEREREQ+7nAKmMWAiIiIiHqYCJiIiIuJhKmAiIiIiHqYC\nJiIiIuJhKmAiIiIiHqYCJiIiIuJhKmAiIiIiHubOAjYdOAJsOWtbJDAfOAB8DUS48foiIiIiXsmd\nBewjoPM52x4mt3zVBg4BD7nx+iIiIiJeyZ0FbBlw4pxtzYEPgSxy75C1cOP1RUSkmFm9ejWDBg3h\noYeGsWXLlovu+8cffzBkyOMMGPAwS5cu9VBCkVzuXoqoGrAAqJ/3+Z/AVUAmYAF2ADHnOU5LEYmI\nSIEkJCRw8819sFqfALIJD5/I8uU/cc011/xn3z179tCoUSvS0+/H5SqLxfIan3zyNj179vR8cPF5\nl7MUUaB7olxQvsPFx8f//XFcXBxxcXFuiCMiIv7ihRfewGp9HbgHgIyMMF59dRKzZ0//z75TprxL\nevpAXK4XAbBar+TZZ8cWSQFLT09n7NhX2bZtDy1bNmTkyOEEBQUV+rziPRISEkhISCjUOTxdwNYB\ndYDEvL+uu9COZxcwERHxHy6Xizlz5pCQsIoaNSozZMgjWCyWQp/XZssCos7aUjpv2/n3dbkq/Wvf\nrKzz71sQOTk5tG3bmR07qpGVdTNLl37K6tUbWbBgzpm7JOIHzr0xNHr06AKfw9MFbA0wEHgy76+r\nPXx9EREx2JNPjmLq1G+wWvsTGrqSTz+dx7p1CQQHBxfqvA8+eCfbtj2B1RoOZGGxPM+gQVPPu+/d\nd/fh//6vF1brlUBZLJah3HffnYW6PsC6dev44480srJmAgHYbL1ZsqQyf/31F5UqVbrk8VJ8uLOA\nzQKuA8oAB4HngWnAJ8AuYCMw0o3XFxERL5OVlcWECW9gtx8AypGZ+T/27r2WxYsXc9NNNxXq3AMH\n9ic7O4dJk57BbDbz3HOv07Vr1/Pu27p1a774YjrPPDMOm83GwIH9ePLJxwt1fQC73Y7JFMo/77gF\nYTIFYbfbC31u8S/eej9Ug/BFRHyM3W4nJSWFcuXKERh4/v+/P336NKVLV8RuT+PMPYDIyO5Mn343\nvXv39mBa97DZbNSp05S//rqFnJwuhIR8TMOGf7J69RKffATpdDo5evQoUVFRhISEGB3Ha13OIHzN\nhC8iIoW2aNEioqKuoEaNhpQtW5lff/31vPtFRkbSokUbgoMfArZhMr1HQMBa2rVr59nAbhIWFsbq\n1Uvo0eMI9euPol+/cH766WufLF/bt2+nSpWrqF69PiVLluWjjz42OpJf8dZ/InQHTETER6SmphIT\nE0tGxlygHbCIyMi7SUr6g8jIyP/sf/LkSR544DFWrFhN5cqV+eCDt6hfv/5/9hNjxcTU4cCBx4H7\ngZ1YLHGsWbOYevXqnXd/h8PBn3/+SWRkJOXKlfNoVqPpDpiIiHjcrl27MJtrkFu+ADphMpVjz549\n592/VKlSfP75DJKSdrJmzWKVLy+Unp7OX3/tBwblbYklIKA9GzduPO/+SUlJxMY2oX7966hcuRaD\nBz+ObqRcnAqYiIgUSuXKlcnO3gMk5W3ZT3Z2EtHR0UbGkkIIDw8nNDQcWJu3JR2Xaz0xMeebOx3u\nuush9u3rhtV6gOzsP5k582c+//xzj+X1RSpgIiJSKFWrVuWFF57GYmlGiRLdCQtryWuvvUz58uWN\njiaXyWQy8dlnH2Gx3EKJEt0ID2/A7bd3uuBYvU2bfsPhGEjuU7hSZGT0YsOG3zya2ddoDJiIiBSJ\nrVu3snv3bmJjY6lbty4Aa9asYfz4qWRn5zB48D107tzZ4JRSEPv37ycxMZHo6GiaN29+wZcJmjS5\njsTEO3G5HgKysVhuZOLEfgwaNOi8+/ubyxkDpgImIiJusXbtWq6//mas1ueBMMLCnmf27Hfo1q2b\n0dH8nsvlYvLkqUyb9n8EBwczZsxwunfv7rbrbd++nXbtbsRur47dnkybNvX59tvPLzgdib9RARMR\nEa/Rt+8gZs26Gvhf3pa5NGv2DmvX/mRkrGJh8uSpPPXUVKzWKcBpLJYHmT9/Jh06dHDbNdPS0tiw\nYQMRERE0bdqUgIDiM8rJFxbjFhGRYsJudwBnLy8UjNPpNCpOsfLuu59itU4E4gCwWv/ko4/muLWA\nlSxZkvbt27vt/P5GBUxERNzikUf68+23fbDZSgEWLJbHeeyxl4yOVSyEhoYAJ//+3GQ6SWho4dba\nlKKlR5AiIuI2P/30E2PHTiQnx86jj95L376FX/BaLu27777jttsGYrWOxGQ6jcUymbVrf/n75Qgp\nWhoDJiIiIgD8+uuvTJ8+i9DQYIYOfVDly41UwEREREQ8TEsRiYiI5JPL5WLnzp0kJiaSnZ1tdBwp\nZi42CH/BRb7mAjSRi4iI+CS73U63bnfwyy+rMZtLUKZMACtWLNLySeIxFytgb5BbtAAqkntrLdnt\niURERNxs6tRpJCSkYbPtAYKx2UYxaNAwvvvuC6OjSTFxsUeQy4AmwEvAOOAV4GWgKbDc/dFERETc\nY+PG7dhs3YEQwITdfhtbt+4o0DmcTidjxrxCzZqNadCgDd9//71bsop/ulgBGw20A+4HqgHVgUFA\nG2Cs25OJiIi4SePGdQkLmw9kAS4CAz+nXr06BTrH6NEvM27c1+zdO5UtW0bQq1d/Vq1a5Za84n8u\nNmL/d+BGYO8522sAi4Ba7gqF3oIUERE3OncMWNmyZpYv/7FAY8CqVq3HwYMzgcZ5W15i6NATTJz4\nulsyi/cq6qWITEDKebanFPQiIiIi3iQwMJCFC79g9+7d2Gw26tatS3BwwWaKDw0NBVL//txsTsVi\nsRRxUvFXFytSs8m9CzbqnO3xQCxwh5syge6AiYiIl/viiy+5996h2GxPEBBwhMjImWzatJqYmBij\no4mHFfVErGWBGcDV5A7Id5E7Jmwb0B84dhkZ80sFTEREvN7SpUuZNWsekZEWhg59mGrVqhkdSQzg\nrpnwI4Gb8j7+DjhdsFiXRQVMREREfIKWIhIRERHxMC1FJCIiIuIDVMBERHyIy+Xi1KlT6CmBiG9T\nARMR8RHff/89JUuWp0yZK6hUqTabNm0yOpKIXCYVMBHxKw6Hg7/++ousrCyjoxSppKQkbrvtHk6f\n/hq7PYPk5NF07NiNnJwco6N5ve+++46QkHKYTIEEBZXmk08+MTqSiAqYiPiPxMRErriiBrVqNaJU\nqfJ8+uksoyMVmc2bNxMY2AhonbelHxkZDpKSkoyM5fVOnTpF1653kJ09GjiF3T6Vu+9+kD///NPo\naFLMqYCJiF9wOp107tyTlJRXsdmOkJm5ggceGMbvv/9udLQiER0dTU7OduBk3pY/cDhOUbZsWSNj\neb2EhASczkhgMGAhdw7xmnzzzTfGBpNiTwVMRPxCSkoKp06lA3fmbalHYGArNm/ebGSsItOwYUMG\nDryT8PDGRETcicXSlgkT3iAiIsLoaF6tRo0awAn+WVkvHTik2erFcJoHTET8Qk5ODiVLlsNm+xlo\nBJzEYmlIQsKXNGvWzOh4RWbFihXs27ePhg0bUr9+faPj+ISWLduzZs1uoBvwI7Vrl2T37o1GxxI/\noolYRaRY+/LLudx770MEBrbAbt/M/ff3ZcKEV42OJV7g1VdfZdWqVTRo0IDRo0cTEKAHQFJ0VMBE\npNjbt28fmzZtomrVqjRu3NjoOCJFxuVy8d57HzJ37g9UqFCaMWOepnr16hw5coTnnnuR/fv/omPH\n1gwfPgyz2Wx03GLFlwrY/cAAIITchb4fO+frKmAiIiJneeGFF3njjS/JyHiKgIBdlCjxLmvWJNC+\n/S0cOXILdntLLJYp3H57PaZPn2p03GLFVwpYaWADUA+wAd8CE4Efz9pHBUxEROQsJUpU4PTpFUAt\nAEJCBtC3r4kvvviL9PQf8vZKw2wuj9V6muDgYMOyFjeXU8AC3RPlomzkhiyZ97mF3FdURERE5AKc\nTgcQ9PfnLlcQdnvmv7ad+bWumxjez4hRiDbgYWA/cBhYAaw1IIeIiIjPePDB+7FY7gC+x2SaQEjI\nfIYPH05Y2BbM5tHAj4SF3catt95BSEgIJ06cYN68eSxcuJDMzEyj48s5jLgDVg6YBtQl987XF8DN\nwMKzd4qPj//747i4OOLi4jwWUERExNuMH/8SZcu+ybx5b1KuXGnGj/+ZunXrsn79r/zvf89x4MAy\nOnRozejRz7Jnzx5atryerKyrgdNcccUo1q79mZIlS17yOnJpCQkJJCQkFOocRowBuxm4m9zpiCH3\nblg1YORZ+2gMmIiIyGW68cZeLF7cHKdzJOAiJGQAw4ZVZty4F42O5pcuZwyYEY8glwFNyR2MHwJ0\nARYZkENExOfY7XZSUlJwOp1GRxEvtn//QZzOtnmfmcjKasvevYcMzST/ZkQBOwW8CMwDlgObgJ8N\nyCEi4lPmz/+GkiXLU6XKVVSoUI0NGzYYHUm8VLt2LQkNnQxkA2lYLB8SF9fC6FhyFk3EKiLiAw4e\nPEhsbGOs1oVAc+ALypR5nOTkvQQFBV3qcK9itVoZPvxZEhJWERNTmbffHkfNmjWNjuVXMjIy6N79\nTn75ZQngpH//Qbz77kStAOAmvjIPWH6ogImInOW7777jzjsncOrUPyM2LJYqbNu2jGrVqhkX7DJ0\n7nwrv/wSSGbmMAICVlC69BR27fqN0qVLGx3N75w6dYrAwEAsFovRUfyar4wBExGRAqpSpQp2+zbg\neN6WnTgcaZQrV87IWAWWnp7OkiXfk5n5CdAap/NJsrLq8/PPhR+J4nK5mDTpbWrVasJVVzVn5sxP\nCh84T0pKCiNGPEO/fvfz2WezfGaerRIlSqh8eSkjpqEQEZECql+/PoMHD2Dq1GsIDGxKTs5ypkyZ\nTHh4uNHRCiQwMDCvvNiAYMAFFM2s7e+99yFPPz0Fq/U9IJuHHx5IiRIR9OjRo1DnTUtL45prWpGS\n0pGcnMbMn/8ie/bsZ9Sopwud+UI2bNjA00+/xMmTp+nbtzvDhj1y5i6L+Alv/WnqEaSIyHmsX7+e\nffv2Ub9+fWJjY42Oc1keeGAon366Hqv1foKDV1C16jo2bVpV6Ds1zZt3ZN26x8id7QhgBrfcsogF\nCz4r1Hk/+ugjHn30GzIy5uVt2U9Y2DVkZJxwSynauXMnTZu2JSNjLFAVi+VZnnqqj1sLnxSOryxF\nJCIil6lp06Y0bdrU6BiF8s47E2jQ4F2WLv2ZmjUr8+yzvxTJYzKLJQxI+ftzkymFiIiwQp83KysL\np7PUWVtKkZOTVejzXshnn83GZhsAPASA1VqZt9++VQXMz6iAiYiIRwUEBDBkyMMMGfJwkZ537NgR\ndO58K1brX0AWFstUnn56SaHPe9NNNzFixPPA+0ADwsLG0qNHX7c9EgwMNGMyZZ+1JZuAALNbriXG\n0SNIERHxG+vXr+ejjz4lMNDMgw8OpG7dukVy3t9++40hQ57myJEUunRpz/jxYwkJCSmSc59r3759\nNGzYkvT0YbhcVbFYxvDyy0MZNmyIW64nhadpKERERPzArl27GDNmPCdPpnPnnd24666+RkeSi1AB\nExEREfEwzQMmIiIi4gNUwEREREQ8TAVMRERExMNUwEREDLR161Y6dbqVa665jlGjxmK3242OJCIe\noHnAREQMcuDAAVq1uoH09Odwuery++9jOXr0GO++O9HoaCLiZnoLUkTEIFOmTGHEiEQyMz/M23KY\n0NCrsNnSDM0lIgWjtyBFRHxIYGAgJpPtrC02zGY9mBApDlTAREQM0qtXL8LDl2M2PwXMxGLpzuOP\nP2Z0rMtis9mIj3+Rnj3v5qWXXiU7O/vSB4kUY3oEKSJioEOHDjFmzDiSk1Pp0aMjAwf2d9sag+7i\ndDpp27YziYmR2GxdCQv7ktatzSxa9LXPfS8FlZKSwrJly7BYLLRv357g4GCjI4kBNBO+iIh43KZN\nm2jd+lYyMnYDZiCbsLDqbN78C7Vq1TI6ntts27aNNm064nA0xuU6Su3aIaxYsYiQkBDmzp3Lnj17\naNiwIV26dDE6qrjZ5RQwDTYQEZFCycnJISAgjH9GtQQREBBS4MeQOTk5fPfdd5w4cYJ27dpRo0aN\nIs9alAYOHEZa2gu4XA8CTnbs6M3kyVPYsGEbCxduJSvrekJCHuPRR1fzyiujC3TurKwsFi5cyOnT\np4mLiyMmJsY934QYRnfARER8iN1u58SJE5QpU4aAAO8Yxpudnc3VVzfnzz9vICenJ8HBs7jyyg0k\nJi4nMDB//5+fnZ1N27ad2b7dhstVE5frRxYu/IK4uDj3hi+E6OirSE6eB9TN2/IGvXtv4LvvVmK1\n7gDCgBSCg2uRlLSHsmXL5uu8NpuNa6/twJ49AUAVYDGLFs3n2muvdcv3IYWntyBFRPzY/PnfULJk\neapUuYqKFauzYcMGoyMBEBwczIoVi+jRI5U6dUbQq1cWv/zyXb7LF8DMmTPZujWQ9PQVZGR8gtX6\nEQMGPOrG1IXXqlULgoMnAg7gGOHh/0ft2tUICooht3wBlCMoKIqTJ0/m+7zvv/8+u3eXJT39V9LT\nPyM9fQqz6Ry3AAAgAElEQVT33eebL2fIhamAiYj4gIMHD9K370Cs1h/IyjpOSsp4bryxBzk5OUZH\nA6B8+fJ8/vkMtm9fxWeffUDp0qULdPzhw4fJzGzCP7+WmpGS8leR5yxK778/kSZN9hIUVIrAwCo8\n8EAXhg8fDuwGPgNOEBDwBlFRIQV6hJiUdBibrSn/3FBpxpEjyUX/DYihVMBERHzAli1bCAxsDDTP\n29IHmw2SkpKMjFVkWrduTWjoZ8AewE5Q0Mtce207j+c4efIkXbveQcmSFalRoyE///zzBfeNiopi\n5cqfOHr0IKdOHefNN1+hTJkyLF26kFq1xhMSUpUGDb7ml1++JygoKN8ZrruuDRbLx8CfQA7Bwa/Q\ntm3bwn9z4lU0BkxExAds2bKFli07Y7VuAUoDuwgJaUZqajLh4eFGxysSkydP5YknnsRuz6Z583Ys\nWDA73+Omikr79l1ZseIKsrNfABKxWAawadMqj7/N+dprbzJq1Cjs9hzatOnA/PmfUapUKY9mkPzT\nNBQiIn7siSeeZdq0/8NsborDsYJJk8Zx3339DU5VtJxOJzk5OYSEhHj82g6Hg+DgUJzOdCD3+hbL\nACZMaMX999/v8TxG/r2QglEBExHxc+vWrWPfvn00aNCA2NhYo+P4FZfLRXh4FDbbWuBKwEVExPV8\n+OFg+vTpY3Q88WJ6C1JExM81a9aMPn36/Kt8ffbZLBo2bEfDhu347LNZbr2+y+Vi0qS3qVevNU2b\n3sD333/v1ut5kslk4vXXx2GxdMBkGoXF0pUaNbLp3r270dHED+kOmIiID/vyy7nce+/jWK3TALBY\nHubjj9+kd+9ebrnehAmTefbZd7FaJwInCAt7hB9//NKvBokvXbqUX375lYoVK9C/f3/CwsIufZAU\na3oEKSJSzNxwQ0+WLu0D3Jm3ZRbt23/OkiXz3HK92NgW7Nr1GnBd3pY3GTDgD6ZPn+qW64n4Ai1F\nJCJSzISEBAOnz9pymtBQ9y0InbvY9D/XM5lOufV64hkHDx5kxoyPycrKpk+f3jRo0MDoSH5Pd8BE\nRHzYihUr6NSpB1brUwBYLK+yaNHXtG7d2i3X+/rrr+nXbzBW6zOYTCcID5/E2rW/UqdOHbdcT9xv\n3759NGrUioyM3jgckYSFvc8PP3zlV4+V3U2PIEVEiqG1a9cyZcqHAAwZch/Nmze/xBGFs3jxYj76\naA5hYSE8/vhg6tate+mDxGs99NAw3n8/Aqfzpbwtn9CixUxWr17Ezp07Wb9+PdHR0Vx//fVnioac\nQ48gRUSKoebNmzNzpntL19k6dOhAhw4dPHY9f+dyuUhLS6NEiRKGLLB+8uRpnM6z72BW4fTpdD7/\n/AsGDHiEgIAbcLk2cfPNLZg9e7pKWBHRNBQiIiL5cOjQIX777TesVmuRnXPLli1Urnwl5ctXJjKy\nLAsWLCiyc+dX3749sFheBVYBW7FYnuSOO7oxYMADWK0/kp4+i4yMDSxcuIalS5d6PJ+/MqqAhQMf\nk7ti6XagpUE5RMTP2O12Dh06RGZmptFR5CIcDgeHDh3CZrMZHSVfRox4jtq1G9Ku3V1UrXoVmzdv\nLvQ5HQ4HHTp046+/niMnJx2r9XvuuGMgf/75ZxEkzr9u3boxceLzVK58HxUq9GDYsC489tgjZGdn\nAtfk7RWGyXQNf/3l3Quk+xKjCtho4ADQIO/PDoNyiIgfSUxMJDq6Jlde2YyoqAp88slnRkeS89i8\neTPR0bW48sqmREVVYMaMmUZHuqglS5YwbdocMjN3c/r0VlJTX6Jnz7sKfd7k5GROn7YB9+ZtaUFg\nYHM2bdpU6HMX1KBBAzl4cDuHD//Byy/HExkZSUzMlZhMEwEXsAmHYzFNmzb1eDZ/ZVQB6wC8DGQC\ndiDNoBwi4iecTiedO/ckJWUcNlsymZkrefDBx/j999+NjlYsJCcn063bndSs2Zhbb72blJSU8+7n\ncrm48caeHD06BpvtMFlZaxk8eAQ7dnjv/4dv27YNu70TUCZvy+3s37+dwr4sVqZMGZxOK7Arb8sp\n7PatVKpUqVDnLSo//DCXGjU+IjAwnLCw65g+fYredi1CRhSwykAoMA1YA4zM+1xE5LKlpKSQlnYa\nuCNvy9UEBrYqkkdFcnFZWVm0atWR77+vzt697/Ltt+Vo27YzOTk5jBv3BlFRlShRojzDhj1Jamoq\nx4+nAHfnHR1LYGA7Q+765FdsbCyBgUuAk3lb5lGlylWFHoweFhbG1KmTCAtrR2Tk7YSHN+Lee3vR\npEmTQmcuCrVq1eKPPzZx/PgR0tOPc8cdWg+zKBnxFmQouaucjgAWA+8CfYB/3YOOj4//++O4uDji\n4uI8FlBEfE/p0qUJCHAAG4HGwAns9kRiYp4zOJn/27x5M6mpAdjtLwEmcnKakpRUmzfeeIuxY2dg\ntS4GLHzwwV2UKvUugYGBZGevA5oBaTidG4iJedzYb+IiOnbsyIABXfnggysJDq6K2ZzMV18VzWD5\ngQP706JFMzZt2kRMzFC3zd9WGJGRkUZH8DoJCQkkJCQU6hxGvUu6AzhzH7MLcA//rKMBmgdMRC7D\n3Llfcc89DxIY2By7fQv339+XCRNeNTrWRW3ZsoX//W8UR4+m0q1bR+LjnyEw0LdmCNq8eTOtWvUg\nI2MXEARkERZWjdatW7J4cXegf96eS2nQIJ7Ro4fTr9+gvJ/TVgYMuI0pU143LH9+7dmzh2PHjlGn\nTh1KlChhdBzxIr40D9jvQAtgHXAzuXfCREQKpVevW2ncuBGbN2+matWqNGrUyOhIF3XgwAFat+5A\nevrzuFx12bNnLEePHuO99yYZHa1A6tWrR7NmdVmzpic2W1cslrnExbUhOroCZvNuHI4ze+6mTJko\nevTozvbt17Bp0yYqV65M48aNjYyfbzVr1qRmzZpGxxA/YdQdsCvJfeQYSm75egHIOOvrugMmIn5v\n8uTJPPnkb2Rmfpi35TChoVdhs/nee0nZ2dm89dZENm3aRZMmVzN06BCSk5PzlrjpiNNpITj4c5Yt\nW+T1xVikoHzpDthuNPeXiBRzQUFBmExnz4NlxWz2rcePZwQHBzNy5Ih/batatSrbtq1nzpw52O12\nevRYXWR3kHJycpgxYwb79x+gZcvmdO3atUjOK+Ip3rqegO6AiYjfO3bsGHXrNuHEidux26/GYnmd\nESPuID7+WaOjeTWHw0Fc3M1s3OjAam1NePhshg+/i9Gj9cKFGEOLcYuI+JikpCTGjn2Nw4dT6d69\nI/3736O19i5hyZIl9OgxnPT0DYAZOExQUA3S0lIJCwszOp4UQ770CFJERIBKlSrxzjsTjY7hU06d\nOkVAQBVyyxdAeUymYKxWqwqY+Awtxi0iIj6lVatW5L5E/ylwkMDAkcTG1qV06dIGJxPJPxUwERHx\nKRUqVGDp0oXUqTOJkiWb067dbhYtmqdHt/mwbds2YmObEhwcTmxsU7Zt22Z0pH9xOp0cOXIEu91u\ndBS389Z/WjUGTEREpAhlZGRQrVodjh17HuiDyfQ5ZcqMYf/+HYSHhxsdjzVr1nDTTb2wWjMxm13M\nmTOTm2++2ehY+XI5Y8B0B0xERKQY2LlzJ9nZUcAgoAQu1yCys6PYuXOn0dHIzMykc+eeHD/+NpmZ\nx8jIWMjtt/cnOTnZ6GhuowImIiL5ZrfbefTRJyhVKpry5Wswdeq7RkeSfCpdujTZ2cnAmYl+08jO\n/ssrxs4dOHAAuz0M6J63pSWBgfW87hFpUVIBExGRfHvuubFMn76BtLQVpKR8yYgR4/jmm2+MjiX5\nUL16dXr37o7Z3AQYitnchN69e1C9enWjo1GhQgVyco6Ru1IhQArZ2TuoUqWKkbHcSgVMRETy7csv\nv8VqfQWoDjTGah3Ol18uNDqW5ENWVhbLlq3E5WoCmHG5mrBs2UqysrKMjkbJkiWZPPktwsLaUKJE\nTyyWxjz++GCuuuoqo6O5jeYBExGRfIuKKgXs5cxqcmbzHsqUKWlopuLC5XIxffoMvvrqR8qXjyI+\n/iliYmLyffz27ds5ftyF0zkbMOF0ujh+/Gp27NiBw+Fg/Pi3yczM5sEH+9GlSxf3fSMXcP/9A2nb\nthVbt26lZs3n/X7NUL0FKSIi+bZ8+XJuvLEnWVn3YDafoESJxWzatJro6GijowFw/PhxlixZgtls\nplOnTkRERBgdqciMGfMK48Z9htU6ErN5FyVLfsS2beupWLFivo7ftWsXjRvfgNW6BwgBsrBYavLJ\nJ5O5664HsFqfASIJC3uBTz6ZzK233urOb8evaCkiERFxu+3btzNv3teEhYXSr18/KlSoYHQkAPbv\n30/z5teRmVkfyKJUqYNs2LCMcuXKGR2tSJQoUYHTp5cDtQEIDb2XceOaMHTo0Hwd73K5uOWWPiQk\nHMdq7Y7FMp/rry9D2bKl+fjj6sDIvD3nc801b5GYmOCOb8MvaSkiERFxu7p161K3bl2jY/zH//73\nHMePD8LhGAVAZuYwnn/+JaZNm2BwsqLhdDrIvXN15vNQHA5Hvo83mUzMnz+L4cOHs27d1zRrVp83\n3niDAQMG/+u8ULDzyuVRARMROUdmZiapqalUrFgRs9l86QPEK/z5ZxIOR/+/P8/Jacn+/fOK7PxW\nq5WTJ09SsWJFAgI8/w7bffcN5IMP+mK1voDJtJOQkHn07Lm2QOd45pnRfPjhtzidN7B587eEhpbk\ngQfu5quvemO1lgNKYLE8zrBhz7jnm5C/6S1IEZGzfPzx/1GqVHlq125CdHQtNm/ebHQkyae2bZsR\nEPAaYANOYTZP4PrrWxTJud98cxJRUeWpWfMaatSox549e4rkvAXL8ApPPXULjRq9TMeOCaxYsYRq\n1arl+/gDBw4wefI0MjLWYrO9R0bGWiZNmkpMTAzz539Kmzaf0azZJKZMeY777hvgvm9EAN0BExH5\n244dO3j44SfIyloD1MFmm0nnzreSlPS71hn0ATk5OcA+IApwAZXIzs4p9HlXrlzJqFGvk529HajK\nwYNv0r17P7ZuXV3ocxeE2Wxm1KinGDXqqcs6/ujRowQHVyEzs0zeljIEB1fh6NGjdOjQgQ4dOhRd\nWLkk3QETEcnz22+/ERh4HVAnb8s9HDt2hLS0tIsdJl5i1apEnM5pwEngFA7HS6xYsbFA50hNTaVT\np55YLFFUrhzLokWLWL9+PQ7HLUBVAJzOR9ixYz2+8LJYeno6mZmZAFx11VUEBBwBPgccwBzM5qPE\nxsYaGbHYUgETEckTExOD07mef5ZqWUdQUBAlSpQo8LlcLpfbfkH7wi9+I9SqVY3AwJ+BUCCY4OCf\nqV07//NkAfTseRcJCdHYbLtJSppEz579CAoKIihoFZCZt9fPlCsX49V3RTMyMrjxxp5ERZUnIqIU\nDz44jPDwcBYvXkClSqMwmYKpVOl5fvrpG7+aqsOXqICJiORp1aoV997bE4ulASVK3IzFchOffvpR\ngQZcu1wuXnxxHOHhpQkJCadfv0FFNtP4ihUrqFTpSgIDg6lbtzm///77pQ8qRiZOfJkKFb4gMrId\nkZHXEhOzjjFjns338Q6HgxUrFpOT8xZQDugEdMVsNtO5cwPCw+tTosRNRETcw+zZH7rr2ygSjz/+\nDL/+GoLdfhKH4zCffLKGd955jyZNmnDo0C5ycrI5dGgXTZo0MTpqseWt9V3zgImIYTZs2EBSUhIN\nGzYs0EzjALNmzWbQoNFYrQuBKMLC7ua+++owefL4QmU6cuQItWrVJz39Q6ATJtO7VKr0Nvv2bSMw\nUMN5z0hPT2f58uWYzWbatm1LaGhovo91uVxERJTGal1J7mNoJxERcUyf/ii9e/dmzZo1pKSk0KRJ\nE6+ZePZCrrqqObt3TwBa5W35kN69l/HFFzMMTOW/LmceMN0BExGv4HK5eOWV16lduykNG7Zl4ULj\n1hds0qQJ3bp1K3D5Avj22yVYrUOBGkAUNtvzfP/90kJn2rhxI2ZzQ6ArEILLNZTjxzNISkoq9Ln9\nSUREBJ07d6Zjx45/l69jx45x6613U6PGNXTpchuHDh0677Emk4mJE9/AYumI2TyS8PAuXHWVie7d\nu2MymWjZsiVdu3b1+vIFEBNTmYCAFXmfuQgOXkGNGpUNzST/pv9tEhGv8PLL43n55TlYrZOAY/Tp\ncx8//vglbdq0MTpagURHlyMoaAs5f798t4Xy5csW+rzlypXDbv8DsAIWIAm7PY2oqKhCn9ufORwO\nrrvuJn7/vRU5OR9y4MACWrXqwK5diYSFhf1n/0GDBhIbeyW//vorFSr04a677iI4ONiA5IUzdepr\ntGx5PdnZCcBpKlY8xTPPvGV0LDmLHkGKiFeoUaMR+/ZN48wiz/AaDz6YxDvvTDQyVoEdO3aMhg2v\n5eTJerhcpTGbF/DLLz/QuHHjQp3X5XLRr98gvvlmPXZ7W8zmb3nuuSE8/fQTRZTcP+3evZvGjTuR\nkbGPM7/yIiOb8MMPk2nVqtXFD/ZxqampJCQkEBwcTIcOHc5bOKVoaCkiEfFZuY+LTvz9eUDACSyW\n/I/f8RZly5Zl27Z1fPXVV2RlZdGly6gCTZZ5ISaTiU8//YBvv/2WvXv30qjRTNq1a1f4wH4uNDQU\nh8NK7huMYYAdp/NUgcaG+aoyZcrQq1cvo2PIBegOmIh4hXnz5nHXXY9gtY4gIOAYEREfkJi4iho1\nahgdTXyYy+Xi1lvvYtGiZKzW2wgL+54mTbL55ZfvDFlOSPzT5dwBUwETEa/x888/83//9yUREWEM\nG/YwNWvWNDqS+AG73c7UqdNYu3YL9evX5rHHhhISEnLpA0XySQVMRERExMM0DYWIiEg+uVwutm3b\nxvr164tsslyR/NIgfA9yuVwkJycDcMUVV3j1MhYiIv4sJyeHrl1vZ9my9ZjNJShd2sWKFYuoVKmS\n0dGkmNAdMA/JzMykU6ce1KzZgJo1G9CpU4+/F0gVERHPmjp1Gr/+mo7V+genT28hKakX9903zOhY\nUoyogHnIqFEvsny5mczMv8jM/Ivly808//yLRscSESmWEhN3YLN1A4IBE3Z7L7Zv32l0LClGVMA8\nZNWqRDIz+5P7L3swmZn9Wbky0eBUIiLFU5MmV2OxzCN3fjAXgYGzqV+/rtGxpBhRAfOQOnVqEBS0\nCHCRuy7Xj9Spo/mNRESM8PDDD3H99WUIC6tBREQsMTHf8cEHvrXqgvg2I0eBm4H1wCFyV5c9m99N\nQ5GamkrLljdw5EgQ4KJCBTurVy+hTJkyRkcT8Qoul4s335zEhx/OxmIJ46WXRnLjjTcaHcvtEhMT\nGTr0WVJSUrnllg688ko8QUFBhT5vWloaQ4Y8yZo1G6lduzrTpr1O1apViyCx/3C5XPzxxx/YbDZi\nY2N9cs1H8Q6+Ng/Y40ATIBLods7X/K6AQe5A/NWrVwPQsmXLYrEUhkh+vfbam4wZM5OMjAnAMSyW\nR1i8+GuuvfZao6O5zb59+2jQoAXp6S8BVxMWFs8dd9Ri+vSphTqvy+WiRYv2bN5ck6ysQZjNiyhf\nfia7d/9GRERE0YQXkb/5UgGrDMwAXiK3iPn9HTARubiaNRuzd+/bwJnCNY6HHkpm2rQJRsZyq8mT\nJzNixGayst7P23KUkJBaZGaeKtR5k5KSqFWrEZmZyeQ+bIASJdowd248HTp0KFxoEfkPX5qI9S1g\nBOA06Poi4mVyH//8UzxMpjRCQ/37kVBwcDABAWeXrTQCAwv/PQcFBeF0ZpE7wBzAidN5Wo/YRLyI\nEROx3gIcBRKBuAvtFB8f//fHcXFxxMVdcFcR8QNjxjxB//4DsVqfxmTKXYx78OAVRsdyq9tuu40X\nXhhHTs5j2O11sVgm8MwzTxb6vOXLl6dHj558++3NWK13ERq6mFq1Iv36ca6IJyUkJJCQkFCocxjx\nCPJl4G7ADoQCJYC5wD1n7aNHkCLF0E8//cSMGZ8TERHG8OFDuPLKK42O5HZHjhzhlVfeIDn5GN27\nd6Rv3zuL5LwOh4PJk6eycuVGYmOrM3LkcMLDw4vk3CLyb740BuyM64An0BgwEbmI3bt3s3fvXmJj\nY6lWrVqBj3c6naxbt45Tp07RtGlToqKiij6kgfSCj4ixLqeAecNakGpaInJB48dP4IUXXiE4uAHZ\n2b/x9ttvMGDAPZc+MI/D4eCWW/qwbNkWzOZKmM27SUj4ngYNGrgxtedoihsR32T0HbAL0R0wEWH/\n/v3UrdsMmy2R3JendxIaei3JyfsoVapUvs4xY8YMhgyZTkbGYnJXoviIevXeZ8uWlW5M7jn33/8o\nH3/sIidnMgDBwUO4554A3n9/ssHJRIoPX3oLUkTkkvbv309wcCy55QsglqCgCiQlJeX7HHv27CUj\n43pyyxdAJw4c2FvESY2zY8decnI6kfvffhPZ2TeyY4f/fH8i/koFTES8VmxsLDk5O4ANeVuWAscL\nNA6sceNGhIfPBVIBF2bz+zRo0KjIsxqlVatGhIbOALKBbEJDZ9Cqlf98fyL+SgVMRDxu79699Op1\nDy1b3siYMa9gt9sBmDFjJq1b30SnTr1YtWoVFStW5JNPPsBi6Uh4eAyRkXcwf/6cAr3N16NHDx54\noCvBwdWxWCpTvfrXzJr1/qUP9BFjxjxHmzYOQkOjCQ2Npk0bB2PGPGd0LHEDp9PJ+PFvce21nenW\n7U527NhhdCQpBI0BExGPSklJITa2ESdPPozT2QiLZTx33nk1TZo04IknXsdqHUfuUkTPsHz5Iho1\naoTVauXw4cNER0df9ht+x48f5/Tp01SuXBmz2Vy035TBXC4XycnJAFxxxRVnxqOInxkx4lmmTl2M\n1ToKk2kXkZGvsWXLOq3x6QV8cRqKC1EBE/FTuYPivyMj4/O8LScIDIwmJqYue/a8BbTL2/4igwen\n8vbbbxmUVMS7REaWIz19HVANgJCQ+3n11at57LHHDM0lGoQvIj4g9z9UjrO2OC64PSBA/4nyBJfL\nxeTJU+nQ4Vbuuut+9u/fb3QkOY/z/Tuiu52+y1t/croDJuKnUlNTqVOnMceP34XD0QiL5U36929J\nkyYNePTRMVitY4FjhIe/yKpVS6lfv77Rkf3eyJGjmDLle6zWkQQEbKNkyffZsWMjFSpUKNB5li9f\nzq5du6hbt65PLHtktVr59ttvsdlsdOjQgUqVKhmSY+PGjSQmJhITE8MNN9xwwVI1atQY3nxzHlbr\nMwQE7KREiWls27ae6OhoDyeWc+kRpIj4hEOHDvHMM2M5dOgIN998Pf/736MEBATw+edf8MEHc4iI\nCOO55/5H48aNjY5aLFgsUdhsm4DcsURhYXfz+uvXMnjw4Hyf48knRzF16qdAO1yupYwcOZjnn3/K\nPYGLQFpaGs2aXUdychmgHCbTzyxbtoiGDRt6NMe0ae8xfPgLBAR0AtZw223tmT797fOWMJfLxbRp\n7zF37g9UqFCaF198lho1ang0r5yfCpiIiBSYxVIKm20bkHsHKDS0P6+/3oxHHnkkX8fv3buXevVa\nYrPtAMoAhwkJqcO+fdu54oor3Ja7MOLjx/LKK7+Tnf0xub8KP6BFizmsXv2TxzJYrVaioiqQnf0b\nUBNIJzy8AUuXzqZ58+YeyyGFpzFgIiJSYIMG3Y/F0gdYSEDAa4SE/EjPnj3zffyRI0cICqpGbvkC\nqEhISCWOHj3qhrRF48CBZLKzm/HP78ymJCcf9miGEydOYDaHk1u+ACIwm+tw+LBnc4gxVMBERIq5\nt956leee60nLlpPp2nUTa9f+UqBxRXXq1AEOAAvIXd73C8zmE9SqVctNiQuvY8e2WCzvAslAJqGh\nr9K+fVuPZqhYsSJRUZGYTO+S+/ftVxyOtTRqpIl0iwM9ghQRkUJbuXIlPXr0JTU1iQoVYliwYA5N\nmjQxOtYFuVwuRo0aw7hxr+J0OujcuTuffz6jQJP8FoUdO3Zw8/+3d+fxUZX3Hsc/k5B1MiQQkLAp\nQiIoVcFLFGRV1oi4tCKL5cJF1LgUFKlQW5XFai3WImrFSIUEroKEK4tgBCqRCCqLQUACkghYpIAK\naZbJZJ37xxlkLAmQkJmTk3zfr1deM+fMYfKd10NmfvOc5zzP0Ls5fHgfDkc0S5cmM3jwYL9mkIun\nMWAiImIql8tV48lyzVBeXk55eTnBwcHnP9iHXC4XISEhmlbCojQGTKSeOHToEL17J9C8eTt6907Q\nvExiGVYqvgACAgL8tjLC3Lmv0abNlbRu3YkXXvgL3h0NoaGhKr4aGBVgInWMy+WiV69BbNnSlx9+\n+IgtW/rSq9dgXC6X2dFE6pWZM58nNNRBSEg4d9wxGqfT6bPflZKymN/9bi7ffZfC0aNvM3PmQubN\nqz9rkkr1qQATqWOysrLIywumomIa0J6Kimnk5QXVq4V3N2zYQLt2v8DhaM4ttwzn1KlTZkeSeqCi\nooKpU5+iSZPWREdfyp/+9CJVDWdZtmwZL7ywiJKSLMrLT/LhhyVMnDjVZ9kWL16B0zkdiAeuw+l8\nlsWLV1R5/DvvLCEmJpbGjVswZsz9FBUV1UqO/fv306VLbyIimtGlS2/2799fK88r1acCTKSOiYiI\noKzsJHD627iTsrKTREREmBmr1hw4cIDbbx/F4cOzKSjYzT/+Ec2dd/7a7FhSD8yePYdXX11Hbm46\nJ09+wKxZb5GcvKjSYz/4YCNOZyLQFrDjcv2B9evTfZYtMjICm+07rz1HiIys/G86IyODe++dzPHj\nKeTn7yA19QQPPfT4RWdwOp306TOEXbtGUFi4l127RtCnz5BaK+6kelSAidQxsbGxDB06ELt9IPA8\ndvtAbr11YJ2+pL86Nm7cCNwKJAAxlJTMJSNjHeXl5ef5lyLn9u67a3A6ZwJxQGeczidZtmxtpce2\nadOC4OBMrz2ZxMRUb+ml6pg+/bfY7X8mIGAKNts07PbpPPts5SsFrFmTRlHR/cCNQBtcrhdZvbry\n135XccwAAA5hSURBVFEdWVlZFBU5cLsfAS7B7X4El8tRr3rXraSR2QFE5OdsNhtLly4gOTmZ3buz\nuPrqCYwdO7beDNCNiooiIOAbjHmPbMA3hIZGWHbh7bS0NPbs2UOnTp0YOnRovWknK4qOjgJyftoO\nCMihWbOoSo+dPHkSixf34sSJBNzuZgQGruP11z/wWbbOnTuTmbmF5ORFuN1u7rlnk2f+tLNFR0cR\nErKL4uLTe3KIjKz8dVRHVFQUZWXHgXzAAeRTWnqcqKiLf26pvrr6TqFpKKTB+Oyzz1i0aCmhocEk\nJk4gLi7O7Eg+VVxcTPfu/fn6awcuV1dCQxfx0kvP8MADE8yOVm2TJ/+OpKT3KCkZQnDwBkaPvpmk\npLlmx2qwdu7cSa9eAykuHonNVozd/j47dnxS5XqJBQUFrFq1iqKiIgYNGkTbtm39nLhyp06d4tpr\ne/D9910oLb2U4OCFpKYu5JZbbrno5x437kFSUz/H6UwgPPwDhg/vzoIFf6uF1A2b5gETsZj169dz\n++33UFT0GDZbPhER89m2LYOOHTuaHc2nXC4XKSkpHD9+nD59+tC3b1+zI1XbkSNHiIu7FpfrANAU\nyCcsLI4vv8yo90V0XZaTk0NqaiqBgYGMGjWK1q1bmx2pRnJzc0lJSaGgoIAhQ4bU2sL0breb1NRU\nsrKyuPLKK7nrrrvUa1sLVICJWMz11w9g27b7gbsBsNmmM2HCjyQlvWJuMDmvXbt20avXKPLzv/pp\nX2RkPGlpr9C9e/cLfp7S0lJWrFjBjz/+SO/evencuXOVx7rdbtLS0jh48CBdu3alR48eF/UaRKR2\n1KQA0xgwERMVFjqB5j9tu92XUFBwxLxAcsHi4uIIDS2koCAJt3s08B6Bgf86ZwH1n0pLS+nTJ4E9\ne1yUl3fCZnuaJUv+zrBhw8461u12M27cgyxf/gkVFT2x2Z5nxozHmTLl0Vp8VSLiL9Yc9SpST4wf\nP4Lw8EeBT4F1hIc/x9ixw82OJRcgLCyM9PS1XHFFEo0aNadDh7/w0UdrcDgcF/wcS5cuZffuMgoK\nNlFUNB+ncznjxz9S6bFffPEFqalpFBZ+TlHRGzidW/j9758iLy+vtl6SiPiResBETDR58kTKysp4\n442HCAkJZvr0l7QQr4VcddVV7Nu3vcb//sSJE5SWXsOZ78JdyM09UeWxQUFxwOnFotvSqFEkJ0+e\npHHjxjXOICLm0BgwERGTbNu2jX79bsPpTAOuJChoKt2772fTprPnfDp27BhxcddQUJACDMBmSyIm\n5iW+/XYfjRr57rt0aWkp+fn5NGnS5GeDtfPy8ggKCiIsLMxnv1vEKrQYt4iIhcTHxzNv3otERPQn\nIMBOt25fsXx5cqXHxsTE8P77y2jePBGbLYTY2DfZuHGNT4uvBQuScTia0rLl5bRvfzXZ2dnk5+dz\n883DiI5uSePGTXn44clVLvcjIlVTD5iISB1QXl5OYGBgrR9bU19++SU9egymqCgd6ITNNofY2GRu\nvLE7S5YUUFz8FlBAePgg5sxJ5L777vVpHpG6TD1gIiIWVZ2CytfFF8D27dsJCBgCdALA7Z5ITs4e\nPv54C8XF1wCjgYk4nQNIT//M53lE6hsVYCIichZjVvjtwOmFmrditzchNDQQmAPcBsQDr2G3+74g\nFKlvVICJiMhZBg4cyNCh3bDbu+Jw/Irw8GG8/fZblJUFAAuAMcBEYCqlpdUbMlJcXMwTTzxFfPwA\nRoz4H44c0dx30vBoGgoRETmLzWZjyZIFbNq0iWPHjtGt25/p0KEDTzzxR8C7xyuIwMCgaj33yJHj\n+fDDf1NUNIXMzM18/HEf9u/PJDIyslZfg0hdpkH4IiJywebPf4tJk57D6ZwNnCQ8fBrp6WuJj4+v\n9Pjs7Gzmz19ASUkpY8aMomPHjkRGRlNWdhIwprBwOAaTnJzInXfe6b8XIlKLtBSRiIj41IQJ4wkO\nDiYpaT5hYSE8/fT/VVl87du3j+uv70th4TgqKuy88cYgVq58x/NoudeRZVoQWhocs/7HtwVSgEuA\n74Ek4G2vx9UDJiJicWPHJrJoURvc7j949qTQs+cSLrusFStWHMTpTCQoaDMxMWns3budiIgIU/OK\n1JSVpqEoBR4DOgN3Ac8CF76AmoiI1Hn5+U7c7hivPS0pLHSSnDyPJ58cTP/+7zBuXAU7dmSo+JIG\np670+a4GXgI2erbVAyYiYnErV65k9OhJOJ3JQDh2+33MmnUvjz32G7OjidSqmvSA1YUCLBZYB1wN\nFHr2qQATEakHFi5MYdasv1JWVkZi4n8zbdoUjfeSeseKBZgDSAdmAiu99rufeeaZnzb69etHv379\n/BpMRBqGgoICQkNDfbqmoojUL+np6aSnp/+0PWPGDLBQARYErAHWYkyr7E09YCLiUydOnCAhYTi7\ndm3FZgvg+eef4/HHJ5kdS0QsyEo9YDYgGfgBmFzJ4yrARMSnbrppGJ98cgVlZbOBfxIe3o9Vq+bT\nv39/s6OJiMVY6SrInsCvgZuBTM/PEJOyiEgDtHXrFsrKpmK8DV6GyzWSTz/91OxYItJAmDXo4RO0\nDqWImKhFi9YcPPgpcDtQTmjoVlq1usfsWCLSQJg9CL8qOgUpIj6VkZFBQsIvCQi4Cbf7INde25SN\nG98nKKh66xqKiFhpDNj5qAATEZ87fPgwmzdvJioqikGDBulKSBGpERVgIiIiIn5mpUH4IiIiIg2W\nCjARERERP1MBJiIiIuJnKsBERERE/EwFmIiIiIifqQATERER8TMVYCIiIiJ+pgJMRERExM9UgImI\niIj4mQowERERET9TASYiIiLiZyrAREQuQl5eHrt37yY3N9fsKCJiISrARERqaNWq1bRseTk9e46g\nVav2vPtuqtmRRMQiqrVytx+53W632RlERKqUm5tL69YdcDrXAjcAOwkL68/Bg3tp0aKF2fFExI9s\nNhtUs6ZSD5iISA0cOnSIRo1aYRRfAF0IDo4lJyfHzFgiYhEqwEREaqBt27aUln4HfOXZc4CSkmza\ntWtnYioRsQoVYCIiNRAdHc2bb75GWFhfIiN7EhbWg5dfnk2rVq3MjiYiFqAxYCIiF+Ho0aNkZ2fT\nvn172rRpY3YcETFBTcaAqQATERERuQgahC8iIiJiASrARERERPxMBZiIiIiIn6kAExEREfEzFWAi\nIiIifqYCTERERMTPVICJiIiI+JkKMBERERE/UwEmIiIi4mcqwERERET8TAWYiIiIiJ+pABMRERHx\nMxVgIiIiIn5mVgHWB8gCDgC/MSmD+Eh6errZEaSG1HbWpvazNrVfw2JWAfYy8AAwAHgYaGZSDvEB\nvYlYl9rO2tR+1qb2a1jMKMAiPbebgMPAOuAGE3KIiIiImMKMAiwe2Oe1vRfobkIOEREREVPYTPid\nA4B7gVGe7USgNfCU1zHZQAc/5xIRERGpiRwg1uwQ5xMJZHptvwIMNSmLiIiISIORiXElZDuM05Ea\nhC8iIiLiY30xpqHIBiaanEVERERERERERMR/ZmP0jH0BzAHCvB6biDFx616gl/+jyQXQBLvW0hbY\nCHwFpAOjPfsdwErgW2AFEGFGOLkggRhDOlZ7ttV21mEHkoGvMT7XbkDtZyX3AVuAHRj1Cli8/QZi\nTI0RALyJcbUkwCUYY8UuxTh9+YUp6eR8To/tuwyN7bOCGKCL534z4BuMN5AnMC6OCQFeBaaYkk4u\nxGTgf4FVnm21nXW8CMwCQoFGGBeoqf2soSlwEKOIDgDWAoOpR+13F5DiuT+MMxUmGB/0Dr8nknP5\nz6tb56KrW61mNXAzkMqZwuw6YJlpieRc2gAbgJs40wOmtrOOnfz8LA+o/awiDDgEtMIowtIxejDr\nTft9CAz33J+FsXTRaUuA/n5PJOcyAHjHazsRo93EGmIxesAiMFaoCPXsD/dsS92zDOiKcVbgdAGm\ntrOGNhhnCRYCnwNTMT7U1X7WkQCUAPnAHz37qtV+ZsyEvx7YXcnPMK9jnsZ4Uaerx8omjHX7MKNI\nQ+IAlgKPAQWYM0GzVM+twAmMXmfv9lLbWUMocAWwHOgHdAbuRu1nFc2B14GrMKbT6oHxN2n59hsH\nbOZMFQlGcfay1/ZOdAqyrtEEu9YUhLEe66Ne+5Zj9KwA/BdGt7rULc8B/8QYh/IvoBBYhNrOSrK8\n7idgnEFQ+1nDUIwzcac9CLxANdvPjB6wcxkC/Ba4DXB57d+KMcDtUoxvCxUYPWRSd/zbc3t6gt2B\nGF3rUnfZgL8De/j5GMvPgfEYp0TGA5/5P5qcx5MYV7FeDowEPgLGoLazkgMY44YCMD7QN6D2s4oM\noBvGYPwQjAJ6HRZvvwMY50wzPT9/83psEsbErXuB3v6PJhdAE+xaSy+MLzM7OfM3NwSLX0rdAPXl\nzFWQajvruALjA3onxhWRdtR+VjIO+BjYhjHeOQC1n4iIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiI\niIiIiIiIiIjIHRjzm3U0O4iISFUCzQ4gIlLLZmAs0dMMSDc3ioiIiEj9FwEcwli2zHutvYnALiAN\neA/4lWd/R4xFdT8HXgOi/RVUREREpL64B5jnub8JuA6I99x3AJ2AIuCXnmNWYaypCPAQMNVvSUWk\nQWtkdgARkVo0Cvir5/4yz3YuRq9XPrCPMwvkXoKxruzpdRQDMXrPRER8TgWYiNQXTYGbgF8AboyC\nyg28WcXxAcCPQFe/pBMRERGph+7HGM/lLR3o47l1YIz5cnLmFOR6jPFgNiAIuMoPOUVECDA7gIhI\nLRmJcarR23LP/hXAFuBljAH3Bz2PP4TRa7YTyAR6+CWpiIiISANg99y2A742MYeICKAxYCLSMKzB\nGCN2FJhkchYREREREREREREREREREREREREREREREREREREREZE66v8BSqjb4+UfAvMAAAAASUVO\nRK5CYII=\n", "text": [ "" ] } ], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
Removing Negative Age Values
\n", "\n", "* The plot allows us to spot the negative age value. Let's remove it.\n", "* `pandas` allows for removal by indexing according to 'truth values'.\n", "* These can also be used to make different colour plots." ] }, { "cell_type": "code", "collapsed": false, "input": [ "df = df[df['Age']>0]\n", "fig, ax = plt.subplots(figsize=(10, 7))\n", "for symbol, sex in zip(['bo', 'ro'], ['Male', 'Female']):\n", " lines= ax.semilogy(df.Age[df.Sex==sex], df.OI[df.Sex==sex], symbol,linewidth=3,markersize=10)\n", "ax.set_xlabel('Age').set_fontsize(18)\n", "ax.set_ylabel('OI').set_fontsize(18)\n", "ax.set_title('Plot of age against OI on a log scale').set_fontsize(20)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAnIAAAHNCAYAAACaWpaoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2YXGV98PHvEoYkC2rFBIXFOL5UsKlalZVs0WWNJkal\n1PZZUvEpiiattLW2Qa0vEQhxaUu1purTCi1gDFpsXKviS2IgcbOiu7goFEUBRZZoQEgKKmYgTJJ5\n/rjPsJPdeZ+dOXN2v5/r2ivZM2fOuefeM3N+c7/8bpAkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIk\nqe2kgUPAJ2MuR61eCnwEuJNQ/i/EW5xZ4xDwjbgLoaYaB+6OuxBtqI9w/V8UczkUgyPiLoBmnUOT\nfvYTPpi3AK8p8Zxcg+c8NzrXmxs8TjU6gP8G3gTcAKwDrmnBeRU0eq3U4lwav65+G/hXYBj4X+Bn\nhPfCB4EnlnhOH7M7aG3l3zhprJtZ6Mi4C6BZKQdcHP3/GOB1wKujnwuAS5p43mZ7GXAC8DfAx1tw\nPk04GcjEcN56r6vzgQEgBYwCnwGOBroJ74U3A6uBbdN8XkmS6nYIOFhk+7nRY1ngqdG2dLTtqgbP\nmT92K1rkPhCd6/QWnEvxOpf6r6s/j577M0LwP9lfAI8QWqyfN+mxvui5O+o4b9KNAz+NuxBtqI9w\nTVwYczkkzQKlArkUsCd6bEW0LU3pQO4Y4MPATcCvgR8RWjSeOWm/IaZ25+Z/FlVZ5lcCXwN2Ebq/\ndgJ/NWmfvjLn6a1w/CcC7ybcmH9OaFW6E/gv4Lllnvc6YDvwEHAroYvuSUy85mJOAT4fHf/XwBih\n+zdVoYyFTiDcML4F/ALYB/wQuBI4vszzziV0N/8yOu8/AEdRfNxTPXVSrLtxHROB9auBawnX2Tjh\nuvqtIsc5EXg7YWzjg4S/+a3AJ4Bjo32GqP+6Oi467kHCeMpS3hsdb3jS9j7qC+ReRHhNdxH+Bt8m\nfPEo1jMzxMQ19A7CNf9r4DZgLbUNy3ku8I+E9+oe4GHgf4B/IVyvtRineCDXQXhP7mSii/qrhPdu\nKedS/fVYzosIPQlDwG+A+4ARigdVRwB/SXjf3kuoix8CHyVcF3m11lkfpQO5udH2GwifFfcQ3gen\nVPsCJalQqUBuLrA3enx5tC1N8UCui3AzOkT4cPowITh5jHCzeVXBvm8m3LwOEcauXVjwU81N5H3R\nc/cCnwI+BtwRbfsq4QYC8AzCQONvFJQ5f55nVDjHEkLLyzZCsPARwgf4AeBRoKfIc/4yOs9DwCbC\njeD7hBvIdylex++MjnkQuBH4N8LN5FD0nLkVypn3BkLw9mVC9/HHCDf4Q4SbYrrIcy6NHr8XuCI6\n913AFwmB2uSbcz11Uiy4WRdt/wTh+rgTuAx4INp+86T9FxKCgCwhUP0YsAH4EuFm+jvRfo1cV/kA\n7YsV9psH3B/tWxjw9VF7IPcmwuv/DWHM5keA70XHGQM6J+0/FD32cSaCyU8TWgkPEa63ar2XcJ1+\nnhCIXM7E+3cXpccCFjPO1GvlCMIXrUOE9+bHgY1MfDF8b5Hj1Ho9ltJLqNdfAVuBf4qONUQI6Aod\nTQieDxGuwyui/T8fPb/wC1+tddZH8UDuBOAH0WP/S6iXbYT30WPAG6t8nZL0uFKB3Fujxx5j4ptp\nmuKB3Kej7f8wafurou0/IXyzzjs32v6mGsv6PMIN/V4OD8bmAddFx1w96TnrqK4VrtATmWjpKbSY\n8IE7Nmn7cYSg4iHgOQXbjyQEHMXq+PmE1/JD4OkF249gIli9oMryLiTclCZbER3nc5O2/25Unrs4\nvNXhGEJAeYipN85a6wTKB3JZDm+dOYowqeAQcGbB9rdE2z5S5PjzCX/7vHOp77raFD3v/Cr2/Wq0\n77kF2/qoLZBbQLhW9gG/V7D9CMKXk0OEsXqFhqLtPyKMPcx7CiEIfjQ6bjVOoHiL73nROT5U5XGg\neCC3OjrOdg7/+ywCdhM+UwpfQz3XYyn5+juzyGOTr98N0b6bmFofnRwenNVaZ30UD+Q2R9v/FphT\nsP1EQqvjHuDJRc4jSSXlu58uItxkP8xEa85BQlCRl2ZqILcw2naAwz+E874fPX5WwbZzqe+G+/+i\n5/1rkcf+T/TYbZO2r6P2QK6cz0XHO6Zg299RPGACOIPigdxno+3nFHnOMYTAsNqbVzljhG/+hf6N\n0jfst1PbjROK1wmUD+S2FjlO/qZYGMD2R9veU0U5zqW+62oset4fVrHvx6J9Ly3Y1kdtgdy7mGhB\nnuwl0WO/nrR9KNperDUrfy29vMrzl3IUIdC6qYbnjDP1Wsm3OJ01Ze+JFsXC9/B0Xo/5z4hTK+z3\nRCZawyePeaxFqTrrY2ogdzITrZQdTLWO+sd4qo04a1Vxyec7OkDoyvg64aa1pcLz8uOjbia0DEz2\nNUKrzUnTUMb8Mb5W5LFthA/lZxNaNg41eK7fA9YQutC6ODxIyRHSVOS7AfPluq7Icb5BqNM5k7Y/\nPyrvcwkf4JPdTwicjyK0YFRyOmHs1AsIrQfzJ5X3twjdrJXKW2pGJtRWJ5VsLrItP/bsxIJtXyJ0\now8Ay4BBQp3eUeV5qlHsptpM5a7j7xLeRwsJLbU/m/R4qXpbyeH1VsmZhMD5ZMI4ysJu/GKtu9U6\ngtAqfYDwGTLZ1whj5wo/D+q9Hou5kjDM4TrgK4TrZydh7Gih50Vl/VH0U41G6+z50b/3UTy/XP7v\nV24crhLAQE5xyDE10KhW/sPnxhKP3zhpv0acSChrsXM9TPhAXkxoGZz8wV2LpYSbUJbQyrKNEAQd\nBF5BCJqeULD/CdG/3y1yrH3A7VG5CqUJdb62TDlyhBnDk2/mk72JMNbm18D1hBvYrwjB7B8BLyS0\nQOQDuROiYxcr752Eupys1jqp5CdFtuVTlRR2YWWj8v8lYWbp0mj7OKEr/z9qOGcpPwJezOHd4qU8\nO/r3hw2cr9J75juEltwuDv/b5yjeMlWs3sq5kPDlYQ+h+/M+wrWTI3RlL6ryOMUsJHz5+CFTWxWh\n+OdBPddjKTcT/pZ/Qxg7+obo2N8htJx/M9ovPzRjZ5XHnY46S0f/9lK6hyAHPK3KMqlNGcgpafI3\nmlJdGUuif38+Def6OeHb+6lM7ZZ6IhNj6O5v8DyXEL6tv5GpA+B/v8j+90b/nsLUm9ExHD4eKO/n\nTG3Vqld+oPgrmNoi9kam5je7l1CPpzC11eQkQkA2uTu21jqZTvcRulvXEWYkLgP+mjDg/FZKB0TV\nynfHvxz45zL7zSd0feZoLJDLvxdOZWqXXAehxTPH9LxnJuskBCUZQgtr4QSAo6hunGA5ewgtyL9N\neE9ODuaKfR7Ucz2WcwshuPorwrX5B4QvAV8kvOceJXwRgOrSEk1XneVf8/s4vGteM4wrOyhp7oz+\nfRHFB1vnV4e4vWBbNvp3HrXJH2NFkceWEVq4fkLjiVlPJqQEmBywLCQM0J98/Hy5lhU51iso/gXt\nFsIN4iX1FxMIMzKfSuhemxzEvYji3TTlyru8yDaovU6a4SAh8PkHJtLN/HnB4/VeV1cRWjDPpHz6\nkXcQWntvoPjkjmqVu45fHJ1jH80J5H6bcJ/ZzNRZnK+jthmrxRwCfkxoHSx2fRX7PKjneqxGhtBC\n/TfAvxMmEeTH7d0elfVkJmY+lzJddXZL9G+xPIWaQQzklDR7gf8kXLvvmPTYKwldij8ljFXJK2zB\nqsW/Esbe/BGHz/ScB7wt+v+/1HjMYr5K6AYp7A5dQFj9olgX9KcIN95XET7081LAn5U4x98TAo//\nINy8JzuaygO2IQQg3472LQykTyTkJCvmE4Sb2B8zkewZQstHsckXUHudTJduDi9jXr4btHASSb3X\n1QNMTOoZBE4rss+fE1oEHyOMk2rERsLf7RVMjJuC8B766+j/H2vwHKV8n9CKvpzDu2JPovHWuLx8\nKpQ/5/BxZE8nvHezhEkJefVcj6W8kuKB/OTr5deEOu4g/O2PmrT/MUwEaNNVZz8iXF+vi55XLC3O\n8wlfjiSpaqXSjxSTpnQeuZ9Gj32T0D31ecIH9q+Z+k37CYQA8GHCjMcLCUFHNd9s10bl3UO4IX6c\niTxyxSZmrKP2WasvY2Lm4H8TAtVHCR/m+RmCk4/311G5CvPI/YAQZH2X4hMW/oqQm+0QobXr36Kf\nrxLqpthg+GL+lIm8VP9JSC6aJdyAtlI8Ke6Hou33EQaI5/N2fYHQEnTnpP3rqZNys1aL/T3STL2+\n/oVQd9cR8s1tILQ+PkoIhl5YsG8j1xWEhMeZ6HzDhGvrSkL3bf5vVKzVqI/a88i9mYk8cv9JeF35\nPHLfpXQeuWLOpbbZuvnVTu4l1HU+1+I3CC2etUwUGmfquL0OJq67OwhB2ycJf5sDwPuLHKfW67GU\nWwjvwS8wke/t5ujY3+fwcZyFeeR+HJ33nwjXzuQ8crXWWR/F048sKCjPfkJX8j8RkqfnZ/iXaxWW\npCmmI5CD8A32nwk3oYcJXRefAZ5V4lgvJ3wI5hPBHqT6QdavIgQ5P2NiZYe3l9j3oujYtaYfOZ3Q\nUnAf4QP2qqh85Y53BuFmnl/Z4f8RunMeINyUijmJEPjdQrip30+4oX+Y4i11xRwBvDY6zl7CjSV/\n7k9Sum7fQkiyW5hJ/6mEIHB7kf1rrZNiwU25+ksz9fp6KeGmfgvhb72H8CXh7UxdNQQau64gdEX/\nK+ELSX5Fgi3AekpP5jid2gM5CH/f/MoOvyIkj76A4l3x36D0+/TN0WPVBnJzCV2MnyesaPFtwqoF\n8yqcp5i7KZ0a5O2EgDhfj1/j8OTgk9V6PRZzFiEwzk+S+CkhQPu/FF81pIPQk5C/Zn5N+AL2EQ5v\nGau1zvooHsjBRMvr9YTUJRnCNfBlQg6+yUG81DJ/SBiHcBV+o5DylhI+0K+PuyBVyieCviLugkh4\nPUqxOI7wrV2aTRYwNf3D0wjf4A8RWhvayVOZOiY3TWgRyc+AlVrF61FqkqsIXT3fn7S9lzBw88dM\nDM7N+2cOX3ZGmg3OI3TNfIYw3meI0N1yiDDLsdVJZyv5R8KYr42E8uZXNzhEWPtTaiWvR6lJXk5I\nYTA5kLuZEMw9gzAW6imEG9U/cfj6idJs8XuE1refEca87CEEcO+iPWejLyWMWbqXMHHgXsJYn7fG\nWSjNWl6PmjHa7Vs7hObtLzMxTf5JhNaGF0W/f4ww8+ZZhEG3Y4SByZdPPtCzn/3s3F13lRrzLUmS\n1FbuorpVXx7Xjt/cJ+vm8GSOPyRk6/44IX/TX1AkiAO46667yOVy/tT5c9FFF8VehiT/WH/WnfWX\nzB/rz/qL64eJZfmqloRATpIkSUUkIZAb4/C1IxcDozGVRZIkqW0kIZD7VfRvL2H83DIaX7RaVejr\n64u7CIlm/dXPumuM9dcY668x1l9rtdtkh2sIWcufQkitcCEhU/zphOVyUoTJDtWuC5iL+pwlSZLa\nWkdHB9QYmxVbmiVOZ5fYvhN4Xj0HXLduHX19fX5DkCRJbWloaIihoaG6nttuLXLTzRY5SZKUCPW0\nyCVhjJwkSZKKMJCTJElKKAM5SZKkhDKQkyRJSigDOUmSpISaE3cBmmxd/j/pdDq+UkiSJJUwNDTE\nxo0b2blzJ8DFtTzX9COSJEltwPQjkiRJs4iBnCRJUkIZyEmSJCWUgZwkSVJCGchJkiQllOlHJEmS\nYmT6kdJMPyJJkhLB9COSJEmziIGcJElSQhnISZIkJZSBnCRJUkIZyEmSJCWUgZwkSVJCGchJkiQl\nlAmBJUmSYmRC4NJMCCxJkhLBhMCSJEmziIGcJElSQhnISZIkJZSBnCRJUkIZyEmSJCWUgZwkSVJC\nGchJkiQllIGcJElSQrmygyRJUoxc2aE0V3aQJEmJ4MoOkiRJs4iBnCRJUkIZyEmSJCWUgZwkSVJC\nGchJkiQllIGcJElSQhnISZIkJZSBnCRJUkIZyEmSJCWUgZwkSVJCGchJkiQl1Jy4C9Bk6/L/SafT\n8ZVCkiSphKGhITZu3MjOnTsBLq7luTUtzJpAuVwuF3cZJEmSKuro6IAaYzO7ViVJkhLKQE6SJCmh\nDOQkSZISykBOkiQpoQzkJEmSEspATpIkKaEM5CRJkhLKQE6SJCmhDOQkSZISykBOkiQpoQzkJEmS\nEspATpIkKaEM5CRJkhLKQE6SJCmh5sRdgCZbl/9POp2OrxSSJEklDA0NsXHjRnbu3AlwcS3P7WhO\nkdpGLpfLxV0GSZKkijo6OqDG2MyuVUmSpIQykJMkSUooAzlJkqSEMpCTJElKKAM5SZKkhDKQkyRJ\nSigDOUmSpIQykJMkSUooAzlJkqSEMpCTJElKKAM5SZKkhDKQkyRJSigDOUmSpIQykJMkSUooAzlJ\nkqSEMpCTJElKKAM5SZKkhDKQkyRJSigDOUmSpIQykJMkSUooAzlJkqSEmhN3AZpsXf4/6XQ6vlJI\nkiSVMDQ0xMaNG9m5cyfAxbU8t6M5RWobuVwuF3cZJEmSKuro6IAaYzO7ViVJkhLKQE6SJCmhDOQk\nSZISykBOkiQpoY6MuwCSJGUyGYaHRxkcHGHXrkdYtGg+/f099PYuobOzM+7iSW3LWauSpFitXHkh\nY2Nz2b27h2x2CdAJZEilRunqGqG7ez+bN6+Pu5htzUB4Zqhn1qqBnCQpNplMhsWLNzA+vrbkPun0\nALfddr4BSQkGwjOH6UckSYkyPDzK7t09ZffZvbuH4eHRFpUoWTKZDGNjcxkfX0s2u5QQxAF0ks0u\nZXx8LWNjR5HJZOIspprIQE6SFJvBwZGoFam0bLaHwcGRFpUoWQyEZSAnSYrNrl2PMNGKVEpntJ8m\nMxCWgZwkKTaLFs0HKnX7ZaL9NJmBsEw/IkmKTX9/D5s2jUbju4pLpUbo7y/ffThbTQTC5YK5ZAfC\nmUyG0eFhRgYHeWTXLuYvWkRPfz9LenudAIOzViVJMXLWamO2bt3BmWdSIRDezrXXdrBiRel92tWF\nK1cyd2yMnt27WZLNRvNxYTSVYqSri/3d3azfvDnuYk4bZ61KkhKls7OT7u79pNMDpFLbmehmzZBK\nbSedHqC7+zGDuBJ6e5fQ1VV+/FtX1wi9veXH0bWjTCbD3LEx1o6PszQK4iC0PS7NZlk7Ps5RY2Oz\nfkauXauSpFht3ryevXv3ctVlV7Dj6r+jY+9D5BY8maXnnMVbzzuPBQsWxF3EtpUPhGEgyiPXw0Qe\nuZEoj1wyA+HR4WF6du8uu0/P7t2MDg+zdMWKFpWq/di1KkmK1YUrV5L6zndY8vOfc9rBg493n31r\nzhxGTzyR7EtfOqO6z5phJq7scMnq1ay58soKo/9gw6pVrL3iilYVq6nq6Vq1RU6SFJtMJsMvv3od\nH8v88rDtncCygwdZds89vGPPL8lkMokNSFqhs7OTFSuWJnIcXCmP7NpVxXzcsN9s5hg5SVJshrZt\n47WZX5fd5zWZhxnatq1FJVK7mL9oURWJacJ+s5mBnCQpNl++7N/p5VDZfU7nEF++7N9bVCK1i57+\nfkZTqbL7jKRS9PT3t6hE7cmuVUlSbH7x/Tur6j77xffvbEVx6maus+m3pLeXDV1dLB0fL7nPSFcX\n5/f2tq5QbchATpIUmz0cU0U627BfuyrMdbamMNfZpk1smIG5zlqls7OT/d3dDBBmp/YU1O1IlEfu\nse7uWR8oG8hJkmLz5N99IcP3/pAVZEvus5MUT/7dF7awVNUrzHVWKJ/rbOn4OAPRfrM94KjH+s2b\nH2/t3DCptfN8WzsB049IkmL0xS9+ja/80Squ4Bcl91nFU/mDL1zF61//2haWrDo7tm6FM89kabZ0\nILo9laLj2mtnda4zVceVHSRJibJ8eR/f61zAKtJsIVWwrgNsIcUq0tzcuYDly/tiLGVpI4ODLCkT\nxAH0ZLOMDA62qESabexalSTFprOzk+e87o+4/kb4z5/P52mHbmch97GH4/nFESdz3ImPcOqptG0X\nWlJync3EhMEKDOQkSbHavHl9QaBxiF27ns4LFs2nv//Utg808rnOKk3WiDPX2cqVFzI2NjdawmsN\n+SW8Nm0apatrA93d+9m8eX1s5VNjHCMnSVKd2n2MXCaTYfHiDYyPry25Tzo9wG23nd/WAfNs4Rg5\nSZJaaElvLyNdXWX3GenqYklMuc6Gh0fZvbun7D67d/cwPDzaohJpuhnISZJUp8dznaXTbE8dPllj\neyrFQDoda66zwcERstklZffJZnsYHBxpUYk03RwjJ0lSA6Yj11mzVobYtesRyo/gA+iM9lMSGchJ\nktSgzs5Olq5YUdc4uGauDLFo0fzoaOWnY4T9lER2rUqSFJPClSGWRkEcTKwMsXZ8nKPGxshkMuUO\nU1J/fw+pVPnxb6nUCP395cfRqX0ZyEmSFJPR4WF6du8uu0/P7t2MDg/Xdfze3iV0dZUf/9bVNUJv\nb/lxdGpfSe5afSawFngScFbMZZGkRGt0jFazxni1qvxxGRkcZE0VK0NsGBysq9u2s7OT7u79wECU\nR66HfB65VGqErq4Rursfa+s6UnkzIY/c5ygdyJlHTlJLNJo5P85ApHCM1pLCMVqpFCNVjNFq9Plx\nlz9OH1i+nIHrrqu837JlDGzbVvd59u7dy2WXbeLqq0fYu3cOCxYc5JxzejjvvDexYMGCuo+r6VVP\nHrl2a5G7Cngd8ADw/ILtvcDlhPJ+DPh464smScU1mjm/mYPdKykco1UoP0Zr6fg4A9F+xQLKRp8f\nd/nj1oqVIQ6/Ps8DOnnwwQzr149y5ZWXV7Wyg0t8ta92C+Q+SQjSNk3a/lHgbcA9wNeBa4C9rS2a\nJE2VyWQYG5tbJHN+J9nsUsbHlwID7N27l1tvumlKi9sLTjkl1kCkljFaxbr2Gn1+o+I+f6N6+vsZ\n3bSJJdkso8AI8AgwH+gBlgAjqRQ9/f11Hb/a67Pc9eUSX+2t3QK5bwLpSdueFP2bH+m5DTiVcL3/\nPfB7wHuAS4sdcN26dY//v6+vj76+vukqqyRVlTl/3j3Xs+F3/5VXPvi/U1rc/uXYY9nzwANln9/M\nQKTRMVrNHuNVSdznb9SS3l7elEoxks3SA0yESTAKbABuTqXYVOfKELWs7LBixdIpj01HIKjShoaG\nGBoaaugY7RbIFdMN3F7w+w8JX1K+CpxX6cmFgZwkTbeQOX9NmT0ynJYb55L7f3HY1sdb3O6/n/WU\nz/TVzEDkkV27qkgXG/ZrxvMbFff5p0OaMHOvUCewNPp5ZwPHrnx95ld22FA0kGs0EFR5kxuYLr74\n4pqPYfoRSWpApcz5KYbp596yxziN0PpSSjMDkfwYrXLKjdFq9PnVymQybN26g9WrL2H58g+wevUl\nbN26g9Txx7fk/M0yOjzMayu0KL42m607/UijKzu4xFf7S0IgNwacXPD7Ysp/5klSy0xkzi+ui0F6\nqdD1RxgrUkozA5Ge/n5GU6my+5Qbo/WiM87gWxXOcUO0X71WrryQxYs3cOaZcOWVa7juugGuvHIN\nZ54JV2yFb3aUv5U1Msas2UYGB1lSRdfwyOBgXcevdH0GpVd2cImv9peEQO5X0b+9hBboZcCNsZVG\nkgpUypy/kCq7/so83sxAZElvLyNdXWX3GenqYkmZMVo7K5yj0uPlFI7RymaXMhFUhDFaP3/gE3zh\niCeWPUal8sep2V3Dja7s0GggqOZrt0DuGuDbwHOBnwFvibb/LSH9yPXAv+GMVSlxMpkMO7Zu5ZLV\nq/nA8uVcsno1O7ZurXvpoXZRKXP+HqrseizzeDMDkc7OTvZ3dzOQTrM9lXq8rBlgeyrFQDrNY93d\nZDIZBgY+wkknncVTnvIGTjrpLAYGPsLI4CCHgAFgOxz+/Gh7Drj5K1+pq3yVx2h1svPQC3n/U59W\ntvztOhC/2V3Tja7s4BJf7a/dJjucXWL7TuB59Rxw3bp1zlbVtElq9vi4xZknrdkqZc4/8OTH+Oae\nI3h17lDJY3yzo4O7Fy5k+0MP0VNQPyNRQttmByLrN29+/NreMOnaPr+3l5NO+lP+/mmXcvDgcsIc\ns5CH7IILvsVL+Rg3cvgsy8L0GedHtfGBOluUqhmsf3vua9y74lI63tBTtPzt/N7Mpx9ZWqZ7tZEW\n2UZXdgiB4IZodmpxIRA8v67yKWhk9upMWNmhHFd20LRJcvb4OGUyGTYsXjwlT1qhgXSa82+7ra1v\nuJWUSph6yikv4PLu7oqv/7yxsaJ55uL+krB3716e9rRLOXjwQ0UfT7Oa27iyYkLbDatWsfaKK2o+\n//LlH+C66wYq7rds2QfYtq3yfs1SKWFuqS+BLzjllKquj0bfH40k9A155I4qGwiaR256zISVHaS2\nlPTs8XFKesLWanV2drJixdKiKRj2d3czQHidpVrcFixYwNIVK9quDi67bFPUElfcbvoZZiMrOFhy\nn0ZalCbGaJUPFeMco1UpYe6zD3yGVxz5SNEW6cu7uvjugQMMpNNlr49GP1fKXZ+VbN68viAQ3DAp\nEDzfz7yY2SInVWHH1q1w5plluz+2p1J0XHtt292I43bJ6tWsubJ5LTZJkdRu+ZNOOos77/wUpQOp\nDKv4Ha7gnpLHaKRFaevWHZx5JtFEh+JSqe1ce21HLHnMMpkMixdvKJIw9/E9eNucLi47+MuSxxhI\np3nTzp18dtN/suPqz9Gx9yFyC57M0nPO4q3nrXYt1FnEFjmpSZKePT5OMyFh63To7Ows2+LWroHe\n3r1zKN8a1sk3eSl/dsQe3jAnO+0tSu0+RqvSZIwUw7z+4MNlj3HKPbs4rfsC7n/ozWSz7wA64cEM\nO9aP8okq10LV7GUgJ1XBYKR+rVgUPOnaYTJIqUDy2GP38+CD5f+Cd7KRQ8/6v5z98bfVPdmg3Biu\nRgbrT5dS9fPZz46Qzb6n5PNCHsHS3c4AvblDHPlAiiyFwapLYKk6BnJSFQxG6veiM85g+MpPsoLS\nszZ3ckRDCWOTrB3GX5YLJF9/9BO4ltdwZ9lscDfw5je/vO4xfpUXZYfbbjs/tjFa5ernmDmdPJcf\ncSfFA+1Jhdn9AAAgAElEQVRq8wguZBfjJR53CSyVM+MDOdOPaDo0O0XATPYYR/EVjmMFvyi5zyAL\n+QOOamGp2kfck0EqBpK/fJAHuYU7y3yVmTNnG+ed9966z1/NouxA3YP1G1GxfrK/YhXfKVk/+TyC\nlb4E7qH0l8Bya6FqZmgk/Ui7JQSedvlATmrEdGS/n62+8pWb+Sa/zyrSbOHwhK1bSLGKNDdwGl/5\nys1xFjM2zV6iqZJqAsl+HmZuxznANg5P+buNOXPexfHH/7TuAfm1LMoeh+rq5+ekKL4WapjVO6fs\n83eSYjflvgS6BNZM19fXx7p16+p67oxvkZOmw+PZ7ymfQmImj2GpdzD+rl2PcCef504yXM0wXQyy\nkF3sYRG76SdLL9DJM3Z9oHUvpo3EPf6ymok8p5Pj3W+cz9yTf8DVV/8He/fOYcGCg5xzTg/nnffe\nhmZVVpPwN84Wqerq5yBdDDLO1BbTLL18cc4TWFFm1uogJ0Tvg1JcAkulGchJVaqU/X4mB3GNDMYv\nzAOWZQXjrCgyFij5N6p6A924x19WG0jmHniAD3z6fD4wzfF2uy/KXm39LJo3xO6D24tOxvjJgeMZ\nOPK3in4J3H7sU/j2A2nIlT6LS2CpHAM5qQaVUkjMRI0Oxu/v72HTptEKecCSfaNqJNCNe/xl3IFk\ntQl/TzhhDju2bi0ZKDcrfUu19bN05Wm87+yOEpMxPlD6S+App3BN9+WUnOmAS2CpPBMCSyqr0WTI\nEwlT15Aq0bWaTn+E225LZob4Rpcgq/b5zVrCK+5k19Uk/D25o48/Pu4OXvng/xZdHu+7Bw7wkiOP\nbMryea2oH5fAUl49CYEN5CSVNR0rM7zq6b/Dc+67j9cffJheDj5+ox1mDl+c8wR+cvzxXP+zHzah\n9M03HTf6C1eu5KioRa/Y+MvvNTFQadVauOXWGu3uvrzulREywEVz5vChg6VztTVS/kwmw4ULF/Lh\nTKbkPu/s7OSDe/bEthaqZg5XdpA07RodjJ/JZHjFkY+wdtKNuBNYwUFWHPwlA0f+VmITnk7Hqh/l\nxl+eV2JR9enKM9eKiTzlup4v7+ri2QfmQzpXtEXqqU/+FH+059cljz0KLC8TxEHj6VvGCQlQeqKf\nx+sn+im9OFn1GlkLVbPbjA/kzCMnNabRMVTV5kkb2raNefPmtd0SVXmlWpR+fffd0zLrtNT4yx1b\ntzY9z1wzJ/JUM8Yyl07z2bHzuOmmW6eMMbvxmg5evql0MukRoPyc18aWzxsdHuYvs1mWEILGDcAj\nwHxCUHc+MJLNNi3Pn2aHRvLIzXQ5SY3ZvmVLbnsqlctByZ/rU6nc9i1bij5/YNWq3L4yz81B7n2Q\ne9cxx+S2p1KP77sPcttTqdxAOp274KyzWvyqD3fBWWflBtLpouV71zHH5N5X4fXtg9zAqlV1nbua\n+mvk+M3W6PWzdtmyss9dW6FuHt9v2bK6yp+U+t+3b19u+5YtuYFVq3Jrly3LDaxaldu+ZUtu3759\nsZZLtQFqHg8241vkJDVmSW8vG7q6WFpmDNVIVxfnl0iGXKlrNgMcDaz9zW8O297KJarKqZzZP8ta\nys+7bGTWadx55hrVaNdzpRbh6ua8Vp512+wW12aqZtb0ezdubMqsXsXPQE6qQbNSHLSzRsdQVboR\njwJLKpShmUtUVVJN1/DphNdRanRTuUA3r9S1NeeEExK9zm+1geiv7767aHqRF51xRtn0LD3At4Bl\nZY5fKZAuFwgdmjuX9wN/X+b4cdZ/NV3X73zgAS593vM4/b77ak6Po/ZnICdVqZFcYUnXyBiqSnnS\nmj3GqVHVtCi9DDjn6KM58OijvOzgxKzcG+bMYfTEEznQ3Q1QMg/aP557bslr6+Cxx7Kmo4PLy8zA\nb+d1fqsZY/l+4NB3vgNnnjnl9X/3hBP4n1Sq5PWzhDBrdVmZCQ/lAum4W1wbVc0XjVdnMhy5a9dh\nXzTapcVbjTOQk6rQaFLc2axS12x1ef3j67qqtkUpdeAAP8jlGAGyQAp4Qi7HwVyOH4yMsGHx4qKB\n2qXHH09m717WT0pv8fi1df/9vGvOHDJRgFhMNS1+cakUyGcIi35/qFTX+j338M7OTi5ctIjT77uv\naIvwXQcOMBClZ6m1xbhVLa7NUu0XjQ2ULn+cLd5qnIGcVIVqZ17O1A/DRlojK3XN3nHUUWT27Wvb\nrsNqZ+0+f/9+puTeP3SIzK5dXDRnDmsntRjlAxV27aJ88gxYdugQ7zjuOM5+6KHErfNbKZAfJQRK\n5bw2m2X/Rz9Kx7x5JVuE620xrjYQevPRR5N77LG2q/+qx1CWeTzOFm81zkBOqsJ05ApLqulojSzX\nNfuWRx9ldOXK2JaoqqSaJbRuIIzVKqZSnrNqupZfnssxumIFHWefXXd6kLgSzlYK5C876ig+tW9f\n2WP0ZLNs+MpXWHvFFSXfX/Uun1dtIPScU0+l493vbrt1lqtOD1Tm8bgna6gxMz6QM4+cpkPSZw42\nYrpaI0vdaDOZTEOzYputmlm7o8C7SjxWKVCrtms5e999da/zG5aAmhsl3F1DPqXtpk2jdHVtoLt7\nf1OXgCoXyD/nQx+ic8eOss9v5nur2kDomGc+sy3XWa7mi8a3Kf1FA9p7ssxs0UgeuVkRyEmNinth\n8Tg1uzWyFSsLNKJS+b42Zw5HP/poyWujUqA2XekzSj43k2FsbG6RJbA6yWaXMj6+FBho+vjOUoH8\nyOBgrO+tagKhdphMUqpF9ZRTTuHyCl80vgGUWgAN4n99szEbwGT5BqeLL7645ufO+EBOmg5J+bBv\nhla0RjZzZYHpsH7zZvbu3cumyy7jsquvZs7evRxcsICec85h7u23897PfKbkcysFaj2EVrtXljl/\nI9fW8PAou3eXa4+B3bt7GB4ejWV5qLjfW43mSWyF8i2ql/PsA/PJpdNFv2h864QTuGPPHjrLrBUb\n5+ubzdkApouBnFSFJHzYN0urWiPrHePUCoU3m/PyN5sHH2R0/Xp2VEgPUinP2RJCjrKygVwD19bg\n4Eh08y8tm+1hcHBDLIFc3O+tdm8RnmhRXUOKYdK8g4XsYg+L2J3tZ3x8DaRzfHbsPG696aYpX4Te\nFaW3GYiu33Z6fWYDmB4GclIV2v3DvpnibjGJW8WbTYX0IJXynHUCt3d2sn7hQk67995pv7Z27apu\nFF7Yr/Xa4b01HS3CzeoeHB4eZd4917GKK+hnN71M1M8wmxiki2/d8wxuuOHFPHHeEU17fc0w27MB\nTJeOuAvQZNHSZdL0mI1jOTKZDBsWL54SyBQaSKc5/7bbZmQd7Ni6Fc48s2wg+/WODj63cGHJ9CDf\nO3CAF1fIc9asJZRWr76EK6/Md8eVkuHNb76UN72hJ7ZrO8nvrcIW2yWF3YPR37eR7sFzz72IIz/1\nKa7gnpL7vI1OUken+OPHMtN+/ma6ZPVq1lx5ZcXW/g2rVrH2iitaVaxYdXR0QI2xmS1yUg2a3f0X\nV4qIctqhxSRO1Uz2qCY9SDUtIs24tvr7e9i0aZRstnS36ckdr+WELbdz8NOXsKZgZYpvbdzIP594\nItmXvrQtA4F20OzuwR/9z51czL1l9/ljMqT2kbiVG2ZzNoDpZCAntYm4U0SU065dM61Q7c2mUnqQ\nuMYA9vYuoatrQzQ7tZgMr+Db/P0DhwerncCygwdZds89vHPPnoYDgXItbuWWKGv3Ae/N7h58cuZO\neqnwRYJkrtwwm7MBTCcDOakNtCpFRCPdV40GIkntOkv6zaazs5Pu7v3AQPQloYf8l4RUaoTjfutK\nXr+nfKDw6kyGoW3beO3rX19XGcrNTKy4RFkbtyjB9KXnKfX+OGHugRm7csNsH387XQzkpDbQihQR\ncU7zT3KKgZlws9m8eX1Bt/2Gw7rtv/jhB3jZ9vLPfxnw7ssuqyuQq9T1WM0SZe3aogTT0z1Y7v3x\npDlzeD9hZnMpSV25Ie4ZyzOFgZxmlKS2+jQ7RUSc0/yTnmJgptxsOjs7WbFi6ZTrZ9Oqc6oKRB68\n9Va2bt1R8/jNSl2P1SxR1q4tStB4i23F90c2y1rK5yJM6soNs3387XSZ8YGcS3TNHklu9Wl2iog4\np/knPcXATL/ZPEJ1K0vcc/9vOPNMah6/WanrsdolytqxRQkab7Gt5v1xOmEZuFJf4dp95YZyZvP4\n20KNLNE10+U0O+zbty83kE7nclDy54PpdG7fvn1xF7WoVasGcrCvXPFzsC+3atVAXccfWLUqt6/8\nwXP7IDewatU0v7J4zz2d9u3bl9u+ZUtuYNWq3Nply3IDq1bltm/Z0rbXVLVWL1+eu77C3+c6yJ3I\n8pK7pNMfLFkPa5ctK3vsgejvn9Tro9HPnmrfH/2Qu76grvZFv18Muf9T4fnt/NmnwwE150yb8S1y\nmh2S3upTTYqIVGqE/v7y4+hKiXOa/0xJMdDOK0804o/+4i/4xrZtZVeW2AH8gr8o+Xi58ZuVuh6b\nvURZszXaYlvt++M5hORiGwitmPMJdfcuYE1HB2uPO46lDz4441qMVZmBnGaEZi/s3myVU0RAV9cI\nvb3n13X8OGdeJn3W50zXt3w5V3V2MpDJ0EMIDh4PBKKfm+nkAMtLHiOb7eGzn72Uo3hsyvjUF51x\nRtmux2YvUdYKjXQPVvv+OIbQtVrsE+Lpz3gG542NFV2iazZ1T85WBnKaEZLe6lMpRURX1wjd3Y/V\n/YEc58zLmTDrcybr7Ozkd173Oh678UZu+fnPueHQIbJACjjmiCP4WeoYfrr/1ZQLNZ7LuRzzX9vg\nP/9hyvjU755wAv+TSpX8+zd7ibJWqbfFtpr3xzc7Orh74UK2l1g55LHubhYsWDAjW4xVmUt0aUaY\nKUu9NGtlhziX2ZrtS3wlRakZ35s+O8KnPvUeSgdyGVbxO2WXkHpnZydHL1jA6ffd1/Ilytpdte+P\nfIvbbKuf2aaeJboM5DQjVLMe5vZUio5rr52131gvXLmSo6JZvaVups3MIxfXudWYrVt3RLNVi3f7\np9jKtZzBijLZ4LanUuzfvJl58+YZiBTh+0N5BnJTGcjNErb6VKeRPHuN5uiL+/mqTyaTYfHiDUVW\nHQnSrOY2kt8aHjevb4GBXDEGcrPIHz/96TzrvvtYfvAgL2NiwPYNwLY5c/jp8cfz3z/7WbyFTKjC\nHH1LCsdARS0Gzc7RF/f5Z7uwDvBRRcdv9sx5GzsfvaviMT6wbBkD27Y1vaxSktUTyDnZQTNCJpPh\nJUceyZqDBxll6hT99QcP8pEjj2zb1QPaWdwrM8R9fpVf4uvGa04js+kuZyVLMbFFTjOCY+SaJ+66\njfv8Ks+/jzR9bJHTrFVtHrkPXXPN4/s3YxzKTBzn0qocfaXq7pvXXMO7E5wjcKabKWvRSkk14wM5\n11qdHarNI3fH5s28/JprmrIWa5LXei2nFTn6ytXdHXPmJDpH4Ew309eilVqhkbVWZ0Ugp5mv2uzo\nix999LDM6NM1zmomj+Nq9soMlepuJJt1ZYg258LnUmPyDU4XX3xxzc91jJwSpVT326OPPsq8lStZ\nks0ySlhWqHCywxLC7NUjKb7EDTQ2jmcmjxNq9murdPwdwEFgWZljJLVuNXPMxGEVaj3HyGlGK9f9\n9t0TTmDo4EGWEgK3NUykH8nPYv0ecHWZ4zcyzirpa72W0+wxUJXqbgnwz5QP5ByDpTjN1GEVSgYD\nOSVCxa7Le+7hoY6OxwO4wx6Pft5Z4RyNjLNK+lqv5TR7DFSluusEssA7583jtQcP1n1+W0zUDDN5\nWIWSwUBOiTA6PEzP7t1l91mRyzFK6a7T10DZxxsZZ9XscWRxa+YYqGrq7r3ApX/yJ3S84Q11nb8V\nLSazPVBM+uuvt/zVfDb17N7N6PBw4lrjpXaQ08wwsGpVbh/kcmV+9kFuoIHHr0+lctu3bKmrfNu3\nbMltT6XKlq+R489kza67ffv25QbS6bLH/2A6ndu3b1/dr+GCs87KDaTTue2p1OPX6T7IbU+lcgPp\ndO6Cs86q+9hJkPTX30j5q/5sWrWqha9ISQXUPLC/Uovcx+s46DtqLYRUSdVdlw083sg4K3Np1a/Z\nddfsFpPZ3rWW9NffaPln8rAKJUOlQO6v6jimgZymXdVdlxUev+Poo9n+2GPTPs7LXFr1a3bdNXsi\nymzvWkv662+0/DN9WIXaX6VA7lktKYVUQU9/P6ObNpVNgXEDYcZqKSOpFG/59KfpmDevKbmuzKVV\nv2bWXbNbTGbyjOVqJP31N1r+aj6bRlIpevr7Gy6rVEylQG68FYWQKqmm+23bnDmsP3iw5OMjXV2c\nv3w5nZ2dTbuh5I/djjesdtesumt2i8ls71pL+utvtPwOq1Dc6p21ugDYO50FmS2SPrMrLtV0v/30\nwAE+cuSRdm3qMM1uMZntXWv51w+UTMad3y9OpT5755xwQkN/P4dVKG7VZg+eB7wH+EPguYT36CPA\nncAXgUuB/c0oYIOiSSDtoTAFwpLCFAjRm92kkZXt3buXqy67gh1Xf46OvQ+RW/Bklp5zFm89bzUL\nFiwwUNYUmUyGDYsXTxnMXmggneb8226r6xqZKat61Pve2bF1K//12teyKJd7PHArTMY9Auzq6OBP\nvva12F5/uc/eHccey54HHuDyMveK7akU+zdvZt68eSXrp9mfPX62zQ71rOxQzc4vBgaBNGEG627g\nXuAEoCs6xt1AP3BzLSdvgbYJ5Jp9M5kNVq68kLGxueze3UM2O3G7SKVG6eoaobt7P5s3r4+7mGpD\nF65cyVHRjbxUi0m9X6Jmwnu7kS+Ze/fu5dKnPY0PlRnW8K45c3jvL37BggULmvMCyqjm7/OuaFhG\nqb9Of2cnL1y4kNPuvTeWL+E2AswezViiqxP4LPAU4ALgo8BvCh5/AmGW6t9F+/0e5TM8zFpJn9kV\nt0wmw9jYXMbH1056pJNsdinj40uBgbZNcaB4NXMyRdK71hpNv3HrTTex7NChsudYdugQt950Uyyf\nbdV89i47dIh3HHccZz/00JS/387jjye9dy8X3HPPYc9pVXqVpKd3UfNVCuTeS5i5ugK4vsjjDwOX\nEFrPvw68D7hwOgs4UyR9ZlfchodH2b273JxU2L27h+HhUVasKLV2g2azZk5ESfKM5Ua/ZI4MDrKm\nQs/Hy3O52D7bqvnsfXkux+iKFXScffaUv9+SRx9l3sqVZZ/fzC/hNgKokkqB3GnA5ygexBXaQeh+\nPW06CjWd1q1bR19fH319fS05X6lxDL++++5Ez+yK2+DgCNnsmrL7ZLM9DA5umNWBnONo4pPUGcuN\nfsls91mr1ZYve999Rf9+l6xeHeuXcBsBZoehoSGGhobqem6lQO4FhBa3anwHeH9dpWiidevWtexc\n5dZzPDR3Lu8H/r7M82fyzLZG7dr1COXnlQF0RvvNTq1YT1QzT6OBWLvP2m20fHEHqnGfX62Rb3C6\n+OKLa35upUBuP2HGajXmAo/VXIIZouI4hmyWtVD2A8WkkaUtWjSf8rUHkIn2m31aMY6m0dY+Wwvb\nU6OBTrsnxG20fHEHqnGfX8m3FfhylfteG+3fTlq20G01C39/HXLbm7hw90y2Zcv2XCq1vVz15lKp\n63NbtmyPu6ixaPbC840uip70RdVnskavnX379uUG0umyz4/zs63R8jX7vVVJ3OdXa1H7+vYcUeHx\n64HXAW+ssN8bgDOoPJZuxhoZHGRJhXEMLwM+cfTRbE+lHk+gmSHkKBpIp9t6ZlvcenuX0NU1Unaf\nrq4RenuXlN1npqrm+uvJZhkZHKz52IWtfUujLluYaO1bOz7OUWNjZDKZpjxfzbWkt5eRrq6y+4x0\ndbGkxMoEj8/aTafb8rOt0fI1Wj+Nivv8an+Vulb/BfgTYCNwOrCBkAT4ECEIfC6wBngL8L1o/1mp\n2nEMzzn1VDre/e7EzWyLW2dnJ93d+4GBKI9cDxN55EaiPHKPzdo6bOY4mkZnzTnrrr1NR/qUdp+1\n20j54k4vE/f51f6qSTr3bODTwKnR7wcIy3MtYCIQ/A7wp8BPpruADYpaKpvvktWrWXPllRXHMWxY\ntYq1V1zRkjLNRJlMhuHhUQYHR9i16xEWLZpPf38Pvb1LZvUHWTOvv0aP7XsjGRzDWF7c9RP3+dUa\nzVrZIb/fucDrgZOBpwM/A24HvgRcVctJW6hlgdxMWaZHydTM6+8Dy5czcN11lfdbtoyBbdum/fmt\nUulG6Y1UUrM1Y2WHvBzwyehHRSzp7WVDVxdLyywDM9LVxfmOY1ATNPP6a3TWXBJm3VVK3fLdAwd4\nyZFHmtpFUtupNNlBVWr3Ab+a2Zp5/fX09zOaSpXdp1z6hkaf32zVTMZ41n33scbJGpLaUE3NdwnU\nsq7VPLtfFKdmXH+NLgrf7ovKV9MtvY0wOHgeYT3CR4D5QA+whBCIOmxCUqOa2bWqKiV1mR7NDM24\n/hqdNdfus+6qWQJpCMgCryFM03+8a5UwlX9fNsvRLpEkKQa2yEmqykxd2aHSZIwMIVhbW+YYA8DD\nS5dy6fbtdZfDGdmSmjlrNakM5CSVVSk9yg7CbK9XljnG9cB/LV/Of3z963WVYeXKCxkbmxvlSFzC\nRI7E0ShH4n42b15f17ElJYddq5JUo0prcY4QulPL+X3gCx3lP3tLtUi+4JRTGBuby/j45Da/TrLZ\npYyPLwUGGlonV9LMZSAnaVarlLrlEcqnTiF6/IllxtmVTW9y7FOY98BJlOu83b27h+HhUVasWFqh\nJJJmG9OPSJrVKqVuuePoo6mUWCQDHPPMZxZ/rEJ6k0vu/wW/nxuPjlJcNtvD4GD5tYYlzU62yEma\n9cqtxfmWRx9ldOXKsulJyuXBq2at2X7u5WqGyVJq1msnu3Y9Uu3LkTSLGMhJEqVTt2QymYZWzagm\nvcnpZOlikPGSgVyGRYvmlz2GpNnJQE6Symg0D94ju3ZVNcZuIbsYL/F4KjVCf39PvS9B0gxmICdJ\nFZTrej2/Qh68atea3UPptWa7ukbo7T2/7vJLmrnMIydJTVTNEmBbO47gzxb+X+5/6M1ksz1M5JEb\nifLIPWYeOWkWMCHwVAZykmJV7Vqz542NcdNNt7qygzSLmRC4iHXr1tHX10dfX1/cRZHadpkqNU+1\nY+wWLFjAihVLzRUnzUJDQ0MMDQ3V9Vxb5KQWKUwKu6QwKWx0M9/f3c36zZvjLqaaxCBeUiV2rU5l\nIKe2UG332vm33eZNXZJmqXoCOVd2kFqgmqSwPbt3Mzo83KISSZJmAgM5qQVGBgdZUiEpbE82y8jg\nYItKJEmaCQzkpBaoNinsI7t2taI4kqQZwkBOaoF8UthyMtF+kiRVy0BOaoGe/n5GU6my+5RbeF2S\npGIM5KQWWNLby0hXV9l9Rrq6WFJi4XVJkoqZ8QmBpXbQ6MLrkiQVYx45qYVMCitJKsWEwFMZyEmS\npEQwIbAkSdIsYiAnSZKUUAZykiRJCWUgJ0mSlFAGcpIkSQllICdJkpRQBnKSJEkJZSAnSZKUUAZy\nkiRJCWUgJ0mSlFAGcpIkSQllICdJkpRQBnKSJEkJZSAnSZKUUAZykiRJCWUgJ0mSlFAGcpIkSQll\nICdJkpRQBnKSJEkJZSAnSZKUUAZykiRJCWUgJ0mSlFAGcpIkSQl1ZNwFaMBc4B+A+cCXgK3xFkeS\nJKm1khzInQaMAdcA/46BnEQmk2F0eJiRwUEe2bWL+YsW0dPfz5LeXjo7O+MuniRpmrVbIHcV8Drg\nAeD5Bdt7gcsJ5f0Y8PHo8ZHo8fktLKPUli5cuZK5Y2P07N7NmmyWTiADjG7axIauLvZ3d7N+8+a4\niylJmkbtNkbuk8CKIts/CrwNeBXwV8AC4FbgWdHjmZaUTmpTmUyGuWNjrB0fZ2kUxAF0AkuzWdaO\nj3PU2BiZjG8VSZpJ2i2Q+ybw0KRtT4r+HQbuAbYBpwLfBk4htM79d6sKKLWj0eFhenbvLrtPz+7d\njA4Pt6hEkqRWaLeu1WK6gdsLfv8hsAT4KvCuSk9et27d4//v6+ujr69veksntYGRwUHWZLNl9+nJ\nZtkwOMjSFcUavSVJrTY0NMTQ0FBDx0hCINeQwkBOmqke2bWLSlMZOqP9JEntYXID08UXX1zzMdqt\na7WYMeDkgt8XA6MxlUVqS/MXLao4UDQT7SdJmjmSEMj9Kvq3F0gDy4AbYyuN1IZ6+vsZTaXK7jOS\nStHT39+iEkmSWqHdArlrCJMYngv8DHhLtP1vCelHrgf+DdgbS+mkNrWkt5eRrq6y+4x0dbGkt7dF\nJZIktUK7jZE7u8T2ncDzWlkQKUk6OzvZ393NAGF2ak9BHrmRVIqRri4e6+42KbAkzTAdcRegyXIX\nXXSRs1U1a7iygyQlT372ajTZoabYbMYHcrlcLu4ySJIkVdTR0QE1xmbtNkZOkiRJVTKQkyRJSigD\nOUmSpIQykJMkSUqoOXEXoMnW5f+TTqfjK4UkSVIJQ0NDbNy4kZ07dwLUtE6Xs1YlSZLagLNWJUmS\nZhEDOUmSpIQykJMkSUooAzlJkqSEMpCTJElKKAM5SZKkhDKPnCRJUozMI1eaeeQkSVIimEdOkiRp\nFjGQkyRJSigDOUmSpIQykJMkSUooAzlJkqSEMpCTJElKKPPISZIkxcg8cqWZR06SJCWCeeQkSZJm\nEQM5SZKkhDKQkyRJSigDOUmSpIQykJMkSUooAzlJkqSEMpCTJElKKAM5SZKkhDKQkyRJSiiX6JIk\nSYqRS3SV5hJdkiQpEVyiS5IkaRYxkJMkSUooAzlJkqSEMpCTJElKKAM5SZKkhDKQkyRJSigDOUmS\npIQykJMkSUooAzlJkqSEMpCTJElKKNdalSRJipFrrZbmWquSJCkRXGtVkiRpFjGQkyRJSigDOUmS\npIQykJMkSUooAzlJkqSEMpCTJElKKAM5SZKkhDKQkyRJSigDOUmSpIQykJMkSUooAzlJkqSEMpCT\nJEgJkgYAAAivSURBVElKKAM5SZKkhDKQkyRJSigDOUmSpISaE3cBmmxd/j/pdDq+UkiSJJUwNDTE\nxo0b2blzJ8DFtTy3ozlFahu5XC4XdxkkSZIq6ujogBpjM7tWJUmSEspATpIkKaEM5CRJkhLKQE6S\nJCmhDOQkSZISykBOkiQpoQzkJEmSEspATpIkKaEM5CRJkhLKQE6SJCmhDOQkSZISykBOkiQpoQzk\nJEmSEspATpIkKaEM5CRJkhLKQE6SJCmhDOQkSZISykBOkiQpoQzkJEmSEspATpIkKaEM5CRJkhJq\nTtwFaLJ1+f+k0+n4SiFJklTC0NAQGzduZOfOnQAX1/LcjuYUqW3kcrlc3GWQJEmqqKOjA2qMzexa\nlSRJSigDOUmSpIQykJMkSUooAzlJkqSEMpCTJElKKAM5SZKkhDKQkyRJSigDOUmSpIQykJMkSUoo\nAzlJkqSEMpCTJElKKAM5SZKkhDKQkyRJSigDOUmSpIQykJMkSUooAzlJkqSEMpCTJElKKAM5SZKk\nhDKQkyRJSigDOUmSpIQykJMkSUooAzlJkqSEMpCTJElKKAM5SZKkhDKQkyRJSigDOUmSpIQykJMk\nSUooAzlJkqSEMpCTJElKKAM5SZKkhDKQkyRJSigDOUmSpIRKciD3TOAK4HNxF0SSJCkOSQ7k7gZW\nx12ImWxoaCjuIiSa9Vc/664x1l9jrL/GWH+t1Q6B3FXA/cD3J23vBX4E/Bj461YXSr4ZG2X91c+6\na4z11xjrrzHWX2u1QyD3SWBFke0fBd4GvAr4K2ABcA6wATihZaWTJElqU+0QyH0TeGjStidF/w4D\n9wDbgFOBq4E1wL3AscBlwO8B72lJSSVJktpIR9wFiKSBLwPPj35/FbAKODv6/TygC7igxuP+BHj2\nNJRPkiSp2e4CnlPLE45sUkHaRU2VIUmSlCTt0LVazBhwcsHvi4HRmMoiSZLUlto1kPtV9G8vodt1\nGXBjbKWRJElSUdcQJi/sB34GvCXafjoh/chPgHfEUzRJkiS1QrF8dE8AvgTsAr4IHBNDuZLg6cA3\ngNuAIeCN0XbrrzrzCC3GtxCGAKyJtlt/tZkD3EyY+ATWXy3GgVsJ9fedaJv1V52jgU8BdwI/JGRI\nsO6qcxLhmsv//IrQ8HIM1l+1/gz4NvBd4F+ibTVdf+3atVqPYvno/oJQEb8N/Jww+1VTZQnBx2Kg\nHxggXEjWX3UeBV5BSIVzOmHG9W9j/dXqbwg30lz0u/VXvRzQB7wIeGm0zfqrzsWEenpB9HM71l21\n7iBccy8CXgJkgC8Af4n1V41jgfcTho91A88FXs0sv/7SHN4iN0i4uQK8GNdlrdaXgaVYf/V4CuFG\nsAjrrxYnAtcTAuJ8i5z1V727CddeIeuvOrcA8ydts+5qt5yQFxasv2rNJ7Smn0BoGR4itAjP6vpL\nc3ggdw+h2wugM/pd5T0H+CmhKdf6q94RwP8AB4C3R9usv+p9jvCt/nQmAjnrr3o/JVx/XwTOjLZZ\nf5WdSPjitZEwPOI9hJurdVe7qwgtcWD91eI1wGPAw8Al0baa6m8mda0W0y4Jj5PiCcB/EbpZf4P1\nV4tDwAsJgfBfEoIS6686ZwAPEMbYFNaZ9Ve90wjX3/uAjwBPw/qrxjxCd9bnCV3Ti4GVWHe1Ogr4\nAyZajqy/6iwEPgH8DqEhqofweVhT/c30QG4MeF70/+dFv6u4FOHD7GrCIEuw/uoxDnyN0Dxu/VXn\n9wmtSHcTZrEvJVyH1l/17ov+/RFwLeGmav1V9hPCOK8vA48Qrr8VWHe1eg1hsP6e6HfrrzovJUyQ\n+wnwv4RA+OXUWH8zPZC7EXgroan8rZhUuJQO4ErgB0zMmgHrr1oLgN+K/v8UwliRL2H9Vev9hJnT\nzwTeAOwAzsH6q1YnoTUdwjf8VwNbsf6q9WPCF68jgNcRxmpad7U5mxAE51l/1fkmcAph0sNcQkC8\njVlcf8Xy0TmFvDovI3QN3sLENPIVWH/Vej7wPcIYpa8Db4q2W3+1O53QogTWX7WeSXjv3gJsJ3zw\ng/VXrecSbpS3AB8mDDq37qp3NLCXiS8TYP3V4lxgJ6HV7YOELxTWnyRJkiRJkiRJkiRJkiRJkiRJ\nkiRJkiRJkiRJkiRJkmasJxOWZzoE/GnMZZEkSVIN3k4I4h4mLA8mSZKkhLiZsLzaB4GDhCWvJEmS\n1OZeTGiNeweQJgRyHyyx798Q1uN8iLCw9QpgXfT8RZP2fQLwkWj/h4EfAhtwzURJkqRp86+E8XFP\njn6/jrBQdcek/T5ECNgeANYDX46et5OpgdxTgXFgH/AFQgD4JWA/cBtw1PS/DEmSpNllHqF17ZqC\nbX9CCMxWFGx7FvAYcAMwv2D7K6J9D3J4IPcpIAucPOl8r4/2/7tpKLskSdKs9kZCYPWqgm1HAXuA\nzQXb/jra791FjnEvhwdyKUIQdy3wFGBBwc/TgL3A16btFUiSJM1S2wldpb8NPKfg5zPAo8Cx0X4b\nCIHcS4ocYxOHd60+L/q93M/t0/9SJM1WR8ZdAEmKwTMJXaMAd5TY5xzgozUeNxX9+yXg4yX2eaTG\nY0pSSQZykmajt0T/rgZ+OemxDmAAeCshkPtptP2VwHcn7fsqIFfw+x2ErtYXYU46SZKkaXcEYWbq\nLWX2uZDQDXoKE5MdvgV0FuyzlOKTHT4dbX99keN2EMbLSZIkqQ4rCIHWBWX2WRzt84no93+Oft9D\nyDM3Of3I0wueeyzw42j7D4CLCBMmNgB3EYJESZIk1eFzhFa0xRX2ux14EJgb/f63hFa8XwLXA68B\nriIEbJ2TnjuPEMDdGO3/EGH1iA1MTUsiSZKkGNxHmPkqSbE4Iu4CSFICzCuybRlhFYfrWlwWSZIk\n1eA84NvA+wldrCOELtVfEXLPSZIkqU11A1uBXxCSBf8Y+CyQjrFMkiRJkiRJkiRJkiRJkiRJkiRJ\nkiRJkiRNg/8PYQGxsmOEUjEAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
Regression Fit
\n", "First example using `statsmodels`.\n", "\n", "* Typically a lot of online resources\n", " * Good one on regression is [Connor Johnson's blog post](http://connor-johnson.com/2014/02/18/linear-regression-with-python/)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import statsmodels.formula.api as sm\n", "df['Eins'] = 1\n", "Y = np.log(df.OI)\n", "Y.name = 'np.log(OI)'\n", "X = df[['Age','Gender', 'Eins']]\n", "result = sm.OLS( Y, X ).fit()\n", "result.summary()\n", "#import statsmodels.graphics.regressionplots as rp\n", "#rp.plot_fit(results, 1)" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: np.log(OI) R-squared: 0.262
Model: OLS Adj. R-squared: 0.247
Method: Least Squares F-statistic: 17.23
Date: Sun, 23 Nov 2014 Prob (F-statistic): 3.96e-07
Time: 12:01:08 Log-Likelihood: -61.671
No. Observations: 100 AIC: 129.3
Df Residuals: 97 BIC: 137.2
Df Model: 2
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [95.0% Conf. Int.]
Age 0.0162 0.004 4.603 0.000 0.009 0.023
Gender 0.3189 0.116 2.757 0.007 0.089 0.548
Eins 0.8292 0.178 4.666 0.000 0.477 1.182
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 3.016 Durbin-Watson: 2.007
Prob(Omnibus): 0.221 Jarque-Bera (JB): 2.772
Skew: -0.408 Prob(JB): 0.250
Kurtosis: 2.975 Cond. No. 200.
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 7, "text": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: np.log(OI) R-squared: 0.262\n", "Model: OLS Adj. R-squared: 0.247\n", "Method: Least Squares F-statistic: 17.23\n", "Date: Sun, 23 Nov 2014 Prob (F-statistic): 3.96e-07\n", "Time: 12:01:08 Log-Likelihood: -61.671\n", "No. Observations: 100 AIC: 129.3\n", "Df Residuals: 97 BIC: 137.2\n", "Df Model: 2 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "------------------------------------------------------------------------------\n", "Age 0.0162 0.004 4.603 0.000 0.009 0.023\n", "Gender 0.3189 0.116 2.757 0.007 0.089 0.548\n", "Eins 0.8292 0.178 4.666 0.000 0.477 1.182\n", "==============================================================================\n", "Omnibus: 3.016 Durbin-Watson: 2.007\n", "Prob(Omnibus): 0.221 Jarque-Bera (JB): 2.772\n", "Skew: -0.408 Prob(JB): 0.250\n", "Kurtosis: 2.975 Cond. No. 200.\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
Load in R Datasets
\n", "\n", "`statsmodels` seems designed for those with experience of statistics and `R`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "import statsmodels.api as sm\n", "import statsmodels.formula.api as smf\n", "\n", "# Fit regression model (using the natural log for the dependent variable)\n", "results = smf.ols('np.log(OI) ~ Age + Sex', data=df).fit()\n", "\n", "# Inspect the results\n", "results.summary()\n" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: np.log(OI) R-squared: 0.262
Model: OLS Adj. R-squared: 0.247
Method: Least Squares F-statistic: 17.23
Date: Sun, 23 Nov 2014 Prob (F-statistic): 3.96e-07
Time: 12:01:16 Log-Likelihood: -61.671
No. Observations: 100 AIC: 129.3
Df Residuals: 97 BIC: 137.2
Df Model: 2
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [95.0% Conf. Int.]
Intercept 0.8292 0.178 4.666 0.000 0.477 1.182
Sex[T.Male] 0.3189 0.116 2.757 0.007 0.089 0.548
Age 0.0162 0.004 4.603 0.000 0.009 0.023
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 3.016 Durbin-Watson: 2.007
Prob(Omnibus): 0.221 Jarque-Bera (JB): 2.772
Skew: -0.408 Prob(JB): 0.250
Kurtosis: 2.975 Cond. No. 200.
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 8, "text": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: np.log(OI) R-squared: 0.262\n", "Model: OLS Adj. R-squared: 0.247\n", "Method: Least Squares F-statistic: 17.23\n", "Date: Sun, 23 Nov 2014 Prob (F-statistic): 3.96e-07\n", "Time: 12:01:16 Log-Likelihood: -61.671\n", "No. Observations: 100 AIC: 129.3\n", "Df Residuals: 97 BIC: 137.2\n", "Df Model: 2 \n", "Covariance Type: nonrobust \n", "===============================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "-------------------------------------------------------------------------------\n", "Intercept 0.8292 0.178 4.666 0.000 0.477 1.182\n", "Sex[T.Male] 0.3189 0.116 2.757 0.007 0.089 0.548\n", "Age 0.0162 0.004 4.603 0.000 0.009 0.023\n", "==============================================================================\n", "Omnibus: 3.016 Durbin-Watson: 2.007\n", "Prob(Omnibus): 0.221 Jarque-Bera (JB): 2.772\n", "Skew: -0.408 Prob(JB): 0.250\n", "Kurtosis: 2.975 Cond. No. 200.\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] } ], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "* There is even provision for loading in `R` datasets." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "import statsmodels.api as sm\n", "import statsmodels.formula.api as smf\n", "\n", "# Load data\n", "dat = sm.datasets.get_rdataset(\"Guerry\", \"HistData\").data\n", "\n", "# Fit regression model (using the natural log of one of the regressors)\n", "results = smf.ols('Lottery ~ Literacy + np.log(Pop1831)', data=dat).fit()\n", "\n", "# Inspect the results\n", "results.summary()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: Lottery R-squared: 0.348
Model: OLS Adj. R-squared: 0.333
Method: Least Squares F-statistic: 22.20
Date: Sun, 23 Nov 2014 Prob (F-statistic): 1.90e-08
Time: 12:01:22 Log-Likelihood: -379.82
No. Observations: 86 AIC: 765.6
Df Residuals: 83 BIC: 773.0
Df Model: 2
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [95.0% Conf. Int.]
Intercept 246.4341 35.233 6.995 0.000 176.358 316.510
Literacy -0.4889 0.128 -3.832 0.000 -0.743 -0.235
np.log(Pop1831) -31.3114 5.977 -5.239 0.000 -43.199 -19.424
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 3.713 Durbin-Watson: 2.019
Prob(Omnibus): 0.156 Jarque-Bera (JB): 3.394
Skew: -0.487 Prob(JB): 0.183
Kurtosis: 3.003 Cond. No. 702.
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 9, "text": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: Lottery R-squared: 0.348\n", "Model: OLS Adj. R-squared: 0.333\n", "Method: Least Squares F-statistic: 22.20\n", "Date: Sun, 23 Nov 2014 Prob (F-statistic): 1.90e-08\n", "Time: 12:01:22 Log-Likelihood: -379.82\n", "No. Observations: 86 AIC: 765.6\n", "Df Residuals: 83 BIC: 773.0\n", "Df Model: 2 \n", "Covariance Type: nonrobust \n", "===================================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "-----------------------------------------------------------------------------------\n", "Intercept 246.4341 35.233 6.995 0.000 176.358 316.510\n", "Literacy -0.4889 0.128 -3.832 0.000 -0.743 -0.235\n", "np.log(Pop1831) -31.3114 5.977 -5.239 0.000 -43.199 -19.424\n", "==============================================================================\n", "Omnibus: 3.713 Durbin-Watson: 2.019\n", "Prob(Omnibus): 0.156 Jarque-Bera (JB): 3.394\n", "Skew: -0.487 Prob(JB): 0.183\n", "Kurtosis: 3.003 Cond. No. 702.\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] } ], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "
`sklearn`
\n", "\n", "* In machine learning the `scikit-learn` would be more familiar." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from sklearn.linear_model import LinearRegression\n", "model = LinearRegression()\n", "\n", "X, y = df[['Age', 'Gender']], np.log(df['OI']).values[:, None]\n", "\n", "model.fit( X, y )\n", "\n", "score = '{0:.3f}'.format( model.score( X, y ) )\n", "print model.coef_" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0.01620833 0.31889832]]\n" ] } ], "prompt_number": 10 }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
My Own Preference
\n", "\n", "* When teaching my own preference is to derive the maximum likelihood solution and have students implement it.\n", "* Solve\n", "$$\n", "\\mathbf{X}^\\top \\mathbf{X} \\boldsymbol{\\beta} = \\mathbf{X}^\\top \\mathbf{y}\n", "$$\n", "for $\\boldsymbol{\\beta}$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "beta = np.linalg.solve(np.dot(X.T, X), np.dot(X.T, y)) \n", "beta = pd.DataFrame(beta, index=X.columns)\n", "beta" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
0
Age 0.016208
Gender 0.318898
Eins 0.829203
\n", "
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 15, "text": [ " 0\n", "Age 0.016208\n", "Gender 0.318898\n", "Eins 0.829203" ] } ], "prompt_number": 15 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Although Darren Wilkinson pointed out to me at the meeting, that this is more correctly done with [QR decomposition](http://en.wikipedia.org/wiki/QR_decomposition) (see [this page](http://www.stat.wisc.edu/~larget/math496/qr.html) for a description).\n", "\n", "$$\\mathbf{X}^\\top \\mathbf{X} \\boldsymbol{\\beta} = \\mathbf{X}^\\top \\mathbf{y}$$\n", "$$(\\mathbf{Q}\\mathbf{R})^\\top (\\mathbf{Q}\\mathbf{R})\\boldsymbol{\\beta} = (\\mathbf{Q}\\mathbf{R})^\\top \\mathbf{y}$$\n", "$$\\mathbf{R}^\\top (\\mathbf{Q}^\\top \\mathbf{Q}) \\mathbf{R} \\boldsymbol{\\beta} = \\mathbf{R}^\\top \\mathbf{Q}^\\top \\mathbf{y}$$\n", "\n", "$$\\mathbf{R}^\\top \\mathbf{R} \\boldsymbol{\\beta} = \\mathbf{R}^\\top \\mathbf{Q}^\\top \\mathbf{y}$$\n", "$$\\mathbf{R} \\boldsymbol{\\beta} = \\mathbf{Q}^\\top \\mathbf{y}$$\n", "So if we let $\\mathbf{z} = \\mathbf{Q}^\\top \\mathbf{y}$,\n", "$$\\mathbf{R} \\boldsymbol{\\beta} =\\mathbf{z}$$\n", "\n" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "import scipy as sp\n", "Q, R = np.linalg.qr(X)\n", "beta = sp.linalg.solve_triangular(R, np.dot(Q.T, y)) \n", "beta = pd.DataFrame(beta, index=X.columns)\n", "beta" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
0
Age 0.016208
Gender 0.318898
Eins 0.829203
\n", "
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 16, "text": [ " 0\n", "Age 0.016208\n", "Gender 0.318898\n", "Eins 0.829203" ] } ], "prompt_number": 16 } ], "metadata": {} } ] }