Created
January 11, 2019 02:41
-
-
Save petigura/956b6295e8f6e993fa05f466105fd323 to your computer and use it in GitHub Desktop.
Correcting Errors in the Labels a Posteriori.ipynb
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Simulate data and fit it\n", | |
| "\n", | |
| "- add errors both in x and y\n", | |
| "- RMS between true and observed x coord is equal to what we put in " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 75, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Populating the interactive namespace from numpy and matplotlib\n", | |
| "RMS(xtrue - xobs) = 0.0947395365753\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.legend.Legend at 0x1818d0cb10>" | |
| ] | |
| }, | |
| "execution_count": 75, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlclNX+wPHPYRtAEAUVTYVBU0szKS1puYlZVyusbLVI\nzVBKU7PS7Fa3uv6iLC2zm2ZoaiZpi6nZdsvc2tTULHfTBMRdQJR9mfP7Y4BAYRhghplhvu/Xi1fN\nPGeecx5hvnPmPOd8j9JaI4QQovHzcHQDhBBCNAwJ+EII4SYk4AshhJuQgC+EEG5CAr4QQrgJCfhC\nCOEmJOALIYSbkIAvhBBuQgK+EEK4CS9HVdyiRQttNBodVb0QQrikLVu2nNJat6zLax0W8I1GI5s3\nb3ZU9UII4ZKUUil1fa0M6QghhJuQgC+EEG5CAr4QQrgJh43hV6WoqIi0tDTy8/Md3RRRC76+vrRr\n1w5vb29HN0UIYYFTBfy0tDQCAwMxGo0opRzdHGEFrTXp6emkpaURERHh6OYIISxwqiGd/Px8QkJC\nJNi7EKUUISEh8q1MCBfgVAEfkGDvguR3JoRrcLqAL4QQwj4k4FeQnp5OZGQkkZGRtG7dmrZt25Y/\nLiwstFk9q1atIigoqPzc/fv3t9m5AbZu3co333xT/njZsmVMnTrVpnUI0RhFR0cTHR3t6GbYjVPd\ntHW0kJAQtm3bBsCLL75IQEAAEyZMqFRGa43WGg+P+n1W9u3bl+XLl9frHNXZunUrO3bsYMCAAQAM\nGjTILvUIIVyL9PCtsH//frp27UpsbCzdunXj0KFDNGvWrPz4kiVLGDFiBADHjx/njjvuoFevXlx5\n5ZVs2LDB6noeeOCBSh8CAQEBgPkbQb9+/bjjjjvo0qULQ4cOLS+zceNGrrrqKnr06EHv3r3Jyclh\n8uTJJCUlERkZyaeffsrcuXMZP348AAcPHqRv375ceuml3HjjjaSlpZXX/dhjj3H11VfToUMHli1b\nVvd/MCGEU3LeHv6S2ZB6wLbnDOsIgx+p00v37NnDwoUL6dWrF8XFxdWWGzduHE899RRRUVEkJycT\nExPDjh07ziu3Zs0aIiMjARg8eDBPP/20xfq3bt3Kzp07CQ0NJSoqig0bNhAZGcngwYNZunQpl19+\nOVlZWfj6+vL888+zY8cO3nzzTQDmzp1bfp7Ro0czYsQIYmNjSUxMZPz48Xz66acAnDhxgp9++ont\n27dzzz33yDcD4VaSkpLYsGEDBQUFGI1GEhISiI2NdXSzbMp5A76T6dixI7169aqx3KpVq9i7d2/5\n48zMTPLy8vDz86tUrrZDOlFRUVxwwQUAREZGkpycjMFgICwsjMsvvxyAoKCgGs+zceNGvvjiCwCG\nDh3Kv//97/Jjt99+O0opLr30Ug4fPmx124RwdUlJScTHx1NQUABASkoK8fHxAI0q6DtvwK9jT9xe\nmjRpUv7/Hh4eaK3LH1ecg661ZtOmTfj4+NS6Di8vL0wmEwAlJSWVvkkYDIby//f09LT4LaOuKtZR\n8fqEcCb2uKla1rOvKDc3l7i4OObMmWPz+tauXWvzc1qjxjF8pVR7pdQapdQupdROpdRjVZRRSqm3\nlFL7lVJ/KKUut09znYOHhwfNmzfnzz//xGQyVRrvvuGGG5g5c2b547KbwNYwGo1s2bIFMM+sKSkp\nsVi+a9eupKamsnXrVgDOnDlDSUkJgYGBnD17tsrXREVF8fHHHwOwaNEirrvuOqvbJ0RjdW6wr+l5\nV2VND78YeFJrvVUpFQhsUUp9p7XeVaHMTUCn0p/ewDul/220Xn31Vfr370+rVq3o2bNn+R/GzJkz\nGTVqFPPnz6e4uJi+fftW+gCw5OGHH+a2227jiy++ICYmplKPuyoGg4HFixczatQo8vPz8fPzY/Xq\n1Vx//fVMnTqVyy67jGeffbbSa2bOnMlDDz3EK6+8QmhoKPPnz6/bP4AQDmKP3rHRaCQl5fw08+Hh\n4Q7rjduDqu1Xd6XUCuBtrfV3FZ57F1irtV5c+ngvEK21PlrdeXr16qXP3QBl9+7dXHzxxbVqj3AO\n8rsTTq+oEI4fhnbn53wqG8PPzc0tf87f35/ExESnG8NXSm3RWtd8Q7EKtZqWqZQyApcBG8851BY4\nVOFxWulzQgjheMfS4JXH4fWnIT/3vMNls9bKvlWHh4c7ZbCvL6tv2iqlAoClwHit9Zm6VKaUigfi\nAcLCwupyCiGEqJ1fVsGit8HLG4Y/Cb7+VRaLjY0tv0HbmIZxKrIq4CulvDEH+ySt9WdVFDkMtK/w\nuF3pc5VorROBRDAP6dS6tUIIYa38PEiaaQ74nbvDiKcg2PLe34010JepMeArcyrE94DdWus3qin2\nOTBGKbUE883aLEvj90IIYVep++HdV+DEUbj1AYi5Dzw8Hd0qh7Omh38NMATYrpQqm2P4DBAGoLWe\nDXwF3AzsB3KB4bZvqhBC1EBr+H4FfPoeBAbBhCnQ5VJHt8pp1BjwtdY/AhYTnmvzVJ9HbdUoIYSo\ntewzMP91+H0jXNobhj9hDvqinCRPO4enp2d52uKyFAabN29m3LhxgHmM7+effy4vv3z5cnbt2lXd\n6apVlhjNmvpt5fTp08yaNav88ZEjR7jrrrtsdn4hHGbfdvjPaNi51bxKf+yLEuyr4LypFRzEz8/v\nvNWxRqOxPI/O2rVrCQgI4OqrrwbMAT8mJoauXbvarX5bKQv4o0ePBuCCCy4oT5wmhEsylcCXS+Dz\nJGjZGv71BoR3cnSrnJb08K2wdu1aYmJiSE5OZvbs2UyfPp3IyEjWrVvH559/zsSJE4mMjOTAgQMc\nOHCAAQMG0LNnT/7xj3+wZ88ewJyW+KqrrqJ79+4899xztap/wYIFjBkzpvxxTExM+WyCgIAAnn32\nWXr06EFUVBTHjx8HzGmaBw0aRI8ePejRowc///wzTz/9NAcOHCAyMpKJEyeSnJzMJZdcApjzAQ0f\nPpzu3btz2WWXsWbNmvK677jjDgYMGECnTp146qmn6vvPKYRtZJ6CaU/Dig/gyj7w/NsS7GvgtD38\ntz/J5ECa7XaZAujYzocxdze3WCYvL688bXFERESlPDlGo5FHHnmk0sYot956KzExMeVDI/369WP2\n7Nl06tSJjRs3Mnr0aFavXs1jjz3GqFGjGDp0qMVUC5bqr0pOTg5RUVEkJCTw1FNPMWfOHJ577jnG\njRtHnz59ynPyZGdnM2XKFHbs2FH+DaLicNHMmTNRSrF9+3b27NnDP//5T/bt2weY8wH99ttvGAwG\nunTpwtixY2nfvn1VzRGiYfy+0TxeX1hgnlt/zY2ObpFLcNqA7yj1GVLJzs7m559/5u677y5/rizH\nzk8//cTSpUsBGDJkCJMmTbJJ/T4+PsTExADQs2dPvvvOnPFi9erVLFy4EDDfFwgKCiIzM7Pa8/z4\n44+MHTsWgIsuuojw8PDygN+vX7/y1Mtdu3YlJSVFAr5wjKJCWDofVi2D9h3g4X9Ba+v+Fv88VMiO\nAwUMig60cyOdl9MG/Jp64s7IZDLRrFmzagO2eUlD7VVMmwyV0zF7e3uXn7ch0ibbqw4hanT8sHlu\nfep+uP5WuHsEeNechrywSPP+l1l8tOoMwU09GRDVBD9f9xzNds+rrodzUw9XfNy0aVMiIiL45JNP\nAHNO+d9//x2Aa665hiVLlgDmRE21YTQa2bZtGyaTiUOHDrFp06YaX9OvXz/eeecdwJxbPysry2La\n5H/84x/l7dq3bx+pqal06dKlVu0Uwm42rIbJY+DUMXj0ebh/tFXBfvv+fEYkHGXxt2fo37sJ7z3X\nxm2DPUjAr7WBAweybNkyIiMj+eGHHxg8eHB5KuIDBw6QlJTEe++9R48ePejWrRsrVqwAYMaMGcyc\nOZPu3bvXejepa665hoiICLp27cq4cePKd7iyZMaMGaxZs4bu3bvTs2dPdu3aRUhICNdccw2XXHIJ\nEydOrFR+9OjRmEwmunfvzr333suCBQtqTM8shN0V5MO812HuaxDWAV6YBZddXePLcvNNzPgog8fe\nOEFxiWbquFZMHBJCoL97h7xap0e2FUmP3LjI707Y3KG/4N2XzUM5t9wHA2PBs+b0CJt25vHG4gxO\nZpYwKDqQuIFBjapXX5/0yE47hi+EcFNaw5qV8PEcaBIIT06Bi3rU+LKs7BLeWXqabzfmEN7ai7ee\nDKVbB/mWWpEEfCGE88g+C+9Ph99+hu5XwENPQmAziy/RWrP+tzze+iiDMzkmhtzUlNgBQfh4122S\nRGPmdAFfa13n2SzCMWTDc2ETf+6AOa9CVibcEw833A4elodi0rNKmLEkgx9/z6NzmA+vjQ2mY7ua\nb+a6K6cK+L6+vqSnpxMSEiJB30VorUlPT8fX19fRTRGuylQCX30Mn38AIaHm9AjGzhZforXmm19y\nmLU0k6JiiL+9GXf3C8TTUxEdHQ00/tz2deFUAb9du3akpaVx8uRJRzdF1IKvry/t2rVzdDOEKzqd\nDnOnwp5tcGU0DBkLfk0svuToqWLe+DCDLXvyufRCA0/GBtM+1Lth2uvinCrge3t7ExFx/gbDQohG\naPuv8N40KMyHB58wp0ew8M2+xKRZvvYs732ehYcHjB/cnJhrA/Dw+Ps1SUlJbNiwgYKCAoxGIwkJ\nCY1uX9r6cKqAL4RwA8VF8NkC+HYptIuA+H/BBZb3uE4+WsS0RensOljIld18eeK+YFoFVw5fSUlJ\nxMfHl6czSUlJIT4+HkCCfimnmocvhGjkTh5l9+PDuNjXg2Wni3knvZhCiyHIE1rEoFrcCqZ89LFF\ncOaXKkuW9ezPZTAYiIqKsk37K3DUPQKZhy+EcH6b1sLCt2jnrXj+aCHrc0yWy/tGoC6IQ/mGobM2\noI99ACVVpwYBqgz2lp53RxLwhRD1UuOsmIJ8WPwO/Pg/6HgxgfFPMzkktNrz5ReaeP+LLD75/izN\nm3oyfnBzrulxD3CPxXYYjUZSUlLOez48PFxm7JSSgC+EsJ+0g+YMl8cOwc2D4dYHwKv6sPP7vnym\nJWVw+GQxt1zThIcHNSfAyvw3CQkJxMfHk5ubW/6cv78/CQkJ9b6MxkICvhCizqqdFaM1rPsKlsyG\nJgHweAJ0rT7pX06eicTlp1n5QzZtWngx7bFWXN6ldms7ym7MxsXFUVBQQHh4uMzSOYfctBVC1EnZ\nrJhze9SJ/32L2IKjsOVH6NYTHpoAQdXvb7Fhex7TF2eQnlXCndcH8mBMEH6Guic7a+wLr+pz01YC\nvhCNXFkAtLXqZsW0a+LLwXuvY056MR+fLqHaCOMZiGodiwq6Gp2fhj46F/L+qrJoYw3edSGzdIQQ\nDa662S+Hc/IZk1bIngILncmmvVGth4CnP/rkMvSplaBlJzV7kx6+EKJOqp0V0749yampVb7m5Oli\n3lycyS/b8+gS7sPEB4Lp0FaSndVGfXr4jWdXACFEg0pISMD/nKR5/v7+JLzyynlltdZ88WM2D00+\nytY9+TxyRzPenhgqwb6ByZCOEI1Mg9y0LC4m1lAAvTvx1OY/OZKTX+2smMMni3g9KYNt+wqI7GxO\ndta2pSQ7c4QaA75Sah4QA5zQWl9SxfEgYBEQVnq+aVrr+bZuqBDCSZw8BnOmwF97iI0byftqKZ30\n+R8wJSbNZ2vOMu/zLLw84Yn7g7nlmiaS+tyBrBnSWQAMsHD8UWCX1roHEA28rpSS72lCOEDZvPh1\n69ZhNBpJSkqybQWb18Pk0XD0EDzyDAwZW2UunINHChk77TjvLD3N5Rf5Mu/5NsRcGyDB3sFq7OFr\nrdcrpYyWigCByvybDAAyALndLkQDs2u2yIJ8+OhdWP81dLgIRj4NLVsDlXv2RcWaD/93hqRvsmji\n58FzD4XQt6e/BHonYYsx/LeBz4EjQCBwr9a6hqxIQrivhpwXn5ubS1xcHHPmzKnzeY0+ihdCvYkw\neJCUWcy8b7fx/TOtzyu3O7mAaYsyOHikiBuu8Gf0Xc1pFuhZ53qF7dki4PcHtgHXAx2B75RSP2it\nz5xbUCkVD8QDhIVZzn8thKgde2SLvKWpJ2NbeJFrggmHC9mcd35fLr/QxPyVWSxdfZaQIE9eHtWS\nqO5+da5T2I9V8/BLh3S+qOam7ZfAFK31D6WPVwNPa603WTqnzMMXwrYsZYtMTk6u3clyc2DhDPOY\nfdfLIW4CBAWfV+y3veZkZ0dPFTPwHwHE396MJn4y29ueHL3SNhXoB/yglAoFugBVr48WQtiNzbJF\n/rUHEl+BjJNw50PQ/y7wqBzEs/NMvPtZJl/+lEPbll5MH9+KHp1lI3tnZ820zMWYZ9+0UEqlAS8A\n3gBa69nA/wELlFLbAQVM0lqfsluLhRBVqne2SJMJ/rcUli+AZiEwaRp07HpesZ/+yOXNxZlknilh\n8I2BDLslCIOP9OpdgaRWEKKRqdPCq6xMmDcNdm6BntfCsPHgH1CpSObZEt7+OJM1W3LpcIE3E4cE\n0yXcYLuGC6s4ekhHCOHKdm2FuVMhLweGjIXrboYK0yi11nz/ay5vf5JJXoGJ4QODGHxjU7y9ZKql\nq5GAL0QjY3XPvrgYPv8Avv4YWreHJ16GdhGVipzIKGb64gw27syna4QPEx4IwdhG0iK4Kgn4Qrij\n9OOQOAUO7IbrboJ7HwbD3zddTSbNyh+zmbP8NCYTjL6rGYOiA/H0kF69K5OAL4S72fIjvP8maBPE\n/wuu7FN+KDo6GnxaEzlgBn/sL6DnRb48cX8wbVpIqGgM5LcohLsoLICP58DaL8DYGR7+F7RsU364\npERDyM2olnfw1+FCJg4JZkCUJDtrTCTgC+EOjqTCuy/D4WTzvPpBw8Dr77H4A2mFvPZBOh6hg9Fn\nNjP/1dsICZK0CI2NBHwhGjOt4adv4cNZ5jH68S/BJX/P6Css0nzwdRZLvj3D6ZQVbP/+WQryztLz\nx1rO4RcuQQK+EI1VXg588F/YtBYuioQRE80Lqkrt/KuAqR+kk3q8mFamr/np+6cpyDOv0rVppk3h\nNGThlRAuxppsm10Miudbe9PaSzEvo5jFmSWUpz1TBlSruyD4RijKQB+dz4bv51SZZM1gMBAVFWXT\n9oOdd+Nq5GThlRACMOc2uauZJw+HeHGqGB47XMiO/AqduibdUG0eQvm0RGd8hz7xCZjy7ZJpUzgf\n6eEL0VicPQ3vvQ47foXLroZhj0NAoPlQrol3lmbyzS85tA/1YkJsMN0v/HvevU0zbQq7kh6+EO5u\n9zaY+xrknIXYMRB9S3l6hB+25TJjSQans03c378pQ28Owse78lRLm2XaFE5NAr4QrqykBFYugi+X\nQGg78yyc9h0AyMgq4a2PM1j/Wx4XtvPm5dGt6BxW9XbT9c60KVyCBHwhXFX6CZgzBfbvgmv/CfeN\nBoMvWmu+3ZjDrE9Pk19oYsStQdxzY1O8PC0voIqNjS3fClFuqjZOEvCFcEW//Qzz3zD38EdOgt59\nATiWbk529uuufC7paGBCbDBhra1PdiaBvnGTgC+EKykqhE/mwurPIbyTORdO6AWYTJoV67OZs+I0\nAGPvac5t1wXgIcnORAUS8IVwFccOwbtT4NABuHGQeftBL29SjxUxLSmDHQcKuKKrL4/fF0zrEHlr\ni/PJX4UQzk5r+HkVfDgTvH1g3H/g0t4Ul2g+/iaL97/KwtfHg0lDg/lnb0l2JqonAV8IZ5afC4ve\nhg2roXN383h98xb8eaiQqR+ksz+tiOsu82PcPcEES7IzUQMJ+EI4q+Q/IfEVOHkMbhsCtwymsMSD\n95ef5qNVZ2gW4MGLI1tw3WX+jm6pcBES8IVwNlrDquXw6XvQtBlMfA06X8L2/flMXZRB2oliBlzV\nhFF3NifQ38PRrRUuRAK+EM7kbJZ5uuUfGyEyCh58glyvAOZ8lMGKddmEBnvy2tiW9LrYz9EtFS5I\nAr4QzmLvHzDnVcg+A/ePhr4D2bgrn+kfHuXk6RLu6BtI3MAg/HylVy/qRgK+EI5mKoGVH8IXi6FV\nGxg3maxgI7PeT+e7TbmEt/birSdD6dbB4OiWChcnAV8IR8o4aU56tm87XH0D+r7RrNsN/515lDM5\nJh64qSkPDDg/2ZkQdSEBX9hE2aYcsjS/FrZtgPmvQ3ERxE0gvWtf3lyYwU+/59E5zIfXxgbTsV3V\nyc6EqAsJ+ELUU60/7IoKYek880ycsI7okU/zzcHmzJp8hKJiiB/UjLuvD8SzhmRnQtRWjQFfKTUP\niAFOaK0vqaZMNPAm4A2c0lr3sWUjhWg0jqWZ59anHoAbbudIn6FM/ySbLXsyuPRCAxMeCKZdK+uT\nnQlRG9b08BcAbwMLqzqolGoGzAIGaK1TlVKtbNc84QqSkpLYsGEDBQUFGI1GyaNenV9WmVfNenlT\nMvoFlp3uxrxX0/HwgPGDmxNzrSQ7E/ZVY8DXWq9XShktFLkf+ExrnVpa/oRtmiZcQVJSEvHx8eV7\nn6akpBAfHw/gFkHfqg+7/DxImmkO+J27k3zrBKat1Ow6eJre3czJzloFy+iqsD9b/JV1BryVUmuB\nQGCG1rrKbwPCscrGmm2pLNhVlJubS1xcXPlmGrbkTDeFrfqwS90P774CJ45SdMsDLPG5hUWzzuJn\n8OCZB0Pod4W/JDsTDcYWAd8L6An0A/yAX5RSG7TW+84tqJSKB+IBwsLCbFC1cLRzg31NzzuCPT7o\noOYPu0FBnoxq4UVWCTyb2ZY/VxpRvmfRWRsoOPYBL/16lpesrMuZPuiE67JFwE8D0rXWOUCOUmo9\n0AM4L+BrrROBRIBevXppG9QtasEeQcNoNJKSknLe8+Hh4Y0+SFn6sHuptTfXBniyPseDBNNACi+4\nCYpPY0qdDtm/NXBLhTCzRcBfAbytlPICfIDewHQbnFe4gISEBOLj48nNzS1/zt/fn4SEBAe2qjJ7\nfPBER0djMBiqDPrhTf25NsjA733HM2d3V4pOFhNzTRMeviOMAD95awjHsWZa5mIgGmihlEoDXsA8\n/RKt9Wyt9W6l1DfAH4AJmKu13mG/JgtnUjZWHRcXR0FBAeHh4W4zSyciIoLU1NTKH3aeHjzf5yqm\n9/wvK9d70qYFTHusFZd38XVgS4UwU1o7ZmSlV69eevPmzQ6pW9ieu620LbvekSNHMjIujvyCAsKa\n+DL81gfZZXyW9CwTd14fyPCBQfj6SLIzYTtKqS1a6151ea3MBRM24S6BHipPxUzet5e3ruxEzIWd\nmXfRi6xKDsLo58mLI1tycYQkOxPORQK+ELVw3lTMo8cYczKd2YEjCW4SxLBbgri/f1O8vWSqpXA+\nMqQjGqWGnIoJYGjSiqhrboSCNJvV5U7fmoT16jOkI4OLQtRCtVMxc07YNNgLYQ8ypCMaJZv3jgvy\n4cNZtNu8mcM5Oecddod1B8L1SQ9fuI3o6Oi6DfUc+ouSyWP5ZIMm5OqX8PSqvJ+ss607EKI60sMX\nojpaw9ovOPjRSqZ6DmOPZxgxg/wYdXcI48eOdLt1B8L1ScAXbqHWKZyzz1I0fwYfbm9Gkte/CPD3\n4t/3BhPd0x+lhrAk6T1AbqwK1yIBXzR6tU7h/OdOds9exNTc20j2uoB+vfwYc08wQQGe5UUk0AtX\nJNMyhdNo8KmUBgNRUVHljz2Ae4L9KGlxJ5953oBnSSaFRxZA9u9W1yUfBMLeZKWtEBZYk8I52BPu\nbdudz5s8yFHVEs/M7yk8/hGY8huqmULYnfTwRaNnKYVzcnIy2Zs3M3vBQb4y9aZtQCETRrSlR2e/\nKs4khOPJwishLEhISMDf37/Sc/7+/iT832R+mrWc4e/58I3pCgZfBXNf6iDBXjRaMqQjGr2qUjg/\nM34Cf21rzXt5F9GhSRYvPdKCLhcGOLilQtiXBHzhFmJjY8v32H324WnM/KEJeRgYfvlp7hveHS9P\nSXYmGj8J+MJt+HgH0z5sOK/80IquPkeYEB+OsWsHRzdLiAYjAV80eiaTZuXKFHw6v8pRE4zpvJfb\nRkfj6ePt6KYJ0aAk4ItG7dCxQl5/Zz9/nAygp0cqTwwLoU3UjY5ulhAOIQFfNEolJZpPvjnFgq/O\n4mNSTAxdw4DxA1HNgh3dNCEcRgK+aHQOpBXy2tzD/HnCk2v1dh77pybktiHgIbOQhXuTgC8ajcIi\nzaKvT7P4f1k0NeXwgv9XXDf6VlTHi+xSn7tt3C5cnwR80Sjs/KuAaQtPkHJCc2PJRkZ3TyXoodHg\nL3PrhSgjAV+4tLx8E++tzGLZmjO05DRT9IdceX8fuG4CKJlbL0RFEvCFy9q8O483kjI4llHC7SVr\nGNFyC/6jnoS2RpvXde7wTa3z6wvhBCTgC5dzNtfEO0sz+eaXHNp7ZTCjaC7d/9EJ7p0KBl+711/r\n/PpCOAnJlilcyg/bcpmxJIPTZ0sYzCqGqv/hM2wMXHHdeWVtmV9/27ZtAERGRlqdX98W5IawOJfk\nwxeNXkZWCW99nMH63/K40P80Lxe+TWejH8TPgJZt7Fr38ePHOXPmDFrraoM9VJ93XwhnUWPAV0rN\nA2KAE1rrSyyUuwL4BRistf7Udk0U7kxrzbcbc5j16WnyC0oY0WQt92R+jNeAO+D2YeBV/Z+wLXrH\nZcM3Zd+ECwoKUEpR1Tfj8PBw6ZELp2ZND38B8DawsLoCSilP4FXgW9s0Swg4ll7M9MUZ/Lorn0ta\n5DAhdzphxWdg/H/gkvO/0dpji8SqevRVBXsPDw/8/Pzq3Qb5wBD2VGPA11qvV0oZayg2FlgKXGGD\nNgk3ZzJpVqzPZs6K0yg049r+yq0H5+JxcSSMmAJBDZcewdIwTVlP32AwEBERQWhoaIO1S4i6qPcY\nvlKqLTAI6IsEfJdX0+pRe68uTT1WxLSkDHYcKOAKYwmPn3qT1il/wh0PwoC7LaZHsEebLG2PaDQa\n7VavEPZgi5u2bwKTtNYmVcNCF6VUPBAPEBYWZoOqRWNRXKL56LszLPwqC18fxaTIA/xzyzRUsxCY\nNA06dnVIuxISEoiPjyc3N7f8OX9/f5l3L1yT1rrGH8AI7Kjm2EEgufQnGzgB3F7TOXv27KmFc1m0\naJE2GAwiPwENAAAWQ0lEQVQa0OHh4XrRokW1Ol5Xe1MK9MiEI7rvqBT94qzDOn3qS1rH9dd61v9p\nnX3GJnXUh72uW4i6ADZrK+J2VT/1DvjnlFsA3GXNOSXg102fPn10nz59bH7eRYsWaX9/fw2U//j7\n+5cHt5qO10V+QYlOXJap+z2aou+cdEivX7lD6ycGa/1wjNarV2ptMtnq8urNXv/uQtRWfQJ+jQuv\nlFKLgWigBXAceAHwLv12MPucsguAL7QV0zJl4VXdREdHs23bNiIjI6stU7ZIqDbK5pmfSylF06ZN\nazxeW5FX3YO6IA5laAOZ67m/8GMealZAapFm8rEi5n27ptbnFMId2HXhldb6PmtPprV+sC6NENap\nmL9lw4YNNp0ZUt0Hf9nzNR23loeXPxdEPoZHxH3owpMEpb3K5MA/6d7cgy+zSvjvqWLyHbP4W4hG\nT1bauohz87cUFBSQmprKc889Z5Obh5ZmoyQnJ9d43BqbdubxxuIMTmaWMCg6kLj2h/D78DiYAmDo\nOG65Mppb6nshQohqScC3MXss/oGqFwDl5uYSFxfHnDlz6n1+Pz8/PDw8MJlM5c+VzUYBy7NVapKV\nXcI7S0/z7cYcwlt78da45nT79X2Y+wUYO0P809DqgnpfgxDCMgn4LsLe+VvKhob27t2L1prw8PBK\nUw/L/hsXF0dBQcF5x6uitWbdb3n896MMzuSYeOCmpjzQ4ww+702Ew8nQ/04Y9CB4edvkGoQQlkm2\nTBdhiyEVa9hq4VV6VglvLsngp9/z6Bzmw8TY5nRMWQsfzgIfX4ibAN1lnZ4QtSXZMt1AfYZUaqO+\nq0a11nzzSw7vLM2ksBjiBzXj7qs88PxwOmxaCxf1gBFPQbMQm7RXCGE9Cfguoi5DKvZg6QPh6Kli\n3vgwgy178rn0QgNPxgbTPu8gJLwCp47D7UPh5nvBw7PhGiyEKCcB34XExsaW36C1Vf4WWwzhlJg0\ny9ee5b3Ps/DwgPGDmxNztT8eq5fD0vkQ1Byeeg06VZtdWwjRACTguxhnS9SVfLSIaYvS2XWwkN7d\nfHn8vmBaeWfDrP/AH5vgsqth2OMQEOjopgrh9iTgi2pZ2qi7qFiz5NszLPomCz+DB888GEK/K/xR\ne/+AOa9Czlm4fzT0HQg1JNUTQjQMCfhuzFJAt7RRd69r72Lqogz+OlxE317+jLm7Oc39gRUfwJeL\nIbQtjH8J2ndw1KUJIaog0zIbgbos9jp+/Dj79u2rtNDKw8ODzp07ExoaWv1G3f7BXDV0CxSfRh99\nH7J/o5UXPBfqw6V+Hnx9ppgZJ2ufHsHZhqqEcFYyLbORqymg2ypZmslkYu/evRw9erT6hV65maQf\nWMbh396gpOgsN7dpzswrOuKlNPGb9vNx6qkqX2cp2ZsQomFIwG8E6hJM161bV+XzWmsiIyOr7+H7\nNiG4cAWtL+nIIyFe3NHMi735JiYfL+JwcDsig9tVeV7pwQvheBLwXUBDb923du1akpKSGDFiJPn5\neeXH/P39SUycTWy/6+DdKXDoANwwiC53DifJ28fmbRRC2Fb1G4SKRi0hIQF/f/9Kz5Wt3M3KLiGl\nsD8R10zB0MScYyc8PJzExERiI0Lh/8ZC5kkY+x8Y/DBIsBfCJUjAd1LR0dF2y7wJ5kVciYmJGAwG\n4O+A3qbL7Tw4+Shrt+Yy6fFhRPXqRp8+fUjes4vYvMMwbxqEd4IXZkGP3k5xLUII68iQjhuruHL3\nk+WreHNxJu/NS+ciow8THwgm4gIfht3yPaT8CZPHwMljcNsQuGWwpEcQwgVJwHdza9as4cufcnho\n8lGKS2DUnc24o28gnh4KtIbvV8Anc6FpM5j4KnTuXqvzW5rrL4RoWBLwnVBDBcnDJ4p4/cMMtu0r\nILKzOdlZ25aluenPZsGCN+D3jRAZBQ8+AQG127vW0uItCfpCNDxZeOVkyoLkuWmQExMTbRYkS0ya\npavPMn9lFl6e8MidzZn6bEz58R6+iuda+xDkCe+eKmZpVkmd6ql2aqfBQFRUVJ3bXx2Z+incgSy8\ncgBX3coQQ1vUBSNQfh25qrsf4+9rTstmXkzFfAd/SHNPhgZ7caRI868jRewvrHuHwN67dAkhakcC\nvpOxW5BUXqgWA6HFQCjJxZQ2k5dmTkGVJjZb+9knMPc12LcdrrqB9rGjmevrX8NJLatprr8QomFJ\nwK8jWwessm8M4eHhNg+Suw8WMHVRBslHi7jhCn8evTucoIBX/y7w+0aY/zoUFcJDE+DqG+pUz7ka\napcuIYR1JOA7GVsGybwCE/NXZrF0zVlaBHny8qiWRHX3+7tAUSEsnQerlkP7jvDw09C6fb2v4dxN\nUxy9S5cQwkwCvhOoOCsnOTmZYcOGMW/evHoFya1783k9KYOjp4oZeG0A8YOa0cSvwjq744fh3Zch\n9QD0uw3uirPLill77NIlhKgbCfgOVtXUxffff5+IiAhCQ0NrHSSzc03MXpbJVz/l0LalF9PHt6JH\nZ9/KhX75Hha9DV5eMOYFiLzKRlcjhHBmEvBrwR4zc6qblXPw4EFCQ0NrV2fA5ag2w8ArCNK/4tDu\nZTy2rqj8sJ+Cr4fEwM+rzPvLjpwEwS1tdCVmVa0hkJ69EM6hxoCvlJoHxAAntNbn7UKtlIoFJgEK\nOAuM0lr/buuGNlY2mZXjGYhqPQQVFIXOT0UfehPyD1YqcqGP4vnW3ube/cBYiLkfPG2bHkEWWgnh\n3GpceKWUug7IBhZWE/CvBnZrrTOVUjcBL2qta8yqJQuvzCxNXUxOTrb4Wq01qzblMvPTTPIKTDxw\nUxCDb2yKt5eqWAjWrISP55hXyo54Ci7q0WDfVkAWWglhS3ZdeKW1Xq+UMlo4/nOFhxuAqnfAEFWy\nNCvn3NkuFZ3IKGb64gw27syna4QPEx4IwdjGu3Kh7LPm9AjbfoFLr4ThT0JgkN2uxZZrCMp28ZKd\nsoSwHVuP4ccBX9v4nI1a2VBHVVMXq1pZazJpvvgxm8TlpzGZ4NG7mnF7dGmys4r+3AFzXoWsTLgn\nHm4cBOrvMo7YVKU2LH3YCSHqxmYBXynVF3PAv9ZCmXggHiAsLMxWVbs8a6cupp0oYtqiDP7YX0DP\ni3x54v5g2rQ451doKoEvl8DnSdAiFP41HYyd7Nj6v9lqDYFk2BTCTrTWNf4ARmCHheOXAgeAztac\nT2tNz549tajeokWLtMFg0IAOCwvTj06ao/uPS9XR8bt0n9v+rU0m0/kvyjyl9WsTtY7rr3XiFK1z\nsx3a7vDwcL1o0aJav97f318D5T/+/v61Po8QjRWwWVsZZ8/9sSpbZukY/he66pu2YcBqYKiuPJ5v\nkdy0rV5VGTM9vPy4dcgbZKZ+A8Wnz/8m8McmmPc6FOZD7KNw9Y2VhnBszdJN3/qMvzfkjV8ZLhKu\nyK43bZVSi4FooIVSKg14AfAG0FrPBp4HQoBZpYm4iuvaGFfUULNdTMV5rFw4BpPJhNYaX19fIiIi\naBsaysgQL+5t7sX+AhOTjxWR+uzLwMtW12frwFefG62SYVMI+7Fmls59NRwfAYywWYtEtcGtpKSk\nUpk/9+3joWAv7m3elmWni3knvZh6ZDOuFXv1jiXDphD2Iytt66mqIFTnXr8yoFrdhSFgPwXZh2ss\nXmIy8d/f95PVrCU/5pjqVidVt9dRwVUybAphPx41FxENokk3VMeXUSH96XDJTXh4WPerOZSTX69g\n72xiY2NJTEzEYDAA5p69LXf7EsKdyRaHdmTNXPKzuSZmL83k619yaB/qxYTYYLpf6EtSUlKlufnZ\n2dmkp6ef93prVuS6IpmHL0TV6nPTVnr4DvTDtlyGTz7C/zbmcH//psx5pg3dLzRntoyNjSUqKoo+\nffqQfPAgMx4ehr9X5dw3tRnqiI6Ottu2jEII1yBj+A6QcaaEtz7KYP1veVzYzptXHm1Fp/bn56Jf\nu3Yt5JyFd14i9vhOuDeGkZ9+Q54bbCYiPXshbE8CfhVsMZxQ1WrR+++/n2835jDr09PkF5owHf+Y\nfbu+ptMz31d9kv27IHEKZKXDXXHE/vNO5qRdX+u2ycpVIQRIwLeLqtIEjxwZz5Jvz5DT5BYu6Whg\nQmwwQwd/UfUJTCXw9SewYiGEhMKkN6BDF5u1RVIWC+Ge3PqmbVVj2sePH2fv3r1orTEYDOU7T9VG\ntatFAy4g6oZHIHM1x48fq1SPh4cHPj4+tDJ4k3jlhUSHBrH00Cke33KQM8UlVdRiVtMiJ1m5KkTj\nYteVtu7k+PHj7Nu3ryw/EAUFBezbtw+gVkG/2tWi2Ucg8/sq6wHoFuTHl/26E+DlydjNB/gg+WR9\nLsdyW2TlqhBuxyUDvr1mmxw8eBCTqfKcdpPJxN69ezl69KjV51HKA63PnxuvlGLbtm2cOXOGqr5Z\nnTyby1nlxcTDRaQ0a0tkZNvaX8Q5DAZDlcFdVq4K4X5kWmYF1fV6azPs5df8Yjr0noCHl995x3x9\nfS2eLzUnn1FphaQU2W6YLSIiAn9//0rPycpVIdxUXdNs1vfHGdMjh4eHV0rLW/YTHh5e42vzC0p0\n4rJM3e/RFH3npEP63wnzqk0TXJ966qK+KYuFEM6DeqRHdskhHXupax6X7fvzmboog7QTxdx0VRMe\nubM5gf7DWf/t+8D5NzMT/vOiuZ7ColrVU1fWbrAihGjcJOBXYGm7wark5puYs+I0K9Zl0zrEk9fG\ntqTXxX8P5VQZXA8nE3vgZ7j6Yh7buI/0vPwGWUQlgV4I4dbTMqtjzcKrTTvzeOPDDE6eLmFQdCBx\nA4Pw87VwS0RrWP8VLHkX/JpA3ESiH32ixnqEEKIimZZpY5YCcFZ2CbM+zeS7TbmEt/birSdD6dbB\nYPmEudmwcAZs/gG6XQ4PTYSg5rZttBBC1EACvpW01qz7LY//fpTBmRwTQ25qSuyAIHy8a9hG8MBu\nmDMFMk/BnQ9B/7ugNPWx9OyFEA1JAr4V0rNKmLEkgx9/z6NzmA+vjQ2mY7vzk51VYjLB/z6F5e9D\n8xYw6XXocFHDNFgIIaogAd8CrTXf/JLDrKWZFBVD/KBm3H19IJ6eNfTqszLhvamwaytrzpbw+p8p\nfDFFgr0QwrEk4Ffj6Kli3vgwgy178rn0QgNPxgbTPtS75hu6O7eag31eDgwZx39eeK3B2iyEEJbI\nStsqHDpeRNxLR9mdXMD4wc15Y3wr2od6W35RcTEsnQfTn4GAQHjuLZLSMtmwYQPr1q3DaDSSlJTU\nMBcghBBVkB5+Fdq18uK+/k3p37sJrYL//ieqNq/8qWPmvPV/7YHrboJ7Hybp06WSllgI4VTcfh6+\ntYnYyjJcVkyu5uHhwa3du/DBZe0BmHayiLXZ5uMNmZYYZMaPEO5C5uGfozbZNLdt22ZVuaoyXJpM\nJn7bf5A9HYKJ2/AnKbl/B3hLaYmtqbOmPPdCCFFbjTLg14a1gXXdunVVPp+Sk8/TGV4073wxFZdS\nWerhW1On9NiFELbWKAO+PYKl0WgkJSXlvOfDw8NZVUV9ZVsLnpuILTExUcbwhRAOIbN0rJTwwvP4\n+1SeqWMpw2VsbCyJiYkYDOa0C+Hh4RLshRAOVWPAV0rNU0qdUErtqOa4Ukq9pZTar5T6Qyl1ue2b\n6WB/7SV273oSr7mYYD/rA3hsbCxRUVH06dOH5ORkCfZCCIeyZkhnAfA2sLCa4zcBnUp/egPvlP7X\n9ZlM8N1n8Nl8CArhZ59gSnx86XNllIyxCyFcTo0BX2u9XilltFDkNmBh6U4sG5RSzZRSbbTW1m8C\n64zOnIZ502DHZrj8Ghg2np23DCQyMrJWwV4+GIQQzsIWN23bAocqPE4rfc51A/6JI/DqBMg5C7Fj\nIPoWUDXkzxFCCCfXoDdtlVLxSqnNSqnNJ0+ebMiqa6dFKFx6JTw7A/rGgFLlq2wlTYIQwlXZIuAf\nBtpXeNyu9LnzaK0Ttda9tNa9WrZsaYOq7cTDE4aNh/YdgL+nWJ6bJkGCvhDCldhiSOdzYIxSagnm\nm7VZjhy/r80qW2tVtYgqNzeXuLi48s3BbUnG/YUQ9lBjwFdKLQaigRZKqTTgBcAbQGs9G/gKuBnY\nD+QCw+3VWEexlCZBCCFchdsnT7OGpVW2ycnJDd8gIYTbqk/yNFlpa4WEhAT8/f0rPWdpla0QQjgj\nCfhWkDQJQojGoFEmT7OH2NjY8hu0clNVCOGKJODXggR6IYQrkyEdIYRwExLwhRDCTUjAF0IINyEB\nXwgh3IQEfCGEcBMS8IUQwk1IwBdCCDchAV8IIdyEBHwhhHATDsuWqZQ6CZyfgtK5tABOOboR9STX\n4BzkGpxDY7iGLlrrwLq80GGpFbTWTrzllZlSanNd05A6C7kG5yDX4BwayzXU9bUypCOEEG5CAr4Q\nQrgJCfiWJTq6ATYg1+Ac5Bqcg1tfg8Nu2gohhGhY0sMXQgg34fYBXyk1QCm1Vym1Xyn1dBXHlVLq\nrdLjfyilLndEOy2x4hpiS9u+XSn1s1KqhyPaaUlN11Ch3BVKqWKl1F0N2T5rWXMdSqlopdQ2pdRO\npdS6hm5jTaz4ewpSSq1USv1eeg3DHdHO6iil5imlTiildlRz3BXe0zVdQ93e01prt/0BPIEDQAfA\nB/gd6HpOmZuBrwEFRAEbHd3uOlzD1UDz0v+/yRWvoUK51cBXwF2ObncdfxfNgF1AWOnjVo5udx2u\n4Rng1dL/bwlkAD6ObnuF9l0HXA7sqOa4U7+nrbyGOr2n3b2HfyWwX2v9l9a6EFgC3HZOmduAhdps\nA9BMKdWmoRtqQY3XoLX+WWudWfpwA9CugdtYE2t+DwBjgaXAiYZsXC1Ycx33A59prVMBtNbOdi3W\nXIMGApVSCgjAHPCLG7aZ1dNar8fcpuo4+3u6xmuo63va3QN+W+BQhcdppc/Vtowj1bZ9cZh7N86k\nxmtQSrUFBgHvNGC7asua30VnoLlSaq1SaotSamiDtc461lzD28DFwBFgO/CY1trUMM2zCWd/T9eW\n1e9p2cTcjSil+mL+47jW0W2pgzeBSVprk7lj6bK8gJ5AP8AP+EUptUFrvc+xzaqV/sA24HqgI/Cd\nUuoHrfUZxzbL/dT2Pe3uAf8w0L7C43alz9W2jCNZ1T6l1KXAXOAmrXV6A7XNWtZcQy9gSWmwbwHc\nrJQq1lovb5gmWsWa60gD0rXWOUCOUmo90ANwloBvzTUMB6Zo8wDyfqXUQeAiYFPDNLHenP09bZW6\nvKfdfUjnV6CTUipCKeUDDAY+P6fM58DQ0jv7UUCW1vpoQzfUghqvQSkVBnwGDHHSnmSN16C1jtBa\nG7XWRuBTYLSTBXuw7u9pBXCtUspLKeUP9AZ2N3A7LbHmGlIxf0NBKRUKdAH+atBW1o+zv6drVNf3\ntFv38LXWxUqpMcD/MM9OmKe13qmUeqT0+GzMM0JuBvYDuZh7N07Dymt4HggBZpX2kIu1EyWQsvIa\nnJ4116G13q2U+gb4AzABc7XWVU69cwQrfxf/ByxQSm3HPNNlktbaaTJQKqUWA9FAC6VUGvAC4A2u\n8Z4Gq66hTu9pWWkrhBBuwt2HdIQQwm1IwBdCCDchAV8IIdyEBHwhhHATEvCFEMJNSMAXQgg3IQFf\nCCHchAR8IYRwE/8Pn9dX2gavZWAAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x1818d87310>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%pylab inline\n", | |
| "ptrue = [1,1]\n", | |
| "n = 30\n", | |
| "xtrue = np.random.random(n) \n", | |
| "ytrue = polyval(ptrue,xtrue)\n", | |
| "xi = linspace(0,1)\n", | |
| "\n", | |
| "xerr = np.ones_like(xtrue) * 0.1\n", | |
| "yerr = np.ones_like(ytrue) * 0.03 \n", | |
| "xobs = xtrue + xerr * np.random.randn(n)\n", | |
| "yobs = ytrue + yerr * np.random.randn(n)\n", | |
| "\n", | |
| "pfit = np.polyfit(xobs,yobs,1)\n", | |
| "\n", | |
| "errorbar(xobs,yobs,yerr=yerr,xerr=xerr,fmt='o')\n", | |
| "plot(xi,polyval(ptrue,xi),label='True Function')\n", | |
| "plot(xi,polyval(pfit,xi),label='Fitted Function')\n", | |
| "\n", | |
| "print \"RMS(xtrue - xobs) = {}\".format(std(xobs - xtrue))\n", | |
| "legend()\n", | |
| "\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### correct x coordinate \n", | |
| "\n", | |
| "- RMS between true and corrected is less than what we put in, because the high precision y coordinate puts the point onto the fitted relationship, and the fitted relationship incorporates many data points. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 76, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "RMS(xcorr - xobs) = 0.0536557067867\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4VNW9xvHvIoRAIkJUFEGYUGuViooahQLVgFW8IGC5\nFAy2WIRysVWriFJrayvgjeppAT0JIKIRFfEQrAKKFTAgYMSoiBpEFChI4oVETAxJZp0/diK5TJhJ\nsjO3vJ/n4SEze8/sn5h5s7L22r9trLWIiEh0aRHqAkRExH0KdxGRKKRwFxGJQgp3EZEopHAXEYlC\nCncRkSikcBcRiUIKdxGRKKRwFxGJQi1DdeATTjjBJiUlherwIiIR6e233/7SWtvB334hC/ekpCSy\ns7NDdXgRkYhkjPk8kP00LSMiEoUU7iIiUUjhLiIShRTuIiJRSOEuIhKFFO4iIlFI4S4iEoUU7iIi\nUUjhLiLiS0qK8ydCKdxFRKKQwl1EJAop3EVEasrIgE2bYN06SEpyHkcYhbuISFUZGTBhApSUOI8/\n/9x5HGEBH7KukCIirnD7pOemTUeCvVJREYwbB+np7hxj7Vp33uco/I7cjTFdjDGvG2O2G2M+MMbc\n5GMfY4z5pzHmE2PMe8aY85qmXBGRJlYz2P09H6YCGbmXAbdaa7caY9oCbxtjXrXWbq+yzxXAaRV/\negGPVvwtItK03B4FJyU5UzE1eTxBGXG7xe/I3Vq731q7teLrb4EPgc41dhsCLLaOTUB7Y8zJrlcr\nIuJDaWkpOTk57rzZjBkQH1/9ufh45/kIUq8TqsaYJOBcYHONTZ2BPVUe76X2DwAREdfl5+dz6aWX\nctFFF5GXl9f4N0xNhbQ0iItzHns8zuPU1Ma/dxAFfELVGHMMsAy42Vpb2JCDGWMmABMAunbt2pC3\nEBH5QU5ODkOHDuWLL75g/vz5nHjiie68cWrqkZOnETQVU1VA4W6MicUJ9gxr7Qs+dvkv0KXK41Mq\nnqvGWpsGpAEkJyfbelcrIlJh6dKljB07lsTERLKyskhOTnb3ABEa6pUCWS1jgAXAh9baf9Sx2wrg\n1xWrZnoDBdba/S7WKSICgNfr5a677mLkyJH07NmT7Oxs94M9CgQycu8LXAe8b4ypPGMxHegKYK19\nDHgZuBL4BCgCrne/VBFp7goLCxkzZgwvvvgi48aNY+7cucRVzo1LNX7D3VqbBRg/+1hgiltFiYjU\ntGPHDoYMGUJubi7/+te/mDJlCs7EgviiK1RFJOytXr2aUaNGERMTw5o1a0iJ4Fa8waLeMiIStqy1\nPPjgg1x55ZV07dqV7OxsBXuAFO4iEpaKi4u57rrruP322/nlL3/Jxo0bSUpKCnVZEUPhLiJhZ+/e\nvfz85z8nIyODe++9l+eee46EhIRQlxVRNOcuImFlw4YNDBs2jKKiIjIzMxk8eLD/F+3fA+9thoHD\nm77ACKGRu4iEjfnz59O/f3/atm3Lpk2b/Ae7tbB+Jfz9Rli5FL49GJxCI4BG7iIScqWlpdxyyy3M\nnTuXyy67jGeeeYbExMSjv+jQt7D4Edi6AbqfC+Nug7btg1NwBFC4i0hI5efnM3LkSNauXcttt93G\nrFmzaNnSTzR99C4seBAKD8KIG+DSX0ILTURUpXAXkZB59913GTJkCF988QVPPvkkY8aMOfoLyspg\nxVOw8lk4sRPc+TAknRacYiOMwl1EQqLejb/y9kH6/bDrY+g3EEZNhNZtglNsBFK4i0hQeb1e7r77\nbmbMmEGfPn1YtmwZHTt2rPsF1sKm1+CpuRATAxOnQ/JFwSs4QincRSRo6t34q+g7yJgDm1+H03rA\nDbfD8S71bI9yCncRCYqqjb/mzJnD5MmTjzT+qmwpULWH+s7tzjTM1/kw5Dq4ahS0iAl22RFL4S4i\nTa5ejb+85fDSM/BiBhzXAaY9BKf+NGi1RguFu4g0GWsts2fPZtq0afTo0YPMzMza/WEyMmDTJigp\ngS5doO/ZcEw59OoPqTdCvNoONIQWhopIk6hs/DV16lSGDRvmu/FXRgZMmOAEO8DevbBsNZxyHoyf\npmBvBI3cRaR+Ami5u7ekhKHbtrH10CFmJCVx54EDmKuuqr1j5Yi9qrJyuO8fsHajO/VG+L1QG0rh\nLiKu2lBQwLAPPqDI62X5mWcy+IQT6t65ZrD7e14CpnAXkfo5ykh4/vz5TJ48GY/Hw+srVtC9e3ff\nK2G8XnhlGby9BQ4V134jj6fZjrjdojl3EWm00tJSbrzxRsaPH8+AAQPYsmWLE+y+HPwKHp4Ozy+A\n4VdBmxpXmcbHw4wZTV90lNPIXUQapWrjr6lTpzJr1ixiYirWo1ddCZOUBBPGQt52OFwCv74Jfn45\n/OJpGDfO2cfjcYI9NTWU/0lRwVhrQ3Lg5ORkm52dHZJji4g7cnJyGDp0KAcOHGD+/PmkVg3lypUw\nRUVHnotpAUMGwJyFcHKXI8/7mroRn4wxb1tr/TTiUbiLSFX1uPn0c3l5XP/xxxwXG8v/nXkmyW3b\nVt/B10oYgLg46N27cXVWaoY/DAINd825i0i9eK3lT7t28asPP6TnMcfw1nnn1Q520EqYENOcu4gc\n4WckXFBQwJgxY/j37t1Hb/xVeNCZY/+moPY2rYQJCo3cRSQgubm59O7dm5UrVzJnzhzS09N9B/u2\nbPjrJDi3mzMFU5VWwgSNwl1E/Fq1ahUXXngh+fn5rFmzhilLl2L696++U+lheDYNHrkLjjkWnngB\nFiw4EvAeD6SlaSVMkPidljHGLAQGAXnW2h4+trcDngK6VrzfQ9bax90uVESCr2bjr+XLl9OtW7fa\nO+7fA2n3wZ6dMGAwDB8HreIgtRukpzv7aComqAKZc18EzAEW17F9CrDdWnu1MaYD8LExJsNae9il\nGkUkBIqLixk/fjwZGRmMGDGCxx9/nISEhOpr1z0e+PWv4KuPoVVr+P09cE6vUJcuBBDu1tr1xpik\no+0CtDVO1/1jgK+BMleqE5GQ2LNnD9dccw1bt25lxowZ3Hnnnc6NNWp2cdy9G2bNhl8NgrnzoP3x\ntd9MI/aQcGO1zBxgBbAPaAv8ylrrdeF9RcQN9Vi7Dk7jr19+8AHFXi+ZP/0pV7/yCrzyirPR19r1\ncq/Tpve/w9ypVz8MXOHGCdWBQA7QCegJzDHGHOtrR2PMBGNMtjEmOz8/34VDi4ib5u/fT/933+XY\nli3ZdO65XF2zo6PWrkcMN0bu1wP3WedS10+MMbuAM4AtNXe01qYBaeBcoerCsUXEnwBGwqWlpdxy\nyy3MXbeOgQMHsmTJEhITE6vvlLcPfnIGFHxb+w20dj3suDFy3w1cAmCMOQk4HfjUhfcVkSDIz8/n\n0ksvZe7cuUydOpWXXnqperBbC2+ugXumQK8zoLXWrkeCQJZCLgFSgBOMMXuBvwCxANbax4C/A4uM\nMe8DBphmrf2yySoWEde8++67DBkyhC+++IInn3ySMWPGVN+h6DvImAObX4efnAV/+19Y9aq6OEaA\nQFbLjPazfR9wmWsViUhQLF26lLFjx5KYmEhWVhbJyTV6UX2yHebfD1/nw9DfwJUjoUWME+Raux72\ndIWqSDPj9Xq56667GDlyJD179iQ7O7t6sJeXw4sZ8MBtgIFps2HQaCfYJWKocZhIM1JYWMiYMWN4\n8cUXfTf++uoAzH8AdnwAvfpD6o0Qn1D7jTRiD3sKd5FmYseOHQwZMoTc3FzmzJnD5MmTnQuTKr21\nHhb/j3MCddxU+NkloStWGk3hLtIMrF69mlGDBhEDvPrqq/Sv2vTr+2JY8ihseAV+dAbccDuc2Clk\ntYo7FO4iUcxay0MPPcQdd9xBjzZtyOzRg6Sqwf5ZLqTf76xhv2oUXD0GWioWooH+L4pEqVqNv/bv\nJ6HyxtVeL6xeBssXwbGJcNv9cPrZIa1X3KVwF4lCtRp/eTyYyrXpXbrAz3tCfCmc3w+uuwmO8XGb\nPIloWgopEmU2bNhAcnIyubm5ZGZmMt3jwVTt5Lh3LyxdCZ3Ohol/UrBHKY3cRcJRPTs5Vkrfv58p\nO3bgad2atWeeSffZs313ciwrhwf+CetrtYBqOC2PDCsKd5EoUOr1csvOnczdt4+BiYks6d6dxNhY\nZ6M6OTZLCneRcFSPUXB+fj4jRoxg3b59TJ06lVmzZhETE+OsV38tE97eAoeKa79QnRyjmsJdJILl\n5OQwdOjQ2o2/Cg/C47Ph/bfgmsvh+VVQXCXg1ckx6umEqkiEWrp0KX379qWsrIysrKwjwb4tG/46\nCT7MgWsnwxPLnEZflW0GPB5IS1MnxyinkbtIhPF6vdx9993MmDGDPn36sGzZMjp27Ailh+GFRfDq\nC9A5CW6d5fwN6uTYDCncRSJInY2/9u+BtPtgz04YMBiGj4NWNW6qoVBvVhTuIhHCZ+MvgPUr4ZnH\noFVr+P09cE6vUJcqYUDhLhIBVq9ezahRo4iJiWHNmjWkpKTAoW9h8SOwdQP89Dz47a3Q/vhQlyph\nQuEuEsastcyePZtp06bRo0cPMjMzSUpKgo/ehQUPOqtiRoyHS6+BFlofIUco3EXCVNXGX8OHD2fR\nokUkxMU5J01XPuu05Z3+MHhOC3WpEoYU7iJhqGrjr3vvvZfp06dj8vfDI/fDro+h30AYPQniWoe6\nVAlTCneRMLNhwwaGDRtGUVERmZmZXD1oEGx6DZ6aCzExTrOv5J+HukwJcwp3kTCSnp7OlClT8Hg8\nvP7663T3dHXuabr5dTitB4yfBsd1CHWZEgEU7iJhoLS0lFtuuYW5c+cycOBAlixZQuLX++Fvk+Hr\nfBj6a7jyV9AiJtSlSoRQuIuE2A+Nv9atcxp/zbiXmFVL4cUMOO5EmDYbTu0e6jIlwijcRUKoVuOv\nKy6Df9wJOz6A3gMgdQq0SQh1mRKBFO4iwVZxI46lU6YwduxYEhMTeeONN7iAYqfhl7Uwbir87JLQ\n1ikRTeEuEkwpKXjfeYc/t2/PzJEj6dOnDy88ncFJ/3kBNrwCPzrDOWna4eRQVyoRzm+4G2MWAoOA\nPGttjzr2SQEeAWKBL621F7tZpEi0KCgrY0xREf8uLOSGG25g7tSbaPW/f4P8L2DQtc6flhpzSeMF\n8l20CJgDLPa10RjTHpgHXG6t3W2MOdG98kSiR25uLmlbtjCnrIwVACuWY3a+Bef3gNvug9PPDnWJ\nEkX8NqOw1q4Hvj7KLtcCL1hrd1fsn+dSbSJRY9WqVazr0YMHSkvxAAYweV/Chu3Qvb+CXVznxu9/\nPwFijTFrgbbA/1hrfY7yRSJKxYnPxrDW8tDeveR8+ilP4mM0dbgUJk+BxU82+liAerbLD9wI95bA\n+cAlQBvgTWPMJmttbs0djTETgAkAXbt2deHQIk0kJQVychr1FsXWMr64mIzSUvI4yq/JJSWNPhYA\nPXs2/j0kargR7nuBr6y13wHfGWPWA+cAtcLdWpsGpAEkJydbF44t0nQaEZZ7vv+eaz74gK2lpaT1\nOJ0Ttn1c985xce4Es0btUoUb4Z4JzDHGtARaAb2Ah114X5HQqW9QVk7hrF1LVlYWw4YNo9gY3nv4\n7/T4aDN8vge+Lar9OmNgwQLdrFpc5/eEqjFmCfAmcLoxZq8xZpwxZqIxZiKAtfZDYBXwHrAFmG+t\n3daURYuEq/T0dAYMGEDS8YnsmzaOHtuynLskPfwIxMdX39kYmDhRwS5Nwu/I3Vo7OoB9HgQedKUi\nkQhU6vVy886dzJswgemDB/L3bu1osecTuHYy9L/aCfLW8TBunDPHHhenEbs0KV0tIdJIhY89RlFW\nFv+yltnHxNO6YB+0S4Kp90PnpCM7pqZCerrztebHpYnpposijfDZzJm0nDyZjtbSAmh9qAje/Ah+\n9LPqwV5p7VoFuwSFRu7SvLiwdr3S0vx8em3fTnzNDYdL4Xe/g8cfd+dA+mEgDaBwF6knr7Xc/dln\nPHpgP1/WtVNJSTBLEqlF4S7NSyNHwQUFBYwZM4ZDpYfYMeZSePZ1OFRce0ePRyNuCSnNuYvUJSWl\n2jRObm4u/X7Wm74HcvnPlReQ2PFkzL1/r73EMT4eZswIaqkiNWnkLhKAVatWcee43/BE759wXmIC\n/PxyGDUR4lrDCR2PLHH0eJxg1xJHCTGFu4gvGRmwaRO2pITCxEQ+i4OsKy6gdcIxcP0f4fx+R/bV\nEkcJQwp3kZoyMmDCBCgpwQDtDh5kYkwLvN/F0GL2/8JxHWq/RqEuYUbhLtHBxSWObNpUe7VLuZcW\nK16HvBHuHEM/DKSJ6YSqSA22rmWMWt4oEUQjd4kOLo2El8x5hGu2vEnr4sO1N2p5o0QQjdxFgNLS\nUhZM/A1XbH4RLjgdGxdXfQctb5QIo3CXZu/LvbtZk3oF48oOUBDfjthnX8YsWOB0bgRnxJ6WpuWN\nElE0LSPN2serX6TVEw8zsG0s73vO4qw7Z0HLllreKBFP4S7Nk9fLew/eTfePt5BvDLlDbuCswS6t\nhBEJAwp3aXa8X+Wx689TOPvwt/znUDln3p/GGd1Orb2jRuwSwRTu0qx8l7WG8oUP0bG8jEWxHRi9\n+HHiWrcOdVkirlO4S/NQ8j0H0x+kfc4Gtn5TyIf9hvCbP96OMSbUlYk0CYW7RL89n3LoH3fR/tuv\n+deOLzj77tmkXnJJqKsSaVIKd4kOle0Hqs6TW4tds5zy59Io/K6Y2/YWceeTmXg8nlBUKBJUCneJ\nToUHKZ//ADHbt/Ly53msOP5U/mflIhISEkJdmUhQKNwl+mzLpjz9AUq/LeAPmz6ic+oE0qdP1/y6\nNCsKd4k8NadgKnqvU1ICHY6HMzuTe2IiYzd+zJ8fm8+gQYNCValIyCjcJbJV6b0OwJdfU/rGNyw6\n8SQW/ed1unfvHtr6REJE4S7B5Ubf9ZycI+/lo/d6rNcy85tviJk0qfHHAl3MJBFJ4S6R5cABKCwE\na33fVKNCjHqvSzPnN9yNMQuBQUCetbbHUfa7AHgTGGWtfd69EiWqNGYUXDkFY63zuKQEC/g8Tare\n69LMBTJyXwTMARbXtYMxJga4H3jFnbIkajTx7e8M4KVG7+oWLaBNm4YfWz8UJAr47edurV0PfO1n\nt98Dy4A8N4oS8amOqRYDULnMMS4OfvITOOmkoJUlEo4aPedujOkMXAP0By5odEUSXeo7Ck5JcU6Y\n9ux55LXWwsY1cPkWOFRc6yXG44GkpIYdTyRKuXFC9RFgmrXW6+8iEWPMBGACQNeuXV04tES9okPw\n1L9gyzpKLuuH/b81tK6cc4cjt7/TXZJEqnHjNnvJwDPGmM+A4cA8Y8xQXztaa9OstcnW2uQOHTq4\ncGiJKhkZsGEDFBTAunWQ2B5+dTlkv0HumX04+bUt3JSQwPexsc7+uv2dSJ0aPXK31nar/NoYswj4\nt7V2eWPfV5qRlBRnieOnn0JZ2ZHnDxZgX9rIhmsGc/Ft93DWWWdx5/LltB471tmuKRiROgWyFHIJ\nkAKcYIzZC/wFiAWw1j7WpNVJeGvf3p33OXQIyst9bjLWcsoLmQyPjWXhrl0k9OzZuOMfPNjAIkUi\ni99wt9aODvTNrLVjG1WNND+HD9cZ7JU8wDNt2qjxl0g96ApVabjGjoIrL0ryw3g88NlnjTuWSDOj\ncJfAuXlBEhy1fUA1uiBJpN4U7hI6dQR75UJHExMDp52mC5JEGkDhLoFzaxT8zZew4CHI3gzffV97\nuzGYiy7SqFukEdxY5y4SuHc2wl8nwacfwh//gLd16+rbW7TAnH66gl2kkRTuEhwl38OT/4S5f4Pj\nT4K75/Lsmefx2/Jy9lIxFePxwOLF8OGHIS5WJPJpWkbc46svDMCeTyHtPti/GwYOp3zwGP58z9+Y\nNWsWffv2JdbrxbRqpdG6iIsU7tJ0rIXXMuH5BZBwDPxxJgWdTyV12HBeeuklxo8fz5w5c2jVqlWo\nKxWJOgp3cU/lXZLWrYOuXeCiZGhdDGf3gutvIXd/HoN79WLnzp3MnTuXSZMm6cIkkSaiOXdxx+TJ\n8NFHR+6StGcvPPsinHgm/P6vrMx6kwsvvJCvvvqKNWvWMHnyZAW7SBMytmr71CBKTk622dnZITm2\nVHDroqQDB5xg98HGxfFgp07csWsXZyckkNmjB56aK2QCpTl5EYwxb1trk/3tp2mZ5qry5KcbCgvr\n3GRLSpi2axcjY2NZGBNDQh0/BPyq2jBMRPxSuDdnbgXmunV1btoN3JuUxPSuXRs3DaNRu0i9KNyb\nKzfC8lAhPPFInVeaeoGvb72VPz30UOOPJSL1ohOq0jAf5jhXmr63BW6c4NzurgovcHDUKM5TsIuE\nhEbuUrfKm2FUbe1bVgrLn4TVS+GkzvCHe6Drjyn76fkU/Pa3JJaXk9+iBQmPPspxAbTzFZGmoXCX\nwB34L6TfB5/tgIuugF/9DuJak5+fz/AFC1hfXs7tXbowc9cuYmJiQl2tSLOmcBffMjKcVTDWOj1f\nxo6GLz+CmJYw6S44vx8AOTk5DBkyhLy8PDIyMrj22mtDXLiIgObcxZeMDPjtb49ckLR7N9z7AHxd\nBn999Idgf+655+jTpw9er5esrCwFu0gY0cg9Wrh5l6QNG6CsrPpzXgvPr4b9Iyi3lj/v2sWsPXvo\nc+yxLOvShY633lq/Y2hpo0iTUrhHAzcvSILawV7l+YJ33iG1qIiXysq4ITaWOUDc9u3uHVtEXKFw\njxZuXZDU0sBra31uskCv0lJ2lpcz98c/ZlKnTuoPIxKmFO7RoCFTHJXTOFVfu2Wtc0ONrFgoKa31\nkq+ArxISePXll0lx+2bZIuIqhXtz4yuUvy+Cp+fBxjVwand45BH4w81QeiTgvwdmd+nCW+vXk5SU\nFKxqRaSBFO7N0YEDsGsXlJRA585w/qlw0jEw6Fq4OhViYqBtO7zXXQfWshtY0asXd732GgkJCaGu\nXkQCoKWQzc2BA5Cb6wQ7wL59sHIjnHExDP21E+zAnosu4gJjaAksmTmT37/5poJdJIKon3skCELf\ndeLioHdvALIKChj2wQcUe7083b07g44/vn7H0TJHkSYTaD93vyN3Y8xCY0yeMWZbHdtTjTHvGWPe\nN8ZsNMac05CCpYlVjtjrUjGST9+/nwHvvku7li3ZfO659Q92EQkLgcy5LwLmAIvr2L4LuNha+40x\n5gogDejlTnnNmNurUXbtAq+3zs02Lo4bd+xg3r59DExMZEn37iTGxtbvGBqxi4QNv+FurV1vjEk6\nyvaNVR5uAk5pfFnSIFlZzt/9+lV/vpU5Msfug23RgntbtGDevn1MPeUUZv3oR8Ro/bpIRHN7tcw4\nYKXL79k8NWQUXNmit/K11sJrmfD8AmgbD98W1XqJbdGCm9q3J72oiKeeeorU1NQGlywi4cO1cDfG\n9McJ935H2WcCMAGga9eubh1aoHoXx6QkuGs6HNoD296Cc3rBw4OdtetFRwK+zBhuiInhtfh43li9\nmuRkv+doRCRCuBLuxpizgfnAFdbar+raz1qbhjMnT3JycmiW6USTynn58eOrd3H8/HP43STofw78\n6R5IGQTGQOt4GDcOW1LCwWOPZUphIZ9fcAFvLVtGx44dQ/afISLua3S4G2O6Ai8A11lrj7IcQ6px\n44RpZbOwsWN9dHH0wvr3oWw23DP7h6cL4uJILS/npcJCxp98Mv+JiSFu1KjAjqcTpiIRw2+4G2OW\nACnACcaYvcBfgFgAa+1jwN3A8cC8iiZSZYGswWz2GtvF8fBhKC4++j6lZdWOk1tezuCiInZ6vcxt\n3ZpJ332Hee+9wI7nVmMyEQmKQFbLjPaz/QbgBtcqai4aE5b+1qz7OM6qr79m1PbtxMbEsOass7i4\n8uRroDRqF4ko6i0TLG6uW/ezZv0HMTFYa3lo717u+PRTzkpIYHmPHiS1bh3YcRToIhFL4R6JjrJm\nvarDp57K9R99xNN5eYzo0IHHTz+dBN24WqRZULgHS31Hwb76rQOUlcLJJ8OXdS5KAo+Hr/74RwYu\nXszWHTuYMWMGd955p26sIdKMKNzDla8Trgf+C+n3wZmd4M1COFzlhhrx8c5SyFatyHrqKYYNG0Zx\ncTErVqxg0KBBwatbRMKCWv5GAmthwyvwtymQtx9mz4OFjzudHAE8HkhLg6Ii0h98kAEDBtCuXTs2\nb96sYBdppjRyDyeVK1iuvRYKCpyvW7aEvufDaYlw+tkwbioc1wHO7wfp6c4+a9dSWlrKzVOmMG/e\nPAYOHMiSJUtITEwMzX+HiIScRu7hpqgIHn30yOPycli/BfZ74dZZTrDXkJ+fzy9+8QvmzZvH7bff\nzksvvaRgF2nmdLMON9R3zbgvRUXV7lnqU7t2tZ7KKS9nyHffkWct89u0IbVVK//H6tlTyxxFIpRr\nN+uQIAgk2H147vBh+hw6hBd4IyEhsGAXkWZBc+4N5cZovVI9g73cWv5cUsKskhL6xsTwfHw8HVv4\n+Tl98GAjChSRSKNwjxQVd0UqsJbUoiJeKitjfGwsc9q0oZXWr4tIDQr3hnJrJPzpx/Dj7kfa9dYU\nEwMTJsC8eeTm5jJ48GB27tzJvHnzmDhxoi5MEhGfFO6h4i2HVc9D5mI458eQs6P2PrGxTvdHYOXK\nlYwePZrY2FjWrFnDxRdfHOSCRSSS6IRqKHzzJfxjOrzwOJzbB7LehkmTjmyPiXEeHz6MtZYHHniA\nq666im7dupGdna1gFxG/NHIPtnc2wqKHofQwjL0F+l7m3CVp3jx4+mlnn4opn+LiYm644Qaefvpp\nRo4cycKFC0lISAhh8SISKRTuwVLyPTyXButeBs9pMH4adDyl+j5Verzv2bOHa665hq1btzJz5kzu\nuOMOza+LSMAU7sGw51NIuw/274aBw+Ga30DL2Nr7VVxYlJWV9UPjr8zMTK6++urg1isiEU/h3pS8\nXngtE5YthIS28MeZ8NPzjvqStLQ0brzxRpKSkli7di3du3cPUrEiEk0U7m7w1Xu94Bt4fDZsy4Zz\nejnz623rvvDp8OHD3HzzzTz66KNcfvnlLFmyhPZuXiglIs2Kwr2xUlKc3utV74n6/luwcDZ8XwSp\nN0LKVc4mNbVjAAAIG0lEQVRJ0zrk5eUxYsQI1q9fz+23387MmTOJ0R2TRKQRFO6NdeAAFBbCunWQ\nlARjr4W9W6FzEtx2n/P3UbzzzjsMHTqUvLw8nnrqKVJTU4NRtYhEOa1zb4zJk+Gjj45cXfr55/DA\nI3BsN7jrn36D/dlnn6Vv3754vV6ysrIU7CLimuY1cq+cG3fDgQNOsNdUXAzzHoetH9b5Uq+1/Pmz\nz5i5ezd9jz2WZV26cNKtt/reWa15RaQBmle4+7ovaUMVFta9raSkzmMVWMuYoiL+Xdn4C2i1fbvv\n96k6jy8iUg/NK9zdDMt16+reFhfn81i5RUUM3raNneXlzP3xj5nUqdPRL0zSqF1EGqh5hbubYZmU\n5Myx12QMLFgANebPf2j8lZDAmpUr1R9GRJqUTqg21IwZEB9f/TljYOLEasFetfFXUlKSGn+JSFD4\nDXdjzEJjTJ4xZlsd240x5p/GmE+MMe8ZY45+CWa0SE2FtDRnCgacv5980mkAVqG4uJgxY8Ywbdo0\nRowYwYYNG/B4PCEqWESak0BG7ouAy4+y/QrgtIo/E4BHG19WiKSk1G9FTWoq9O7t3Li6d+9qI/Y9\ne/bQr18/lixZwsyZM3nmmWfU0VFEgsbvnLu1dr0xJukouwwBFltrLbDJGNPeGHOytXa/SzWGNx/z\n+FUbf61YsYJBgwYFvy4RadbcmHPvDOyp8nhvxXPN0vLlyxkwYADt2rVj8+bNCnYRCYmgnlA1xkww\nxmQbY7Lz8/ODeWj/MjJg06YjbQQyMhr0NhdeeCGjR49my5Yt6ugoIiHjRrj/F+hS5fEpFc/VYq1N\ns9YmW2uTO3To4MKhXZKR4dyEuqTEefz5587jBgR8p06deOKJJ9TRUURCyo117iuAG40xzwC9gIIm\nn293s40AOCP2ymCvVFQE48ZBero7x9AFSSISRH7D3RizBEgBTjDG7AX+AsQCWGsfA14GrgQ+AYqA\n65uq2CZTM9j9PS8iEuYCWS0z2s92C0xxraJAuD0KrutqU49HI24RiUi6QhV8X20aH+88LyISgRTu\nUPtqU4/Heaz+6iISoZpX47CjSU09cvJUUzEiEuEU7lUp1EUkSmhaRkQkCincRUSikMJdRCQKKdxF\nRKKQwl1EJAop3EVEopDCXUQkCincRUSikMJdRCQKGaepYwgObEw+4KMVo08nAF82YTlui7R6IfJq\njrR6IfJqVr1NryE1e6y1fu92FLJwrw9jTLa1NjnUdQQq0uqFyKs50uqFyKtZ9Ta9pqxZ0zIiIlFI\n4S4iEoUiJdzTQl1APUVavRB5NUdavRB5NaveptdkNUfEnLuIiNRPpIzcRUSkHsIy3I0xxxljXjXG\n7Kj4O9HHPl2MMa8bY7YbYz4wxtwUgjovN8Z8bIz5xBhzh4/txhjzz4rt7xljzgt2jTXq8VdvakWd\n7xtjNhpjzglFnTVqOmrNVfa7wBhTZowZHsz6fNTht15jTIoxJqfi+3ZdsGusUYu/74l2xpgXjTHv\nVtR7fSjqrFLPQmNMnjFmWx3bw+ozV1GTv5qb5nNnrQ27P8ADwB0VX98B3O9jn5OB8yq+bgvkAj8N\nYo0xwE7gR0Ar4N2axweuBFYCBugNbA7hv2kg9fYBEiu+viKU9QZac5X9/gO8DAwP53qB9sB2oGvF\n4xPDvN7plZ8/oAPwNdAqhDVfBJwHbKtje9h85upRc5N87sJy5A4MAZ6o+PoJYGjNHay1+621Wyu+\n/hb4EOgctArhQuATa+2n1trDwDM4dVc1BFhsHZuA9saYk4NYY1V+67XWbrTWflPxcBNwSpBrrCmQ\nf2OA3wPLgLxgFudDIPVeC7xgrd0NYK0NZc2B1GuBtsYYAxyDE+5lwS2zSjHWrq+ooS7h9JkD/Nfc\nVJ+7cA33k6y1+yu+/gI46Wg7G2OSgHOBzU1bVjWdgT1VHu+l9g+XQPYJlvrWMg5nBBRKfms2xnQG\nrgEeDWJddQnk3/gnQKIxZq0x5m1jzK+DVl1tgdQ7B+gO7APeB26y1nqDU16DhNNnriFc+9yF7AbZ\nxpg1QEcfm/5U9YG11hpj6lzSY4w5BmfUdrO1ttDdKpsnY0x/nG+yfqGuJQCPANOstV5ncBn2WgLn\nA5cAbYA3jTGbrLW5oS2rTgOBHGAAcCrwqjHmDX3W3Of25y5k4W6t/UVd24wxB4wxJ1tr91f8SuXz\nV1djTCxOsGdYa19oolLr8l+gS5XHp1Q8V999giWgWowxZwPzgSustV8Fqba6BFJzMvBMRbCfAFxp\njCmz1i4PTonVBFLvXuAra+13wHfGmPXAOTjnjIItkHqvB+6zzoTwJ8aYXcAZwJbglFhv4fSZC1hT\nfO7CdVpmBfCbiq9/A2TW3KFiDnAB8KG19h9BrK3SW8BpxphuxphWwCicuqtaAfy64gx+b6CgynRT\nsPmt1xjTFXgBuC5MRpJ+a7bWdrPWJllrk4DngckhCnYI7HsiE+hnjGlpjIkHeuGcLwqFQOrdjfNb\nBsaYk4DTgU+DWmX9hNNnLiBN9rkL9ZnkOs4eHw+8BuwA1gDHVTzfCXi54ut+OCd73sP5tTEHuDLI\ndV6JM+LaCfyp4rmJwMSKrw0wt2L7+0ByiP9d/dU7H/imyr9ndhh8Lxy15hr7LiKEq2UCrReYirNi\nZhvOdGLY1lvxmXul4vt3GzAmxPUuAfYDpTi/BY0L589cgDU3yedOV6iKiEShcJ2WERGRRlC4i4hE\nIYW7iEgUUriLiEQhhbuISBRSuIuIRCGFu4hIFFK4i4hEof8HY5DHC2gmWkAAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x1818a91c50>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "xcorr = (yobs - pfit[1]) / pfit[0] \n", | |
| "plot(xi,polyval(ptrue,xi),label='True Function')\n", | |
| "plot(xi,polyval(pfit,xi),label='Fitted Function')\n", | |
| "errorbar(xcorr,yobs,yerr=yerr,xerr=xerr,fmt='or')\n", | |
| "print \"RMS(xcorr - xobs) = {}\".format(std(xcorr - xtrue))" | |
| ] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 2", | |
| "language": "python", | |
| "name": "python2" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 2 | |
| }, | |
| "file_extension": ".py", | |
| "mimetype": "text/x-python", | |
| "name": "python", | |
| "nbconvert_exporter": "python", | |
| "pygments_lexer": "ipython2", | |
| "version": "2.7.15" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 2 | |
| } |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Simulate data and fit it\n", | |
| "\n", | |
| "- add errors both in x and y\n", | |
| "- RMS between true and observed x coord is equal to what we put in " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 75, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Populating the interactive namespace from numpy and matplotlib\n", | |
| "RMS(xtrue - xobs) = 0.0947395365753\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.legend.Legend at 0x1818d0cb10>" | |
| ] | |
| }, | |
| "execution_count": 75, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlclNX+wPHPYRtAEAUVTYVBU0szKS1puYlZVyusbLVI\nzVBKU7PS7Fa3uv6iLC2zm2ZoaiZpi6nZdsvc2tTULHfTBMRdQJR9mfP7Y4BAYRhghplhvu/Xi1fN\nPGeecx5hvnPmPOd8j9JaI4QQovHzcHQDhBBCNAwJ+EII4SYk4AshhJuQgC+EEG5CAr4QQrgJCfhC\nCOEmJOALIYSbkIAvhBBuQgK+EEK4CS9HVdyiRQttNBodVb0QQrikLVu2nNJat6zLax0W8I1GI5s3\nb3ZU9UII4ZKUUil1fa0M6QghhJuQgC+EEG5CAr4QQrgJh43hV6WoqIi0tDTy8/Md3RRRC76+vrRr\n1w5vb29HN0UIYYFTBfy0tDQCAwMxGo0opRzdHGEFrTXp6emkpaURERHh6OYIISxwqiGd/Px8QkJC\nJNi7EKUUISEh8q1MCBfgVAEfkGDvguR3JoRrcLqAL4QQwj4k4FeQnp5OZGQkkZGRtG7dmrZt25Y/\nLiwstFk9q1atIigoqPzc/fv3t9m5AbZu3co333xT/njZsmVMnTrVpnUI0RhFR0cTHR3t6GbYjVPd\ntHW0kJAQtm3bBsCLL75IQEAAEyZMqFRGa43WGg+P+n1W9u3bl+XLl9frHNXZunUrO3bsYMCAAQAM\nGjTILvUIIVyL9PCtsH//frp27UpsbCzdunXj0KFDNGvWrPz4kiVLGDFiBADHjx/njjvuoFevXlx5\n5ZVs2LDB6noeeOCBSh8CAQEBgPkbQb9+/bjjjjvo0qULQ4cOLS+zceNGrrrqKnr06EHv3r3Jyclh\n8uTJJCUlERkZyaeffsrcuXMZP348AAcPHqRv375ceuml3HjjjaSlpZXX/dhjj3H11VfToUMHli1b\nVvd/MCGEU3LeHv6S2ZB6wLbnDOsIgx+p00v37NnDwoUL6dWrF8XFxdWWGzduHE899RRRUVEkJycT\nExPDjh07ziu3Zs0aIiMjARg8eDBPP/20xfq3bt3Kzp07CQ0NJSoqig0bNhAZGcngwYNZunQpl19+\nOVlZWfj6+vL888+zY8cO3nzzTQDmzp1bfp7Ro0czYsQIYmNjSUxMZPz48Xz66acAnDhxgp9++ont\n27dzzz33yDcD4VaSkpLYsGEDBQUFGI1GEhISiI2NdXSzbMp5A76T6dixI7169aqx3KpVq9i7d2/5\n48zMTPLy8vDz86tUrrZDOlFRUVxwwQUAREZGkpycjMFgICwsjMsvvxyAoKCgGs+zceNGvvjiCwCG\nDh3Kv//97/Jjt99+O0opLr30Ug4fPmx124RwdUlJScTHx1NQUABASkoK8fHxAI0q6DtvwK9jT9xe\nmjRpUv7/Hh4eaK3LH1ecg661ZtOmTfj4+NS6Di8vL0wmEwAlJSWVvkkYDIby//f09LT4LaOuKtZR\n8fqEcCb2uKla1rOvKDc3l7i4OObMmWPz+tauXWvzc1qjxjF8pVR7pdQapdQupdROpdRjVZRRSqm3\nlFL7lVJ/KKUut09znYOHhwfNmzfnzz//xGQyVRrvvuGGG5g5c2b547KbwNYwGo1s2bIFMM+sKSkp\nsVi+a9eupKamsnXrVgDOnDlDSUkJgYGBnD17tsrXREVF8fHHHwOwaNEirrvuOqvbJ0RjdW6wr+l5\nV2VND78YeFJrvVUpFQhsUUp9p7XeVaHMTUCn0p/ewDul/220Xn31Vfr370+rVq3o2bNn+R/GzJkz\nGTVqFPPnz6e4uJi+fftW+gCw5OGHH+a2227jiy++ICYmplKPuyoGg4HFixczatQo8vPz8fPzY/Xq\n1Vx//fVMnTqVyy67jGeffbbSa2bOnMlDDz3EK6+8QmhoKPPnz6/bP4AQDmKP3rHRaCQl5fw08+Hh\n4Q7rjduDqu1Xd6XUCuBtrfV3FZ57F1irtV5c+ngvEK21PlrdeXr16qXP3QBl9+7dXHzxxbVqj3AO\n8rsTTq+oEI4fhnbn53wqG8PPzc0tf87f35/ExESnG8NXSm3RWtd8Q7EKtZqWqZQyApcBG8851BY4\nVOFxWulzQgjheMfS4JXH4fWnIT/3vMNls9bKvlWHh4c7ZbCvL6tv2iqlAoClwHit9Zm6VKaUigfi\nAcLCwupyCiGEqJ1fVsGit8HLG4Y/Cb7+VRaLjY0tv0HbmIZxKrIq4CulvDEH+ySt9WdVFDkMtK/w\nuF3pc5VorROBRDAP6dS6tUIIYa38PEiaaQ74nbvDiKcg2PLe34010JepMeArcyrE94DdWus3qin2\nOTBGKbUE883aLEvj90IIYVep++HdV+DEUbj1AYi5Dzw8Hd0qh7Omh38NMATYrpQqm2P4DBAGoLWe\nDXwF3AzsB3KB4bZvqhBC1EBr+H4FfPoeBAbBhCnQ5VJHt8pp1BjwtdY/AhYTnmvzVJ9HbdUoIYSo\ntewzMP91+H0jXNobhj9hDvqinCRPO4enp2d52uKyFAabN29m3LhxgHmM7+effy4vv3z5cnbt2lXd\n6apVlhjNmvpt5fTp08yaNav88ZEjR7jrrrtsdn4hHGbfdvjPaNi51bxKf+yLEuyr4LypFRzEz8/v\nvNWxRqOxPI/O2rVrCQgI4OqrrwbMAT8mJoauXbvarX5bKQv4o0ePBuCCCy4oT5wmhEsylcCXS+Dz\nJGjZGv71BoR3cnSrnJb08K2wdu1aYmJiSE5OZvbs2UyfPp3IyEjWrVvH559/zsSJE4mMjOTAgQMc\nOHCAAQMG0LNnT/7xj3+wZ88ewJyW+KqrrqJ79+4899xztap/wYIFjBkzpvxxTExM+WyCgIAAnn32\nWXr06EFUVBTHjx8HzGmaBw0aRI8ePejRowc///wzTz/9NAcOHCAyMpKJEyeSnJzMJZdcApjzAQ0f\nPpzu3btz2WWXsWbNmvK677jjDgYMGECnTp146qmn6vvPKYRtZJ6CaU/Dig/gyj7w/NsS7GvgtD38\ntz/J5ECa7XaZAujYzocxdze3WCYvL688bXFERESlPDlGo5FHHnmk0sYot956KzExMeVDI/369WP2\n7Nl06tSJjRs3Mnr0aFavXs1jjz3GqFGjGDp0qMVUC5bqr0pOTg5RUVEkJCTw1FNPMWfOHJ577jnG\njRtHnz59ynPyZGdnM2XKFHbs2FH+DaLicNHMmTNRSrF9+3b27NnDP//5T/bt2weY8wH99ttvGAwG\nunTpwtixY2nfvn1VzRGiYfy+0TxeX1hgnlt/zY2ObpFLcNqA7yj1GVLJzs7m559/5u677y5/rizH\nzk8//cTSpUsBGDJkCJMmTbJJ/T4+PsTExADQs2dPvvvOnPFi9erVLFy4EDDfFwgKCiIzM7Pa8/z4\n44+MHTsWgIsuuojw8PDygN+vX7/y1Mtdu3YlJSVFAr5wjKJCWDofVi2D9h3g4X9Ba+v+Fv88VMiO\nAwUMig60cyOdl9MG/Jp64s7IZDLRrFmzagO2eUlD7VVMmwyV0zF7e3uXn7ch0ibbqw4hanT8sHlu\nfep+uP5WuHsEeNechrywSPP+l1l8tOoMwU09GRDVBD9f9xzNds+rrodzUw9XfNy0aVMiIiL45JNP\nAHNO+d9//x2Aa665hiVLlgDmRE21YTQa2bZtGyaTiUOHDrFp06YaX9OvXz/eeecdwJxbPysry2La\n5H/84x/l7dq3bx+pqal06dKlVu0Uwm42rIbJY+DUMXj0ebh/tFXBfvv+fEYkHGXxt2fo37sJ7z3X\nxm2DPUjAr7WBAweybNkyIiMj+eGHHxg8eHB5KuIDBw6QlJTEe++9R48ePejWrRsrVqwAYMaMGcyc\nOZPu3bvXejepa665hoiICLp27cq4cePKd7iyZMaMGaxZs4bu3bvTs2dPdu3aRUhICNdccw2XXHIJ\nEydOrFR+9OjRmEwmunfvzr333suCBQtqTM8shN0V5MO812HuaxDWAV6YBZddXePLcvNNzPgog8fe\nOEFxiWbquFZMHBJCoL97h7xap0e2FUmP3LjI707Y3KG/4N2XzUM5t9wHA2PBs+b0CJt25vHG4gxO\nZpYwKDqQuIFBjapXX5/0yE47hi+EcFNaw5qV8PEcaBIIT06Bi3rU+LKs7BLeWXqabzfmEN7ai7ee\nDKVbB/mWWpEEfCGE88g+C+9Ph99+hu5XwENPQmAziy/RWrP+tzze+iiDMzkmhtzUlNgBQfh4122S\nRGPmdAFfa13n2SzCMWTDc2ETf+6AOa9CVibcEw833A4elodi0rNKmLEkgx9/z6NzmA+vjQ2mY7ua\nb+a6K6cK+L6+vqSnpxMSEiJB30VorUlPT8fX19fRTRGuylQCX30Mn38AIaHm9AjGzhZforXmm19y\nmLU0k6JiiL+9GXf3C8TTUxEdHQ00/tz2deFUAb9du3akpaVx8uRJRzdF1IKvry/t2rVzdDOEKzqd\nDnOnwp5tcGU0DBkLfk0svuToqWLe+DCDLXvyufRCA0/GBtM+1Lth2uvinCrge3t7ExFx/gbDQohG\naPuv8N40KMyHB58wp0ew8M2+xKRZvvYs732ehYcHjB/cnJhrA/Dw+Ps1SUlJbNiwgYKCAoxGIwkJ\nCY1uX9r6cKqAL4RwA8VF8NkC+HYptIuA+H/BBZb3uE4+WsS0RensOljIld18eeK+YFoFVw5fSUlJ\nxMfHl6czSUlJIT4+HkCCfimnmocvhGjkTh5l9+PDuNjXg2Wni3knvZhCiyHIE1rEoFrcCqZ89LFF\ncOaXKkuW9ezPZTAYiIqKsk37K3DUPQKZhy+EcH6b1sLCt2jnrXj+aCHrc0yWy/tGoC6IQ/mGobM2\noI99ACVVpwYBqgz2lp53RxLwhRD1UuOsmIJ8WPwO/Pg/6HgxgfFPMzkktNrz5ReaeP+LLD75/izN\nm3oyfnBzrulxD3CPxXYYjUZSUlLOez48PFxm7JSSgC+EsJ+0g+YMl8cOwc2D4dYHwKv6sPP7vnym\nJWVw+GQxt1zThIcHNSfAyvw3CQkJxMfHk5ubW/6cv78/CQkJ9b6MxkICvhCizqqdFaM1rPsKlsyG\nJgHweAJ0rT7pX06eicTlp1n5QzZtWngx7bFWXN6ldms7ym7MxsXFUVBQQHh4uMzSOYfctBVC1EnZ\nrJhze9SJ/32L2IKjsOVH6NYTHpoAQdXvb7Fhex7TF2eQnlXCndcH8mBMEH6Guic7a+wLr+pz01YC\nvhCNXFkAtLXqZsW0a+LLwXuvY056MR+fLqHaCOMZiGodiwq6Gp2fhj46F/L+qrJoYw3edSGzdIQQ\nDa662S+Hc/IZk1bIngILncmmvVGth4CnP/rkMvSplaBlJzV7kx6+EKJOqp0V0749yampVb7m5Oli\n3lycyS/b8+gS7sPEB4Lp0FaSndVGfXr4jWdXACFEg0pISMD/nKR5/v7+JLzyynlltdZ88WM2D00+\nytY9+TxyRzPenhgqwb6ByZCOEI1Mg9y0LC4m1lAAvTvx1OY/OZKTX+2smMMni3g9KYNt+wqI7GxO\ndta2pSQ7c4QaA75Sah4QA5zQWl9SxfEgYBEQVnq+aVrr+bZuqBDCSZw8BnOmwF97iI0byftqKZ30\n+R8wJSbNZ2vOMu/zLLw84Yn7g7nlmiaS+tyBrBnSWQAMsHD8UWCX1roHEA28rpSS72lCOEDZvPh1\n69ZhNBpJSkqybQWb18Pk0XD0EDzyDAwZW2UunINHChk77TjvLD3N5Rf5Mu/5NsRcGyDB3sFq7OFr\nrdcrpYyWigCByvybDAAyALndLkQDs2u2yIJ8+OhdWP81dLgIRj4NLVsDlXv2RcWaD/93hqRvsmji\n58FzD4XQt6e/BHonYYsx/LeBz4EjQCBwr9a6hqxIQrivhpwXn5ubS1xcHHPmzKnzeY0+ihdCvYkw\neJCUWcy8b7fx/TOtzyu3O7mAaYsyOHikiBuu8Gf0Xc1pFuhZ53qF7dki4PcHtgHXAx2B75RSP2it\nz5xbUCkVD8QDhIVZzn8thKgde2SLvKWpJ2NbeJFrggmHC9mcd35fLr/QxPyVWSxdfZaQIE9eHtWS\nqO5+da5T2I9V8/BLh3S+qOam7ZfAFK31D6WPVwNPa603WTqnzMMXwrYsZYtMTk6u3clyc2DhDPOY\nfdfLIW4CBAWfV+y3veZkZ0dPFTPwHwHE396MJn4y29ueHL3SNhXoB/yglAoFugBVr48WQtiNzbJF\n/rUHEl+BjJNw50PQ/y7wqBzEs/NMvPtZJl/+lEPbll5MH9+KHp1lI3tnZ820zMWYZ9+0UEqlAS8A\n3gBa69nA/wELlFLbAQVM0lqfsluLhRBVqne2SJMJ/rcUli+AZiEwaRp07HpesZ/+yOXNxZlknilh\n8I2BDLslCIOP9OpdgaRWEKKRqdPCq6xMmDcNdm6BntfCsPHgH1CpSObZEt7+OJM1W3LpcIE3E4cE\n0yXcYLuGC6s4ekhHCOHKdm2FuVMhLweGjIXrboYK0yi11nz/ay5vf5JJXoGJ4QODGHxjU7y9ZKql\nq5GAL0QjY3XPvrgYPv8Avv4YWreHJ16GdhGVipzIKGb64gw27syna4QPEx4IwdhG0iK4Kgn4Qrij\n9OOQOAUO7IbrboJ7HwbD3zddTSbNyh+zmbP8NCYTjL6rGYOiA/H0kF69K5OAL4S72fIjvP8maBPE\n/wuu7FN+KDo6GnxaEzlgBn/sL6DnRb48cX8wbVpIqGgM5LcohLsoLICP58DaL8DYGR7+F7RsU364\npERDyM2olnfw1+FCJg4JZkCUJDtrTCTgC+EOjqTCuy/D4WTzvPpBw8Dr77H4A2mFvPZBOh6hg9Fn\nNjP/1dsICZK0CI2NBHwhGjOt4adv4cNZ5jH68S/BJX/P6Css0nzwdRZLvj3D6ZQVbP/+WQryztLz\nx1rO4RcuQQK+EI1VXg588F/YtBYuioQRE80Lqkrt/KuAqR+kk3q8mFamr/np+6cpyDOv0rVppk3h\nNGThlRAuxppsm10Miudbe9PaSzEvo5jFmSWUpz1TBlSruyD4RijKQB+dz4bv51SZZM1gMBAVFWXT\n9oOdd+Nq5GThlRACMOc2uauZJw+HeHGqGB47XMiO/AqduibdUG0eQvm0RGd8hz7xCZjy7ZJpUzgf\n6eEL0VicPQ3vvQ47foXLroZhj0NAoPlQrol3lmbyzS85tA/1YkJsMN0v/HvevU0zbQq7kh6+EO5u\n9zaY+xrknIXYMRB9S3l6hB+25TJjSQans03c378pQ28Owse78lRLm2XaFE5NAr4QrqykBFYugi+X\nQGg78yyc9h0AyMgq4a2PM1j/Wx4XtvPm5dGt6BxW9XbT9c60KVyCBHwhXFX6CZgzBfbvgmv/CfeN\nBoMvWmu+3ZjDrE9Pk19oYsStQdxzY1O8PC0voIqNjS3fClFuqjZOEvCFcEW//Qzz3zD38EdOgt59\nATiWbk529uuufC7paGBCbDBhra1PdiaBvnGTgC+EKykqhE/mwurPIbyTORdO6AWYTJoV67OZs+I0\nAGPvac5t1wXgIcnORAUS8IVwFccOwbtT4NABuHGQeftBL29SjxUxLSmDHQcKuKKrL4/fF0zrEHlr\ni/PJX4UQzk5r+HkVfDgTvH1g3H/g0t4Ul2g+/iaL97/KwtfHg0lDg/lnb0l2JqonAV8IZ5afC4ve\nhg2roXN383h98xb8eaiQqR+ksz+tiOsu82PcPcEES7IzUQMJ+EI4q+Q/IfEVOHkMbhsCtwymsMSD\n95ef5qNVZ2gW4MGLI1tw3WX+jm6pcBES8IVwNlrDquXw6XvQtBlMfA06X8L2/flMXZRB2oliBlzV\nhFF3NifQ38PRrRUuRAK+EM7kbJZ5uuUfGyEyCh58glyvAOZ8lMGKddmEBnvy2tiW9LrYz9EtFS5I\nAr4QzmLvHzDnVcg+A/ePhr4D2bgrn+kfHuXk6RLu6BtI3MAg/HylVy/qRgK+EI5mKoGVH8IXi6FV\nGxg3maxgI7PeT+e7TbmEt/birSdD6dbB4OiWChcnAV8IR8o4aU56tm87XH0D+r7RrNsN/515lDM5\nJh64qSkPDDg/2ZkQdSEBX9hE2aYcsjS/FrZtgPmvQ3ERxE0gvWtf3lyYwU+/59E5zIfXxgbTsV3V\nyc6EqAsJ+ELUU60/7IoKYek880ycsI7okU/zzcHmzJp8hKJiiB/UjLuvD8SzhmRnQtRWjQFfKTUP\niAFOaK0vqaZMNPAm4A2c0lr3sWUjhWg0jqWZ59anHoAbbudIn6FM/ySbLXsyuPRCAxMeCKZdK+uT\nnQlRG9b08BcAbwMLqzqolGoGzAIGaK1TlVKtbNc84QqSkpLYsGEDBQUFGI1GyaNenV9WmVfNenlT\nMvoFlp3uxrxX0/HwgPGDmxNzrSQ7E/ZVY8DXWq9XShktFLkf+ExrnVpa/oRtmiZcQVJSEvHx8eV7\nn6akpBAfHw/gFkHfqg+7/DxImmkO+J27k3zrBKat1Ow6eJre3czJzloFy+iqsD9b/JV1BryVUmuB\nQGCG1rrKbwPCscrGmm2pLNhVlJubS1xcXPlmGrbkTDeFrfqwS90P774CJ45SdMsDLPG5hUWzzuJn\n8OCZB0Pod4W/JDsTDcYWAd8L6An0A/yAX5RSG7TW+84tqJSKB+IBwsLCbFC1cLRzg31NzzuCPT7o\noOYPu0FBnoxq4UVWCTyb2ZY/VxpRvmfRWRsoOPYBL/16lpesrMuZPuiE67JFwE8D0rXWOUCOUmo9\n0AM4L+BrrROBRIBevXppG9QtasEeQcNoNJKSknLe8+Hh4Y0+SFn6sHuptTfXBniyPseDBNNACi+4\nCYpPY0qdDtm/NXBLhTCzRcBfAbytlPICfIDewHQbnFe4gISEBOLj48nNzS1/zt/fn4SEBAe2qjJ7\nfPBER0djMBiqDPrhTf25NsjA733HM2d3V4pOFhNzTRMeviOMAD95awjHsWZa5mIgGmihlEoDXsA8\n/RKt9Wyt9W6l1DfAH4AJmKu13mG/JgtnUjZWHRcXR0FBAeHh4W4zSyciIoLU1NTKH3aeHjzf5yqm\n9/wvK9d70qYFTHusFZd38XVgS4UwU1o7ZmSlV69eevPmzQ6pW9ieu620LbvekSNHMjIujvyCAsKa\n+DL81gfZZXyW9CwTd14fyPCBQfj6SLIzYTtKqS1a6151ea3MBRM24S6BHipPxUzet5e3ruxEzIWd\nmXfRi6xKDsLo58mLI1tycYQkOxPORQK+ELVw3lTMo8cYczKd2YEjCW4SxLBbgri/f1O8vWSqpXA+\nMqQjGqWGnIoJYGjSiqhrboSCNJvV5U7fmoT16jOkI4OLQtRCtVMxc07YNNgLYQ8ypCMaJZv3jgvy\n4cNZtNu8mcM5Oecddod1B8L1SQ9fuI3o6Oi6DfUc+ouSyWP5ZIMm5OqX8PSqvJ+ss607EKI60sMX\nojpaw9ovOPjRSqZ6DmOPZxgxg/wYdXcI48eOdLt1B8L1ScAXbqHWKZyzz1I0fwYfbm9Gkte/CPD3\n4t/3BhPd0x+lhrAk6T1AbqwK1yIBXzR6tU7h/OdOds9exNTc20j2uoB+vfwYc08wQQGe5UUk0AtX\nJNMyhdNo8KmUBgNRUVHljz2Ae4L9KGlxJ5953oBnSSaFRxZA9u9W1yUfBMLeZKWtEBZYk8I52BPu\nbdudz5s8yFHVEs/M7yk8/hGY8huqmULYnfTwRaNnKYVzcnIy2Zs3M3vBQb4y9aZtQCETRrSlR2e/\nKs4khOPJwishLEhISMDf37/Sc/7+/iT832R+mrWc4e/58I3pCgZfBXNf6iDBXjRaMqQjGr2qUjg/\nM34Cf21rzXt5F9GhSRYvPdKCLhcGOLilQtiXBHzhFmJjY8v32H324WnM/KEJeRgYfvlp7hveHS9P\nSXYmGj8J+MJt+HgH0z5sOK/80IquPkeYEB+OsWsHRzdLiAYjAV80eiaTZuXKFHw6v8pRE4zpvJfb\nRkfj6ePt6KYJ0aAk4ItG7dCxQl5/Zz9/nAygp0cqTwwLoU3UjY5ulhAOIQFfNEolJZpPvjnFgq/O\n4mNSTAxdw4DxA1HNgh3dNCEcRgK+aHQOpBXy2tzD/HnCk2v1dh77pybktiHgIbOQhXuTgC8ajcIi\nzaKvT7P4f1k0NeXwgv9XXDf6VlTHi+xSn7tt3C5cnwR80Sjs/KuAaQtPkHJCc2PJRkZ3TyXoodHg\nL3PrhSgjAV+4tLx8E++tzGLZmjO05DRT9IdceX8fuG4CKJlbL0RFEvCFy9q8O483kjI4llHC7SVr\nGNFyC/6jnoS2RpvXde7wTa3z6wvhBCTgC5dzNtfEO0sz+eaXHNp7ZTCjaC7d/9EJ7p0KBl+711/r\n/PpCOAnJlilcyg/bcpmxJIPTZ0sYzCqGqv/hM2wMXHHdeWVtmV9/27ZtAERGRlqdX98W5IawOJfk\nwxeNXkZWCW99nMH63/K40P80Lxe+TWejH8TPgJZt7Fr38ePHOXPmDFrraoM9VJ93XwhnUWPAV0rN\nA2KAE1rrSyyUuwL4BRistf7Udk0U7kxrzbcbc5j16WnyC0oY0WQt92R+jNeAO+D2YeBV/Z+wLXrH\nZcM3Zd+ECwoKUEpR1Tfj8PBw6ZELp2ZND38B8DawsLoCSilP4FXgW9s0Swg4ll7M9MUZ/Lorn0ta\n5DAhdzphxWdg/H/gkvO/0dpji8SqevRVBXsPDw/8/Pzq3Qb5wBD2VGPA11qvV0oZayg2FlgKXGGD\nNgk3ZzJpVqzPZs6K0yg049r+yq0H5+JxcSSMmAJBDZcewdIwTVlP32AwEBERQWhoaIO1S4i6qPcY\nvlKqLTAI6IsEfJdX0+pRe68uTT1WxLSkDHYcKOAKYwmPn3qT1il/wh0PwoC7LaZHsEebLG2PaDQa\n7VavEPZgi5u2bwKTtNYmVcNCF6VUPBAPEBYWZoOqRWNRXKL56LszLPwqC18fxaTIA/xzyzRUsxCY\nNA06dnVIuxISEoiPjyc3N7f8OX9/f5l3L1yT1rrGH8AI7Kjm2EEgufQnGzgB3F7TOXv27KmFc1m0\naJE2GAwiPwENAAAWQ0lEQVQa0OHh4XrRokW1Ol5Xe1MK9MiEI7rvqBT94qzDOn3qS1rH9dd61v9p\nnX3GJnXUh72uW4i6ADZrK+J2VT/1DvjnlFsA3GXNOSXg102fPn10nz59bH7eRYsWaX9/fw2U//j7\n+5cHt5qO10V+QYlOXJap+z2aou+cdEivX7lD6ycGa/1wjNarV2ptMtnq8urNXv/uQtRWfQJ+jQuv\nlFKLgWigBXAceAHwLv12MPucsguAL7QV0zJl4VXdREdHs23bNiIjI6stU7ZIqDbK5pmfSylF06ZN\nazxeW5FX3YO6IA5laAOZ67m/8GMealZAapFm8rEi5n27ptbnFMId2HXhldb6PmtPprV+sC6NENap\nmL9lw4YNNp0ZUt0Hf9nzNR23loeXPxdEPoZHxH3owpMEpb3K5MA/6d7cgy+zSvjvqWLyHbP4W4hG\nT1bauohz87cUFBSQmprKc889Z5Obh5ZmoyQnJ9d43BqbdubxxuIMTmaWMCg6kLj2h/D78DiYAmDo\nOG65Mppb6nshQohqScC3MXss/oGqFwDl5uYSFxfHnDlz6n1+Pz8/PDw8MJlM5c+VzUYBy7NVapKV\nXcI7S0/z7cYcwlt78da45nT79X2Y+wUYO0P809DqgnpfgxDCMgn4LsLe+VvKhob27t2L1prw8PBK\nUw/L/hsXF0dBQcF5x6uitWbdb3n896MMzuSYeOCmpjzQ4ww+702Ew8nQ/04Y9CB4edvkGoQQlkm2\nTBdhiyEVa9hq4VV6VglvLsngp9/z6Bzmw8TY5nRMWQsfzgIfX4ibAN1lnZ4QtSXZMt1AfYZUaqO+\nq0a11nzzSw7vLM2ksBjiBzXj7qs88PxwOmxaCxf1gBFPQbMQm7RXCGE9Cfguoi5DKvZg6QPh6Kli\n3vgwgy178rn0QgNPxgbTPu8gJLwCp47D7UPh5nvBw7PhGiyEKCcB34XExsaW36C1Vf4WWwzhlJg0\ny9ee5b3Ps/DwgPGDmxNztT8eq5fD0vkQ1Byeeg06VZtdWwjRACTguxhnS9SVfLSIaYvS2XWwkN7d\nfHn8vmBaeWfDrP/AH5vgsqth2OMQEOjopgrh9iTgi2pZ2qi7qFiz5NszLPomCz+DB888GEK/K/xR\ne/+AOa9Czlm4fzT0HQg1JNUTQjQMCfhuzFJAt7RRd69r72Lqogz+OlxE317+jLm7Oc39gRUfwJeL\nIbQtjH8J2ndw1KUJIaog0zIbgbos9jp+/Dj79u2rtNDKw8ODzp07ExoaWv1G3f7BXDV0CxSfRh99\nH7J/o5UXPBfqw6V+Hnx9ppgZJ2ufHsHZhqqEcFYyLbORqymg2ypZmslkYu/evRw9erT6hV65maQf\nWMbh396gpOgsN7dpzswrOuKlNPGb9vNx6qkqX2cp2ZsQomFIwG8E6hJM161bV+XzWmsiIyOr7+H7\nNiG4cAWtL+nIIyFe3NHMi735JiYfL+JwcDsig9tVeV7pwQvheBLwXUBDb923du1akpKSGDFiJPn5\neeXH/P39SUycTWy/6+DdKXDoANwwiC53DifJ28fmbRRC2Fb1G4SKRi0hIQF/f/9Kz5Wt3M3KLiGl\nsD8R10zB0MScYyc8PJzExERiI0Lh/8ZC5kkY+x8Y/DBIsBfCJUjAd1LR0dF2y7wJ5kVciYmJGAwG\n4O+A3qbL7Tw4+Shrt+Yy6fFhRPXqRp8+fUjes4vYvMMwbxqEd4IXZkGP3k5xLUII68iQjhuruHL3\nk+WreHNxJu/NS+ciow8THwgm4gIfht3yPaT8CZPHwMljcNsQuGWwpEcQwgVJwHdza9as4cufcnho\n8lGKS2DUnc24o28gnh4KtIbvV8Anc6FpM5j4KnTuXqvzW5rrL4RoWBLwnVBDBcnDJ4p4/cMMtu0r\nILKzOdlZ25aluenPZsGCN+D3jRAZBQ8+AQG127vW0uItCfpCNDxZeOVkyoLkuWmQExMTbRYkS0ya\npavPMn9lFl6e8MidzZn6bEz58R6+iuda+xDkCe+eKmZpVkmd6ql2aqfBQFRUVJ3bXx2Z+incgSy8\ncgBX3coQQ1vUBSNQfh25qrsf4+9rTstmXkzFfAd/SHNPhgZ7caRI868jRewvrHuHwN67dAkhakcC\nvpOxW5BUXqgWA6HFQCjJxZQ2k5dmTkGVJjZb+9knMPc12LcdrrqB9rGjmevrX8NJLatprr8QomFJ\nwK8jWwessm8M4eHhNg+Suw8WMHVRBslHi7jhCn8evTucoIBX/y7w+0aY/zoUFcJDE+DqG+pUz7ka\napcuIYR1JOA7GVsGybwCE/NXZrF0zVlaBHny8qiWRHX3+7tAUSEsnQerlkP7jvDw09C6fb2v4dxN\nUxy9S5cQwkwCvhOoOCsnOTmZYcOGMW/evHoFya1783k9KYOjp4oZeG0A8YOa0cSvwjq744fh3Zch\n9QD0uw3uirPLill77NIlhKgbCfgOVtXUxffff5+IiAhCQ0NrHSSzc03MXpbJVz/l0LalF9PHt6JH\nZ9/KhX75Hha9DV5eMOYFiLzKRlcjhHBmEvBrwR4zc6qblXPw4EFCQ0NrV2fA5ag2w8ArCNK/4tDu\nZTy2rqj8sJ+Cr4fEwM+rzPvLjpwEwS1tdCVmVa0hkJ69EM6hxoCvlJoHxAAntNbn7UKtlIoFJgEK\nOAuM0lr/buuGNlY2mZXjGYhqPQQVFIXOT0UfehPyD1YqcqGP4vnW3ube/cBYiLkfPG2bHkEWWgnh\n3GpceKWUug7IBhZWE/CvBnZrrTOVUjcBL2qta8yqJQuvzCxNXUxOTrb4Wq01qzblMvPTTPIKTDxw\nUxCDb2yKt5eqWAjWrISP55hXyo54Ci7q0WDfVkAWWglhS3ZdeKW1Xq+UMlo4/nOFhxuAqnfAEFWy\nNCvn3NkuFZ3IKGb64gw27syna4QPEx4IwdjGu3Kh7LPm9AjbfoFLr4ThT0JgkN2uxZZrCMp28ZKd\nsoSwHVuP4ccBX9v4nI1a2VBHVVMXq1pZazJpvvgxm8TlpzGZ4NG7mnF7dGmys4r+3AFzXoWsTLgn\nHm4cBOrvMo7YVKU2LH3YCSHqxmYBXynVF3PAv9ZCmXggHiAsLMxWVbs8a6cupp0oYtqiDP7YX0DP\ni3x54v5g2rQ451doKoEvl8DnSdAiFP41HYyd7Nj6v9lqDYFk2BTCTrTWNf4ARmCHheOXAgeAztac\nT2tNz549tajeokWLtMFg0IAOCwvTj06ao/uPS9XR8bt0n9v+rU0m0/kvyjyl9WsTtY7rr3XiFK1z\nsx3a7vDwcL1o0aJav97f318D5T/+/v61Po8QjRWwWVsZZ8/9sSpbZukY/he66pu2YcBqYKiuPJ5v\nkdy0rV5VGTM9vPy4dcgbZKZ+A8Wnz/8m8McmmPc6FOZD7KNw9Y2VhnBszdJN3/qMvzfkjV8ZLhKu\nyK43bZVSi4FooIVSKg14AfAG0FrPBp4HQoBZpYm4iuvaGFfUULNdTMV5rFw4BpPJhNYaX19fIiIi\naBsaysgQL+5t7sX+AhOTjxWR+uzLwMtW12frwFefG62SYVMI+7Fmls59NRwfAYywWYtEtcGtpKSk\nUpk/9+3joWAv7m3elmWni3knvZh6ZDOuFXv1jiXDphD2Iytt66mqIFTnXr8yoFrdhSFgPwXZh2ss\nXmIy8d/f95PVrCU/5pjqVidVt9dRwVUybAphPx41FxENokk3VMeXUSH96XDJTXh4WPerOZSTX69g\n72xiY2NJTEzEYDAA5p69LXf7EsKdyRaHdmTNXPKzuSZmL83k619yaB/qxYTYYLpf6EtSUlKlufnZ\n2dmkp6ef93prVuS6IpmHL0TV6nPTVnr4DvTDtlyGTz7C/zbmcH//psx5pg3dLzRntoyNjSUqKoo+\nffqQfPAgMx4ehr9X5dw3tRnqiI6Ottu2jEII1yBj+A6QcaaEtz7KYP1veVzYzptXHm1Fp/bn56Jf\nu3Yt5JyFd14i9vhOuDeGkZ9+Q54bbCYiPXshbE8CfhVsMZxQ1WrR+++/n2835jDr09PkF5owHf+Y\nfbu+ptMz31d9kv27IHEKZKXDXXHE/vNO5qRdX+u2ycpVIQRIwLeLqtIEjxwZz5Jvz5DT5BYu6Whg\nQmwwQwd/UfUJTCXw9SewYiGEhMKkN6BDF5u1RVIWC+Ge3PqmbVVj2sePH2fv3r1orTEYDOU7T9VG\ntatFAy4g6oZHIHM1x48fq1SPh4cHPj4+tDJ4k3jlhUSHBrH00Cke33KQM8UlVdRiVtMiJ1m5KkTj\nYteVtu7k+PHj7Nu3ryw/EAUFBezbtw+gVkG/2tWi2Ucg8/sq6wHoFuTHl/26E+DlydjNB/gg+WR9\nLsdyW2TlqhBuxyUDvr1mmxw8eBCTqfKcdpPJxN69ezl69KjV51HKA63PnxuvlGLbtm2cOXOGqr5Z\nnTyby1nlxcTDRaQ0a0tkZNvaX8Q5DAZDlcFdVq4K4X5kWmYF1fV6azPs5df8Yjr0noCHl995x3x9\nfS2eLzUnn1FphaQU2W6YLSIiAn9//0rPycpVIdxUXdNs1vfHGdMjh4eHV0rLW/YTHh5e42vzC0p0\n4rJM3e/RFH3npEP63wnzqk0TXJ966qK+KYuFEM6DeqRHdskhHXupax6X7fvzmboog7QTxdx0VRMe\nubM5gf7DWf/t+8D5NzMT/vOiuZ7ColrVU1fWbrAihGjcJOBXYGm7wark5puYs+I0K9Zl0zrEk9fG\ntqTXxX8P5VQZXA8nE3vgZ7j6Yh7buI/0vPwGWUQlgV4I4dbTMqtjzcKrTTvzeOPDDE6eLmFQdCBx\nA4Pw87VwS0RrWP8VLHkX/JpA3ESiH32ixnqEEKIimZZpY5YCcFZ2CbM+zeS7TbmEt/birSdD6dbB\nYPmEudmwcAZs/gG6XQ4PTYSg5rZttBBC1EACvpW01qz7LY//fpTBmRwTQ25qSuyAIHy8a9hG8MBu\nmDMFMk/BnQ9B/7ugNPWx9OyFEA1JAr4V0rNKmLEkgx9/z6NzmA+vjQ2mY7vzk51VYjLB/z6F5e9D\n8xYw6XXocFHDNFgIIaogAd8CrTXf/JLDrKWZFBVD/KBm3H19IJ6eNfTqszLhvamwaytrzpbw+p8p\nfDFFgr0QwrEk4Ffj6Kli3vgwgy178rn0QgNPxgbTPtS75hu6O7eag31eDgwZx39eeK3B2iyEEJbI\nStsqHDpeRNxLR9mdXMD4wc15Y3wr2od6W35RcTEsnQfTn4GAQHjuLZLSMtmwYQPr1q3DaDSSlJTU\nMBcghBBVkB5+Fdq18uK+/k3p37sJrYL//ieqNq/8qWPmvPV/7YHrboJ7Hybp06WSllgI4VTcfh6+\ntYnYyjJcVkyu5uHhwa3du/DBZe0BmHayiLXZ5uMNmZYYZMaPEO5C5uGfozbZNLdt22ZVuaoyXJpM\nJn7bf5A9HYKJ2/AnKbl/B3hLaYmtqbOmPPdCCFFbjTLg14a1gXXdunVVPp+Sk8/TGV4073wxFZdS\nWerhW1On9NiFELbWKAO+PYKl0WgkJSXlvOfDw8NZVUV9ZVsLnpuILTExUcbwhRAOIbN0rJTwwvP4\n+1SeqWMpw2VsbCyJiYkYDOa0C+Hh4RLshRAOVWPAV0rNU0qdUErtqOa4Ukq9pZTar5T6Qyl1ue2b\n6WB/7SV273oSr7mYYD/rA3hsbCxRUVH06dOH5ORkCfZCCIeyZkhnAfA2sLCa4zcBnUp/egPvlP7X\n9ZlM8N1n8Nl8CArhZ59gSnx86XNllIyxCyFcTo0BX2u9XilltFDkNmBh6U4sG5RSzZRSbbTW1m8C\n64zOnIZ502DHZrj8Ghg2np23DCQyMrJWwV4+GIQQzsIWN23bAocqPE4rfc51A/6JI/DqBMg5C7Fj\nIPoWUDXkzxFCCCfXoDdtlVLxSqnNSqnNJ0+ebMiqa6dFKFx6JTw7A/rGgFLlq2wlTYIQwlXZIuAf\nBtpXeNyu9LnzaK0Ttda9tNa9WrZsaYOq7cTDE4aNh/YdgL+nWJ6bJkGCvhDCldhiSOdzYIxSagnm\nm7VZjhy/r80qW2tVtYgqNzeXuLi48s3BbUnG/YUQ9lBjwFdKLQaigRZKqTTgBcAbQGs9G/gKuBnY\nD+QCw+3VWEexlCZBCCFchdsnT7OGpVW2ycnJDd8gIYTbqk/yNFlpa4WEhAT8/f0rPWdpla0QQjgj\nCfhWkDQJQojGoFEmT7OH2NjY8hu0clNVCOGKJODXggR6IYQrkyEdIYRwExLwhRDCTUjAF0IINyEB\nXwgh3IQEfCGEcBMS8IUQwk1IwBdCCDchAV8IIdyEBHwhhHATDsuWqZQ6CZyfgtK5tABOOboR9STX\n4BzkGpxDY7iGLlrrwLq80GGpFbTWTrzllZlSanNd05A6C7kG5yDX4BwayzXU9bUypCOEEG5CAr4Q\nQrgJCfiWJTq6ATYg1+Ac5Bqcg1tfg8Nu2gohhGhY0sMXQgg34fYBXyk1QCm1Vym1Xyn1dBXHlVLq\nrdLjfyilLndEOy2x4hpiS9u+XSn1s1KqhyPaaUlN11Ch3BVKqWKl1F0N2T5rWXMdSqlopdQ2pdRO\npdS6hm5jTaz4ewpSSq1USv1eeg3DHdHO6iil5imlTiildlRz3BXe0zVdQ93e01prt/0BPIEDQAfA\nB/gd6HpOmZuBrwEFRAEbHd3uOlzD1UDz0v+/yRWvoUK51cBXwF2ObncdfxfNgF1AWOnjVo5udx2u\n4Rng1dL/bwlkAD6ObnuF9l0HXA7sqOa4U7+nrbyGOr2n3b2HfyWwX2v9l9a6EFgC3HZOmduAhdps\nA9BMKdWmoRtqQY3XoLX+WWudWfpwA9CugdtYE2t+DwBjgaXAiYZsXC1Ycx33A59prVMBtNbOdi3W\nXIMGApVSCgjAHPCLG7aZ1dNar8fcpuo4+3u6xmuo63va3QN+W+BQhcdppc/Vtowj1bZ9cZh7N86k\nxmtQSrUFBgHvNGC7asua30VnoLlSaq1SaotSamiDtc461lzD28DFwBFgO/CY1trUMM2zCWd/T9eW\n1e9p2cTcjSil+mL+47jW0W2pgzeBSVprk7lj6bK8gJ5AP8AP+EUptUFrvc+xzaqV/sA24HqgI/Cd\nUuoHrfUZxzbL/dT2Pe3uAf8w0L7C43alz9W2jCNZ1T6l1KXAXOAmrXV6A7XNWtZcQy9gSWmwbwHc\nrJQq1lovb5gmWsWa60gD0rXWOUCOUmo90ANwloBvzTUMB6Zo8wDyfqXUQeAiYFPDNLHenP09bZW6\nvKfdfUjnV6CTUipCKeUDDAY+P6fM58DQ0jv7UUCW1vpoQzfUghqvQSkVBnwGDHHSnmSN16C1jtBa\nG7XWRuBTYLSTBXuw7u9pBXCtUspLKeUP9AZ2N3A7LbHmGlIxf0NBKRUKdAH+atBW1o+zv6drVNf3\ntFv38LXWxUqpMcD/MM9OmKe13qmUeqT0+GzMM0JuBvYDuZh7N07Dymt4HggBZpX2kIu1EyWQsvIa\nnJ4116G13q2U+gb4AzABc7XWVU69cwQrfxf/ByxQSm3HPNNlktbaaTJQKqUWA9FAC6VUGvAC4A2u\n8Z4Gq66hTu9pWWkrhBBuwt2HdIQQwm1IwBdCCDchAV8IIdyEBHwhhHATEvCFEMJNSMAXQgg3IQFf\nCCHchAR8IYRwE/8Pn9dX2gavZWAAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x1818d87310>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%pylab inline\n", | |
| "ptrue = [1,1]\n", | |
| "n = 30\n", | |
| "xtrue = np.random.random(n) \n", | |
| "ytrue = polyval(ptrue,xtrue)\n", | |
| "xi = linspace(0,1)\n", | |
| "\n", | |
| "xerr = np.ones_like(xtrue) * 0.1\n", | |
| "yerr = np.ones_like(ytrue) * 0.03 \n", | |
| "xobs = xtrue + xerr * np.random.randn(n)\n", | |
| "yobs = ytrue + yerr * np.random.randn(n)\n", | |
| "\n", | |
| "pfit = np.polyfit(xobs,yobs,1)\n", | |
| "\n", | |
| "errorbar(xobs,yobs,yerr=yerr,xerr=xerr,fmt='o')\n", | |
| "plot(xi,polyval(ptrue,xi),label='True Function')\n", | |
| "plot(xi,polyval(pfit,xi),label='Fitted Function')\n", | |
| "\n", | |
| "print \"RMS(xtrue - xobs) = {}\".format(std(xobs - xtrue))\n", | |
| "legend()\n", | |
| "\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### correct x coordinate \n", | |
| "\n", | |
| "- RMS between true and corrected is less than what we put in, because the high precision y coordinate puts the point onto the fitted relationship, and the fitted relationship incorporates many data points. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 76, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "RMS(xcorr - xobs) = 0.0536557067867\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4VNW9xvHvIoRAIkJUFEGYUGuViooahQLVgFW8IGC5\nFAy2WIRysVWriFJrayvgjeppAT0JIKIRFfEQrAKKFTAgYMSoiBpEFChI4oVETAxJZp0/diK5TJhJ\nsjO3vJ/n4SEze8/sn5h5s7L22r9trLWIiEh0aRHqAkRExH0KdxGRKKRwFxGJQgp3EZEopHAXEYlC\nCncRkSikcBcRiUIKdxGRKKRwFxGJQi1DdeATTjjBJiUlherwIiIR6e233/7SWtvB334hC/ekpCSy\ns7NDdXgRkYhkjPk8kP00LSMiEoUU7iIiUUjhLiIShRTuIiJRSOEuIhKFFO4iIlFI4S4iEoUU7iIi\nUUjhLiLiS0qK8ydCKdxFRKKQwl1EJAop3EVEasrIgE2bYN06SEpyHkcYhbuISFUZGTBhApSUOI8/\n/9x5HGEBH7KukCIirnD7pOemTUeCvVJREYwbB+np7hxj7Vp33uco/I7cjTFdjDGvG2O2G2M+MMbc\n5GMfY4z5pzHmE2PMe8aY85qmXBGRJlYz2P09H6YCGbmXAbdaa7caY9oCbxtjXrXWbq+yzxXAaRV/\negGPVvwtItK03B4FJyU5UzE1eTxBGXG7xe/I3Vq731q7teLrb4EPgc41dhsCLLaOTUB7Y8zJrlcr\nIuJDaWkpOTk57rzZjBkQH1/9ufh45/kIUq8TqsaYJOBcYHONTZ2BPVUe76X2DwAREdfl5+dz6aWX\nctFFF5GXl9f4N0xNhbQ0iItzHns8zuPU1Ma/dxAFfELVGHMMsAy42Vpb2JCDGWMmABMAunbt2pC3\nEBH5QU5ODkOHDuWLL75g/vz5nHjiie68cWrqkZOnETQVU1VA4W6MicUJ9gxr7Qs+dvkv0KXK41Mq\nnqvGWpsGpAEkJyfbelcrIlJh6dKljB07lsTERLKyskhOTnb3ABEa6pUCWS1jgAXAh9baf9Sx2wrg\n1xWrZnoDBdba/S7WKSICgNfr5a677mLkyJH07NmT7Oxs94M9CgQycu8LXAe8b4ypPGMxHegKYK19\nDHgZuBL4BCgCrne/VBFp7goLCxkzZgwvvvgi48aNY+7cucRVzo1LNX7D3VqbBRg/+1hgiltFiYjU\ntGPHDoYMGUJubi7/+te/mDJlCs7EgviiK1RFJOytXr2aUaNGERMTw5o1a0iJ4Fa8waLeMiIStqy1\nPPjgg1x55ZV07dqV7OxsBXuAFO4iEpaKi4u57rrruP322/nlL3/Jxo0bSUpKCnVZEUPhLiJhZ+/e\nvfz85z8nIyODe++9l+eee46EhIRQlxVRNOcuImFlw4YNDBs2jKKiIjIzMxk8eLD/F+3fA+9thoHD\nm77ACKGRu4iEjfnz59O/f3/atm3Lpk2b/Ae7tbB+Jfz9Rli5FL49GJxCI4BG7iIScqWlpdxyyy3M\nnTuXyy67jGeeeYbExMSjv+jQt7D4Edi6AbqfC+Nug7btg1NwBFC4i0hI5efnM3LkSNauXcttt93G\nrFmzaNnSTzR99C4seBAKD8KIG+DSX0ILTURUpXAXkZB59913GTJkCF988QVPPvkkY8aMOfoLyspg\nxVOw8lk4sRPc+TAknRacYiOMwl1EQqLejb/y9kH6/bDrY+g3EEZNhNZtglNsBFK4i0hQeb1e7r77\nbmbMmEGfPn1YtmwZHTt2rPsF1sKm1+CpuRATAxOnQ/JFwSs4QincRSRo6t34q+g7yJgDm1+H03rA\nDbfD8S71bI9yCncRCYqqjb/mzJnD5MmTjzT+qmwpULWH+s7tzjTM1/kw5Dq4ahS0iAl22RFL4S4i\nTa5ejb+85fDSM/BiBhzXAaY9BKf+NGi1RguFu4g0GWsts2fPZtq0afTo0YPMzMza/WEyMmDTJigp\ngS5doO/ZcEw59OoPqTdCvNoONIQWhopIk6hs/DV16lSGDRvmu/FXRgZMmOAEO8DevbBsNZxyHoyf\npmBvBI3cRaR+Ami5u7ekhKHbtrH10CFmJCVx54EDmKuuqr1j5Yi9qrJyuO8fsHajO/VG+L1QG0rh\nLiKu2lBQwLAPPqDI62X5mWcy+IQT6t65ZrD7e14CpnAXkfo5ykh4/vz5TJ48GY/Hw+srVtC9e3ff\nK2G8XnhlGby9BQ4V134jj6fZjrjdojl3EWm00tJSbrzxRsaPH8+AAQPYsmWLE+y+HPwKHp4Ozy+A\n4VdBmxpXmcbHw4wZTV90lNPIXUQapWrjr6lTpzJr1ixiYirWo1ddCZOUBBPGQt52OFwCv74Jfn45\n/OJpGDfO2cfjcYI9NTWU/0lRwVhrQ3Lg5ORkm52dHZJji4g7cnJyGDp0KAcOHGD+/PmkVg3lypUw\nRUVHnotpAUMGwJyFcHKXI8/7mroRn4wxb1tr/TTiUbiLSFX1uPn0c3l5XP/xxxwXG8v/nXkmyW3b\nVt/B10oYgLg46N27cXVWaoY/DAINd825i0i9eK3lT7t28asPP6TnMcfw1nnn1Q520EqYENOcu4gc\n4WckXFBQwJgxY/j37t1Hb/xVeNCZY/+moPY2rYQJCo3cRSQgubm59O7dm5UrVzJnzhzS09N9B/u2\nbPjrJDi3mzMFU5VWwgSNwl1E/Fq1ahUXXngh+fn5rFmzhilLl2L696++U+lheDYNHrkLjjkWnngB\nFiw4EvAeD6SlaSVMkPidljHGLAQGAXnW2h4+trcDngK6VrzfQ9bax90uVESCr2bjr+XLl9OtW7fa\nO+7fA2n3wZ6dMGAwDB8HreIgtRukpzv7aComqAKZc18EzAEW17F9CrDdWnu1MaYD8LExJsNae9il\nGkUkBIqLixk/fjwZGRmMGDGCxx9/nISEhOpr1z0e+PWv4KuPoVVr+P09cE6vUJcuBBDu1tr1xpik\no+0CtDVO1/1jgK+BMleqE5GQ2LNnD9dccw1bt25lxowZ3Hnnnc6NNWp2cdy9G2bNhl8NgrnzoP3x\ntd9MI/aQcGO1zBxgBbAPaAv8ylrrdeF9RcQN9Vi7Dk7jr19+8AHFXi+ZP/0pV7/yCrzyirPR19r1\ncq/Tpve/w9ypVz8MXOHGCdWBQA7QCegJzDHGHOtrR2PMBGNMtjEmOz8/34VDi4ib5u/fT/933+XY\nli3ZdO65XF2zo6PWrkcMN0bu1wP3WedS10+MMbuAM4AtNXe01qYBaeBcoerCsUXEnwBGwqWlpdxy\nyy3MXbeOgQMHsmTJEhITE6vvlLcPfnIGFHxb+w20dj3suDFy3w1cAmCMOQk4HfjUhfcVkSDIz8/n\n0ksvZe7cuUydOpWXXnqperBbC2+ugXumQK8zoLXWrkeCQJZCLgFSgBOMMXuBvwCxANbax4C/A4uM\nMe8DBphmrf2yySoWEde8++67DBkyhC+++IInn3ySMWPGVN+h6DvImAObX4efnAV/+19Y9aq6OEaA\nQFbLjPazfR9wmWsViUhQLF26lLFjx5KYmEhWVhbJyTV6UX2yHebfD1/nw9DfwJUjoUWME+Raux72\ndIWqSDPj9Xq56667GDlyJD179iQ7O7t6sJeXw4sZ8MBtgIFps2HQaCfYJWKocZhIM1JYWMiYMWN4\n8cUXfTf++uoAzH8AdnwAvfpD6o0Qn1D7jTRiD3sKd5FmYseOHQwZMoTc3FzmzJnD5MmTnQuTKr21\nHhb/j3MCddxU+NkloStWGk3hLtIMrF69mlGDBhEDvPrqq/Sv2vTr+2JY8ihseAV+dAbccDuc2Clk\ntYo7FO4iUcxay0MPPcQdd9xBjzZtyOzRg6Sqwf5ZLqTf76xhv2oUXD0GWioWooH+L4pEqVqNv/bv\nJ6HyxtVeL6xeBssXwbGJcNv9cPrZIa1X3KVwF4lCtRp/eTyYyrXpXbrAz3tCfCmc3w+uuwmO8XGb\nPIloWgopEmU2bNhAcnIyubm5ZGZmMt3jwVTt5Lh3LyxdCZ3Ohol/UrBHKY3cRcJRPTs5Vkrfv58p\nO3bgad2atWeeSffZs313ciwrhwf+CetrtYBqOC2PDCsKd5EoUOr1csvOnczdt4+BiYks6d6dxNhY\nZ6M6OTZLCneRcFSPUXB+fj4jRoxg3b59TJ06lVmzZhETE+OsV38tE97eAoeKa79QnRyjmsJdJILl\n5OQwdOjQ2o2/Cg/C47Ph/bfgmsvh+VVQXCXg1ckx6umEqkiEWrp0KX379qWsrIysrKwjwb4tG/46\nCT7MgWsnwxPLnEZflW0GPB5IS1MnxyinkbtIhPF6vdx9993MmDGDPn36sGzZMjp27Ailh+GFRfDq\nC9A5CW6d5fwN6uTYDCncRSJInY2/9u+BtPtgz04YMBiGj4NWNW6qoVBvVhTuIhHCZ+MvgPUr4ZnH\noFVr+P09cE6vUJcqYUDhLhIBVq9ezahRo4iJiWHNmjWkpKTAoW9h8SOwdQP89Dz47a3Q/vhQlyph\nQuEuEsastcyePZtp06bRo0cPMjMzSUpKgo/ehQUPOqtiRoyHS6+BFlofIUco3EXCVNXGX8OHD2fR\nokUkxMU5J01XPuu05Z3+MHhOC3WpEoYU7iJhqGrjr3vvvZfp06dj8vfDI/fDro+h30AYPQniWoe6\nVAlTCneRMLNhwwaGDRtGUVERmZmZXD1oEGx6DZ6aCzExTrOv5J+HukwJcwp3kTCSnp7OlClT8Hg8\nvP7663T3dHXuabr5dTitB4yfBsd1CHWZEgEU7iJhoLS0lFtuuYW5c+cycOBAlixZQuLX++Fvk+Hr\nfBj6a7jyV9AiJtSlSoRQuIuE2A+Nv9atcxp/zbiXmFVL4cUMOO5EmDYbTu0e6jIlwijcRUKoVuOv\nKy6Df9wJOz6A3gMgdQq0SQh1mRKBFO4iwVZxI46lU6YwduxYEhMTeeONN7iAYqfhl7Uwbir87JLQ\n1ikRTeEuEkwpKXjfeYc/t2/PzJEj6dOnDy88ncFJ/3kBNrwCPzrDOWna4eRQVyoRzm+4G2MWAoOA\nPGttjzr2SQEeAWKBL621F7tZpEi0KCgrY0xREf8uLOSGG25g7tSbaPW/f4P8L2DQtc6flhpzSeMF\n8l20CJgDLPa10RjTHpgHXG6t3W2MOdG98kSiR25uLmlbtjCnrIwVACuWY3a+Bef3gNvug9PPDnWJ\nEkX8NqOw1q4Hvj7KLtcCL1hrd1fsn+dSbSJRY9WqVazr0YMHSkvxAAYweV/Chu3Qvb+CXVznxu9/\nPwFijTFrgbbA/1hrfY7yRSJKxYnPxrDW8tDeveR8+ilP4mM0dbgUJk+BxU82+liAerbLD9wI95bA\n+cAlQBvgTWPMJmttbs0djTETgAkAXbt2deHQIk0kJQVychr1FsXWMr64mIzSUvI4yq/JJSWNPhYA\nPXs2/j0kargR7nuBr6y13wHfGWPWA+cAtcLdWpsGpAEkJydbF44t0nQaEZZ7vv+eaz74gK2lpaT1\nOJ0Ttn1c985xce4Es0btUoUb4Z4JzDHGtARaAb2Ah114X5HQqW9QVk7hrF1LVlYWw4YNo9gY3nv4\n7/T4aDN8vge+Lar9OmNgwQLdrFpc5/eEqjFmCfAmcLoxZq8xZpwxZqIxZiKAtfZDYBXwHrAFmG+t\n3daURYuEq/T0dAYMGEDS8YnsmzaOHtuynLskPfwIxMdX39kYmDhRwS5Nwu/I3Vo7OoB9HgQedKUi\nkQhU6vVy886dzJswgemDB/L3bu1osecTuHYy9L/aCfLW8TBunDPHHhenEbs0KV0tIdJIhY89RlFW\nFv+yltnHxNO6YB+0S4Kp90PnpCM7pqZCerrztebHpYnpposijfDZzJm0nDyZjtbSAmh9qAje/Ah+\n9LPqwV5p7VoFuwSFRu7SvLiwdr3S0vx8em3fTnzNDYdL4Xe/g8cfd+dA+mEgDaBwF6knr7Xc/dln\nPHpgP1/WtVNJSTBLEqlF4S7NSyNHwQUFBYwZM4ZDpYfYMeZSePZ1OFRce0ePRyNuCSnNuYvUJSWl\n2jRObm4u/X7Wm74HcvnPlReQ2PFkzL1/r73EMT4eZswIaqkiNWnkLhKAVatWcee43/BE759wXmIC\n/PxyGDUR4lrDCR2PLHH0eJxg1xJHCTGFu4gvGRmwaRO2pITCxEQ+i4OsKy6gdcIxcP0f4fx+R/bV\nEkcJQwp3kZoyMmDCBCgpwQDtDh5kYkwLvN/F0GL2/8JxHWq/RqEuYUbhLtHBxSWObNpUe7VLuZcW\nK16HvBHuHEM/DKSJ6YSqSA22rmWMWt4oEUQjd4kOLo2El8x5hGu2vEnr4sO1N2p5o0QQjdxFgNLS\nUhZM/A1XbH4RLjgdGxdXfQctb5QIo3CXZu/LvbtZk3oF48oOUBDfjthnX8YsWOB0bgRnxJ6WpuWN\nElE0LSPN2serX6TVEw8zsG0s73vO4qw7Z0HLllreKBFP4S7Nk9fLew/eTfePt5BvDLlDbuCswS6t\nhBEJAwp3aXa8X+Wx689TOPvwt/znUDln3p/GGd1Orb2jRuwSwRTu0qx8l7WG8oUP0bG8jEWxHRi9\n+HHiWrcOdVkirlO4S/NQ8j0H0x+kfc4Gtn5TyIf9hvCbP96OMSbUlYk0CYW7RL89n3LoH3fR/tuv\n+deOLzj77tmkXnJJqKsSaVIKd4kOle0Hqs6TW4tds5zy59Io/K6Y2/YWceeTmXg8nlBUKBJUCneJ\nToUHKZ//ADHbt/Ly53msOP5U/mflIhISEkJdmUhQKNwl+mzLpjz9AUq/LeAPmz6ic+oE0qdP1/y6\nNCsKd4k8NadgKnqvU1ICHY6HMzuTe2IiYzd+zJ8fm8+gQYNCValIyCjcJbJV6b0OwJdfU/rGNyw6\n8SQW/ed1unfvHtr6REJE4S7B5Ubf9ZycI+/lo/d6rNcy85tviJk0qfHHAl3MJBFJ4S6R5cABKCwE\na33fVKNCjHqvSzPnN9yNMQuBQUCetbbHUfa7AHgTGGWtfd69EiWqNGYUXDkFY63zuKQEC/g8Tare\n69LMBTJyXwTMARbXtYMxJga4H3jFnbIkajTx7e8M4KVG7+oWLaBNm4YfWz8UJAr47edurV0PfO1n\nt98Dy4A8N4oS8amOqRYDULnMMS4OfvITOOmkoJUlEo4aPedujOkMXAP0By5odEUSXeo7Ck5JcU6Y\n9ux55LXWwsY1cPkWOFRc6yXG44GkpIYdTyRKuXFC9RFgmrXW6+8iEWPMBGACQNeuXV04tES9okPw\n1L9gyzpKLuuH/b81tK6cc4cjt7/TXZJEqnHjNnvJwDPGmM+A4cA8Y8xQXztaa9OstcnW2uQOHTq4\ncGiJKhkZsGEDFBTAunWQ2B5+dTlkv0HumX04+bUt3JSQwPexsc7+uv2dSJ0aPXK31nar/NoYswj4\nt7V2eWPfV5qRlBRnieOnn0JZ2ZHnDxZgX9rIhmsGc/Ft93DWWWdx5/LltB471tmuKRiROgWyFHIJ\nkAKcYIzZC/wFiAWw1j7WpNVJeGvf3p33OXQIyst9bjLWcsoLmQyPjWXhrl0k9OzZuOMfPNjAIkUi\ni99wt9aODvTNrLVjG1WNND+HD9cZ7JU8wDNt2qjxl0g96ApVabjGjoIrL0ryw3g88NlnjTuWSDOj\ncJfAuXlBEhy1fUA1uiBJpN4U7hI6dQR75UJHExMDp52mC5JEGkDhLoFzaxT8zZew4CHI3gzffV97\nuzGYiy7SqFukEdxY5y4SuHc2wl8nwacfwh//gLd16+rbW7TAnH66gl2kkRTuEhwl38OT/4S5f4Pj\nT4K75/Lsmefx2/Jy9lIxFePxwOLF8OGHIS5WJPJpWkbc46svDMCeTyHtPti/GwYOp3zwGP58z9+Y\nNWsWffv2JdbrxbRqpdG6iIsU7tJ0rIXXMuH5BZBwDPxxJgWdTyV12HBeeuklxo8fz5w5c2jVqlWo\nKxWJOgp3cU/lXZLWrYOuXeCiZGhdDGf3gutvIXd/HoN79WLnzp3MnTuXSZMm6cIkkSaiOXdxx+TJ\n8NFHR+6StGcvPPsinHgm/P6vrMx6kwsvvJCvvvqKNWvWMHnyZAW7SBMytmr71CBKTk622dnZITm2\nVHDroqQDB5xg98HGxfFgp07csWsXZyckkNmjB56aK2QCpTl5EYwxb1trk/3tp2mZ5qry5KcbCgvr\n3GRLSpi2axcjY2NZGBNDQh0/BPyq2jBMRPxSuDdnbgXmunV1btoN3JuUxPSuXRs3DaNRu0i9KNyb\nKzfC8lAhPPFInVeaeoGvb72VPz30UOOPJSL1ohOq0jAf5jhXmr63BW6c4NzurgovcHDUKM5TsIuE\nhEbuUrfKm2FUbe1bVgrLn4TVS+GkzvCHe6Drjyn76fkU/Pa3JJaXk9+iBQmPPspxAbTzFZGmoXCX\nwB34L6TfB5/tgIuugF/9DuJak5+fz/AFC1hfXs7tXbowc9cuYmJiQl2tSLOmcBffMjKcVTDWOj1f\nxo6GLz+CmJYw6S44vx8AOTk5DBkyhLy8PDIyMrj22mtDXLiIgObcxZeMDPjtb49ckLR7N9z7AHxd\nBn999Idgf+655+jTpw9er5esrCwFu0gY0cg9Wrh5l6QNG6CsrPpzXgvPr4b9Iyi3lj/v2sWsPXvo\nc+yxLOvShY633lq/Y2hpo0iTUrhHAzcvSILawV7l+YJ33iG1qIiXysq4ITaWOUDc9u3uHVtEXKFw\njxZuXZDU0sBra31uskCv0lJ2lpcz98c/ZlKnTuoPIxKmFO7RoCFTHJXTOFVfu2Wtc0ONrFgoKa31\nkq+ArxISePXll0lx+2bZIuIqhXtz4yuUvy+Cp+fBxjVwand45BH4w81QeiTgvwdmd+nCW+vXk5SU\nFKxqRaSBFO7N0YEDsGsXlJRA585w/qlw0jEw6Fq4OhViYqBtO7zXXQfWshtY0asXd732GgkJCaGu\nXkQCoKWQzc2BA5Cb6wQ7wL59sHIjnHExDP21E+zAnosu4gJjaAksmTmT37/5poJdJIKon3skCELf\ndeLioHdvALIKChj2wQcUe7083b07g44/vn7H0TJHkSYTaD93vyN3Y8xCY0yeMWZbHdtTjTHvGWPe\nN8ZsNMac05CCpYlVjtjrUjGST9+/nwHvvku7li3ZfO659Q92EQkLgcy5LwLmAIvr2L4LuNha+40x\n5gogDejlTnnNmNurUXbtAq+3zs02Lo4bd+xg3r59DExMZEn37iTGxtbvGBqxi4QNv+FurV1vjEk6\nyvaNVR5uAk5pfFnSIFlZzt/9+lV/vpU5Msfug23RgntbtGDevn1MPeUUZv3oR8Ro/bpIRHN7tcw4\nYKXL79k8NWQUXNmit/K11sJrmfD8AmgbD98W1XqJbdGCm9q3J72oiKeeeorU1NQGlywi4cO1cDfG\n9McJ935H2WcCMAGga9eubh1aoHoXx6QkuGs6HNoD296Cc3rBw4OdtetFRwK+zBhuiInhtfh43li9\nmuRkv+doRCRCuBLuxpizgfnAFdbar+raz1qbhjMnT3JycmiW6USTynn58eOrd3H8/HP43STofw78\n6R5IGQTGQOt4GDcOW1LCwWOPZUphIZ9fcAFvLVtGx44dQ/afISLua3S4G2O6Ai8A11lrj7IcQ6px\n44RpZbOwsWN9dHH0wvr3oWw23DP7h6cL4uJILS/npcJCxp98Mv+JiSFu1KjAjqcTpiIRw2+4G2OW\nACnACcaYvcBfgFgAa+1jwN3A8cC8iiZSZYGswWz2GtvF8fBhKC4++j6lZdWOk1tezuCiInZ6vcxt\n3ZpJ332Hee+9wI7nVmMyEQmKQFbLjPaz/QbgBtcqai4aE5b+1qz7OM6qr79m1PbtxMbEsOass7i4\n8uRroDRqF4ko6i0TLG6uW/ezZv0HMTFYa3lo717u+PRTzkpIYHmPHiS1bh3YcRToIhFL4R6JjrJm\nvarDp57K9R99xNN5eYzo0IHHTz+dBN24WqRZULgHS31Hwb76rQOUlcLJJ8OXdS5KAo+Hr/74RwYu\nXszWHTuYMWMGd955p26sIdKMKNzDla8Trgf+C+n3wZmd4M1COFzlhhrx8c5SyFatyHrqKYYNG0Zx\ncTErVqxg0KBBwatbRMKCWv5GAmthwyvwtymQtx9mz4OFjzudHAE8HkhLg6Ii0h98kAEDBtCuXTs2\nb96sYBdppjRyDyeVK1iuvRYKCpyvW7aEvufDaYlw+tkwbioc1wHO7wfp6c4+a9dSWlrKzVOmMG/e\nPAYOHMiSJUtITEwMzX+HiIScRu7hpqgIHn30yOPycli/BfZ74dZZTrDXkJ+fzy9+8QvmzZvH7bff\nzksvvaRgF2nmdLMON9R3zbgvRUXV7lnqU7t2tZ7KKS9nyHffkWct89u0IbVVK//H6tlTyxxFIpRr\nN+uQIAgk2H147vBh+hw6hBd4IyEhsGAXkWZBc+4N5cZovVI9g73cWv5cUsKskhL6xsTwfHw8HVv4\n+Tl98GAjChSRSKNwjxQVd0UqsJbUoiJeKitjfGwsc9q0oZXWr4tIDQr3hnJrJPzpx/Dj7kfa9dYU\nEwMTJsC8eeTm5jJ48GB27tzJvHnzmDhxoi5MEhGfFO6h4i2HVc9D5mI458eQs6P2PrGxTvdHYOXK\nlYwePZrY2FjWrFnDxRdfHOSCRSSS6IRqKHzzJfxjOrzwOJzbB7LehkmTjmyPiXEeHz6MtZYHHniA\nq666im7dupGdna1gFxG/NHIPtnc2wqKHofQwjL0F+l7m3CVp3jx4+mlnn4opn+LiYm644Qaefvpp\nRo4cycKFC0lISAhh8SISKRTuwVLyPTyXButeBs9pMH4adDyl+j5Verzv2bOHa665hq1btzJz5kzu\nuOMOza+LSMAU7sGw51NIuw/274aBw+Ga30DL2Nr7VVxYlJWV9UPjr8zMTK6++urg1isiEU/h3pS8\nXngtE5YthIS28MeZ8NPzjvqStLQ0brzxRpKSkli7di3du3cPUrEiEk0U7m7w1Xu94Bt4fDZsy4Zz\nejnz623rvvDp8OHD3HzzzTz66KNcfvnlLFmyhPZuXiglIs2Kwr2xUlKc3utV74n6/luwcDZ8XwSp\nN0LKVc4mNbVjAAAIG0lEQVRJ0zrk5eUxYsQI1q9fz+23387MmTOJ0R2TRKQRFO6NdeAAFBbCunWQ\nlARjr4W9W6FzEtx2n/P3UbzzzjsMHTqUvLw8nnrqKVJTU4NRtYhEOa1zb4zJk+Gjj45cXfr55/DA\nI3BsN7jrn36D/dlnn6Vv3754vV6ysrIU7CLimuY1cq+cG3fDgQNOsNdUXAzzHoetH9b5Uq+1/Pmz\nz5i5ezd9jz2WZV26cNKtt/reWa15RaQBmle4+7ovaUMVFta9raSkzmMVWMuYoiL+Xdn4C2i1fbvv\n96k6jy8iUg/NK9zdDMt16+reFhfn81i5RUUM3raNneXlzP3xj5nUqdPRL0zSqF1EGqh5hbubYZmU\n5Myx12QMLFgANebPf2j8lZDAmpUr1R9GRJqUTqg21IwZEB9f/TljYOLEasFetfFXUlKSGn+JSFD4\nDXdjzEJjTJ4xZlsd240x5p/GmE+MMe8ZY45+CWa0SE2FtDRnCgacv5980mkAVqG4uJgxY8Ywbdo0\nRowYwYYNG/B4PCEqWESak0BG7ouAy4+y/QrgtIo/E4BHG19WiKSk1G9FTWoq9O7t3Li6d+9qI/Y9\ne/bQr18/lixZwsyZM3nmmWfU0VFEgsbvnLu1dr0xJukouwwBFltrLbDJGNPeGHOytXa/SzWGNx/z\n+FUbf61YsYJBgwYFvy4RadbcmHPvDOyp8nhvxXPN0vLlyxkwYADt2rVj8+bNCnYRCYmgnlA1xkww\nxmQbY7Lz8/ODeWj/MjJg06YjbQQyMhr0NhdeeCGjR49my5Yt6ugoIiHjRrj/F+hS5fEpFc/VYq1N\ns9YmW2uTO3To4MKhXZKR4dyEuqTEefz5587jBgR8p06deOKJJ9TRUURCyo117iuAG40xzwC9gIIm\nn293s40AOCP2ymCvVFQE48ZBero7x9AFSSISRH7D3RizBEgBTjDG7AX+AsQCWGsfA14GrgQ+AYqA\n65uq2CZTM9j9PS8iEuYCWS0z2s92C0xxraJAuD0KrutqU49HI24RiUi6QhV8X20aH+88LyISgRTu\nUPtqU4/Heaz+6iISoZpX47CjSU09cvJUUzEiEuEU7lUp1EUkSmhaRkQkCincRUSikMJdRCQKKdxF\nRKKQwl1EJAop3EVEopDCXUQkCincRUSikMJdRCQKGaepYwgObEw+4KMVo08nAF82YTlui7R6IfJq\njrR6IfJqVr1NryE1e6y1fu92FLJwrw9jTLa1NjnUdQQq0uqFyKs50uqFyKtZ9Ta9pqxZ0zIiIlFI\n4S4iEoUiJdzTQl1APUVavRB5NUdavRB5NaveptdkNUfEnLuIiNRPpIzcRUSkHsIy3I0xxxljXjXG\n7Kj4O9HHPl2MMa8bY7YbYz4wxtwUgjovN8Z8bIz5xBhzh4/txhjzz4rt7xljzgt2jTXq8VdvakWd\n7xtjNhpjzglFnTVqOmrNVfa7wBhTZowZHsz6fNTht15jTIoxJqfi+3ZdsGusUYu/74l2xpgXjTHv\nVtR7fSjqrFLPQmNMnjFmWx3bw+ozV1GTv5qb5nNnrQ27P8ADwB0VX98B3O9jn5OB8yq+bgvkAj8N\nYo0xwE7gR0Ar4N2axweuBFYCBugNbA7hv2kg9fYBEiu+viKU9QZac5X9/gO8DAwP53qB9sB2oGvF\n4xPDvN7plZ8/oAPwNdAqhDVfBJwHbKtje9h85upRc5N87sJy5A4MAZ6o+PoJYGjNHay1+621Wyu+\n/hb4EOgctArhQuATa+2n1trDwDM4dVc1BFhsHZuA9saYk4NYY1V+67XWbrTWflPxcBNwSpBrrCmQ\nf2OA3wPLgLxgFudDIPVeC7xgrd0NYK0NZc2B1GuBtsYYAxyDE+5lwS2zSjHWrq+ooS7h9JkD/Nfc\nVJ+7cA33k6y1+yu+/gI46Wg7G2OSgHOBzU1bVjWdgT1VHu+l9g+XQPYJlvrWMg5nBBRKfms2xnQG\nrgEeDWJddQnk3/gnQKIxZq0x5m1jzK+DVl1tgdQ7B+gO7APeB26y1nqDU16DhNNnriFc+9yF7AbZ\nxpg1QEcfm/5U9YG11hpj6lzSY4w5BmfUdrO1ttDdKpsnY0x/nG+yfqGuJQCPANOstV5ncBn2WgLn\nA5cAbYA3jTGbrLW5oS2rTgOBHGAAcCrwqjHmDX3W3Of25y5k4W6t/UVd24wxB4wxJ1tr91f8SuXz\nV1djTCxOsGdYa19oolLr8l+gS5XHp1Q8V999giWgWowxZwPzgSustV8Fqba6BFJzMvBMRbCfAFxp\njCmz1i4PTonVBFLvXuAra+13wHfGmPXAOTjnjIItkHqvB+6zzoTwJ8aYXcAZwJbglFhv4fSZC1hT\nfO7CdVpmBfCbiq9/A2TW3KFiDnAB8KG19h9BrK3SW8BpxphuxphWwCicuqtaAfy64gx+b6CgynRT\nsPmt1xjTFXgBuC5MRpJ+a7bWdrPWJllrk4DngckhCnYI7HsiE+hnjGlpjIkHeuGcLwqFQOrdjfNb\nBsaYk4DTgU+DWmX9hNNnLiBN9rkL9ZnkOs4eHw+8BuwA1gDHVTzfCXi54ut+OCd73sP5tTEHuDLI\ndV6JM+LaCfyp4rmJwMSKrw0wt2L7+0ByiP9d/dU7H/imyr9ndhh8Lxy15hr7LiKEq2UCrReYirNi\nZhvOdGLY1lvxmXul4vt3GzAmxPUuAfYDpTi/BY0L589cgDU3yedOV6iKiEShcJ2WERGRRlC4i4hE\nIYW7iEgUUriLiEQhhbuISBRSuIuIRCGFu4hIFFK4i4hEof8HY5DHC2gmWkAAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x1818a91c50>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "xcorr = (yobs - pfit[1]) / pfit[0] \n", | |
| "plot(xi,polyval(ptrue,xi),label='True Function')\n", | |
| "plot(xi,polyval(pfit,xi),label='Fitted Function')\n", | |
| "errorbar(xcorr,yobs,yerr=yerr,xerr=xerr,fmt='or')\n", | |
| "print \"RMS(xcorr - xobs) = {}\".format(std(xcorr - xtrue))" | |
| ] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 2", | |
| "language": "python", | |
| "name": "python2" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 2 | |
| }, | |
| "file_extension": ".py", | |
| "mimetype": "text/x-python", | |
| "name": "python", | |
| "nbconvert_exporter": "python", | |
| "pygments_lexer": "ipython2", | |
| "version": "2.7.15" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 2 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment