{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Nonlinear Least Squares\n", "====\n", "\n", "## Unit 12, Lecture 2\n", "\n", "*Numerical Methods and Statistics*\n", "\n", "----\n", "\n", "#### Prof. Andrew White, April 17 2018" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Goals:\n", "---\n", "\n", "1. Be able to apply the same analysis from 1D/ND OLS to NLS\n", "2. Compute the F-matrix and understand its use in error analysis\n", "3. Distinguish between when linearized OLS is possible and NLS is required" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "%matplotlib inline\n", "import random\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.stats\n", "import scipy.linalg as linalg\n", "import matplotlib" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Linearizing An Exponential\n", "---\n", "\n", "We previously learned how to linearize polynomials into OLS-ND. What about a more complex equation, like:\n", "\n", "$$y = \\beta_0 e^{-\\beta_1 x^2} + \\epsilon$$\n", "\n", "Well, you could of course just do this:\n", "\n", "$$\\ln y = \\ln\\beta_0 - \\beta_1 x^2 + \\epsilon$$\n", "\n", "and then choose our $x$ matrix to be $[1,-x^2]$\n", "\n", "What is wrong with this?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "**We changed our assumption of noise!** To do the math I just above, it should have been that our original equation was:\n", "$$y = \\beta_0 e^{-\\beta_1 x^2}\\times e^{\\epsilon}$$ so that after taking the log, we ended up with that above. But that equation doesn't our assumption that we have normally distributed 0-centered noise." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Can we neglect our assumption of linear normal noise?\n", "===" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt4FOX5//H3nQMQCBKRABJEUqEcBAkYDpqKKCCiVhBoxUPVfquoqK21paK11VqrVPxVa6VSFEtb8VRBRKGggIiiVECgokAFxJKAcpKTBAnJ8/tjEkhCIBuyu7Oz+3ld1167Ozu7cy/RTybPPHOPOecQEZH4kuR3ASIiEn4KdxGROKRwFxGJQwp3EZE4pHAXEYlDCncRkTikcBcRiUMKdxGROKRwFxGJQyl+bbhJkyaudevWfm1eRCSQli5dus05l1nder6Fe+vWrVmyZIlfmxcRCSQz+zyU9TQsIyIShxTuIiJxSOEuIhKHfBtzFxH/FBUVkZ+fz/79+/0uRY6iXr16tGzZktTU1ON6v8JdJAHl5+fTsGFDWrdujZn5XY5U4pxj+/bt5Ofnk52dfVyfEahwn7asgLGz17BpZyEtMtIYNaAdg7tm+V2WSODs379fwR7DzIyTTjqJrVu3HvdnBCbcpy0r4K6pH1FYVAxAwc5C7pr6EYACXuQ4KNhjW21/PoE5oDp29ppDwV6msKiYsbPX+FSRiEjsCky4b9pZWKPlIhIc9913H4888shRX582bRqffPJJFCsKvsCEe4uMtBotF5HwmbasgLwx88gePYO8MfOYtqwguttXuNdYYMJ91IB2pKUmV1iWlprMqAHtfKpIJDGUHe8q2FmI4/DxrtoG/O9+9zvatWtHv379WLPGG1596qmn6N69O126dGHo0KHs27eP9957j+nTpzNq1ChycnJYt25dletJRYEJ98Fds3hoSGeyMtIwICsjjYeGdNbBVJEIi8TxrqVLl/LCCy+wbNkypk6dyuLFiwEYMmQIixcvZsWKFXTo0IGJEydy9tlnc+mllzJ27FiWL1/OaaedVuV6UlFgZsuAF/AKc5HoisTxrnfeeYfLLruM+vXrA3DppZcCsHLlSu655x527tzJ3r17GTBgQJXvD3W9RBaocK8tzZMXqbkWGWkUVBHktT3eVdVUv+uuu45p06bRpUsXJk2axPz586t8b6jrJbLADMvUVqTGDUXiXSSOd/Xu3ZtXXnmFwsJC9uzZw2uvvQbAnj17OPnkkykqKmLy5MmH1m/YsCF79uw59Pxo68lhCRPumicvcnwicbyrW7duXH755eTk5DB06FDOOeccAH7729/Ss2dP+vfvT/v27Q+tP3z4cMaOHUvXrl1Zt27dUdeTw8w5d+wVzE4B/g40B0qACc65P1Zax4A/AhcB+4DrnHMfHutzc3NzXTQv1pE9egZVfVMDPhtzcdTqEIkFq1atokOHDn6XIdWo6udkZkudc7nVvTeUPfeDwM+ccx2AXsAtZtax0joDgbaltxHAk6EUHk2aJy8iiaTacHfObS7bC3fO7QFWAZX/HhsE/N15FgEZZnZy2KutBc2TF5FEUqMxdzNrDXQF/l3ppSxgY7nn+Rz5CyA8tqyCaSNh5RTYtyPkt2mevIgkkpCnQppZOjAFuN05t7vyy1W85YghbjMbgTdsQ6tWrWpQZjk71sPqGbB8srfZrDOhTV9o0w9adIPko38lzZMXkUQRUribWSpesE92zk2tYpV84JRyz1sCmyqv5JybAEwA74BqjasFaH8x/GI9FHwIa+fAurmwYCy8/Xuo1wi+1ccL+tP6QiMFuYgkpmrDvXQmzERglXPuD0dZbTpwq5m9APQEdjnnNoevzEqSkuGU7t7tvLu84Zn1872gXzsPPnnVWy+zQ+lefV9odTak1otYSSIisSSUMfc84AfA+Wa2vPR2kZndZGY3la4zE1gPrAWeAkZGptyjqN8YOg2BQePgjk/g5vfhggegYTP4YAL84zL4fWt4digsehK2r4tqeSISXfPnz+e9996r9eds2LCBTp06Vbvegw8+WOH52WefXett11a1e+7OuXepeky9/DoOuCVcRdWKGTTr6N3Ovg0O7IPPF3pDOGvnwqzR3q33L6DPXZCUMOdxiSSM+fPnk56eHrWQffDBB7n77rsPPQ/HL5baiv9kq1Mf2vaHgb+H25bAT/4DOVfDgofh+cuh8Cu/KxRJSM8++yw9evQgJyeHG2+8keLiYj7//HPatm3Ltm3bKCkp4ZxzzuGNN94AYPDgwZx55pmcfvrpTJgw4dDnzJo1i27dutGlSxf69u3Lhg0bGD9+PI8++ig5OTm88847Fbb79ttvk5OTQ05ODl27dmXPnj045xg1ahSdOnWic+fOvPjii0fUO2nSJG699dZDzy+55BLmz5/P6NGjKSwsJCcnh6uuugqA9PR0gKN+7vz58+nTpw/Dhg2jffv2XHXVVVR3QmlNJVTjMABOPBUGPQFZXeFfo2HCeTB8MjQ73e/KRPzxr9HwxUfh/czmnWHgmKO+vGrVKl588UUWLlxIamoqI0eOZPLkyVxzzTXceeed3HTTTfTs2ZOOHTtywQUXAPDMM8/QuHFjCgsL6d69O0OHDqWkpIQbbriBBQsWkJ2dzY4dO2jcuDE33XQT6enp/PznPz9i24888gjjxo0jLy+PvXv3Uq9ePaZOncry5ctZsWIF27Zto3v37vTu3TukrzpmzBieeOIJli9ffsRrx/rcZcuW8fHHH9OiRQvy8vJYuHAh3/nOd0LaZijif8+9KmbQ/Xq4bgYUFcLT/eCjl/2uSiRhzJ07l6VLl9K9e3dycnKYO3cu69evB+D6669nz549jB8/vsKl9x5//HG6dOlCr1692LhxI59++imLFi2id+/eZGdnA9C4ceNqt52Xl8cdd9zB448/zs6dO0lJSeHdd9/liiuuIDk5mWbNmnHuuece6jFfG8f63B49etCyZUuSkpLIyclhw4YNtd5eeYm3515eq55w49vw0rUw5UewaRn0+80x58qLxJ1j7GFHinOOa6+9loceeuiI1/bt20d+fj4Ae/fupWHDhsyfP585c+bw/vvvU79+ffr06cP+/ftxzlXZOvhYRo8ezcUXX8zMmTPp1asXc+bMCWlIJCUlhZKSkkPP9+/fX+17jvW5devWPfQ4OTmZgwcPVvt5NZGYe+7lNWwO174GPUbA+0/APwbD19uqXNXv60iKxIu+ffvy8ssvs2XLFgB27NjB559/DsCdd97JVVddxf33388NN9wAwK5duzjxxBOpX78+q1evZtGiRQCcddZZvP3223z22WeHPgeObBFc3rp16+jcuTN33nknubm5rF69mt69e/Piiy9SXFzM1q1bWbBgAT169KjwvtatW7N8+XJKSkrYuHEjH3zwwaHXUlNTKSoqOmJboXxupGgXFSClDlw0Flp0hdd/Cn85Fy7/B2R1O7RKWT/4srbBZf3gAZ31KlJDHTt25IEHHuCCCy6gpKSE1NRUxo0bx4YNG1i8eDELFy4kOTmZKVOm8Ne//pUrr7yS8ePHc8YZZ9CuXTt69eoFQGZmJhMmTGDIkCGUlJTQtGlT3nzzTb773e8ybNgwXn31Vf70pz8daikM8Nhjj/HWW2+RnJxMx44dGThwIHXq1OH999+nS5cumBkPP/wwzZs3rzBUkpeXR3Z2Np07d6ZTp05063Y4H0aMGMEZZ5xBt27dKvSXv+yyy6r83NWrV0f837jalr+REu2WvyHbtBxevBr2boFL/gBdrwYgb8y8Kq9Gk5WRxsLR50e7SpFaUcvfYIh0y9/E0iIHRrwNrXrBq7fA63fAwQMRuY6kiEikKNyr0uAkuHoq5P0ElkyESRdzRqOqQ1z94EUkFincjyY5BfrfD9+bBF9+zIvcSV7qpxVWUT94CTK/hmQlNLX9+Sjcq3P6ZXD9HOo1OIF/pPyW29LfwnDH1Q9es20kVtSrV4/t27cr4GOUc47t27dTr97xNzvUAdVQFe6EV26E/86CLlfCd//ozbIJUeXZNuDt+euCIeKHoqIi8vPzQ5qrLf6oV68eLVu2JDU1tcLyUA+oaipkqNIyYPjzXt/4t8d4bYcv/ZN3tmsIxs5eUyHYAQqLihk7e43CXaIuNTX10FmdEp8U7jWRlOT1j3fF3gVCmnwb8n4c0ls120ZEoklj7sejz93eWPybv4ZVr4f0lqPNqtFsGxGJBIX78UhKgsFPemewTr0BNq+o9i2jBrQjLTW5wjLNthGRSFG4H6/UNG8MPq0xPDccdh9xydgKBnfN4qEhncnKSMPguGbbiIiESrNlauuLlfDMADjpNPjhv6BOA78rEpE4pvYD0dK8Ewx7xrvYwSs3QrmWoCIiflG4h8O3B8AFv4NVr8G8+/2uRkREUyHDptfNsP1TePdROKnNoW6SIiJ+0J57uJjBwIfhW+fBa7fDhnf9rkhEEpjCPZySU71GY42zvZ7w29f5XZGIJCiFe7ilZcCVLwIGz30fCr/yuyIRSUAK90ho/C0Y/hzs/B+8dA0UH3ltRRGRSFK4R8qpZ8F3H4fPFsCMO0CtVUUkijRbJpJyroDta+GdR7wmY2ff5ndFIpIgFO6Rdt4vvSmSb/wKGp8G7S/yuyIRSQAalom0pCQYPB5adIUp18Pm//hdkYgkAIV7NNSpD1c8D2knwnOXw+7NflckInFO4R4tDZvDlS/A/l3w/HAo0kU6RCRyFO7R1LwzDJsIm5fDG/f4XY2IxLFqw93MnjGzLWa28iiv9zGzXWa2vPT26/CXGUfaDYSzboXFT8OaWX5XIyJxKpQ990nAhdWs845zLqf0praI1en7a2jWGV4dCXu+9LsaEYlD1Ya7c24BsCMKtSSOlLow9Gk48DVMu1k94EUk7MI15n6Wma0ws3+Z2elh+sz41rQ9DPgdrJsLH/zF72pEJM6EI9w/BE51znUB/gRMO9qKZjbCzJaY2ZKtW7eGYdMBl/sj+PZAePPX3uX6RETCpNbh7pzb7ZzbW/p4JpBqZk2Osu4E51yucy43MzOztpsOPjMY9IQ3/33K9ZoeKSJhU+twN7PmZmalj3uUfub22n5uwmjQBAb/Gbaugjfv9bsaEYkT1faWMbPngT5AEzPLB+4FUgGcc+OBYcDNZnYQKASGO6cWiDXSph/0GgmL/uw9/vYFflckIgFnfuVwbm6uW7JkiS/bjklF++HpvrD3S7j5PUhv6ndFIhKDzGypcy63uvV0hmqsSK3nTY/8Zg+8eov6v4tIrSjcY0nTDnDBA/DpG/DBU35XIyIBpnCPNd2vh7YDvN4zX37idzUiElAK91hjBoPGQb0TSqdH7ve7IhEJIIV7LErPhMFPwpaPYc59flcjIgGkcI9VbftDz5vg30/Cp3P8rkZEAkbhHsv6/QaadvSai+1VuwYRCZ3CPZaVTY/cvwum36rpkSISMoV7rGt2OvS/H/47y7vAh4hICBTuQdDzRq8twRv3wJbVflcjIgGgcA8CM2/2TJ10b3rkwW/8rkhEYpzCPSjSm3rdI7/8CObqSoYicmwK9yD59gDofgO8/wSse8vvakQkhincg+aC30KTdjBtJBR+5Xc1IhKjFO5Bk5oGQ/4CX2+BmaP8rkZEYpTCPYhadIXev4CP/gkrp/pdjYjEIIV7UJ3zM8g6E17/Keze7Hc1IhJjFO5BlZwCl03wpkXq4h4iUonCPciatPEOsK6bC0sm+l2NiMQQhXvQdb8eTjsf3vgVbF/ndzUiEiMU7kFXdnGP5Drwyo1QfNDvikQkBijc48EJLeDi/wf5i2Hho35XIyIxIMXvAiR005YVMHb2GjbtLKRFRhqjBrRjcNcs78XOw2DNTJg/Btr0hxY5/hYrIr7SnntATFtWwF1TP6JgZyEOKNhZyF1TP2LasoLDK130CDTI9IZndO1VkYSmcA+IsbPXUFhUXGFZYVExY2evObygfmNv/H3rajUXE0lwCveA2LSzMLTlbfp6zcUWjYPPFkShMhGJRQr3gGiRkRb68v73w0lt4JWbvUv0iUjCUbgHxKgB7UhLTa6wLC01mVED2h25cp363tmrezbDv+6MUoUiEksU7gExuGsWDw3pTFZGGgZkZaTx0JDOh2fLVNbyTOj9c1jxPHwyPaq1ioj/zPnUkyQ3N9ctWbLEl20njOIimNgfvvocRi6Chs38rkhEasnMljrncqtbT3vu8Sw5FS77CxTtg+m3qbmYSAJRuMe7zHbQ7z74dDZ8+De/qxGRKKk23M3sGTPbYmYrj/K6mdnjZrbWzP5jZt3CX6bUSo8bIftcmHU37FjvdzUiEgWh7LlPAi48xusDgbaltxHAk7UvS8IqKQkG/xmSUrzpkSXF1b9HRAKt2nB3zi0AdhxjlUHA351nEZBhZieHq0AJk0Yt4aKxsHERLPyj39WISISFY8w9C9hY7nl+6TKJNWd8HzoOgrcehE3L/a5GRCIoHOFuVSyrclqGmY0wsyVmtmTr1q1h2LTUiBlc8pjXXOzlH8L+3X5XJCIREo5wzwdOKfe8JbCpqhWdcxOcc7nOudzMzMwwbFpqrH5jGDbRm/v++u2aHikSp8IR7tOBa0pnzfQCdjnnNofhcyVSTj0bzrsbVk6BpZNq9NZpywrIGzOP7NEzyBszr2LLYRGJGdVerMPMngf6AE3MLB+4F0gFcM6NB2YCFwFrgX3ADyNVrITRd+6AzxfCrNHQsjs071TtW8p6ype1Hi7rKQ8cvQ2CiPhC7QcS2d6tMD4P6p4AI+ZD3fRjrp43Zh4FVbQezspIY+Ho8yNTo4hUoPYDUr30TBj6NOxYBzN+Vu34e8g95UXEdwr3RJfdG869E/7zAiyffMxVa9RTXkR8pQtkJ5gqL7Lde5Q3/j7j55CVC03bV/neUQPaVRhzh2P0lK/J9jVeLxJ22nNPIEe9yPaKL2DIU96Y+z+vgwP7qnx/jXvKh7p9zbgRCTsdUE0g1R4QXTcP/jEEul7lXWg72tsXkWrpgKocodoDoqedD+f8DJY9CytejP72RSRsFO4JJKQDon3uglZnw+s/hW2fRn/7IhIWCvcEEtJFtpNTvPYEKXW98fei8O1V1+gi3yJSKwr3BBLyAdETWsCQCfDlSph1V/S3LyK1pgOqcnRv/trr/T7sGeg01O9qRAQdUJVwOP9XcEpPmP4T2L7O72pEpAYU7nJ0yakwdCIkJXv93w9+43dFIhIihbscW8YpMPhJ2LwC3rjH72pEJEQKd6le+4ug1y3wwQT4ZLrf1YhICBTuEpp+90HWmfDqrfDVBp+LEZHqKNwlNCl1vFkzAP/8IRw84G89InJMCncJ3YmtYdATsOlDeO3Huv6qSAxTuEvNdLwUzvslrHge3nrQ72pE5CjUz11qrvco2Pk/WPAwNGoJZ17rd0UiUonCXWrODC55FHZv8hqMndAC2vb3uyoRKUfDMnJ8klPh+3+DZqfDS9fCpuV+VyQi5Sjc5fjVbQhXvgT1G8Nz3/eGakQkJijcpXZOOBmu+icU7Ydnh0HhV35XJCIo3CUcmnaA4ZPhq8/ghavVg0YkBijcJTyyz/F60Hz+Lky7GUpK/K5IJKFptoyET+dhsGsjzLkPGp0C/X/jd0UiCUvhLuGVdzvs3AgLH/M6Sna/3u+KRBKSwl3CywwGPgy7C2DmKDghC9oN9LsqkYSjMXcJv+QUr8nYyV3g5f+DgqV+VySScBTuEhl1Gnhz4BtkwnOXw47P/K5IJKEo3CVy0pvC1VOg5CBMHgb7dvhdkUjCULhLZDVpC8Of9w6yPj8cigr9rkgkIYQU7mZ2oZmtMbO1Zja6itf7mNkuM1teevt1+EuVwDr1LBjyF9j4Abxyo+bAi0RBtbNlzCwZGAf0B/KBxWY23Tn3SaVV33HOXRKBGiUenH4Z7CqAN37pXWj7QvWCF4mkUKZC9gDWOufWA5jZC8AgoHK4ixzbWbd4JzktGuddtq/vvd7USREJu1DCPQvYWO55PtCzivXOMrMVwCbg5865j8NQn8QTMxjwIBzcD+8+Cnu+gEv/5LUPFpGwCiXcq9q1qnzxzA+BU51ze83sImAa0PaIDzIbAYwAaNWqVQ1LlbiQlAyXPOad3PTW72DvFq8vfN2GflcmEldCOaCaD5xS7nlLvL3zQ5xzu51ze0sfzwRSzaxJ5Q9yzk1wzuU653IzMzNrUbYEmhmc+wv47uOwfj5MusQLeREJm1DCfTHQ1syyzawOMByYXn4FM2tu5g2emlmP0s/dHu5iJc6ceS0Mfw62roGJ/WH7Or8rEokb1Ya7c+4gcCswG1gFvOSc+9jMbjKzm0pXGwasLB1zfxwY7pyrPHQjcqR2F8J1r8P+3V7A50e2VcG0ZQXkjZlH9ugZ5I2Zx7RlBRHdnohfzK8Mzs3NdUuWLPFl2xKDtq2FZ4fA11vhe3+Db18Q9k1MW1bAXVM/orCo+NCytNRkHhrSmcFds8K+PZFIMLOlzrnc6tbTGaoSG5q0gR+9CSe18c5kXfZs2DcxdvaaCsEOUFhUzNjZa8K+LRG/KdwldjRsBj+cCdm94dVb4O2xEMa/LDftrLr1wdGWiwSZwl1iS92GXjfJMy6Htx6AGT+DkuLq3xeCFhlpNVouEmQKd4k9KXVg8Hjvqk5LJsJL14Sl4dioAe1IS02usCwtNZlRA9rV+rNFYo3CXWJTUpJ3DdaBD8PqGfD3QbVuGTy4axYPDelMVkYaBmRlpOlgqsQtzZaR2PfxKzB1BJzY2usPn6GzmyVxabaMxI/TL4MfvAJ7voSJF8AXK/2uKGHpPIHgULhLMLT+DvzfLMDgrwPho5fDOpNGqld2nkDBzkIcULCzkLumfqSAj1EKdwmOZh3h+tK58FN+BC9cCbs3+11VwtB5AsESSldIkdjRqKV3stOiP3tdJcf1hAEPQNcfRKU3/LRlBYydvYZNOwtpkZHGqAHtAnVAtjb16zyBYNGeuwRPcgrk/Rhufg+ad4Lpt3mzab7aENHNBn1Yorb16zyBYFG4S3CddBpc+zpc/Aco+BD+fBYsejJsJz1VFvRhidrWr/MEgkXhLsGWlATdfwS3LIJT82DWaHjmQq+NcJgFfViitvXrPIFg0Zi7BMpRx4wbtYSr/gn/eQlm3Qnjv+NdECTv9rBdxq9FRhoFVQRhUIYlwlH/4K5ZCvOA0J67BEa1Y8Zm0OVyuGUxtL8Y5j0AT50Hm1eEZftBH5YIev1SMwp3CYyQx4zTM+F7k+Dyyd7l+yacB3N+A0X7a7X9oA9LBL1+qRm1H5DAyB4944grs4N3BffPxlxc9ZsKv4LZ98DyZ+GktjDoCWjVK5JlShyLhamwaj8gcee4puKlnQiDx3ntCw5+4x1snfmLWjchk8QTtKmwCncJjFqNGZ92Pox8H3qMgA8mwB86wmu3w9b/RqhaiTdBmwqrcJfAqPWYcd10uOhhL+Q7D4Plz8G47jD5e7B+vnrVyDEFbSqspkJKoIRlKl7TDt7Ye997vYuBLH7aO8O1WSfoNdIL/pS64SlY4kbQpsJqz10SV3om9BkNt6+ES58AVwKvjoRHO8HbD8PX2/yuUGJI0KaSaraMSBnnYP1b8P6fYe2bkFLPu5Zrr5HQtL3f1UkMCNJsGYW7SFW2rvE6T654AQ7uhzb9vJA/7fyodJ88mlgIF/GXwl0kHL7eBkv+6s2w+XoLZHaAs0ZCh0shLSOqpZRNxSs/YyMtNVknIiUYhbtIOB38BlZOgffHwZcrwZLhlB7Qpq+3V9+8i9fELILyxsyr8oBeVkYaC0efH9FtS+wINdw1W0YkFCl1IedK6HIF5C+G/86GtXO8/jXzHoAGmXBaadCfdj40OCnsJQRtKp74S+EuUhNm3h77KT2g76+83jXr5nlB/+kb8J8XAIOsbl7Qt+kHWWdCUnK1H12doE3FE39pWEYkXEqKYdNyL+jXzoGCJd70ynoZ3t58m37eME7D5sf18RpzD75wHBDXmLuI3/bt8KZWrp3rhf3eL73lzTpBi67eyVSZ7b37hieHNAtHs2WCK1y/nBXuIrHEOfjiIy/k18+HLz+GfeVOkqrbyJtLn9kemnYsfdwB0pv6OvUy3vj5yzFcB8R1QFUklpjByWd4t3Pu8JZ9vQ22rIKtq2HLJ7BlNayaDh/+7fD70k70Qr4s7MvuGzRR6NdQ5T3nsq6OQFQCPtoHxEMKdzO7EPgjkAw87ZwbU+l1K339ImAfcJ1z7sMw1yoSXxo0gexzvFsZ57yDtFtXeWFfdr9yCuzfdXi95LqQ3szbs09vBg2bVXyeXu65+uQAx+7qGI1wj/YB8WrD3cySgXFAfyAfWGxm051zn5RbbSDQtvTWE3iy9F5EasLMC+qGzeBbfQ4vdw72bPb29Lf9F3Zv8n4J7P0SvtoAGxfBvu1Vf2a9jIq/ABo0hXqNvC6ZddKhTgOo29B7fGhZ6ePU+nHzF4LfU0lHDWhX5Zh7pHrThLLn3gNY65xbD2BmLwCDgPLhPgj4u/MG8BeZWYaZneyc2xz2ikUSkRmc0MK7telb9TrFRfD1Vi/w926BPV8c/gVQtix/sXdftC/UDVcM/bL7lLreXw/JqZBcB1LqePfJdUqX1T38OKVupeWp3tRQS650nxT6ckvy/k3MvBoP3SdVXGZJhx53bPQNm3Z9g8NweL+wHNCiURoU7jz871z62uFfanaM16i4TuWfWTmDO2diJe34wxv/JX/XAZpnNIjomH8o4Z4FbCz3PJ8j98qrWicLULiLREty6uFfANUpPggH9nq3b/bCga/hwJ7Sx+WXl93v8dYpW7bva++XSfE3UHyg9PEBOHig9PkBcMXV1xFFMwDqVfHCN8Dvo1PDoNIbfW+H/r+J6LZCCfeq/iarPMUmlHUwsxHACIBWrVqFsGkRiYjkFK83TiT745QUHw768qFffMB7zRWXuy+p9PxYy0tKbw5wlR6XPi97XOF1WLHxK+at+oLdhUU0Skvh/PbNOCPrhNKCHYcv2FJ6X/555dfKVDnjsJp1Wnav2b/lcQgl3POBU8o9bwlsOo51cM5NACaANxWyRpWKSLAkJUNSGqQePmDoTUX8zLd5+l16QJeobc1foXQ6WgxnmCBIAAAEyElEQVS0NbNsM6sDDAemV1pnOnCNeXoBuzTeLiLlBe0C00FXbbg75w4CtwKzgVXAS865j83sJjO7qXS1mcB6YC3wFDAyQvWKSEAF7QLTQRfSPHfn3Ey8AC+/bHy5xw64JbyliUg88XsqYqLRNVRFJCqOdrKOulpGhsJdJIqmLSsgb8w8skfPIG/MvIQabw7aBaaDTr1lRKLE794mfiv7jupqGR0Kd5Eo8bu3SSwY3DUrYb6r3zQsIxIlOqAo0aRwF4kSHVCUaFK4i0SJDihKNGnMXSRKdEBRoknhLhJFOqAo0aJhGRGROKRwFxGJQwp3EZE4pHAXEYlDCncRkTik2TIiEhjelZw0lTQUCncRCYREb7xWUxqWEZFA0JWcakbhLiKBoMZrNaNwF5FAUOO1mlG4i0ggqPFazeiAqogEghqv1YzCXUQCQ43XQqdhGRGROKRwFxGJQwp3EZE4pHAXEYlDCncRkThkzjl/Nmy2Ffj8ON/eBNgWxnKCQN85Meg7J4bafOdTnXOZ1a3kW7jXhpktcc7l+l1HNOk7JwZ958QQje+sYRkRkTikcBcRiUNBDfcJfhfgA33nxKDvnBgi/p0DOeYuIiLHFtQ9dxEROYbAhbuZXWhma8xsrZmN9rueSDOzU8zsLTNbZWYfm9lP/K4pGsws2cyWmdnrftcSLWaWYWYvm9nq0p/3WX7XFElm9tPS/6ZXmtnzZlbP75oiwcyeMbMtZray3LLGZvammX1aen9iuLcbqHA3s2RgHDAQ6AhcYWYd/a0q4g4CP3POdQB6AbckwHcG+Amwyu8iouyPwCznXHugC3H8/c0sC/gxkOuc6wQkA8P9rSpiJgEXVlo2GpjrnGsLzC19HlaBCnegB7DWObfeOXcAeAEY5HNNEeWc2+yc+7D08R68/+HjuuepmbUELgae9ruWaDGzE4DewEQA59wB59xOf6uKuBQgzcxSgPrAJp/riQjn3AJgR6XFg4C/lT7+GzA43NsNWrhnARvLPc8nzoOuPDNrDXQF/u1vJRH3GPALoMTvQqLoW8BW4K+lw1FPm1kDv4uKFOdcAfAI8D9gM7DLOfeGv1VFVTPn3GbwduCApuHeQNDC3apYlhDTfcwsHZgC3O6c2+13PZFiZpcAW5xzS/2uJcpSgG7Ak865rsDXROBP9VhROsY8CMgGWgANzOxqf6uKL0EL93zglHLPWxKnf8qVZ2apeME+2Tk31e96IiwPuNTMNuANu51vZs/6W1JU5AP5zrmyv8pexgv7eNUP+Mw5t9U5VwRMBc72uaZo+tLMTgYovd8S7g0ELdwXA23NLNvM6uAdgJnuc00RZWaGNw67yjn3B7/riTTn3F3OuZbOudZ4P995zrm436Nzzn0BbDSzsqs99wU+8bGkSPsf0MvM6pf+N96XOD6AXIXpwLWlj68FXg33BgJ1DVXn3EEzuxWYjXd0/Rnn3Mc+lxVpecAPgI/MbHnpsrudczN9rEki4zZgcumOy3rghz7XEzHOuX+b2cvAh3gzwpYRp2eqmtnzQB+giZnlA/cCY4CXzOxHeL/ovhf27eoMVRGR+BO0YRkREQmBwl1EJA4p3EVE4pDCXUQkDincRUTikMJdRCQOKdxFROKQwl1EJA79f1Z3DVp+jN1yAAAAAElFTkSuQmCC\n", "text/plain": [ "