{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Error propagation \n", "\n", "## Motivation\n", "\n", "Say we have some function $f(x)$. We have measurements of the variable $x$ with some error $\\sigma_x$. We want to know the error in $f(x)$.\n", "\n", "The uncertainty in $f(x)$ depends on $\\sigma_x$ and the derivative $df/dx$." ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "tags": [ "remove-input" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAggAAAGKCAYAAABpbLktAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCK0lEQVR4nO3dd3hUddrG8e/MpBdCCb0jJaF3xAKiqCsggqKClNCs7Iq4qGtZYVXAxVUXF10VIVRRXBSV8rorggWFQKihSQ0tEGp6m5nz/hEdiCkmYZIzydyf6+JyOHlm5uEImTunPD+LYRgGIiIiIlewmt2AiIiIeB4FBBEREclHAUFERETyUUAQERGRfBQQREREJB8FBBEREclHAUFERETyUUAQERGRfHxK+0Sn08mpU6cIDQ3FYrG4sycREREpI4ZhkJKSQr169bBaCz9OUOqAcOrUKRo2bFjap4uIiIiJjh8/ToMGDQr9eqkDQmhoqOsNqlSpUtqXERERkXKUnJxMw4YNXZ/jhSl1QPj1tEKVKlUUEERERCqY37s8QBcpioiISD4KCCIiIpKPAoKIiIjko4AgIiIi+SggiIiISD4KCCIiIpKPAoKIiIjkU+o5CCIiInJ1HE6DtGw7qZl2UrPspPzy39zf5+T5vdVqoVqQHzWC/age7Ef1kNzH1YL9CPX3cfuyBwoIIiIiZSQpPYe9p5PZm5DMvoQUDiSmcCkjxxUI0rMdbnkfP5uVasG+VA/2d4WGX4NEeIg/betXIbJuFXxtxT9xoIAgIiJylRxOg6Pn09ibcDkM7E1I5lRSZrGe72uzEBrgS4i/D8H+PoT6+xAS4EPIL/8N9ffB4TS4kJ7NhbTcX+dTs7mYnk16toNsh5MzyVmcSc4q9D0CfK20b1CVNuG+xepJAUFERKQE0rPt7DqR9EsYSGHf6WT2n0khM8dZYH39qoFE1q1CZN1QWtUJJTzEnxB/H0KvCAD+PrZS95OR7cgNDqnZnE/L4mJ6bnj4NUgkJGWy/fglkjJyiDlygY370ov1ugoIIiIiv+PUpQzW7j3D13sT+enwebLt+cNAgK+VVnWq0LpuKBF1cg/pt6oTSlhg8X5iL61APxv1/QKpXzWw0Bqn0+DwuTS2xl/kx73HmVWM17UYhmGUpqHk5GTCwsJISkrSYk0iIlKpOJ0Gu04muULBnoTkPF+vGxZAm3q5ISA3DITSuEYwNqt7LxQsC8X9/NYRBBEREXIP1W84eI61+86wdm8iiSmXz+dbLNC5UTX6Rtamb2QtmtcKcftdA55GAUFERLzWmeRM1u5NZO3eM/xw8BxZV5w6CPaz0atlTW6JrE2fVjWpEeJvYqflTwFBRES8SlqWnc+3n+LjzcfYcSIpz9fqVw3klsha9I2sTY9m1a/q4sGKTgFBRES8woEzKSzeGM+nW0+SkmV3be/YsCp9I2txS2RtIuqEVvpTB8WlgCAiIpVWtt3Jf/ecZtFP8Ww6csG1vUmNIEZc25iBHetRKzTAxA49lwKCiIhUOqcuZbA05hgfbT7O2V8uNrRaoG9kbUb2bMz114RjrQB3HJhJAUFERCoFp9Pgh4PnWLwxnq/3nsH5y038NUP9GdatIUO7N6JeEbMCJC8FBBERqdAupWfzyZYTLNkUz9Hzl6cE9mxWgxHXNua2NrVLtAaB5FJAEBGRCikhKYPZ3xzkP7EnXLcnhvr7cE+XBgzv0YgWtUNN7rBiU0AQEZEK5WxKFu+sP8iSTcdcI49b163CyJ6NuatjPYL89NHmDtqLIiJSIVxKz+a97w4zf8NRMnJyl0nu3rQ6f761Jd2bVtftiW6mgCAiIh4tJTOHuT8cYe73R1zzCzo0rMpTt7Xi+uY1FAzKiAKCiIh4pPRsOwt+jOe97w5xKT0HgMi6VZh8W0tujqilYFDGFBBERMSjZOY4+HDTMd5Zf5BzqdkAXFMzmCdvbcUdbetofkE5UUAQERGPkG138knscf619iCnkzMBaFQ9iCf6tuCujvUrxFLKlYkCgoiImMrpNPhs20n+ufZnjl/IAKBuWACP39KCIV0aaIaBSRQQRETENLtPJfHCiji2HbsEQHiIPxP6XMOw7o0I8PXelRQ9gQKCiIiUu5TMHN74388s+PEoTgOC/WxMuLk5o69rojkGHkL/F0REpNwYhsHKnQm8vHIPib8sotS/fV3+2r81dcK0qqInUUAQEZFycfhsKlO+2M33B84BuUsuv3RXW3q1rGlyZ1IQBQQRESlTmTkO3ll3kHe/PUy2w4mfj5XHbrqGR3pfo+sMPJgCgoiIlJl1+xOZ8vlujl3IXWWxV8uavDSwDU3Cg03uTH6PAoKIiLjdqUsZvLxyD2viTgNQp0oAL97Zmjva1tEExApCAUFERNwmx+EkesMR/vn1AdKzHdisFsZc14Qnbm1JiL8+cioS/d8SERG32HrsIs8u38X+MykAdGlcjVcGtSWybhWTO5PSUEAQEZGrkuNw8tbaA7y97iBOA6oF+fLsHZEM6dJA6yZUYAoIIiJSagcTU5n08XZ2nUwCYFDHeky5sw3Vgv1M7kyulgKCiIiUmGEYLNoYz/TVe8nMcRIW6Mu0wW0Z0L6e2a2JmyggiIhIiZxJzuSp/+zku5/PAnBji3BeG9JBkxArGQUEEREptjW7Enj2s11cSs/B38fKX+6IIKpnE11rUAkpIIiIyO9Kzsxh6he7+XTrSQDa1KvCP+/vSIvaoSZ3JmVFAUFERIoUc+QCkz7ezslLGVgt8Ejva3iib0v8fKxmtyZlSAFBREQKlGV38Ob/DvDed4cwDGhQLZA37+9ItybVzW5NyoECgoiI5PPzmRQmfrSdvQnJANzbpQEv3tma0ABfkzuT8qKAICIiLoZhEL3hKK/+3z6y7U6qBfky4+52/KFtXbNbk3KmgCAiIkDuhYiTl+3gv3vOAHBTq5rMHNKeWqG6fdEbKSCIiAj7Tifz6OKtHDmXhp/NygsDIhl5bWOtvOjFFBBERLzcim0nefbTXWTkOKgXFsA7I7rQsWFVs9sSkykgiIh4qWy7k1dW7WHhT/FA7kTEWUM7UV3rKAgKCCIiXikhKYPHlmxl27FLAPzp5uY80bclNk1ElF8oIIiIeJkNB8/xp6XbuJCWTZUAH968vyO3RNY2uy3xMAoIIiJewuk0+Pe3h3j9v/txGtC6bhXeHdGFRjWCzG5NPJACgoiIF0jKyOHPy3bw9d7cWxjv7dKAlwe1JcDXZnJn4qkUEEREKrm9Cck8sjiW+PPp+Nms/O2uNgzt1lC3MEqRFBBERCqx5bEneH7FLjJznNSvGsi/R3SmfYOqZrclFYACgohIJZRld/DSl3tYsukYAL1a1mTW/R2pplsYpZgUEEREKpnzqVk8vCiWLfEXsVjg8Ztb8PgtLXQLo5SIAoKISCVy4EwKYxds5viFDEIDfHhraCf6RNQyuy2pgBQQREQqie8PnOWxxVtJybLTqHoQ80Z3pXmtULPbkgpKAUFEpBJYvDGeKV/sxuE06NakGu+N7KqRyXJVFBBERCowh9Ng2qq9zNtwBIDBnerz6j3t8PfRfAO5OgoIIiIVVGqWnYlLt7F2XyIAk29ryYQ+zTXfQNxCAUFEpAI6dSmDsfM3s+90Cv4+Vl6/rwMD2tczuy2pRBQQREQqmB3HLzF+4RbOpmQRHuLPnFFd6NSomtltSSWjgCAiUoGs3pXAk8u2k5njJKJOKB9EdaVBNS22JO6ngCAiUgEYhsE76w/x2lf7AejTqiZvDetEaICvyZ1JZaWAICLi4bLsDp77NI7lW08AMPq6JrzQPxIfm9XkzqQyU0AQEfFgF9KyeWRRLDFHL2CzWph6Z2tG9mxidlviBRQQREQ81LHz6Yyat4mj59MJ9fdh9vDO9G5Z0+y2xEsoIIiIeKDdp5KImreZc6lZ1K8aSPSYbrSsrbHJUn4UEEREPMxPh87z4MItpGbZiaxbhQVjulGrSoDZbYmXUUAQEfEga3YlMPGj7WQ7nPRoWp05UV2pojsVxAQKCCIiHmLJpnheWBGHYcDtbWoza2gnAny1poKYQwFBRMRkhmHw1tqDvPn1zwAM696IVwa1xWbVmgpiHgUEERETOZwGf/tyNwt/igfg8ZubM+nWllpwSUyngCAiYpIsu4MnP97Bql0JWCww9c42RF3XxOy2RAAFBBERU6Rk5vDwolh+PHQeX5uFN+7ryJ0dtBqjeA4FBBGRcnY2JYsx82OIO5lMsJ+N90Z25YYW4Wa3JZKHAoKISDk6fiGdkXNzpyPWCPYjekw32jeoanZbIvkoIIiIlJM9p5KJio7hbEoWDaoFsmhcD5qGB5vdlkiBFBBERMrBxsPneXDBFlKy7ETUCWXB2O7U1nRE8WAKCCIiZex/e84w4cOtZNuddG9anTmjuhIWqOmI4tkUEEREytCXO04x6ePt2J0Gt7WuzVvDNB1RKgYFBBGRMvKf2BM8/Z8dOA0Y3Kk+rw1pj4/NanZbIsWigCAiUgaWbIrn+c/iABjWvSHTBrXDqtHJUoEoIIiIuNncH47w8so9AIy+rglT7myt0clS4SggiIi40dvrDvLaV/sBeLh3M/7yhwiFA6mQFBBERNzAMAze/N/PvPXNQQCe6NuCibe0UDiQCksBQUTkKhmGwYw1+3j/u8MAPPOHCB696RqTuxK5OgoIIiJXwek0mHrFcs1T7mzNmOubmtyVyNVTQBARKSWH0+C5T3fx8ZbjWCwwbVA7HujRyOy2RNxCAUFEpBTsDieTP9nBiu2nsFrgH/d24O7ODcxuS8RtFBBEREoo2+5k4kfbWBN3Gh+rhVlDO9G/fV2z2xJxKwUEEZESyMxxMGHJVtbuS8TPZuXt4Z25tXVts9sScTsFBBGRYsrIdvDQoi18f+Ac/j5W3h/Vld4ta5rdlkiZUEAQESmGtCw7Y+dvZtORCwT52Zgb1Y2e19Qwuy2RMqOAICLyO9KzL4eDUH8f5o/tRpfG1c1uS6RMKSCIiBQhI9vBuPlbXOFg4bjudGpUzey2RMqc1h0VESlEZo6D8Qs389Ph84T4+7BA4UC8iAKCiEgBMnMcPLhwCxsOnifYz8aCsd3orHAgXkQBQUTkN34NB98fOEeQn435Y7vrmgPxOgoIIiJXyLI7eGRxLN8fOEegr43o0d3o1kThQLyPAoKIyC+y7A4eXbyV9fvP5oaDMd3o0Uy3Mop3UkAQESF3fPKEJVv5Zl8iAb5W5o7uyrUKB+LFFBBExOtl251M+HArX+9NxN/Hytyoblx3TbjZbYmYSgFBRLxajsPJn5Zu5X97zuDnY+WDqK5c31zhQEQBQUS8Vo7DyeNLt/HV7txwMGdUV25sobUVREABQUS8lN3h5ImPtrMm7jR+NivvjeyihZdErqCAICJex+5w8sTH21m1KwFfm4V3R3amT6taZrcl4lEUEETEq9gdTp5ctoOVO3PDwb+Hd+HmiNpmtyXicRQQRMRrOJ0GzyzfxRc7TuFjtfD2A53p21rhQKQgCggi4hUMw+DFL+JYvvUENquF2Q904rY2dcxuS8RjuT0gnD59mqioKOrWrYvNZsNisRAbG0uHDh3o16+fu9/Oox08eBAfHx/eeecds1sR8WqGYTBjzT4WbzyGxQJv3NeBP7Sta3ZbIh7Nx50vZhgGAwYMYPv27dx///00b94ci8XCli1b2LlzJ3PmzHHn23m85s2bM3z4cKZOncqIESOoUqWK2S2JeKW31h7k/e8OAzBjcDvu6ljf5I5EPJ9bA8L69euJjY1l4sSJ/POf/wTA4XDQrFkzevfuTffu3d35dhXCU089xcKFC3nrrbd44YUXzG5HxOvM+e4wb379MwAvDmjN0O6NTO5IpGJw6ymGb775BoDBgwe7tq1evZpjx44xcuRId75VoW666SaaNGlSLu9VHG3btqVDhw7MmTMHp9NpdjsiXmXxxnimrd4LwFO3t2LsDU1N7kik4nBLQIiOjsZisfDKK68AuR/SFouF6tWrM3/+fCwWC/fcc0+Bz+3ZsycWi4XNmzfn2X7x4kXatGlDQEAA3377rTvaLJazZ8/y1FNPERkZSVBQEBaLJd+v22+/vUSved9993Hs2DHWrl1bRl2LyG99uvUEf/08DoDHbrqGCX2am9yRSMXiloDQrFkzpkyZQnBwMLVq1WLKlClMmTKFGTNmsH79eiIiIqhatWqBz3311VcBePHFF13bMjMzGThwIPv27WPJkiX07t3bHW3+rqNHj9KpUyf+8Y9/ULNmTSZOnMiYMWMICAgAoEqVKjRu3JgbbrihRK/bs2dP4PIRFhEpW2t2JTD5kx0YBoy+rglP3d7K7JZEKh6jlJKSkgzASEpKcv3eYrEYgwYNctXs3r3bAIzhw4cX+Vp33HGHARgbNmwwHA6HcffddxuA8c4775S4r969exuNGzcu8fMcDofRvXt3AzDeeuutPF9btmyZARh9+vQp8esahmEkJycbgNGrV69SPV9Eiu+bvWeM5s+tMho/s9J46pPthsPhNLslEY/y28/vwrjtGoTt27djGAadOnVybTtx4gQAtWsXPYhkxowZWCwWXnzxRR5//HE+/fRTXnzxRR599FF3tfe7vvjiC2JiYrjnnnv405/+lOdrAwcOJCgoiJiYmFK9dmhoKAEBAa79ISJl48dD53hkcSw5DoM7O9Rjxt3tsVotZrclUiG57S6Gbdu2AdCxY0fXtvPnzwNQrVq1Ip/boUMHHnjgAZYsWcLatWt56KGH+Nvf/va772mxFP4Pv6CvHTlypNALGD/88EMAJk2alO9r/v7++Pv7k5GR8bs9FaZ69eqcO3eu1M8XkaLFxl9k/IItZNmd9I2szRv3dcCmcCBSam4LCNu3bwfyBoTAwECAYn2whofnrr8eFhbGv/71r2K955QpU/Jtmz9/PpcuXeKJJ57I97XCroMA+OGHH6hSpYrreoErpaWlkZycTERERLH6KkhGRgZBQUGlfr6IFC7uZBKjo2NIz3ZwY4twZj/QCV+bBsWKXA23HkGoXr06jRpdvse4Zs3cpVMvXLhQ5HNnzZrFrFmzqF27NmfOnGHx4sWMHTv2d99z6tSp+batX7+eo0ePFvi1wmRlZZGQkECLFi2wWvN/U1mzZg0Oh4Obb745z/Z58+YxadIkjh49SrVq1bDb7fTr1w+bzcaXX36Jj0/u7nU6nSQlJdGmTZti9yQixXPgTAqj5sWQkmmnW5NqvDeyCwG+NrPbEqnw3BKxs7Oz2bNnT56jBwBt2rTBarVy4MCBQp/70UcfMWnSJPr27cvWrVsJDQ1l6tSpZGZmuqO1EklKSsIwjDzbnE4nb7zxBhaLhXHjxuX5WlRUFLVr12bWrFkAPProo5w9e5Zly5a5wgHAgQMHcDqdtGvXruz/ECJeJP58GsM/2MSFtGzaNwhj3uhuBPm5df6biNdyS0CIi4sjJycnzwWKkHtIv3379mzZsiXfBy/A119/TVRUFB07duTTTz+lXr16TJw4kePHj/P222+7o7Vi8ff3JyIigsTERL7++mvXdsMweO655/jpp58YP348HTp0yPM8m83G1KlTmTVrFs899xxfffUVq1atIjQ0NE/dpk2bAMrtdk0Rb3DqUgYPzNlEYkoWEXVCWTi2O6EBvma3JVJpuCUgFHSB4q8GDRpEUlJSvkFIW7du5e6776ZBgwasWbPG9aE6efJkqlWrxowZM0hOTnZHe8Xy6xyGQYMG8dBDD/HMM8/QpUsX/v73vzNw4EBmz55d4POGDh1KeHg4s2fPZvXq1dSrVy9fzf/+9z9sNhsDBgwo0z+DiLc4l5rFiA82cfJSBs3Cg1k0rgdVg/zMbkukUinzgDB+/HhsNhuLFy92bTt06BD9+vUjICCAr776Ks9tkGFhYUyePJnz58/z2muvuaO9Yhk2bBhLly4lIiKCxYsX8+9//5ugoCDmz5/PihUr8PMr+JvPihUrOH78OE6nk1q1auX7enp6OitWrODOO+8sMDyISMkkZ+YQNS+Gw+fSqF81kMXje1Az1N/stkQqHYtR0LH/YkhOTiYsLIykpKTfXaXwgQce4L///S/x8fEEBweXqlFPtHHjRm699VaWLFnCyy+/zA033MCbb76Zp2bevHmMGzeOb7/9ll69epnUqUjlkJnjYNTcGGKOXiA8xI9PHrmOpuGV53uKSHko7ud3udwHNG3aNFJTU8v1uoKydujQIQYOHMiMGTMYOHAgf/3rX3n33Xc5deqUq8ZutzN9+nQGDhyocCBylXIcTiYs2UrM0QuE+vuwYGx3hQORMlQuAaFp06YsWLCg0hw9uHDhAv369WPEiBH88Y9/BHKnLUZGRroWrILcSZIjRozgjTfeMKtVkUrB6TR46pMdrN2XSICvlXljutGmXpjZbYlUald9iuHUqVMFHqKw2WyuRY4gd9hQYaxWq2uoUklr09PTC7xDAnKnKV45nKgktRkZGUUuz3xl2ClJbWZmJg6Hwy21v642CbmzHOx2u1tqAwMDXfMgsrOzycnJcUttQEAANputxLU5OTlkZ2cXWuvv7++6rbQktXa7naysrEJr/fz88PX1LXGtw+Eo8jZdX19f1zUtJal1Op1FDh0rSa2Pjw/+/rnn7Q3DID093S21Jfl3X9xawzD4xzfxzP/xKD5WC3NGdaVPRP7rfUSkeIp9icDVLvZQ2K9+/frlqQ8KCiq0tnfv3nlqw8PDC63t2rVrntrGjRsXWtu6des8ta1bty609rcLPHXt2rXQ2vDw8Dy1vXv3LrQ2KCgoT22/fv2K3G9XGjJkSJG1qamprtqoqKgiaxMTE121jz32WJG1R44ccdVOnjy5yNq4uDhX7ZQpU4qsjYmJcdXOnDmzyNp169a5amfPnl1k7cqVK1210dHRRdYuW7bMVfvrIlyF/YqOjnbVrly5ssja2bNnu2rXrVtXZO3MmTNdtTExMUXWTpkyxVUbFxdXZO3kyZNdtUeOHCmy9rHHHnPVJiYmFlkbFRXlqk1NTS2ydsiQIXn+DhdVW5LvEY2fWWk0+ctKY8W2E4aIXJ1yX6xJRKQs/W1gG+7qWN/sNkS8hk4xFFKrUww6xaBTDCWvdecphs+3n+Qvy3cBMLl/ex6/pUWhryUixVfcUwzlcpujiEhJfL3nDA8vjsWek02DQ1/Qs1l1pk+fXug8EhEpPo+6zVFEpLg2Hj7PYx9uxeE0uKtdHX78LJrXX3+9yCNOIuJ+Cggi4jHiTiYxfsEWsu1O+kbW5uVBWgFVxCwKCCLiEQ6dTWXUvBhSs+xc26w6sx/ohI9N36JEzKJ/fSJiulOXMhj5y7LN7eqHMWdUVwJ8bWa3JeLVFBBExFTnU7MYOXcTp5IyaVYzmPljumnZZhEPoIAgIqZJycxhdPRmDp1No15YAIvH9aBGiFZmFPEECggiYoosu4OHF8Wy62QSNYL9WDS+B/WqBv7+E0WkXPiY3YCIeB+H0+DJj3fw46HzBPvZmD+mO9fUDMlXFxgYSFxcnOuxiJQfBQQRKVeGYTD1i92s2pWAn83K+6O60q5BWIG1VquVNm10q6OIGXSKQUTK1VtrD7JoYzwWC7xxfweubx5udksiUgAdQRCRcrN4Yzxvfv0zkLv40oD29Yqsz87OZvr06QA899xzGrUsUo60FoOIlIs1uxJ47MOtGAY8fnNznryt1e8+Jy0tjZCQ3GsTUlNT8yxmJiKlo7UYRMRj/HToPBM/2o5hwLDujZh0a0uzWxKR36GAICJlavepJB5auIVsh5Pb29TmlUFtXUuPi4jnUkAQkTITfz6NqHmbScmy06NpdWYN7YTNqnAgUhEoIIhImTibksWoeTGcS80ism4V5kRpfQWRikQBQUTcLneEcgzx59NpWD2QBWO6UUXrK4hUKAoIIuJWv45Q3n0qmRrBfiwc24NaVQLMbktESkhzEETEbRxOg0kfb88zQrlpeOlvTQwICCAmJsb1WETKjwKCiLiFYRhM+SKO1btO/+4I5eKy2Wx069bNTR2KSEnoFIOIuMWstQdYvPGYRiiLVBI6giAiV23xxnj++fUBoHgjlIsrOzubWbNmATBx4kSNWhYpRxq1LCJX5f/iEnh0SclGKBeXRi2LuJ9GLYtImYs5coHHXSOUG2qEskglooAgIqXy85kUxi/YTLbdSd/I2rx8l0Yoi1QmCggiUmKnLmUQNS+G5Ew7XRpX41/DOuFj07cTkcpE/6JFpESS0nOImhdDQlImzWuFMDeqK4F+GqEsUtkoIIhIsWXmOBi/cDMHElOpXcWfBWO7UzVIdxaIVEYKCCJSLA6nweNLt7H56EVCA3xYMLY79asGmt2WiJQRzUEQkd9lGAZ//TyO/+45g5+PlTmjuhJRp+xvbw4ICGDdunWuxyJSfhQQROR3/eubg3y4KXdK4qz7O3Jtsxrl8r42m42bbrqpXN5LRPLSKQYRKdJHMcd4438/AzD1zjbc0a6uyR2JSHnQEQQRKdTXe87w3Ge7AJjQ5xqirmtSru+fk5PD+++/D8BDDz2Er69vub6/iDfTqGURKVBs/EWGf7CRzBwnQ7o04LUh7ct9EJJGLYu4n0Yti0ipHUxMZdyCzWTmOOnTqiYz7m6nKYkiXkYBQUTyOJOcSdS8GC6l59ChYVXeHt4ZX01JFPE6+lcvIi5JGblTEk9eyqBZeDDRo7sR5KdLlUS8kQKCiACQZXfw8KIt7DudQs3Q3CmJ1YM1JVHEWykgiAhOp8GTy3aw8fAFQvx9mD+mGw2rB5ndloiYSAFBxMsZhsHLq/awamcCvjYL743sQpt6YWa3JSIm08lFES/3/neHid5wFIDX7+vI9c3DzW3oCv7+/qxcudL1WETKjwKCiBf7bNsJZqzZB8AL/SMZ2KGeyR3l5ePjQ//+/c1uQ8Qr6RSDiJf6/sBZnvpkJwDjbmjK+BubmdyRiHgSHUEQ8UJxJ5N4ZFEsdqfBnR3q8Xy/SLNbKlBOTg5LliwBYPjw4Rq1LFKONGpZxMscv5DO4Hd+5FxqFj2b1WD+2G74+9jMbqtAGrUs4n4atSwi+VxIy2bUvBjOpWYRUSeU90Z18dhwICLmUkAQ8RLp2XbGzt/MkXNp1K8ayIKx3akSoEP2IlIwBQQRL2B3OPnTh9vYfvwSYYG+LBjbjdpVAsxuS0Q8mAKCSCVnGAYvrIhj7b5E/H2szBvdlea1Qs1uS0Q8nAKCSCX3z68P8NHm41gt8K9hnejSuLrZLYlIBaCAIFKJfbjpGLPWHgDg5UFtua1NHZM7EpGKQnMQRCqp/+05wwsrdgHw+M3NGd6jsckdlZy/vz/Lli1zPRaR8qOAIFIJxcZf5E9Lt+I04L6uDZh0a0uzWyoVHx8f7r33XrPbEPFKOsUgUskcOpvK+AWbycxx0qdVTaYNbofFYjG7LRGpYHQEQaQSSUzOZNTcGC6m59ChYVXeHt4ZX1vF/TnAbrfz2WefATB48GB8fPQtS6S86F+bSCWRnJlDVPRmTl7KoEmNIOZFdSXIr2L/E8/KyuK+++4DckctKyCIlJ+K+6OFiLhk2R08siiWvQnJhIf4s3BsD2qE6KI+ESk9BQSRCs7pNJj8yU5+PHSeYD8b88d0o1GNILPbEpEKTgFBpIKbtnovX+44hY/Vwrsju9C2fpjZLYlIJaCAIFKBzfnuMHN/OALAP+7twI0taprckYhUFm4NCKmpqTz44IPUr18fHx8fmjVrxksvvUTr1q1xOp0lfr25c+dSv3590tLS3NmmSKWwYttJpq3eC8Bz/SIY1Km+yR2JSGXi1oDw5JNPsnz5cl5//XW+//57Fi9ezMyZM3nppZewWkv+VlFRUQQHBzNz5kx3tilS4X1/4CxP/WcHAONuaMqDNzYzuSMRqWzcFhCys7NZunQp48ePZ+jQofTs2ZPPP/+cqlWrcvfdd5fqNX18fHj44YeZNWsW6enp7mpVpEKLO5nEI4tiyXEY3NmhHs/3i6y0g5D8/PyIjo4mOjoaPz8/s9sR8SpuCQhjxozB39+f1NRUXnvtNSwWC507d2bu3Lk88MAD+Y4eJCQkEBISwtChQ/NsX7lyJb6+vjz//POubcOHDyc5OZmPPvrIHa2KVGjHzqczOnozadkOejarwT/ubY/VWjnDAYCvry+jR49m9OjR+Pr6mt2OiFdxS0B45plnePbZZwH44osv+Omnn5g1axbnz5+nT58++err1q3L008/zbJly4iNjQVg/fr13HvvvTz66KNMmzbNVVunTh0iIiJYtWqVO1oVqbDOp2YRFR3DudQsIutW4b1RXfD3sZndlohUUm4JCBEREaSmplKtWjXuvPNOrr32Wn766ScAOnfuXOBzJk+eTN26dXnmmWfYvHkzAwcOZNiwYcyaNStfbefOndmwYYM7WhWpkNKz7Yydv5kj59KoXzWQBWO6USWg8v9EbbfbWbVqFatWrcJut5vdjohXcdvc0tjYWLp06eL6/alTp7BYLISHhxdYHxQUxCuvvMLYsWPp06cP/fv3Z86cOQWeS61VqxaJiYnY7XaNWhWvk+NwMmHJVnacSKJakC8Lx3WnVpUAs9sqF1lZWQwYMADQqGWR8uaWIwgOh4Pt27fnCQgZGRn4+vpisxV+CLRly9wlaC0WC/Pnzy+0NiAgAMMwyMzMdEe7IhWGYRg8++ku1u0/S4Cvlbmju3FNzRCz2xIRL+CWgLB3717S09PzBITw8HCys7MLnWGwfft2BgwYwPXXX09qairz5s0r9PUvXLiAv78/ISH6xije5fX//sx/Yk9gs1p4+4HOdG5UzeyWRMRLuCUgbNmyBSBPQIiIiADg0KFD+er379/P7bffTs+ePVm3bh133XUXU6dOJSkpqcDXP3z4MK1bt3ZHqyIVxqKfjjJ73UEApg1qyy2RtU3uSES8iVsCQmxsLFWrVqVZs8vDWm666SYANm7cmKf26NGj9O3bl1atWrF8+XJ8fX159dVXuXjxItOnT8/32k6nk5iYmALvhhCprP4vLoEXv9gNwKS+LRnavZHJHYmIt3FbQPjt3QoNGzbkxhtv5PPPP3dtS0hIoG/fvtSqVYuVK1cSGBgI5B5tGDt2LLNmzeLo0aN5Xmf9+vUkJSUxfPhwd7Qq4vF+OnSex5duxzDggR6NePyW5ma3JCJeyGIYhlGaJyYnJxMWFkZSUhJVqlQpsGb58uXcf//9xMfHU79+6ebEjxw5ksOHD+s2R/EKe04lc/97P5GSZef2NrV5Z3gXbJV4ENLvSUtLc117lJqaSnBwsMkdiVR8xfn8Bjfe5liQu+++m27dujFjxgxmz55d4ucfOnSIjz/+mG+++aYMuhPxLMcvpBMVHUNKlp3uTasza2gnrw4HkDtq+dfvHRq1LFK+rjogpKWlFXh7os1mIyAggDlz5vDFF1+QkpJS6IJNVqvVdbrh19eE3IsZX3/9dTp16uTa9tva9PR0CjsIYrFYCAoKKlVtRkZGkStQXvmTTElqMzMzcTgcbqkNCgpyzY3IysoqcpBMSWoDAwNd/6+ys7PJyclxS21AQIDr70pJanNycsjOzi601t/f33V/fElq7XY7WVlZhdb6+fm5xvuWpNbhcBR5S66vr6/rw+7X2vOpWUQt2snZlCwi6oQyZ1RXAnw1JdHX15cJEyaY3YaIdzJKKSkpyQAK/dWvX7889UFBQYXW9u7dO09teHh4obVdu3bNU9u4ceNCa1u3bp2ntnXr1oXWNm7cOE9t165dC60NDw/PU9u7d+9Ca4OCgvLU9uvXr8j9dqUhQ4YUWZuamuqqjYqKKrI2MTHRVfvYY48VWXvkyBFX7eTJk4usjYuLc9VOmTKlyNqYmBhX7cyZM4usXbdunat29uzZRdauXLnSVRsdHV1k7bJly1y1y5YtK7I2OjraVbty5coia2fPnu2qXbduXZG1M2fOdNXGxMRc/jv4zErjuhlrjdNJGYaISFn59fM7KSmpyDq3LvcsIiWTY7985MnnzB6iR3ehtpdMSSwOh8PB+vXrWb9+fZFH00TE/a76IsVTp04VeJHDr6cYflXYwCQo/BRDcWp1ikGnGCrqKQan0+CJj2L59Kf9nJw9AtCFeL+lixRF3K/cLlIMDg4u1j/akvzDLkntlR/q7qy9MoS4s/bK0OTOWn9/f/z9/d1e6+fnV+yLw8qq1tfXt9hL/Zak1sfHp9iz/UtSa7PZivV3eMaavXyx8wx+/sX/+yMiUl50ikHEBO9/d4g53x8B4JXBbU3uRkQkPy2NJlLOPt16gumr9wHwXL8I7upYx+SORETy0xEEkXK0bn8iT/9nJwAP3tiUh3pdY3JHIiIFU0AQKSfbjl3kscVbsTsNBnWsx7N3RJrdkohIoRQQRMrBobOpjJ2/mYwcB71a1mTmkA5YvXxKooh4Nl2DIFLGziRnMmpuDBfTc+jQIIx/D++Mn8/lbO7r68vMmTNdj+Uy7RsR85TpYk0i3i4pI4f73/uJfadTaBoezH8e6UmNkOLdYioiUhaK+/mtUwwiZSQzx8GDC7aw73QKtUL9WTi2u8KBiFQYOsUgUgZyHE7++OFWYo5eIDTAh/ljutOwesGDuhwOB1u3bgWgc+fOBS5+5q20b0TMo4Ag4mZOp8Ezy3fy9d5E/H2szI3qRut6hR/Gy8zMpHv37oDGCf+W9o2IeXSKQcSNDMNg2uq9fLr1JDarhXeGd6Z70+pmtyUiUmIKCCJu9M76Q8z9IXeE8mtD2nNLZG2TOxIRKR0FBBE3+XDTMV77aj8Afx3Qmrs7NzC5IxGR0lNAEHGD1bsSeH7FLgAm9LmGcTc0NbkjEZGro4AgcpV+OHCOJz7ajmHAsO6NmHxbK7NbEhG5agoIIldh+/FLPLRoC9kOJ/3a1eGVQW2xWDRCWUQqPt3mKFJKBxNTGBMdQ3q2gxuah/Pm/R2xlWJ9BV9fX6ZMmeJ6LJdp34iYR6OWRUrh5KUMhvz7RxKSMunQIIwlD15LiL/ytoh4Po1aFikj51OzGDl3EwlJmVxTM5joMd0VDkSk0tF3NZESSM2yM2b+Zg6fTaNeWACLxvWgerDfVb2m0+lk7969AERGRmK1Krf/SvtGxDwKCCLFlGV38PCiLew8kUS1IF8WjutBvaqBV/26GRkZtG3bFtA44d/SvhExj+K4SDE4nAZPfLSdDQfPE+xnY/6Y7jSvFWJ2WyIiZUYBQeR3GIbBCyviWBN3Gj+blfdHdaVDw6pmtyUiUqYUEESKYBgGr67Zx9KYY1gtMGtoR65vHm52WyIiZU4BQaQIb687yHvfHQZg+uB23NGurskdiYiUDwUEkUJEbzjCP/77MwAv9I9kaPdGJnckIlJ+FBBECvDJluP87cs9AEy8pQXjb2xmckciIuVLtzmK/MbqXQk8s3wnAONuaMoTfVuU6fv5+voyefJk12O5TPtGxDwatSxyhfX7E3lw4RZyHAZDuzVkxt3ttPiSiFQqGrUsUkKbDp/nkcWx5DgMBrSvy7TBCgci4r10ikEE2HniEuMWbCEzx8nNEbVKvTJjaTidTo4dOwZAo0aNNE74Cto3IuZRQBCv9/OZFEbNiyE1y861zarzzvDO+NrK74MoIyODpk2bAhon/FvaNyLmURwXrxZ/Po0RH2ziUnoOHRpW5YOobgT42sxuS0TEdAoI4rUSkjIY/sEmElOyiKgTyoIx3bRss4jILxQQxCudS81ixAebOHExgyY1glg4rjtVg65u2WYRkcpEAUG8TlJGDqPmxnDobBr1wgJYPL4HtUIDzG5LRMSjKCCIV0nPtjN2/mb2JCQTHuLH4vE9aFAtyOy2REQ8jgKCeI3MHAcPLYwlNv4iVQJ8WDSuB81qhpjdloiIR9IVWeIVsu1O/vjhNn44eI4gPxvzx3Ynsq5nTAD18fHhsccecz2Wy7RvRMyjUctS6dkdTv60dBtr4k7j72MlenQ3rmsebnZbIiKm0KhlEXLDwaRlO1gTdxo/m5X3RnZROBARKQYds5NKy+E0ePo/O/lyxyl8bRb+PaIzN7WqZXZb+RiGwblz5wAIDw/X+g9X0L4RMY8CglRKTqfBs5/u5NNtJ7FZLfxrWGduiaxtdlsFSk9Pp1at3OCiccJ5ad+ImEenGKTSMQyDFz6PY9mWE1gtMGtoR/7Qto7ZbYmIVCgKCFKpGIbB377cw4ebjmGxwBv3dWRA+3pmtyUiUuEoIEilYRgG01btZf6PRwGYeU97BnWqb25TIiIVlAKCVAqGYTDzq/188MMRAKYPbse9XRua3JWISMWlgCCVwj+/PsC/1x8C4KW72vBAj0YmdyQiUrEpIEiFN/ubA8xaewCAvw5ozaieTcxtSESkEtBtjlKhvfftIf7x358BePaOCMbd0NTkjkrOx8eHqKgo12O5TPtGxDwatSwV1rwfjvDSyj0ATL6tJX+8uYXJHYmIeD6NWpZKbdFPR13h4PFbWigciIi4mY7ZSYXzUcwx/vr5bgAevekaJvWt2OHAMAzS09MBCAoK0jjhK2jfiJhHRxCkQlkac4y/fLoLgPE3NOXp21tV+A+N9PR0QkJCCAkJcX0YSi7tGxHz6AiCVBiLfjrqOnIw+romPN8/ssKHAxERT6WAIBVC9IYj/O3L3GsOHryxKc/1UzgQESlLCgji8eZ8d5hpq/cCudccVIbTCiIink4BQTzaO+sPMvP/9gPw+M3NmXRrS4UDEZFyoIAgHmvW1wd48+vcIUiT+rZkYgW/W0FEpCJRQBCPYxgGb/7vZ9765iAAT93eigl9mpvclYiId1FAEI/y66qMvy689Fy/CB7qdY3JXZUtm83GkCFDXI/lMu0bEfNo1LJ4DMMwmL56L3O+z12y+cUBrRlbAddWEBHxZMX9/NYRBPEIhmHw0so9RG84CsDLd7VhpFZlFBExjQKCmM7pNHjxizgWbzwGwPTB7XigRyOTuxIR8W4atSymcjoNnvtsF4s3HsNigZlD2ntdOEhLS8NisWCxWEhLSzO7HY+ifSNiHh1BENM4nAbPLN/Jf2JPYLXAP+7twN2dG5jdloiIoIAgJrE7nDz1n518tu0kNquFN+7rwF0d65vdloiI/EIBQcpdlt3BpI+3s3rXaXysFmYN7UT/9nXNbktERK6ggCDlKi3LziOLY/n+wDn8bFbeGtaJP7StY3ZbIiLyGwoIUm4upWczZv5mth27RJCfjfdHduWGFuFmtyUiIgVQQJBykZicyci5Mew/k0JYoC/zx3SjU6NqZrclIiKFUECQMhd/Po0Rczdx/EIGtUL9WTSuB63qhJrdlsew2Wz069fP9Vgu074RMY9GLUuZ2nc6mZFzYzibkkXjGkEsHteDhtWDzG5LRMRradSymC42/iJjomNIzrQTUSeUhWO7U6tKgNltiYhIMSggSJn47uezPLwolowcB50bVSV6dHfCgnzNbktERIpJo5bF7VbtTGDcgs1k5Djo1bImi8f3UDgoQlpaGsHBwQQHB2uc8G9o34iYR0cQxK2Wxhzj+c924TSgf/u6vHlfR/x8lEN/T3p6utkteCztGxFzKCCI27z77SFeXbMPgGHdG/HKoLbYrBaTuxIRkdJQQJCrZhgGf/+//bz77SEAHrvpGp66vRUWi8KBiEhFpYAgV8XhNHhhxS6WxhwH4Nk7Ini49zUmdyUiIldLAUFKLTPHwZPLchddslpg+uB2DO3eyOy2RETEDRQQpFQupmXz4MItbIm/iK8td0XGfu20IqOISGWhgCAldvRcGmPmb+bIuTRCA3x4b0QXrmuuRZdKy2q10rt3b9djuUz7RsQ8GrUsJRIbf5EHF27hQlo29asGMn9MN1rU1roKIiIVhUYti9ut2ZXAEx9vJ8vupF39MOaO7kqtUI1OFhGpjBQQ5HcZhsEH3x9h+pq9GAb0jazFW8M6EeSnvz4iIpWVTupJkewOJy9+vptpq3PDQVTPxrw3sqvCgRulpaVRs2ZNatasqXHCv6F9I2IefZeXQqVl2Xl86TbW7kvEYoHn+0Uy7oamGoBUBs6dO2d2Cx5L+0bEHAoIUqDE5EzGLthM3Mlk/H2szBrakT+01W2MIiLeQgFB8vn5TApjojdz8lIG1YP9+CCqK50bVTO7LRERKUcKCJLHjwfP8fDiWFIy7TQLDyZ6TDca1wg2uy0RESlnCgji8p/YE/xl+U7sToNuTarx/siuVAv2M7stERExgQKCYBgGs9Ye4J9fHwDgzg71eG1IewJ8bSZ3JiIiZlFA8HIZ2Q6eXr6TL3ecAnKXap58WyusVt2pUF6sVitdu3Z1PZbLtG9EzKNRy17s+IV0HloUy96EZHysFl4e1JZhWo1RRKRS06hlKdIPB87xx6VbuZSeQ3iIH+8M70L3ptXNbktERDyEAoKX+XVs8ow1e3Ea0KFBGO+O7ELdsECzWxMREQ+igOBFMrIdPLN8J1/8cr3BkC4NeGVQW12MaLL09HRat24NwJ49ewgKCjK5I8+hfSNiHgUEL3H8QjoPL4plzy/XG7x4Z2tGXttYY5M9gGEYxMfHux7LZdo3IuZRQPACGw6e448fbuVieg41gv14Z3hnejSrYXZbIiLiwRQQKjHDMJj7wxGmr8693qB9gzDeHdGFelV1vYGIiBRNAaGSysh28OynO1mxPfd6g3s6N2DaYF1vICIixaOAUAldeb2BzWrhr/0jibquia43EBGRYlNAqGR+PHiOCVdcb/D28M5cq+sNRESkhBQQKolf5xu8+n/7cDgN2tXPnW9QX9cbeDyLxeK6lU9HefLSvhExjwJCJXAuNYunPtnBuv1nAbi7c32mD26n6w0qiKCgIHbv3m12Gx5J+0bEPAoIFdz3B87y5LIdnE3Jws/Hyl/7RzJC8w1EROQqKSBUUNl2J6//bz/vfXsYgJa1Q3hrWCci6mjhLBERuXoKCBVQ/Pk0Hl+6jR0nkgAY3qMRL/RvTaCfTilUROnp6XTr1g2AzZs3a5zwFbRvRMyjgFDBfLbtBC98FkdatoOwQF/+fk97/tC2jtltyVUwDIM9e/a4Hstl2jci5lFAqCBSs+y8uCKOT7edBKB70+r88/6OmoooIiJlQgGhAthx/BKPf7SN+PPpWC3wRN+WTOjTHJtVFyKKiEjZUEDwYE6nwZzvD/PaV/uxOw3qVw1k1tCOdG1S3ezWRESkklNA8FCJyZn8+ZMdfH/gHAD929Vl+uB2hAX5mtyZiIh4AwUED7RuXyKTP9nB+bRsAnytTL2zDfd3a6jZBiIiUm4UEDxIUkYOr67Zy9KY4wBE1q3Cv4Z1pHmtUJM7k7JksVho3Lix67Fcpn0jYh4FBA/xf3EJ/PXz3ZxNyQJg9HVN+MsdERqX7AWCgoI4evSo2W14JO0bEfMoIJjsTHImL34ex1e7zwDQrGYwr97dnu5NdSGiiIiYRwHBJE6nwUebjzNj9V5Ssuz4WC08etM1TOjTXEcNRETEdAoIJjh0NpVnP91FzJELAHRoWJW/39NO6yh4qYyMDHr16gXAd999R2Cghl/9SvtGxDwKCOUox+Hk/e8OM2vtAbLtToL8bEy+rRVR1zXR0CMv5nQ62bJli+uxXKZ9I2IeBYRysuP4JZ5ZvpN9p1MA6NWyJtMGtaVhdS0+IyIinkcBoYylZ9t5478/M2/DEZwGVAvyZcqdbbirYz3dtiUiIh5LAaEMfffzWZ77bBcnLmYAMLhTfV7oH0mNEH+TOxMRESmaAkIZOH4hnde+2s8XO04BUL9qINMGt+WmVrVM7kxERKR4FBDcKCk9h7fXH2T+hqNkO5xYLDDmuqb8+baWBPtrV4uISMWhTy03yLI7WPhjPLPXHSQpIweA65vX4Nk7ImlbP8zk7qQiCA8PN7sFj6V9I2IOBYSr4HQafLnzFK99td91nUFEnVD+ckcEvVvW1EWIUizBwcGcPXvW7DY8kvaNiHkUEErpx0PnmLF6H7tOJgFQu4o/f76tFfd0bqCZBiIiUuEpIJTQz2dSeHXNPr7ZlwhAiL8Pj950DWOvb0qgn0Yki4hI5aCAUExnkjN5478/80nscZwG+FgtDO/RiMdvaaHbFuWqZGRkcMcddwCwZs0ajRO+gvaNiHkUEH5HSmYO7393mDnfHyYzJ3fU6x1t6/D0HyJoGh5scndSGTidTr799lvXY7lM+0bEPAoIhbiQls3ijfEs+PEo59OyAejSuBrP9YukS+NqJncnIiJSthQQfuPQ2VTm/nCE5bEnyLLn/sTSNDyYZ/4Qwe1tauvOBBER8QoKCIBhGGw8fIG5Pxzm672Jru3t6ocx/sam9GtXF1+b1cQORUREypdXB4Qch5NVOxP44IfDxJ1MBsBigVsiajP+xqb0aFpdRwxEKqGbbrqJJk2aMH/+fLNbEfFYXvljcVJGDu99e4heM9fxxMfbiTuZTICvleE9GrH2yd58ENWVa5vVUDgQ8TAJCQmEhIQwdOjQPNtXrlyJr68vzz//fJm9t91u57XXXqN9+/YEBgZisVjy/GrUqFGZvbeIGbzqCMLxC+lEbzjKx5uPkZbtACA8xJ+ono0Zfm1jqgf7mdyheKugoCCzW/BYV+6bunXr8vTTTzN16lSeeuopunTpwvr167n33nt59NFHmTZtWpn1MW7cOJYsWcLjjz/OzJkzOXPmDE8//TSJiYlMmjSJdu3aldl7i5ih0gcEh9Ng0+HzLNl0jDVxCTiN3O2taocy7samDOxQjwBfDTgS8wQHB5OWlub6fUJCAi1atGDAgAF89NFHru0rV65k8ODBPP3002X6QehJfrtvACZPnsx7773HM888w4wZMxg4cCDDhg1j1qxZBb6GYRg4HI582wzDwG6359nu41Pwt8QPP/yQhQsX8v777/Pggw/meZ0xY8Zw8803M2DAgNL8EUU8lsUwDKM0T0xOTiYsLIykpCSqVKni7r6uitNpsCX+Iit3nmL1rtOcS81yfe3GFuGMv7EZvVqE6xSCeKyXXnqJqVOnsnnzZtdPyXfccQcPPvggb731ltntmS46OpqxY8cSHBxM//79+fDDD7HZCg7669evp0+fPsV63SNHjtCkSZN823v37s3Zs2fZs2dPnu2bNm3i2muvZcGCBYwaNarEfw4RMxT387vSHEEwDINtxy+xckcCq3clcDo50/W1sEBf+rWrQ9R1TYio41lhRqQgJf0p2du0bNkSAIvFwvz58wsNBwBdunRh8+bNebY9/PDD1KtXjylTpuTZXq9evXzPv3jxIt9//z1PPvlkvq+dOHECgAYNGri2TZgwgczMTObOnYvT6WTw4ME0btxYwU4qnAodEAzDYNfJJFbtTGDlzgROXspwfS3U34fb2tRhQIe63NA8XLcpisfKzMzknnvuAWD58uUEBAQQFBTEK6+8wtixY+nTpw/9+/dnzpw5XnfUq6B9s337dgYMGMD111/Phg0bmDdvHhMmTCj0NUJDQ+natWu+bTVq1Mi3vSDHjx/HMIwCw8OKFSsIDw+nV69erm3PPvsskZGRTJkyhdmzZ+NwOHjzzTeL+0cW8RgVLiAYhsHehBRW7jzFql0JxJ9Pd30t2M9G39a1GdC+Hr1ahuPvo2sLxPM5HA5Wr17tevyrkvyUXFn9dt/s37+f22+/nZ49e/L5559z7733MnXqVEaMGEFYWFiZ9FC1alUA9u7dm2f7jz/+yNKlS5kyZUqeaxcaNGjAqFGjGDhwIAA//PCDV/6/k4qvQgSEpPQcth2/yJajF1kdl8Dhs5cvWgrwtXJLRG0GtK9Ln4hauuBQKoWS/pTsDeLj47n99ttp1aoVy5cvx9fXl1dffZW2bdsyffp0/v73v5fJ+zZq1IhevXoxf/58mjZtSvfu3YmJiWH69OncdtttBd5a2bFjR9555x1+/PFHQkJCyqQvkbLmcRcpOp0GB8+msjX+IluPXWTrsUscTEzNU+PnY+WmljUZ0KEet0TUIti/QuQckQKlpaW5PkRSU1M5ceIEvXr1okuXLq6fkjds2MDBgwfL7KdkT3XlvmnWrBlVq1Zl3bp1eb7nPPTQQyxcuJB9+/YVeIFhQUo6KCkxMZEnnniCNWvWkJ6eTsuWLRk9ejQTJ07Md+dDTEwMd999N9dddx3Vq1fn3XffLdZ7iJSX4n5+mx4QkjJy2H78kisQbD9+iZRMe766JjWC6NyoGje0COfW1rUJDfAt9XuKeJIrPwR3797N7bffTtOmTfnqq68IDAxk3759tG3blj//+c9l9lOyp/pteAoO9uwVVI8dO8YNN9zAwoULadWqFa1atWLXrl00btzY7NZEXDwuIKRm2UlMziQxJYv482lsjb/E1mMXOXg2ld92EOhro0PDMDo3qkbnRtXo1KgqNUL8S9OmiMcrq5+SK4OKFBBSUlK4/vrrmThxIuPGjQPgscceIycnhzlz5pjcnchl5RYQjp46S6bFn8SUTBKTs0hMycp9nJLF2eTLj9OzHYW+VuNfjg50blSVTo2qEVEnFB/ddSBeoiJ9CJY37RsR9yvzOQi/5orrX16F1b94Y2KD/KzUDA2gTpUA2tQPo2PDqnRoEJbv6EB6WmohryBS+Vw5KTA5OTnf1D9vpn0j4n7JybmLE/7e8YFSH0E4ceIEDRs2LM1TRURExGTHjx/PM+Trt0odEJxOJ6dOnSI0NNTrhreIiIhUVIZhkJKSQr169bBaCz+dX+qAICIiIpWXrgQUERGRfBQQREREJB8FBBEREclHAUFERETyUUAQERGRfBQQREREJB8FBBEREclHAUFERETyUUAQERGRfBQQREREJB8FBBEREclHAUFERETy+X/8P5T9EKyUZAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "x = np.linspace(0, 10)\n", "x1 = 3\n", "sigma1 = 1\n", "\n", "def f(x):\n", " f = -x**3 + 10*x**2 + x\n", " return f\n", "\n", "plt.figure()\n", "plt.plot(x, f(x))\n", "\n", "plt.plot([x1, x1],[0, f(x1)], 'k--')\n", "plt.plot([x1 + sigma1, x1 + sigma1],[0, f(x1 + sigma1)], 'k--')\n", "\n", "plt.plot([0, x1],[f(x1), f(x1)], 'k--')\n", "plt.plot([0, x1 + sigma1],[f(x1 + sigma1), f(x1 + sigma1)], 'k--')\n", "\n", "\n", "plt.xlim([0, 7])\n", "plt.ylim([0, 160])\n", "\n", "plt.text(0, f(x1), '$f(x)$', fontsize = 12, verticalalignment = 'bottom')\n", "plt.text(0, f(x1 + sigma1), '$f(x + \\sigma_x)$', fontsize = 14, verticalalignment = 'bottom')\n", "\n", "plt.text(x1, 0, '$x$', fontsize = 12, verticalalignment = 'bottom')\n", "plt.text(x1 + sigma1, 0, '$x + \\sigma_x$', fontsize = 12, verticalalignment = 'bottom')\n", "\n", "plt.xticks([]) \n", "plt.yticks([]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\\frac{df}{dx} \\approx \\frac{f(x + \\sigma_x) - f(x)}{\\sigma_x} $$\n", "\n", "$$ \\sigma_x \\frac{df}{dx} \\approx f(x + \\sigma_x) - f(x)$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## How error propagates through calculations and anlysis\n", "\n", "Now imagine you have three variables:
\n", "$a$, $b$, $c$
\n", "Each of those variables has a measure of uncertainty:
\n", "$\\sigma_a,\\sigma_b,\\sigma_c$\n", "\n", "We want to know the uncertainty of a calculated variable $y = f(a,b,c)$ \n", "\n", "### Oceanographic example\n", "\n", "An oceanographic example would be density, which is usually expressed as a function of temperature, salinity and pressure. Oceanographers measure these parameters with a CTD (conductivity, temperature and depth) instrument. This instrument does measure salinity directly. Salinity is calculated as a function of conductivity, temperature and pressure. This means that errors in salinity are not independent of errors in temperature. When estimating the uncertainty in density, we have to think about how the errors in salinity and temperature co-vary.\n", "\n", "### General equation for uncertainty\n", "\n", "What is the uncertainty of y when measurements are dependant/inter-related?\n", "\n", "$\\sigma_y^2 \\approx \\sigma_a^2 \\left(\\frac{\\partial y}{\\partial a}\\right)^2 +\\sigma_b^2 \\left(\\frac{\\partial y}{\\partial b}\\right)^2 + \\sigma_c^2\\left(\\frac{\\partial y}{\\partial c}\\right)^2 + 2\\sigma_{ab} ^2 \\left(\\frac{\\partial y}{\\partial a}\\right)\\left(\\frac{\\partial y}{\\partial b}\\right) + 2\\sigma_{bc}^2 \\left(\\frac{\\partial y}{\\partial b}\\right)\\left(\\frac{\\partial y}{\\partial c}\\right) + 2\\sigma_{ac}^2 \\left(\\frac{\\partial y}{\\partial a}\\right) \\left(\\frac{\\partial y}{\\partial c}\\right)$\n", "\n", "This equation shows that the uncertainty in the calculated variable $y$ depends on:\n", "\n", "* The uncertainty in the measurement variables ($\\sigma_a,\\sigma_b,\\sigma_c$)\n", "\n", "* The covariance among measurment variables ($\\sigma_{ab}^2,\\sigma_{bc}^2,\\sigma_{ac}^2$)\n", "\n", "* The relationship between magnitudes of the measurement variables and calculated variable ($\\partial y$/$\\partial a$, $\\partial y$/$\\partial b$, $\\partial y$/$\\partial c$)\n", "\n", "### Simplified equation for uncertainty, if errors are uncorrelated\n", "\n", "If the errors are uncorrelated, the covariance terms are zero and the above equation simplifies to\n", "\n", "$\\sigma_y^2 \\approx \\sigma_a^2 \\left(\\frac{\\partial y}{\\partial a}\\right)^2 +\\sigma_b^2 \\left(\\frac{\\partial y}{\\partial b}\\right)^2 + \\sigma_c^2\\left(\\frac{\\partial y}{\\partial c}\\right)^2$\n", "\n", "## Special cases
\n", "\n", "Rules for error propagation in common operations can be derived from the equation above. If there all measurement variables are independent, then the terms involving covariance become zero. The examples listed below assume no correlation between measurement variables.\n", "\n", "A more complete table, including covariance terms, can be found on Wikipedia:\n", "https://en.wikipedia.org/wiki/Propagation_of_uncertainty#Example_formulas\n", "\n", "\n", "### Addition, subtractions \n", "\n", "$y = a + b - c$\n", "\n", "_error:_\n", "\n", "$\\sigma_y ^2 = \\sigma_a ^2(1)^2 + \\sigma_b ^2(1)^2 + \\sigma_c ^2(-1)^2$\n", "\n", "$\\sigma_y ^2 = \\sigma_a ^2 + \\sigma_b ^2 + \\sigma_c ^2$\n", "\n", "### Multiplication, division\n", "\n", "$y = \\frac{ab}{c}$\n", "\n", "_error_:\n", "\n", "$\\sigma_y ^2 = \\sigma_a ^2\\left(\\frac{b}{c}\\right)^2 + \\sigma_b ^2\\left(\\frac{a}{c}\\right)^2 + \\sigma_c ^2\\left(\\frac{-ab}{c^2}\\right)^2$\n", "\n", "_relative error:_\n", "\n", "$\\left(\\frac{\\sigma_y }{y}\\right)^2 = \\left(\\frac{\\sigma_a }{a}\\right)^2 + \\left(\\frac{\\sigma_b }{b}\\right)^2 + \\left(\\frac{\\sigma_c }{c}\\right)^2$\n", "\n", "### Multiplying by an exact number\n", "\n", "$y=ma$\n", "\n", "_error:_\n", "\n", "$(\\sigma_y) = m(\\sigma_a)$\n", "\n", "_relative error:_\n", "\n", "$\\left(\\frac{\\sigma_y}{y}\\right)^2 = \\left(\\frac{\\sigma_a}{a}\\right)^2$\n", "\n", "### Power\n", "\n", "$y = a^m$\n", "\n", "$(\\frac{\\sigma_y}{y})^2 = (m \\frac{\\sigma_a}{a})^2$\n", "\n", "### Natural log \n", "\n", "$y= \\ln{a}$\n", "\n", "$\\sigma_y = |\\frac{\\sigma_a}{a}|$\n", "\n", "### Log, base 10\n", "\n", "$y = \\log_{10}{a}$\n", "\n", "$\\sigma_y = |0.434 \\frac{\\sigma_a}{a}|$\n", "\n", "### Exponential\n", "\n", "$y = e^a$\n", "\n", "$\\sigma_y = |y \\sigma_a|$\n", "\n", "### Exponential, base 10\n", "\n", "$y = 10^a$ \n", "\n", "$\\sigma_y = |2.303 y \\sigma_a|$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Python uncertainties package\n", "\n", "The [uncertainties](https://pythonhosted.org/uncertainties/) package in Python makes it easy to keep track of uncertainties in complicated equations and account for correlations, which are ignored in the special cases above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Example__: Measure diameter of a sphere:\n", "\n", "relative error = 10%\n", "\n", "Error in calculating volume?\n", "\n", "$V = \\frac{4}{3} \\pi r^3$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This problem can be approached by breaking it into several parts. First, we have the relative error in measuring the diameter, what is the error in calculating the radius, \n", "\n", "$r = \\frac{1}{2}d$?\n", "\n", "The _absolute error_ changes when calculating the radius,\n", "\n", "$$\\sigma_r = \\frac{1}{2} \\sigma_d$$\n", "\n", "but the _relative error_ does not,\n", "\n", "$$\\frac{\\sigma_r}{r} = \\frac{\\sigma_d}{d}$$\n", "\n", "Similarly, the relative error of the volume, $V$, is equal to the relative error in calculating $r^3$,\n", "\n", "$$\\frac{\\sigma_V}{V} = \\frac{\\sigma_{r^3}}{r^3}$$\n", "\n", "The error in $r^3$ can be found by using the rule for powers,\n", "\n", "$$\\left(\\frac{\\sigma_{r^3}}{r^3}\\right) = 3 \\left(\\frac{\\sigma_r}{r}\\right)$$\n", "\n", "which is equivalent to \n", "\n", "$$\\left(\\frac{\\sigma_{V}}{V}\\right) = 3 \\left(\\frac{\\sigma_d}{d}\\right)$$\n", "\n", "$$\\left(\\frac{\\sigma_{V}}{V}\\right) = 3 \\left(0.1\\right) = 0.30$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This result can be checked with the uncertainties package." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.029999999999999995" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from uncertainties import ufloat\n", "import numpy as np\n", "\n", "d = ufloat(10,0.1)\n", "V = (4/3)*np.pi*(0.5*d)**3\n", "V.std_dev/V.nominal_value" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Application: Marine Carbonate Chemistry\n", "\n", "![images/carbonate_example.png](images/carbonate_example.png)\n", "\n", "[source](http://acmg.seas.harvard.edu/people/faculty/djj/book/bookchap6.html)\n", "\n", "Relative concentrations of dissolved CO2 and carbonate ions in seawater are determined by pH pf seawater and dissociation constants $K_1'$ and $K_2'$\n", "\n", "DIC = dissolved inorganic carbon\n", "\n", "DIC = [CO2] + [HCO3-] + [CO3-2]\n", "\n", "\n", "$CO_2 + H_2 O \\iff HCO_{3}^- + H^+$\n", "\n", "$k_1'= \\frac{[HCO3-][H^+]} {[CO_2]}$\n", "\n", "$HCO_3 ^- \\iff CO_3 ^2- +H^+$\n", "\n", "$k_2' =\\frac{ [CO3- ^2-][H^+]}{[HCO_3 ^-]}$\n", "\n", "\n", "partial pressure of CO2\n", "\n", "$pCO2 = \\frac{[CO2]}{K_{h,CO_2}}$\n", "Henry's law constant = 3.24X10$^-2$ $\\frac{mol}{kg atm}$ at 20 degrees C, s=35 PSU\n", "\n", "$pK_1' = 5.847$\n", "\n", "$pK_2' = 8.966$\n", "\n", "$pH = -\\log_{10}{[H^+]} \\to [H^+] = 10^{-pH}$\n", "\n", "$pK_1' = -\\log_{10}{K_1'} \\to K_1' = 10^{-pK}$\n", "\n", "We can measure pCO_2 and pH well, but not DIC, but not alkalinity\n", "\n", "$[CO_2] = \\frac{DIC} {1 + \\frac{k_1'}{H^+} + \\frac{k_1' k_2'}{[H^+]^2}}$\n", "\n", "\n", "Analytical error in pH = 0.0020
\n", "Analytical error in DIC = umol kg-1\n", "\n", "What is the $\\sigma_{H^+}$?\n", "\n", "$\\sigma_{H^+} = 2.303[H^+]\\sigma_{pH}$\n", "\n", "To get the order of magnitude assume a pH around 8.2 (typical for seawater)\n", "\n", "$[H^+] \\approx 10^{-8.2}$\n", "\n", "\n", "Chemical oceanographers may be interested in knowing the error in calculating the partial pressure of CO2 in seawater from these measurements. The PYCO2SYS package provides [tools for propagating uncertainty](https://pyco2sys.readthedocs.io/en/latest/uncertainty/) in its carbonate system calculations." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.6" }, "vscode": { "interpreter": { "hash": "0ef88d3abb6b62f34a20525ce337090c4512fe8aecf32c74604482b944e1c3bd" } } }, "nbformat": 4, "nbformat_minor": 4 }