{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Example applications of F and chi-square tests of variance\n", "Elizabeth Maroon, 7/16/2021
\n", "Adapted from AOS 575 materials\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.stats import f, norm, chi2\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### F-statistic\n", "\n", "The F-statistic can be used to compare two sample standard deviations $s_1$ and $s_2$ from samples of size $N_1$ and $N_2$. It assumes that the two samples are taken from underlying normal populations. If the populations are not normal, the F-test will likely fail because it is *very sensitive to normality*. Wilks states that a *permutation test* is a more reliable way to test sample differences for a complicated statistic like variance and isn't any harder to compute than an f-test (Wilks textbook, pp169-172). \n", "\n", "Random variable $F$,\n", "\n", "\\begin{equation}\n", "F = \\frac{s_1^2}{s_2^2},\n", "\\end{equation}\n", "\n", "is compared to the F-distribution using the parameters $\\nu_1$ and $\\nu_2$, the degrees of freedom for the two samples." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "#constructing two test time series of the same length\n", "x1 = np.random.normal(loc=0.0, scale=1.0, size=100)\n", "x2 = np.random.normal(loc=0.0, scale=2, size=100) \n", "#x2 has standard deviation 2x that of x1 (scale) by design\n", "#time series are not required to be same length - dof below accounts for that\n", "#loc = mean for distribution, here centered at 0." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "sample 1 mean: 0.16733089871925952\n", "sample 1 variance: 0.829555502045668\n", "\n", "sample 2 mean: -0.0022216687028138614\n", "sample 2 variance: 4.42306495474155\n" ] } ], "source": [ "#calculate means and variances\n", "mu1 = np.mean(x1)\n", "mu2 = np.mean(x2)\n", "\n", "var1 = np.var(x1, ddof = 1) #setting DOF to N-1 b/c sample variance not pop variance\n", "var2 = np.var(x2, ddof = 1) #setting DOF to N-1 b/c sample variance not pop variance\n", "\n", "print('sample 1 mean:', mu1)\n", "print('sample 1 variance:', var1)\n", "print()\n", "print('sample 2 mean:', mu2)\n", "print('sample 2 variance:', var2)\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAArIElEQVR4nO3de5yMdf/H8ddnDrubc6GQw0rkGMmhIomckkOhSEoHUtGvu4huktKBkuq+U1IhlSQkISREIqcohw4SWRQWsdbaOXx/f8zmXmyZXTPznZn9PB+PHreZuWav97h33679Xtf1/YoxBqWUUrHPYTuAUkqp0NBCV0qpOKGFrpRScUILXSml4oQWulJKxQmXrR2XKFHCJCcn29q9UkrFpLVr1+43xpTM6TVrhZ6cnMyaNWts7V4ppWKSiOz4u9d0yEUppeKEFrpSSsUJLXSllIoT1sbQlVKxw+PxkJKSQkZGhu0o+UZSUhJly5bF7XYH/R4tdKXUGaWkpFC4cGGSk5MREdtx4p4xhtTUVFJSUqhYsWLQ7wtqyEVEWovIjyKyVUQG5fB6UxH5U0TWZ/03NBfZlVJRLiMjg+LFi2uZR4iIULx48Vz/RnTGI3QRcQJjgBZACrBaRGYZYzafsukyY8wNudq7UipmaJlHVl7+voM5Qm8AbDXGbDPGZAJTgA653pNSSqmwCqbQLwR2ZnuckvXcqa4UkQ0i8pmI1AhJOqWUsmjo0KEsXLgwJF9rwIAB1KhRgwEDBjB27FgmTZoEwMSJE9m9e3dI9hHMSdGcjvtPXRVjHVDBGJMmItcDM4HKp30hkd5Ab4Dy5cvnLqlScSh50Bwr+90+oq2V/caap556KmRf64033mDfvn0kJiae9PzEiROpWbMmZcqUOet9BHOEngKUy/a4LHDSPyfGmMPGmLSsP88F3CJS4tQvZIwZZ4ypZ4ypV7JkjlMRKKVUjoYPH07VqlVp0aIF3bp1Y9SoUQC8+eab1K9fn9q1a9OpUyfS09MB6NmzJ9OmTTvx/kKFCgGwZ88emjRpQp06dahZsybLli3D5/PRs2dPatasSa1atXjppZdO+xpPPfUU9evXp2bNmvTu3Zu/Vntr2rQpAwcOpEGDBlSpUoVly5adlr19+/YcPXqUhg0b8uGHHzJs2DBGjRrFtGnTWLNmDd27d6dOnTocO3bsrP6OgjlCXw1UFpGKwC6gK3Br9g1EpBTwhzHGiEgDAv9QpJ5VMqVUdPpsEPz+fWi/Zqla0GbE3768Zs0apk+fzrfffovX66Vu3bpcfvnlANx000306tULgCFDhvD222/Tr1+/v/1akydPplWrVgwePBifz0d6ejrr169n165dbNy4EYBDhw6d9r6+ffsydGjgAr4ePXowe/Zs2rVrB4DX62XVqlXMnTuXJ5988rRhmlmzZlGoUCHWr18PwLBhwwDo3Lkzr776KqNGjaJevXpn/ns6gzMWujHGKyJ9gfmAExhvjNkkIn2yXh8LdAbuExEvcAzoanSxUqVUiHz11Vd06NCBc845B+BEkQJs3LiRIUOGcOjQIdLS0mjVqtU/fq369etz11134fF46NixI3Xq1OGiiy5i27Zt9OvXj7Zt29KyZcvT3rd48WKef/550tPTOXDgADVq1DiR46abbgLg8ssvZ/v27SH61LkX1I1FWcMoc095bmy2P78KvBraaEqpqPQPR9Lh8k/Hhz179mTmzJnUrl2biRMnsmTJEgBcLhd+v//E+zMzMwFo0qQJS5cuZc6cOfTo0YMBAwZw++23s2HDBubPn8+YMWOYOnUq48ePP7GPjIwM7r//ftasWUO5cuUYNmzYSdeI/zUu7nQ68Xq9of74QdO5XJRSUa9x48Z8+umnZGRkkJaWxpw5/zuZfOTIEUqXLo3H4+H9998/8XxycjJr164F4JNPPsHj8QCwY8cOzj//fHr16sXdd9/NunXr2L9/P36/n06dOjF8+HDWrVt30v7/Ku8SJUqQlpZ20tj82SpcuDBHjhwJydfSW/+VUlGvfv36tG/fntq1a1OhQgXq1atH0aJFgcDJ0oYNG1KhQgVq1ap1ohx79epFhw4daNCgAc2bN6dgwYIALFmyhBdeeAG3202hQoWYNGkSu3bt4s477zxxRP/cc8+dtP9ixYrRq1cvatWqRXJyMvXr1w/ZZ+vZsyd9+vThnHPOYcWKFSeGlfJCbA1116tXz+gCFyq/i5XLFrds2UK1atXClCY4aWlpFCpUiPT0dJo0acK4ceOoW7eu1UzhltPfu4isNcbkeAZVj9CVUjGhd+/ebN68mYyMDO644464L/O80EJXSsWEyZMn244Q9fSkqFJKxQktdKWUihNa6EopFSe00JVSKk7oSVGlVK6F+nLLaJv9sWnTprmaX+Wjjz5i2LBhbNmyhVWrVoVkXpa80CN0pZQ6SzVr1mTGjBk0adLEag4tdKVU1Dt69Cht27aldu3a1KxZkw8//BD45ylt//Wvf9GkSROqVavG6tWruemmm6hcuTJDhgwBYPv27VStWpU77riDSy+9lM6dO5+Yeje7BQsWcOWVV1K3bl26dOlCWlraadtUq1aNSy65JIx/A8HRQldKRb158+ZRpkwZNmzYwMaNG2ndujUQmNJ29erVbNy4kWPHjjF79uwT70lISGDp0qX06dOHDh06MGbMGDZu3MjEiRNJTQ3M7v3jjz/Su3dvvvvuO4oUKcJrr7120n7379/P008/zcKFC1m3bh316tVj9OjRkfvguaSFrpSKerVq1WLhwoUMHDiQZcuWnZjHZfHixTRs2JBatWqxaNEiNm3adOI97du3P/HeGjVqULp0aRITE7nooovYuTOwqma5cuVo1KgRALfddhtfffXVSftduXIlmzdvplGjRtSpU4d33nmHHTt2ROIj54meFFVKRb0qVaqwdu1a5s6dy2OPPUbLli159NFHg5rS1uFwnLTsm8PhODHFrcjJK2ye+tgYQ4sWLfjggw/C9dFCSo/QlVJRb/fu3RQoUIDbbruN/v37s27dupBMafvbb7+xYsUKAD744AMaN2580utXXHEFy5cvZ+vWrQCkp6fz008/neWnCR89QldK5VqkLzP8/vvvGTBgAA6HA7fbzeuvvx6SKW2rVavGO++8w7333kvlypW57777Tnq9ZMmSTJw4kW7dunH8+HEAnn76aapUqXLSdh9//DH9+vVj3759tG3bljp16jB//vy8f+A80ulzlbJIp8+1Z/v27dxwww0n1hGNRrmdPleHXJRSKk7okItSkeI5Bhl/QsZhOH4YMg7R1rGSRDL51lTmV1PadsJ8JTk5OaqPzvNCC12pcPJmwg+zYc142L7stJfHJPzvz7/6L2Cx/zIW++vwjb8ambgjGPTMjDGnXQWiwicvw+Fa6EqFw8EdsHYifPsuHN0HxcrD1f2hSGlILApJRSCxCK3GrsePgyscm2nm+JZbnV9wl2se6SaR5f6ajPe1ZoW/hu1PQ1JSEqmpqRQvXlxLPQKMMaSmppKUlJSr92mhKxVKO1fD0ufh589BBKq0hnp3Q6Vm4Dj9lNWP5gAAP/vK8q6vJUkczyr39bR0ruED5zN85G3CM97uHKJwpD/NCWXLliUlJYV9+/ZZy5DfJCUlUbZs2Vy9RwtdqVDw+2H5S7DoGShYAq55FOreDkVz9wOZQSJL/JexxH8Zz3i708/1Mfc6Z9PM+S3DPT2Y6W8ERP4I2e12U7FixYjvV+WOFrpSZyttL3x8L/yyCGrcBO1eCQypnKXjJDDKewuf+q7kOfdbvJzwGjf5ljHYexc7zQUhCK7ijV62qNTZ2LYExjaGHV8Hirzz+JCUeXY/mvJ0zhzGEM+dXObYyoKEgdzm/Dyk+1DxQQtdqbzweQPDK5M6QlIx6LUILu8ZGDcPAz8O3vO14LrjL7DCX52n3RN40DkDsHNjoIpOOuSiVG55jsEH3WDbYqhzG1z/PCQUjMiu/+A8enkeYSRv8rB7GgXlGM95b8XGuLqKPlroSuWGNxM+7BEYamn3H7j8johH8OFkgKc3R00i97rmUJAMHvfeidFfuPM9LXSlguXzwvS7YevngfFyC2X+F4ODJ7w9Oco53O+aRQE5zgDPvfhwWsuk7NNCVyoYfj988gBsmQWtnguMl1snPO/tSppJ4lH3VApwnAc9faPuDlMVOfo7mlJnYgzMfQS+mwLXDoEr77ed6CSv+TrypKcHrZ2recM9Gic+25GUJVroSv0TY2DBkMBcLI0egib9bSfK0QRfGwZ77uJa5wYGu963HUdZElShi0hrEflRRLaKyKB/2K6+iPhEpHPoIipl0ZcjYcWrUL8XXDcsbJclhsL7vusY723NXa55dHEusR1HWXDGQhcRJzAGaANUB7qJSPW/2W4kEPllOpQKh00fw5LnoPat0Ob5qC7zvzzj7c5yXw2edo2njmy1HUdFWDBH6A2ArcaYbcaYTGAK0CGH7foB04G9IcynlB2pv8An/aBs/cAVLTlMrBWNfDh5wPMgf5hzeSNhNOdz0HYkFUHBfJdeCOzM9jgl67kTRORC4EZg7D99IRHpLSJrRGSNztqmopbnGEy9A5wu6DwBXAlnfk8UOURhenseoTDHGJvwEgl4bEdSERJMoef0e+ap9xu/DAw0xvzj6XVjzDhjTD1jTL2SJUsGGVGpCPtsIPzxPdw4DoqVs50mT34w5XnYcx91HVsZ7pqAThGQPwRzHXoKkP27uiyw+5Rt6gFTsia+LwFcLyJeY8zMUIRUKmI2fAjr3oHGD0OVlrbTnJV5/gb8x9uRB10z2WQqMMnXynYkFWbBFPpqoLKIVAR2AV2BW7NvYIw5MVGyiEwEZmuZq5iz9weY/RBUaATXDradJiRe8nammvzGUNe7bPBXYoO52HYkFUZnHHIxxniBvgSuXtkCTDXGbBKRPiLSJ9wBlYqIzKPw0R3gLgCd3g6Mn8cBg4OHPffzB+fyonssiWTajqTCKKhT98aYucaYKsaYSsaYZ7KeG2uMOe0kqDGmpzFmWqiDKhU2xsDsh2Hfj9DprcC6n3HkCAUY6OnNxY7d9HdNtR1HhVFsXIulVDhtnB64rb/pIKh0re00YfGVvxbveZtzt/Mz6ssPtuOoMNFCV/nb0VT47FG48HJoMsB2mrB61tudFFOCUe6xcDzNdhwVBlroKn+bNwgyDkP7V8ER31PPppNEf08fysk+WPiE7TgqDLTQVf710wL4fipc/QhccNpsFnFplanGBF9rWP1WYJEOFVfi41S+Url1/AjM/heUrApXP2w7TUQ9772Fpo71JL5zN22Oj+AIBSK27+0j2kZsX/mRHqGr/Gnhk3B4F7T/L7gSbaeJqOMk8IjnPkqTyhDXe7bjqBDSQlf5z44VgSGHhn2gXAPbaaxYby7mDV87bnEt4RrHBttxVIhooav8xZMBs/pB0XLQbIjtNFa97O3EL/7SDHNN1Am84oQWuspflr4AqT9Du5chsZDtNFZl4maY9w4qOv6gl3OO7TgqBLTQVf7x+0ZY/nJgwYqLm9tOExWW+S/lM199+rpmciE6pXWs00JX+YMxMLc/JBWFVs/YThNVhnt6ADDErSdIY50WusofNs2A31ZA86FQ4DzbaaLKbkrwX29H2jhX00RPkMY0LXQV/zLTYcFQKHUpXNbDdpqo9JavLdv8pRjmekdPkMYwLXQV/5a/AodToM3IuL+9P68ycfOk9w4ucvzOPc65tuOoPNJCV/Ht0M7AidAaN0GFq2yniWpf+msz31ePvq6ZlGG/7TgqD7TQVXz7fCgg0OIp20liwnBvDxz4GawnSGOSFrqKX9uXB06GNn4oZhd7jrQUU5JXvR1p61xFY8f3tuOoXNJCV/HJ74N5A6FIWbjqQdtpYsqbvrbs8J/PENd7OPDbjqNyQQtdxadv34Xfv4eWT0FC5GYTjAfHSWCktytVHTvp5FxqO47KBS10FX+OHYIvnoLyVwVOhqpcm+tvyDr/xTzi+ohzyLAdRwVJC13Fn6UvQPoBaDMCRGyniVHCM57ulJKD3O38zHYYFSQtdBVfDm6Hb96Ay7pD6dq208S0teYS5vnq08f1KSX403YcFQRdsUhFneRBeZ/5b7T7Na53QNMVDfl9Re6+jq6mc7qR3q40T1jH/7mm87j3Lttx1BnoEbqKG9VlOx0dyxnva8PvFLcdJy78akoz2deMbs5FVJJdtuOoM9BCV3FjoGsKhynAWG8721HiyiveThwjkYGuKbajqDPQQldx4SrHRq5xfser3o4cpqDtOHHlAEV43duOls611JcfbMdR/0ALXcU8wc8g1wekmBK862thO05cGu9rwx5zHoPd7wPGdhz1N7TQVcxr6/iGSx2/MtrTmeMk2I4TlzJI5EVvF+o4fuEGx0rbcdTf0EJXMc2NlwGuD9niL8dMf2PbceLaDN/VbPGXZ4DrQ9x4bcdROdBCVzGtm/MLKjj2MtLbDb9+O4eVHwcjvbdQwbGXm51LbMdROdDr0FXMKsgxHnR9zApfdZb4z/4morO5/j2/WOKvw2p/FR50zWC672oySLQdSWWjhzQqZvV2zaGEHGaEtyugt/hHhvC8pysXyCFudy6wHUadIqhCF5HWIvKjiGwVkUE5vN5BRL4TkfUiskZEdDBThVVx/uQe5xzm+BqwwVxsO06+stpUZYmvNve5PqUw6bbjqGzOWOgi4gTGAG2A6kA3Eal+ymZfALWNMXWAu4C3QpxTqZPc55pFEpmM9naxHSVfesF7M+dKGve4dJgqmgRzhN4A2GqM2WaMyQSmAB2yb2CMSTPG/HVxakH0QlUVRqVIpYdzITN8V/OLudB2nHxpk6nIbF9D7nHOpbhO3BU1gin0C4Gd2R6nZD13EhG5UUR+AOYQOEo/jYj0zhqSWbNv37685FWKfq6ZCH5e8XWyHSVfG+3tQhKZ3O+aZTuKyhJMoed0tum0I3BjzMfGmKpAR2B4Tl/IGDPOGFPPGFOvZMmSuQqqFEB5+YObnUuY4mtGitHvIZu2mTJM813Dbc7PKcN+23EUwRV6CpB9hd2ywO6/29gYsxSoJCIlzjKbUqd5yDUdL07+6+1oO4oCXvEGVoR60DXDchIFwRX6aqCyiFQUkQSgK3DS71gicrFIYGkYEakLJACpoQ6r8rfKkkJHx3Le8bVkH+fajqOA3ZTgfd91dHYu5SL52+M8FSFnLHRjjBfoC8wHtgBTjTGbRKSPiPTJ2qwTsFFE1hO4IuaWbCdJlQqJh10fcZQknR43yozxduA4bh52TbMdJd8L6k5RY8xcYO4pz43N9ueRwMjQRlPqf2rKNto4V/Oy9yYOUdh2HJVNKkV529eGB10zec3bns0m2XakfEvvFFUxob/rIw6aQrzlvd52FJWDt7xt+dMU4F+u6baj5Gta6Crq1ZMfaOrcwOvedqRRwHYclYPDFGSc9wZaONdSW7bajpNvaaGrKGcY4J7KXlOMSb6WtsOofzDR14pUU5hHXB/ZjpJvaaGrqNbYsZGGjh/4r7ejzuwX5Y5yDmO97Wji/F6XqrNEC11FMUN/11RSTAk+9F1rO4wKwru+Fuw1xejvnorOABJ5WugqajV3rKOO4xf+472RTNy246ggZJDIq94ONHT8QCPHRttx8h0tdBWVBD8Pu6bxq/8CZviuth1H5cIUXzN2meL0d32EHqVHlha6ikqtHaup4djBK95OeHVhrZiSiZv/em/kMsdWmjm+tR0nX9FCV1HHkXV0/rP/Qmb5r7IdR+XBNF8TdvjP5xHXRwh+23HyDS10FXXaOb6msmMXL3k76cLPMcqLi5e9najh2EErxxrbcfIN/WlR0cXn4SHXdDb7K/CZv4HtNOosfOJvxFZ/GR52fYRDj9IjQgtdRZcNH1DR8Qcvejtj9NszpvlxMNrbmSqOXbR3fG07Tr6gPzEqeniPw5fPs95fiS/8dW2nUSHwmb8Bm/0VeMg1HRde23Hinha6ih7rJsGfO3nR24WcF8pSscbg4EVvZ5Idf9DJucx2nLinha6ig+cYLB0F5a9imb+W7TQqhL7w1+Vb/8WBVY28x23HiWta6Co6rH4b0n6HZoPRo/N4I4zyduFCSYW1E22HiWta6Mq+40fgq9Fw0bWQ3Nh2GhUGy/01WeGrHvgtLPOo7ThxSwtd2bfydUhPhWaP206iwiZwlM7RvbDqTdth4pYWurIr/QB8/V+oegOUvdx2GhVGa80lcHELWP4yZBy2HScuaaEru77+T2DI5drBtpOoSGg2GI4dhJWv2U4Sl7TQlT1H/oCVY6FWZ7iguu00KhLKXAbV2sGKMYHfzlRIaaEre5a9CL5MaPqY7SQqkq4dHPitbPkrtpPEHS10Zceh32DNeLjsNiheyXYaFUnnV4NaXeCbNwK/pamQ0UJXdnw5EsQB1zxqO4myoemgwG9ny0bZThJXtNBV5O3/GdZ/APXvhqJlbadRNhSvBHVvhzUT4OB222nihha6irzFz4IrCRo/bDuJsumageBwwuLnbCeJG1roKrJ+/x42zYAr7oNCJW2nUTYVKQ0N+8B3H8LvuqB0KGihq8j6YjgkFYWr+tlOoqJB44cgqQgsGm47SVzQQleRs/0r+Hl+YKjlnGK206hocM650Ogh+Gke7FhhO03M00JXkWEMfP4EFLkQGt5rO42KJg37QKFSsHBY4PtE5ZkWuoqMLZ/CrjVw7b/BfY7tNCqaJBQIXL66cyX8vMB2mpimha7Cz+eFL56EktWgdjfbaVQ0qns7nFsRFj4Jfl1QOq+00FX4ffsupG6F654IXKam1Kmcbmg2BPZugo3TbKeJWUEVuoi0FpEfRWSriAzK4fXuIvJd1n9fi0jt0EdVMSnzKCwZAeWvhCqtbadR0azGTVCqFix6GryZttPEpDMWuog4gTFAG6A60E1ETp0a71fgGmPMpcBwYFyog6oYtfL1wNJy1z0JokvLqX/gcEDzYXBohy5Vl0fBHKE3ALYaY7YZYzKBKUCH7BsYY742xhzMergS0Pu5FRxNDcyoV/UGKN/QdhoVCy5uDslXw5cjIONP22liTjCFfiGwM9vjlKzn/s7dwGc5vSAivUVkjYis2bdvX/ApVWxa9iJkpkHzobaTqFghAi2HB5YkXDbadpqYE0yh5/R7co4Xi4rItQQKfWBOrxtjxhlj6hlj6pUsqbd9x7WDO2D1m4HpcUteYjuNiiVlLoNLuwaG6w7usJ0mpgRT6ClAuWyPywK7T91IRC4F3gI6GGNSQxNPxazFzwSmx9XFK1ReNH88cLT+xVO2k8SUYAp9NVBZRCqKSALQFZiVfQMRKQ/MAHoYY34KfUwVU1LWBiZcuvIBKFLGdhoVi4qWhSv7Bi5hTFlrO03MOGOhG2O8QF9gPrAFmGqM2SQifUSkT9ZmQ4HiwGsisl5E1oQtsYpuxsC8QVDoAmj8L9tpVCxr/BAUPB/m/1unBAiSK5iNjDFzgbmnPDc225/vAe4JbTQVkzZOh5RV0P5VSCxsO42KMsmD5uRq+27Odjx39G3uHTyM+f4Ged7v9hFt8/zeWKJ3iqrQ8RwLTMBV6lKoc6vtNCoOTPU15Ud/WQa5PsCN13acqKeFrkLn61fhcAq0fk5v8Vch4cPJs97uVHT8QQ/n57bjRD0tdBUah/fAV6OhWntIbmw7jYojX/prs9RXiwddMyhKmu04UU0LXYXGF0+B3wst9DIzFXrPeLtTmHT6uT62HSWqaaGrs7drHWyYHFgn9LyKttOoOPSjKc+Hvqbc4VzAxZJiO07U0kJXZ8cYmPcYFCwJV/e3nUbFsVHeW0gnkWGud/ibm9XzPS10dXY2zwysNNNsSGCxX6XC5ABFeMF7C42dm7je8Y3tOFFJC13lXWY6LBgKF9SEy3rYTqPygcm+5mzyV2CI+z0KkGE7TtQJ6sYilT+d6SaQAa4pPOD6jVv2Ps43/54XoVQqP/Pj4HHPncxIHEZf10ye93a1HSmq6BG6ypNKsotezjlM913NN6aa7TgqH1lnqjDN14R7nHO4SE6bJzBf00JXeWB42jWBdJJ41qN3hKrIG+HpRoaeID2NFrrKtY6O5Vzp3Mzz3q6kUtR2HJUP7acoo72daeL8nlaO1bbjRA0tdJUrRUhjsPs9vvVfzAe+a23HUfnYu74WbPGX43H3eyRx3HacqKCFrnJlgGsq53GEIZ67MPrtoyzy4eQJT0/Kyn7ud31iO05U0J9IFbTaspXuzi94x9eKTSbZdhylWGWq8bGvEfc6Z1NJdtmOY50WugqKAz9Pu8ezl2KM9na2HUepE57x3MYxEhnpfhPBbzuOVVroKig9nJ9Ty7Gd4Z4epFHAdhylTthPUZ7y9KCe4yduz+dT7GqhqzO6gAM84prKUl8t5vgb2o6j1Glm+K9mia82j7qmUFb22Y5jjRa6OgPDc+63cOPjce+dgNgOpFQOhH977sYgPOt6i/x6bboWuvpHXZxf0sy5nhHebuwwpWzHUepv7aYEI7zdaOL8ni7OL23HsUILXf2t0qTyuOtdVvqrMcnXwnYcpc7ofV9zvvFXZYjrPUpy0HaciNNCVzkzhpHucTjxM8DTW685VzHB4GCQpxeJeBjunkh+G3rRn1KVs7UTaeL8nue8t7LTXGA7jVJB+9WU5iVvZ1o7V9PGscp2nIjSQlenO7gDFgzhK18N3vc1t51GqVx7y3c93/kr8pR7AsU4YjtOxGihq5P5/TCrLwADdahFxSgfTh713EtRjjLS/WZgqcR8QH9a1cnWvA2/LoVWz7CLkrbTKJVnP5jyjPR2pZVzDaydYDtORGihq/85sA0+HwqVmkHdO2ynUeqsjfe1YamvFsz7N+z9wXacsNNCVwHe4/DRneB0Q/v/gugNRCr2GRw84rkPEgrC9LvBE9/rkGqhq4DPh8Ke9dDhNSha1nYapUJmH8Wg4+vwx0ZY+ITtOGGlha5g8yz4Ziw0vA+q3WA7jVKhV6UlNOwT+D7/ab7tNGGjhZ7fHdwOn/SFMnWhxVO20ygVPtc9CRfUhJn3w5HfbacJCy30/MybGRg3B+gyAVwJdvMoFU7uJOj0NmQehZn3BS7RjTNBFbqItBaRH0Vkq4gMyuH1qiKyQkSOi0j/0MdUYfH5UNi9DjqOgXOTbadRKvzOrwqtn4VfFsHyl2ynCbkzFrqIOIExQBugOtBNRKqfstkB4EFgVMgTqvDYMhu+eT0wrlitne00SkXO5XdCzU7wxXD4aYHtNCEVzBF6A2CrMWabMSYTmAJ0yL6BMWavMWY14AlDRhVqB3fAJ/dDmct03FzlPyLQ/lUoVStwKeP+n20nCplgCv1CYGe2xylZz+WaiPQWkTUismbfvvy7qohVx9NgSvfAJHSdJ4Ar0XYipSIvoQB0nQzOBPigKxw7ZDtRSART6DndYZKniRGMMeOMMfWMMfVKltTbyiPO7wsckezdDF3Gw3kVbSdSyp5i5eDmSYErvWb0Cvx8xLhgCj0FKJftcVlgd3jiqLBa8Dj8NA/ajISLr7OdRin7khtBm+fh5wWwaLjtNGctmEJfDVQWkYoikgB0BWaFN5YKudVvw8oxgZOgDXrZTqNU9Kh/d+BE6VcvwcbpttOcFdeZNjDGeEWkLzAfcALjjTGbRKRP1utjRaQUsAYoAvhF5CGgujHmcPiiq6Bt/QLmDoDKLaHVs7bTKBV92jwP+36AmQ/AeZWgTB3bifLkjIUOYIyZC8w95bmx2f78O4GhGBVt9v4AH/WEklWh83hwOG0nUir6uBIC4+njroXJt8Bd82LyHJPeKRrPju6HyTeDKwlu/RASC9tOpFT0KnQ+3DYNfMfh3Y4xOT1AUEfoyq7kQXNy/Z6CHOO9hOeoJru5JfNxNoz4Dvgu9OGUiifnV4Pu0+GddvDujdBzDhQ4z3aqoOkRehw6hwwmJDxPTfmVfp5+bDAX246kVOwoezl0mwypWwO/4WYetZ0oaFrocSaJ47ztHsXl8hP/53mAz/31bEdSKvZc1DRwzmnX2sCNeN7jthMFRQs9jiSSyTj3aK5wbOFfnvuZ67/CdiSlYle1doHVu7Ytjpkbj7TQ40QCHl53v0xjx0YGeO5llr+R7UhKxb7Lbgtc6rv5E/j0wagvdT0pGgfceBnj/g/NnOsZ5LmH6f4mtiMpFT+ufAAy/oQvR0JmOtz4RtSuHaCFHuMS8PCyewwtnGsZ4rmTKb5mtiMpFX+u/XdgoenPhwbK/ZZ3A4+jjA65xLAipDEpYQTXO1fxpKcH7/la2I6kVPxq9H+BaXe3LYZJHSH9gO1Ep9FCj1EXso/pCU9SV37iwcwHmOBrYzuSUvGvbg/o8g7sWQ8T28LhPbYTnUSHXGJQDfmVCQkvkEgmt3seY6X/1AWklFLZ5eXmvL/n5CpHf8b9MZoDo66mh+cxdphSOW65fUTbEO73zPQIPcY0dXzL1ISnyMRFp8xhWuZKWfC1vya3Zg6mkBxjWsKTXC4/2o4EaKHHlK7ORbzlfpFfTWluPP4kW43Oh6aULd+ZSnTJfII0k8SUhKfp6ZxHHtf+CRkt9BiQSCbPut5ihPstlvlrcXPmUPZxru1YSuV7v5gL6ZD5NEv8tRnmnsTL7jGcQ4a1PFro0W7vD3yS8Di3uhbxurcd93j6k06S7VRKqSyHKUhvz8O84LmZ9o4VzEh4ggpiZ6ZGLfRoZQysexfGNaWE/MntmQMZ6e2GD53PXKloY3AwxteRnp5HKSUH+TRhCM0dayOeQws9GmUcDswdMasvlKtPm+PPsdRf23YqpdQZLPXXpl3mM+ww5/N2wouBlcKOp0Vs/1ro0WbnKhh3TWBtw2uHQI+ZOl6uVAxJMSXpnDmMCd5WsOpNeO1K+GVRRPathR4t0g/ArH7wdgvwZsIds+GaAbpknFIx6DgJPOm9I7CUnSsxsFjGzAfg2MGw7jcmbywK7U0CwQvLTQJ+P2yYDAseD8wRcWVfaPoYJBYK/b6UUpFV/gro81VgYq/lr8DWhXDDaKganhuO9Ajdpj82w8Tr4ZMHoEQV6LMMWj2jZa5UPHEnwXVPQK9FULAkTLkVFj4Zll3F5BF6zDu8B756Cda8DYlFAhP+1OkODv33Vam4VaYO9F4My1+Gi64Nyy600CPpyO9ZRT4B/N7ARD/NhkLB4raTKaUiwemGJgPC9uW10CPhyO/w1cuwdgL4PFCnG1zdH86raDuZUiqOaKGH094fAsMq6yZpkSulwk4LPdQ8GbBlVmBY5bevwZkAtW6GJo/AeRfZTqeUimNa6KGyf2tgSGX9ZDh2IFDeLYYHTnbqGLlSKgK00HPh5OvfDZVlF20cq2jtXE11xw48xskC/+VM9jXn6901MLsd8OlKa3mVUvmLFnouCH5qynZaO1fR2rGaSo49+I2w1lRmuKc7s3yN2Ecx2zGVUvmUFvo/EPxcIilc4djMlY7NNHRsoZgcxWscrPBXZ7ynDQt8l+tcK0qpqKCFnk0R0qjh2EFN+ZW6jp9p6NjCeRKYKe03f0nm++qzwl+dxf46/InezamUii75stDdeCkne6koe7hEdlLTsZ2a8ivlHftObLPTX5JF/rqs8FVnpb8auyhpMbFSSp1ZnBa64TyOUFoOUFpSKS2pJMsfVJQ9VJQ9lJN9uMR/Yutf/RfwnanEZE9zNpqKbPJX4CBFLOZXSqncC6rQRaQ18ArgBN4yxow45XXJev16IB3oaYxZF+KsAekHaChbOE8OU1wOU5zDWX8+QnEOU0pSKS0HSBLPSW87ZhL41ZRmk0lmtv9KfvWXYpspwy+mDEcoEJaoSikVSWcsdBFxAmOAFkAKsFpEZhljNmfbrA1QOeu/hsDrWf8betsW82Hi8JOeOmQKkmqKkEoRNpqKLPDX43dzHrtNcfaY4uwx57GfohidXFIpFceCOUJvAGw1xmwDEJEpQAcge6F3ACYZYwywUkSKiUhpY8yekCdOvppbM//NAVOEVFOYgxTGG68jR0oplQvBNOGFwM5sj1M4/eg7p20uBE4qdBHpDfTOepgmIj/mKm3sKgHstx3CAv3c+Yt+7lPIyLDsr8LfvRBMoUsOz5k8bIMxZhwwLoh9xhURWWOMqWc7R6Tp585f9HPbF8ygcgpQLtvjssDuPGyjlFIqjIIp9NVAZRGpKCIJQFdg1inbzAJul4ArgD/DMn6ulFLqb51xyMUY4xWRvsB8ApctjjfGbBKRPlmvjwXmErhkcSuByxbvDF/kmJTvhpmy6OfOX/RzWyaBC1OUUkrFOr0wWyml4oQWulJKxQkt9AgRkRdE5AcR+U5EPhaRYrYzRYKIdBGRTSLiF5GouLQrnESktYj8KCJbRWSQ7TyRICLjRWSviGy0nSWSRKSciCwWkS1Z3+P/ZzuTFnrkfA7UNMZcCvwEPGY5T6RsBG4CltoOEm7ZpsloA1QHuolIdbupImIi0Np2CAu8wCPGmGrAFcADtv//1kKPEGPMAmOMN+vhSgLX6sc9Y8wWY0x+uSP4xDQZxphM4K9pMuKaMWYpcMB2jkgzxuz5axJCY8wRYAuBO+St0UK34y7gM9shVMj93RQYKs6JSDJwGfCNzRw6q1UIichCoFQOLw02xnyStc1gAr+qvR/JbOEUzOfOJ4KaAkPFFxEpBEwHHjLGHLaZRQs9hIwx1/3T6yJyB3AD0NzE0Q0AZ/rc+YhOgZHPiIibQJm/b4yZYTuPDrlESNYiIQOB9saYdNt5VFgEM02GihNZC/u8DWwxxoy2nQe00CPpVaAw8LmIrBeRsbYDRYKI3CgiKcCVwBwRmW87U7hknfT+a5qMLcBUY8wmu6nCT0Q+AFYAl4hIiojcbTtThDQCegDNsn6m14vI9TYD6a3/SikVJ/QIXSml4oQWulJKxQktdKWUihNa6EopFSe00JVSKk5ooSulVJzQQldKqTjx/1xjrn1XoMyJAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0eklEQVR4nO3deZyN9fvH8dd1zplB9mUsWRpE2cVY2kRFtFC08CVUX74q7RRpES0qpU1JJSoVUZKUSgsqZcguJSmDmFCWwcw55/r9cQ6/MYY5w5y5z3I9H495zLn3920517k/574/H1FVjDHGxB+X0wGMMcY4wwqAMcbEKSsAxhgTp6wAGGNMnLICYIwxccrjdID8qFChgiYnJzsdwxhjosrixYv/VtWknPOjqgAkJyeTmprqdAxjjIkqIvJHbvOtCcgYY+KUFQBjjIlTVgCMMSZORdV3AMaY6JGVlUVaWhr79+93OkrcKFq0KNWqVSMhISGk9a0AGGPCIi0tjZIlS5KcnIyIOB0n5qkq27dvJy0tjZo1a4a0jTUBGWPCYv/+/ZQvX97e/AuJiFC+fPl8XXGFVABEpKOIrBWRdSIyJJflPUVkefDnOxFpkte2IlJORD4XkV+Dv8uGnNoYExXszb9w5ffPO88CICJuYCzQCagP9BCR+jlW+x04T1UbAyOB8SFsOwSYq6p1gLnBaWOMMYUklCuAlsA6VV2vqpnAu0CX7Cuo6nequjM4uRCoFsK2XYBJwdeTgMuP+yyMMSYCPPDAA3zxxRcFsq/BgwfToEEDBg8ezLhx43jjjTcAmDhxIps3by6QY4TyJXBVYGO26TSg1THWvwH4JIRtK6nqFgBV3SIiFXPbmYj0B/oD1KhRI4S4xhS+5CEfO3LcDaMuceS4JncjRowosH29/PLLpKenU6RIkcPmT5w4kYYNG3LyySef8DFCuQLIrVEp12HERKQdgQJwT363PRpVHa+qKaqakpR0RFcWxhiTq5EjR3L66afTvn17evTowejRowF45ZVXaNGiBU2aNKFbt25kZGQA0LdvX6ZNm3Zo+xIlSgCwZcsW2rRpQ9OmTWnYsCHz58/H5/PRt29fGjZsSKNGjRgzZswR+xgxYgQtWrSgYcOG9O/fn4OjL7Zt25Z77rmHli1bUrduXebPn39E9s6dO7N3715atWrFlClTGD58OKNHj2batGmkpqbSs2dPmjZtyr59+07ozyiUK4A0oHq26WrAEdcfItIYeBXopKrbQ9h2q4hUCX76rwJsy294Y0yU+GQI/LWiYPdZuRF0GpXrotTUVKZPn85PP/2E1+ulWbNmNG/eHICuXbvSr18/AO677z5ee+01brnllqMe5u233+aiiy5i2LBh+Hw+MjIyWLp0KZs2bWLlypUA/PPPP0dsN3DgQB544AEArr32WmbNmsVll10GgNfr5ccff2T27Nk89NBDRzQbzZw5kxIlSrB06VIAhg8fDsCVV17JCy+8wOjRo0lJSQntz+kYQrkCWATUEZGaIpIIdAdmZl9BRGoA7wPXquovIW47E+gTfN0H+PD4T8MYY/7fggUL6NKlC8WKFaNkyZKH3ngBVq5cybnnnkujRo2YPHkyq1atOua+WrRoweuvv87w4cNZsWIFJUuWpFatWqxfv55bbrmFTz/9lFKlSh2x3VdffUWrVq1o1KgRX3755WHH6dq1KwDNmzdnw4YNBXPSxyHPKwBV9YrIQGAO4AYmqOoqERkQXD4OeAAoD7wYvA3JG2y2yXXb4K5HAVNF5AbgT+CqAj43Y0ykOMon9XA52NySm759+zJjxgyaNGnCxIkT+frrrwHweDz4/f5D22dmZgLQpk0b5s2bx8cff8y1117L4MGD6d27N8uWLWPOnDmMHTuWqVOnMmHChEPH2L9/PzfddBOpqalUr16d4cOHH3Z//sF2fbfbjdfrLejTD1lIzwGo6mxVrauqtVX1keC8ccE3f1T1v6paVlWbBn9SjrVtcP52Vb1AVesEf+8o6JMzxsSnc845h48++oj9+/ezZ88ePv74/7+k3717N1WqVCErK4vJkycfmp+cnMzixYsB+PDDD8nKygLgjz/+oGLFivTr148bbriBJUuW8Pfff+P3++nWrRsjR45kyZIlhx3/4Jt9hQoV2LNnz2HfLZyokiVLsnv37gLZl3UFYYyJOS1atKBz5840adKEU045hZSUFEqXLg0Evhxu1aoVp5xyCo0aNTr0ZtqvXz+6dOlCy5YtueCCCyhevDgAX3/9NU8++SQJCQmUKFGCN954g02bNnHdddcdumJ47LHHDjt+mTJl6NevH40aNSI5OZkWLVoU2Ln17duXAQMGUKxYMb7//nuKFSt23PuSY10qRZqUlBS1AWFMJLLbQI+0Zs0a6tWr59jx9+zZQ4kSJcjIyKBNmzaMHz+eZs2aOZansOT25y4ii7O3zBxkVwDGmJjUv39/Vq9ezf79++nTp09cvPnnlxUAY0xMevvtt52OEPGsABgTITx4OUPW0di1nmqSTlX5m+LsJxMPB0hkm5ZhjdbgZ38NVmkymSQ41vQEkd38ZEJjBcAYB3nw0sGVyhXubznTtYoSErh7ZI8WJU2T2MVJFGcfRcniLNdKekvgyc9dWowv/M352NeKef4mZNl/ZXMc7F+NMQ4oxV76uOfQ0zOXyrKTNK3ADN/ZzPc35kf/aeykJEf2pKJUk3QayAbOdy3lIvciuroXsEnLM857GVN9bTlAohOnY6KUFQBjCpEbHz3cX3Kn5z3KyR6+8TVmqO+/fONvgj/Px3KENK1ImlZkjr8l93mv5zzXMv7n+YiRCRO51fMBT3qvZqqvLbl3w2XM4awAGFNIGstvPJEwntNdG/neV5+R3l6s1uTj3l8WHr7wN+eLzGa0kp+5K2EqTyS8Qjf3fO7NuoHftGrBhS8ABf19RSR9B9G2bdt89c8zePBgPvroIxITE6lduzavv/46ZcqUCW/IXNiQkMaEmQs/N7k/ZHricEpKBv/LvIMeWcNO6M3/cMIPWo9rMu/n7qx+nCYbmZ04lL7uT8ln57umkLRv356VK1eyfPly6tate8SDZIXFCoAxYVSef5mc8Ch3J0xhjr8FnQ48xhx/C8LRRKO4mOprxwUHRjPP34ThCW/wQsLzFOfEugyORnv37uWSSy6hSZMmNGzYkClTpgDH7qL5jjvuoE2bNtSrV49FixbRtWtX6tSpw3333QfAhg0bOP300+nTpw+NGzfmyiuvPNSVdHafffYZZ555Js2aNeOqq65iz549R6zToUMHPJ5AA0zr1q1JS0sL1x/FMVkBMCZM6kgaMxIfoKlrHYOy/sfArFvYRYmwH3c7pemfdQejsrrTyfUDMxPvo4ZsDftxI8mnn37KySefzLJly1i5ciUdO3YEAl00L1q0iJUrV7Jv3z5mzZp1aJvExETmzZvHgAED6NKlC2PHjmXlypVMnDiR7dsDPdyvXbuW/v37s3z5ckqVKsWLL7542HH//vtvHn74Yb744guWLFlCSkoKTz/99DGzTpgwgU6dOhXwn0BorAAYEwbnuFYwPfFBikgW12TezzTfeRTmF7OKi3G+zvTMGkZZ2c30xAdpIL8X2vGd1qhRI7744gvuuece5s+ff6gfoGN10dy5c+dD2zZo0IAqVapQpEgRatWqxcaNgYENq1evztlnnw1Ar169WLBgwWHHXbhwIatXr+bss8+madOmTJo0iT/++OOoOR955BE8Hg89e/Ys0PMPlRUAYwpYR9ePvJ7wBJu0ApcfGMFyre1YloX++lyV+SAHSGRK4kjOcq10LEthqlu3LosXL6ZRo0YMHTqUESNGHOqiedq0aaxYsYJ+/frl2kWzy+U6bBhGl8t1qMvmYHf3h+ScVlXat2/P0qVLWbp0KatXr+a1117LNeOkSZOYNWsWkydPPmI/hcUKgDEF6BLXQl5IeI7lWourMx9kMxWcjsRvWpVuB4aTpkm8nvAEbV1LnY4Udps3b+akk06iV69eDBo0iCVLlhRIF81//vkn33//PQDvvPMO55xzzmHLW7duzbfffsu6desAyMjI4JdffjliP59++imPP/44M2fO5KSTTsp3joJit4EaU0A6u75jTMJYFmtdrsu8m70cfze9BW0r5bg6837eTnyUcQljuC7rbr73NyjUDIV52+aKFSsYPHgwLpeLhIQEXnrppQLporlevXpMmjSJ//3vf9SpU4cbb7zxsOVJSUlMnDiRHj16cODAAQAefvhh6tate9h6AwcO5MCBA7Rv3x4IFI5x48Yd59keP+sO2pgCcP29Ixmf8DSpehrXZw4mg6JOR8pVWXbxTuIj1JBt9M68h1Q9/bj3ldcbutPdQRe0DRs2cOmllx4aBzhS5ac76JCagESko4isFZF1IjIkl+Wni8j3InJARAZlm3+aiCzN9rNLRG4PLhsuIpuyLbs4vydqTETYuIixCc+xWk/hhsxBEfvmD7CTUvTKvJctWo7XE5/kNPnT6UjGQXkWABFxA2OBTkB9oIeI1M+x2g7gVmB09pmquvbgMJFAcyAD+CDbKmOyDSM5+/hPwxiHpP8Cb1/FVi0bcc0+R/M3pemVeS97KcqExCdJYqfTkaJCcnJyxH/6z69QrgBaAutUdb2qZgLvAl2yr6Cq21R1EZB1jP1cAPymqke/J8qYaLJ3O0zuBi4PvbOGsJ3STicK2RbKc0PmYMqwh9cSR1OM/XlvdByiqYk5FuT3zzuUAlAV2JhtOi04L7+6A+/kmDdQRJaLyAQRKZvbRiLSX0RSRSQ1PT39OA5rTBj4suC9PrB7K/SYwp9ayelE+bZKkxmYdSsNZAPPJoxF8Bfo/osWLcr27dutCBQSVWX79u0ULRp6E2QodwHldoNqvv5GRSQR6AwMzTb7JWBkcF8jgaeA6484kOp4YDwEvgTOz3GNCZs5w2DDfLh8HFRrDjg3MMuJ+Mp/BiO8vXkoYRID/TN43te1wPZdrVo10tLSsA9uhado0aJUq1Yt5PVDKQBpQPVs09WAzfnM1QlYoqqHnkfP/lpEXgFm5bahMRFnyZvw48tw5kBo2sPpNCdskq8DjV2/cYdnOiu0Fl/7mxbIfhMSEqhZs2aB7MuERyhNQIuAOiJSM/hJvjswM5/H6UGO5h8RqZJt8gogtr5dMbHprxXw8V1Qqy1c+JDTaQqIMCzrBn7WGjyTMJZqss3pQKaQ5FkAVNULDATmAGuAqaq6SkQGiMgAABGpLCJpwJ3AfSKSJiKlgstOAtoD7+fY9RMiskJElgPtgDsK7KyMCYcDe+C9vlCsLHR9Fdyx8xzlfoowIOt2BOWlhGdIPOb9HCZWhPQvOHiL5uwc88Zle/0Xgaah3LbNAMrnMv/afCU1xkmq8PGdsGM99J4JJZKcTlTg/tRKDMoawCuJTzPYM4VHvL2cjmTCzPoCMiYUSyfD8ilw3hCoea7TacLmc38Kb3jb088zmzauZU7HMWFmBcCYvOxYD7PvhuRzoc2gvNePco94e/KzvzpPJYyjPP86HceEkRUAY47F74MPbgSXB64YBy6304nC7gCJ3Jo1kJJk8ETCeGxYydhlBcCYY/n2Wdi4EC5+EkqHfn91tPtFqzPK24ML3D9xlfsbp+OYMLECYMzR/LUCvnoU6neBxlc7nabQTfJ1YKG/Hvd73qQK252OY8LACoAxufFlwYwbA7d8XjIGHBqxyUmKi8FZ/XHj53FrCopJVgCMyc13zwWuAC59GoofcRdz3NiolXjM+x/auFfwH/eXTscxBcwKgDE5pf8CXz8eaPqpd5nTaRz3lu9CvvU1YIjnbSqxw+k4pgBZATAmO78fZt4CCcWg05NOp4kQwr3eG0jEy0MJk5wOYwqQFQBjskt9LXDXT8dRUDL6ungOlz+0Ms94u9HRvYiLXIucjmMKiBUAYw7avRXmjgh09Naku9NpIs6rvotZ7T+FhxImUpIMp+OYAmAFwJiD5twL3v1w8VNxeddPXrx4GJL1X5L4h7s8U52OYwqAFQBjAH77ClZOg3PuhAqnOp0mYi3X2rzpa8+17s9h81Kn45gTZAXAmKz9gT7+y9WCc6xX8rw87b2KHZSE2YMCX5qbqGUFwJjvnocdv8ElT0FC6OOpxqtdFOexrP9A2iJY+pbTccwJsAJg4ts/G2H+U1CvM9Q+3+k0UeN9/7lQ4yz4/EHIsGcDopUVABPfPrsv8PuiR5zNEXUELhkN+/+FLx92Oow5TiGNCCYiHYFnATfwqqqOyrH8dOB1oBkwTFVHZ1u2AdgN+ACvqqYE55cDpgDJwAbgalXdeWKnY+Jd8pCPQ173TNcq3kmcwdNZV/LcqBXAivAFi0WVGkCL/8KiVyDleqjc0OlEJp/yvAIQETcwFugE1Ad6iEj9HKvtAG4FRpO7dqra9OCbf9AQYK6q1gHmBqeNKRRufDzoeYON/iRe9l3qdJzo1XYIFC0Nnw4JDJtpokooTUAtgXWqul5VM4F3gS7ZV1DVbaq6CPI1knQX4OBz5ZOAy/OxrTEn5D/uuZzu2sjD3l4cINHpONHrpHLQbhhsmA9rZjqdxuRTKAWgKrAx23RacF6oFPhMRBaLSP9s8yup6haA4O+KuW0sIv1FJFVEUtPT0/NxWGNyV4o93OGZxne++szxp+S9gTm25tdBxQaB71Oy9jmdxuRDKAUgt0ci83Otd7aqNiPQhHSziLTJx7ao6nhVTVHVlKSkpPxsakyubvV8QBn28rC3F7n/8zb54vZAp1Hwz5+w8EWn05h8CKUApAHVs01XAzaHegBV3Rz8vQ34gECTEsBWEakCEPy9LdR9GnO8kmULvd2fMdV3Hqs12ek4saNmGzjtYpg/BvbYlXq0CKUALALqiEhNEUkEugMhNfaJSHERKXnwNdABWBlcPBPoE3zdB/gwP8GNOR73et4mkwSe8sbfEI9h134EePfB1486ncSEKM8CoKpeYCAwB1gDTFXVVSIyQEQGAIhIZRFJA+4E7hORNBEpBVQCFojIMuBH4GNV/TS461FAexH5FWgfnDYmbFrJGjq4F/OitzPplHE6TuypUAdSboDFE2HbGqfTmBCE9ByAqs4GZueYNy7b678INA3ltAtocpR9bgcuCDmpMSdA8HNvwmQ2azle813sdJzYdd49sOxd+Ox+6DXN6TQmD/YksIkLl7kW0sS1nqeyrrbbPsOpeHloMwjWfR7oYdVENCsAJuYlksVgzxRW+0/hA/85TseJfa3+B6VrwOcPWG+hEc4KgIl517o/o7ornUe9/8Fv/+TDz1MEzr8P/loOK6c7ncYcg/1vMDGtFHu4xTODb3yNWeBv5HSc+NHoKqjcCL4cAd4DTqcxR2EFwMS0Gz0fUYoMRnl7OB0lvrhcgdtC//kTFr3qdBpzFFYATMyqxA6uc3/Kh/6zWKOnOB0n/tQ+H2q1g3lPwr5/nE5jcmEFwMSs2z3TceHnKe9VTkeJX+0fgn07A6OumYhjBcDEpNqyiavdXzPZdyFpmms/g6YwVGkCDbsF+gjavdXpNCYHKwAmJg3yTGUfRXjBe7nTUUy7YeDLDDQFmYhiBcDEnCayjk7uRbzivYTtlHY6jilfG5r1gcWvw471Tqcx2VgBMDFnkGcq27Ukr1qXD5HjvLvBlQBfWUdxkSSkvoCMiRZnulZxrnslI7N6sZdiTscxB5WsDK1vhAVPw9m3BZ4ROE75Gfe5IG0YdYkjxw0nuwIwsUOVezzvsknL85bvQqfTmJzOvjUwfvCXjzidxARZATCxY+1smrp+41lvV+vwLRIVKxv49P/LJ7DxR6fTGKwAmFjh98GXD/ObvwrTffkaddQUplYDoHgSzB0Bmp+RZU04WAEwsWHldNi2mjHeK/HhdjqNOZrE4tBmMGyYD+u/djpN3LMCYKKfLytwd0mlRnzsb+V0GpOX5n2hdHW7CogAIRUAEekoImtFZJ2IDMll+eki8r2IHBCRQdnmVxeRr0RkjYisEpHbsi0bLiKbRGRp8Mfu2TPHZ+lk2Pk7nD8Mtc80kc9TJDBy2OYlsHZ23uubsMnzf4uIuIGxQCegPtBDROrnWG0HcCswOsd8L3CXqtYDWgM359h2jKo2Df7YvwSTf1n74ZsnoGoK1O3odBoTqiY9oFztwB1BNmiMY0L5uNQSWKeq61U1E3gX6JJ9BVXdpqqLgKwc87eo6pLg690EBpWvWiDJjYHA06W7NsEF94OI02lMqNweaHcvbFsFqz9wOk3cCqUAVAU2ZptO4zjexEUkGTgD+CHb7IEislxEJohI2aNs119EUkUkNT09Pb+HNbEscy/MfwqSz4VabZ1OY/KrQVeoWD/w/Y3P63SauBRKAcjtY1W+vrkRkRLAdOB2Vd0VnP0SUBtoCmwBnsptW1Udr6opqpqSlJSUn8OaWPfjeNibDuff73QSczxcrkBHcdvXwfIpTqeJS6EUgDSgerbpasDmUA8gIgkE3vwnq+r7B+er6lZV9amqH3iFQFOTMaHZ/y8seAbqdIAadudP1Dr9Ejj5DPhmFHgznU4Td0IpAIuAOiJSU0QSge7AzFB2LiICvAasUdWncyyrkm3yCmBlaJGNARa+BPv/CbQjm+glAu3uCwwd+dObTqeJO3l2BqeqXhEZCMwB3MAEVV0lIgOCy8eJSGUgFSgF+EXkdgJ3DDUGrgVWiMjS4C7vDd7x84SINCXQnLQB+F8BnpeJZRk74PuxUO+ywKdHE91OvQCqt4Z5o6FpT0go6nSiuBFSb6DBN+zZOeaNy/b6LwJNQzktIPfvEFDVa0OPaUw23z0PB3ZDW/v0HxNE4PxhMOmywF1drW90OlHcsKdmTHTZkw4/jAsMM1gp5+MoJmrVbBP4mf9U4O4uUyisAJjo8u0z4N0PbYc6ncQUtHb3Be7q+nG800nihhUAEz12bYFFrwaeIq1wqtNpTEGr0QpObQ/fPgv7d+W9vjlhVgBM9Jj/FPi9geEFTWw6fxjs2xm4y8uEnQ0JaQpcOIbsq0o6XxWZwHu+tgx7fBWwqsCPYSLAyWfA6ZcG7vJq1T8wiIwJG7sCMFFhoGcGivCC93Kno5hwa3cvHNgF373gdJKYZwXARLxT5C+ucn/D274L2EJ5p+OYcKvUABpcEWgG2vu302limhUAE/Fu9bxPFh5e9HbJe2UTG9oOBe8+WDDG6SQxzQqAiWi1ZRNXuL5lkq8D6ZRxOo4pLEl1ofE1gbu+dv/ldJqYZQXARLQ7PNPJoAgvey91OoopbOfdHbjra36uHQWbAmAFwESsevIHl7oXMsHXkZ2UcjqOKWzlasEZvSD19UBncabAWQEwEetOzzR26Um86rXhouNWm8GBvoLmPel0kphkBcBEpCayjvbuxbzsvZRdlHA6jnFK6WrQ/Dr4aTJs/83pNDHHCoCJSHd53mO7lmSi7yKnoxinnXsnuBPhm8edThJzrACYiNNK1tDGvYIXvZ3ZSzGn4xinlawMLfvB8qmwbY3TaWKKFQATYZS7EqayVcvwlq+902FMpDjnDkgsAV894nSSmBJSARCRjiKyVkTWiciQXJafLiLfi8gBERkUyrYiUk5EPheRX4O/rdMPQxvXclq61vK89woOkOh0HBMpTioHZ94Maz6ioax3Ok3MyLMAiIgbGAt0IjDMYw8RyTkSxw7gVmB0PrYdAsxV1TrA3OC0iWvKIM9UNvqTmOJr53QYE2nOvAmKluEuz3tOJ4kZoVwBtATWqep6Vc0E3gUOeyZfVbep6iIgKx/bdgEmBV9PAi4/vlMwseIi1yIau37nGW83sqyjWpNT0dJwzu20cy8jRX52Ok1MCKUAVAU2ZptOC84LxbG2raSqWwCCvyuGuE8Tg1z4GeR5j3X+k/nAf47TcUykatmfbVqGwQlTAXU6TdQLpQDkNqh7qH/yJ7JtYAci/UUkVURS09PT87OpiSJdXN9Sx7WJp7xX4bd7E8zRJBbnee/ltHL9TBvXcqfTRL1Q/qelAdWzTVcDNoe4/2Ntu1VEqgAEf2/LbQeqOl5VU1Q1JSkpKcTDmmiSgJc7PNNY4U/mU38Lp+OYCPeu73w2+pMY7JmC4Hc6TlQLpQAsAuqISE0RSQS6AzND3P+xtp0J9Am+7gN8GHpsE0uucX9FDVc6T3mvRu3Tv8lDFh6e8XajkWsDHV2LnI4T1fL836aqXmAgMAdYA0xV1VUiMkBEBgCISGURSQPuBO4TkTQRKXW0bYO7HgW0F5FfgfbBaRNnirGfWz0f8IP/dL72N3E6jokSH/jP4Vd/Ve7yvIcbn9NxolZIt1qo6mxgdo5547K9/otA805I2wbnbwcuyE9YE3v6uj+jovzDgMzbyf0rI2OO5MfFaO9VvJz4DF3d83nP19bpSFHJrreNY0qxhwGemXzhO4MlWtfpOCbKzPG3YKm/Frd7plOETKfjRCUrAMYxAzyzKMk+RnuvcTqKiUrC494eVJXt9HJ/7nSYqGQFwDiiIju5zv0pH/rP4met4XQcE6W+9zdgnq8RN3s+pCQZTseJOlYAjCNu90zHjY8x3iudjmKi3OPe7pSTPfzX87HTUaKOFQBT6GrJZq52f81k34X8qZWcjmOi3CqtyUe+1vzXPZsK/Ot0nKhiBcAUukGeqewnkRe8lzsdxcSIp7xXkYiXWz3vOx0lqlgBMIXqDPmVi90/Mt57Kdsp7XQcEyM2aBXe8Z1PD/eXJMsWp+NEDSsAphAp9yS8S7qW4lWfDfRuCtZz3q4cIIHBnilOR4kaVgBMoTnf9ROtXWt4ztuVDIo6HcfEmL8pzSveS7jE/SNNZZ3TcaKCFQBTKNz4GOp5h9/8gUt1Y8LhFd8lpGtphia8jXUXnTcrAKZQXO3+mjquTTzh7Y7XBnsxYZJBUZ7xdqOV62cudC1xOk7EswJgwu4k9nOnZxqL/HWZ409xOo6Jce/62rHOfzJDPW/jwet0nIhmBcCEXX/PLJLkXx7N6ol1+GbCzYebR73/obZrC93dXzkdJ6JZATBhVZGd9Hd/zCxfK37SOk7HMXHiS/8ZfOerzx2eadZFxDFYATBhNcgzFTc+Hvd2dzqKiSvCI96elJfd3OgJdfyq+GMFwIRNA/mdK93zeN3XkY3W5YMpZKu0JtN953CD+xOqiY0nnhsrACZMlPs8k9lJCV70dnE6jIlTo7OuwY9wj+cdp6NEpJAKgIh0FJG1IrJORIbkslxE5Lng8uUi0iw4/zQRWZrtZ5eI3B5cNlxENmVbZo+GxpAOrlTOdK9mjPdKdlHc6TgmTm2hPOO8l3GZeyHNZa3TcSJOngVARNzAWKATUB/oISL1c6zWCagT/OkPvASgqmtVtamqNgWaAxnAB9m2G3NweXDoSBMDEsliqOdtfvVXtYe+jONe9l3KZi3HAwlvIvidjhNRQnkipyWwTlXXA4jIu0AXYHW2dboAb6iqAgtFpIyIVFHV7L0yXQD8pqp/FFB2E6Guc39KTddWemfegw+303FMmCQPiY7+9/dThMezuvNs4otc4VrA+/42TkeKGKE0AVUFNmabTgvOy+863YGcDXEDg01GE0SkbG4HF5H+IpIqIqnp6fZFTsTb/Re3eD7gc18z5vmbOJ3GGABm+s/iJ/+p3JPwLsXZ53SciBFKAcjtyZ2cnWwccx0RSQQ6A+9lW/4SUBtoCmwBnsrt4Ko6XlVTVDUlKSkphLjGUXNHkEgWD3t7OZ3EmEMUFw9l9aaS/MNAzwyn40SMUApAGlA923Q1YHM+1+kELFHVrQdnqOpWVfWpqh94hUBTk4lmaYth6WQm+C7mD63sdBpjDrNUT+U9bxtucM+mpo0ZAIRWABYBdUSkZvCTfHcg55MVM4HewbuBWgP/5mj/70GO5h8RqZJt8gpgZb7Tm8jh98Mnd0PxijxvI32ZCPW4twf7SeR+z5tOR4kIeRYAVfUCA4E5wBpgqqquEpEBIjIguNpsYD2wjsCn+ZsObi8iJwHtgZxjtT0hIitEZDnQDrjjRE/GOGjpW7ApFTqMZC/FnE5jTK7+pjTPertyvnsp51tvoUjgxp3okJKSoqmpqU7HMDll7IDnm0PSaXDdJyQPtTt6TeRKwMsniUNIwEuHzCc4QGJI220YdUmYk4WPiCxW1SO64rUngc2J+3Ik7P8XLh4NYr19msiWhYcHvH05xbUt7vsJsgJgTszmnyD1dWjZHyo3dDqNMSH5zt+Qmb4zudH9EafIX07HcYwVAHP8/D6YdScUT4J2Q51OY0y+PJzVi0w8POSZRLwOH2kFwBy/Ra/B5iXQ8TEoWtrpNMbkyzbK8rT3Stq6l9HRtcjpOI6wAmCOz64tMHcE1GoHDbs5ncaY4/KGrwOr/afwYMIblIjDgWOsAJjjM2co+DLhkqfsi18TtXy4GZp1A5XYySDPVKfjFDorACb/fv0cVn0AbQZD+dpOpzHmhCzTU5nk60Bv9+ecIb86HadQWQEw+XNgN8y6AyqcBmff6nQaYwrEaO/V/EVZHk14FQ9ep+MUGisAJn/mjoR/06Dz8+Ap4nQaYwrEXorxYFZf6rk20t8dHd1cFwQrACZ0G3+EH8cH7vmv0crpNMYUqM/9Kcz2teQ2z/vUlk1OxykUVgBMaLwH4MOBULoaXPCA02mMCYsHs/qSQRGeSBiPKw5GD7MCYELzzRPw91q49BkoUsLpNMaERTpleCirN81dv9LXPcfpOGFnBcDkbdNiWDAGmvaEOhc6ncaYsJrhP5u5vjMY7JlCDdma9wZRzAqAObas/fDBjVCyMlz0qNNpjCkEwrCs68nCzZMJL8d0U5AVAHNsXz8aaPrp/BwUK+N0GmMKxV+UZ3hWH1q5fuYGd+x2b24FwBzdnz/Ad89D875wqjX9mPjyvv9c5vhSGOSZSl3Z6HScsLACYHJ3YDe83w9KV4f2I51OY4wDhHuzbmA3J/F0wkvgzXQ6UIELqQCISEcRWSsi60RkSC7LRUSeCy5fLiLNsi3bEBz6camIpGabX05EPheRX4O/yxbMKZkC8ck98O9G6DoeipZyOo0xjthOae7NuoGGrg3w9WNOxylweRYAEXEDY4FOQH2gh4jUz7FaJ6BO8Kc/8FKO5e1UtWmOIcmGAHNVtQ4wNzhtIsGqGbB0Mpx7F9Ro7XQaYxz1mb8F73rbBu6E+32+03EKVChXAC2Bdaq6XlUzgXeBLjnW6QK8oQELgTIiUiWP/XYBJgVfTwIuDz22CZt/N8FHt8HJzeC8e5xOY0xEGOHtHej48P3+gTGwY0QoBaAqkP0bkLTgvFDXUeAzEVksIv2zrVNJVbcABH9XzO3gItJfRFJFJDU9PT2EuOa4+bww/b/gy4Kur4A7welExkSEDIpCt1dhbzp8dCtobIwgFkoByK2z95xnf6x1zlbVZgSaiW4WkTb5yIeqjlfVFFVNSUpKys+mJr++eRz+/A4uHQMVTnU6jTGR5eQzAt2grPkIUl9zOk2BCKUApAHVs01XAzaHuo6qHvy9DfiAQJMSwNaDzUTB39vyG94UoPXfwLwnA0/7NrnG6TTGRKYzB8Kp7eHTobB5qdNpTlgoBWARUEdEaopIItAdmJljnZlA7+DdQK2Bf1V1i4gUF5GSACJSHOgArMy2TZ/g6z7Ahyd4LuZ47d4auOWzQh24+Emn0xgTuVwuuOJlKJ4E7/WB/f86neiE5FkAVNULDATmAGuAqaq6SkQGiMiA4GqzgfXAOuAV4Kbg/ErAAhFZBvwIfKyqnwaXjQLai8ivQPvgtClsviyYdh3s3wVXvg6JxZ1OZExkK14erpwA/2wM9JAbxd8HeEJZSVVnE3iTzz5vXLbXCtycy3brgSZH2ed24IL8hDVh8MVw+ONbuGI8VG7odBpjokON1nDhg/D5A/D9WDhroNOJjos9CRzPVs2A71+AFv+1dn9j8uusW6FeZ/j8/sB3aFHICkC82roaPrwZqrWAi2LvCUdjwk4ELn8RKtSF9/rCP386nSjfrADEo4wd8E73QHv/1W+AJ9HpRMZEpyIlofvb4PfBlF6QmeF0onyxAhBvfFkwtTfs/ivwD7fUyU4nMia6la8N3V6BLcthxgDwR8/4AVYA4olqoJO3DfMD/ftXS8l7G2NM3upeBB1GwuoPA2NoRImQ7gIyMeL7sYEnGM+6FZp0dzqNMbHlzIHw9y+BByrL14mKGyusAMSL1R/CZ/dB/S5w4UNOpzEm9ojAxU/Bjt9h5sBA82rNc51OdUzWBBQPNv4Y6MWwWovAU4wu+2s3Jiw8iXDNm1CuFrzbE7aucjrRMdk7Qazb9jO8fXXg00iPdyGhmNOJjIltxcpCz2mQeBK8dSX8m+Z0oqOyAhDL/vkT3rwC3InQ6/3AI+zGmPArUz1QBDL3wFvdYO92pxPlygpArNqzDd64HLL2wrUfQLmaTicyJr5Ubhi41XrnBnjriojsOM4KQCzK2AFvdoVdm+E/70GlBk4nMiY+1TwXrn4z8OT95Kshc6/TiQ5jBSDWZOyAN7oEbkfr8TbUaOV0ImPiW90OgdHE0n6Ed3pE1NPCVgBiyb6dgTb/9J8Dl561z3c6kTEGoMHlcPlL8Pu8wE0ZEXIlYAUgVuzdHvjkv3UVXDMZ6lzodCJjTHZNugduw/7j20Bz0IE9TieyAhATdm2BiRdD+trAJ/+6HZxOZIzJTZNrAmNv/PkdvNU1cNXuINEoGs0mJSVFU1NTnY4RWXb+Efjkvzc9cJ9/ticPk4d87GAwY8zRdHL9wLMJL/CbVqV35j2kU/aY628YdckJHU9EFqvqEZ1/hXQFICIdRWStiKwTkSG5LBcReS64fLmINAvOry4iX4nIGhFZJSK3ZdtmuIhsEpGlwZ+LT+QE49KWZfBae9i3A3p/GPGPnRtjAj7xt+L6rLupIVuZlvgQNWSrIznyLAAi4gbGAp2A+kAPEamfY7VOQJ3gT3/gpeB8L3CXqtYDWgM359h2jKo2Df4cNuSkycO6ufD6xeBKgOs/s549jYkyC/yN6Jk5jFKSwfuJD3KG/FroGUK5AmgJrFPV9aqaCbwLdMmxThfgDQ1YCJQRkSqqukVVlwCo6m4Cg8pXLcD88WnJG4E7Ccomw3+/gIqnO53IGHMcluqpdMsczl4tyjuJD3Oxa2GhHj+UAlAV2JhtOo0j38TzXEdEkoEzgB+yzR4YbDKaICK5NoKJSH8RSRWR1PT09BDixjCfFz4ZAjNvgZpt4LpPoFQVp1MZY07Aej2ZKzJHsEqTeTHxOW52zwAK57vZUAqA5DIvZ7pjriMiJYDpwO2quis4+yWgNtAU2AI8ldvBVXW8qqaoakpSUlIIcWNUxg6YfCX88BK0vinwhG/RUk6nMsYUgB2U4j+Zw5jhO4vBCVN5KeEZirMv7McNpQCkAdWzTVcDNoe6jogkEHjzn6yq7x9cQVW3qqpPVf3AKwSamkxuNi2Gl8+DDQug8wvQ8TFw21AOxsSSAyRye9bNjMzqSQdXKh8kPkBN2RLWY4ZSABYBdUSkpogkAt2BmTnWmQn0Dt4N1Br4V1W3iIgArwFrVPXp7BuISPa2iyuAlcd9FrFKFX58BSZ0BBSunwPNrnU6lTEmbITXfJdwbdZQKsi/zEq8l8tdC8J2tDwLgKp6gYHAHAJf4k5V1VUiMkBEBgRXmw2sB9YR+DR/U3D+2cC1wPm53O75hIisEJHlQDvgjgI7q1iwdztM6QWzB0GttvC/eVCtudOpjDGF4Dt/Qy4+8BgrtSbPJL4IM24OS/cR9iBYJFo3F2bcFLi//4IHoPXNxzWKlz0IZkx0c+PjVs/73OaZAde8BfUuPa79HO1BsLhpSHbyzTDkp/j2/wuf3Q9LJkGF06Dne1ClcXjDGWMilg83Y7xXcdttQyHptALff9wUgIj3yxz46HbY8xecdQu0G2bDNxpjAsLw5g9WAJy38w+Ycy/8PAsq1g9c5llbvzGmEFgBcEpmBnz/Asx/GkTgwuGBtn5PotPJjDFxwgpAYfP7YNk78OUjsHsz1OscuK+/dDWnkxlj4owVgMLi98OamfDN47BtNVRtDle+Bqec5XQyY0ycsgIQZi78sGoGzHsStq6E8nXgygnQoGug6ccYYxxiBSBMipBJN/d8+rlnwXtboVztwEhAja4El9vpeMYYYwWgoFWTbfRyf8E17q8pK3tY5q8FV02CepfZG78xJqJYASgAiWRxoWsxV7u/oY1rOX6Ez/3NmeS7iIX+emxocHxP7xljTDhZAThOLvy0kLV0dn/Hxe4fKCt72KzleMHXhbe9F/AX5Q+ta10yGGMikRWAfEjASyvXGjq4UungTqWy7CRDi/C5vznTfG341t8Qf2jDLBtjjOOsAOShCttp417Oea5lnONaSSnJYJ8m8o2/CbN8rZnrP4N9FHU6pjHG5JsVgMMo1eRvmstaWrvW0Nq1mpqurQBs1nLM9rXkC39zFvgbsp8iDmc1xpgTE9cFoBy7aOj6nYbyO41dv3OG61cqyT8A7NKT+MFfj7ey2jPP35hftSq5j3xpjDHRKT4KwL6dNJe1nOrazKmyibqSRj3Xn1QMvtkDrPdX5jt/Axb767LEX4eftYa15xtjYlp8FIDP7md6kTcB2K8J/KYnM8/fmDX+6qzWZFb6a7KbkxwOaYwxhSukAiAiHYFnATfwqqqOyrFcgssvBjKAvqq65Fjbikg5YAqQDGwArlbVnSd+SrlIuZ7rfqjMOj2ZTZpkn+yNMYYQxgQWETcwFugE1Ad6iEj9HKt1AuoEf/oDL4Ww7RBgrqrWAeYGp8OjajO+8p/BRq1kb/7GGBMUyrthS2Cdqq5X1UzgXaBLjnW6AG9owEKgjIhUyWPbLsCk4OtJwOUndirGGGPyI5QmoKrAxmzTaUCrENapmse2lVR1C4CqbhGRirkdXET6E7iqANgjImtDyHw0FYC/T2D7aGHnGVvsPGNLvs9THj/hY56S28xQCkBu9z5qiOuEsu0xqep4YHx+tjkaEUlV1ZSC2Fcks/OMLXaesSWSzjOUJqA0oHq26WrA5hDXOda2W4PNRAR/bws9tjHGmBMVSgFYBNQRkZoikgh0B2bmWGcm0FsCWgP/Bpt3jrXtTKBP8HUf4MMTPBdjjDH5kGcTkKp6RWQgMIfArZwTVHWViAwILh8HzCZwC+g6AreBXnesbYO7HgVMFZEbgD+Bqwr0zHJXIE1JUcDOM7bYecaWiDlPUc1Xk7wxxpgYYTfFG2NMnLICYIwxcSouC4CI3CIia0VklYg84XSecBKRQSKiIlLB6SzhICJPisjPIrJcRD4QkTJOZyooItIx+O90nYiE70l5h4lIdRH5SkTWBP9P3uZ0pnAREbeI/CQis5zOAnFYAESkHYGnkBuragNgtMORwkZEqgPtCXzJHqs+BxqqamPgF2Cow3kKRIhdsMQKL3CXqtYDWgM3x/C53gascTrEQXFXAIAbgVGqegBAVWP5+YMxwN3k8+G7aKKqn6mqNzi5kMCzJrEglC5YYoKqbjnYeaSq7ibwBlnV2VQFT0SqAZcArzqd5aB4LAB1gXNF5AcR+UZEWjgdKBxEpDOwSVWXOZ2lEF0PfOJ0iAJytO5VYpqIJANnAD84HCUcniHwgczvcI5DYnI8ABH5Aqicy6JhBM65LIFLzRYEnkWopVF4P2we53kv0KFwE4XHsc5TVT8MrjOMQFPC5MLMFkYn3I1KtBGREsB04HZV3eV0noIkIpcC21R1sYi0dTjOITFZAFT1wqMtE5EbgfeDb/g/ioifQOdM6YWVr6Ac7TxFpBFQE1gWGKqBasASEWmpqn8VYsQCcay/TwAR6QNcClwQjYX8KELpgiVmiEgCgTf/yar6vtN5wuBsoLOIXAwUBUqJyFuq2svJUHH3IFjwCeaTVfUBEalLYCyCGjH0xnEEEdkApKhqzPW0GBxw6GngPFWNuiJ+NCLiIfCl9gXAJgLdqvwn25P0MSM4oNQkYIeq3u5wnLALXgEMUtVLHY4Sl98BTABqichKAl+s9YnlN/848AJQEvhcRJaKyDinAxWE4BfbB7tRWQNMjcU3/6CzgWuB84N/h0uDn5RNmMXdFYAxxpiAeLwCMMYYgxUAY4yJW1YAjDEmTlkBMMaYOGUFwBhj4pQVAGOMiVNWAIwxJk79H8bolHxVyapTAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#check visually that the two samples at least appear Gaussian. \n", "#Doing something like a Kolmogorov-Smirnov test would be a more objective way to do this...\n", "plt.figure()\n", "plt.hist(x1, density = True, label='sample 1')\n", "xrange = np.arange(np.min(x1), np.max(x1)+0.1, 0.1)\n", "plt.plot(xrange, norm.pdf(xrange, mu1, np.sqrt(var1)), label='gaussian fit')\n", "plt.legend()\n", "\n", "plt.figure()\n", "xrange = np.arange(np.min(x2), np.max(x2)+0.1, 0.1)\n", "plt.hist(x2, density = True, label='sample 2')\n", "plt.plot(xrange, norm.pdf(xrange, mu2, np.sqrt(var2)), label='gaussian fit')\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x1: Max. Difference = 0.1361100216491875, p-value = 0.044587556750829975\n", "x2: Max. Difference = 0.21001978411137656, p-value = 0.00023906536465364013\n" ] } ], "source": [ "# Use Kolmogorov-Smirnov normality test:\n", "# compare our data sample with a TRUE normal\n", "import scipy.stats as stats\n", "D, p = stats.kstest(x1,'norm')\n", "print('x1: Max. Difference = ' + str(D) + ', p-value = ' + str(p))\n", "D, p = stats.kstest(x2,'norm')\n", "print('x2: Max. Difference = ' + str(D) + ', p-value = ' + str(p))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "F-statistic: 0.24626867737245473\n", "p value: 1.0577726440922956e-11\n", "null hypothesis rejected: variances are statistically different to 95% confidence.\n" ] } ], "source": [ "#calculate f-stat for the two samples\n", "fstat = var1/var2\n", "print('F-statistic:',fstat)\n", "\n", "#dof for both samples\n", "dof1 = len(x1)-1 \n", "dof2 = len(x2)-1 \n", "\n", "#p value of this f-statistic\n", "p = f.cdf(fstat, dof1, dof2) #p-value of F statistic \n", "\n", "#set significance\n", "alpha = 0.025 #two tailed, 95%\n", "print('p value:', p)\n", "\n", "#Do test\n", "if (1-p" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXyV1bno8d+TnZGMZGBKGMIgEiCEGBlVwFIFRahWj1CrV3t6uPSItnrrLfece21P7WCtRz1UTynHCrWnTq0TbXECRRFFGQwoo0lIIEwJCQmZw85e94+9d9yEneydZCd7er6fTz7Z+13rfffz5iUPK+td71pijEEppVToivB3AEoppfqWJnqllApxmuiVUirEaaJXSqkQp4leKaVCXKS/A3AnPT3djBo1yt9huHXo0CEAxo8f7+dIvHDOHitJQRBroNOfpQpwu3btOmOMyXBXFpCJftSoUezcudPfYbg1d+5cALZs2eLXOLyyaa79+/wt/owiNOjPUgU4ESnrrEy7bpRSKsRpoldKqRCniV4ppUJcQPbRK6X85/z585SXl9Pc3OzvUJQbsbGxZGVlERUV5fU+muiVUhcoLy8nMTGRUaNGISL+Dke5MMZQVVVFeXk52dnZXu+nXTdKqQs0NzeTlpamST4AiQhpaWnd/mtLE71S6iKa5ANXT66Nx0QvIsNF5D0ROSAi+0Tk+27qiIisFpEiEdkrIvkuZQtE5JCjbFW3I1RBwxjDp0eq+d37xbx3sAKbTafAVioQeNOitwL/yxgzAZgB3C0iOR3qLATGOb6WA78FEBEL8JSjPAdY5mZfFQJarTYe+Mte/uF3H/PLNw5y1/od/I91n1LbdN7foakgZLFYyMvLa/8qLS3tsn5paSmTJk0CYOfOndx7771d1n3uuec6LYuLiyMvL4+cnBxWrFiBzWZr3z516lQmTJjAtGnT+MMf/tC+3/r168nIyGiP94477uj+SfchjzdjjTEngZOO13UicgDIBPa7VFsCPGvsq5hsF5EUERkKjAKKjDElACLygqOu674qBDzy5kH+squce68eyx2zRvHmF6f4t7/u457nP2PdnZdjidCuAOW9uLg4CgsLe7RvQUEBBQUFnZY7E/23vvUtt+VjxoyhsLAQq9XK1VdfzWuvvUZ+fj5jxozhs88+A6CkpISbbroJm83GXXfdBcCtt97Kk08+2aOY+1q3+uhFZBQwFfikQ1EmcMzlfbljW2fb3R17uYjsFJGdlZWV3QlL+dn2kiqe/vAId8wcyf3XjCc9IYZvzxjJQ0sm8cHhSn7/YYm/Q1QhaNeuXUyZMoWZM2fy1FNPtW/fsmULixYtAuD9999vb2VPnTqVuro6Vq1axdatW8nLy+Pxxx/v9PiRkZHMmjWLoqKii8pGjx7NY489xurVq31/Yn3A6+GVIpIAvAz8wBhzrmOxm11MF9sv3mjMWmAtQEFBgXbuBgljDI+8eZAhSbH8y3UTLihbOm0E7+w/zX9s+pJv5GUyKCnWT1Gqnvq3v+5j/4mOv+69kzMsiR/fMLHLOk1NTeTl5QGQnZ3Nq6++elGdu+66i9/85jfMmTOHBx54wO1xHn30UZ566ilmz55NfX09sbGxPPzwwzz66KP87W9/6zKGxsZGNm/ezE9/+lO35fn5+Rw8eLD9/YsvvsiHH34IwPe///32ln4g8KpFLyJR2JP8n4wxr7ipUg4Md3mfBZzoYrsKEVu/PMPuozV8f/44YqMsF5X/v0U5tFhtrHlfW/XKe86um8LCQrdJvra2lpqaGubMmQPA7bff7vY4s2fP5v7772f16tXU1NQQGem5bVtcXExeXh6zZ8/m+uuvZ+HChW7rdVxv+9Zbb22POZCSPHjRohf7WJ7fAweMMY91Um0DsNLRBz8dqDXGnBSRSmCciGQDx4GlgPuOMRWU/ri9jPSEaL6Zn+W2fFR6PIvzhvH8p0dZefVYUuOj+zlC1RueWt796a677uKzzz5j2LBhPPfcc14NM1y1ahXXX389GzduZMaMGWzatMnjPs4+ek8+++wzJkyY4LFeIPCm62Y2cDvwuYg4z/5fgBEAxpg1wEbgOqAIaATucpRZRWQl8BZgAZ4xxuzz6RkovzlZ28TmA6dZMWcM0ZGd/3H4vTljeGX3cZ7/9Ch3zxvbjxGqULJu3boL3icnJ/Phhx9yxRVX8Kc//cntPsXFxUyePJnJkyfz8ccfc/DgQYYPH05dXV2vYiktLeWHP/wh99xzT6+O01+8GXXzIe772l3rGODuTso2Yv+PQIWY1wtPYDOw9PIRXdYbNziR6dmpvLjjGN+bM4YIHYGjfGDdunV85zvfYcCAAVx77bVu6zzxxBO89957WCwWcnJyWLhwIREREURGRjJlyhTuvPNO7rvvPq8+r7i4mKlTp9Lc3ExiYiL33HNPwHXRdEY69jMFgoKCAqMLj/hAHy+WsfjJDxHg9ZVXeKz7euFxvv9CIf/9j9O5Ylx6n8TTp8Jo4ZEDBw4ETZdEuHJ3jURklzHG7bhSnQJB9cix6kb2ltdy3eShXtW/duIQEmIieb3weB9HppTqSBO96pF39p8GYMGkIV7Vj42ycM3Ewby57xQt1ra+DE0p1YEmetUj7x+uZHRGPCPT4r3eZ/GUYdQ1W3n/kD4Qp1R/0kSvuq35fBvbS6qYc4nbBec7NXtsOomxkWw+UNFHkSml3NFEr7rtkyPVtFht3U70UZYIrrokg/cO6cyWSvUnTfSq2z4uriLKIkzPTuv2vlePH0RFXQv7fPxYvVKqc5roVbd9cqSK3KwU4qIvnvLAk7njMxCBdw9q943qXHenKXbasmULH330UbfrrVmzhmeffdZn9QONrhmruqWx1crn5bUsv2p0j/ZPS4ghb3gK7x6q4Pvzx/k4OhUqejpN8ZYtW0hISGDWrFndqrdixQqf1g802qJX3bK7rAarzTB9dPe7bZyuHj+IPcdqqKxr8WFkKtysXr2anJwccnNzWbp0KaWlpaxZs4bHH3+cvLw8tm7dyl//+lemT5/O1KlTmT9/PqdPn3Zb7yc/+QmPPvqo18d1rV9UVMT8+fOZMmUK+fn5FBcX+/PH4pa26FW37D56FhHIH5HS42NcPWEQ//7OYd4/XMnNl7mfDE0Fhh/84Ac9XgCkM3l5eTzxxBNd1vFmmuKHH36YI0eOEBMTQ01NDSkpKaxYsYKEhAR++MMfAnD27Fm2b9+OiPD000/zyCOP8O///u8X1du8eXO3juta/7bbbmPVqlXceOONNDc3Y7PZevcD6gOa6FW37C2vYXR6PImxUT0+Rs7QJNLio/mo+IwmeuWWN103ubm53HbbbXzjG9/gG9/4hts65eXl3HrrrZw8eZLW1lays7M9frY3x3Wqq6vj+PHj3HjjjQDExgbmmgua6JXXjDHsKa/lyrG9m6tGRJgxOo3txVUYY3q0qr3qH55a3v3JdZrijRs38ve//50PPviADRs28NBDD7Fv38UT495zzz3cf//9LF68mC1btvCTn/zE4+d4c1ynQJwrzB3to1deO3Wumcq6FnKzknt9rBlj0jhR28zR6kYfRKbCwbp16ygsLGTjxo3YbDaOHTvGvHnzeOSRR6ipqaG+vp7ExMQLpiCura0lM9O+eqnrYt4d6zl5e1ynpKQksrKyeO211wBoaWmhsTHw/k1rolde23OsBoDc4T3vn3eaOToVsI/JV6q72tra+Pa3v83kyZOZOnUq9913HykpKdxwww28+uqrF9w0veWWW7jyyitJT//qL9GO9bp7XFd//OMfWb16Nbm5ucyaNYtTp07128/BW9p1o7y2p7yWyAghZ2hSr481JiOBjMQYtpdUsXRa1/PZq/BTX1/fZXlUVFT7+qyuLrnkEvbu3XvBtiVLlnisd+WVV7a/9ua4rvXHjRvHu+++22W8/ubNUoLPAIuACmPMJDflDwC3uRxvApBhjKkWkVKgDmgDrJ3NlayCw97yGsYPSXS7Nmx3OfvpPy7Rfnql+po3XTfrgQWdFRpjfm2MyTPG5AH/B3jfGFPtUmWeo1yTfBCz2Qx7y2vJzep9t43TjNGpnD7XwpEzDT47plLqYh4TvTHmA6DaUz2HZcDzvYpIBaTSqgbqmq1M8cGNWKeZjoeutpd4+89LKdUTPrsZKyIDsLf8X3bZbIC3RWSXiCz3sP9yEdkpIjsrK3W+8kCzt7wWwKct+uz0eNITYthZpoleqb7ky1E3NwDbOnTbzDbG5AMLgbtF5KrOdjbGrDXGFBhjCjIyujf9rep7B06eI8oijBuc4LNjigiXjUxhd9lZnx1TKXUxXyb6pXTotjHGnHB8rwBeBab58PNUPzp4qo4xGQlEWXw7IveykQMprWrUeW+U6kM++a0VkWRgDvC6y7Z4EUl0vgauAb7wxeep/nf4dB2XDkn0+XEvGzkQsM+ho1R3bNiwgYcffhiA1157jf3797eXPfjgg2zatKnbxywtLWXSpIsGF/rtOL7izfDK54G5QLqIlAM/BqIAjDFrHNVuBN42xrgOnxgMvOoYNhcJPGeMedN3oav+Utt4npO1zYwf0vvx8x1NHJZMtCWC3WVnuXaidwuNK2W1Wlm8eDGLFy8G7Il+0aJF5OTkAPDTn/7Un+EFHG9G3Swzxgw1xkQZY7KMMb83xqxxSfIYY9YbY5Z22K/EGDPF8TXRGPPzvjgB1fcOnbY/+t0XLfrYKAuTs5LZpf30ysWzzz5Lbm4uU6ZM4fbbbwfgzjvv5P7772fevHn86Ec/Yv369axcuZKPPvqIDRs28MADD5CXl0dxcTF33nknf/nLXwDYsWMHs2bNYsqUKUybNo26ujpKS0u58soryc/PJz8/3+NiJbfeeisbN25sf3/nnXfy8ssve3UcZ5xOixYtYsuWLQC8/fbbzJw5k/z8fG655Zb2B8VWrVrVPlWyc8bM3tAnY5VHh07Zl/0b3weJHuzdN+s/KqXF2kZMZO8fxlI+tOsHcNa30xQzMA8u63yytH379vHzn/+cbdu2kZ6eTnX1V+M7Dh8+zKZNm7BYLKxfvx6AWbNmsXjxYhYtWsTNN998wbFaW1u59dZbefHFF7n88ss5d+4ccXFxDBo0iHfeeYfY2Fi+/PJLli1bxs6dOzuNaenSpbz44otcd911tLa2snnzZn77299ijOnWcVydOXOGn/3sZ2zatIn4+Hh+9atf8dhjj7Fy5UpeffVVDh48iIhQU1Pj1fG6ooleeXTodB2JsZEMTe6bKVjzRwxk7QclfHH8XHufvQpf7777LjfffHP73DSpqantZbfccgsWi/eNgUOHDjF06FAuv/xywD4JGUBDQwMrV66ksLAQi8XC4cOHuzzOwoULuffee2lpaeHNN9/kqquuIi4ujtra2m4dx9X27dvZv38/s2fPBuz/Kc2cOZOkpCRiY2P57ne/y/XXX8+iRYu8PmZnNNErjw6dqmP84MQ+m6ag/YZs2VlN9IGmi5Z3X+lqSoz4+HifHOvxxx9n8ODB7NmzB5vN5nEe+djYWObOnctbb73Fiy++yLJly7w+TmRk5AWLkTQ3N7fH9vWvf53nn7/4GdNPP/2UzZs388ILL/Dkk0/2ei4dnb1SdckYw8FTdX3WbQOQkRhD1sA4Cst7/yeqCn5f+9rXeOmll6iqss9s6tp105nOphG+9NJLOXHiBDt27ADsC4VYrVZqa2sZOnQoERER/PGPf6Strc3jZyxdupR169axdetWrr32WgCvjjNq1CgKCwvbp0D+9NNPAZgxYwbbtm2jqKgIgMbGRg4fPkx9fT21tbVcd911PPHEEz5Z4Utb9KpLJ2ubqWu29smNWFdTslLap0FW4W3ixIn867/+K3PmzMFisTB16tT2/vjOLF26lH/6p39i9erV7TdhAaKjo3nxxRe55557aGpqIi4ujk2bNvHP//zPfPOb3+TPf/4z8+bN8+ovhWuuuYY77riDxYsXEx0dDeDVcWbPnk12djaTJ09m0qRJ5OfnA5CRkcH69etZtmwZLS3250h+9rOfkZiYyJIlS2hubsYYw+OPP+7tj65TEogrpBQUFBhvb2j0t7lz5wK03zUPaJvm2r/P39LjQ7x3qIK71u3gpf85k2nZqZ536KG1HxTzi40H2fV/55OWENNnn9NjPvhZBosDBw4wYcIEf4ehuuDuGonIrs4mj9SuG9Wlw6fsfw5f4sOpD9yZ4phDxzmnjlLKdzTRqy4VV9aTnhBDyoDoPv2cSZnJRAgUaveNUj6niV51qbiygTEZ3Rvp0BPxMZGMG5TIXr0hGxACsUtX2fXk2miiV10qqaxndEbfdts4TRmezJ7yWk0yfhYbG0tVVZVehwBkjKGqqsrjcNCOdNSN6lR1QytnG8/3S4se7HPdv7SznPKzTQxPHdAvn6kulpWVRXl5ObouRGCKjY0lKyurW/tooledKqm0z7sxpp9a9HnD7Tdk95TXaKL3o6ioKLKzs/0dhvIh7bpRnSru50Q/fkgi0ZEROp5eKR/TRK86VVLZQHRkBJkD4/rl86IsEUwclsSeYzrEUilf0kSvOlVcWU92WjyWiL6Z48adKVkpfH68FmubzXNlpZRXPCZ6EXlGRCpExO3qUCIyV0RqRaTQ8fWgS9kCETkkIkUissqXgau+V1LZwOh+uhHrlDc8habzbRQ5uo2UUr3nTYt+PbDAQ52txpg8x9dPAUTEAjyFfWHwHGCZiOT0JljVf1qtNsqqG/utf94pNysZQPvplfIhb1aY+gDwPH3cxaYBRY6VplqBF4AlPTiO8oOj1Y202Uy/t+hHpcWTFBtJofbTK+Uzvuqjnykie0TkDRGZ6NiWCRxzqVPu2OaWiCwXkZ0islPH7/pff4+4cYqIECZnJfPFcU30SvmKLxL9bmCkMWYK8BvgNcd2d3fwOn3Uzhiz1hhTYIwpyMjI8EFYqjdKKu3rvPd3ix5gcmYKB0+do8XqeY5wpZRnvU70xphzxph6x+uNQJSIpGNvwQ93qZoFnOjt56n+UVbVQHpCNImxUf3+2blZyZxvMxw6dfFCEkqp7ut1oheRIeJYq0tEpjmOWQXsAMaJSLaIRANLgQ29/TzVP8qqGhnhp6dTJ2fab8jqlMVK+YbHKRBE5HlgLpAuIuXAj4EoAGPMGuBm4HsiYgWagKXGPhuSVURWAm8BFuAZY8y+PjkL5XNHqxv7dKGRrmQNjGPggCg+10SvlE94TPTGmGUeyp8EnuykbCOwsWehKX9psbZxorbJby16EWFyVgp79YasUj6hT8aqi5SfbcIY/JboAXIzkzl8uo7m83pDVqne0kSvLnK0qhGAkWn+S/STs5Jpsxn2nzzntxiUChWa6NVFyqrsQytH+DHRO5+Q1X56pXpPE726SFl1IwOiLWQkxPgthiFJsaQnxOjIG6V8QBO9ushRx9BKx6hZvxARcrOS+fy4znmjVG9polcXKav23xh6V5MzkymqqKehxervUJQKapro1QVsNsPR6ka/3oh1ys1KxmbQG7JK9ZImenWB03XNtFptjEjr/zluOnI+Ias3ZJXqHU306gJlzqGVAdB1MygpliFJsXyuD04p1Sua6NUFAmEMvavJWcnsLdcbskr1hiZ6dYGy6gYsEcKwlP5ZENyT3MxkSs40UNd83t+hKBW0NNGrC5RVNZKZEkeUJTD+aUzOSsYY2HdCb8gq1VOB8dusAkagjLhx0huySvWeJnp1AX/OQ+9OWkIMmSlxOpOlUr2giV61q208T23T+YBq0YN9PP3nekNWqR7zmOhF5BkRqRCRLzopv01E9jq+PhKRKS5lpSLyuYgUishOXwaufO9otX3EzYhU/4+hdzU5K5nSqkZqG/WGrFI94U2Lfj2woIvyI8AcY0wu8BCwtkP5PGNMnjGmoGchqv5SVm2ftTLgWvSZKQB8cUK7b5TqCY+J3hjzAVDdRflHxpizjrfbsS8CroKQ82GpQOqjB11DVqne8nUf/T8Cb7i8N8DbIrJLRJb7+LOUjx2taiQ9IYb4GI8rTPar5AFRjEwboDNZKtVDPvuNFpF52BP9FS6bZxtjTojIIOAdETno+AvB3f7LgeUAI0aM8FVYqhvKqhsCrtvGaXJmMoXHNNEr1RM+adGLSC7wNLDEGFPl3G6MOeH4XgG8Ckzr7BjGmLXGmAJjTEFGRoYvwlLddLSqMSDmuHEnNyuZ8rNNVDe0+jsUpYJOrxO9iIwAXgFuN8YcdtkeLyKJztfANYDbkTvK/1qsbZw81+zX5QO7MtlxQ1YnOFOq+zx23YjI88BcIF1EyoEfA1EAxpg1wINAGvCfjhWJrI4RNoOBVx3bIoHnjDFv9sE5KB84Vt2EMYE34sZpUmYSAJ+X1zDnEv2LT6nu8JjojTHLPJR/F/ium+0lwJSL91CB6KhjaGWgjaF3SoyNYnRGvI68UaoH9MlYBbjMQx+gLXqwz2SpXTdKdZ8megXYE318tIW0+Gh/h9KpyVkpnKxtpqKu2d+hKBVUNNErwD79wYi0eBz3VAJSbpb9wakvtFWvVLdoolcAlFU1BOzQSqecoUlEiD4hq1R3aaJX2GyGY2ebArp/HiA+JpKxgxJ0bnqlukkTveLUuWZarbaAHUPvanJmCnuP12KM8XcoSgUNTfTqqxE3ATq00lVuVjKVdS2cPtfi71CUChqa6FX7GPpA77oB+9z0AHt1IRKlvKaJXlFW1UhkhDA0OdbfoXiUMzQJS4ToeHqlukETvaKsupGsgXFEWgL/n0NslIVLBifqyBuluiHwf7NVnztaZR9DHyycT8jqDVmlvKOJXgXFGHpXk7OSqW5o5XhNk79DUSooaKIPczWNrZxrtgbFjVgn5xOyOp5eKe9oog9zgbpObFfGD0kkyiLs1RuySnlFE32YK6t2JPogatHHRFq4dEiStuiV8pIm+jB3tMo5D33wJHqw99PvLa/RG7JKecFjoheRZ0SkQkTcLgModqtFpEhE9opIvkvZAhE55Chb5cvAlW+UVTWSkRjDgGifrRPfL3IzkznXbOWo4y8SpVTnvGnRrwcWdFG+EBjn+FoO/BZARCzAU47yHGCZiOT0Jljle2XVgbsgeFe+ekJWu2+U8sRjojfGfABUd1FlCfCssdsOpIjIUGAaUGSMKTHGtAIvOOqqAGIfQx98if6SwYlER0boE7JKecEXffSZwDGX9+WObZ1td0tElovIThHZWVlZ6YOwlCfN59s4da45KCYz6yjKEkHO0CQKj+mcN0p54otE725JItPFdreMMWuNMQXGmIKMjAwfhKU8OVYd+OvEdmXqiBQ+L6/lfJvN36EoFdB8kejLgeEu77OAE11sVwGifQx9kCb6/BEDaTrfxsGTdf4ORamA5otEvwG4wzH6ZgZQa4w5CewAxolItohEA0sddVWAcI6hD8absQD5IwcCsPvoWT9HolRg82Z45fPAx8B4ESkXkX8UkRUissJRZSNQAhQB/wX8M4AxxgqsBN4CDgAvGWP29cE5qB46Vt1IQkwkqfHR/g6lR4YlxzIkKVYTvVIeeBw8bYxZ5qHcAHd3UrYR+38EKgCVVTUwInUAIu5upwQ+ESF/ZAq7yjTRK9UVfTI2jJVVNwbtjVin/BEDKT/bRMW5Zn+HolTA0kQfptpshvLqpqC9Eeuk/fRKeaaJPkydOtdMa5stKMfQu5o4LIloSwS7j+p4eqU6o4k+TJVVBc+C4F2JibQwKTOJ3dpPr1SnNNGHqaNBOA99Zy4bOZC9x2tpteqDU0q5o4k+TJVVNxJlEYalxPk7lF7LHzGQVquNfSd03hul3NFEH6bKqhrIGjgAS0RwDq109dUNWe2nV8odTfRhqvRMI6OCvH/eaXBSLJkpcdpPr1QnNNGHIWMMpVUNjEoP7hE3rgpGDWRHabWuOKWUG5row1BlXQuNrW1kh1Cin5adSkVdC6VVuuKUUh1pog9DR844h1aGTqKfnp0GwCclVX6ORKnAo4k+DJU6xtBnh1CiH5MRT3pCNJ8e6WoxNKXCkyb6MHTkjHNoZay/Q/EZEWFadiqfaKJX6iKa6MNQWVUDw1MHEGkJrcs/PTuN4zVN7StnKaXsQus3XXnlyJmGkOq2cZo+OhVAW/VKdeBVoheRBSJySESKRGSVm/IHRKTQ8fWFiLSJSKqjrFREPneU7fT1CajuMcZQVtUYUjdinS4ZlEjKgCg+PaI3ZJVy5XHhERGxAE8BX8e+DuwOEdlgjNnvrGOM+TXwa0f9G4D7jDGuzap5xpgzPo1c9cjpcy00nW8jOz00HpZyFREhXD5K++mV6sibFv00oMgYU2KMaQVeAJZ0UX8Z8LwvglO+5xxaGUoPS7manp1KWVUjp2p1IRKlnLxJ9JnAMZf35Y5tFxGRAcAC4GWXzQZ4W0R2icjyngaqfMM5PfGoEOy6AZgx2jGeXrtvlGrnTaJ3N+tVZ8+Z3wBs69BtM9sYkw8sBO4WkavcfojIchHZKSI7KysrvQhL9cSRqgaiLREhMWulOxOGJpEYE6ndN0q58CbRlwPDXd5nASc6qbuUDt02xpgTju8VwKvYu4IuYoxZa4wpMMYUZGRkeBGW6onSMw0MT40LiVkr3bFECAWjBrK9WFv0Sjl5k+h3AONEJFtEorEn8w0dK4lIMjAHeN1lW7yIJDpfA9cAX/gicNUzpWcaQ2qOG3dmj02n5EwDx2ua/B2KUgHBY6I3xliBlcBbwAHgJWPMPhFZISIrXKreCLxtjGlw2TYY+FBE9gCfAn83xrzpu/BVd9hsjlkrQ7R/3umqS+x/EX74pXYBKgVeDK8EMMZsBDZ22Lamw/v1wPoO20qAKb2KUPnM6bpmWqy2kB1x4zRuUAKDk2L44Msz3Hr5CH+Ho5Tf6ZOxYaR9aGWIt+hFhCvHZbCt6AxtNp2fXilN9GGk9Ix9DphRIfiwVEdXjkunpvE8XxzXdWSV0kQfRkqrGoiOjGBYcmgOrXR1xdh0ALZqP71SmujDyZEzDYxMHUBEiA6tdJWWEMOkzCQ++FJn3lBKE30YKa6oZ0xGgr/D6DdXjM1gd9lZ6lus/g5FKb/SRB8mWq02yqobGTsofBL9VePSsdqMLi+owp4m+jBRVtVAm80wZlBoj7hxdZ6Q74cAABFeSURBVNmogcRGRbBVu29UmNNEHyaKKuoBGJuR6OdI+k9MpIUZo9P44LDekFXhTRN9mHAm+nBq0QPMGz+IkjMNlFTW+zsUpfxGE32YKKqsJzMljgHRXj0MHTK+NmEQAJsPVPg5EqX8RxN9mCiurGdMGN2IdcoaOIBLhySy6cBpf4eilN9oog8DNpuhuKKBsWE0tNLV/AmD2Vl2lrMNrf4ORSm/0EQfBk7UNtF0vi2shla6mp8zmDabYcth7b5R4UkTfRhovxGbEV43Yp1yM5PJSIxhk/bTqzCliT4MFFfaZ60M1xZ9RITwtUsH8f6hSlqsbf4OR6l+p4k+DBRV1DNwQBRpCTH+DsVvFkwaQn2Lla2H9eEpFX68SvQiskBEDolIkYisclM+V0RqRaTQ8fWgt/uqvldUUce4QeHzoJQ7s8emkxwXxcbPT/o7FKX6ncdELyIW4ClgIZADLBORHDdVtxpj8hxfP+3mvqqPGODgyTouHRreiT7KEsE1OYN558Bp7b5RYcebFv00oMgYU2KMaQVeAJZ4efze7Kt8oNVqo67FyqVDkvwdit9dN3kodc1WthVp940KL94k+kzgmMv7cse2jmaKyB4ReUNEJnZzX0RkuYjsFJGdlZU6N4mvNDqm6A33Fj3Yu2+SYiP5+95T/g5FqX7lTaJ3t0pFx4U4dwMjjTFTgN8Ar3VjX/tGY9YaYwqMMQUZGRlehKW80dhq76YYP1gTfXRkBF/PGcLb+09p940KK94k+nJguMv7LOCEawVjzDljTL3j9UYgSkTSvdlX9a2G822MTBtAfEx4zXHTmSV5w6hrtvKujqlXYcSbRL8DGCci2SISDSwFNrhWEJEhIiKO19Mcx63yZl/VtxpbrFw6RFvzTrPHpjMoMYZXPjvu71CU6jcem3nGGKuIrATeAizAM8aYfSKywlG+BrgZ+J6IWIEmYKkxxgBu9+2jc1EdtNkMzefbuDRLb8Q6WSKEb0zNZN22I1Q3tJIaH+3vkJTqc179Pe/ojtnYYdsal9dPAk96u6/qH03n7f3QE/RG7AVunJrJ2g9K+NveE9wxc5S/w1Gqz+mTsSHMeSNWh1ZeaMLQJC4dksgru7X7RoUHTfQhrLHVSoQII1IH+DuUgHPzZVkUHqvh8Ok6f4eiVJ/TRB/CGlvbiIu2EBHhbpRreLspP4toSwTPfXLU36Eo1ec00Ycom83Q0GIlQYdVupUaH83CyUN4eXc5Ta06pl6FNk30IarkTD1tNqPj57vwrWkjqGu28te9+miHCm2a6EPUnmO1ANqi78K07FTGDkrgT9p9o0KcJvoQtae8hogIIS7K4u9QApaIcNv0Eew5VkPhsRp/h6NUn9FEH6L2HKshPjoS0fuwXbqlYDiJsZH819YSf4eiVJ/RRB+CWqxtHDhZp902XkiIieRb00bwxucnOVbd6O9wlOoTmuhD0MGTdbS22TTRe+nO2aOIEGH9R6X+DkWpPqGJPgTtLbf3NyfEaqL3xtDkOK7PHcoLnx6ltvG8v8NRyuc00YegwmO1pMVHEx2pl9dbK+aMoaG1jd9/qH31KvRoJghBnx07y5ThKW5XfVHuTRiaxMJJQ1i3rZSaxlZ/h6OUT2miDzEV55opqWxgenaqv0MJOvd+bRx1LVae+fCIv0NRyqc00YeYj0uqAJg5Js3PkQSfCUOTuG7yEJ7ZVkp1g7bqVejwKtGLyAIROSQiRSKyyk35bSKy1/H1kYhMcSkrFZHPRaRQRHb6Mnh1se0l1STGRJIzVKcm7on75l9C0/k2/mPTYX+HopTPeEz0ImIBngIWAjnAMhHJ6VDtCDDHGJMLPASs7VA+zxiTZ4wp8EHMqguflFRxeXYqkRb9Y60nxg1OZNm04fz3J0cpqqj3dzhK+YQ32WAaUGSMKTHGtAIvAEtcKxhjPjLGnHW83Y59EXDVz06fa6bkTAMzR2u3TW/cN/8SBkRZ+MXGA/4ORSmf8CbRZwLHXN6XO7Z15h+BN1zeG+BtEdklIss720lElovIThHZWVlZ6UVYqqPtjv75GZroeyUtIYaVV4/l3YMVbNp/2t/hKNVr3iR6d6P0jNuKIvOwJ/ofuWyebYzJx971c7eIXOVuX2PMWmNMgTGmICMjw4uwVEfbS6pIjI0kZ5j2z/fWXbOzuWRwAg++/gX1LVZ/h6NUr3iT6MuB4S7vs4CLJvAWkVzgaWCJMabKud0Yc8LxvQJ4FXtXkPIxYwzbiqqYnp2KRVeU6rXoyAh+eVMuJ8818+hbh/wdjlK94k2i3wGME5FsEYkGlgIbXCuIyAjgFeB2Y8xhl+3xIpLofA1cA3zhq+DVV4orGzha3cic8YP8HUrIuGzkQG6fMZI/fFzKuWZt1avg5THRG2OswErgLeAA8JIxZp+IrBCRFY5qDwJpwH92GEY5GPhQRPYAnwJ/N8a86fOzULx70N6XfPWlmuh96X8vuJQRqQMoqqjD2ua2x1KpgOfVrFfGmI3Axg7b1ri8/i7wXTf7lQBTOm5Xvrf5QAWXDkkkMyXO36GElISYSP5j6VRa3jSUVDUwzhhEJ/lXQUYHW4eAM/Ut7Cit5us5g/0dSkjKG55C1sA4qutb+INOZayCkCb6EPDmF6ewGbg+d6i/QwlZw5LjGDggmof+foBtRWf8HY5S3aKJPgT8fe9JRmfEM35wor9DCVkiMGZQAmMy4rn7ud2UVOpTsyp4aKIPcsdrmth+pIobcodp33Efi4wQ/uuOAiwifPvpTzhR0+TvkJTyiib6IPfyrnKMgZsv01kn+sPItHj+8J1p1DVb+fbTn1BR1+zvkJTySBN9EGuzGf686xgzR6cxPHWAv8MJG5Myk3nmrss5WdvMP6z5mPKzuqi4Cmya6IPYuwcrOFbdxLdnjPR3KGHn8lGp/Pd3p1Hd0Motaz7m0Kk6f4ekVKc00Qex339YwrDkWK6dqMMq/eGykam8sHwmbTbDTf+5jbf3nfJ3SEq5pYk+SO0orWZ7STV3zc7Wuef9KGdYEhtWXsHYQQks/+MufvnGAVqtNn+HpdQFNEMEIWMMj719mPSEGO22CQBDkmN58X/O5FvTR/C790u46bfbtCtHBRRN9EHorX2n+bikipXzxhAXbfF3OAqIjbLwixsn87vbL+NETTPXr97KLzce0CmOVUDQRB9kzjWf56G/7Wf84ERtzQegaycOYfP9c/hmfha/+6CEK3/1Lms/KKaptc3foakwpok+iBhj+Mnr+zhZ28QvbpqsffMBamB8NL+6OZfX757N5KwUfrHxIFf9+j2eeq+IyroWf4enwpBmiiDyzLZSXvnsOPd+bRyXjRzo73CUB1OGp/Dsd6bx5xUzGT84kV+/dYhZD2/mnuc/Y8uhCr1pq/qNV9MUK/977pOjPPS3/SyYOIR7rx7n73BUN9jH3E+nqKKeP31Sxsu7yvnrnhMkxUby9ZwhfG3CIGaMTiM1PtrfoaoQpYk+wDWfb+NXbx5k3bZS5o3P4ImleUToUoFBaeygBH58w0R+tOBSPvzyDBu/OMk7+0/x8u5yAHKGJlEwaiCTMpOZNCyZcYMTiNLuOeUDXiV6EVkA/AdgAZ42xjzcoVwc5dcBjcCdxpjd3uyr3GtqbeMvu46xdmsJx6qbuHPWKP71+gn6ix8CYqMszM8ZzPycwZxvs7G3vJaPi8/wUXEVr+w+zrMflwH2dWvHZCSQnT6AkWnxZKfFk5Uax6DEWAYnxZAQE6kT2SmveEz0ImIBngK+jn2h8B0issEYs9+l2kJgnONrOvBbYLqX+4Y1Ywz1LVZqGs9TVtXIlxV1fFxcxcfFVdS1WJmSlczDN+Uye2y6v0NVfSDKEsFlIwdy2ciBrLx6HDabobSqgc+P1/LF8VqKKuo5eLKOt/edxmq7cCnDuCgLg5NiSEuIITkuiqTYSJLiokiKjSIpLpKEmChioyKIjbLYv0daiHG+jrIQExlBlCUCS4RgEcFiESIjBEuEEBkRQYSg/5GECG9a9NOAIseygIjIC8ASwDVZLwGeNcYYYLuIpIjIUGCUF/v6zKLfbKX5vA17GND+a2Eu+HZRuWkvNxe+77BEqDGGfUfPYoCZv9zc+X503N99uc0Y6pqttHX4Bc5MieP63KF887IsCkYO1F+2MBIRIYzOSGB0RgJL8jLbt1vbbByvaeJ4TROVdS2cPtdMxbkWKupaqGpoobKuheLKes41neecm39TPWVpT/z274I9+Yvw1Wvs8/XbuZaBcGFdnNs7K+ttwL08gC9+03rz+5o6IJqXVsz0QRQX8ibRZwLHXN6XY2+1e6qT6eW+AIjIcmA5wIgRI7wI62JjMxI471zAWS74dtE/JPFQ/tX+ckH9Y3FRAFw5Lt1tuXQ4wMWf89VxRSAxNpLkuChS4qIZnjqA0RnxDEqM0eSuLhBpiWBkWjwj0+I91jXG0NjaRl2zlebzbTRb22g+b7O/Pm9/3WJto+W8DavN0GZzfrd/WS/4bi+z2Uz775Yx9qaNMfZGjmsDxv7avs213LU+7e+Nyz4uDbMeMh1bZt3dv5ef74uDJMb2zW1Tb47qLuN0PJ3O6nizr32jMWuBtQAFBQU9+nE9sXRqT3brlp1PJQDwyM265rkKTCJCfEwk8TE61kLZefMvoRwY7vI+CzjhZZ1oL/ZVSinVh7wZwrEDGCci2SISDSwFNnSoswG4Q+xmALXGmJNe7quUUqoPeWzRG2OsIrISeAv7EMlnjDH7RGSFo3wNsBH70Moi7MMr7+pq3z45E6WUUm551YlnjNmIPZm7blvj8toAd3u7r1JKqf6jT98opVSI00SvlFIhThO9UkqFOE30SikV4qS3T5P1BRGpBMp6uHs6cMaH4fiTnkvgCZXzAD2XQNXTcxlpjMlwVxCQib43RGSnMabA33H4gp5L4AmV8wA9l0DVF+eiXTdKKRXiNNErpVSIC8VEv9bfAfiQnkvgCZXzAD2XQOXzcwm5PnqllFIXCsUWvVJKKRea6JVSKsQFZaIXkQUickhEikRklZtyEZHVjvK9IpLvjzi94cW5zBWRWhEpdHw96I84vSEiz4hIhYh80Ul5MF0XT+cSFNdFRIaLyHsickBE9onI993UCYrr4uW5BMt1iRWRT0Vkj+Nc/s1NHd9dF2NMUH1hn+64GBiNfWGTPUBOhzrXAW9gX+FqBvCJv+PuxbnMBf7m71i9PJ+rgHzgi07Kg+K6eHkuQXFdgKFAvuN1InA4iH9fvDmXYLkuAiQ4XkcBnwAz+uq6BGOLvn2xcmNMK+BccNxV+2LlxpjtgHOx8kDjzbkEDWPMB0B1F1WC5bp4cy5BwRhz0hiz2/G6DjiAfS1nV0FxXbw8l6Dg+FnXO95GOb46jozx2XUJxkTf2ULk3a0TCLyNc6bjT7w3RGRi/4TWJ4LlungrqK6LiIwCpmJvPboKuuvSxblAkFwXEbGISCFQAbxjjOmz6xKMqwf3ZrHyQONNnLuxz2FRLyLXAa8B4/o8sr4RLNfFG0F1XUQkAXgZ+IEx5lzHYje7BOx18XAuQXNdjDFtQJ6IpACvisgkY4zrPSGfXZdgbNH3ZrHyQOMxTmPMOeefeMa+WleUiKT3X4g+FSzXxaNgui4iEoU9Mf7JGPOKmypBc108nUswXRcnY0wNsAVY0KHIZ9clGBN9bxYrDzQez0VEhoiIOF5Pw37Nqvo9Ut8IluviUbBcF0eMvwcOGGMe66RaUFwXb84liK5LhqMlj4jEAfOBgx2q+ey6BF3XjenFYuWBxstzuRn4nohYgSZgqXHckg80IvI89lEP6SJSDvwY+02moLou4NW5BMt1mQ3cDnzu6A8G+BdgBATddfHmXILlugwF/iAiFuz/Gb1kjPlbX+UxnQJBKaVCXDB23SillOoGTfRKKRXiNNErpVSI00SvlFIhThO9UkqFOE30SikV4jTRK6VUiPv/QthQ5zmV5EwAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "xrange = np.arange(0,3,0.01)\n", "fdist = f.pdf(xrange,dof1,dof2)\n", "lowerbound = f.ppf(alpha,dof1,dof2)\n", "upperbound = f.ppf(1-alpha,dof1,dof2)\n", "\n", "plt.plot(xrange, fdist, label = 'F-dist PDF')\n", "plt.axvline(fstat, label = 'F-statistic', color='black')\n", "plt.axvline(lowerbound, label = 'critical values', color = 'orange')\n", "plt.axvline(upperbound, color = 'orange')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the two standard deviations are statistically distinct and $s_1$ (numerator) < $s_2$ (denominator), then $F$ will appear in the LHS tail of the F-distribution. If $s_1$ > $s_2$, then $F$ will be in RHS tail." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Adjusting DOF for filtering\n", "I assume that the DOF of the samples need to be adjusted if there's autocorrelation and/or filtering has been applied. The effective sample size ($N^*$) of Leith (1973) is\n", "\n", "\\begin{equation}\n", "N^* = N \\frac{\\ln(\\rho(1))}{-2},\n", "\\end{equation}\n", "\n", "where $\\rho(1)$ is the lag-one autocorrelation. See Hartmann's notes for other effective sample approaches (https://atmos.uw.edu/~dennis/552_Notes_6a.pdf). I suspect that low pass filtering may break the normal assumptions for the f-test. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Original DOFs: 99 99\n", "New DOFs: 7.520915488338549 6.536532399300415\n", "\n", "F statistic 0.2153572393716421\n", "p value: 0.027085159636082483\n", "cannot reject null hypothesis: variances are not statistically different to 95% confidence.\n" ] } ], "source": [ "#quick running mean of x1 and x2 to give them some autocorrelation\n", "n = 5 #5-sample running mean\n", "x1_lp = np.convolve(x1, np.ones(n)/n, mode='same')\n", "x2_lp = np.convolve(x2, np.ones(n)/n, mode='same')\n", "#low pass filtering may break normal assumption for the f-test\n", "\n", "#new stats\n", "mu1_lp = np.mean(x1_lp)\n", "mu2_lp = np.mean(x2_lp)\n", "\n", "var1_lp = np.var(x1_lp, ddof = 1)\n", "var2_lp = np.var(x2_lp, ddof = 1)\n", "\n", "#lag-one autocorrelation\n", "autocorr1 = np.corrcoef(x1_lp[0:-1],x1_lp[1:])[0,1]\n", "autocorr2 = np.corrcoef(x2_lp[0:-1],x2_lp[1:])[0,1]\n", "\n", "#effective samples\n", "nstar1 = len(x1)*np.log(autocorr1)/-2\n", "nstar2 = len(x2)*np.log(autocorr2)/-2\n", "\n", "#new DOF\n", "dof1_lp = nstar1-1 #does 1 need to be subtracted? I think so, but not 100% sure\n", "dof2_lp = nstar2-1\n", "print('Original DOFs:', dof1, dof2)\n", "print('New DOFs:', dof1_lp, dof2_lp)\n", "print()\n", "\n", "#recompute f-stat\n", "fstat_lp = var1_lp/var2_lp\n", "print('F statistic', fstat_lp)\n", "\n", "#p value of this f-statistic\n", "p_lp = f.cdf(fstat_lp, dof1_lp, dof2_lp) #p-value of F statistic \n", "\n", "print('p value:', p_lp)\n", "\n", "#Do test\n", "if (1-p_lp" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD5CAYAAAA3Os7hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3xU9Z3/8ddnJpfJPeTCNUDCTbklgAERpGitLVQK2lWBWlvc7bpsi7Zldcv+drfbte6uWrdaqlvq2opaK1itlrasF1C8oChBLpV7uIdrLiTkQjKZzPf3x0zCEHKZhJmcmZPP8/HgMTNnvnPmM0l455vvOef7FWMMSimlop/D6gKUUkqFhga6UkrZhAa6UkrZhAa6UkrZhAa6UkrZhAa6UkrZREwwjURkFvAzwAk8bYx5qNXz9wN3BOxzNJBtjKlob59ZWVkmNze3OzW3a++urQBcMWZiSPerVItze323qVdYW4fqtbZs2VJmjMlu6znp7Dx0EXEC+4AbgRJgM7DQGLOrnfZfAb5vjPl8R/stLCw0RUVFQZQfvOsK0gHYsL0ypPtVqsW663y3X9hgZRWqFxORLcaYwraeC2bIZQpQbIw5aIxxA6uAeR20Xwi82PUylVJKXY5gAn0QcCzgcYl/2yVEJBGYBbxy+aUppZTqimACXdrY1t44zVeAje2NnYvI3SJSJCJFpaWlwdaolFIqCMEcFC0BBgc8zgFOtNN2AR0MtxhjngKeAt8YepA1KqXCoLGxkZKSEurr660uRbXB5XKRk5NDbGxs0K8JJtA3AyNFJA84ji+0v9a6kYikATOBrwf97kopy5SUlJCSkkJubi4ibf0hrqxijKG8vJySkhLy8vKCfl2nQy7GGA+wBHgD2A28ZIzZKSKLRWRxQNNbgDeNMbVdrF0pZYH6+noyMzM1zCOQiJCZmdnlv56COg/dGLMWWNtq24pWj1cCK7v07kopS2mYR67ufG965ZWih8pqefr9g+w8UWV1KUopFTK9LtAPlNbwlZ9/wIN/3s3cJzayYe8Zq0tSqldyOp1MmDCh5d/hw4c7bH/48GHGjRsHQFFREffee2+HbX/729+2+1xCQgITJkxgzJgxLF68GK/X27J94sSJjB49milTpvDss8+2vG7lypVkZ2e31PuNb3yj6x86zIIacrGTB/64C4fAH74znR+8soP7freDd+6bSYor+CPJSqnLl5CQwLZt27r12sLCQgoL27xYErgQ6F/72iXnbwAwfPhwtm3bhsfj4fOf/zyvvfYakyZNYvjw4Wzd6ptC5ODBg3z1q1/F6/Vy1113ATB//nyeeOKJbtXcE3pVD/1QWS3v7ivlb2cMo2BwOg//VT5lNQ0899ERq0tTSrVhy5YtFBQUcM011/Dkk0+2bN+wYQNz5swB4N13323pNU+cOJHq6mqWLVvG+++/z4QJE3jsscfa3X9MTAzTpk2juLj4kueGDRvGT3/6U5YvXx76DxYmvaqH/sftJxCB+ZN9p9UXDE7n2hFZPP/REe7+3DBinb3q95tSLf79jzvZdeJcSPc5ZmAq//aVse0+f/78eSZMmABAXl4er7766iVt7rrrLn7+858zc+ZM7r///jb38+ijj/Lkk08yffp0ampqcLlcPPTQQzz66KP86U9/6rDGuro61q9fzwMPPNDm85MmTWLPnj0tj1evXs0HH3wAwHe/+92Wnnuk6FUJ9sH+MsYOTKVvqqtl213Tczl1rp639+hYulI9qXnIZdu2bW2GeVVVFZWVlcycOROAO++8s839TJ8+naVLl7J8+XIqKyuJiem8n3rgwAEmTJjA9OnTuemmm5g9e3ab7VpPXjh//vyWmiMtzKEX9dBrGjx8evQsf/u5YRdtnzkqm8ykOP64/QRfGtvfouqUslZHPemedNddd7F161YGDhzIb3/726BO3Vu2bBk33XQTa9euZerUqaxbt67T1zSPoXdm69atjB49OqjaI0GvCfTNhyvweA3Xjsi6aHuM08Hs8f15eUsJtQ0ekuJ7zZdEqYjzzDPPXPQ4LS2NDz74gGuvvZYXXnihzdccOHCA8ePHM378eD766CP27NnD4MGDqa6uvqxaDh8+zH333cc999xzWfvpSb1myOWzkipEfOPmrd00fiD1jV4+KC6zoDKlVHueeeYZvvOd73DNNdeQkJDQZpvHH3+ccePGUVBQQEJCArNnzyY/P5+YmBgKCgo6PCja2oEDB1pOW7z99tu55557InJopT2dLnARLj29wMXi57ew93Q179x33SXPNTZ5mfTAW8wpGMB/fTU/pDUpm7HRAhe7d++OquGE3qit79HlLnBhC7tOnmPMgNQ2n4t1Opg+IosNe0svOQiilFLRolcE+rn6Ro5W1DFmYNuBDnD9ldmcrKpn7+nLG3dTSimr9IpA33vKF9KjB6S022bmqL4AvLdPF95QSkWnXhHoh8p8M/oOz05ut03/NBfDspL4+GCbiy0ppVTE6xWBfqS8lhiHMCi97aPkza4elsEnhyto8uo4ulIq+vSKQD9cXkdOnwRiOrm0f+qwTKrrPew+GdpLoJVSqif0ikA/Ul7L0MykTttdnZcJwKaD5eEuSSlF16fQbbZhwwY+/PDDLrdbsWIFzz33XMjaRxrbXxZpjOFIWR1XDenTadv+aS5yMxPZdLCCb80Y1ml7pdTl6e4Uuhs2bCA5OZlp06Z1qd3ixYtD2j7S2D7QK2rdVDd4guqhg6+X/vrOU3i9BodDl+dSymrLly9nxYoVxMTEMGbMGB566CFWrFiB0+nkN7/5DT//+c+prKzkwQcfxO12k5mZyQsvvMD58+cvabd+/XqSk5O57777gtpvYPvi4mIWL15MaWkpTqeT3/3udwwfPtzqL89FbB/oJWfPA5DTp+MDos2uyu3D6qJjHCyrZUTf9s+KUcpOvve973V7sYn2TJgwgccff7zDNsFMofvQQw9x6NAh4uPjqaysJD09ncWLF7cELcDZs2fZtGkTIsLTTz/NI488wn//939f0m79+vVd2m9g+zvuuINly5Zxyy23UF9fj9frvbwvUBgEFegiMgv4GeAEnjbGPNRGm+uAx4FYoMwYMzOEdXbbySrfqtkDOznDpdmkIb7pAz49elYDXakwC2bIJT8/nzvuuIObb76Zm2++uc02JSUlzJ8/n5MnT+J2u8nLy+v0vYPZb7Pq6mqOHz/OLbfcAoDL5eqwvVU6DXQRcQJPAjcCJcBmEVljjNkV0CYd+B9gljHmqIj0DVfBXXWqytdDH5AW3DdgWFYyqa4Yth6t5PbCweEsTamI0VlPuicFTqG7du1a/vznP/Pee++xZs0afvzjH7Nz585LXnPPPfewdOlS5s6dy4YNG/jRj37U6fsEs99m0TIlSDBnuUwBio0xB40xbmAVMK9Vm68BvzfGHAUwxkTMahEnz9UT53SQkRQXVHuHQ5gwpA9bj54Nc2VKqbY888wzbNu2jbVr1+L1ejl27BjXX389jzzyCJWVldTU1JCSknLR9LhVVVUMGjQI4KKFnVu3axbsfpulpqaSk5PDa6+9BkBDQwN1dXWh/uiXLZhAHwQcC3hc4t8WaBTQR0Q2iMgWEYmY5bBPVtbTP80V1ET5zSYNSWfv6WpqGjxhrEwp1Zmmpia+/vWvM378eCZOnMj3v/990tPT+cpXvsKrr77KhAkTeP/99/nRj37EbbfdxowZM8jKurDmQet2Xd1voOeff57ly5eTn5/PtGnTOHXqVI99HYIVzBh6W0nY+u+PGOAq4AYgAfhIRDYZY/ZdtCORu4G7AYYMGdL1arvhVJUv0Lti0pA+GAPbj1UyvdWCGEqp0Kmpqenw+djY2JY1PAONGjWKHTt2XLRt3rzWAweXtpsxY0bL/WD2G9h+5MiRvP322x3Wa7VgeuglQOBgcg5woo02rxtjao0xZcB7QEHrHRljnjLGFBpjCrOzs7tbc5ecPHeegV0M9OZFMD49osMuSqnoEUygbwZGikieiMQBC4A1rdr8AZghIjEikghcDewObald5/UaTlc10D8tuDNcmqUlxDKybzJbj126UIZSSkWqTodcjDEeEVkCvIHvtMVfG2N2ishi//MrjDG7ReR1YAfgxXdq42fhLDwYZ+vcuJu89EuN7/JrJw5J561dpzHGdGn8XSmlrBLUeejGmLXA2lbbVrR6/BPgJ6Er7fKV17oByE7pTqD34aWiEo6U15GbFdxVpkopZSVbT85VVtMAQGZS1wM9PycNgB3Hq0Jak1JKhYvNA93XQ89KDu4c9ECj+qUQH+Ngh46jK6WihK0Dvby5h57c9R56rNPB2IGp7CjRHrpSVlmzZg0PPeSbaeS1115j166WC9T54Q9/yLp167q8z8OHDzNu3LjLri1U+wklW0/OVV7jxukQ0hNiu/X6/Jx0Xio6RpPX4NSZF5XqUR6Ph7lz5zJ37lzAF+hz5sxhzJgxADzwwANWlheRbN1DL6tpICMprtvT4ObnpFHnbuJAaccXPyiluu65554jPz+fgoIC7rzzTgAWLVrE0qVLuf766/nBD37AypUrWbJkCR9++CFr1qzh/vvvZ8KECRw4cIBFixbx8ssvA7B582amTZtGQUEBU6ZMobq6msOHDzNjxgwmTZrEpEmTOl0QY/78+axde+Hcj0WLFvHKK68EtZ/mOpvNmTOHDRs2APDmm29yzTXXMGnSJG677baWi6mWLVvGmDFjyM/Pb5nd8XLZuodeVuMmM8g5XNqSn+O7wGj7sUpG9UsJVVlKRZ4t34OzoZ0+lz4T4Kq2J/3auXMn//Ef/8HGjRvJysqiouLC4uz79u1j3bp1OJ1OVq5cCcC0adOYO3cuc+bM4dZbb71oX263m/nz57N69WomT57MuXPnSEhIoG/fvrz11lu4XC7279/PwoULKSoqarfcBQsWsHr1ar785S/jdrtZv349v/jFLzDGdGk/gcrKynjwwQdZt24dSUlJPPzww/z0pz9lyZIlvPrqq+zZswcRobIyNMfqbB3o5bUN3TplsdmwrCSS42PYUVLFbTrzolIh8/bbb3Prrbe2zLuSkZHR8txtt92G0+kMel979+5lwIABTJ48GfBNpAVQW1vLkiVL2LZtG06nk3379nW0G2bPns29995LQ0MDr7/+Op/73OdISEigqqqqS/sJtGnTJnbt2sX06dMB3y+fa665htTUVFwuF9/61re46aabmDNnTtD77IitA72spoGhGYndfr3DIYwblMqOEj3TRdlcOz3pcOnogr2kpK5d99Hevh577DH69evH9u3b8Xq9nc5h7nK5uO6663jjjTdYvXo1CxcuDHo/MTExFy14UV9f31LbjTfeyIsvvnjJaz755BPWr1/PqlWreOKJJ0IyT4ytx9DLa9zdOsMlUEFOOrtPVuP2RN7qJEpFqxtuuIGXXnqJ8nLfguyBQy7taW9q2yuvvJITJ06wefNmwLcYhcfjoaqqigEDBuBwOHj++edpamrq9D0WLFjAM888w/vvv8+XvvQlgKD2k5uby7Zt21qm5f3kk08AmDp1Khs3bqS4uBiAuro69u3bR01NDVVVVXz5y1/m8ccfD9lqUbYN9Dq3hzp3E5ndOAc9UH5OOu4mL3tPXfqDpJTqnrFjx/LP//zPzJw5k4KCApYuXdrpaxYsWMBPfvITJk6cyIEDB1q2x8XFsXr1au655x4KCgq48cYbqa+v59vf/jbPPvssU6dOZd++fUH1/L/4xS/y3nvv8YUvfIG4OF92BLOf6dOnk5eXx/jx47nvvvuYNGkSANnZ2axcuZKFCxeSn5/P1KlT2bNnD9XV1cyZM4f8/HxmzpzJY489FuyXrkNi1UochYWFJtgDC8G6rsB3EHPD9kqOVdQx45F3eOTW/Mtaeah5Pw/ePI6vTx0aqlJVtFp3ne/2CxusrCIkdu/ezejRo60uQ3Wgre+RiGwxxhS21d62PfTmy/67c5VooJw+CfRJjNVxdKVUxLNtoFf4J+bK6MY8LoFEhPycdL1iVCkV8Wwb6JV1jQDdvko0UEFOGvtOV3Pe3flBFaWiSbQsftwbded7Y99AP+8P9MTLD/TxOel4Dew8ob10ZR8ul4vy8nIN9QhkjKG8vLzTUy1bs+156FXnGxGBFFdoeugA20uqKMzN6KS1UtEhJyeHkpISSktLrS5FtcHlcpGTk9Ol19g30OvcpLpiQzKpVt9UF/1TXXpgVNlKbGwseXl5VpehQsjWQy5pIRg/b5afk8Zf9MCoUiqC2TfQ6xpDMn7eLD8njYNltVT5x+aVUirS2DbQq0LcQy8Y7LtoSXvpSqlIpYEepPxB/ql0dRxdKRWhggp0EZklIntFpFhElrXx/HUiUiUi2/z/fhj6Urumss4d0iGXtMRY8rKS2K5rjCqlIlSnZ7mIiBN4ErgRKAE2i8gaY8yuVk3fN8aEZlLfy+T1GqrON5KecHmX/bdWkJPGRwfLQ7pPpZQKlWB66FOAYmPMQWOMG1gFzAtvWZenxu3BawjpkAv4xtFPn2vgVFV9SPerlFKhEEygDwKOBTwu8W9r7RoR2S4i/yciY9vakYjcLSJFIlIUzosZqvyX/aeFcMgFApak03F0pVQECibQ27oyp/W1wp8CQ40xBcDPgdfa2pEx5iljTKExpjA7O7trlXZBKOdxCTR2YCoxDtELjJRSESmYQC8BAicUzwFOBDYwxpwzxtT4768FYkUkK2RVdlFVyzwuoR1Dd8U6uaJ/CtuP6amLSqnIE0ygbwZGikieiMQBC4A1gQ1EpL/4F/UTkSn+/Vp29LDyvG/q3FCPoYNvHH17SSVer05opJSKLJ0GujHGAywB3gB2Ay8ZY3aKyGIRWexvdivwmYhsB5YDC4yFU7i1DLmEeAwdYEJOOtX1Hg6X14Z830opdTmCmpzLP4yyttW2FQH3nwCeCG1p3dc85BKOHnr+4OaZFysZlp0c8v0rpVR32fJK0XP1jcTFOHDFOkO+75F9U0iMc+o4ulIq4tgy0KvrPaTEh2dmYKdDGDcwTU9dVEpFHPsGuit8U70XDE5j54lzuD3esL2HUkp1lS0Dvaa+keSwBno6bo+Xfaerw/YeSinVVbYMdN+QS+gPiDYr8F8xuk0n6lJKRRBbBnpNgyesPfScPglkJMXpFaNKqYhiy0AP9xi6iJCfk6ZnuiilIopNA70xbGe5NCvISWf/mWpqGzxhfR+llAqWLQO9psFDiit8Y+gAEwan4zXwl+PaS1dKRQbbBboBvIawDrnAhTVGtx7VcXSlVGSwXaA3T+wbzoOiABlJcQzLSmLLkbNhfR+llAqW7QLd+BM93EMuAJOG9uHTo2excB4ypZRqYb9A92druA+KAlw1tA8VtW6OlNeF/b2UUqoz9gt0/224x9DBF+iADrsopSKC/QK9h8bQAUZkJ5PiimHLUQ10pZT1bBfo9OAYusMhTBzSh0+1h66UigC2C/SWHnoPjKEDXDWkD3tPV1Nd39gj76eUUu2xX6D7b3ss0If2wRidqEspZT37BbqBpDgnTof0yPsVDE7DIXpgVCllPfsFOqZHxs+bpbhiGdUvRQNdKWW5oAJdRGaJyF4RKRaRZR20mywiTSJya+hK7KIeuOy/tauG9mHb0Uq8Xr3ASCllnU4DXUScwJPAbGAMsFBExrTT7mHgjVAX2RWGnjllMdBVQ/tQ3eBh/5maHn1fpZQKFEwPfQpQbIw5aIxxA6uAeW20uwd4BTgTwvq6zJieHXKBCxcYbT5c0aPvq5RSgYIJ9EHAsYDHJf5tLURkEHALsKKjHYnI3SJSJCJFpaWlXa01KIaeuew/0JCMRPqlxvPxIQ10pZR1ggn0tk4XaT1Y/DjwA2NMU0c7MsY8ZYwpNMYUZmdnB1tjlxgLxtBFhKvzMvnkULlO1KWUskwwgV4CDA54nAOcaNWmEFglIoeBW4H/EZGbQ1JhN/TUOeiBrh6WwelzDTpRl1LKMsEE+mZgpIjkiUgcsABYE9jAGJNnjMk1xuQCLwPfNsa8FvJqg2AwPX5QFODqvAwAPj5U3uPvrZRSEESgG2M8wBJ8Z6/sBl4yxuwUkcUisjjcBXZFT18lGmh4djJZyXF8fFDH0ZVS1ggq+Ywxa4G1rba1eQDUGLPo8svqnubh6yQLAl1EmJKXoQdGlVKWsdmVor5EtyLQAabkZnC88jwlZ3UcXSnV82wV6C099DinJe9/9bBMAB12UUpZwl6B7r9NjLOmh35FvxTSEmL5RIddlFIWsFWg0zKGbk0P3eEQJudm6JkuSilL2CrQjcVj6ABTh2VwuLyOk1XnLatBKdU72SzQfZIsGnIBmD4iC4AP9pdZVoNSqneyVaA3J3qiRUMu4BtHz0qOY2OxBrpSqmfZKtAjoYfucAjTR2TxQbHO66KU6ln2CnQDgvTY8nPtuXZEFmU1Dew5VW1pHUqp3sVWgX7pJJDWuHakbxxdh12UUj3JVoFuALG2cw7AgLQEhmcn8b4eGFVK9SB7Bbppe/J2K8wYmc3Hh8pp8HQ4RbxSSoWMrQIdiIwuOr7TF+sbvXx6pNLqUpRSvYStAt0YEzE99KnDMnA6hA+Kw7PUnlJKtWavQCdyhlxSXLFMHJyu4+hKqR5jq0AHIifRgZmjstlRUkVpdYPVpSilegFbBXokHRQF+PzovgBs2HvG4kqUUr2BvQId38pBkWLMgFT6p7p4RwNdKdUDbBPovtMDI+PComYiwvVXZvPevjLcHq/V5SilbM42gV7X4DvfO3L65z7XX9GXmgYPRYd10QulVHgFFegiMktE9opIsYgsa+P5eSKyQ0S2iUiRiFwb+lI7Vuv2+Gvp6Xfu2PQRWcQ5Hby9R4ddlFLh1Wmgi4gTeBKYDYwBForImFbN1gMFxpgJwF8DT4e60M7UuZuvyIysRE+Kj2Hq8Eze1nF0pVSYBdNDnwIUG2MOGmPcwCpgXmADY0yNuTBXbBIWDGbXNERmDx3g81dkc7C0lkNltVaXopSysWACfRBwLOBxiX/bRUTkFhHZA/wZXy/9EiJyt39Ipqi0NLRXUEbqGDrADaP7AfDWrlMWV6KUsrNgAr2tjLykB26MedUYcyVwM/DjtnZkjHnKGFNojCnMzs7uWqWdaB5Dj0SDMxIZNyiV1z/TQFdKhU8wgV4CDA54nAOcaK+xMeY9YLiIZF1mbV1SG8FDLgCzxvbn06OVnKqqt7oUpZRNBRPom4GRIpInInHAAmBNYAMRGSH+K3pEZBIQB5SHutiO1Lqbh1wiM9FnjesPwJs67KKUCpNOA90Y4wGWAG8Au4GXjDE7RWSxiCz2N/sr4DMR2YbvjJj5pocX1Kzz99AjNM8Z0TeFEX2TddhFKRU2Qa2mbIxZC6xttW1FwP2HgYdDW1rXtAy5WFlEJ2aN7c8v3j1ARa2bjKQ4q8tRStmMba4UrXU3RexwS7NZ4/rT5DWs23Xa6lKUUjZkm0Cvi+CzXJqNHZhKTp8E/vyXk1aXopSyIdsEek1DU8Se4dJMRJiTP5APissor9E50pVSoWWbQK9r8ET4gIvPzRMH0uQ12ktXSoWcbQK91u2J3JPQA1zZP5Ur+qXwh23tnsqvlFLdYp9Ab2iKih46wLyJA9ly5CzHKuqsLkUpZSP2CXS3Jxo66ADMLRgIwB+2Hbe4EqWUndgm0Jsn54oGOX0SmZzbh9e2naCHr79SStmYbQK9NkoOijabN2EQxWdq+Oz4OatLUUrZhC0C3RjjH3KJnkj/SsFA4mMcrC46anUpSimbsEWgN3i8eKNs5CItIZYvjx/AH7ad4Lw7eoaLlFKRyxaBXhMF87i05fbCwVTXe/i/z/ScdKXU5bNFoLesVhRliT51WAa5mYms2nys88ZKKdUJWwT6hdWKoivRRYTbJw/mk0MVHCytsbocpVSUs0WgN0/MFW09dIBbJ+XgdAirtZeulLpMtgj0mgheILozfVNdfHFMP1ZtPqYHR5VSl8UWgR7pqxV1ZtG0XKrON+qVo0qpy2KLQI/09UQ7MyUvg9EDUln54WG9clQp1W32CPQoPW2xmYiwaNpQ9pyqZtPBCqvLUUpFKXsEuju6h1zANxVAemIsz3542OpSlFJRKqhAF5FZIrJXRIpFZFkbz98hIjv8/z4UkYLQl9q+2gYPDonqPMcV62ThlCG8uesUh8tqrS5HKRWFOg10EXECTwKzgTHAQhEZ06rZIWCmMSYf+DHwVKgL7UhtQxNJ8TE9+ZZhcde0XGKcDn753kGrS1FKRaFgeuhTgGJjzEFjjBtYBcwLbGCM+dAYc9b/cBOQE9oyO1bb4CEpLvoDvW+qi1uvyuGVLSWcPldvdTlKqSgTTKAPAgKveinxb2vP3wD/dzlFdVWt20NSvLMn3zJs/u5zw/B4vfzqg0NWl6KUijLBBHpbQ9NtnlsnItfjC/QftPP83SJSJCJFpaWlwVfZidqGJpJtMOQCMDQziTn5A3lh0xGq6hqtLkcpFUWCCfQSYHDA4xzgkhWORSQfeBqYZ4wpb2tHxpinjDGFxpjC7Ozs7tTbptoGD4k2GHJp9vfXDafW3cSvNmovXSkVvGACfTMwUkTyRCQOWACsCWwgIkOA3wN3GmP2hb7MjtW67XFQtNnoAanMHtefX71/kIpat9XlKKWiRKeBbozxAEuAN4DdwEvGmJ0islhEFvub/RDIBP5HRLaJSFHYKm5DbYN9xtCbLb1xFOcbm1jx7gGrS1FKRYmgurXGmLXA2lbbVgTc/xbwrdCWFrw6t8dWPXSAkf1SuHniIJ798DB/c20e/VJdVpeklIpwtrhStKbBQ1KcvXroAN//wii8xrB8/X6rS1FKRYGoD3RPk5f6Rq/teugAgzMSWThlCKs2H6P4TLXV5SilIlzUB3pdo2+mRbucttjad28YSWKckwf+tFtnYlRKdSjqA715pkU7nbYYKDM5nu/eMJL39pXyzt4zVpejlIpgNgh0Xw/dbme5BPrGNbkMy07iwT/txu3xWl2OUipC2SDQfT10O8zl0p64GAf/etMYDpbV8vQHOnGXUqpt9gl0m46hN7v+yr7MGtufn63br9PrKqXaFP2B7rb3QdFA/z5vLHExDv7p93/RA6RKqUtEf6A3HyOdMsIAAA8mSURBVBS18Rh6s36pLv5p9mg+OljO74pKrC5HKRVhoj/Q/cvP9YYeOsCCyYOZkpvBg3/excmq81aXo5SKINEf6C2nLdq/hw7gcAgP35pPY5PhH17ajterQy9KKR8bBLpvDN2u56G3JS8riR/NHcOHB8r1rBelVAsbBLqHhFgnTkc0LxHddbcXDmbW2P785I29fHa8yupylFIRIPoD3YYzLQZDRPivr44nIymOe17cyrl6Xd1Iqd4u+gO9oYnkXnCGS1v6JMWxfMFEjlXUsXS1jqcr1dvZINDttfxcV109LJN/uWk063af5ol3iq0uRylloegPdLen15yy2J5vTsvlqxMH8di6fazffdrqcpRSFon+QG9o6hUXFXVERPjPr45n3MA0lvx2KztKKq0uSSllgegP9F56ULQ1V6yTXy0qJDM5jr9euZljFXVWl6SU6mHRH+g2XX6uO/qmuFh51xQamwzffOYTKmrdVpeklOpBQQW6iMwSkb0iUiwiy9p4/koR+UhEGkTkvtCX2b7ahibtoQcY0TeZp79ZyPGz57nj6Y+prNNQV6q36DTQRcQJPAnMBsYAC0VkTKtmFcC9wKMhr7ADXq/Rg6JtmJybwf9+o5ADpTV8/VcfU1Wn56gr1RsE00OfAhQbYw4aY9zAKmBeYANjzBljzGagR5Oj1u3BGEhxaaC39rlR2fzyzqvYd6qGO3/9MWd1+EUp2wsm0AcBxwIel/i3Wa663jcxV4or1uJKItP1V/TlF1+fxJ5T1dz2y484UamzMyplZ8EEeluTpHTrkkQRuVtEikSkqLS0tDu7uEjz5e6pGujtumF0P5776ymcrqrnr37xIcVnqq0uSSkVJsEEegkwOOBxDnCiO29mjHnKGFNojCnMzs7uzi4ucqGHrkMuHZk6LJNVfzeVxibDrSs+4sPiMqtLUkqFQTCBvhkYKSJ5IhIHLADWhLes4Jw77+uha6B3buzANH7/99PISo7nzl9/wsqNh3QZO6VsptNAN8Z4gCXAG8Bu4CVjzE4RWSwiiwFEpL+IlABLgX8RkRIRSQ1n4XChh56aoEMuwRiSmcir357G9Vf05Ud/3MU/vryD+sYmq8tSSoVIUF1bY8xaYG2rbSsC7p/CNxTTo6rrtYfeVSmuWJ668yoeX7eP5W8Xs72kkuULJ3Jl/7D//lVKhVlUXyl6rrmHrgdFu8ThEJZ+8Qqe/espVNQ2MveJjTz30WEdglEqykV5oDcS53QQHxPVH8MyM0dl8/r3ZjBteCY//MNOvvHrT3QOGKWiWFQnYXW9hxRXDCK9a/m5UMpKjueZRZP58byxfHrkLF987D2efv8gniav1aUppbooqgP93PlGHT8PARHhzmtyeWvpTKYNz+TBP+9m3pMb+ehAudWlKaW6IKoDvbreo2e4hNDA9ASe/mYhT35tEmdr3Sz8303c/VwRh8pqrS5NKRWEKA907aGHmohwU/4A3r7vOu7/0hVsLC7ji4+9y7+89hedOkCpCBfVgX6u3kNKvPbQw8EV6+Q714/gnfuv47bCwazefIyZP3mHf371LxzXYFcqIkV1oGsPPfz6prj4z1vGs+H+65k/eTC/Kyph5iPvcO+LW9l69KzV5SmlAkRtGhpjOFvXSEZSnNWl9AqD0hN48ObxfPu6Efzqg0O8tPkYa7afYMLgdO6ansuXxvbHFasrRyllpagN9Fp3E26PVwO9hw1MT+Bf54zh+zeO4pUtJTyz8RDfXbWNtIRY5hYM5NarcsjPSdNTSZWyQNQGevOCDX000C2RHB/DN6flcufUoWw8UMbvikp4qegYz286wqh+ycwtGMiscQMY0TfZ6lKV6jWiNtCbF0DOSNRAt5LDIcwYmc2MkdlUnW/kTztO8PtPj/Pom/t49M19jOqXzOxxA/jS2P6MHpCiPXelwih6A71Oe+iRJi0hljuuHsodVw/lZNV53vjsFGs/O8Xyt/fzs/X76Zcaz4yR2cwclc21I7L0e6dUiEVtoDcPuegYemQakJbAoul5LJqex5nqejbsKeXd/aW8tes0L28pQQTyc9K5Oi+DwqF9mJyboQGv1GWK2kCv0ECPGn1TXNw+eTC3Tx5Mk9ewvaSSd/eWsrG4jJUbD/PUewcBGNk3mcLcDCYOSWf8oDRG9E0m1hnVZ9Yq1aOiNtDLatzEOoVUPQ89qjgdwqQhfZg0pA/fv3EU9Y1N7CipYvPhCjYfruBP20/w4idHAYiPcXDlgFTGD0pl3MA0Rg9IZUTfZJLi9XuuVFui9n/Gyarz9E9z6UG2KOeKdTIlL4MpeRkAeL2GQ+W1fHa8is+OV/GX41X8YesJfrPpaMtrBqUnMLJfMiP7JjOybwrD+yaTl5VEn8RY/XlQvVr0BnplPQPSEqwuQ4WYwyEMz05meHYy8yYMAnwhf6Sijr2nqik+U83+MzXsP13DRwfKafBcmOY3OT6GwRmJDM1IZEhmIoMzEhmSkcjgPgn0T3ORGBe1P+5KBSVqf8KPV55v6dUpe3M4hLysJPKykoD+LdubvIaSs3XsP13DkYo6jlXUcbSijuLSGt7eewa35+I53VNdMQxIS6BfmosBqS7fbZqL/qkuslPiyUiKIyMpTq94VVErKgO9yWs4fa6eAWkuq0tRFnI6hKGZSQzNTLrkOa/XcKa6gaMVdZScrePUuXpOVfn/natnz8lzlNY00Naqe8nxMWQm+8I9MymOzKR4MpJ99+dVN+B0CPsPlpOaEEtaQiypCbEkxTl1uEdZLqhAF5FZwM8AJ/C0MeahVs+L//kvA3XAImPMpyGutcWZ6no8XsOAdB1yUW1zOIT+aS76p7na/UuusclLaXUDJ6vqKatpoKLWTUWtu+V+eY2b45X17CipoqLWjcdrGDesBoAFH226aF9Oh+8AfXPAp7piSU2IISkuhqT4GBLjnBdu42JIjPffBm733ybGxeB06C8H1XWdBrqIOIEngRuBEmCziKwxxuwKaDYbGOn/dzXwC/9tWOw9VQ34TnNTqrtinQ4GpicwMIiOgTGGc+c9xL/7MB6vlxc+fzVV5xs5d76Rc/WN/vuegPuNnDpXT12Dh1p3E7UNHjze4BfhjnUK8TFOXLEO4mOcxMc6cAU8dsU6cMU6iY/x3Tbfj4/1PRfndBAX4yDG4SDWKcTFOIh1OohxCLExvudjnQ5inNJyP9Yp/lv//RgHsf7XOx2if4FEgWB66FOAYmPMQQARWQXMAwIDfR7wnPEtG79JRNJFZIAx5mTIKwY+PXIWh8DoAanh2L1SlxAR0hJjIdYBOJg+IqvL+3B7vNS5PdQ0eKjzh/xFt24PdQ2+2waPl4ZGL/WeJuobm/yPfbf1jU2U1Xho8DRR3+i95DYcRCDG4Qv2GIcDh0CM04HTITjFv9154b6zpa3g8N9e2O77xeIQ/3b/6y5u43veIb6vffN9h0MQoeWxU6TD5x0tz9P2/iSgvSP49iKCEPAY/zb/fYdD/NsutA18TXZKPP3DMGQcTKAPAo4FPC7h0t53W20GASEP9Lf3nOaX7x1kSl4Gabr8nIoicTEO4mLiSA/j/EPGGBo8XhqbvDQ2GRqbvLg9XjzeC/ebn/M0eXEHtGtsdd/X9kI7j9fQ1Oqfb5uXJi80eX1tvMbgafK3Mf52Tb77jY1emrxNLa/1eg0er/dC2ybT8j4G8BpfG2P891tuL9xv6zhIpFs8czjLZl8Z8v0GE+ht/Z3V+ksYTBtE5G7gboAhQ4YE8daX6p+awPVX9OUfZ13RrdcrZWci0jIE01sYc2ngm1bBb1r9Mriovbdr7Zu8BvBtM/hfD77H5sJ9b8v95rb+WwNDMxPD8rUIJtBLgMEBj3OAE91ogzHmKeApgMLCwm79Xh0zMJUVd17VnZcqpWyoeajD0Wa/sncJZqKMzcBIEckTkThgAbCmVZs1wDfEZypQFa7xc6WUUm3rtIdujPGIyBLgDXynLf7aGLNTRBb7n18BrMV3ymIxvtMW7wpfyUoppdoS1Hnoxpi1+EI7cNuKgPsG+E5oS1NKKdUVOjepUkrZhAa6UkrZhAa6UkrZhAa6UkrZhAa6UkrZhBiLrpsVkVLgSDdfngWUhbCcaKCfuXfQz9w7XM5nHmqMyW7rCcsC/XKISJExptDqOnqSfubeQT9z7xCuz6xDLkopZRMa6EopZRPRGuhPWV2ABfQz9w76mXuHsHzmqBxDV0opdalo7aErpZRqJeoCXURmicheESkWkWVW1xNuIjJYRN4Rkd0islNEvmt1TT1FRJwislVE/mR1LT3Bv3TjyyKyx//9vsbqmsJNRL7v/7n+TEReFJHQr8tmMRH5tYicEZHPArZliMhbIrLff9snFO8VVYEesGD1bGAMsFBExlhbVdh5gH8wxowGpgLf6QWfudl3gd1WF9GDfga8boy5EijA5p9dRAYB9wKFxphx+KbnXmBtVWGxEpjVatsyYL0xZiSw3v/4skVVoBOwYLUxxg00L1htW8aYk8aYT/33q/H9Jx9kbVXhJyI5wE3A01bX0hNEJBX4HPArAGOM2xhTaW1VPSIGSBCRGCCRNlY6i3bGmPeAilab5wHP+u8/C9wciveKtkBvbzHqXkFEcoGJwMfWVtIjHgf+EQjPMvaRZxhQCjzjH2Z6WkSSrC4qnIwxx4FHgaP4FpSvMsa8aW1VPaZf86pu/tu+odhptAV6UItR25GIJAOvAN8zxpyzup5wEpE5wBljzBara+lBMcAk4BfGmIlALSH6MzxS+ceN5wF5wEAgSUS+bm1V0S3aAj2oxajtRkRi8YX5C8aY31tdTw+YDswVkcP4htU+LyK/sbaksCsBSowxzX99vYwv4O3sC8AhY0ypMaYR+D0wzeKaesppERkA4L89E4qdRlugB7Ngta2IiOAbV91tjPmp1fX0BGPMPxljcowxufi+x28bY2zdczPGnAKOicgV/k03ALssLKknHAWmikii/+f8Bmx+IDjAGuCb/vvfBP4Qip0GtaZopGhvwWqLywq36cCdwF9EZJt/2//zr/Oq7OUe4AV/Z+UgNl9s3RjzsYi8DHyK72yurdjwqlEReRG4DsgSkRLg34CHgJdE5G/w/WK7LSTvpVeKKqWUPUTbkItSSql2aKArpZRNaKArpZRNaKArpZRNaKArpZRNaKArpZRNaKArpZRNaKArpZRN/H97eRbkS24drgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "xrange = np.arange(0,10,0.01)\n", "fdist = f.pdf(xrange,dof1_lp,dof2_lp)\n", "lowerbound = f.ppf(alpha,dof1_lp,dof2_lp)\n", "upperbound = f.ppf(1-alpha,dof1_lp,dof2_lp)\n", "\n", "plt.plot(xrange, fdist, label = 'F-dist PDF')\n", "plt.axvline(fstat_lp, label = 'F-statistic', color='black', zorder=10)\n", "plt.axvline(lowerbound, label = 'critical values', color = 'orange')\n", "plt.axvline(upperbound, color = 'orange')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Chi-square distribution\n", "To determine the bounds of the true population variance $\\sigma$, the chi-square distribution can be used. **Again, the samples must be from an underlying normal distribution, otherwise, the test is not valid, and some sort of non-parametric test should be used instead.**
\n", "\n", "95% confidence bounds on $\\sigma^2$:\n", "\\begin{equation}\n", "s^2\\frac{N-1}{\\chi^2_{0.975}} \\le \\sigma^2 \\le s^2\\frac{N-1}{\\chi^2_{0.025}}\n", "\\end{equation}\n", "\n", "where $\\chi^2$ is the chi-square distribution, which depends on degrees of freedom $\\nu$." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "variance of x1: 0.27176954758552535\n", "95% confidence bounds on variance of x1: 0.12159093658237595 1.0531255173713834\n" ] } ], "source": [ "#applying to the low pass filtered version of x1\n", "#chi2 2.5 and 97.5 percdentiles\n", "chi2_lower = chi2.ppf(1-alpha, dof1_lp)\n", "chi2_upper = chi2.ppf(alpha, dof1_lp)\n", "\n", "#confidence bounds on variance\n", "var1_lower = var1_lp * dof1_lp / chi2_lower\n", "var1_upper = var1_lp * dof1_lp / chi2_upper\n", "\n", "print('variance of x1:', var1_lp)\n", "print('95% confidence bounds on variance of x1:', var1_lower, var1_upper)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU1fn48c+TyR4SlrBmDxgNS0jYwiouKIKoiLRVrAvKV8SqtVqt1PanttXaVmuV1g2Kgq0VN1S0uKFYQVEWWUMgbAFCIAmE7AlJZs7vjxliCFmGZMIsed4v55WZe8+581wmPjlz7rnniDEGpZRSvsvP3QEopZRqX5rolVLKx2miV0opH6eJXimlfJwmeqWU8nH+7g6gMd27dzcJCQnuDgNKdtp/Rpzn3jiU79HfLeViGzZsOGqM6dHYPqcSvYhMAp4FLMA/jTF/arA/GXgFGAr8xhjzlGN7LPAq0BuwAfONMc+29H4JCQmsX7/emdDa14oL7T8v+dKdUShfpL9bysVEZH9T+1pM9CJiAZ4DLgVygHUisswYs71esULg58DVDarXAr80xnwvIuHABhH5rEFdpZRS7ciZPvp0YLcxZq8xphpYAkytX8AYk2+MWQfUNNh+2BjzveN5KZAJRLskcqWUUk5xJtFHAwfrvc6hFclaRBKAIcB3TeyfLSLrRWR9QUHBmR5eKaVUE5zpo5dGtp3RvAki0gl4B/iFMaaksTLGmPnAfIDhw4frvAzK69TU1JCTk0NVVVXLhSMfsf/MzGzfoJTPCQ4OJiYmhoCAAKfrOJPoc4DYeq9jgFxn30BEArAn+deMMUudjkwpL5OTk0N4eDgJCQmINNY+qqfE8WVaR92oM2CM4dixY+Tk5JCYmOh0PWe6btYBSSKSKCKBwHXAMmcOLvbf9oVApjHmaaejUsoLVVVVERkZ2XKSV6qVRITIyEjnvjXW02KL3hhTKyJ3AZ9gH175sjEmQ0TmOPa/KCK9gfVABGATkV8AA4DBwI3AVhHZ5DjkQ8aY5WcUpVJeQpO8am+t+R1zahy9IzEvb7DtxXrPj2Dv0mloNY338aszUHailtW7CthTUA5Avx5hjDmnOxHBzvfRKaU6Lo+8M1bZVdVY+fsXu3jl62wqqq2n7Avy9+PaEbHcMyGJyE5BbopQKeUNNNF7qIOFFcxavI6svDKuTI3ihpFxDI7pgsGQkVvCOxtyeO27A/x3y2GevjaNC85t9M5npZTSSc080f5j5Ux/4RuOFFex+NZ0/j5jCCP7RhISaCE00J8RCd340/TBLP/5+fQID+LWRet4Z0OOu8NWHuq9997jtttuY+rUqXz66afuDke5gSZ6D1NUUc1NL6+l2mrjrTljmm2pn9c7nLfmjGZU32788q3NLP4m++wFqrzG1VdfzYIFC1i0aBFvvPFGi+UfffRRnnrqqbrXY8aMabJsUVERzz//fKtjO3nshsfJzs5m0KBBLda3WCykpaUxcOBAUlNTefrpp7HZbHX7c3JymDp1KklJSfTr14977rmH6urq0+qffGRnZ7f4njt37jylTkREBM8888xp5RISEkhJSSEtLY3hw4fXbT948CAXXXQR/fv3Z+DAgTz7bIvTf7WdMcbjHsOGDTMe4bML7I+zxGazmVmL1ppzHvqvWZ99zOl6VTW15rbF60zC3A/Nx9sOt2OEqjnbt293vnDxDvvjLLrvvvvMhg0bWiz3yCOPmCeffNKpY+7bt88MHDiwraGddhxnjxsWFlb3PC8vz0yYMME8/PDDxhj7/08jRowwL7/8sjHGmNraWnPrrbea+++/v9H6rVFbW2t69eplsrOzT9sXHx9vCgoKTtuem5tb9zmUlJSYpKQkk5GRcUbv29jvGrDeNJFTtUXvQd5Yd5AVmfk8dHl/hsV3c7pekL+FeTOGMDimC79YsoktOUXtGKXyVOXl5URHR/PQQw8BsG7dOtLS0qisrOTBBx9k8uTJDB06tNG6jz/+OOeddx6XXHIJO3fuPGVfp06dKC8vZ8qUKaSmpjJo0KC6bwZz585lz549pKWl8cADD5xS7y9/+Qvz5s0D4N577+Xiiy8G4PPPP+eGG26oO3ZTx7Fardx2220MHDiQiRMnUllZ2ez59+zZk/nz5/OPf/wDYwxffPEFwcHB3HLLLYC99f63v/2Nl19+mYqKCuf+UVvw+eef069fP+Lj452u06dPn7rPITw8nP79+3Po0CGXxNMUTfQeoqD0BI//N5PRfSO5eXTCGdcPDrDwz5uG0y0skDv+/T0lVTUtV1I+JSwsjC1btvD6669TUVHBrbfeyqJFi1iwYAErVqzg7bff5sUXXzyt3oYNG1iyZAkbN25k6dKlrFu37rQyH3/8MVFRUWzevJlt27YxadIkAP70pz/Rr18/Nm3axJNPPnlKnfHjx7Nq1SoA1q9fT1lZGTU1NaxevZrzzz//lLKNHWfXrl3ceeedZGRk0KVLF955550W/w369u2LzWYjPz+fjIwMhg0bdsr+iIgI4uLi2L17NwCVlZV1XTDTpk0D4Pzzzz+la+bkY8WKFae935IlS5gxY0ajsYgIEydOZNiwYcyfP7/RMtnZ2WzcuJGRI0e2eG5toaNuPMQ/vthFRY2Vx6YNws+vdbce9AgP4u/XD+HHL67h4fe28cx1Q1wcpXLW7z7IYHtuo9M62VkdLUpLodPHHBAVwSNXDmy2TGRkJKGhocyaNYsbb7yxLkn9/Oc/b7LOqlWrmDZtGqGhoQBcddVVp5VJSUnh/vvv58EHH+SKK644LVE3ZtiwYWzYsIHS0lKCgoIYOnQo69evZ9WqVXUt/eYkJiaSlpZWdyxn+s/B3h198mdjNxfV3x4SEsKmTZtO2X/yj1NLqqurWbZsGU888USj+7/++muioqLIz8/n0ksvJTk5mfHjx9ftLysrY/r06TzzzDNEREQ49Z6tpS16D7D/WDmvfXeA60bE0q9HpzYda2hcV+6ZkMR7m3J5b2P7fh1Unmnw4MEcPnyY+++/3+k6Ld1tee6557JhwwZSUlL49a9/ze9///sWjxkQEEBCQgKvvPIKY8aM4fzzz2flypXs2bOH/v37t1g/KOiH+0MsFgu1tbUt1tm7dy8Wi4WePXsycODA0xYwKikp4eDBg/Tr16/JYzjbov/oo48YOnQovXr1avQ4UVFRgL1Ladq0aaxdu7ZuX01NDdOnT+enP/0p11xzTYvn1VbaovcAf/00iwCLH/dMSHLJ8X52YT/+l1XAox9kcH5Sd72hyg1aanm311KCBQUFrFy5kt/+9rf4+TnXjhs/fjwzZ85k7ty51NbW8sEHH3D77befUiY3N5du3bpxww030KlTJxYtWgTY+5hLS0ubPfZTTz3Fyy+/TEpKCvfddx/Dhg077Q9LS8dxRkFBAXPmzOGuu+5CRJgwYQJz587l1Vdf5aabbsJqtfLLX/6SmTNn1n17aYyzLfrXX3+9yW6b8vJybDYb4eHhlJeX8+mnn/Lwww8D9m8Us2bNon///tx3331nfqKtoC16N8vKK2XZ5lxuHZdAz4hglxzT3+LHE9ekUFZVyxMf7XDJMZV3mDVrFhdffDGbN292us7QoUO59tprSUtLY/r06Y12y2zdupX09HTS0tJ4/PHH+e1vfwvYu4rGjh3LoEGDTrsYC/bW8eHDhxk9ejS9evUiODi40eO3dJymnOxjHzhwIJdccgkTJ07kkUfsU0CLCO+++y5vvfUWSUlJnHvuuQQHB/PHP/7R6eM3paKigs8+++y01vjll19Obm4ueXl5jBs3jtTUVNLT05kyZUrddY2vv/6af/3rX3zxxRd13xaWL2/f6b/kZH+WJxk+fLjpKGvG/nrpVpZ+n8O3v55A17BAlx77Lx/v4Pkv97Bk9ihG9Y106bHV6TIzM53qkgDapUX/0ksv8cUXX/DUU08xefJktm3b5rJjK8/S2O+aiGwwxgxvrLy26N2oqKKadzfmMG1ItMuTPMDdFycR2y2Eh9/fRq3V1nIF5bV27drF008/zYsvvkhsbCx9+vRh4sSJ7g5LeQhN9G60ZN1BqmpszByb0C7HDwm08JvL+5OVV8ZbOkWCT0tKSmLnzp107doVgM8++0ynO1B1NNG7Sa3VxqvfZDO6byTJvdtvaNVlA3szPL4rf/00i/ITLY9aUEr5Hk30bvLVrgJyi6u4eUxCu76PiPCbKf05WnaCl/63p13fSynlmTTRu8m7G3PpGhrAxck92/29hsR1ZcrgPixcvY/C8uqWKyilfIomejcorarh04wjXDE4ikD/s/MR3HtJEhU1Vhas2ntW3k8p5Tk00bvBJxl5nKi1cfWQ6LP2nuf0DOfKwVEs/iZbW/VKdTCa6N3g3Y05xHULZWhcl7P6vj+fcA6V2qpXqsPRRH+WHSmu4ps9x7h6SHSrVnNvC23VK9UxaaI/yz7bfgRj4MrBfdzy/idb9fO/0la9r3n22WcZNGgQAwcOPGXFo0cffZTo6OjTbrf/+uuvGTx4MCNGjKibtreoqIjLLruM9r5j/q233qJ///5cdNFFrF+/vsnZNRMSEjh69Gi7xuIJ7/vll19yxRVXtNvxdVKzs+zT7Xkkdg/jnJ5tm6Wytc7pGc4Vg6P415ps7riwH51DAtwSh3Ktbdu2sWDBAtauXUtgYCCTJk1iypQpJCXZJ8q79957T5vN8q9//SvvvPMO2dnZvPDCC/z1r3/lD3/4Aw899FC7f9tcuHAhzz//PBdddBHAKUvtKdfTFv1ZVFJVw7d7j3HpgF5nvdumvtvH96W82sp/vjvgthiUa2VmZjJq1ChCQ0Px9/fnggsu4N133222TkBAAJWVlVRUVBAQEMCePXs4dOgQF1xwQZN11q1bx5gxY+om6yotLaWqqopbbrmFlJQUhgwZwsqVKwFYtGgR11xzDZMmTSIpKYlf/epXAPz+979n9erVzJkzhwceeOCU1uyxY8eYOHEiQ4YM4fbbbz/lm8W///3vuonVbr/9dqxWK2Bfpeo3v/kNqampjBo1iry8PADy8vKYNm0aqamppKam8s033zR7nIaefPJJ0tPTSU9Pr/vGs3//fiZMmMDgwYOZMGECBw7Y/x+aOXMmb7/9dl3dkytnffnll1x44YX86Ec/Ijk5mZ/+9Kd15/Txxx+TnJzMuHHjWLp0aV3d//3vf3XfvoYMGdLmWT0BXTO2WS5eM3bZpkMm/sEPzdp9zq8H216uX7DGjHjsM1NVU+vuUHzGKet4rr/nh9+fxh4fjbA/mivT8LH+nmbfOykpyRw9etSUl5ebUaNGmbvuussYY18DNj4+3qSkpJhbbrnFFBYWGmOM2bhxoxk5cqS58MILzcGDB821115rsrKymnyPEydOmMTERLN27VpjjDHFxcWmpqbGPPXUU2bmzJnGGGMyMzNNbGysqaysNK+88opJTEw0RUVFprKy0sTFxZkDBw4YY4y54IILzLp164wxxqxcudJMmTLFGGPM3XffbX73u98ZY4z58MMPDWAKCgrM9u3bzRVXXGGqq6uNMcbccccdZvHixcYYYwCzbNkyY4wxDzzwgPnDH/5gjDHmJz/5ifnb3/5mjLGv7VpUVNTsceqLj483jz32mDHGmMWLF9fFd8UVV5hFixYZY4xZuHChmTp1qjHGmJtvvtm89dZbdfVPrkW7cuVKExERYQ4ePGisVqsZNWqUWbVqlamsrDQxMTEmKyvL2Gw28+Mf//iU91i9erUxxpjS0lJTU1PT6OfdEG1dM1ZEJonIThHZLSJzG9mfLCJrROSEiNx/JnU7ks+25xEZFsjQuK7uDoXZ4/uRX3qC9zflujsU5QL9+/fnwQcf5NJLL2XSpEmkpqbi72/vmb3jjjvYs2cPmzZtok+fPvzyl78EIC0tjW+//ZaVK1eyd+9eoqKiMMZw7bXXcsMNN9S1jE/auXMnffr0YcSIEYB9WT5/f39Wr17NjTfeCEBycjLx8fFkZWUBMGHCBDp37kxwcDADBgxg//79zZ7HV199Vbee7JQpU+rm7vn888/ZsGEDI0aMIC0tjc8//5y9e+3XmQIDA+u+EdRfieqLL77gjjvuAOwLl3Tu3LnZ4zR0cq75GTNmsGbNGgDWrFnD9ddfD8CNN97I6tWrm/9ggPT0dGJiYvDz8yMtLY3s7Gx27NhBYmIiSUlJiEjdOQOMHTuW++67j3nz5lFUVFT3ObZFi0cQEQvwHHApkAOsE5Flxpjt9YoVAj8Hrm5F3Q6hxmpj5c58Jg3sjaWVSwW60vik7iT3DmfBV3v50dCYVi9fqJow7Jnm97fDNMWzZs1i1qxZADz00EPExMQAnLIC0m233XbaRT9jDI899hhvvPEGd911F7/73e/Izs5m3rx5PP7446eUa2ppvqa0ZpWopt7j5ptvbnTZvoCAgLo6Lb1Hc8dpLo6mulpPbvf398dms9W9R3X1D6Pamvo3aOqYc+fOZcqUKSxfvpxRo0axYsUKkpOTW4y3Oc606NOB3caYvcaYamAJMLV+AWNMvjFmHdBwReoW63YU67ILKa2q5ZIBjS87draJCLPH92VXfhlfZuW7OxzlAvn59s/xwIEDLF26tK5Fevjw4boy7777LoMGDTql3uLFi+tazxUVFfj5+eHn50dFRcUp5ZKTk8nNza1bPLy0tJTa2lrGjx/Pa6+9BkBWVhYHDhzgvPNa9wes/rE++ugjjh8/Dti/Gbz99tt151hYWNjit4MJEybwwgsvAGC1WikpKTmj47zxxht1P0ePHg3AmDFjWLJkCQCvvfYa48aNA+yjdDZs2ADA+++/T01Nw1R4quTkZPbt28eePfb5p15//fW6fXv27CElJYUHH3yQ4cOHs2NH2xcPcuY7QTRwsN7rHMDZJcudrisis4HZAHFxcU4e3nus2nUUfz9h7Dnd3R1KnStTo3jyk5289L+9XJzsGX+AVOtNnz6dY8eOERAQwHPPPVfX7fGrX/2KTZs2ISIkJCTw0ksv1dWpqKhg8eLFdVMa33fffUyfPp3AwMBTkg/Yu0jeeOMN7r77biorKwkJCWHFihX87Gc/Y86cOaSkpODv78+iRYtOacWeiUceeYQZM2YwdOhQLrjggrpcMGDAAB577DEmTpyIzWarO8f4+Pgmj/Xss88ye/ZsFi5ciMVi4YUXXmD06NFOH+fEiROMHDkSm81W928xb948br31Vp588kl69OjBK6+8Ati/KU2dOpX09HQmTJhAWFhYs+cZHBzM/PnzmTJlCt27d2fcuHF1C8U888wzrFy5EovFwoABA5g8eXKr/i3ra3GFKRH5MXCZMeb/HK9vBNKNMXc3UvZRoMwY89SZ1q3PF1eYuvLvqwkJsPDmnNFtPpYrvfS/PTzx0Q4+uud8+vdp35XofZ27V5hSHUd7rDCVA8TWex0DOHsFry11fUZheTXbcosZl+Q5rfmTrh0RS3CAH4u/yXZ3KEqpduJMol8HJIlIoogEAtcBy5w8flvq+oyvdx/FGDwy0XcJDeTqtGje23SIogqdFkEpX9RiojfG1AJ3AZ8AmcCbxpgMEZkjInMARKS3iOQA9wG/FZEcEYloqm57nYynWr3rKOHB/gyO7uzuUBp185gEqmpsvLn+YMuFVbNa6gpVqq1a8zvm1ABNY8xyYHmDbS/We34Ee7eMU3U7EmMMq3cfZUy/SPwtnnkjcv8+EaQnduPVNfuZNa6vRwz/9EbBwcEcO3aMyMhIt975rHyXMYZjx44RHBx8RvV0rpt2tu9oOYeKKrnjwn7uDqVZM8ck8LPXvueLHflc6iFDQL1NTEwMOTk5FBQUtFy46oj9Z7CtfYNSPic4OLjuHglnaaJvZ2v2HgPwqGGVjbl0QC96RwSz+JtsTfStFBAQQGJionOFV9jv2HTFiC6lWuKZfQk+ZN2+QnqEB5EQGeruUJoVYPHjhlFxrN59lN35LphESSnlMTTRt7O1+wpJT+zmFX2216XHEWjx49U1zd9xqJTyLpro21HO8Qpyi6tIT+jm7lCc0r1TEFek9uGdDTmUnWh5ThKllHfQRN+O1u4rBCA90TsSPcANo+Ipr7by/qZD7g5FKeUimujb0brsQiKC/TmvV7i7Q3HakNguJPcO5z/fHdAx4Ur5CE307ei7fYWMSOjmVVMAiwg/HRlHRm4JW3KK3R2OUsoFNNG3k4LSE+wtKGeEF3XbnDR1SDQhARZdalApH6GJvp1s2G/vnx/hJRdi64sIDuCq1CiWbc6lpKr5ebWVUp5PE3072XigiECLH4OivXPq35+OiqOyxsp7G/WirFLeThN9O9l4oIiB0REE+VvcHUqrDI7pwqDoCL0oq5QP0ETfDmqsNrYcKmJIrPsXAW+L69Pj2XGklO8PFLk7FKVUG2iibwc7j5RSVWNjSFwXd4fSJlelRREWqBdllfJ2mujbwcYD9gWNvT3RdwryZ+qQaD7ckktxhV6UVcpbaaJvBxsPFNG9UxDRXULcHUqbXZ8ex4laG0s35rg7FKVUK2mibwebDhYxJK6LV0xk1pJB0Z1JjemsF2WV8mKa6F3seHk1e4+We323TX3Xj4xjV34ZG/Yfd3coSqlW0ETvYpty7CNUvH3ETX1XpkbRKcif/6zVi7JKeSNN9C62NacYEUiJ8cyFwFsjNNCfqWlR/HfLYb0oq5QX0kTvYlsPFZPYPYxOQb61SuMMx0XZd/WirFJeRxO9i2UcKmZQlO+05k8aFN2ZwTGdeX3tQb0oq5SX0UTvQsfKTpBbXEVKtO8lerC36nfm6Z2ySnkbTfQutC23BICBXjqRWUuuTLXfKbtEL8oq5VU00bvQtkP2hToG+mDXDdjvlL0qLZoPtuj0xUp5E6cSvYhMEpGdIrJbROY2sl9EZJ5j/xYRGVpv370ikiEi20TkdREJduUJeJJth4qJjwylc0iAu0NpNzPSY6mqsfG+Tl+slNdoMdGLiAV4DpgMDABmiMiABsUmA0mOx2zgBUfdaODnwHBjzCDAAlznsug9zLZc37wQW19KdGcGRkXwmt4pq5TXcKZFnw7sNsbsNcZUA0uAqQ3KTAVeNXbfAl1EpI9jnz8QIiL+QCiQ66LYPUpRRTUHCysZ5KMXYk8SEWakx7HjSCmbdU1ZpbyCM4k+GjhY73WOY1uLZYwxh4CngAPAYaDYGPNpY28iIrNFZL2IrC8oKHA2fo+R4bgQ660rSp2JqWlRhARYeF2nL1bKKziT6Bubmavhd/ZGy4hIV+yt/UQgCggTkRsaexNjzHxjzHBjzPAePXo4EZZn2eq4EOvrXTcA4fXWlC3Vi7JKeTxnEn0OEFvvdQynd780VeYSYJ8xpsAYUwMsBca0PlzPte1QMdFdQugaFujuUM6KGSPta8ou2+yTPXFK+RRnEv06IElEEkUkEPvF1GUNyiwDbnKMvhmFvYvmMPYum1EiEir2OXsnAJkujN9jbDtU3CG6bU5KjelMcu9wXtcx9Up5vBYTvTGmFrgL+AR7kn7TGJMhInNEZI6j2HJgL7AbWAD8zFH3O+Bt4Htgq+P95rv6JNytpKqG7GMVPntHbGNEhOtHxrHtUAlb9aKsUh7NqZm3jDHLsSfz+tterPfcAHc2UfcR4JE2xOjxMg6dvCO24yR6gKlp0fxxeSb/WXuAJ2JS3B2OUqoJemesC2TkdpwLsfV1DgngisFRLNt0iLITte4ORynVBE30LrDjSCndOwXRIzzI3aGcdTPS4yivtvKBXpRVymNponeBHUdK6N8n3N1huMXQuC6c1ytcJzpTyoNpom+jWquNrLwyknt3zERvv1M2ls05xXWTuimlPIsm+jbKPlZOda2N5N4dZ2hlQ9OGxBDk78eSddqqV8oTaaJvo8zDpQAkd9CuG4DOoQFMSenDextzqajWi7JKeRpN9G2040gJFj/hnJ6d3B2KW80YGUfZiVo+3HzY3aEopRrQRN9GOw6X0q9HGEH+FneH4lbD47tyTs9O/EcvyirlcTTRt9GOI6Udun/+pJPTF286WETm4RJ3h6OUqkcTfRsUV9ZwqKiyQ/fP13fNkGgC/f10qKVSHkYTfRvsPGK/ENtfW/QAdA0LZPKg3izdeIjKaqu7w1FKOWiib4MdR+xdFNqi/8GM9DhKq2r571a9KKuUp9BE3waZh0vpHBJA7wifXe/8jI1M7Ebf7mE6fbFSHkQTfRvsPFJCcu9w7FPtK/jhouyG/cfruraUUu6lib6VbDbDziOl9O+j/fMNTR8WQ6DFT1v1SnkITfStlHO8kvJqa4ed46Y53cICuWxQb5Z+n0NVjV6UVcrdNNG3UmbdhVht0TdmRnosJVW1fLhFL8oq5W6a6Ftpx+FSRODcXh176oOmjO4bSb8eYfxrTba7Q1Gqw9NE30pZ+aXEdg0lNNCp1Rg7HBHh5jEJbM4pZtPBIneHo1SHpom+lXbllWprvgXXDI2hU5A/r36T7e5QlOrQNNG3Qo3Vxr6j5ST10guxzekU5M81Q6P5cMthjpWdcHc4SnVYmuhbIftoOTVWoy16J9w0Op5qq40l6w66OxSlOixN9K2QlVcGQFJPbdG35Jye4Yw9J5LXvt1PrdXm7nCU6pA00bdCVl4pfkKHX2zEWTeNTiC3uIoVmfnuDkWpDkkTfStk5ZUS1y2U4ICOvdiIsyYk9yS6Swj/+jbb3aEo1SE5lehFZJKI7BSR3SIyt5H9IiLzHPu3iMjQevu6iMjbIrJDRDJFZLQrT8AdsvJK9ULsGfC3+PHTUXF8vfsYu/N1/hulzrYWE72IWIDngMnAAGCGiAxoUGwykOR4zAZeqLfvWeBjY0wykApkuiButzlRayX7WIVeiD1D1w6PJdDfj1fX7Hd3KEp1OM606NOB3caYvcaYamAJMLVBmanAq8buW6CLiPQRkQhgPLAQwBhTbYzx6rtn9h0tx2oznKst+jMS2SmIKwdH8c6GHEqratwdjlIdijOJPhqoPzYux7HNmTJ9gQLgFRHZKCL/FJGwxt5ERGaLyHoRWV9QUOD0CZxtOuKm9W4aHU95tZW3N+S4OxSlOhRnEn1jk60bJ8v4A0OBF4wxQ4By4LQ+fgBjzHxjzHBjzPAePXo4EZZ77HKMuOnbo9G/V6oZqZjv9tQAABiNSURBVLFdGBbflVe+zsZqa/grpJRqL84k+hwgtt7rGCDXyTI5QI4x5jvH9rexJ36vlZVXSkJkmI64aaX/G5fIgcIKVmTmuTsUpToMZxL9OiBJRBJFJBC4DljWoMwy4CbH6JtRQLEx5rAx5ghwUETOc5SbAGx3VfDusCuvjCS9ENtqEwf2JqZrCAtX7XN3KEp1GC0memNMLXAX8An2ETNvGmMyRGSOiMxxFFsO7AV2AwuAn9U7xN3AayKyBUgD/ujC+M+qqhor2cfK9UJsG1j8hFvGJrI2u5DNOqulUmeFU3PsGmOWY0/m9be9WO+5Ae5sou4mYHgbYvQYewvKsRl0DH0b/WR4DM98lsXC1fuYN2OIu8NRyufpnbFnYJfjZh8dQ9824cEBXDsiluVbD5NbVOnucJTyeZroz0BWXikWPyGxu464aauZYxOwGcPiNdnuDkUpn6eJ/gxk5ZWREBlKkL+OuGmrmK6hTE7pw3++O0D5iVp3h6OUT9NEfwbsq0pp/7yrzBqXSGlVLW+u17nqlWpPmuidVFVjZX9hhSZ6Fxoa15URCV3556p91Ohc9Uq1G030TtqdX4YxaKJ3sTsu7MehokqWbWp4D55SylU00TtJR9y0j4vO60ly73Be/N8ebDotglLtQhO9k7LyyvD3ExJ0xI1LiQh3XNiPXfllOi2CUu1EE72TduWVktg9jACL/pO52pSUPsR2C+H5L/dgv/dOKeVKmrWctCu/TPvn24m/xY/Z4/ux6WAR3+4tdHc4SvkcTfROqKy2cqCwQicza0c/HhZD906BvPC/Pe4ORSmfo4neCXsK7CNudLGR9hMcYOGWsYl8lVXAtkPF7g5HKZ+iid4JOuLm7LhxdDzhQf7844vd7g5FKZ+iid4JuxwjbuIjdcRNe4oIDuCWsQl8nHGE7bkl7g5HKZ+hid4JWXllJHYPI9Bf/7na26xxfQkP8mfe57vcHYpSPkMzlxN255fqhdizpHNoALeMS9RWvVIupIm+BTZj2F9YoRdiz6JZYxMJD/bn2c+z3B2KUj5BE30LKqut9hE32qI/azqHBnDr2EQ+ycgjI1dH4CjVVproW1BZYwV0MrOz7dZxjlb9Cu2rV6qtNNG3oKLaap/jRkfcnFWdQwKYNS6RT7fn6bh6pdpIE30LKmusJOiIG7e4ZWwiEcH+/O0z7atXqi00e7WgstpKUk/tn3eHziEBzLmwH5/vyGftPp0DR6nW0kTfDJsxVNVYSdL+ebe5ZUwivSKC+NNHmTqzpVKtpIm+GZXV9gux2qJ3n5BAC/deci7fHyji0+06X71SraGJvhk64sYz/GhYDP16hPGXj3dQq2vLKnXGnEr0IjJJRHaKyG4RmdvIfhGReY79W0RkaIP9FhHZKCIfuirws6Gi2goiJHQPdXcoHZq/xY8HLktmT0E5b2/IcXc4SnmdFhO9iFiA54DJwABghogMaFBsMpDkeMwGXmiw/x4gs83RnmWVNVaC/f0I8re4O5QO77KBvRgS14W/rciq61JTSjnHmRZ9OrDbGLPXGFMNLAGmNigzFXjV2H0LdBGRPgAiEgNMAf7pwrjPispqK6GBmuQ9gYgwd1IyeSUnePnrfe4ORymv4kyijwYO1nud49jmbJlngF8BzXauishsEVkvIusLCgqcCKt9VdVYqaqxEhKgid5TjOwbySX9e/L8yt3kl1S5OxylvIYziV4a2dZwnFujZUTkCiDfGLOhpTcxxsw3xgw3xgzv0aOHE2G1r70F5QCEBPq7ORJV32+mDKDaauPPH+90dyhKeQ1nEn0OEFvvdQyQ62SZscBVIpKNvcvnYhH5d6ujPYtOriqlXTeeJbF7GLPG9eWd73PYeOC4u8NRyis4k+jXAUkikigigcB1wLIGZZYBNzlG34wCio0xh40xvzbGxBhjEhz1vjDG3ODKE2gvu/LKQITgAB2B6mnuuvgceoYH8egH27HZ9CYqpVrSYhYzxtQCdwGfYB8586YxJkNE5ojIHEex5cBeYDewAPhZO8V71uzKLyXY3w8/aaxXSrlTpyB/HpyUzOaDRSzdeMjd4Sjl8ZzqgDbGLMeezOtve7HecwPc2cIxvgS+POMI3WRXXhkhfbTbxlNNGxLNv77dz58/3sGkQb3pFKTXUpRqivZLNOJErZXsY+WE6ogbj+XnJzx61UAKSk/wd11fVqlmaaJvxN6CcmxGR9x4urTYLlw7PJaFq/eReVjXl1WqKZroG7ErvwywT6ilPNuvL0+mc0gAc5duxaoXZpVqlCb6RuzKK8VPIERH3Hi8LqGBPHzlADYfLOJfa7LdHY5SHkkzWSOy8kpJiAzTETde4qrUKC44twdPfrKT3KJKd4ejlMfRRN+IHUdKSe6jUxN7CxHhsasHYTPw8PvbdIESpRrQRN9A+YlaDhRWcF6vCHeHos5AbLdQ7rv0XFZk5vPRtiPuDkcpj6KJvoGsvFKMQVv0XuiWsQkMjIrg4fczKCyvdnc4SnkMTfQN7Dhin+Omf29t0Xsbf4sfT/04leLKan773lbtwlHKQRN9AzuPlBIaaCGma4i7Q1Gt0L9PBL+45FyWbz3Css0N595TqmPSRN9A5uESzusdjp+fjrjxVreP78uQuC78v/e2caRY561XShN9PcYYduaVktxb++e9mb/Fj6d/kka11cYDb2/WGS5Vh6eJvp68khMUVdSQrP3zXi+xexi/mTKAVbuOsnC1Lj2oOjZN9PVkHrHPl6Itet9ww8g4LhvYiz9/vIPNB4vcHY5SbqOJvp6djhE32qL3DSLCX6an0isimLtf30hJVY27Q1LKLTTR17PjcAl9OgfTOTTA3aEoF+kcGsCz16VxqKiSh5bqkEvVMWmir2fHkVLO024bnzM8oRv3XXouH245zKJvst0djlJnnSZ6h+paG3sKyrTbxkfdcUE/Lunfi8f/m8nafYXuDkeps0oTvcPeo2XUWA39deoDn+TnJzx9bSqx3UL52Wvfk1ei4+tVx6GJ3uHkhVjtuvFdEcEBvHTjMCqqa7nj3xuoqrG6OySlzgpN9A6Zh0sJsAh9u3dydyiqHZ3bK5ynfpzK9weK+LVenFUdhCZ6hx1HSujXoxOB/vpP4usuT+nD/RPP5d2Nh/j7F7vdHY5S7U6zGvapD7YdKmZgVGd3h6LOkjsvOodrhkbz9GdZfKCTnykf5+/uADxBXskJjpZVkxKtI246ChHhiWtSyCms5JdvbaZHeBCj+ka6Oyyl2oW26IFth4oBGBStLfqOJMjfwks3DiOuWyi3LV5f93uglK9xKtGLyCQR2Skiu0VkbiP7RUTmOfZvEZGhju2xIrJSRDJFJENE7nH1CbjC1kPFiMCAKG3RdzRdwwJ59dZ0woP9mfnKWrKPlrs7JKVcrsVELyIW4DlgMjAAmCEiAxoUmwwkOR6zgRcc22uBXxpj+gOjgDsbqet2GbnF9OvRidBA7cnqiKK6hPDqrJFYbYYbX/6Ow8WV7g5JKZdypkWfDuw2xuw1xlQDS4CpDcpMBV41dt8CXUSkjzHmsDHmewBjTCmQCUS7MH6X2HqomEHamu/QzunZiUW3pFNUXsOM+d/qgiXKpziT6KOBg/Ve53B6sm6xjIgkAEOA7xp7ExGZLSLrRWR9QUGBE2G5Rn5pFXklJ7R/XpEa24VFt6ZztKya6xd8q3fPKp/hTKJvbE29hneZNFtGRDoB7wC/MMaUNPYmxpj5xpjhxpjhPXr0cCIs18jItYejiV4BDIvvyuJbR5BXUsWMBd9qN47yCc4k+hwgtt7rGKDhwOMmy4hIAPYk/5oxZmnrQ20f23LsIy0GateNchgW341Ft6aTX3KCH72whj0FZe4OSak2cSbRrwOSRCRRRAKB64BlDcosA25yjL4ZBRQbYw6LiAALgUxjzNMujdxFth4qJrF7GOHBOge9+sGIhG4smT2KqhorP35xDVtzdOil8l4tJnpjTC1wF/AJ9oupbxpjMkRkjojMcRRbDuwFdgMLgJ85to8FbgQuFpFNjsflrj6J1jLGsPFgEakx2m2jTjcoujNv3zGGkAAL181fwzd7jro7JKVaxanxhMaY5diTef1tL9Z7boA7G6m3msb77z1CzvFKCkpPMDS+q7tDUR4qsXsY79wxhpte/o6ZL6/jzz9KYdqQGHeHpdQZ6dB3xn5/4DgAQ+M00aum9e4czJu3j2ZofBfufWMzTyzPxGrTWS+V9+jQiX7jgSJCAiwk6xz0qgVdQgP516yR3DAqjpe+2sttr66nVBcbV16iQyf67w8cJzW2M/6WDv3PoJwUYPHjsatT+MPVg/gqq4Bpz3/D7nwdkaM8X4fNcFU1Vrbnlmi3jTpjN46K59VZ6RSWV3Pl31fz1vqDuoCJ8mgdNtFvySmm1mY00atWGdOvO8t/fj6psZ154O0t3PfmZspO1Lo7LKUa1WET/ckLsUPiurg5EuWtencO5rX/G8W9l5zL+5sOccW8VWx0/F4p5Uk6bqLff5yEyFAiOwW5OxTlxSx+wj2XJPGf20ZxotbG9Be+4fH/bqeyWhceV56jQyZ6q83w3b5CRiR0c3coykeM6hvJp/eO57r0OBas2sfkZ7/iu73H3B2WUkAHTfQZucUUV9YwLqm7u0NRPiQ8OIA/TkvhP7eNxGbg2vnf8uDbWzhadsLdoakOrkMm+tW77beyj+mniV653ph+3fn4F+cze3xflm7M4aInv2TBV3uprrW5OzTVQXXMRL/rKMm9w+kRrv3zqn2EBvrz0OX9+fgX4xmW0JXHl2cy6dmv+CTjiA7FVGddh0v0ldVW1mcfZ9w52ppX7a9fD/vKVa/MHAEGbv/XBqY+9zVFFTWnLeqgVHvpcIukrt9fSLXVxljtn1dn0UXJPTk/qTtLNx7i2RW72HGkhPDgAI7vOsrYcyKxz+itVPvocC361buPEmAR0nXEjTrL/C1+/GR4LCvvv5CE7mFU1Vq5YeF3TJm3mvc3HaLGqn34qn10uET/VdZRhsR1JSyow32ZUR4i0N+P3hHBDIntwl+mD6baauOeJZu48Mkveel/ezimo3SUi3WoRH+wsILMwyVc0r+nu0NRCj8RfjIilk9/MZ6FNw8npmsIT3y0g9FPfMHdr29kzZ5jeuFWuUSHatZ+knEEgMsG9nZzJEr9wM9PmNC/FxP69yIrr5T/fHeApd/n8MHmXOIjQ7lycBRT06JI6qXTaavW6VCJ/uNtR0juHU58ZJi7Q1GqUef2CufRqwYyd3Iy/91ymPc2HeL5L3fzj5W7Se4dzlVpUVyREkVcZKi7Q1VepMMk+gPHKli//zgPXHaeu0NRqkXBARamD4th+rAY8kurWL7lMMs25/KXj3fyl493ck7PTkxI7snFyT0ZFt9V11RQzeowiX7pxhxEYNqQaHeHotQZ6RkezMyxicwcm8jBwgo+3Z7HFzvyWLh6Hy99tZfOIQGcn9Sd0f0iGdU3kr7dw3S4pjpFh0j0Vpvh7Q05jO4bSVSXEHeHo1SrxXYLZda4RGaNS6S0qobVu46yIjOfVbsK+HDLYQB6hgcxqm8k6YndSIvtwnm9wwnQFn+H1iES/acZR8g5Xslvp/R3dyhKuUx4cACTU/owOaUPxhiyj1WwZs8xvt1rfyzbnAtAkL8fA6MiSI3tQlpsF1KiOxMfGYbFT1v9HYXPJ3pjDPNX7SW2WwiXDtDRNso3iQiJ3cNI7B7G9SPjMMZwsLCSzTlFbD5YxOacIpasPcgrX2cD9uR/Ts9OnNcrnKRe4ZzXuxNJPcOJ6hKifwB8kM8n+k+357HxQBGPTxukv8CqwxAR4iJDiYsM5crUKABqrTZ25Zex9VAxu/JK2ZlXxjd7jrF046G6egEWIbarvV5cN/sjPjKM+MhQencOJjzIX/v/vZBPJ/qyE7U8/t9MzunZiWuHx7o7HKXcyt/iR/8+EfTvE3HK9uLKGnbllZKVV8aBwgoOFJaz/1gFG7KPU9pgHdzQQAu9I4LpFRFMr4ggenUOrnsdGRZIt7BAuoYF0iUkQEcCeRCnEr2ITAKeBSzAP40xf2qwXxz7LwcqgJnGmO+dqdterDbDr5duJed4Ba/fNkp/6ZRqQueQAIYndGN4g/mfjDEcr6hh/7FyDhRWkFdSxZHiE/afJVWs33+c/JITVDcxR0/nkAC6hgbQNSyQbqH2PwARwQF0CvYnPMifTsH+dHL8rP86PCiA0CCLXkB2oRYTvYhYgOeAS4EcYJ2ILDPGbK9XbDKQ5HiMBF4ARjpZ1+Vyiyr53QcZfJKRx68mncfIvpHt+XZK+SQRoZujlT4krmujZWw2w/GKao6UVHG8vIbCimqOl1dTWF5NUUU1hRU1HC+3799+uITSqlrKGnxLaIrFTwgJsBAc4EeQv4WQQPvzYMfzIH/765AAC0EBfgRYTj4Efz8/Av1/eB7g70eAnxBg8cPfIgRa/PB3lD1Zz+JnP2eLCBY/QcQeg0XEvr3uuWP7yTKOfX5+gp+jvp8f9Z67v6vLmRZ9OrDbGLMXQESWAFOB+sl6KvCqsU/M8a2IdBGRPkCCE3Vd5oq/r+J4eQ2HiioJtPjx2yn9+b/z+7bHWymlsE/fENkpiMhOzi/iY7MZyqvtCb+sqpZSx8/6r8tP1FJVY6WqxkZVrZWqaqv9Z42NqhorZSdqOVpWzYkaK5U1VqpqrNRaDTU2GzVWg9XmeXME+Yn9D4lg/yOA/T/8HH88BOgRHsSXD1zk8vd2JtFHAwfrvc7B3mpvqUy0k3UBEJHZwGyAuLg4J8I6XVJP+1wg5/TsxFWpUcR209vElfI0fn5CeHAA4cEB0Ll93sNmsyf9WquhxmpP/jVW++tqq41am42aWscfhlr7fpsxWI3BZjPYjL3712Yc208+t1FXxmrs5Ww206AsjrI/lMEYDPbtxlD3nJPPbfb9YYGWdvn3cCbRN/a9o+Gfy6bKOFPXvtGY+cB8gOHDh7fqz/Hfrk1rTTWllI/x8xOC/CzobOR2zvwz5AD1h6zEALlOlgl0oq5SSql25Mxl7XVAkogkikggcB2wrEGZZcBNYjcKKDbGHHayrlJKqXbUYoveGFMrIncBn2AfIvmyMSZDROY49r8ILMc+tHI39uGVtzRXt13ORCmlVKOc6sEyxizHnszrb3ux3nMD3OlsXaWUUmeP3pGglFI+ThO9Ukr5OE30Sinl4zTRK6WUjxP7dVTPIiIFwP5WVu8OHHVhOO6k5+J5fOU8QM/FU7X2XOKNMT0a2+GRib4tRGS9MWa4u+NwBT0Xz+Mr5wF6Lp6qPc5Fu26UUsrHaaJXSikf54uJfr67A3AhPRfP4yvnAXounsrl5+JzffRKKaVO5YsteqWUUvVooldKKR/nM4leRCaJyE4R2S0ic90dT1uISLaIbBWRTSKy3t3xnAkReVlE8kVkW71t3UTkMxHZ5fjZ+AKkHqaJc3lURA45PptNInK5O2N0lojEishKEckUkQwRucex3es+m2bOxes+GxEJFpG1IrLZcS6/c2x36efiE330jkXIs6i3CDkwo70XIW8vIpINDDfGeN0NICIyHijDvobwIMe2vwCFxpg/Of4IdzXGPOjOOJ3RxLk8CpQZY55yZ2xnyrGGcx9jzPciEg5sAK4GZuJln00z5/ITvOyzEREBwowxZSISAKwG7gGuwYWfi6+06OsWMDfGVAMnFyFXZ5kx5iugsMHmqcBix/PF2P+n9HhNnItXMsYcNsZ873heCmRiX9PZ6z6bZs7F6xi7MsfLAMfD4OLPxVcSfVOLk3srA3wqIhsci6Z7u16OFcdw/Ozp5nja6i4R2eLo2vH4ro6GRCQBGAJ8h5d/Ng3OBbzwsxERi4hsAvKBz4wxLv9cfCXRO70IuZcYa4wZCkwG7nR0ISjP8ALQD0gDDgN/dW84Z0ZEOgHvAL8wxpS4O562aORcvPKzMcZYjTFp2NfUTheRQa5+D19J9M4sYO41jDG5jp/5wLvYu6a8WZ6jX/Vk/2q+m+NpNWNMnuN/TBuwAC/6bBx9wO8Arxljljo2e+Vn09i5ePNnA2CMKQK+BCbh4s/FVxK9zyxCLiJhjgtMiEgYMBHY1nwtj7cMuNnx/GbgfTfG0iYn/+dzmIaXfDaOi34LgUxjzNP1dnndZ9PUuXjjZyMiPUSki+N5CHAJsAMXfy4+MeoGwDGU6hl+WIT8cTeH1Coi0hd7Kx7sa/r+x5vORUReBy7EPtVqHvAI8B7wJhAHHAB+bIzx+IucTZzLhdi7BgyQDdx+si/Vk4nIOGAVsBWwOTY/hL1v26s+m2bOZQZe9tmIyGDsF1st2Bvebxpjfi8ikbjwc/GZRK+UUqpxvtJ1o5RSqgma6JVSysdpoldKKR+niV4ppXycJnqllPJxmuiVUsrHaaJXSikf9/8BDn0u6EpEKNYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#visual demonstration of chi2 distribution with dof1_lp\n", "xrange = np.arange(0,30,0.1)\n", "plt.plot(xrange, chi2.pdf(xrange,dof1_lp), label='$\\chi^2$ dist with DOF='+str(dof1_lp)[0:4])\n", "plt.axvline(chi2_lower,color='orange', label='95% confidence bounds')\n", "plt.axvline(chi2_upper,color='orange')\n", "plt.legend()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Note for Yeager et al. (in revision):\n", "\n", "If the LSW formation time series in Table 2 look normally distributed, then: \n", "\n", "1) the f-test portion of this notebook can be used to determine if the standard deviations from two of formation time series are statistically different from each other (e.g., LS+SPG-west versus SPG-east).\n", "\n", "2) the chi-square portion of this notebook can be used to put confidence bounds on a single sample standard deviation.\n", "\n", "If the time series do not look normally distributed, then I would just ignore trying to put confidence bounds on the standard deviations. Figuring out a more appropriate method does not seem worth it since I don't think the confidence bounds are that important for R1." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Footnote\n", "The F-statistic and F-test was named for Ronald Fisher (as was the Fisher-Z transform), which are now seen as controversial choices. See: https://nautil.us/issue/92/frontiers/how-eugenics-shaped-statistics " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:miniconda3-rapcdi-analysis]", "language": "python", "name": "conda-env-miniconda3-rapcdi-analysis-py" }, "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.8.0" } }, "nbformat": 4, "nbformat_minor": 4 }