Skip to content

Instantly share code, notes, and snippets.

@ljmartin
Created April 27, 2020 05:28
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ljmartin/ca6d05e211db66fbe87e0c862860184f to your computer and use it in GitHub Desktop.
Save ljmartin/ca6d05e211db66fbe87e0c862860184f to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import pymc3 as pm\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from seaborn import kdeplot\n",
"\n",
"import statsmodels.api as sm"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-2.7112278332640463\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.lines.Line2D at 0x123779110>"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2dd5zc1LXHf0czs9W79tped+O1jXHFjTVgGwjFNhgSQyC0kABphIRUEhLzCASSQEiBJISEUENIoYTy4GGqscEYDC649467171sm3LfH5qrudJIGs2MpJnV3u/ns5+d0Wikq9HVueeeewoxxiCRSCSSYKIUugESiUQi8Q4p5CUSiSTASCEvkUgkAUYKeYlEIgkwUshLJBJJgAkXugEiXbt2ZXV1dYVuhkQikbQpFi1atI8xVmv2WVEJ+bq6OixcuLDQzZBIJJI2BRFttfpMmmskEokkwEghL5FIJAFGCnmJRCIJMFLISyQSSYCRQl4ikUgCjBTyEolEEmCkkJdIJJIAI4W8RCKRuMiKHYexZNuhQjdDo6iCoSQSiaQtM3f9Pnzp8Y8BAFvuvajArVGRmrxEIpG4xE3/+aTQTUhDCnmJRCJxicNN0UI3IQ0p5CUSiSTASCEvkUgkAUYKeYlEIgkwUshLJBKJS4QU0l7vPNRUwJakkEJeIpFIAowU8hKJROIBr6/YXegmAJBCXiKRSFyDhNfF4k4phbxEIpG4RCSkoHNlCQDgpO4dCtwaFVeEPBE9QUR7iWiFsO1OItpBREuSfxe6cS6JJAh8sGEf4glW6GZIXCaeYJgwsAsAYP2eYwVujYpbmvyTAC4w2f4Hxtjo5N9rLp1LImnTvL++Adc89jEeendDoZsicZlYIoHKEjUl2J/eWV/g1qi4IuQZY3MAHHDjWBJJ0NlzpAUAsKnheIFbInGTRIIhwYCK0lChm6LDa5v8d4hoWdKcU+PxuSSSNgFfnJPGmmBxtCUGAKgoaT9C/iEAAwGMBrALwH1mOxHRDUS0kIgWNjQ0eNgciaQ4oKSUZ0yK+SDx85fVJcnVu44WuCV6PBPyjLE9jLE4YywB4FEAp1rs9whjrJ4xVl9bW+tVcySSooEo8z6Stgf3i2+NJQrcEj2eCXki6im8/TyAFVb7SiTtEanHB4uWpHDvUFpctZhcaQ0RPQ3gbABdiWg7gJ8DOJuIRkPty1sAfNONc0kkbZ1DjWqQjPSgDCYdygIo5BljV5tsftyNY0skQeOu/1sFQNrkg0r36lJMHtYdb6/aA8YYqMD2ORnxKpFIJC4wrGc1AOCC4T3x9qo9AICNDYUPiCqueYVE0o6QenywGNKjCkeaozi5T8dCN0WH1OQlkkIhpXygiCUYworeNFNoUw0ghbxEUjBmLN+FRVsPFroZEpeIJxjCIVWkXjqmNwAgFi/8SC6FvERSQC576MNCN0HiErFEQtPkp56sepBH44X3mZdCXiKRSFwgnmBa+b9ISP3fUgSBUVLISyQ+wwWAJFiINvmSsCpajzQXvnCIFPISic+EFCnkg4ioyfNc8r98dVUhmwRACnmJxHeKYTFO4j6xOENYUUXq6QPUwiGVJWEs234ID7+3ETsONRWkXdJPXiLxEcYYYjKfQSARNfluVaUAgOU7DmPagx8AAF78ZAfe/OFZvrdLavISiY+YCfjN+2TxkCAQSyQQTq63hEzWXQ40tvrdJABSyEskvmJmqnllyc4CtETiNjrvGiVdtBaqpq8U8hKJj7Sa+E0nZKKyQCB615gtrhdqvV0KeYnER2ImQl5mowwGsXhKkzemN1ApjJSXQl4i8REzm7xchw0GasSrKlIVqclLJO0THuZ+y/mDtW3SbT4YqLlrrG+mUqBkZVLISyQ+whdee3Ys07YVQ6ZCSf7EhIVXM6QmL5G0A7gmHwmlHr2ySKhQzZG4SNwk1bBIoQZzKeQlEh+JJjX5SIjwty+NBQB07VBSyCZJXGD+5gPYdbjZNqr1SFNh8ti4IuSJ6Aki2ktEK4RtnYnobSJan/xf48a5JJK2TCyhavJhRcHwXsVVQUiSO9c9MR8A8MGG/Zb7HG2J+dUcHW5p8k8CuMCwbTqAdxhjgwC8k3xfVDS1xjFvo/VNkUjchmvy4RBpHhjSg7Lt0xSNF7oJlrgi5BljcwAcMGy+GMA/kq//AeASN87lJre+uAxXP/oRth1oLHRTJO2EmGCT5+ZbGQwl8RIvbfLdGWO7ACD5v5uH58qJNbuPAgCOFWgaJWl/cD/5sEKaS50U8e2H5dsP+x78VvCFVyK6gYgWEtHChoYG3847f/MBTchLJH7B0xpEwooW/yg1+eBw/xWjbD//3INz0f/W13xqjYqXQn4PEfUEgOT/vWY7McYeYYzVM8bqa2trPWyOnisenqe9lm7KEr/gfvIRRdFc6mTEa3A4f3gP7fUL35qAi0f3KmBrVLwU8q8AuC75+joAL3t4rqxoLuJFEkmw4X2vLJKyycvcNcFBDIY6pV8NvnBKH9P9/LznbrlQPg1gHoDBRLSdiL4G4F4Ak4loPYDJyfdFgbG4LhUocZAkuGxsOIa66TMwe41+ApsS8iFNk3fzef/3x1vx5Aeb3TugxBGj+nYCUXpgmxj0JuJn2mFXKkMxxq62+Og8N47vNq1GIS9lvMRlFn96CADwf8t24pwhKZ+D5mTfKxU0eTfd7257SQ1VuX5if9eOKclMaUjBaf07p223KtoeSzCEfQp0LvjCayEw5vSWs2WJ22jukQaNrblV0OSTM8h7X1/ja9sk7tMST6DERGqHTYqHAP5q8u1TyBs0+QfeWY9zfv9uYRojCSQ8R42xEBQ315RHQiiNtMvHL5AcaYqixMQ0Y5WV0s86v+2ykLdRyM9YvqtALZEElZ++sBxA+gJbcyyOkEKIhBREQkBtVSkmDuxSiCZK8mTPkWZMe3Au/nrNKdi877hprd5isMm3S1XCKOQlEjd5ZM5G7bUxh3hzNIGycOqxqy4L+6rVSbLnww37cPC4WoT7uBA4+fryXdhzpAX3v73W8rtWSSl5DiM/aJ+afNxfF8qXFm/HD59dite/fyaWbT+EK8ed4Ov5Jf5yz2spG7vxIW+OxnUeGGFFKViBZ0lmWmJxfPGxjzGqT0cQEZZsO4Tnvjkep/bvrOUeysU7r81517Q1jC6UHMaYJzmfH5y1AQAw9U/vAwAuHdvHchonCRZGTb41ltDd+5BCnmjyhxuj6FgRcf247Q0+61+6/bC2bc66Bpzav7MmK3IRGTHjYo2HtEtJY2WuicYZ5m8+gD1Hml09n3GFvbFVBmO1F4xKQ5zpS8SFFPJEq7v+yfmuH7M9svtwuiw40qzmhed30U4xLC8x16OlTd5jrIT8pn3HcMXD8/DZP8919XxRg/1NRtwGgz1HmjMmtzOaa4zVg9zU5A81tmqvuZ++JD/eWrUnbRu/X+QgYrl3p3LT7UY3bi9pn0Le4gf+v6U7AQANR1tcPd+Og/pqMXKhLRicds87mPagXiEwamhGc42xDmhYIcRdWoRrjkqHArcZ2Se9sAuX6fxe8/e/uexk02OcOagragymsyl/mONeIzPQPoW8hSZ/vMUbDdtY99EYICNpu2xq0LvNGYX8swu36T+PM535zk1zDZNJi10naqoQqr/zHS+vTL5T35eEzcXpP792GhbfMcWT9jlBCnmB4x7llTc+etKbIrhkErRpmnzIPSEfjcl+5TbmssKwzpK8f8ZZW7HQPoW8hblmaM9qT85nzBcel3kUAkumWxtPJHQLrwq5Z5OfvdY0m7ckD8w88Ywzc17SUQr5IsJKk//Fq6s8Od/UET1175duk4tiQcUo5Mee0En33twm746Qf3D2Bu319RPqXDlme8dsfS6UJuRVeZJJyE8a2g03fmYghveqxri6GvcamYH2KeR9XNkG0jvFLc8v8/X8kvyZvWYvnl3wacb9jLO2zpWluvfxBEOIRO8axTWfaXGtZ/GnB1E3fQa2H5T1i/PhVzNWp20z3mOuNFpFt3Ieu24cpk8dgqqysK/pzdunkPc5rYFx8YZ3hqfnf4pVO4/42hZJbnzlyQX46QvLMxZ7MH5qVO681ORF4cODdxZtPejKsSUpjOY1vuDqNJAyrCi+pjWQQt6CQ42t2HesBWN/+TZW7jyccX87jJ2CRzze+uJyXPjA+3kdW+Iv05OJx6wwanlG7S6RMARDhci1B97MXGCcRUrypyXpqnrZWLXq094jqkmnJOzstw4plJad1EukkLfgo037MWddAw4cb8WjczblfC7GGGYs02e5/NzIwtd9lOTGswu3oW76DMsgKC7jv3x6P5RFFCzYotekVU1ecKEkcq3G6zfOGpC2rVgXA9syL3yyHUCqIMjuZIT8ih3OZuVuxkY4oX0KeQc2+dJwCOGkxu3U+6E1lkhzw/zk0/TpcucOJdJXvo2zzyJgjptz+netRHM0gQPHW3WfGyNew4p7mnxFiZr4TNTepZD3hl2Hm7DCMMOfNsqZ8hZSSOau8RonmnxpRMHPX1ZLqTm1mV75yDwM//mbum0Hj0fT9ovGEr4v/krcxUwwb9l3XHO5s7KSROMJnRAOKYS4Sw8876diyTkZk5Efg7p1MN0+/tez0jT3mooSR8d0MzbCCe1WyBt9XY2UhhUcbFQFtFNN3ixfyM3PLUk/fzwhUxu0cYwl+461xHD279/VPKcUi/4VSzBdBaFwyD0/eX4Y8fh3z/DGLbi90LvGPPeMGVZVoIyEfE4v7bmQJ6ItRLSciJYQ0UKvz+eElngCnSv1o24fw83ceSiVfS6fG3KkOWW+ef8n56BbVSlaYwnXtDeJu6zdfRQX/+WDjInHZq7WBx7xpHNz1jUAUGMiv35Gf1SW6Ot+RuMJnaadbVqDeIJZJrjjpiIxvH6nSRZFiXOyuTdOhXzYo/TSVvilyZ/DGBvNGKv36Xy2tMZUIT+ga6W2TSzkAADffXqx9to8f0X29O1cgYqSEBpb42iSmSiLkt++sQZLtx3CvI37s/qe0ca6+0gzFIXSopu37m/UDSChLCNeb/zXIgy5/Q0AwG/eWKMNKkBKIJnVGpVkz4cb9uH99fsc7x+xKNptxKv00la0y97QGkugNKzgT1eN0baVWiQXAqyLjFghPnhGKkvDONYSw/l/9C8LncQ5XBuLCQO7k0Vy4zrPI3M2QTF4zvBjirOAbKfubwupbx96dyOufSKVN54fRlRY+nZ2bm6Q6Mk2aNHKRGfEzcV2J/gh5BmAt4hoERHdYPyQiG4gooVEtLChwVo4uklrLIGSsILBPapQXaYm9Q/baD+NrdklLhMfPCORkIJZa/bicFP6gqyk8PB+EBUE75Mfbsn4PWNJyViCQSF9rnGusX/p9FT5x/ISBU3ReMYgKyOi186WZAFpPohUlKaEvDGlhtcs2nowMPUSjmf53DsliJr8RMbYWABTAdxERGeJHzLGHmGM1TPG6mtraz1vDGMM8zbtRzTOUBJW8NTXTuMfWH7Hzex+MjiluFm3+ygAvXB+OVlnwIr9x1rScrkzprovRuMMddNnYM66Bk3I9+ucMhN2Ki9BPMFwPMtqYc8JKYzP/v27eH7Rds2sWC5o8k0+ViHbdbgJlz30IYbc/ob06rEhcDZ5xtjO5P+9AF4CcKrX57SDuz0tSSYJ4yLX7jd3sqASc2i3lzK+uFm/91jatkwJ5e57e12aS+wPJg3STd//9dFWrN2t9j1xoC+NqI9gS5bar9G758f/XYq/vacG7YkLr68u25n1LCFXRGVo2fa2n4TPq59NcdFt1tH5vDw4EVUSURV/DWAKgBVenjMTxsT+PFjELg+4kxTEj76/WffemPb1/OHddeeTBIeIQmk2+YkndtUN6G+t2oPLHpoHQK808EVSN+Im+DHEhdeDjVG8tnx33sd2Qki4rvvfXoeNDekD5s5DTTjc2DZMlWKKil9cPBxXjeuLm84ZmPdxg6bJdwcwl4iWApgPYAZj7A2Pz2kLn0Z2LFfLcXGZy9dBzBZgnQRPvbdOL9Tvf2ud7v2wnmoZMWmuaRtEsvBQKQkraYvzCpEu26TI/M0H0s7jZtI8Y9u3Hjhusae7iAvU76/fh68+uSBtnwn3zsJ597/rS3vyRVzAvnZ8He69bKTOFJYrIUXxtaaEp0KeMbaJMTYq+TecMXa3l+dzAhfy910+CoAg5BnDml9egJe/MzHtO69ksMkC6Q/W8h36kGc+U5CafHHDb49CwNm/m4266TMyfoexdCGtkLW3hZj6gs8s3agrXFtVilP7d06brfpVMcoqBa+RfcdaTbcXG1eN65u2LeqCmSWskGPzrhu0OxdK7rrEp5ZiXueySAil4dxG6kyan50bnsxjUzx0ryoDoHrCbNmv5mK/sj79YRd5bO7mNIEWUigtzTDn5smDtdc88voLf5uXa5M1WqJxlIaVNCHvl7ueMYBslxCI9dzCbbj52fTo72KGK2RlkdTv+bUz++v2mTqiR/bHVVTXWr/WStqdkOeaPH+4+iT9iL8ysU63PVsymWH49GzuhvTgir8IFX0kxYHoHWIsxm1GS0y/cGpnrikXomCd+lY7oYkLeYPC4Yb26YQ/vL3e8rOfPL8MLy7e4Us73IIL4dk/PlvbVl0WwcybUw6CNZXO8tWI8H7hl27X7oQ87/BcKFeXRbDl3otw5bgTdNtFenfKLaBk+faUyaZHdZnlfnPW+xMfIHFOtlkC//SOXsCFFLI0zRmLhrgFdws2avL5prWNJ1haNk0zZq7ek3GftkSCqfeqZ0f9819eEtZe53L3+Bjsl5tpuxPyfOpqZV4xE/JO0hqYLdg2HEtNV685rZ/ld405xyWFJ5so5ynDumPrfn2ZPYXIUksXN2e7EN+1Q6nt56XhUNrgkq8nx2Pvb8LYX76NPUes8+AE0S8+wZipy7OYe4gIuHZ8P3z33BMdH5fXE/DrNwtn3iU4HDzeiv3JRR8rDSpXIX9C54q0beLD5ua0XOI9//OSfQUoEatn1eqWi/0i7DDfCae6PIx9x6wXac2UjXyFyavJojdz1+9Dz05lGHtCTVquJ7NztHUfgwQzL+knmsMIhF9cPCKr42qavE82+XYl5Mf88m3ttaUmb3JTnUzdjzSn+/5KT5r2gZk5hMhauIrPdraafKa9zeRGvjb5tXvUKOAf/XcpAODMQV3xTx4pnuTjzekJ3cpydGIoFuKJhIUmLwj5HB7x3YfVQfpocxQdSr0Xwe3OXMOxFPLCVOyNH5yJy0/pg2gGm2ZzNI5/ffSpbtvA2kop5NsguWhXVuYQq5QCoquh0/S0nEzFoo+3xtIWgfN11zN6Dr2/fl/atd32UnqMY1M0jqMmyk9bIRpnpnJCJ+RzOO4TH6iBk6/7FKTWboW81cMlavKDu1ehtqo0oyZvLPkHqMFWfmaak7hDLu6sZuY8AvDSEnNvEvEMbisCzdEE3lypXwD1wijw7X8vcrTf9Bedm72KjVgiYSHkRZt89vevawfVI6c6GZDpNe1GyG8w5CSxyv0sTp+JCOGQgliC2fq0mk3LYwmGW9twB2+v5PLgWSkBjS16bbcqOTUXNXm3zTUzV+/B2BNqsjpmLsxeq/cIs5J12ebkKSZicWa6dpeLYBf57RdGAgB6dbL2uHOTdiPkJ93/nu69lSZv9DHmnVSs8GTEzBMjGmdaMMi4Ou8fOok75DL7slqYH1BbqXt/2Sl9AKhuu5xsF0WdyBc/1/gXbT2IIbe/7sjFsq3xzIJt2JshEjkXed+xXNXk3UxlYUe7WngVsRLyRi8Ybj97Yu5m/HDySabfEZNLfe+8QVi67RB2HmrStt1wVv5JjST+YJUdcGBtJTY2mOeAWbr9sOn22y4aijteXonHrq1HazyBLpUl+NoZ/VFblXKDzNZeTjlYgcVBxW0efm8jmqOJtFTLnJmr9+LNlf7YngtBLveDK5LZxmLkSrvR5I04vTncZmpXOEQckW+efBKqyyO6xTjpPdl2sFpEPdwUw5Z7L8Ifrxzt6DgdyyMY3qsjXvjWBNRUlqB7dRnCIQV9Da622fqwZ9Ic59xyDvonZxB/v36co+94zY+TXjmcSJaLzYWid6dyXDq2t+lnP7toKIDcfltuKfYrE2W7FfJiPgo7bjpHDXII2fgzp2cgBDbvO66lKB7eq2OOrZT4jZX5hPumXzKmt86OPqRHlen+3WwinEXcdKFbf/dUnNClAj+aPBhPXF+Pc4Z0Q1VZuOCBSsblrLYSN9WSLBNqBr+mXIarsM/BUO1GyBtDvascTmE/O1Itn/a39zZa7mO0rb28RM1aeVL3DggphB4dUw+8k9z0ksLhRLsSH85PDzTa7JmZUX07YURvtU84SVhlJRieueF0zROkJKzg3CFq/YKQQmnZIf3GeP5CDzpOaInFse9YC/YeMbfJ86yyuWjyXEnwy/uu3Qh5q2RRmXASkWhlV315yc608z5/43jc+blhObVF4j3ZCqBGE1/4iSd2yeoYU4b1cHxuKz/+0weYn1MtJp6bUG2NJbBgywHTz/p3VU1Cb63KnK/G7PxLMlTbKjQfJBMJvrNmr+nng3uoA/PIPp2yPjb32PFr8G03C6/ZJIL601Wj0ZJcSKoqy/wTGbW/7557Iv48S80saRwjKkvD6O5wKi/xHycpLM4d0g2zLB7+FXedbznFt4I7AcQSDJmCRLMdhBQi5BoL9bV/LMD769OzpgLZeYaYLTAes/FWKwYyrdl95qRavHfL2ejXpdJ2PzM0TV4uvLqL6DVz/YQ6230vHt0bVyQLBtRUlqAkpOCikfqq9x9t2o9H56g1NY0P3kndU3Zas5nAiN7SRu82K3YcdqUQgxMh+tdrxppuv+/yUehQGs6qqhSQitnIZCr6dH+jLhGak8lpSMm9XoGVgC+LpFfCGpDU7FfcdX7a/n6WunMLJ/0gFwEPpAZ1aZN3GXGx7M5pw7P6bv+ulWkC5KpHPsLdr60GkN6JxVmD2QSib+cKbLn3Iq3uqyQ/1u85is/+eS5+9+bavI7DGHMkkMoiIbzynYl47pvjddu5H3y28Ic+mkE7Put3s3XvnbhGhohcT4SlEKHVmD9fISjkfCHZC4+f5xdtx0uLt7tyLL6GV9/P/RiXlE1eCnlXKctyCq37ronmwpm9Zi9mrdHbJcVZg11EozGTnyQ3DiYLQy/aml/KZv7Mmd2yE7t10L0f2acTTu3fWXtfWZL7veR1Q5scRodWlYYx5oROeOL6+oz7Kgq5XnmssTWeVng8kWBaX//hJPN4Eq/58X+X4ofPLs28owP4tdxy/uAMe2ZP4LxriOgCIlpLRBuIaLrX57NC9HDJltJwCM3CAyh6QXzlyQV4er6+clDYoZC/83OpGYVfpcCCCPe7dmJPt4N7O5gNvndZzP56JvvVcYtkZE7glaKcCvmjLTG89O2JOKVfZ8y8+Sy89r0zLffN1bsm08BgtMlHEwkt3P/7kwahS4aKScXuYcMHMaNXnhsESpMnohCAvwCYCmAYgKuJqCCuJfyH7VZlX3TBjNKIgo82HdBCtzOlbg3pzDXWQr6msgQ/SkbRFnunz8S0B+fijzPXFeTc3AbemudCFr8H5SZC3mox9fbP5t+dNU0+h4HixG5VGNbL2i1XIUIuP4tRUzeSYHqvsm0HmrJajPVCeLrJkSZ1dpjt+ooTuBKYb8Uup5CXGiQRjQdwJ2Ps/OT7WwGAMfZrs/079xvKJv/PE560ZeXOwzjWEseYvp2y7mArdx7BsZYYqkrDGNarGgnGTKs5jezTEeWREA43RbFmt5qDuyREGGOTMGrnoSZsO9iEcXU1bTo18cebVVe70wQThl8cb4lhxc4jeZ8/nmBYuPUgSkJKmpAb3qva1N584Hgr1u89hk4VEQzubh4YlQneX4b1rLb15uK/ce9O5ehT46wk5dLth1BREsYgg7kpE/y3sGNcXU3ac8B//zW7j+Jwk3Wa4cE9qtDJ5SyMbvZBfiz+TLsJlx81FREM6tYh74RnAPDcjRMWMcZM7XdeD6e9AYi2jO3JbRpEdAMRLSSihdGod7mnGVNDzXPRIPgtaE4uNlmNi9wnXnfPMtzANizXiwa31BR+nCyLNQHILfKRE8rSbzobxYyy3J/jpC12k89B3e0HlbVJJcgtvPI59+Lx5EL9YGPUMh+Sm3jtJ2/2G+nuBmPsEQCPAEB9fT171uCx4AaLth7EZQ99iMnDuuLRazMvVhm5+pGPMG/TftRUlODZb47HvmMtqP/VzLT9Hv7yKehWXYb5mw/giofnAVBttnbX9Pjczfjlq6vw2HXj0NGn/NJeUDd9BgDYXqtXLNhyAJf/bR4Uyu/8DUdbMO7umRjQtQNW7TqC4b2qkWDA6l1HcPfnT8bovumBLzOW7cJN//kE4wd2wV+vOSWn867ceRgXPTAXP5h0Ei4Y0cNyP/4bf25UL9x64VBHx66bPgNN0QT+/MUxWLD5YJorsBW7Djdh/K9npW0/uXdHnDe0G/44cz0eumYsTr3nHd3n4u/P22uFm33laHMUJ9/5lmvH5W3/65dOQe9OzmZNTmGMof+trwEADja2utLe5260/sxrTX47gL7C+z4Adnp8zjQue+hDAOlphJ1izFhpFcQQTh5fPE2mSNuUfa54bPIb9h41LYRSrPD7kW1udiP8HvD7fXLvjpqnjdV95BpkPlPuihJV12qKOvvNc+krVz38EW76zyemZSrNsLKvP3LtKeiVFHrZFDv3Gjefn3kbU6UMc42Ut0PsK3489l4L+QUABhFRfyIqAXAVgFc8PqcluWa/MwoPKy8Ovp94EzMV8PY7j4UTJt0/B9c+Mb/QzXAMf8DzfWD4PdAXjrH/Dhfy+aynpBZerfuA6N2Vi1fGpn2qWYA57GavrzBPDxxSSFuEzrQ46ydueqocakzlxg+COdVTIc8YiwH4DoA3AawG8BxjbKWX57Rjj0WyoUwYH2CrDsW1cnHvTLbCYtTkgex8zgvt/slr8OZrl+X3QKH0wZpZWP7HJ3PGfGViXc7n5S6UdumsxcpmefUVh0Lr3tfX6N6L3kV8RtxikUPeKVv2uWePFn+TYXe8gbrpM/D++gabb1gjPu5dO2TvjZcNfqQh99yPiTH2GpYekYsAACAASURBVGPsJMbYQMbY3V6fz455m9IryjvBeCOy0eQzPQirdqleIe+uza1Duk0uAvtYgU07vNBHvmMNH7z53QspqRwmVsfuVl2GLfdelFfJPa7JN9v4yYsDWD4RrLkOyL+6ZATCCqG6LOW8cOED79t+J9PAd/bv382pLWaIihdPGvflx+fnOCA6i3NxA2N9AS8obmfVosGZuYb71IrBUHYPLgDMT7pqvbPaPOGV3+TyTNz3VmH84zniA75DqMiVLVwg8DE6RISfXDAYXSpL0iJe3aQkrCCskGlGS44om62qVzkh11nAZWP7YMM9F6IsEnLsofbzzw33rfSl1W+Sb4Cc12zd34iVO80ri7mFFPIOMA7mVguvfD/RvJNpccoNH1k3ETXGfcdaMg5SAPDkh1s8bFFmxPWMiffOykpb3dhwDIeTaRH4feX3RFEIZw6qxaLbJ6PSxeIeZsQSzNZE9tzCbbp9cyXXr4prS9k4MHitCXOsZje5TFz8eCSvHd9Pe+11rVcp5B3Abzp/QKwWSblwEDv2cRs7KwBUJO2x2aRC9hJRyNf/aia+8dRCx9/164E2YtROX1vuvKboefe9h88++L7uOJq5xucBmAfgmPHvjz/VXueT2M4Nf/JIFrEmTuox5MqB46147P1NYIxZRo/mYtryY4npFxeP0F57EVUr0q6EfK6LKFwz5x0pU1oDMU9Opg5zW7JW5JQiyUhpbK9VulkzCiHiDx5v1ULQOYeaWi32NmfbAdXEwwdvzVxTJAOvkSnDrX3pM+HGAn82mryb4+Thxii+/8xi7DnSDAC45b9L8asZq7Fk2yHL2U0u1+t3JS2v+1nghbx4k62STGWCC3k+nc+U7L9jeQSzf3y2o2PzRE7Fkp8snw5eCMvTmF++jdtf1jtsNbbkliwspcmnL6AHBTcEmJW7qJmwclNgzly9By8v2Yk/vbMeADSf/2icWT6TuWTg7FShBiX+6SpnRdvzxetuFnghLy685Dor4jeB+wVzl70Km/SyTs0voSJzoczPO684hGKtIQnd4aYodh5qwh9nrkPd9Bnag//E3M26/bg2eMagrgDU6j9BI5dwjMnD9LNMKwuMWTpvJ/EDJ//8TVzwxzkZ9+MupnuTmnzK84lZenjlYq7hv1HPju5Gulrh9bMfyPJ/rbEE5m8+gDMGddVqNQLA3qP5+clzIcC1hrumDcfzi7ab2lIzBUFxuD0uWiTBUPloXsUSHGP0/pj6xznYebhZG3h3H2lGr07l+MWrq3T78YdtXF1nbLh7qhbB3JYZUFuJTUJ+lGzv79+vH4eJJ3bVbbMS3M9/a0LaNiezoaMtMS2hnx3/XaQWBJnJPdGShz7eGsNXnzRfO8pFk+e/kV+33+tHv+33YhN++8YafOnxj7H404NYtycVRNKlMjebPO+ncU3Iq3dlRO+OePab402z3jn18Cg2Td5pRCTH+BD5GRhl5ZVg3L7zsKr5ceE04d5ZaImlm3T4IB5SKBACHgD+96aJuvfTX1yWleA7tX/ntEGz1mJta2jP9JTHH+UYm2LGsu0pV8OjzVHN/fioTb3YXDyR3EhVkQ1eR7sHoycb2NCgCvaDja06e9eFJ+e+YCXCOw5Pk/DE9ePSaltyWde92n5g4fU9My3m+oVxeptpQmJ8iNbtOYZXlvqTnujvH2w23W4mwAG97dPMtZUvrBeLp5MbGEsEfrTpQMZYAnEQMHMdrakswW+/MNLR+Y0DrlvJvl5dtsvRfrkoT/wR8Cv1t9cLvYEU8vw3IxAG91BzfN81bXjOI3NNhbo4yrNEpnKcqD9fZWk4Ldc4X7y5sr4v7AglB4pcQ7DdxtjhMq38Gx+iLz3+Mb739GJLQesmvOwf571bzgZgHZsgPrRmMxa3Ep3lwsg+anF3My3b7dlRJrdYJ6bDTM4HVoQUwpdPT/mIbz/YaLO3NeJPYickc6qKpeUjyvqrOeG1lTOQQl6DUm59I3p3zPkwP71gCADVVgukOridxldVFsHqX1yAH2Sod8mP8e7aBry3rvCCPlshb5xqNiTXPf4j+HV7wbGWGJZtP6Tbxgdhq1QS4hhvtiBnzELpJ1OSi5tm7RJnS24Inkz2bycCXBzELz+lD7599kDT/TobygAayxGe8ZvZGc8FpA90YuCYXXtz0eSNOYy8RpprckC8rev2qB06nyl4eUkII/t01DpnzKEwKC8JZVyAFdvFvQYKiVHGZApmsXqIHn5vk1tNMuXGfy7Chxv19l6e6MvJAnDMZJ9oIvPg7RV8VmgmsESTxwfTz/W8LU6EfLMwkF44sid+klSEjBiVBKLcPLiM/eyFT7Zrr40mw1F9O+GqceoMOjdNXv3vm7lGavLZw0d9AnDPa2o2vXy1M4Uo5V2jCYP8fz7xIZi9tvD5azJp8it3Hsb6PSlN0GphqyzibdeauyE9SKskpIAIaLFIxSBeW9Sk3XGDGc5P+MBiptWJbsA9qnMvSO8Ubq755cXWcSVi0XG74KiRhhl0iAgvL9mRdZvsFlCNnykEnDlIdX/NxRTC5Ydf3SCfhHNOCKSQ5+jyuuc5KocV0uylXAt0Q+MT2/ja8t2YtWYPALWjzVy1J6MnxP1vq77fbnnnGA9jrNN50QNzMfkPc7R2WZ3XS++UT/eb23GJCIwBM5abL8qJ2qdZQisnZjiv4EqI2e8pavJ+eHxsTf6+TTZ5i04foJouzx/eHRMGdrHc755LT0aP6jI8+MUxAFSlwW5/K+w0cuOsjJByf8wt4lX9758mL4V8VjDGPFmtDimkaVRe2m73JnPev7R4B77+1EL86+Ottvv/7b2NANxLcuS0w/F6t/w3qe+nzzbopaDctO9Y2jZx4dtJ3UyzxcVC2uT572XmZcXNT049WvLl/rfXAgDeWrnHcp8JA7ti/d1T8fCX620Hnu7VZfjof87D1BE9cfWpJ+DBL47F8F7Zr4/ZafJGQa4QabOxtrDw6mbBEzMCJ+TH3f0OPtig2mrdvEclYUUTaFFN43P/5+MP9M6km9vuw/Z2ei4cnGSLdILTZ4IPKk3J9LjGtQfjgpubmAkVJ948Jwi5u800PNFP3m/4zMdOk8+1fOWphjiOMwd1tdhTf75MibOySawVUgi/vvTknFM2izMvY98yCsn6us6aJp+Pn7zXmjxPm+B1jEzghPy+Y6moVvEe5avcl4ZD2nRf86f2QOO7I5mHhc9AMwkcnoP83XXu2POdaj58MOLnn2+I+jUreu0WZr+IVZyBKPzFh8ksz3jcxbWWbLErA8mvzWkedyP//NqpmPvTc3D5KX0wpEeVdo/fX9+AGSb+5gu2qJ4r2WSb9BpRWB84rk9AJ85iH7u2Hj+ecpKQVLB4hTx375ZCPg/EXCqVpdZ5ZpxQGlE0gRH1wXYbd9jRuDAtCeV3fYCa5c9pkQX+YHFhbxyMvOy42RQZP3g8taYgXlvRafJcyNt41+SakrY0HEKfmgr87vJRqC6PaOf48uPzcdN/PrH8XlWZN1lPvjqxf9bfsetP4kA+pGcVwiFFu4c5mWuS3cRrIc/v+fq9mVM65EOghfz8zSkXu35dKvM6VmlY0YJs4gmGkEKeLoIlHAqcrh3UqWu+FXCi8QRG/eIt3PbSCst9eNY/QBU8R5uj+O5/FgPQ1wAFvBXyv5qxOm3btFG9TPfVeYEIbTRq/mJO8kIsvHKB8kAyw6JIa1y9hlw1eZEQUUbBd0V9HwDAj6cMzvt8ZnSsiOCik3tm9R2u9Izqk27PFwugVydjJUIuaPKeZ4dMnuCPM9e7Zm41PY9XByaiO4loBxEtSf5d6NW5rHhg1gbUdanAxaPNBUA2cHMNYwwPzt7gqhAz2khfX75LE06ZhDz/PN/kYFxbnL8lPdkaH3AOHEtNk1vjCTw1byt2J337jQLIiS00kWBoONqCg8db8dC7Gx1Hdl4yJv1+WmUEFafyI4QFP+P9S7DCavLPJ5Nvvbg43b2wNaZPo5EPIYUy9t2qsgg6lIbRv2t+ipEdThP4cbhNvjScfp9FDzCexoEf32jacYKW1sDjfiCaBfMtim6H15r8Hxhjo5N/r3l8rjSG9axGUzTuSuWV0rCClmjck5XwnxoCSb7170/weDINbiZtQhPyeXrX2F0XH0BEG21rLKGbPVx+Sh/dd5wMgg/MWo9xd8/EV/+xAL95Yw0+3nwAG/ame84YsUqQJTJv43589ckFOg2JCWFyd7y8QreQmWBMEySF0OQnnGjtVsh/f+NsKRdCCiFTrFMsnvDcw8jsJ7Yb5PlahZkSIgZGcQ4lU158+9/W5igrtCyUHqvy4tKPlxlcA22uOdwUxZ4jLZZ+1dlQElYQTSQ8cc+0O2YmuyB3FctbyNt0Mi4oE4aFS1GQjxd8n7t2KHUU4PFm0kVv8z7V5fGhdzdi0v3vpaUrSGurgwHkmsc+wqw1ezFHSBUhtnfN7qO6ByueYAXV5E+tS89kyonmaZMXaYrGsXTbIWw7YP1MRBPM84HOrF/bKQbZzpzzmWnHfXKhFB+RP89KN9O5hddC/jtEtIyIniAi07LtRHQDES0kooUNDe7mbuHZ9txY2FCIkEh4U8GpX2fraXEmbSLlX+2dJv/MArWItPjgxBN6L5deQnbBsEKmwUZGVu86AiAlvLhw50LfCicPMN/lvrfXadvsMn0y5s9aixWdKqxdTvlg5IZNnntBnfnb2Zb7xOPMcw8jMyFv6wvv4METKznlo4zxZnjdD8Q2PjXPPh4mH/K6k0Q0k4hWmPxdDOAhAAMBjAawC8B9ZsdgjD3CGKtnjNXX1uZXicdqumfMVpgLIUXtaF5o8h0rIthy70Wmn2XqZ/xhyVeTtxskeKbAmE7IMzz1UapjhhUFM753Bj669bykScD578TNJtxXPNO1GIXBpKHd8N3zBgGw9wG3GxzW7D6CjzfvL1hef+4/blaNit+bXP3ksyWa8N5c88mnB9O22f323CPox1NSCf/OG9JNt484UNo9p4s/PYhDjda2euaTJu9XJbW8fKQYY5Oc7EdEjwJ4NZ9zOYEXZPaCEKkLVlc+/JFn5zA9b4aeZixNmCt2Wi4vgyY+OIwxze4JqIuCJ3ZTFzatFvfiCYZ4gllqpHxWcsgwKB843orG1hj61KjBTE9+uAUAsOzOKYgoipaYDADOOLEr3l+/Dxee3AOvLd+tO45dtr9vPLVIF2NRKN5b14AjzVEcbY6hc0UJyktCmleXG+YaJ8TizPNzmc3W/vruBtxyvnmiM96fqoT8+OUlIYSVVE6piYLJ0Kp0H2MMn//rhxjWsxp//8o4dKksSUvBwc2SXrtQ9u2sb2PD0Za00pVu4KV3jegj9XkA1r55LtAcjeOs31lPQfOF276X7zicYU93ydTReIfkQj6RYLjvrbXYe1QfKTt3/T5bDdnOJv+7N9cm9xE0eYOmJD4oVkJ+yO2v46SfvY61hlS3fPDgA83dr+ldJMf/+h1dSlqeUriqNKwT8PzcQMojRcROUywGAc95dekuTLx3Fr7y5HwAgibvU3BSLJEoyOLzX2ZvtPyMC3Ix8V1ZJKQtRvfsWKbrg6f274zhvarTomN5H1i16whOu+edtBKQgJC7xuPfwGgOuvFfizw5j5e95rdEtJyIlgE4B8APPTwX3lntbQbHQlWDy9TRuLDlAnzxtoP486wN+MEzS7R9Pvn0IL70+Mf4/VtrLY/jZCYgavJGgSkKBSshz4X4Wyv1GjY/rlVqAmMRkNF9O6FXxzJTmykX8sZBDlAFhdfZMd1gaE81EvKjTar9PBWJ64/gjcZZ0ZU/jGtCPjWoR0KKNvCZKUPj6jqnKS/GpGsvmbis+p27htOQYw3qTHh2JxljX2aMncwYG8kYm8YYc1avK0fsIvcmDe1u+ZlTvB7VAWCiiRtdptPyzs+F/GUPzQMAzX8dSHWeTTaJu5zkEBdt4Uabp07Ik70vtqKQbv1kTzIpmyjM7UoIJhJMq6hlhLdDrAfKiSdYmp/11BH6kpA/ucCbAKBsMK45uFnE4jeXnZz5/PGEKz75bsJ/g3JByCuU8ps3WydWkllJRYz9ysw/3a+0BgBw8+TUGoNXg3hxDdce0acm/7qSXvvMAsADV41J25ZpHTBl6tB3VrG4sZOFJCfVaYzeNSJGc00mN0ezNLaiOel7Ty+2TF8QSzDL+2GXC161Neu/d/Zg/UJnpAB5a4yIv8Pv3lyj5QdyQ9FwYvON+eBCmS2akBfMc0QpE5bZepxCemVk24FGLNyiX/A1m8H6mWpYfKa8UiS9SU5RZIjZB3PFD99pM7e1TFGgXAM32tsTOq1b/W93DWY2bCN2Cb5E4amaa6wHDYUIzybdMu3Ye7QF/U0KSccZs7wWOytDLJFI+57xNy8GDVYM4PrL7I2oSv4GbnRBJwuq0XiiIEna7OBKiKjJE8g2QExRSKck2bmNiviV1kA8F+CdIllcdzIPeA4XM66fUJf38Y3C4Rc2VXNyPoeJgMnk0sc7iVEj2S+EczuZfmaryRuFvCgUjFGVt764HHXTZ+jaY1eQgrN1v7l5KR63E/LWXXrdnmOaaUhrt+E3L4bMi6t2HtG9P5qc0bihaDgR8rE4K0hOfcBaqeF9T2z/Z06qtV2MJnLmX29cC2IOlCK3EIW8V5p84Xu0S0wb1dt0e+fKEld+POMNdyPE3IjZFNlKyH/hoQ9x3n3vCjZ5686cCu6wPrdRaD96bT2++ZkBlm0xzhzsNPmn5+uLepdFQo6Kk/z8lZVp56+bPgNvrNxtKcyzNTMY76Nfbop2iAFcIm6YD5x46EQT/i68isnKjIMwh/c98TmcNKy71reNayuAqhk7yYf08SZ9qgS/XCgBvTnWqzGl8D3aJRgy38x8MN5wcZXfq3MA1hG2C7cexMaG45rG3BpPWHZoJ51W9JMPKYTJw7rj1qlD0btTOT6XzPAoakXGmYPo6ZIpCVYiwdDdQa3SrYZ0FOLAYiWDsh3QxxoqWvkVcJQLbggds0Rua3brZw6xeAIRj7XYGd87A2GF8OH0c/GXa8ZqVa+OtZgHLlpV7eJR02ZZHBUiROMsY0GZa5+Yrz+Xj941YjqLMSd4U4OheHt0llhphrlkoTPDqCGaZcNz+xxA5ukmv+5oLGEpWJ24hOl84IXjlJeEhHquKSEbtfG5D2cQ8i2xOLrYmNesEEv2ZaPJ2wnujuUR/HBSysOhGDR5K9wQOuUmysnVj3ykRYDOWrMHK3ce8dxcM7xXR2y450ItHUbnZLTq8ZaUQJ6/+QBW7lS9pPhCflghXDq2N/5w5Sjd8WavTU+Jwvv9tD9/kFXb/EprAADnDOmGpXdMARFQY5PaIh+Kt0dniSgMzxzUFTedM9DV4xs1RGMQjhfnADLn4NDMNfGEpUeLk+AOq7QGISLNXi/uYudXzzX5yfe/h0v+kv6AtcQSeNWkIpEZPP8QoNfWrGSQmR2V5xg38sntk1EaDukEWiGSk3FEdzoz3Gib2TEONka1HP1ffXIhAG8LsZvB3Wd/88YabdsVD8/DRQ/MBaB3I73/itH4/Jg+6Qcx8Nd31eCqtXsy564SlRLGmK8+8h0rIqbunm4RGO8aUeac0q9Gi4p0C+PKt5lG5AWZbNdiMJTVgOBk4dVKyCsKab+tfuHVul1cyK+3SBv81so9jh48QK1WxRF9mq28P8w0+Y7lYdOIVi7cFSoOIT9+YBfgbevP3dAsrfpA2hqLz78Dzwn/4cb9pp+nNHlvBp9oPIGQoj7TCRvvLa8wunu6emxPjloARHt0PMHyTthlxHjTvYqcvP+KUboF0kzrkzyroDH1b6+OKZu3Mz958xOFhUVUsRPa/b6q9m/dcKcCHtAv6IpJrayedbPZilWGR+4Tr4/Wddw017ETK24JHavjGGW/35o8N99ZXaZWKc0jM5Ko5CSYP6Ya/fmZZylTAiPk9YE6TJv+XTY287TOCUbh0bcmf997My4d2wd1QqlCp1kRWw02ef4yFk/gpy8sB5CHJs/S2/LacmtzS2s8gZUGN8BcESMUvy+kashGk6+2qFWqafLCd/zwqLDCrmCKW4ql1eUt3HJQd36/g6F49s3vn2dustJy/VtcgF32USeIM9OEz+Yazvvr93ly3OAIeUHL7FQRwdwN6g/mVpFcY6evqfRmkQRIz/bohJZYQhflyoW2uCBlp51Y5c0IK6RpUWK7jPlkRNzsrH+etcF0u1VNTDNN1WqRnN9TY96dQmG3ZuLW4GMlJHccatKZtPz2k+f3wOon4M+31f35+eeGpW17+hunOz6/mOMmkWAFHezdJjBCXrRdf2Vif3RLhm9/dmR2BYOt8POmi3LdaV72aDyBqX96X3vPF0bFWXeJzYMrml90Pu/Cwqso5Os8rP/phIVb0/ORA+YaviiwqoQIWj7oicLVj/QVVtid22tzDaDv435HvPJ2RRMMj8zZqFuLAYB43LzI+oa7p+K1752JE7tVpR3ztP7W1bYANdEdR3Qk+HDjfi2VhJ8M7VntyXGDs/CalD8fTj8XkZCidZpuVZn9sZ2wwscUw6LPv5m1Zs+R9AyLrfEEjrWka/IlISGhk80DLtrQjdGrZt41blWXt/KpnzKsO95atSeH45mfg2PmFVQsmrzdgO4kQtgJdrM58bfz+3cgIigEfLBhHxZtPZiWYI73T2MfDocUDOtlLhwzxUxMHtYdS7ap1chEc41bpsZsmHPLOeicg1uxEwKjyR9viWFA10rN79ZJlGc2HLNIluUF4rNu5l3z2zfSUwYb9+OdVtTK7XzARUErar6lEQXN0VSueo6dkDfLpmmFlTkq19zpov/8NaedAECvoZqtPWxqSNmiCynkpyWDzsxwy/HC7vrsUkL4QVhR0JTUoLm3DYfPInNdK/it4JrJEQMa8y2fmS8ndKlAB5M8TW4QCCGfSDDMWrMXm4RqM6kkQ+48tH4VbAD0D/THmw9g24FGrBc8Usw8e4xjAa/CJF6/nVupKOTF4KGKkpCmRYqa5rGk/f+68f3wzA1626fdYHLt+H6691ayK9eHWfyeloZW+A3unJaec2ipSVriQuBFFLURcVAd0kNv4vDKhc8pimLdhnyLrHOfeRHxOXLbG6+YCIaQN+sYWrpQd87hZ7i7+CAu2XYIk//wHib/YY62zWzcEn8DHrq+81CTLvGYnaeOKMDPG5qqnVkeCWvalfh9njTr5smDcfoAveb+E4sSbgBQXaYfaKzkSq5apSgE+EMs3rprx9elfUfMX1NYMQe8+YOzPD0+V1amjuiB31+ujxp9/P3Nnp47E2FF0fqYsV/wwvDZCvlJQ7tZflYWLh5N3ksCIeTNbJmaJu9Ssdywzuzh7ZTeKIubDYUNjJ93ry7VCWC+aHTPa6s1jVv9nrkI+9qTC/CckPpXrLNZURLCjkNNmPbgXPzsf1MVHHmud7PIXysbKWA96Bq9I3LV5EUhwGcUmRbNRe+bQvtUDO6RvoDoJqXhEJbdOQUPfnFs2oxwhuAWWwil/lhLTAugMwpdzbsmy5m53aAgzpwy1T9oywRCyJt1yJF91JVzNwqGAMCV4/pqr1f94gJXjmlFpmRrxrwxCWb+G7y+YjdueX6Z9r4lltC0cpF31uzVdXLRNMWFuHEhrLE1jrBCWZuxrMxnxtnAZwzFPBwfX3jNB4pMAsuLjKLFTHVZBCGFitpN8OPN+syQqukx+wR0dl5C4iAnPlMDayt1mTHbOoHo3WZmiG+eNQBv/OBMjOrrTmY3sfq710msRmdoc7NByIeIdLOZr5/RX3stLmA99O5GDL3jjYznN6YNtiIXG7LxO1zOhBXS2fZP6FyBIT2qUFORe3oK7st/rNV+0dzP9ZZiwu+oznzItVqV3aAgpiZpNUS8FnIB3m3y6t1EdDkRrSSiBBHVGz67lYg2ENFaIjo/v2baY2aGUBTCkB7e+J16zR+uHI2Xb5po+fm63foAr0vG9Nb9Bl06ZC7xZofOhdJGEORitrJ6dkIK4fQBXYQqSIQ1u4/iYKN56lmr+AexJzw1bwsAYEYyGZqVMB8/0Lk3UJBwGmhXCEQngYajLWpd3xwEr10XLdV516R+C7MKYm2ZfFWYFQAuBTBH3EhEwwBcBWA4gAsA/JWIPHMd4Iq8Vy5IflNREracgRxqbNXlfikJKSgNKzqTRDbmBzMXTVF422lC2cxouDujEX50bWBJbsj0kE0YmF0Y+ye3T8bCn00y/ezyU9xJfeEl3RzUZs0WOxHvdX2GTJw/vLv2+oqH59nW9bXDrj69zlwjaPKtsURR1xXIlryuhDG2mjGW7rQNXAzgGcZYC2NsM4ANAE7N51x2cEH1oyn2qVrbOokEw4NCmP/Xz+iPF789IU0gZpOr3bhoff7w7rppvJ2szUbImxWrAFImA2PiqUyHdvK8i79L58qSNM8eYxuKlT415fjLNWNdP6546wcYIpgLreSL59+87zjiOWrydllcSy28a5pa456kEi8UXg1XvQGIlZq3J7elQUQ3ENFCIlrY0JCe+N8J3FQRpCkWAPyvwWQTZwyrhSo+t144FCN6d9QEMU/SdN7Q7nCKWCxkSI8qPPxlndXNNXPNFfV9oRAwdYS5mYXbW/kRjYuComYHAP0t0iqc0LkCJ3SuwBPX1xf1wmI2zP3puRhXZx+inwvFbK4xyuYjzdGMGVnNj2Mn5M395JtjCZR6lGW2EGS0bxDRTADpBRSB2xhjL1t9zWSb6a/NGHsEwCMAUF9fn1Ov49posWtk2WKcMsYTTLeND2rcpBJWCCd175Ax133d9BlYedf5qCwN6/zozQZJM3PNOYNrMXttQ1bpaAd1r8KmX19k+Tk/t6bZG8PXhXWCx6+rT/PG4ZRFQpjzk3Ms254J7pVVDCy/c4qnfVpMwVzoQCgjxpJ9L36yI6fjWF3XpWN7o1Iw78YSTEvX3RpL+FYvwg8yCnnGmLkh057tAPoK7/sA2GmxWPyRrQAAEMBJREFUb95o1dUDJuSN3l+xBDMtdMw11lgye56TGc2nBxoxtGc1PtiQKtJg5r1gPFa/LhWamcYuNa4dlSUhHDe4ckYMF2vUwvt2TqV2tvPD1x/DfPusH33GUssvJk+bKgvzklvUVpXiw+nn4t7X12DxNvOEb4XCLQ82q7Hr/itG60qDRuMJDLrtde29H9HHfuFVj34FwFVEVEpE/QEMAjA/w3dyJlUazKszFAbjoBWPM6zalZ48ie+3dX8jjjSle6OYCW+u4cxZnzKRmWm+xjbc8dlhOSUOE+kuFDThR+c2eX464+By4cmpyaRTs5yVIB9Q26HgWTSLhV6dyhEJKUgYAj4Lbcpx6/x2Jh7xHMa0BmVFNNjnS74ulJ8nou0AxgOYQURvAgBjbCWA5wCsAvAGgJsYY57l7tTK2wVMyhun6l9/aoHpfvy6Pz3QiJ2H0zNUmuUG5/07LtjkzWZC+dS27dvZPBDt0WtTdn/RTx5ICWajIBcFdmnIWRu4ScfrCGW3ubK+L07pV+Pb+Yyl5y46uSe+e94g385vhl0N4WywGyysXCgBb2o4F4q8fA4ZYy8BeMnis7sB3J3P8Z3CtZCgLLRxjPbEBVvMp9SZxraIouAb5w7QFeDgx45msMkb5aPTPOMb77kQBGD/8da0EPWBtR201x3LS7DvWIsm7PlimDHMXLy3VRaVnox8ZWIdXlq8w1GqgMqSkG6GUUh+84WRvp4vpJCur3nhyZMtYj3fblWl2GtR1CYTdvOBDqVhzLz5M5h0/3tpfVSaa4qMlHdNgRviMk4z45kJ5xe/PUF7HQ4RfjRlsO7z11fsBqD3rjEV8oZtjUL0qJ2wDSkERSHUVpVq6Z/NeP7G8bhr2nDNna0+6UVi1L7FscWLGdvyO8/HzB9+xvXjtgWIKCfPFa/o2bEMzcLC66DuHWz2tifTgvLA2koQqak9RKSQLzK4d43Xmnx9vxqM6O1fFK3dlFWsXWvmgTFGCKYy84J5KJl6tbNQxtCJd83W/Y3aa6uyetlQ17US102o097/7gsj8b83TUS3qjL8VdAo4wmGx6+rx08vsM5waSSqVRPK3M2V5KDUHlFI9ScnAr537okFa0fvpDKgEOk0ead1js3IlF2SiBAJKdhgKC4fJCEfiBBR5pOQf/5bEzLv5CKDu1ubGX520VDttZktnYhApNreIzbCa2jP1DnMarMaj91JyCWTT14ZK8oiIS13jzgAtcYSOG9o96xiAHj2TrP8+5IUIUXNfcRYYd2QX/3uGdh7tAXfeGqhroawcVE4G0TvMStKQkqay2a7cqFsC8QDapOvLA2jpiJimr9FvFZRhl8vaMVKMnGZMZpUpNUu7hvptvFJgpC9TRhosuXD6edmTL8gmmyMC2NOGNRNneZfZ5JDXpJCoVQJxkI+QjWVJaipLEFIIV3lMad1ju2ws+tHQoRjLfpzBEkxCMSVBNUmDwAvfnsizjZJu0sWNmoxfQDfbPRBF4llmM6uNrhsip46+Uxpe3Uqz5hITSwckktRh27VZdhy70WYGqC0sV5AJLohF15RUgg6TT4fcw3nfy4cir9fP870MzOf/CCZawIhFlNaSOE7qNv071qJa07rl7ZdrGojTi3FDst/DzMXSk6mqvRvrtytex9RFJzaX10ctcsL4gaif399nX8uhe2NEFGqUHYRPEIhhTTzCVF+0bg8sdvpA7rgnCHmVaLMhHyQzDWBEPJBjXjlmAUziZGZYoZCXQZJzQfd+jbvP5aK+psyLN3e3a1K71aoKKSdw2uPDD441VREXFnklZijKKR5chWDoqQQaea5kpCi8wDLFq6QVJRm13+ClLsmEFeiedcE4mrSMXq9fNGQtldcoBS1Ej715kL5S6enp/sVzSClJtrLFfWqF88vLxmB//vOGQCAwd1VD6NOHiy8ivDBLWiJ54oNcUZWBDIeIYW0Sk0lYSUvTf73l4/CWz88yzIDKQDsONSUti1ImnwgFl4H1lbiP18/DUN6ts0iIZkwavLG96L9MGwi5Pm2oSa/z+y1Kf9gM/PLdRPqcM6QbujXJZUGYPrUIZg0rBtG9O6YzWVkDbfJF4OdOMjsF3K4uFUTOR9CCmnuwyWhVHHv7+cQhVsWCeEkGy81u+8FhUAI+aqyCCacmF0RibaE0c/dqNmKIdglpuYa9YW4gMVr324/mNJizBa4iEgn4AFVu8q2aEcu8HbnWtRb4oxjLakAt2JIAaEQaQuvkZAirBf41zavS3z6SXCuJMDY5XEBbDR5hZtr1G2iG6Io3DnFVrGez9Lba5CSX4iR1cWQhVPs7yVhRWtfgOSurwRCkw86Rk124VZ9DhsxY56ogRxK+tfzBcz6DEmvii2neFxzjZVC3kvEdZliKHsnOlBEQimfeTcWhf989RgMqLXPQDppqLkXTlul8HdUkpE0IWcQxqL2bjbd5t41o/p2wvq7p2KgRScvNiHPyxh++fR0F1KJe4iafDGYKUQHCtFc48Zg/7lRvTC8l/Va0rRRvfDYdeb+9G0Vqcm3AYx+7nZTajNXw5DhoTm1fxccboql7edG0ImbVJdFsOVe62pSEncQcyQVo7mGe9r44SL9wNVjPD+H3xT+jkoyku5dY33bzPJgf7TpgO79ocZW7DvWoltwA4pPk5f4Q52wsF4UmrzOXKNo6bClk1VuFP6OShxgnXbXiJl/72FDtSieZvjFT7brtucTdCJpu3xXyDxZDN41oiYfCZGw8Fr4trVFpJBvE1gX0DBSaRLZd4JQH1XkjpdX6t5LTb59MkjwIy8GTV40y4QVRYuslkI+Nwp/RyUZyUb2VpakL7NkmuaedZKaAK3ITPKSAlAMQl50mRVfOy2iky1WjghBofB3VJIRo+y10+TFHB2ZUvlqx+euitLo2e4pBmVZbINoPXpmwTZPzlcMA5uX5FvI+3IiWklECSKqF7bXEVETES1J/v0t/6a2X/p31WsadrJY1OSfv1EtcpLpueVVeIKa+0eSGZ57vxgC4kSzjKjQNGXImCoxJ9/HegWASwHMMflsI2NsdPLvxjzP064xahp2mry48Gplw7zlfH29V17AmpttJO2PHsk+YFdy0i/E/v3hxlRlJ+bRmlExZN70kryEPGNsNWNsrVuNkTjDbkqtGDwTzJg6oofufd+acsy79Vx86zMDXWmfpO3RvVoV8sUQ8SrWERYznXo1yZgwsIs3By4SvLyj/YloMRG9R0RnWu1ERDcQ0UIiWtjQ0OBhc9o2j15bL+SNd6Z5cA3FqKkYZwLxBEPPjuWB12gk1tw5bTjumja8KAQed/mdNqoXugrVw7zy/rr61PQU3EEio5AnoplEtMLk72Kbr+0CcAJjbAyAmwH8h4hM8wAzxh5hjNUzxupra6W5wIrJw7rjwS+OBQD062LuEmmEy2yj6DaacYrBDispLB1Kw7huQl1RDPSXjO4FAOhYHtHNWsW6Ce4S7P6fMa0BY2xStgdljLUAaEm+XkREGwGcBGBh1i2UaJzavzMevbYeZ53kLM2vleJjfI6LLZ2BpH0j1vYVB52BycVht2nxyDWzWPAkdw0R1QI4wBiLE9EAAIMAbPLiXO2NySYl+gBg2Z1TTIw4SeFt+CBdkw92J5e0LZqSWScPNrbqNPnbLhzqyfm88r8vFvJ1ofw8EW0HMB7ADCJ6M/nRWQCWEdFSAM8DuJExdsDqOJL8qS6LoMpQ4qymQp3eThmmX2g1+sNLTV5STDyz4FMAwKvLdukUkl6dyj05X21yravSJO9TEMhLk2eMvQTgJZPtLwB4IZ9jS/KnS4dSLLhtUpot02h3jcqcNZIiQlQ6/ChH2KemAvP/5zzdIm+QKLy/lMRTaqtK08wzxvdSk5cUEzyOo6rMv0zo3arLAluBTAr5dog010iKmQuScRyl4ZDmNnn7Z4cVskltGink2yFkuOveuaZJJNlTGlJt4/FEQnNuHNXHupqTxB5ZGaodImryv/vCSHx2ZK8CtkYi0VNdHsbXz+iPS8b0xoOzNgDQF6uXZIcU8u0Q0SZ/eX3fArZEIkmHiPCzpHnmN18YiXOHdMOI3lKTzxVprmmHFEFQo0TiiI7lEVwxTioi+SCFfDtE5o2XSNoPUsi3Q+xSFUskkmAhhXw7hPsDhwPqFyyRSFLIhdd2ys8uGoozB8msnxJJ0JFCvp3y9TMHFLoJEonEB6S5RiKRSAKMFPISiUQSYKSQl0gkkgAjhbxEIpEEGCnkJRKJJMBIIS+RSCQBRgp5iUQiCTBSyEskEkmAIcaKpyoQETUA2JrHIboC2OdScwpJUK4DkNdSrATlWoJyHUB+19KPMWYawl5UQj5fiGghY6y+0O3Il6BcByCvpVgJyrUE5ToA765FmmskEokkwEghL5FIJAEmaEL+kUI3wCWCch2AvJZiJSjXEpTrADy6lkDZ5CUSiUSiJ2iavEQikUgEpJCXSCSSABMIIU9EFxDRWiLaQETTC90eJxDRFiJaTkRLiGhhcltnInqbiNYn/9cI+9+avL61RHR+Adv9BBHtJaIVwras201EpySvfwMRPUDkf+FZi2u5k4h2JO/LEiK6sI1cS18imk1Eq4loJRF9P7m9Td0bm+toc/eFiMqIaD4RLU1ey13J7f7eE8ZYm/4DEAKwEcAAACUAlgIYVuh2OWj3FgBdDdt+C2B68vV0AL9Jvh6WvK5SAP2T1xsqULvPAjAWwIp82g1gPoDxAAjA6wCmFsm13Angxyb7Fvu19AQwNvm6CsC6ZJvb1L2xuY42d1+S5+2QfB0B8DGA0/2+J0HQ5E8FsIExtokx1grgGQAXF7hNuXIxgH8kX/8DwCXC9mcYYy2Msc0ANkC9bt9hjM0BcMCwOat2E1FPANWMsXlM7cFPCd/xDYtrsaLYr2UXY+yT5OujAFYD6I02dm9srsOKorwOAGAqx5JvI8k/Bp/vSRCEfG8A24T322HfKYoFBuAtIlpERDckt3VnjO0C1M4OoFtye7FfY7bt7p18bdxeLHyHiJYlzTl8Kt1mroWI6gCMgao5ttl7Y7gOoA3eFyIKEdESAHsBvM0Y8/2eBEHIm9mm2oJf6ETG2FgAUwHcRERn2ezbVq/Rqt3FfD0PARgIYDSAXQDuS25vE9dCRB0AvADgB4yxI3a7mmwrmusxuY42eV8YY3HG2GgAfaBq5SNsdvfkWoIg5LcD6Cu87wNgZ4Ha4hjG2M7k/70AXoJqftmTnJoh+X9vcvdiv8Zs2709+dq4veAwxvYkH8wEgEeRMosV/bUQUQSqYPw3Y+zF5OY2d2/MrqMt3xcAYIwdAvAugAvg8z0JgpBfAGAQEfUnohIAVwF4pcBtsoWIKomoir8GMAXACqjtvi6523UAXk6+fgXAVURUSkT9AQyCuhBTLGTV7uQU9SgRnZ70ErhW+E5B4Q9fks9DvS9AkV9L8tyPA1jNGLtf+KhN3Rur62iL94WIaomoU/J1OYBJANbA73vi52qzV38ALoS6Cr8RwG2Fbo+D9g6Auoq+FMBK3mYAXQC8A2B98n9n4Tu3Ja9vLQrgvSG042mo0+UoVA3ja7m0G0A91Ad1I4AHkYy+LoJr+SeA5QCWJR+6nm3kWs6AOoVfBmBJ8u/CtnZvbK6jzd0XACMBLE62eQWAO5Lbfb0nMq2BRCKRBJggmGskEolEYoEU8hKJRBJgpJCXSCSSACOFvEQikQQYKeQlEokkwEghL5FIJAFGCnmJRCIJMP8PhhGibfhgcI8AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"####Generating correlated curves:\n",
"np.random.seed(2) \n",
"def gen_correlated_curve(rho, num=100):\n",
" num_with_runup = num+5000\n",
" y = np.zeros((num_with_runup,))\n",
" for i in range(1,num_with_runup):\n",
" y[i] = rho * y[i-1] + np.random.normal()\n",
" return y[-num:]\n",
"\n",
"timeseries = gen_correlated_curve(0.99, 3000)\n",
"print(timeseries.mean())\n",
"plt.plot(timeseries)\n",
"plt.axhline(0)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x129857410>"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxV1b338c8vJxMhCSETIRMJIcwzIQJSUKFlqEqrrYqdJ2odqn06XGvvffq09/b23j5Prdpap2qrrYpaJ2pxQAQcGEyYA0kgBELmhMzzcM56/sjBG0MwJ+Qk++Sc3/v1ysucvdc5+5ctfNlZe621xRiDUkop7+VndQFKKaWGlwa9Ukp5OQ16pZTychr0Sinl5TTolVLKy/lbXUB/oqOjTUpKitVlKKXUqLF///5zxpiY/vZ5ZNCnpKSQnZ1tdRlKKTVqiEjRxfZp141SSnk5DXqllPJyLgW9iKwVkXwRKRCRu/vZLyLygHP/ERFZ2GtfhIj8XUTyRCRXRJa68wdQSin1yQYMehGxAQ8C64CZwEYRmdmn2Tog3fm1CXio1777gTeMMdOBeUCuG+pWSinlIleu6DOBAmNMoTGmE9gMbOjTZgPwlOmxF4gQkYkiEg6sAB4HMMZ0GmPq3Vi/UkqpAbgS9AlAca/XJc5trrSZDFQDfxaRgyLyJxEZO4R6lVJKDZIrQS/9bOu75OXF2vgDC4GHjDELgBbggj5+ABHZJCLZIpJdXV3tQllKKaVc4UrQlwBJvV4nAmUutikBSowx+5zb/05P8F/AGPOoMSbDGJMRE9PvmH81Cuky2EpZz5UJU1lAuoikAqXATcDNfdpsAW4Xkc3AZUCDMaYcQESKRWSaMSYfWAUcd1v1yiPtOVXDC/uLOXS2nuK6Vq6YFsuNGUlcMS0Gf5uO6FVqpA0Y9MaYbhG5HXgTsAFPGGOOicgtzv0PA1uB9UAB0Ap8o9dH3AE8LSKBQGGffcrLPJd1lnteziE82J9FkyJZNiWKN3Iq2Xa8kstSI3nsaxmEBwdYXaZSPkU88VfrjIwMo0sgjD4PbD/JvdtOsGJqDH/80kJCg3quI7rsDl46UMK/vpLDlNgwnvzmYmLDgi2uVinvIiL7jTEZ/e3T36OVW7x1rIJ7t53guoUJPP61jI9CHiDA5seNi5N5/GuLKapp4cZH9tLQ1mVhtUr5Fg16NWTNHd38fMsxpseF8d/XzyXgIv3wK6bG8OQ3MymubeWnLx3RG7VKjRANejVkv9t2gorGdn71+TkXDfnzFqdE8uM109h6tIKn950doQqV8m0a9GpIckob+PMHp7k5M5lFk8a79J7vfGoyK6bG8MvXjpNf0TTMFSqlNOjVkPz+nZOEjwngJ2unu/wePz/h3hvmERrkz7+9kqNdOEoNMw16dclK6lrZdrySjZnJjBszuCGT0aFB/HjNND48U8uWw33n3yml3EmDXl2yv+4tQkT48pJJl/T+GzKSmJ0Qzq+35tHS0e3m6pRS52nQq0vS1mnnuaxiPjNzAgkRYy7pM2x+wi+unU1FYzsP7ihwc4VKqfM06NUlefVQKfWtXXx9WcqQPmfRpPF8fkECj79/mqrGdvcUp5T6GA16dUme2lPE9LgwMlMjh/xZP1g9FbvD6FW9UsNEg14N2qnqZo6XN3JDRhIi/a1QPTjJUSHcsDiJZz48S0ldqxsqVEr1pkGvBu2NnAoA1s2Jc9tn3nHVFESE32/Xq3ql3E2DXg3a1qPlLEyOYOK4S7sJ25+J48bwpcuS+fuBEs6ca3Hb5yqlNOjVIBXVtHCsrJH1cya6/bO/d0Ua/n7CQztPuf2zlfJlGvRqUF53dtusne2+bpvzYsOC2ZiZzIsHSiitb3P75yvlqzTo1aC8nlPB3MRxJI4PGZbP37RiMiLwyC69qlfKXTTolctK69s4XFzPutnu77Y5Lz5iDF9YlMjmrGIdV6+Um2jQK5e9k1cFwJpZE4b1ON9bOQW7w/Dou4XDehylfIUGvXLZ+yerSYgYQ2r02GE9TnJUCBvmxfP0vrPUNHcM67GU8gUa9Moldodh96kalk+JdsskqYHceuUU2rvtPPHB6WE/llLeToNeueRIST1N7d0sT48ekeNNiQ1l/ZyJPLm7iIZWfb6sUkOhQa9c8v7Jc4jA5VNGJugBbr9yCs0d3fxl95kRO6ZS3kiDXrnkvYJzzIoPJ3Js4Igdc8bEcFbPmMATH5ymWderV+qSadCrAbV0dHPwbN2IXs2fd8dVU2ho6+Jve4tG/NhKeQsNejWgD0/X0mU3fGpKzIgfe15SBCumxvCn9wpp67SP+PGV8gYa9GpA7508R5C/Hxkp4y05/h1XTeFccyebs85acnylRjuXgl5E1opIvogUiMjd/ewXEXnAuf+IiCzste+MiBwVkUMiku3O4tXI2He6hoXJ4wkOsFly/MUpkVyWGskjuwrp6NareqUGa8CgFxEb8CCwDpgJbBSRmX2arQPSnV+bgIf67L/SGDPfGJMx9JLVSGru6Ca3vJHFbniS1FDccVU6FY3tPJ9VbGkdSo1GrlzRZwIFxphCY0wnsBnY0KfNBuAp02MvECEiw7cgihoxB8/W4TCw2KJum/MunxJFZkokD7xTQGunjsBRajBcCfoEoPdlVIlzm6ttDPCWiOwXkU2XWqiyRtaZOvwEFiRbG/Qiwo/XTqO6qUPH1Ss1SK4EfX/z3c0g2lxujFlIT/fObSKyot+DiGwSkWwRya6urnahLDUSss/UMmNiOKFB/laXwuKUSK6cFsPDO0/R0KazZZVylStBXwIk9XqdCJS52sYYc/6/VcDL9HQFXcAY86gxJsMYkxETM/LD+NSFuuwODhXXszjF2v753n60ZhqN7d26Xr1Sg+BK0GcB6SKSKiKBwE3Alj5ttgBfdY6+WQI0GGPKRWSsiIQBiMhY4DNAjhvrV8Mot7yR1k47iyZZ223T26z4cWyYH8/j75+mpK7V6nKUGhUGDHpjTDdwO/AmkAs8b4w5JiK3iMgtzmZbgUKgAHgMuNW5fQLwvogcBj4E/mmMecPNP4MaJlln6gAsGz9/MT9ZOx0R+PXreVaXotSo4FLHqzFmKz1h3nvbw72+N8Bt/byvEJg3xBqVRfYX1ZI4fgwTx42xupSPSYgYw3dXpHH/9pN8bWktmRYP/VTK0+nMWNUvYwxZZ+rI8KBum95uWZnGxHHB/OIfx7A7+o4NUEr1pkGv+lVS10Z1UweLPOhGbG9jAm3cs34Gx8oadcEzpQagQa/6dai4HoAFSREWV3JxV8+dyIqpMfzmjTzK6tusLkcpj6VBr/p1uLieQH8/psWFWV3KRYkIv/rcbBwG/u2VHHpuFSml+tKgV/06XFLP7PhwAmye/UckKTKEH35mKtvzqnjtSLnV5SjlkTz7b7GyRLfdwdHSBuYneeaN2L6+viyFeUkR/NurOVQ2tltdjlIeR4NeXSC/son2LgfzksZZXYpL/G1+/O6GeXR0OfjRC4dx6CgcpT5Gg15d4HBxAwDzPfhGbF+TY0L516tn8N7Jczy554zV5SjlUTTo1QUOF9cTERJAcmSI1aUMys2ZyayaHsuvX88jp7TB6nJcUtnYzt/3l3D3i0d4PruYbrvD6pKUF9KgVxc4XFLPvMQIRPpblNRziQi/+cJcIkMCue2ZAx69wqUxhnvfyuey/9zOj144zKuHyvjJ34/wmfve5Z28SqvLU15Gg159TEtHNycqm0ZVt01vUaFBPPilBZTWtfHjFw575JBLh8Pwi38c54F3CrhuYQJbv/8pjv1iDY98ZRF+Imx6aj8Hz9ZZXabyIhr06mNyShtwmNHVP9/XokmR3L1uOm8dr+Tx909bXc4Ffr7lGH/ZfYZvL0/lt1+cx8z4cPz8hDWz4njxe8uYEB7MHc8epLHdc38jUaOLBr36mMMlPTNi5yaOjhE3F/Ot5amsmTWB/3o9j/1FtVaX85E3csr5694ivr08lZ99dsYF3WPjxgTwwMYFlDe0c89LRz3yNxI1+mjQq485WtpIQsQYokKDrC5lSESE//vFeSSMH8NtTx+kprnD6pKoamrnpy8dZXZCOP+ybvpF74EsmjSe//Xpqbx2pJy3c6tGuErljTTo1cfklDYwOyHc6jLcIjw4gD9+aSG1rZ3c9dwhS1e5NMbw0xeP0tJp53c3zB9wxvF3V0wmcfwYHtYnaSk30KBXH2lq7+L0uRbmJIzubpveZsWP45fXzuK9k+e47+0TltWx9WgF2/Oq+Je100mfMPD6Qf42P769PJX9RXUe1fWkRicNevWRY2WNAMzyoqAHuCkzmRsyEvn9OwVszx35oYvtXXb+641cpseF8fVlKS6/74bFSUSEBPDIrsLhK075BA169ZHzk4xmx3tX0AP8csNsZieEc9dzhyiqaRnRY/9l9xmKa9v418/OxObn+tyEkEB/vrJkEttyKymsbh7GCpW306BXH8kpbSAuPJiYsNF9I7Y/wQE2HvpSzzj1W/52gLZO+4gc91xzBw++U8Cq6bEsT48e9Pu/ujSFAJsfT3zgecNE1eihQa8+crS0gdle1m3TW1JkCPfdNJ+8ikZ+9srIDF287+0TtHXZ+en6GZf0/piwINbPjuO1I+V0duvyCOrSaNArAJo7uik81+I1I24u5sppsdy5Kp2XDpTy9L6zw3qsE5VNPLPvLF9eMokpsaGX/DnXzIunvrWLDwrOubE65Us06BUAueWNGINXjbi5mO9flc4V02L4xT+ODetSA7/6Zy6hQf7cuSp9SJ/zqfQYxo0JYMvhMjdVpnyNBr0C4GhJz41YXwh6Pz/hvhvnMyE8mFufPjAsk6l25lex60Q131+VzvixgUP6rEB/P9bNjuOtYxW0d43MvQXlXTToFQA5ZQ3EhAURGx5sdSkjIiIkkIe/vIialk6+v/mgWydTddsd/OfWXCZFhfCVpZPc8pnXzIunpdPOjjydKasGT4NeAT0jbnzhar632Qnj+I8Ns/mgoIbfvpXvts99ck8RJyqbuWf9DIL8bW75zCWTo4gODeIfR7T7Rg2eBr2itbObgqpmZsd7943Y/tywOImNmUn8cecp3jpWMeTPq2ps53fbTrByagyfmTnBDRX2sPkJn50Tx/bcKlo6ut32uco3uBT0IrJWRPJFpEBE7u5nv4jIA879R0RkYZ/9NhE5KCKvuatw5T655U04DF49tPKT/PyaWcxJGMcPnz/M6XNDm0z169fz6Ox28H+uneX2B7esmR1HR7eD3adq3Pq5yvsNGPQiYgMeBNYBM4GNIjKzT7N1QLrzaxPwUJ/9dwK5Q65WDYvzM2LnjPKliS9VcICNh768EJtN+N7f9tPaeWlXzLsLzvHywVI2rZhMavRYN1cJGZMiCQ3yZ0e+9tOrwXHlij4TKDDGFBpjOoHNwIY+bTYAT5kee4EIEZkIICKJwGeBP7mxbuVGOaUNRI0NJM5HbsT2J3F8CPfftID8yiZ+9nLOoCdTnWvu4K7nDpEaPZbbrpwyLDUG+vtx+ZQoduZV6Tr1alBcCfoEoLjX6xLnNlfb3Af8BNBpfR7q/IzY0faMWHdbOTWGH6yeyssHS/nzB2dcfp/DYfjBc4eob+viwZsXMibQPTdg+3PltFjKGto5Ualr3yjXuRL0/f3t73s50W8bEbkaqDLG7B/wICKbRCRbRLKrq6tdKEu5Q3uXnZNVzV4/I9ZVt185hU/PnMAvXzvOn95zbdXIB945yXsnz/Hza2Yyc5hvaF8xLRZAu2/UoLgS9CVAUq/XiUDfMV4Xa3M5cK2InKGny+cqEflbfwcxxjxqjMkwxmTExMS4WL4aqryKJuwO43NDKy/Gz0948OaFrJ8Tx3/8M5d7t53AcZEx9g6H4ddbc7nv7ZN8fkECN2cmD3t9ceOCmTExXMfTq0FxJeizgHQRSRWRQOAmYEufNluArzpH3ywBGowx5caYnxpjEo0xKc73vWOM+bI7fwA1NEedN2JneeHSxJcq0N+PB25awBcXJfLA9pNc99BujjifpXteSV0rdz53iEfeLeTLS5L5f1+cN2JdX1dOiyG7qE4fHq5c5j9QA2NMt4jcDrwJ2IAnjDHHROQW5/6Hga3AeqAAaAW+MXwlK3fKKWkgIiSAxPFjrC7Fo/jb/PjNF+ayZHIUv349jw0PfsD0uHBSokKob+1iT2ENIvDjNdO49Yq0Eb2/ceX0WP648xTvnTjHZ+dOHLHjqtFrwKAHMMZspSfMe297uNf3BrhtgM/YCewcdIVqWOWU9cyI9fUbsf0REa5flMinZ03gifdPc6SkgfyKJkTgB6unct3CBJIiQ0a8rgVJEYQF+/N+QbUGvXKJS0GvvFNHt50TlU18a/lkq0vxaOHBAdy1eqrVZXzE3+bHZamROnFKuUyXQPBhJyqa6bLrjdjRaGlaNEU1rZTWt1ldihoFNOh92JHSnhuMOrRy9FmWFgXAHr2qVy7QoPdhOaUNjBsTQLIF/cxqaKZNCCNybCC7T+lTp9TANOh92JESvRE7Wvn5CUsnR7HnVI0uh6AGpEHvo9q7em7E+uqKld5gaVoU5Q3tnKlptboU5eE06H1UfkUTXXbDXB9dsdIbnO+n1+4bNRANeh91fkasjrgZvVKjxxIXHqzDLNWANOh91FGdETvqiQhL06LYq/30agAa9D7qaKneiPUGS9OiqGnp1GWL1SfSoPdB52/EarfN6Kf99MoVGvQ+KK+iiW6H3oj1BonjQ0iODNF+evWJNOh90NGS8zNiNei9wbK0KPYW1mC/yLr5SmnQ+6CjpQ1Ejg0kIUJvxHqDpWlRNLV3c6yswepSlIfSoPdBR0r0GbHeZOlH/fTafaP6p0HvY84/I3audtt4jdiwYNJjQ3WBM3VRGvQ+5nh5I3aH0f55L7MsLYqsM7V0djusLkV5IA16H5PjnBGrI268y9K0KFo77Rzu82xbpUCD3uccKWkgamwgE8cFW12KcqMlk6MQgd0F2n2jLqRB72NyShuYk6g3Yr1NREggs+LDdeKU6pcGvQ9p69QZsd5sWVo0B8/W09Zpt7oU5WE06H3I8fJGHEZXrPRWy9Ki6LQ7yC6qtboU5WE06H3I+RmxcxMjLK5EDYfFKZH4+wkfaD+96kOD3occLW0kOjSICeFBVpeihsHYIH8WJEewR/vpVR8a9D7kaGk9c/VGrFdbmhbN0dIGGtq6rC5FeRANeh/R2tlNQVWzTpTycpenReEwsK9Qu2/U/9Cg9xFHShpwGFiQpP3z3mx+cgTBAX667o36GJeCXkTWiki+iBSIyN397BcRecC5/4iILHRuDxaRD0XksIgcE5FfuPsHUK45eLbnRux8DXqvFuRvY3FKpI6nVx8zYNCLiA14EFgHzAQ2isjMPs3WAenOr03AQ87tHcBVxph5wHxgrYgscVPtahAOnq0jNXos48cGWl2KGmbL0qI5UdlMdVOH1aUoD+HKFX0mUGCMKTTGdAKbgQ192mwAnjI99gIRIjLR+fr8wywDnF/6dIQRZozhYHG9dtv4CH28oOrLlaBPAIp7vS5xbnOpjYjYROQQUAVsM8bs6+8gIrJJRLJFJLu6utrV+pULSuraqG7qYEGyBr0vmJ0wjrBgf122WH3ElaDvbyxe36vyi7YxxtiNMfOBRCBTRGb3dxBjzKPGmAxjTEZMTIwLZSlXHSzu6Z9fkDze4krUSLD5CUsmR/GBXtErJ1eCvgRI6vU6ESgbbBtjTD2wE1g76CrVkBw8W0dwgB/T48KsLkWNkMvToiiubaO4ttXqUpQHcCXos4B0EUkVkUDgJmBLnzZbgK86R98sARqMMeUiEiMiEQAiMgZYDeS5sX7lgoNn65mbGIG/TUfT+oplU6IB7adXPQb8m2+M6QZuB94EcoHnjTHHROQWEbnF2WwrUAgUAI8Btzq3TwR2iMgRev7B2GaMec3NP4P6BB3ddo6XNWr/vI9Jjw0lOjRIx9MrAPxdaWSM2UpPmPfe9nCv7w1wWz/vOwIsGGKNagiOlTXSaXewIEn7532JiLAsLYrdp2owxuiyFz5Of5f3cgeK6gD0it4HLZ8STXVTB/mVTVaXoiymQe/lss7UkhQ5hgnh+uhAX7Nias/otV35OlzZ12nQezFjDNln6shMibK6FGWBuHHBTI8LY6cGvc/ToPdip6pbqGnpJDNV++d91cqpMWQX1dLc0W11KcpCGvRe7MPTPY+UW5wSaXElyiorp8XQZTc6S9bHadB7sawztUSHBpIaPdbqUpRFMiZFEhJoY2d+ldWlKAtp0HuxD0/XsjglUofW+bBAfz+WpUWz60Q1PaOglS/SoPdSZfVtlNa3abeNYuW0GErq2ig812J1KcoiGvReKutMT/98ZqoGva+7wjnMckeedt/4Kg16L7XvdC2hQf7MmBhudSnKYkmRIaTHhvKOBr3P0qD3Ulmna1k0aTw2P+2fV7B65gQ+PF1LQ1uX1aUoC2jQe6GqxnZOVjWzNE0nSqkeq2fE0u0w7Dqhk6d8kQa9Fzr/wInlzqVqlZqfNJ6osYG8fbzS6lKUBTTovdD7J2uICAlgpvbPKyebn3Dl9Fh25lfRZXdYXY4aYRr0XsYYwwcF57g8LRo/7Z9XvayeMYHG9u6PRmQp36FB72VOVbdQ0djO5dpto/r4VHo0gf5+bM/V0Te+RoPey3xQoP3zqn9jg/xZlhbFtuOVOkvWx2jQe5kPCs6RFDmG5KgQq0tRHmjNrDjO1raSW64PI/ElGvRepNvuYE9hjV7Nq4v69MwJ+Am8kVNudSlqBGnQe5EjpQ00tXezLE2DXvUvOjSIxSmRvJ5TYXUpagRp0HuRHXlV+EnPTTelLmbd7DhOVjVTUNVsdSlqhGjQe5HtuVVkTIokIiTQ6lKUB1s7eyIAbx7Tq3pfoUHvJcob2jhe3shVM2KtLkV5uLhxwSxIjuB17af3GRr0XuL82OjVGvTKBWtnxZFT2khxbavVpagRoEHvJd7JqyI5MoS0mFCrS1GjwPo5Pd03rx3Rq3pfoEHvBdo67XxQcI6rpsfqYwOVS5IiQ1iYHMGrh0qtLkWNAJeCXkTWiki+iBSIyN397BcRecC5/4iILHRuTxKRHSKSKyLHROROd/8ACnafOkdHt4NV2m2jBmHD/ATyKprIr9DJU95uwKAXERvwILAOmAlsFJGZfZqtA9KdX5uAh5zbu4EfGmNmAEuA2/p5rxqit3OrGBto08cGqkFZP2ciNj9hy2G9qvd2rlzRZwIFxphCY0wnsBnY0KfNBuAp02MvECEiE40x5caYAwDGmCYgF0hwY/0+r9vu4K1jFVw5PZYgf5vV5ahRJCYsiMunRPPqoTJd+8bLuRL0CUBxr9clXBjWA7YRkRRgAbCvv4OIyCYRyRaR7OpqfQqOq/adrqWmpZOr5060uhQ1Cm2YF09JXRsHztZbXYoaRq4EfX939/r+8/+JbUQkFHgRuMsY09jfQYwxjxpjMowxGTExMS6UpaBn1ERIoI0rpmn/vBq8z8yaQJC/H1v0pqxXcyXoS4CkXq8TgTJX24hIAD0h/7Qx5qVLL1X11W138EZOOatnTCA4QLtt1OCFBQewesYEXjtSTrc+ecpruRL0WUC6iKSKSCBwE7ClT5stwFedo2+WAA3GmHLpGev3OJBrjLnXrZUr9hTWUNfaxWe120YNwbXz46lp6eSDUzVWl6KGyYBBb4zpBm4H3qTnZurzxphjInKLiNzibLYVKAQKgMeAW53bLwe+AlwlIoecX+vd/UP4qn8eKWdsoI2VU7WrS126K6bFEBbsr2PqvZi/K42MMVvpCfPe2x7u9b0Bbuvnfe/Tf/+9GqIuu4M3jlXw6ZnabaOGJsjfxrrZcfzzSDntn7frnycvpDNjR6l38qqob+3i2vnxVpeivMCG+Qm0dNp5J0+fJ+uNNOhHqb/vLyEmLIgV6dpto4ZuyeQoYsKCtPvGS2nQj0LnmjvYkVfFdQsS8Lfp/0I1dDY/4Zq58ezIq6a+tdPqcpSbaUqMQq8eKqPbYbh+UaLVpSgvct3CBDrtDv5xuO/oaTXaadCPMsYYXsguZl7iOKZOCLO6HOVFZieMY8bEcF7YX2J1KcrNNOhHmWNljeRVNPGFjKSBGys1SF9YlMiRkgZd0dLLaNCPMpuzzhLo78c1OklKDYPPzY/H3094Ibt44MZq1NCgH0Wa2rt4+UAp18yN1weAq2ERFRrEqhmxvHKolC5dEsFraNCPIq8cLKWl085Xlk6yuhTlxb64KIlzzZ3s0DH1XkODfpQwxvDUniLmJo5jflKE1eUoL7ZyWgyxYUFsztLuG2+hQT9K7Dtdy8mqZr68RK/m1fAKsPlxQ0YSO/OrKK1vs7oc5QYa9KPEX/cUMW5MANfM1SUP1PC7cXESBnhOr+q9ggb9KFBc28obxyq4cXESYwJ1wSk1/JIiQ1iRHsNzWWd1nXovoEE/Cjz+/mkE+PqyFKtLUT7k5suSqWzs0IXOvIAGvYerb+3kuaxirp0XT3zEGKvLUT5k1fRYYsOCeHrfWatLUUOkQe/h/ra3iLYuO99ZMdnqUpSP8bf5sTEzmV0nqjl9rsXqctQQaNB7sPYuO3/ZXcSKqTHMmBhudTnKB31pSTIBNuHJ3WesLkUNgQa9B3v5YCnnmjv4rl7NK4vEhgVz9dx4XsguprG9y+py1CXSoPdQDofhsfcKmRUfzrK0KKvLUT7sm5en0tJp54VsXdVytNKg91Bv51ZSWN3Cd1emIaKP3VXWmZM4joxJ43ly9xnsDmN1OeoSaNB7qEffLSRx/BjWz46zuhSl+MblqZytbWXb8QqrS1GXQIPeA+0vqiW7qI5vLU/VRwUqj7B2dhyp0WP5w44CjNGr+tFGU8QDPbyrkHFjArhBHy6iPITNT/jeyjRyShvZeaLa6nLUIGnQe5jc8ka2Ha/kG5enMDbI3+pylPrI5xYkkBAxhj+8o1f1o40GvYf5w44CQoP8dbkD5XEC/f347srJ7C+qY29hrdXlqEHQoPcgBVXNbD1azleWTtInSCmPdENGEjFhQdy7LV+v6kcRl4JeRNaKSL6IFIjI3f3sFxF5wLn/iIgs7LXvCRGpEpEcdxbujf64o4AgfwXHb5UAAAzzSURBVD++vTzV6lKU6ldwgI07V6WTdaaOt45XWl2OctGAQS8iNuBBYB0wE9goIjP7NFsHpDu/NgEP9dr3F2CtO4r1ZkU1Lbx6uIwvXTaJqNAgq8tR6qJuWpxEWsxY/vv1PH2u7CjhyhV9JlBgjCk0xnQCm4ENfdpsAJ4yPfYCESIyEcAY8y6gHXoDeGjnKWx+wiZd7kB5OH+bH3evm0HhuRY2f6grW44GrgR9AtD7MTMlzm2DbfOJRGSTiGSLSHZ1tW8N3yqtb+PFAyXcmJHEhPBgq8tRakCrZ8SSmRrJfW+fpL610+py1ABcCfr+5t/3vQvjSptPZIx51BiTYYzJiImJGcxbR71Hd53CGPjuSr2aV6ODiPDza2ZS39bFv7+Wa3U5agCuBH0J0HvmTiJQdgltVD+qGtt5NquY6xcmkjg+xOpylHLZrPhx3LJyMi8eKGGXTqLyaK4EfRaQLiKpIhII3ARs6dNmC/BV5+ibJUCDMabczbV6pYd3FdJtd/C9K9KsLkWpQbvjqnQmx4zlnpeO0tzRbXU56iIGDHpjTDdwO/AmkAs8b4w5JiK3iMgtzmZbgUKgAHgMuPX8+0XkWWAPME1ESkTkW27+GUatioZ2/raviOsXJpISPdbqcpQatOAAG7+5fi5lDW382ys5OrbeQ7k0x94Ys5WeMO+97eFe3xvgtou8d+NQCvRmf9hxEofD8P1V6VaXotQly0iJ5Aerp3LvthPMT4rgazqr2+PozFiLlNS18lxWMTcsTiIpUvvm1eh2+5VTWD0jln9/7TjZZ3Q0tafRoLfI77cXICLccdUUq0tRasj8/IR7b5xPUmQIm/66nxOVTVaXpHrRoLfAmXMt/P1ACTdnJjNx3Biry1HKLcKDA/jz1xfj7yfc/Ng+CqqarS5JOWnQW+D+7ScJsAm3XqkjbZR3SYkeyzPfWQLAzY/t5XhZo8UVKdCgH3EFVU28cqiUry5NITZMZ8Eq7zMlNpRnvnMZfiJc/9BuXj+qI62tpkE/wn739klCAmx8V9e0UV5s6oQwttx+OdMnhvG9pw/w6625tHfZrS7LZ2nQj6BjZQ3880g537g8VVeoVF4vNjyYZ7+zhJsvS+aRdwtZf/97fHhaR+RYQYN+hBhj+NU/cxkfEsB39Gpe+YjgABv/+fk5/O1bl9Fpd3DDI3u49en9FFbrjdqRpEE/QrbnVrH7VA13rZ7KuDEBVpej1Ihanh7Nm3et4M5V6ezMr+bTv3uXuzYfJKe0werSfIJ44pTljIwMk52dbXUZbtNld7DmvncBePOuFQTY9N9X5buqmzp4aOcpnss6S0unnctSI/nW8lRWzZiAza+/hXCVK0RkvzEmo799mjgj4Om9RRRWt3DPuhka8srnxYQF8b+vmcmee1bxs/UzKKlrY9Nf97Pqtzt57N1Calt0fXt30yv6YVbV2M6q3+5iXlIEf/1WJiJ6xaJUb912B28cq+DPH5xhf1EdgTY/1s2J4+bMZDJTI/XvjIs+6YrepUXN1KX75WvH6bA7+PfPzdY/sEr1w9/mx9Vz47l6bjz5FU08s6+Ilw6W8uqhMqbEhnJzZjLXL0xkXIje27pUekU/jHbmV/H1P2fxg9VTuXO1rlCplKtaO7t57Ug5T+87y+HieoL8/fjs3Il86bJkFiaP14umfnzSFb0G/TBp7uhm3f3vEmDz4/U7P0WQv83qkpQalY6VNfDMvrO8crCUlk476bGhXLcwkc8tiNe1onrRoLfAj144zIsHSnhu01IyUyOtLkepUa+lo5sth8v4+/4S9hfVIQJLJ0fx+QUJrJkdR3iwb3ftaNCPsH8cLuOOZw9y+5VT+NGaaVaXo5TXOXOuhZcPlvLKoVKKalrx9xOWTI5i9YxYVs2Y4JPPeNCgH0Elda2su/890mJCeeGWpTqcUqlhZIzhwNl63jpewdvHKzlV3QLA9LgwrpgWy8qpMSyaNJ5Af+//e6hBP0Ka2rv44sN7KK1r47XvL2dSlD4HVqmRVFjdzPbcKt7OrWR/UR3dDsPYQBvLpkSzcmoMK6fGeO3Vvg6vHAHddge3P3OQk1XN/OUbizXklbLA5JhQJseE8p0Vk2lq72L3qRp2nahmV341245X9rSJHsunZ07gmnnxzIoP94kRPHpF7wYOh+Gel4+yOauYX183h42ZyVaXpJTqxRhD4bkWduVXsyO/ij2nauh2GCbHjOWaufFcOz+etJhQq8scEu26GUZddgc/euEwrx4q05uvSo0StS2dvJFTwZbDpew7XYsxsCA5go2ZyVw9dyIhgaOvs0ODfpg0d3Rz57MH2Z5XxY/XTOPWK9J84tdApbxJZWM7Ww6VsTnrLKeqWwgL8ufa+fFszExmdsI4q8tzmQb9MDhUXM+dmw9SXNvKLzfM5stLJlldklJqCIwxZBfV8ey+s/zzaDkd3Q7mJIzjpswkrp0XT5iHj9PXoHejhrYuHtl1ikffLWRCeDC/u3G+TohSyss0tHbx8sESNmcVk1fRREigjWvmxnNTZhLzkyI88jd3DXo3qGpq54XsEh59t5CGti6uW5DAz6+ZpQstKeXFjDEcKq7n2Q/P8o/D5bR12ZkeF8ZNi5NYN2ciE8KDrS7xI0MOehFZC9wP2IA/GWP+q89+ce5fD7QCXzfGHHDlvf3xhKA3xlBU08ruUzVsz61k54lq7A7DFdNi+PGaacyKHz19d0qpoWtq72LL4TKe/fAsOaWNAMyKD2d5ejQZkyJZkBxBtIXPgh5S0IuIDTgBfBooAbKAjcaY473arAfuoCfoLwPuN8Zc5sp7++OuoDfGYHcY7MZgDB99b7cbmju6aensprm9m+aObupbu6hobKesvo2Tlc2cqGyixvkAhInjgvncggSuX5jIlNjRPQRLKTV0eRWNvJNXxY68Kg4V19Nl78nR6NBApk4II3H8GGLDgokNDyI2LIiYsCDGBvkzJsDGmAAbwYE2gvz9sIlg8xO3dAUNdcJUJlBgjCl0fthmYAPQO6w3AE+Znn819opIhIhMBFJceK/bzP/lW7R0dONwhvqlCAvyJy02lFUzYpmTGMGytCgmR4/1yD45pZQ1pseFMz0unFuvmEJ7l52jpQ0cLq4nv6KJE1XN7Myv5lxzB67GkAjYRIgNC2L3T1e5vV5Xgj4BKO71uoSeq/aB2iS4+F4ARGQTsMn5sllE8l2ozV2igXPnX+SM4IE93MfOiwL0nFyMnpcLDfqcnALknks+3kWH/rkS9P1dyvb9d+pibVx5b89GYx4FHnWhHrcTkeyL/crjy/S8XEjPSf/0vFzIk86JK0FfAiT1ep0IlLnYJtCF9yqllBpGrqzdmQWki0iqiAQCNwFb+rTZAnxVeiwBGowx5S6+Vyml1DAa8IreGNMtIrcDb9IzRPIJY8wxEbnFuf9hYCs9I24K6Ble+Y1Peu+w/CRDY0mX0Sig5+VCek76p+flQh5zTjxywpRSSin38f7HriillI/ToFdKKS/n00EvIv9XRPJE5IiIvCwiEb32/VRECkQkX0TWWFnnSBKRL4rIMRFxiEhGn30+eU7OE5G1zp+9QETutroeq4jIEyJSJSI5vbZFisg2ETnp/O94K2scaSKSJCI7RCTX+ffnTud2jzgvPh30wDZgtjFmLj1LNfwUQERm0jNCaBawFvijczkHX5ADXAe823ujj5+T80uBPAisA2YCG53nxBf9hZ4/A73dDWw3xqQD252vfUk38ENjzAxgCXCb88+HR5wXnw56Y8xbxphu58u99Izzh55lGjYbYzqMMafpGU2UaUWNI80Yk2uM6W9Wss+eE6ePlgIxxnQC55fz8DnGmHeB2j6bNwBPOr9/EvjciBZlMWNM+fmFHI0xTUAuPSsDeMR58emg7+ObwOvO7y+2pIMv8/Vz4us//0AmOOfO4PxvrMX1WEZEUoAFwD485LyMvgcjDpKIvA3E9bPrZ8aYV51tfkbPr15Pn39bP+29ZhyqK+ekv7f1s81rzokLfP3nVy4QkVDgReAuY0yjpyyG6PVBb4xZ/Un7ReRrwNXAKvM/kwpcWfZh1BronFyEV58TF/j6zz+QShGZaIwpd65cW2V1QSNNRALoCfmnjTEvOTd7xHnx6a4b50NR/gW41hjT2mvXFuAmEQkSkVQgHfjQiho9iK+fE13O45NtAb7m/P5rwMV+M/RKzocvPQ7kGmPu7bXLI86LT8+MFZECIAiocW7aa4y5xbnvZ/T023fT82vY6/1/incRkc8DvwdigHrgkDFmjXOfT56T85wP2LmP/1nO41cWl2QJEXkWuIKeZXgrgZ8DrwDPA8nAWeCLxpi+N2y9logsB94DjgIO5+Z76Omnt/y8+HTQK6WUL/DprhullPIFGvRKKeXlNOiVUsrLadArpZSX06BXSikvp0GvlFJeToNeKaW83P8Hc97cRi8quEsAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"kdeplot(timeseries)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"now fit an AR(1) model to the timeseries. We do not assume the time series is centred on zero (because in production you won't know the centre). Hence the timeseries is shifted to zero by the 'center' parameter. Fitting the AR(1) model has the consequence of also finding where the timeseries is centred, which gives us an HPD for the timeseries mean. "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [center, tau, rho]\n",
"Sampling 2 chains, 0 divergences: 100%|██████████| 2000/2000 [00:01<00:00, 1144.14draws/s]\n"
]
}
],
"source": [
"with pm.Model() as ar1:\n",
" rho_ = pm.Uniform('rho',-1,1) #we assume process is stationary, so -1<k_<1\n",
" tau_ = pm.Gamma('tau',mu=1,sd=1)\n",
" center = pm.Normal('center', mu=timeseries.mean(), sigma=5) #set the prior for the true mean initially centred on the population mean\n",
" likelihood = pm.AR1('likelihood', k=rho_, tau_e=tau_, observed=timeseries-center)\n",
" trace = pm.sample(target_accept=0.9)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Analysis:\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-13.55096279, -9.67315923])"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#HPD is way off:\n",
"pm.stats.hpd(trace['center'])"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x129984690>]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2dd7wdRdnHf8+5Lb3fNEJIIbQQauhNINIVaYrYRXl9xYrlDYKKhSYqqCBFrKggCggYapTeA0lMQk+BNEgjPbm59555/9idPbOzM7Oz7Zx7T+b7+eSTe3ZnZ2Z3Z5955plnniHGGBwOh8NRn5RqXQGHw+FwFIcT8g6Hw1HHOCHvcDgcdYwT8g6Hw1HHOCHvcDgcdUxjrSsgMmTIEDZmzJhaV8PhcDi6FS+++OIqxlir6lyXEvJjxozBjBkzal0Nh8Ph6FYQ0Vu6c85c43A4HHWME/IOh8NRxxQm5IloHyJ6lohmEdEMIjqwqLIcDofDoaZITf4nAH7AGNsHwPf83w6Hw+GoIkUKeQagn/93fwDLCizL4XA4HAqK9K75GoAHiein8DqTQ1WJiOg8AOcBwOjRowusjsPhcGx/ZBLyRDQdwHDFqYsAHAvg64yxO4jowwB+C2CKnJAxdhOAmwBg8uTJLiSmw+Fw5EgmIc8YiwhtDhH9CcBX/Z9/B3BzlrIcDkd1mLV4LRpLhD136F/rqjhyoEib/DIAR/l/HwPgjQLLcjgcOfGh657CKb96stbVcOREkTb5zwP4BRE1AtgK3+7ucDgcjupRmJBnjD0JYP+i8nc4HA5HPG7Fq8PhcNQxTsg7HA5HHeOEvMPhcNQxTsg7HA5HHeOEvMPhcNQxTsg7HA5HHeOEvMPhcNQxTsjnzIoNW/HA3OW1robD4XAAcEI+dz5+83P4wp9fwtb2zlpXxeFwOJyQz5vFa7YAAMrMBdR0OBy1xwn5guiOMn75ui049mePYvm6LbWuisPhyAkn5HOGqNY1SM+tz72N+Ss34fYXliS+dt2Wdlz98Ovo6CwXUDOHw5EWJ+QLohsq8kGd03RUP3voNfzi32/gvrnv5Fonh8ORDSfkc6YbK/KZTEwdZe/i9Vvac6qNw+HIAyfkC4J1Q6M883X5NB1Vz6YGAKiZV9HUO/6Lo656pCZlOxxdmSI3DdkuoW5slOf9Uppb6NXsCfnN22oj5G97YXFNynU4ujpOky+I7qfHizb55FK+h6/Jb3HrAxyOLoUT8o6ALBYmLuTb2p13jcPRlXBC3hEQ2ORTmGv4HESp+1qrasrazdvwoeuewturN9e6Ko46wwn5guiG866BvYZSTL1y75qGBifl0zBtznLMWrwW1z82v9ZVcdQZTsjnTCDiuqOQ90mjyXf6Qr7RqfKpcM/PURSFCXki2puIniGiOUR0LxH1K6qsatPRWcYV97+KdZv1PuFMI+UfeW0FLrlnXlFVy0SWfokLqYaS0xvS0NHJn58T8o58KfKLvBnAVMbYJAB3AfhWgWVVlQfnvYsbHpuPH017OXrS/0Z15prP/P4F/OHpRYXVLQvlcnq7ekeGax1iJ+keoCNfihTyuwJ43P/7YQBnFFhWVen0JbjJXbA7WmsCF8oUNvnOsvOqyQJvU85c48ibIoX8XAAf9P8+C8COqkREdB4RzSCiGStXriywOvnR4ButueYrwj/RbrniNcNiKK7Jy7f94ltr8POHXstYs/rHafKOosgk5IloOhHNVfw7FcBnAZxPRC8C6AtgmyoPxthNjLHJjLHJra2tWapTNRr8p9apEPKc7ifi9fMINnR2qq894/pn8Mv/vJk63+0FZ5N3FEWmsAaMsSkxSY4DACLaBcDJWcrqSpS4Jm+Qid1QkRc0+fQulLrb7ugso7HBTcrq4OYaJ+QdeVOkd81Q//8SgIsB3FBUWdWmIuQV5hr/XBatuNYkETOrN7bhgbnL0cFt8prebZtlnPkxU6fh8vtfSVCD+oDPaTibvCNvilStPkpErwN4FcAyAL8vsKyqwrUtk7kmTsZ3ZZt9EkX+3D/OwBf+/BJWbVBa49DkL45KEu7gxscW2FegTuhwLqiOgigsCiVj7BcAflFU/rWkVDJp8t7/cSKcsa63i1SajmfJe94y/K0dnqeRnENzQwntnZ1o63DeNyY6A5t8jSviqDtck0pBg8FcY0tX3Oi74kKZnE6Nd01zo9fE2jpcdEoTFZu8+yQd+eJaVAq42dToXRMjwzu7opDPMPHarrG5V4R8vCbflU1YRePCGjiKwgn5FATmGoPcipt47YryjNc5jZwJNHnpvgMhb2GTN01x1Dsdzk/eURBOyKcgmHg1SOo4Id4VzTWBkM3iQimbaxrszTXGiew6p9P5yTsKwgn5FHBBppx45Wli8shLni1ftwVjpk7DvbOXZc4rMNekuLZDsxiquZHv/WqjyW/HQt7F43cUhBPyKeDCSBXWgBNnX85LoL36zgYAwN9fXJJDbuk3DdEthqp0evH3G/dM1m1pxzvrtiavnMOxHeOEfAq4LOpkDOUyCwn0YDFUnJ+8QrFdunYLTr32SazYYC/Ikspjxhienr/K2AllCVCmy9amT5PNNYde/m9ccf+rwe/3XfUIDr7834nr1pW5+J9zMP3ld4Pf2/FgxlEQTsingGulm7d1Ytx37sPV099InIfKnn/9o29i9pJ1+Nfs5cnrZCkd7nxpKc75zXO486WlijwSFxugM9ckyVIeGC1btxU3CDslvWeI399d+fOzb+Nzf5oR/DaZ8To6y/jIjc/g6fmrqlAzR73ghHwa/A9x49YOAMAtzywKTlWiUJqzUJkm3vL399xpcC/rqsS5O765YgPGTJ2Gl95+DwDwznpvlPDGio2RtFmiUC72F0XpzDI2wl5n/lq/tR2n/OqJ5JXqhpjMWis2tOG5hWvwjdtnV7FGju6OE/Ip4LKIx7BRabFxNmiVkF/ra6ot/mRlHjz6mhe+mY8OejV7eW/Z1hFJG2zknaKcdv4MpNtK4vuu81Z6/PWVmLt0fYpa1Z4FKzcmegbOXOPIGyfkU8CFIf+/Q6GBxtrkFedV+cRhO3LgcCG/aVvUpTGLJh/koTtuUUGdJt/Q1eI/WDJnyToc87PHcPMTC5XnVc/EyXhH3jghnwL+bXINnk8YvrV6E1Zv8gJ16T5WLq9UmnwweSlc/d6mbXh3vX4iNqn869nshSvaohLyPM9Uunx2dH2c7DveYRnRsta8vcYzYc1c/J7yvOp+bTpDp+07kuCEfArkNUPtvnA+6qpHK2k0XyI38agW/qgWFB1w6XQcdFl+HiW9fU1+s8pck4cmr7lvG7mkM9fIQr5egp11KJZMmwR4Nx3QOGqME/IpkLXwJJoVl1eqa0yCPw7b+PV84w61uYb7yaeXJvJ9mZ5NZ5mFRhQ6c02pToW86n1vz/F7HMXghLzA9Y/Ox7T/WrgvWnyHWnMN9BEsufknyWfO87OVDVyIZAkhsHJDGx5/PeF+vIriLrh9Fnb/3gPBb91iqBLJQr57RbTUvRuxA+d3aPNabDv0re2duOBvs4zmvqJYtbENv31yoeu0ugBOyPswxnDlA6/i/L++FJ+WT7wa2q/uXMUmHz1XCddr/2GolO7X392Are3JBaFtqR+7+Vl88nfPKzVv7cSr4szds8KhGHQdj3yLSTYgqSVxAyJxX1wm/a/ML+FcyYPz3sGdM5fix9Oqv9PWl/86Ez/618t47d0NVS/bEcYJeZ+VG9qs09rZrjUCyzTxypJr8n99/u3Q7xffeg/HXf04dvvuA8r0prx55xLXySxYuQlAZS4inIdcnv3d6DR52VavsmV3R9ReWdWdeP3075/H2Aun5Zehz7otnjuwbpGco3oUtjNUd2NDW3QiUofVkFqTJtgfVmV/T+E1ws1LZcYwZuo0jGvtHVOxyp8LVm7E8wvX4OwDR4dOxd1eYwOho8zQ3skSmIni04hBzERh1ykJinqxAHSFqJt8HYUjnvbOMhqIInNEXR2nyafAStsCcN+c5ZERQmUT8Og1gWaX4tvnAoNr2bH1Ywwf+NWTmHrnHOEY/8N8baO/e1G7YgKUa+5fvnUmTvnVExGBvKmtI2JT58/z079/IVoXRDXeWovGJ95YiblL12XOh49ISiTa5OPvThxBPvb6Sjy/cI0xfS3s4qoSJ/94Oj77hxcUZ9LzwqI1eHlZdRbKTbjofnz/nnlVKStPnJD3SfIdBLLQcM26Le344l9ewmf+8HzouOlj5tp9EhMHx1YrFPPmHjZMMhPFlc83527vLEdMVvy27p29LLRKlR+f+P0Hceb1z4Su4VVftbFNOMaUf3cFPvHb53HKr57MnA9/ZyUiqzYFRZpP/e55fPjGZ5TpsnhJFcGqjW34z6srcs3zrBuewUm/LD7kBf9Gbnn2rcLLyhsn5APsBYmNTZ5vh7f0vS2h4yabfHuG4bvtpSohEgwgApu8OQ/uhinW94gJQ9CvR9T6p8pqjqQFqzoo8VBEk+9aMj81/L7E4b9x4rXGMntreyfGTJ2GO1+KD2vdtbqX9Jzzm2dx1g1PK9vcBX+bhS9ZOGrUGifkfRJp8lbqlvqwKRSxbjNsG2y03Q1b2zF78dpI9eSJzLicmnyhJJtrTJqjKU/lqEa0yUfqF3+vm9o6MGbqNPzuyYWxaYuGMWD6y+9i3x8+FPJ64u9bDNuQtQPbsq0TS/xgcXmzYr030vr5w6/Hpq11P7yto5yLq+3T81fjhUXqFct3zlyKf9m4XNeYTEKeiM4ionlEVCaiydK5C4noTSJ6jYiOz1bN4hEbpWkzEDlt0jSmTcCTCnlRgNl4MZz3pxfxy/+8GTnOZahtuYEm31kOXUMU7QDTeouEbPIpJl75XMgfnl4Un7ggxC7vsvtfwXub27FEGNl1CFv+qcx4Nz42HzMWRe3tptv/7B9ewOFXPmKdPgmVPYDt9fRajT4OveLf2PVitYdZGmrdaWUhqyY/F8DpAB4XDxLRHgDOBjARwAkAfk1E+YVWLADTRJ8ubRqFvjLxmr3Z/PBfLwd/2+Q3e8laqS7e/9xFsWxtrvEu3JbAG8gk7OM1+eTPil9RaxMHpxJILnpfOmeNy+9/FWfeULG382QrN7RhzhL1xO8zC1ZnraoW/hp0z3Tzto5gpFhrVm3clmt+3XlRVyYhzxh7hTH2muLUqQBuY4y1McYWAngTwIFZyioa0QQQDVvAAhu7mNYkQNIshqrknxwrrwyxDBaNoxMI+biJV+5dI2nYhGjdbe4lTshnscnbyPgf/+tlXPC3WfaZpoA/a5WZTIzNc9WDr2GbJmyDeO0HrrWf+M3az7V1dOL7d8/F2s3bjPn93x1zcOp1T2HlhrZChOKfn30L10yPNxUVQfcV8cXZ5HcAsFj4vcQ/FoGIziOiGUQ0Y+XK2vnsim1S1hx/+e83MeGi+7HR96XPsvS8YpPXZ5LmA7GJgyIPs2WffdsRSqPgXSNCRPrt/wz5qdY2ibcT6XQtPrmjf/pobBrOzU8uxJ0zoztl5YlqBBfY5CVV/m8vhBe4peFNYVOYrALq7lnL8Mdn3gq2YtSZa3iZRe3De/E/5+KaFLuw5UFefVZ7Zxljpk7DTx9U6cbFECvkiWg6Ec1V/DvVdJnimPIxMcZuYoxNZoxNbm1tta137oSEvPRGb5/h9Vdck8myCQR/MKYgW+k0eYtjwlshAnyFPNCUbd30mxoqfvKx5hCLm+GCb0CvpsploiavmG94aN47eN1iyXxXcSMMRnDCa+fPXa6jrm0kETRTfv5Ybh4u/F0Enbom4yF9mgEAqze1dZnnnhdp3JpV8Hf7+6eq5xAQK+QZY1MYY3sq/t1tuGwJgB2F36MALNOkrTmPv74y5GsrT7zy9pokwqKOPi2em2Hes/Lq+PRRcwqHsYpXR8UWb6fKNwk2+dDEKzxfaN4Zipiy5OWfMHE4AODAsYNCHZR8H4wB593yIo67OjQVpKSaoqato1OrAAQjOIVZUK6jzvSWl6BRsWFrO/73zy8qN5GvBNXjv9UM6u0J+TWb7Ozhp177JC6/r/pxddLQjU3yhZlr7gFwNhG1ENFYABMAPB9zTc246fEFod+yDZgPTxe/txmzF69NtGBF/ugnjeoPAOjZpJ+HTtOgVNquLCxk7Upvk48iCm5uXojY5Am4f+472OeHDwfHxBQ6r6Wy1Lcwxow2+URIEumddVuD/W7zpK2jE7te/ACufKAyDA+UAzB1iGlhIjN0uAYC5Y4Xl+D+ue/gWoX3VUXJMVesh79tZVtH2Wq0O3vJOtwofXv1Ti0W9mV1oTyNiJYAOATANCJ6EAAYY/MA3A7gZQAPADifMVZYfNiNbR2YsWiNUoPkPDN/NcZ/577QqkodsjDiQu2c3zyHU697qqKBmSZe48owvuzouT89syhxfhFTEal/B66bmusWr9mMfX74MP7gDzEr4Y2TNFgWEtYn/kIYOUkTvmVmXvGaZeL1mJ89itN//bR9BpbwydKwX35lTYRqEVxFMw7PZej6tKTyodPW/oaKp1RzQ1QkcOWAmzF1NnnViLdWu4x1VeLcs4sgq3fNXYyxUYyxFsbYMMbY8cK5Sxlj4xljuzLG7s9eVT1vrtiIM294xqih/fm5t9BZZlZx0GWbvOziZmPViOuxk06u3viYWeOJ81ABwh8nQ2WlJU8nr3zl8IiCss+5mMoTBIbFUJLgfmV5JeSB7LrJWDjwmVyfM65PL6Q3KzZLyQNua9e5lQbeNQw45VdP4K6ZS0LKgsm7i5NUPCQJN807qeZGhZDn9eImec1rFkcunCJNTNVEfiVpg8vVIiZdXax4rfgg69OMHtQLALB4zRZ9Ih/5Bcqai817kuPBVI6H8/jyrTMxZuo0ZZokqCdezZ0Vt8lXth0M13nd5nZ8959zg/oslkI0iBoqY0z58YcmUOPMNcJvse7y84jzz6+FT3Nc+GNRG567dD2+/rfZ4bZh6NTSsiUvIS+NQvTaebINbLoTcmcle5bZUovIo/Uh5DUToyItfuNV7SMqC6dYIW9hrokLec4/mHtnR+ej0zQDZfyXSB3CFeaCpyzZ5DlXT38dtzz7VuBdFC0jQahhREMGV+oZ1eTFopI+j5C5IEcvjy3bOnHXTHXcFrFdqUJGV7RhtcZussmPmToN67e2Jxb+NuEHOG2dUSE/Zuo0fOeuORH3z7hHGhrh1cBcY2sSmb14rXVa8dF/4/bZyv0aymWG9VvbzXULnmH1nkt9CHmuQRjSxHkGiES8UjTmGhNlWWVH+KfctrJqb+ogX4b7YCy64lUKb8Dz1G2wLQvTOL/ZuHxCNnnhfkyP5tV3wmFm7529TBIyHk+9uQr/fuVdfUYW/Hjay/j632Yrz4kd6svL10dWfvJn3RG6L64Zq231IkvWbEmsIW/Y6q3rmDZnOU76hTlSY6DJSzb5vz73dmTuRieguorX5E1PxE/mPrtgNU697inc/KQ67VurN4Xai/jo79AEaLv5yQXY65KHjNstptn9LSt1sWmIzey/jfbNkYWjvFjFavu/mDJMNr505hobmzyk33beNbb1UT3bt1ZvDvLQmTQ6pQ6xLHnXmOy6J//yScy/7KTg95dvnYkrH3g1UqeP3fyczS0YWWHYPUy8tw9e+xQA4IaP7y/UwzeNdarfc3jiNX8B8PJyc8x1LuRblOYaswvli2+twZwl6yrHa2yveSXmXgFg2dotflr1Woujrno09NtGKD80z+sUXli0BnOXrseH9h2J3Yb3C6Xpdt41XQ3T41u4yttMw2Z4Jo+2o+Yai7potVZ1Tx7S8BR3EtfIVPclH5J1bS70ZQ3eVstI2lx19ki5c/EmaYVyDAWp8lzyXvy8Sxoeflk/Eogzz/FnLW6ZWIkFQ5K5pvqCwGSTl4PqlaQkZ1z/DC6592Vh4rVCEu1+zaZt2sniJM+kwaJQU8jvtGz1o16u2bQNNzw2H2ff9GwkTWXy2plrEmFjk+eLj2xitkdt8uHz/KzRJq8pRhRkIrFB0Yxn7fYLlevLvWu4FsobPN/4Oc4zQpx4Bcz2VwamFfJyxxfV5NNTLZuwaeKVQa3JhyeXzeYaL03lb5UwzgKfzG5SuVAiPOLTPVNSTLyqvsnVG9uUoQ/2+9HDOPMGtedUEllssz2fqq4mbJKt8YOimaLJ8lHrRj8UdjWoD3NNxb8mNq3NPqoRIS81Gpve37SwSJWHWK80yoWpQXHkps/NUN/+x38xuE+z4elphLM08RqnnOg1+XApjEkdVAZtq1oKk6pNhDtAj3bNe05qrhniry7NC65Bq57Xms1h4aXDdtHU/j+erj0n7iYmkqQFJNHkbfO1aYK8ozQ9p+jqbVa4Vr/daPKc9k4vouSZ1z+Ni+6ao0wjC0e50VTMGvpydD7nFbtztF5SkkSoJjUj5hrNitc3VmzEswvW4MW31OsMdPeZRPYyZmGuCdnkhWvti7FG5dXEufOlJdjrkge10SBVxI3E+LMOC3lfM474yavzsEmTli3+JuqqAcl3/zkXgLCTlc5PPkPd4lYhJzHXyOYkFTaBAsMVsEkU9kJSZiOdq4ZHZV0J+W/8fTbO+U3UDhb21S7j6odfx4y33sNfnntb+ZLT+Mmv2xx2nYprPLIpJDzxmvzNq10o9eYaBnstV66OTguKy04n5OUVtxFzTYYPQaclffnWmdprLrh9NtZv7cCWBAun4rTcIBhcZ1RQe5utVNJ6IySV+a3yd94TeFv9ezXl2xmzGqoyQcuEY3blx61CTnK3Npua2KytCZcfn7Ayd6FPIytj1ZiIrQ8h77+yzds68fT86KYJof1COxneXlPZHu29zVG/VqPrIRC0DPH4cdc8ZlVXnWdO3OKKPEw4EZu85Reoa+ALVm4Mwi97+Zts8nptV57wZSysUdZy1eTHbn4WR131SGy6uIlX3kZF231oZahkuol733lrgJvbOyL1kKnY5M0UIbfEPDe2dWDxGv0WhzbtuhLf366yNvfEixVNr+Uyw3MLVgdtW1YGnJC3JO6dysGuxMfa0VnGE2+sCqWXg33JjUb1gb27vk1KYxZoEZt87MRr8sYQ6aw03jVx6NrhNdPfMIZMjquPfDywyUNvu05KFmsnA8NTb64O3EBVrNrYhodfflc58aqar2hXafKgsCmmrH7b4rGId1ZnGVN+bqdoqOCjFlM7i9sZqgiPFY5Yr7NvegZH/ETf8couzypkz7L48uMJTHLCt3zXzKX4yE3P4h7fPCiXVw1HqvoQ8jHnxQcpT7zahAOQG43NNnnyyzzw0un463OVzSDka0X7b14v3jQiYSxe48mzAcrxaETKonT30154Z2W+JEs1ssxp8cVEvE4qPn7zc/j8n2YoY+KI1/Bn3aGxyctzEGpzTeVYJwubdNZs3hbaKMTEmKnT8MN7Xw4d28pt8oaHzTsybYAyKSRxnoiPQzU5O2bqNFzmhy22MtcU0CHxUsV3vHmb14Yef32Vsjwn5C2R3+ki3yeeIz7Y9jILSQ3VhGXcitcfSB+IClE7ZYxhxYY2fOeuOVpN3jYGdxJMAcqSYNsOk4yoVMdFbU1cvJPNJp/+2hOuqcSr1wkuvv5CNRLrFJSBYCivWfQm29vl3OQjaze34xnBNGnjUeKV4+XzO2nTirYOG5u82VxTpCZvAw8ZrvACxdurN+P6R+cLHSM319hhM0/GzZV8FOptzOMd4+EOnE0+NeFm9z5p6zfxOa7b3B58mIB6EZF8RDf8M4YaZpX/1VvzhX+v3lQx91yt2McyTVuILoZKR8RBKGW71GrycuwaAH17VLx7i7LJx3l0bBK087gOSnWer6oERO8a0VwjCs2w8Fc9c7mEc1Ks4tW9Ax5p1KSFV8IaxJVhH9PIliT5qZSZT//+eVz5wKtY6YcaT+KRB9h1BiqT3EZ/NPjqO+tx42PzI/KmzBiO/umjOO7q9Ka2OOpCyCfRIJ98c1VIS1THYZfMNSnUQVEwqVazyoJL1ORNNuAkGM01CQSnbdr0mrxfjvD/Wft7G4sN69eSyV4zd+l6rRYme3RsbOvA0rXq1bJxdVdt2nLZfZXwCsGKV2Eov35LZdI/qslHhYFxUtRSWunSccFk0ljFOQQVlVFq9NwTb+hDfNt4MSVpr6rFUFyTDhyEhJxtSDLxKip13DFh8ZotuPz+V7FaGrHPXboeC1dtwuvv2pna0lAfQj7mvM0QVEQ+ovMaMdrkpY82cl6y2W9qMzf0NHJOLkO8j7b2Ml59J36PVFXh2om3mHjyunsoM4av3DoTD8x7xy+u4kLZIC35T8Miy07z9F8/hcOu+I/ynO5d8/YTv8guapO/RDD7yZukyOXFuWjaTiDG73NgkYnmNVcCzIVHJWs2bcMnfqvfGC5uQxzbevERt2rgza/ni754drbzB+b4SV7wN5UHlTiv49UjnM9HFS7feVMfQj5GhbQZgoaQDtl6oYTLjGryJRK1neJtc1HvmgpvWE7SAdkmPit5qH2/Aa+e9wiLkzxzBRcYdvZQE0vesxPyJm0qrgqPvLbCeL4SuyaakRy7Ru0MAJjehK0mH5fMph3qPgduirhm+htYu8Xf9B4ssPfrsNne0ebuuJBXjbz5ffEY+7oYUmkqMG/ZeixbuyVwzeZ7Vqzd3B7ZbL4W0xX1IeTjEhgerDLmi3RBmglLMVtumxU/5jID/vj0okr6AlwS8pp4tSXeXKM+rrI/B5OWOXQxK9bHb/kYh8ncAAC3z1CHn+WovGs4BHkSNqrJv7J8vVFZsW0/cULcJhtdOxI7GtGlOE6w2WykYSOMt3WUsWpjm9Jcw6/mQp6/hrw+u0OFEeCTb1ZcsuV1O7adcZ7Uh5C3tAWrgjrZTIraLJOW2Sr4HZ9wDR/OVZi1eC2+f8+84Hesn3yKthENa5A8D69sS5u8MY94k0f4mPd/2XCdLXG7Ntlw3i0vRo7NXboucizOh1ze/JwjHvVMW+F0F/9zbkgpkLHdcSgumZ0Xia4OqvzitfC8NHkAmPzj6cpOSDbXmCbMs5Qfh2rupmjqQ8grxMvT81fh9he8HY34ixzatyWSTrXS9OGX38Xds5Zi3rJ1GDN1Gu6b807iOl3q++yKbahElQ2bN7WFbXVxjW3VxjbcqdmsQGLuazMAACAASURBVEdWE1CSq236D119VEGbxOF0XD24TVRHyp3aYjnlV09GjqmiOAIsGLYrOxwKC1c5dg9n3rKKw8AZ+40KnbPVEPOwyeuFvPpBx3UcNkEDkzRl+R7XbW4PvnNZyMtsbOuwUvzSYtOh5U19RKFUNLpzfuO5l334gB2DD0a1IYLqhd72wmLc9sJifPrQMXlWMyQJ5Y/9jhfjBfgFt8/G6f7H/dC8+I4nGmq4OD95FpO/SZNXrZqtbGQS/4GJwk9FNYfITSWCvOJhS3sn5q/03HbbOxQ2ecjeNfGCMeJ9Yyk8WIw8fW5hNCyIjD4GkabMmKpZeZOpps40GV8z/Y3Q771/+FDw95Zt4bDaoedeZtjz+w/iowfuqCg+nzZ06bT4NTZ5k0mTJ6KziGgeEZWJaLJwfDARPUJEG4no2uzVzAZvDC2NDZFzpp41bsPopIj+0PKwfZkivrYJlflARh4a2uyYw2GMaaMK6r1r9Hzj77O1H6X8nBmELQkV7oRJucIfVVWDRoUmL77rdoW2683VhDV51R2bHADy0uSnv2KeQAbUnbIub5tqTZuzPDaN6omkUYq3dfqaPN8oR8iXt8M7XlwaLT8nPUEVK6tosppr5gI4HcDj0vGtAL4L4JsZ87fCdmPhHk12mjwnSahZG5JOfF73yJuZyhM//LdWbzKkjCItDI7FzlyjPt4uPWfGKpppuZx9Yc2mBNEks6JaOCe2MZ1NNnSPmtGLaaRka5O32f80Dt13YePEkBbV80jjdcUvkfc1Biqm28YG/cRtdySTkGeMvcIYe01xfBNj7El4wr5w4l0o02ny4mIVdbnxdRNzl0PKmpi/ciOuejDyaAP6tMRb2sQGvEGaA4hDF/rYuJoxtrNVX/jCojWRdPICqe6CSsiLbUxlk4+aa+J72LGDe4d+284tX//ofLuEBsTVyCIqm7zJTJcEVRY2/ZpsxpI30ha/Yd55Naq8c7pbQxSok4lXM4FNXqnJ67+Ohwx7elojT7xaXiZPzMp8cJ+RsXmEwtombKRptt+zfQ8yt/kT5EF5TDLXKOq+46CelrWqLiofbXFi8dHXoq6YXoAy0VwTrwF/8eidQ79tTHF5bTc3vrWP8rh2e8cc9GBlwDaLfGUzFmPe+1jghzYRR9fcrKaePO++xKqDRDQdwHDFqYsYY3dnrQARnQfgPAAYPXp0yjzM53lvrpp4zeLSlLRzT2Kt0fuUe9uF2ZQtCg7b4XylHPGH5nhCkmwQXhYnXtMXWRgDejVhrcK+GmeuUUWrBML3qAthIObcUCIcNHYQnlvojYK+fcd/reqdlSGGbSJVbexrf5sVihUl0qu5Qfs8ZFRl2jQnVfz2Kx94Fb95YiEAjSavMtd0xUZoSWyXxRibwhjbU/Evs4D387+JMTaZMTa5tbU1VR5xmzXzF6Qy1yQVflnwhuV25X3ouqeUx2959i3r8sQO7KW311pfB4Q1oNDmFoZr4sxm1osLmehdoxN4+Szs2rA13USYakgPqAVE3IYwACIdqepRVWu/WjN6BUP1LekEPKCfBP7yrTNDETYBnU0eaFW4RZuuKzPg2QUV86DYZvnEq7w3RHenLsYl8ROv3ptWTbxW0281jw17v3c3X0AVX+8sfvK67fdMnVSWGELRtH55UA/L8xJ4x10t+wzYoWs3SnNNTBuLbBqiMVHl1bFxeKzzJHi359Xt8vtfCcX9Tzoo1llK7529DB/9zbOhcAhq75r4SXnZXPPDf80LtZ2SQpNXUdeavAkiOo2IlgA4BMA0InpQOLcIwM8BfJqIlhDRHplqaqpHzPmKn3y+mrzNlWLjLFF+pgebRpelAxP9qecIKzsZsqyctU4p+DHHmy6ysDyh6yoATNl9qLbdKCdeLaSfHNBOeYWUdda2ZNKyVUw9cbfQJPGNjy3Arc9XNsIxzW+piFuJHHIzVmnyFmXI72lrexn/XVJpz+IjNY24arkFZVYyLYZijN0F4C7NuTFZ8k6E4Yv/75K1uOExz6NApclbDaU12GimYhJbW3ocW9s7U9kjkyBqQEveq4TfNZWbNtSwTNhcozNd1M52QUTaxUcqN9m4NkbSiledR0ok54xtqXdzss//qF1a8funFiYy15hIklyV9NePvIlVG82mlbhFYiWFuUZFNc26eVMfK14NUv6D11Zs27lr8gknP/MSSxu2dth5FmS4N91Q2FRu7NyIZdkM4dgiRWryaShRslW0VgG4hL83tnVgpmJDk7z7taTmPKKoaSmUX87hI8JuvNHzv7ZwB417T6KyIK/XEPnU7/Whkrs69SHkLRv/sH7RSZpMJg2Lj0TMXg4pm5b1lpOFWe4tbiPyNFh71zDRXAMou4caSvkSkVagqQRgnCYvx/V/ev7qSPRCIDxhqCsrCYmFPAhEeg+hRQkX3MUhauFp7zXuHkXrWpvhPfHwwd2R+ph4tUw3uE8LPnHwTqFjWSYnTVcO7NUUOZaXJrZ+S7tR2HK7cJbwxWm0srShhmUzGkNYi1NpwrXV5ClRZMtaBKWyIWm1iLw5jH/9d3kQ/E9kxYZ8vVI6mVmTtyHuNYlt1qTJd2fqQ8hbSs8SAQeMHRQ6pgv9aoOp4alMQ3EulLabk/zthcXBtmIquJDvepp89NiEoX0iH6IX1qDyW6Vg1dYmn0xAdlUhn9ScJz7x++fGx5vJSshckzaPWE3eu6sbH5uPS4TQ3/VEfZhrbNMRoUThl57UI0DEvB9m9Fxc7JqGEqFs0enIK0RleCTEvFwoRXTDZiL78BIiDSVSBu2St8OLlGcsqViSxiCyCaWbhqyT+L99cmHqa/nmG0WiC62RBJuJ13+/8i4uv/9VY7ruTJ1o8pbpEPVjzmSTN5xTZRtXz7x2bgo0+QyjFK3NWZOlzTeoStLUUFLsDMVCGphSyNfUJp8sfRdV5PEPi/DWIuIzl+cHiiBu4jVpHkoIOPePM9Jl3k2oDyFvqdeViCJbgxUV1kC9mMV8TV5CnsfeyKrJiz7QHFOOsesVFB+cLuJfOdYmn/1ZNaeMUZL2PaUtT0f1+478e9bzjx4f2QCFk8dEbuw30EU74DypCyFv2/ZKpegHWoTdGlDbAuPMGaqFNGnggrMYm7zBhTKm+qr6NKn2VmSoiiafdjMR1R6iNgzrb16C3xUR22RRo6fLTt8T44b0jhz/xG8rbos2r0p2qgAsdsLaDqR8XQh5a3ONahf3osw1Ku2TzI0qJxmPRl9wZvWTV5Hlk1A9E1XH5mnyld9FLURJm6/pPelkys5D+2BIn3yFfDXC34r3WoSMJxBaGhswYZg6siXHRhirRoVFbf3YndiuJl5LRNjaHvZKufHxDPG1jeaa6LFYm3zOmnwW4aiNlqjJ0pt4Neep0uTVEf9YqJPpct41KcRdU0MpN3NcNfGeMxP+zjt//3/DM/35Q6/hl/+J30BHZQ6L+wa6c0waW+pEk7drfISoV0AmF0qDlNd51xj923OeeE0j5PmGEOu3aFw0mfqDZJrjIsvWRheUqGJ3M3Rt7xqTCUB3pkT5jdSqSdGafJC3IXMbAQ+o21IeG5d3d+pDyFumKxGhLUfXL5MMVXrXxOSXlyZfIvKW3qcQ8v16eIu41ml2xdrWWU5tx7xWsZ2hLmSv+PFlmUAugrS2/Lw14Wo8FXH0UcRAJM8sUwn5nJ9imk/4YweNxuDezbnWQ6Q+zDW2E6/kRaHLi6R+8nEByvLS9AieXT6NMOrfswlL127Rhk444NLp6jItzDUqlJo8CwswVWelmq+tFmn7nG6oyIdGl3mHOvYyJfG/TKht8tVVEBpLJWOgMxVJdoxLQ51o8pbmGiLst9PA3Mo1+8lHzy5ctQmbDDG887LZEgG9Wxq02riJAX44hvc2b0tebuIrdC6U4XjqanNN7USmSXCYNMfuaJMvFexdk2eWKpu86X3sNrxv7uaaNB5yJSp2tFofQj6BJr9/nkLefy87DIjuN6qTA2Isa5nchDwIY4b0xsKVyf2MezU3YEifZrz0VrKdpNKi9K6Rnp1qYVYt5aVphGTy1sp79FENK1a15hHy6LSbFAqDbh0Mbz95P8I0Qt4UujoP6kLI21KUR8bw/j0ix0zaXu/maFwbID8hQASMHNATi9/bnGpieUifFkx/JYdNzC1Q+ckzSOaarjbxani3prUJeY8+qm2TL4Ig+xyKaVLs4az7DktEue3vICIK+QPG2CmUnmt1cdSFkE+iyRdBUq+YLBOsp+w1wipdS2MJS97bgisfSB6T45jdhia+BkjXUFXmGlnKm7ZlqwUmQZ6XDThu79JqQYVPvOaXqUph0I26CLxf0b+vNO8gtHjM1oxs2Dc3D+pDyCewyRdB0iGazqPE5j56NqlHAaF8iDItoR+pMD8VhdqFMuzzMGuxwnRUQ3uNSZCbRk5JqtxV3C1DLpQFPvM8clYpDDpzjafJm/Mb0DMaLjyONCuES1Tswrb6EPLWK16LKT+pmaVBc4FN/Ww6FIJaeNqSV3gFG/QulOZGX0shaBLySaOaHqsZNVmNDquy4lXUTPMnz29S1ea1oy6/XNPAa2Cv5G6NDSlGPqUSFRrErj6EvGU6ZZyUXMrPR5O3wcbUQwQ0K+yT1mWkrF4amdOgXPEab/rpqouhkga8u/S0STh43KDI8Vqu6BUJa/L55x+Y5HPIXDXxqtuVq0Tequ7/vLpCm19/xcY/cajMNQeMGYirztxLe017Zxlb2juViwXzoD6EvGUDUdp/cyk/WXqdpmyTjY2GR5RNk0872ZZmU3RVh2cjJmspBE2CXBUbH9DH2y9RJdZQ6HgX+TJDNvkCu9Y8cla1eV2bJBAWrjJ7n/XPyVzTUCKMGthLe80Tb6wCAHz3n3MTl2dDF2lK2bDW5HMO9RqUn1Dg6IS8zUYMduYaQnPqDo1Sm2vkvUptUHUozN+82/RYa6nnmnbl0nUAXtiHKESkVD5sOtpqrwMuRJMv2Fxz96xlqcvt05J8rWijQsiXiIwLpEhKnzeZpB4RnUVE84ioTESThePvJ6IXiWiO//8x2atqqoddOtVwLpfyLdPxKIQ6c8276+P3yLQx9WQ311RPhOqiUDKw3GL55M3ydVu155KGdy6R5wkl8vTUY+yEfBWkvFiNYmzyFCknLSohrzPH2DzfHhZODpF8FeYaIrOHWNHfW1bVdi6A0wE8Lh1fBeADjLFJAD4F4JaM5Rix1aSL0uRtFV8en0IWbKrFVDqqMfGaVwwdHaJrmkqQM+b9M9WjlvK/X8/kGt6wfj2UdSaiyH7Aw6W0XaavK6AepknspK68SZQ4m5Ry52viZ2ftHS2DeFmEIyYM0del4PebSeoxxl5hjL2mOD6TMcbHSfMA9CCimjv+8mHxbz81OSZlMmx74p7+IihZUHck8ciwKYooo00++TVJLhkoTGjpBDlj5nrU0iZ/y7kHaT/M848erzz+3VP2UB4nRIUJkZ1XS9rgWl+fsou6LspOKFUR1vBJbLmYNKPuxgRt3ua+bDX5X5y9D/Ya1d9YVo+mhmDydczgXrjxE/tb5Z0H1bDJnwFgJmNMaYsgovOIaAYRzVi5cmWhFeG+47sM61toOTq4cJdtsElWpR41oTVWCBOymWvSmEmSiBtRgJnMT6Z62NTwwDGD8L1T9sC3jt81Qe3i2WFAT5yp2bLuW8fvhkVXnBw61rdHI8Yqdj7itDTJQp6s/NPTuPgBwFenTMAH9x4ZOa5SVsKdTYETr1LZ3gKhZJ2YqS29ePEUY3kqbDV5ceKc13lgr6bI8+RlNjWUQjKo5uYaIppORHMV/061uHYigCsB/I8uDWPsJsbYZMbY5NbW1mS1Twjv6fN+prZaJW+Dsp98hz8pc/HJuxuvP2TcYBy68xCr8rIshiraXCPWX2d+YoxlNtd8/shx+OzhYyuda4739Z2TzO9K5MGvHQlA3zE1N0Q1RhtN/pqP7GNdB5mPHjhaUWY0nXioCFmk3Rg+xTjFZMqUQx7YNAW58zWVu+OgXujXoxFTT9wdt37+YNz/1SNDE69ymaryi5pjiTUuMsamxKVRQUSjANwF4JOMsQzbL+UH/8jzHurbyo5gkkk6zifr4jQz3lBLBMT54TQ1pr/HNJpFkivCjV0j5FPWQ1UnPiIY1LsZKzbET27bMDBB/G++SlnX7lTCxCacwOA+LRjat0V5T15kQ32dVALR09TDF5FFZ5MFXpp8j3xeJkngLlMnLq+RsdPk7cw1jSVCj6YG/PeS48NlBGWF/wfCbZubrIqaRy/EXENEAwBMA3AhY+ypIspIA7dT591Yk8bOkbVsLuTjTCz8+rhhMxHQpyW5jy+/tqD56QCxgevig5QZw46D9BPSNqYD0U+5lpjaR8/mBqVZILwjkz4D/U5U5ntWPpJYTb6456h0pYV+3YEK0+Il2USa58Srbv2NrNSJ71G8Xf79FxXaIKsL5WlEtATAIQCmEdGD/qkvAdgZwHeJaJb/L13UqxxpCDT5fPONEzh//8Ih+M83jgoacktTCU/+39HB+U5LIU+RP/Tp+vVItx8MofhJTf4e/vGFQ7S+yGUG7DGiH847cpzyvE0VZSFfq/2lAlc66fhL330/ejQ1KDXGzKOYmMtVprA4wV9Iq/AFm1w2gyf0kqwgHto3Gg2W01giTPvK4RjpR4wtyiYvUtHko3JHfL/8+++Smjxj7C7G2CjGWAtjbBhj7Hj/+I8ZY70ZY/sI//Trh6tMlgkk1Sq4uNWJQ/q0YFxrn+DFEnkr4LhGz19yXKNS2fZ09Ekp5IH89prVIQ5fTTZ5Amk9HGyEIH/PRc8xWFTE+0+qRu8W795iNfkU1Y9r46p33LvZ3GYKscn7/8vvk2u1ScNEvHHpicrjRISJI/tjx0G9/N/xedk6L2hXsAcj70od5HNAZVVuUfFr6mLFq8ykHfTuTECyxvqXzx0U+n3beQcr8rP7oHTJzjnImwTbfUQ/Yz6i320cvWI+WBOFxxD3/2dMP9Qt+ytedR1OkirWYlHVHz97YPC3rnheL1UcdNsQv9oRfswtq0LwqkxnIXNNFb1rPE0+mbkGiF8fUpkMjc9LNcI6aGw0zpDe3bOi1FV+8fIrv3hH1iXNNV2VWBdDy7a6/04DcdjO4UUMKs0yLjsehKskvWzuP/DBvUdi0RUnBytidSRZHTi4T/qNgYuOmyJ+0LqhLmMMRMniDf3u0+H1D3xCq+g5BhVH7dIaaIK6qUveHlQThrY2ed0gP+6pbd0WnbofoLBph4RvEZq8X335Eew5sj8YkplrbKh8gxbmmgTeNSq2+mFKVm9qC5UNhL9h/q0WpVzVpZCPk4K2GomqZ1W6mVlOcgUTp8Tz9/4P5goU14o+3pV84uvfr0eTcZWdiWppvgwGcw2vSwJTyyHjwvdbDgSIb5OvslG+8r7VnbNpYnhEwpj+Pzlzr5BbZNwr3KwQ8jsNivryFx2Fkis6Ypvee1R/tDSWPE0+RdA7E0k0eVs3ZN3o4ck3vcBjM99eGyrbK7/y43efPgAA0C9FQDQb6lLIn7HfDsbzto1VJRPUC0bM+XCbu044lwwa+nlHjsOBY7whom4CT0drzMhAR+F+8sLfOrc3z3WO9BusKB6WfKiiydfGJi8rE3ItyKDJqxYrqeAd1747DsBxewwLjscpAuOH9gn9/sXZ+6hdOQsOBafT5Im8c6pYQHHmWBMkmVBMZLXJR8uuIL6fYf164A+fOQCfPWyMVT5JqUshv9Pg3vjJGfr4zVmaraptxOXHzYo6P3lZ4xNpIMIg3yebWzbi2mdFc0x3p9UKUCbb5N+89MRg5MKjUPK6jGvtjfu+ckSQVlVD3eSd6namfeXwjLW3J9acpxASx08cjj985gBl+v87YbdoGWRRkMDYIb3x/Q+oQy2IiNETCflv1lLxkw8PGQjkm2uimnyWTjtJ07YNDWK7yE68R/mS9+06FPuOttsTNil1KeQJwOZt+nCwtsJPNbxXCUBTfnuP6o9RA3v610r5+/+bGi2RKNyTCe+0srpoxbdirmIhm3xjQyl4vmXGQKh0AuUyC9mMVXWMPF/JXMOf+PB+PTBxZHpt0JbgPmPTqR+4LnCduMmIUkjCzivFZivJTUJYZVUwtbyIrPD1NXlVyI8s7ZM/p7h5p/OPHo+hiolo1VO1jZkjepUVPUISqUshXyJCmyG0p+3jVb1Q5UpBQ4Y/OHXPwPxRkoQ0Y2YzzuePGAvytRqx3rbCO1WgMSrevCE2cLksfm+Bd41/vtOfiK2kU5nNwsf4SJ8LszWbtgGo2IFVE415EjXPqNPpHrdO+CtNVVJ5pvjlQbkW71m03ROyxURSwTtiSZEP1t6qgvdlaZ+mBYVi1MtvHb8berc0Rs1mklA4YsKQQImLI1A1GEBVlLx1KuTNEzZZrBGqa01aYcgOV4oeA/SNduqJu4cu4BNBtuaUtGaXqplrEB3q8l9l30+en+/sZKEPc+SA6MIXudpcmB88fjAARDylHv/20fjRh/bMcAd2xO9XqxPm6vTi4YpJikLCX6WFyoT3I1UXJo6IiYAbPl5M9ESVTR4aTV6u61ePnYAFl51kVQ6/Uizvgvd7UTlV72l4f/0CK8CLSGpr1hHfczX3bKhLIQ8CtrabNHnLB6x46SrPk4PHDcKfzz0ochyQNRR1uTohL3vdcNOFbfNIb64p2l5T+TOy3FwycXDNcVtnWJP/wlHjsdPgXtK1ak2+X48mvPLDE3D+0Tt7eQvH5Um8739gD1xx+qSEN6RmmLS6UtfudIpp4BVksRZS1uT79GiMRMOUsdGIxXlPAuGQ8YNx6Wn5dYwq7xoCsKmtE88vWoO1m7dFrlF1CLbOAvJoGgB6+6uuVU85zy9BNNdU0xegLoU8gdDWYQjhZSvjFcd09jddjy9+2LIXTWW1n7keYohS8bc2vXRdUrKaaz55yE5W6RjT+8nDN9dwG7A8MmsslXDiniNi8q+8wZ7NDcrOS+60P3PYWOwx0rwoTeawnQcrj//1cwfjZ2ftHYRu0GrmOk1eU56YXFwxmvR1mwTjlN2HYcruQ3GGEFJZNKXlRmCuCddlztJ1AIAbH18QuSSLPVsUtBz+GJQDrhxH/ToXyqKpSyFfitPkNc/3pe++P/Rb9dL1Ln129VIftxPaTQlDJaeV1WmvO3KXVkwc2c/6fhhY1Cbvn+UTr9xEta2jHF59SfHPQX5/wccs5ROtX7IHwNMPkRagDe/fA2fsr447L6Jbl6A14yjql3UyT76ytW8Lbv7UAer2XsCCA138fJXZVR7ZJBshRUfDRYlb+b3armLOm7oU8mOH9DZq8vLH87nDx+JLR+8cuCqaUK3A1G3SDOh6bwquA+I1Z35ZYynaQFXwHajSfvRp/eT/9NkDMU1wc9QhOrtEbPKBtsg8bw7fd7u9sxy6cdk8oaIsCSPVh5V0Il0FT/+RA3a0SiejG8zY1EO8RTG9HFpXRbvknBB+XN6PdVvaI/kXoMhrOzS+kEhEHknYKm1eOfwaYYQtBLA7cpdWjBM2eIl8QwnaRlSBEf+unpRPH+Cki/LbT03G0H49EnnXfPSg0Rjf2ieSTiXQVeaFzjKzC5iVUpPn5yvx5M3ph/fr4aeLrZKSwgOUofJRmWLXABVNPrIoxkKTjwqDcOcK5Dtsjg0BrTmvPa4176iPiYd58DORq87cKxQ6w+ScwJ1a1m8VhDyizy8rqrUMsW/EsnzdiAeQw0ZU+JMQc0iuV5KyAU+BESP9m/zki6TuhDznzP1H4e5Zy5TnIrYyRZr/OWocPnFw1Las0vxk9z5dWbqVrXFhXwNzjWWo5BEJwqnKECjdpiEJhqLiebnTlPNpEXy5Q/G4QbFCNaLJK9KopliSa/KWFyQQ2ony5WmF5L0VIZzPmhweachulmFbv/fsRCHPkZ9rFpjCJh9323JwNe0oWvFuVbFrZJfmvJC/a5vNcoqg7sw1/NkdMcF+K0HVA7/wxN0xamAvReooPCyusj6hiVd+LIxKcxaFH2/UTZYulH17eP7fcrKBln7haQKU2XwgPzljL/z6Y/sJ15iHtEA4fkjYUyl5neRFUV4+2T82W61MFwtFH+4iPs/KWotwW7NZ6LRNGO0ShTV0/vfPztonlEY8lwcq54OknbcOVS6qidfdhnv7rZ6pmD/J0joie7wG8wHJJ8mzUH9CPsVrydqrdpZtJ141H7PiLYjCj5sqbFfWNUpRLzkzv3ec1fW2nVtSPnzAjjhp0gjBtsu0NnnAe5diPJWwTTP5xCtHjCSq6mCTCjHb1qMLSauNR25hxmFCWvG4SpOXMW0gz8/sOryy4TTPPk9NnpPkG5SLT2IGVa0aHzGgJxZdcTJO3Sca88pUrQ9PNk+qR5UM/zhYLsqFLXUn5HVf3GDDpGrW591p8CkLa59qc4tSkxcEAo/fwYVEXH3z2OawV3Oy5eup4+TEhHTQab/eSmAz8mvhIV93FoJzqQRBWiFm+15sr0tit5Vt8jbvb5th3kr1CFRzGjJfmzLBqt77jh4Qyit0TWznLU+ox3eGctY6m7wtp++7A35y5t7GNJHHVEXtXaTuhLzqOV7ygT0w/YKjtNeIgqZfj0b8j2bLOR1lxiLCSuUbrWv8Kk1O1HB5J2JrrglWxtZ6RyQNpgk8+VsPafKi3dY7EPxWLUaTXe0OHjcYJ+81ApedVlnsJI6iTpg4HEByP3DbDk4XEkD7PrWHFRqqVI8vHzMhtj4fPkDwgZdGAqoFWGQ4x/nalF2w4HLzIiyg8n3wvJJsGG77fsyavHhMn4f8rHsmUX4ivgLR+YBqUH9CXvHGJgzri4GCJi9rh+IV/73keFx40u6Jyuwss8hrU7kx6jYsUEafFKQPH1Y32mryjV6CLYqY4bYUMSTnGD8qaQJuSG/PG2R4vx6hp1YqVZ7i2QfsiMMVsfNlYdCruRHXnbMfRgqBv/g7GdKnBTd8Yn//umT3TtL/OpLb5NWdsN8GpQAAGWhJREFUYSi5MHHJj+8xol/scnwAGNG/J07cc3glK7EcpSavSJcAcTQt31sSc02WtqmKXWMSunK1TtnLW4BnUwM5Tb+eXse2amNbNHGB1J+QVxyTNeXGhhIWXl6JdZHVJl9WeNfwia+0mrxov43T5GVbL0+3ckP6xhS4MKYISGW9KYvqWsm8VSoRnrnwGNzzpcNC5xpKFaGmG7HYTAYHm3wLaVN7WcS0I525Rqw+H00AdqP7ShRKobNJ0Zzla8QnsOuwvn4aL5GNJv3vb4RHzrsO64sXLpoS/FbGqlGUrcLWT16pyfv/i99ckueVZDW43I74c9x/p2JCCuuoOxdK1QtTTXjl6bPKg2mJcJuoGPI1iU2+QWmu4SOBMAePG4wn3lgV/A42CM+g8fAG2tJQMtpu80b1Kkb09zRvMY6J6E2i8+u3CRzFLxW1w6SbEWU11/A5guvO2Q8n71UJ1RC3o1TomOZ4WkQBdet5B+O1dzZUzinE8JjBvbBo9ebgt2pCvRQSrOHzYvK43aBsO2F1OGrvoPguTI9NPCdOtCZRBP73feMBeMrlc985Fv16FBv9VKZuhDx3AVNpkQ0xPoG2H6muF999RL/IB8Y9OLYKK2/1fvLRPEM2+WCHo3AYS37P8vCVL5oyTQjrqAg97//mxhJQ0OhS+aEYbLMq0xegfi9fOXaC0iVOpkGhnSY211gKVy5Y9h7VP/jwAWBAr2ZlMDFtCGLlIh8Cf2Jiff7vhN1w5QOv2lVQQHwCg3o345Dxlfg8qh3H/nn+YeifYPs6+d7E92kKSQKo1j+oH5Tyu/YPhUxnFu/vf44ahwtP3B13zVwCIJm55sPC+oRh/eLNaHmTyVxDRGcR0TwiKhPRZOH4gUQ0y/83m4hOy15VM6q4LvzvuJ1bbD/SGcJwk3PH/x6CEf17as014obJWnNNjCbPWws/xD+IYIMN6ZvQrhJNABfAXDAl9bYxwTtA03A6+iP8WwzIpXp/F7x/FytNnj9nUXCktcnHweuz7+iBOCEmuJqXb3zOIT95hb1Z7EySYGo6Z+w3Cteds1/gIQPwOQE7O7dX37CZjEJC3jyXZD/xGj3Gv42wJm+oqzSPlmTVr25rw2qT1SY/F8DpAB5XHJ/MGNsHwAkAbiSiQkcNLYoPmn9UuqXznDib/BEThuDwnYeEJm8rhF8+ANz4if0xZojnay4Kaz5clRuVqnxRQAV7lfrp+G49XFOXzTL82nIGIS/b5D+TYP/JPool9SKXnz4JXzhqfCS+u0x0grryd4Mw8ZolDEMgXBQLgezz8P+PSTfBN8vsMqxvTEo/v5K6Plo/ecX5OHT3ajJHlEqEk/caEZ4ItzApqc7zUsT9aeMcBmxiEnnHoyf4aEO3yC6Sh0UaHaowyrUgk5BnjL3CGHtNcXwzY4zvNtAD+cY0UsJNFOLj5C8yTpOP62lvOfcg/Plz6njxlYBH3v8tjSUcP3E4vv+BifjZWXvjwLGVrdrkd803dlBNHIqdQ2VyzTu2fN3WUFpZmHPbvWpXnaSkGRV88eid8a3jd9Uu/hnSpwVTT9wtNjiY6dsQNdcsrqL80myavF35R+7SimlfORwfPdAcyKySb4LjZG92FDnnoNEAPL/14yZWBG3SD9a2bHk0yh/1joN64fmLjgUAbNFo8nw0KYcbSXLX/f1V36JiZHO9LOytzDWB55N19QqhMO8aIjqIiOYBmAPgC4LQl9OdR0QziGjGypUrU5cXCBThgXINPi+bPOdj/ochXivn0KOpAWfsP0qa4A3X8a7zD8O15+yrLEPsmFhgk1fXUxRKh+08OIiHnsYmL3PVWXvjmN2GYo8R9jHWezQ14Pyjd7beMUck7NomnxP+FswDWTSlysYcFcYpgtWZSFL8xJH9EwhEG3NNpQ6B6cO+Ojhyl1YsuuJkjOjfE0dMaMWvPuq3RxtzRKiuduWZ3lnvZm+wrzPXDOjZhEVXnIzPHDbWrjBlHt5oXBwtmN6HfGpPf5MZ0fVURyVkQxfX5IloOhHNVfw71XQdY+w5xthEAAcAuJCIlDMOjLGbGGOTGWOTW1vt483I8FgvopBostTkk76DS0+bFNjcI5qJ4To5ds0OA3rilL1GKtOKIQxU8T1EuANPa98W/OVzBwcbbWSZeOVMHNkPv/v0AYnzAdI1blsfekDs/BIXEyBuHM7ZYUBPzLfcTs6rV/ry0+QbNtdUTAJMlSBlmTa7UYlE4rTo6m4oh39Th2vMeA2akWGS2+X7+m4StzU0pJdHaeNb+2DBZSfhpEnxcyqVubTaCvlYOzljLDrbmADG2CtEtAnAngBmZMnLBLcdi89Tt/hEJs1LkO1tQRaGb8OmnOMnDsOD894Nm2tiGgtjDPd/9YjIvp5ZNPmTJg3HfXPeCTpIU+hmHbaP9YlvHx3kH9bWzddxT7ssH5HsTcRpKBHOOWg0Dhs/BOf/9SUAwGWnTcKClRujeSi8WvJAbwZSmLmMZ+3RLcBSkcZD18ufCesThHMlwuPfOhqtfVuw+/ceiFyr3UUsAQN8m3xYk4+/Tkxjax6syAj7+hVBIeYaIhrLJ1qJaCcAuwJYVERZHJUPOTfXxPndpnkJkckwi09LFRxJhmv2KnONrnGVGcPuI/phsOTalsa75pvH7QoAuOYj++KFi6YEdU3jK2/7WHcc1CvwFQ/Z5OWJV+m6ThbuaNPAR3tfUISyuOy0SSG/9XMOGo2LT9kjmkmVNXld2jwWKQe6SsK85Hegmw6Km6QePbhXJHTAT8/yYsRo4xgleAH9VZp8AnNNGqoZjExFVhfK04hoCYBDAEwjogf9U4cDmE1EswDcBeCLjLFVunzyQLX/6dQTdkNTAwWLaXSk0+Slay2GudwrRqUNcgK3T4O5Zqywcw2gX7yTVJPv37MpsEc3N5bQKowM0mjyaSZEw3G+pXPSbz7hnGVP2oYSYdEVJ+MCv3NLQzU+YT4pCUjmGm6Th6chy+eTksRcI6aQy9QtxJPdf21a6KiB3vcrxjESSWauUdnk469LE2/GxoXykHHq/YHzJJNbI2PsLnhCXD5+C4BbsuSdFJWf/ImTRuANG9tZGqSZc5uG8sirKwAA85at16YJVqsKqpDsQvnHzxyII696pHJeI8yTCnmTsOSa/EmThmPskN647pH5ofOqK7MKv+jEq6QtxkxI58V9XzkCzY0mbY+PIvOtB7+vE/YcjqF9K1NaymdNosCPcrzgOWPG3lwjJpIVJV3bq3wvyc1CPRqzr9XgLpQtjeJmNPGkc6H00CmRr/34hFxMUHHUzYrXinBMPmbNwyZvY8uMW+QBAEP9FXHvbarsyCPvniO3C53LX9JnYZKVh+88BFcCOPfwsdh/p0ERIa8qKZUZxWCU12mL4vEbPr5/EAgqL7i3ko6iupimhhKev+jYwCNEheq5y+aBl394vPX8VBIXwdB10m9d2+NtolKd+JL4DlaiJv+Zw8bg908tSlhLL/rl907ZA+/btRXH/OwxAMWZU1iMObElh07LhroR8jzyYpz9XUUWm7xsYzQ12REDemKZ5OMuM6SP90Gv2VSJ0yJv+C1rrrqhscom/6MP7Yl7Zy3D84vWRM6ZhPKkUf2VS+9NpNt+UI/sJRWYa4RyTrBwbcubIgcSogbP0T1XXdvr1Wz/mVds8vbmmgPHDop2wHGafIKusc1XjlqEVarf/8BElMsMf3zmLet8OJ89fGy4Toa0WeY5gvUtNQ4DWTdRKLm5psOw242ObDZ573+bLK4Xtr7TwT/qzwsTgXJZcn11DVGl0X7i4J1w+xcOUabP5KWiOPbHzx4QWlNglY9QBzlPuXPjcqTWLmrVJnS3ipW6WZ5GEjMK56KTdo90PLrRJU+V5JXx+SCd5ptVE7ezyScnTRjlIqg7IS9vTmyCx4ZOZW+TYm7YaCbcFDNEEeCJ09xYwqIrTsa5grYRKUuegNR8UL/55GTlcR1JbNu3nXcwLj7ZHHd/4sj+uFTYoMMGg7Um8jFzbbGrbI5SrW9ZfE8XnrQbAM9cyRTmq6TYjEg5po5g9xH9cPYBOwZtJGpytO9MuJlTnnjVLURMStGbeNS6edaNuYbbHJOYa37+4X1wyQcnptIEKjFD/P+5LTOm1T499ZjEwb4i5hpLrSnOq0gmyWM4eNxgHDxuMEYP6oXzbnkxUTm2dYj7+ILFUDX+iEo5CRtb+gsbsn/msLGRFaBZhJZtO1ZdI9JQIlxxxl545LUVoeO9mxuwbkt7YGpLMmIoKkSvqd0nXRSmwmnyOWHrEy/S3FgyatUm5KGY7YTVyAE9Azcu67KkRRW2PslJSeOlwoO27TAwWYeiIyTkY6pTCcFc24+o2n7QfTQ29hzc5JPFZrFIJT8Z7qo6yG83Nnmcus8OOP/o8fjm8Wo31zzuO5YcFtzVirrR5Lm5pr2jKq88IM1EUlJk27Ms5D8uBWxKSxqN44Axg3DDx/fD+3YdmksdTLFrZIIVrzUW8nyyvFofs34nLP+PTOaa5K6N5u3zwufO3H8UTtlrBO6etdS6nObGEr51/G6KvNXp7/3S4Vi6drP6pLKO+nN5LDCrtSZfN0J+hL+nZaKNdnMgoskX0MfIrliiC2VSjxcTH9hbHUcnDpvY6EXAh/w1/4hqbS/yCTbFzpJJChdK0+Ov2PgrOfZoaihUKZo0qj8mjepvnd5Ulxz6zdq3z5qWniNfPmYCRg3siZOLWvwUQ5GvkTc0bpYootHMueS4IApgLQlvEWdOKy8SqzYfnjwKr7+7EcfsNhTXPTIfR0xIH2DPhvu/egTeE1xrdeQy8WrjQmnRE3C3V92SjTx0otR78vokjV2TFDfxmhPNjSV85IBk7np5EGwEUqSgkZZHF2GD7lvlfSd18EiEQPwzVS2GqiY/OXPv4O88R1Q6do8L95xH7JqcHyZvqxG/+RxGvnmNBoy55DA079axaxziZGhxZZRjXCizMGmH/vjmcbvkl2FGRCEfRx6xa+qJPNYNVDT5+LQ24q8x2HQm7B2gMuMkxfY2dxiQPXZVmg5lyu62oSSKpW40+VpRCTPr/b/nDvaba9hSpLnm3i8fnlteedCz2W5bNmD7XQylI7DJZzItUCgvG0zl8eiwndIixTy1W1OH9Nx3jo11WS6q+Vz3sX2xdnN7fMKCcUI+I6ISecf/Horxrb31iVMiR7OrlQ26GvQQzTUx2lNXcaHsauTjJx+f1sYWzgNwacNeZzLXxDOsn3KvonA+hu8pi7GmpbEBw/pV1xFEhRPyWRHax/47DSykiMBcU9DmFF0JZ65Jz/47DcS+owfgOyeZVyKbSGKuqVyjf/6BJi8J+Yq5Jjt5LFiKozt/c07IZ6QapoJIMLQcypz53fdXZxFJQsRgWrGLoQIXyiJr1H3o1dyIu754WLZMAhfKfMw1jZqN4G2b8F4GV8juLHiriRPyGemu7YyvVO1qhLxrYtI6m3z+JFkMxcMMmPZQbtR51/jEmXzu+VL8nFERa1PkvLtzC3PeNRmpjiaffkJt39EDAAB/10Se7Go0NyaZeO0ai6HqiSQ2+Ws/ti8uPnn3YOtGFU2afR7Sxq0P58EniYsjj8nsWuM0+YxUQ8Bk0VhvOfcgvLt+K8a36j/EroS4mldl6338W0cHE67V2hlqeyKJa+PQvj3wuSOie+OKNGps8pN28Mwwx+2RPv5/Nd96kaOFonFCPitVaGlZtIk+LY3o000EPIDQdmiq+x09uFfwNxcc3VnL6mqkiSdvoqmk9pPfeWhfvHHpiYGmn4UiBXCDX3/dxjzdAWeuyUg1lMiusvlANUjyPJ0mnz95mFFEuCavMslnFvBVeO3Nfv23pdjIvqvghHxGqrFkuVwHkz+2JHmenYrt/xzZSBK7xoauErwtLXyOKM22ol2FTEKeiM4ionlEVCaiyDZERDSaiDYS0TezlNOVqYYSyfe2rHUMjGoTd7/lLhJquJ7IezKzSd51vgCK9JMPQpin2Fa0q5D1DcwFcDqAxzXnrwZwf8YyujRFbx0GALeceyC+NmVCELd8eyHuyX7uiLEAgF2H9S2+MtsJeYfMTtsB89DhJvq2eFOK6woMHZBmW9GuRqaJV8bYK4Ba4yKiDwFYAGBTljK6OtVQrse19sHXpnSdIGLVIu7ZHjdxeFWiP25P5LkSVeSkScm8aP715cOxbO1WY5pxvkPB22vsNwhJSrCtaDe2yRfiXUNEvQH8H4D3AzCaaojoPADnAcDo0dUPFZyV7cyCUlXco60+lLcqD2DuD45Hj8ZkRoPBfVowOGZrzmN2G4rT990B5x+zc5bqGWlq9J5Hd7bJxwp5IpoOQNUNX8QYu1tz2Q8AXM0Y2xhnV2WM3QTgJgCYPHlytzN8bQ8eL9WGqHv7JXdnitDk+7QU46ndo6kBP//IPoXkzQls8roAa92A2KfPGJuSIt+DAJxJRD8BMABAmYi2MsauTZFXl6ZWQv66c/bD6EG94hN2Q0pE6GRsu5to7grwR152vSwAce/oOtbk08AYO4L/TUSXANhYjwIeqJ1J4eS9arPNYTVoIEInmDOF1YAksWu2B5rrYOI1qwvlaUS0BMAhAKYR0YP5VKv74ASRHYeOH2ydNoi2WVBdHHqK3JC+OzLS31Vq4sj8NwOqFlm9a+4CcFdMmkuylNHVcSaFeOZcchxaGu3jxAcrWN2zrRlOxnvsOrwv7vvKEdhlWPcJDSLjVrw6Cqdvj6ZQdMk4+ApWJ+KrT0WTr38xv8+OA6zS7TGyXxAXvzviApQ5uhxuBWvt2J68xW4772BsauuodTUKxwl5R5dj3Zbab368vbI92eR7NDWE9hSuV7rvGMRR97yzzrzi0ZE/gXeNs8rXDU7IO7os25HloMtQ8ZOvbT0c+eGEvKPL0p1jeHdX8g417Kg9Tsg7uixb2ztrXYXtjrw3DXHUHifkHV2WNqfJ1wBnI6s3nJB3dFmcuab6UBERyhw1xQn5lAyJCYPqSM+vProvAKfJ1wIu412AsvrB+cmn5P6vHuFc/AqCb3fY1uFs8tVmp8G9cfq+O+Bcf9ctR/fHCfmUtPZtQWtfp80XQa9mr1m6uEDVp6FEhcdod1QXJ+QdXY5Dxw/Gl47eGZ88dKdaV8Xh6PY4Ie/ocpRKhG8ev2utq+Fw1AVu4tXhcDjqGCfkHQ6Ho45xQt7hcDjqGCfkHQ6Ho45xQt7hcDjqGCfkHQ6Ho45xQt7hcDjqGCfkHQ6Ho46hrrQ5ABGtBPBWhiyGAFiVU3W6A9vb/QLunrcX3D0nYyfGWKvqRJcS8lkhohmMscm1rke12N7uF3D3vL3g7jk/nLnG4XA46hgn5B0Oh6OOqTchf1OtK1Bltrf7Bdw9by+4e86JurLJOxwOhyNMvWnyDofD4RBwQt7hcDjqmLoQ8kR0AhG9RkRvEtHUWtcnL4hoRyJ6hIheIaJ5RPRV//ggInqYiN7w/x8oXHOh/xxeI6Lja1f79BBRAxHNJKJ/+b/r+n4BgIgGENE/iOhV/30fUs/3TURf99v0XCK6lYh61OP9EtHviGgFEc0VjiW+TyLan4jm+Od+SUn2xmSMdet/ABoAzAcwDkAzgNkA9qh1vXK6txEA9vP/7gvgdQB7APgJgKn+8akArvT/3sO//xYAY/3n0lDr+0hx3xcA+CuAf/m/6/p+/Xv5I4DP+X83AxhQr/cNYAcACwH09H/fDuDT9Xi/AI4EsB+AucKxxPcJ4HkAhwAgAPcDONG2DvWgyR8I4E3G2ALG2DYAtwE4tcZ1ygXG2HLG2Ev+3xsAvALvAzkVnlCA//+H/L9PBXAbY6yNMbYQwJvwnk+3gYhGATgZwM3C4bq9XwAgon7whMFvAYAxto0xthb1fd+NAHoSUSOAXgCWoQ7vlzH2OIA10uFE90lEIwD0Y4w9wzyJ/yfhmljqQcjvAGCx8HuJf6yuIKIxAPYF8ByAYYyx5YDXEQAY6ierh2dxDYBvAygLx+r5fgFvFLoSwO99M9XNRNQbdXrfjLGlAH4K4G0AywGsY4w9hDq9XwVJ73MH/2/5uBX1IORVtqm68gsloj4A7gDwNcbYelNSxbFu8yyI6BQAKxhjL9peojjWbe5XoBHekP56xti+ADbBG8br6Nb37dugT4VnkhgJoDcRfdx0ieJYt7nfBOjuM9P914OQXwJgR+H3KHhDv7qAiJrgCfi/MMbu9A+/6w/h4P+/wj/e3Z/FYQA+SESL4JndjiGiP6N+75ezBMASxthz/u9/wBP69XrfUwAsZIytZIy1A7gTwKGo3/uVSXqfS/y/5eNW1IOQfwHABCIaS0TNAM4GcE+N65QL/gz6bwG8whj7uXDqHgCf8v/+FIC7heNnE1ELEY0FMAHehE23gDF2IWNsFGNsDLz3+B/G2MdRp/fLYYy9A2AxEe3qHzoWwMuo3/t+G8DBRNTLb+PHwptvqtf7lUl0n75JZwMRHew/r08K18RT69nnnGawT4LneTIfwEW1rk+O93U4vGHZfwHM8v+dBGAwgH8DeMP/f5BwzUX+c3gNCWbgu9o/AO9Dxbtme7jffQDM8N/1PwEMrOf7BvADAK8CmAvgFngeJXV3vwBuhTfv0A5PIz83zX0CmOw/q/kAroUfrcDmnwtr4HA4HHVMPZhrHA6Hw6HBCXmHw+GoY5yQdzgcjjrGCXmHw+GoY5yQdzgcjjrGCXmHw+GoY5yQdzgcjjrm/wEDSZxbRo2ADAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#Trajectory looks OK:\n",
"plt.plot(trace['center'])"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x1393edb10>"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxV9Z3/8dcnN3vIAkkIISsJYd8Nu+ICIrihdUOlaq2lVpm2v047tfOb6W86ne6dn7/p1Eqt40xbbKm2WBGQtSqyE7awhEDCkn0jhASy3/v9/ZGoMQZyEpKcu3yejwcPcu85J3nnSt6efO85368YY1BKKeW9/OwOoJRSqn9p0SullJfToldKKS+nRa+UUl5Oi14ppbycv90BuhITE2NSU1PtjqGUUh7jwIEDVcaY2K62uWXRp6amkpWVZXcMpZTyGCJy/mrbLA3diMgiEckVkTwReaGL7UtEJFtEDotIlojc2GHbORE5+tG23n0LSimleqvbM3oRcQAvAbcDRcB+EVlrjDnRYbdtwFpjjBGRScAbwJgO2281xlT1YW6llFIWWTmjnwHkGWPOGGOagdXAko47GGMum09usQ0D9HZbpZRyE1aKPgEo7PC4qP25TxGR+0XkJLAeeLrDJgNsFpEDIrL8al9ERJa3D/tkVVZWWkuvlFKqW1aKXrp47jNn7MaYt4wxY4D7gO932DTXGDMNWAw8LyLzuvoixphXjDGZxpjM2Ngu3zhWSinVC1aKvghI6vA4ESi52s7GmO1AuojEtD8uaf+7AniLtqEgpZRSA8RK0e8HMkRkhIgEAkuBtR13EJGRIiLtH08DAoELIhImIuHtz4cBC4FjffkNKKWUurZur7oxxrSKyApgE+AAXjPGHBeRZ9u3rwQeAJ4QkRagAXik/QqcOOCt9v8H+AN/MMZs7KfvRal+Z4yhqdVFoMMPP7+uRjWVcj/ijvPRZ2ZmGr1hSrmLwup61hwsJut8NceKL3GxvgWAQH8/piZFcfPoWO6eOJzk6FCbkypfJiIHjDGZXW1zyztjlXIHO05X8fIHeezMu4AIjBkWwe3j4kiJDqPF6aK2oZXdZy7w0425vLjlFF+YO4IVt40kIjjA7uhKfYoWvVKdHCu+xI/fPcmOvCriI4P5xu2jePCGRIZHhXS5f0lNAy9uOcVvPjzDmoPFvPZUJpMSowY4tVJXp0M3SrW70tTKi1tO8drOs0SGBLDitgyWzUomyN9h6fjsohqee/0gF68085snMpkzMqafEyv1iWsN3eg0xUoBu/MvsPDF7by64yyPzkjm/W/dyhdvHGG55AEmJUbxl6/MIWFwCE/9937ez63ox8RKWadFr3xaY4uTf1t3gsde3UOgvx9vPjubH9w/kciQ3o2zx0UE88aXZ5M+dBBfW32Y4pqGPk6sVM9p0SufdbbqCvf/ahev7jjLspkprP/qjUxPHXLdnzcqNJCXH5+G02X4uz8cpMXp6oO0SvWeFr3ySe8cKeHuX3xI6aUGXnsqk+/fN4HQwL67NiE1JowffW4iBwtq+Pmm3D77vEr1hl51o3xKY4uTf1t/glV7CpiWHMUvH5t21atprtc9k4ezK7+K33x4hgduSGRUXHi/fB2luqNn9MpnFFbX88DLu1i1p4Avz0vjT1+e3W8l/5F/uGMMYUH+/GhDTr9+HaWuRYte+YTDhTXc/6udFFbX8+oTmXznzrEEOPr/n//gsEBW3DqS93Ir2ZWna+8oe2jRK6+3+XgZS1/ZTUiggzXPzWXBuLgB/fpPzkklISqEH2zIweVyv/tWlPfToldebeOxUr7y+kFGD4vgrefmMnLooAHPEBzg4Ft3jOZ4SS0bj5cN+NdXSoteea1Nx8tY8YdDTE6M5PVnZhIzKMi2LPdMHk7SkBD+Z+c52zIo36VFr7zS3jMXWPGHg0xIiOS3T89gUJC9F5g5/IQnZ6ey71w1x0su2ZpF+R4teuV1ii7W89zrB0kaEspvn55BuJvMJvlQZhIhAQ5+u+uc3VGUj9GiV16lvrmV5b87QLPTxW+eyOz1VAb9ITIkgM9NS+Cvh0uovtJsdxzlQ7TolVf53toT5JTV8otHp5IeO/BvvHbnqTmpNLe6+OO+ArujKB+iRa+8xvu5Ffwpq5Avz0vn1tFD7Y7TpYy4cGanRfNmViHuOEW48k5a9Mor1Da28J01Rxk5dBBfX5Bhd5xrum/qcM5dqOdosb4pqwaGFr3yCj/akEN5bSM/e3ASwQHW55C3w6Lx8QQ6/Hj7cIndUZSP0KJXHi+7qIY/7ivkizeOYGryYLvjdCsyNIBbRseyLrsEp94pqwaApaIXkUUikisieSLyQhfbl4hItogcFpEsEbnR6rFKXQ9jDD9Yn0N0WCBfne/eQzYd3TtlOOW1Tew9e8HuKMoHdFv0IuIAXgIWA+OAR0VkXKfdtgGTjTFTgKeBV3twrFK9ti2ngr1nq/naggy3uV7eivlj4ggLdLBWh2/UALByRj8DyDPGnDHGNAOrgSUddzDGXDafXEIQBhirxyrVW61OFz96N4e0mDAenZFsd5weCQl0cMf4YWw4Wkpzq65ApfqXlaJPAAo7PC5qf+5TROR+ETkJrKftrN7yse3HL28f9smqrKy0kl35uDUHi8mvvMK3F48ZkCmH+9pdk+KpbWzV4RvV76z8dEgXz33mHSRjzFvGmDHAfcD3e3Js+/GvGGMyjTGZsbGxFmIpX+Z0GV7+IJ8JCREsHOBph/vK3JExBAf4sfVEud1RlJezUvRFQFKHx4nAVQcWjTHbgXQRienpsUpZteFoKWerrvD8LSMR6ep8wv0FBzi4cWQMW3Mq9OYp1a+sFP1+IENERohIILAUWNtxBxEZKe0/bSIyDQgELlg5VqmeMsbw0nt5pMeGccf4YXbHuS7zx8ZRXNNAbnmd3VGUF+t27lZjTKuIrAA2AQ7gNWPMcRF5tn37SuAB4AkRaQEagEfa35zt8th++l6Uj3gvt4KTZXX8/KHJ+Pl55tn8R+aPaZuqYVtOBWOGRdicRnkrccdfGTMzM01WVpbdMZSbevjXuym+2MD737rFI9+E7ezeX+7A4Se89dxcu6MoDyYiB4wxmV1t8/yfEuVTcsvq2He2midmp3hFyUPbNfWHC2uorGuyO4ryUt7xk6J8xut7zxPo78dDmUnd7+wh5o8dijHw3skKu6MoL6VFrzzGlaZW1hws5u6J8QwJC7Q7Tp8ZPzyCoeFBfHBa7x9R/UOLXnmMtw+XcLmplcdnpdgdpU+JCDdmxLAzr0onOVP9QoteeQRjDKv2nGdsfATTkqPsjtPn5mXEUlPfoguHq36hRa88wqHCGk6U1rJsVrLH3iB1LXNHxgDw4ekqm5Mob6RFrzzCqt3nGRTkz31TupwqyePFhgcxNj6CD3WcXvUDLXrl9i5eaWbd0VLun5pAWFC39/h5rHkZMRw4f5ErTa12R1FeRoteub03DxTS3OpimZe9CdvZTRmxtDgN+85W2x1FeRkteuXWXC7D63sLmJE6hNHDwu2O068yUwcT5O/Hdh2+UX1Mi165tR15VZy/UM/jszxrYZHeCA5wMGPEEHboG7Kqj2nRK7e2as95osMCWTTBs2eptOrGkTGcrrhMRV2j3VGUF9GiV26r9FIDW3PKeXh6EkH+DrvjDIg56W2XWe7O11WnVN/Roldu64/7CjHAYx62Huz1GDc8gohgfy161ae06JVbanG6WL2vgFtGxZI0JNTuOAPG4SfMSotmlxa96kNa9MotbTlRTkVdk9dfUtmV2enRFFTXU1hdb3cU5SW06JVbWrXnPAlRIdwyeqjdUQbcx+P0Z/SsXvUNLXrldvIrL7Mr/wKPzUzG4eFLBfbGqLhBRIcF6ji96jNa9MrtvL6ngACH8Mh071lcpCdEhNnp0ezKr8Idl/pUnkeLXrmVhmYnfz5QyKIJ8cQMCrI7jm3mpMdQXtvE2aordkdRXkCLXrmVd7JLqG1sZdlM37mksitz0qMB9Oob1Se06JVbeX3PeUbFDWLGiCF2R7FVSnQowyODdZxe9QlLRS8ii0QkV0TyROSFLrY/LiLZ7X92icjkDtvOichRETksIll9GV55l+yiGo4UXeLxmSleubhIT7SN08ew+8wFXLq8oLpO3Ra9iDiAl4DFwDjgUREZ12m3s8DNxphJwPeBVzptv9UYM8UYk9kHmZWX+p9d5wgLdHD/NO9cXKSn5qRHU32lmdzyOrujKA9n5Yx+BpBnjDljjGkGVgNLOu5gjNlljLnY/nAPkNi3MZW3q6xrYt2RUh7KTCIiOMDuOG5hto7Tqz5ipegTgMIOj4van7uaLwLvdnhsgM0ickBEll/tIBFZLiJZIpJVWanzcfua1/eep9np4sk5qXZHcRvDo0JIjQ5ld75OW6yuj5Wi72qwtMtBQxG5lbai/3aHp+caY6bRNvTzvIjM6+pYY8wrxphMY0xmbGyshVjKWzS1Olm1p4BbR8cyIibM7jhuZXZ6DHvPVNPqdNkdRXkwK0VfBHS8cyURKOm8k4hMAl4FlhhjPv5d0xhT0v53BfAWbUNBSn1sfXYpVZeb+MLcEXZHcTtz0qOpa2rlWEmt3VGUB7NS9PuBDBEZISKBwFJgbccdRCQZWAN83hhzqsPzYSIS/tHHwELgWF+FV57PGMOrH55l5NBB3JQRY3cctzMrrW2cXi+zVNej26I3xrQCK4BNQA7whjHmuIg8KyLPtu/2XSAa+FWnyyjjgB0icgTYB6w3xmzs8+9Ceaztp6s4UVrL8nlpPn9JZVdiw4MYHRfOLh2nV9fB38pOxpgNwIZOz63s8PEzwDNdHHcGmNz5eaU+8vL7eQyLCOa+KXpJ5dXMTo9m9f4CmltdBPrrPY6q5/RfjbLNoYKL7DlTzTM3jdACu4Y56dE0trg4XFhjdxTlofSnS9lm5Qf5RIYE8KgPLRXYGzPTovETdPhG9ZoWvbLFybJaNp8o58nZKYQFWRpB9FmRIQGMHx6pN06pXtOiV7b46cZcwoP8+eKNaXZH8Qhz0qM5VHCRhman3VGUB9KiVwNu39lq/naygq/cMpLIUJ3uwIrZ6dG0OA1Z56vtjqI8kBa9GlDGGH78bg5xEUE8pdMdWDY9dQj+fqLDN6pXtOjVgNpyopyDBTV8fcEoQgIddsfxGGFB/kxJitKiV72iRa8GTKvTxc825ZIWG8ZDN+gEpz01Jz2ao0U1XGposTuK8jBa9GrArDlYzOmKy3xr4Wj8HfpPr6fmjozBZXQ6BNVz+tOmBkRji5MXt55iclIUiyYMszuOR5qWMphBQf5sP63TeKue0aJXA+J3u89ReqmRby8arXPa9FKAw4/Z6dFsP1WJMbq8oLJOi171u0sNLbz0Xj7zRsUyJ11nqLwe80bFUnSxgbNVV+yOojyIFr3qd7/+IJ9LDS38wx2j7Y7i8W7OaFuUZ/spHb5R1mnRq35VXtvIazvPcu/k4UxIiLQ7jsdLjg4lJTqU7ad13htlnRa96lf/se00Tpfhmwv1bL6vzMuIZXf+BZpadToEZY0Wveo3+ZWX+dP+Qh6bkUxydKjdcbzGvFGxNLQ4OXD+ot1RlIfQolf95t835xLk78eK2zLsjuJVZqdHE+AQ3s/VcXpljRa96hdHCmvYcLSMZ25KIzY8yO44XmVQkD8zR0SzLafc7ijKQ2jRqz5njOEnG08yJCyQL900wu44Xmn+2KHkV17hnF5mqSzQold97sPTVezKv8CKW0cSHqzTEPeHBWPjANiqZ/XKAi161aeMMfxsUy4JUSE8PkuXCOwvSUNCGRU3iG05FXZHUR7AUtGLyCIRyRWRPBF5oYvtj4tIdvufXSIy2eqxyrtszangaPElvrYggyB/nYa4P80fG8f+c9U6m6XqVrdFLyIO4CVgMTAOeFRExnXa7SxwszFmEvB94JUeHKu8hMtl+L9bTpESHcrnpibYHcfrLRg7lFaX4QO9S1Z1w8oZ/QwgzxhzxhjTDKwGlnTcwRizyxjz0UW9e4BEq8cq77H5RBk5pbV8bX6GTkM8AKYkDWZIWKBefaO6ZeWnMQEo7PC4qP25q/ki8G5PjxWR5SKSJSJZlZV6huJpXC7Di1tOkxYTxr2Th9sdxyc4/ITbxgzlbycraG512R1HuTErRd/VnLJdzpEqIrfSVvTf7umxxphXjDGZxpjM2NhYC7GUO9maU05ueR1/N3+kns0PoDsnDqOusZWd+Tr3jbo6Kz+RRUBSh8eJQEnnnURkEvAqsMQYc6EnxyrP9+vtZ0gcHMI9k/RsfiDNHRlDeJA/G7JL7Y6i3JiVot8PZIjICBEJBJYCazvuICLJwBrg88aYUz05Vnm+rHPVHDh/kS/dlKZn8wMsyN/B7ePi2HyinBanDt+ornX7U2mMaQVWAJuAHOANY8xxEXlWRJ5t3+27QDTwKxE5LCJZ1zq2H74PZaOVH5xhcGgAD2Xqgt92WDwxnksNLezStWTVVfhb2ckYswHY0Om5lR0+fgZ4xuqxynvkVdSxNaecry/IIDTQ0j8n1cduyohhUPvwzc2j9P0t9Vn6e7a6Lv+14xzBAX48MTvV7ig+KzjAwfyxQ9l0okyHb1SXtOhVr11qaOGvh4pZMjmBIWGBdsfxaXdOjKemXodvVNe06FWv/flAEQ0tTj4/O8XuKD7v5lGxRAT789dDxXZHUW5Ii171istlWLXnPNOSo3QtWDcQHODgrknxbDxWxpWmVrvjKDejRa96ZUdeFWerrujYvBu5f2oiDS1ONh0vszuKcjNa9KpXfrf7PNFhgSyeOMzuKKpdZspgEgeH8JYO36hOtOhVj1XUNvK3k+U8PD1JpyJ2I35+wv1TE9iZV0VFbaPdcZQb0aJXPbbmUDEuAw/doDdIuZv7pybgMvD2YZ1pRH1Ci171iDGGN7MKuSFlMGmxg+yOozpJix3ElKQo3jxQiDFdzh+ofJAWveqRQ4U15Fde0bN5N7Z0ehKnyi9zqLDG7ijKTWjRqx55M6uI4AA/7poUb3cUdRV3Tx5OaKCDP+0r7H5n5RO06JVlDc1O1h0p4c4J8YQHB9gdR13FoCB/7pk0nHeyS7is19QrtOhVD2zJKaeuqZUHddjG7S2dkUR9s5N3juibskqLXvXA2sPFDIsIZlZatN1RVDemJEUxOi6c1ft1+EZp0SuLauqb+eBUJXdPisfPr6sVIpU7EREemZ7EkcIajpdcsjuOspkWvbJk47EyWpyGe6foUoGe4oFpiQQH+LFqT4HdUZTNtOiVJWuPlJAaHcpEncDMY0SGBnDv5OG8fbiYusYWu+MoG2nRq25V1Day+8wF7p08HBEdtvEky2alUN/s1PlvfJwWverWuuxSjEGHbTzQpMQoJiVGsmrPeb1T1odp0aturcsuYWx8BCOHhtsdRfXCspkpnCq/zP5zF+2OomyiRa+uqfRSAwcLarhLpyP2WPdMHk5EsD+/33Pe7ijKJpaKXkQWiUiuiOSJyAtdbB8jIrtFpElEvtlp2zkROSoih0Ukq6+Cq4Gx6VjbIhaLJ+qUB54qJNDBAzcksvFYKZV1TXbHUTbotuhFxAG8BCwGxgGPisi4TrtVA18Ffn6VT3OrMWaKMSbzesKqgbfhWBmj48JJ15kqPdrjM1NocRreyNIbqHyRlTP6GUCeMeaMMaYZWA0s6biDMabCGLMf0Gu4vEhFXSP7z1WzaIIO23i6kUMHMTstmj/sLcDp0jdlfY2Vok8AOp4GFLU/Z5UBNovIARFZ3pNwyl6bj5djDNypwzZeYdmsFIprGvjgVIXdUdQAs1L0XV043ZNTgrnGmGm0Df08LyLzuvwiIstFJEtEsiorK3vw6VV/efdYKWmxYYyK02Ebb7BwfByx4UF6p6wPslL0RUBSh8eJgOUp8YwxJe1/VwBv0TYU1NV+rxhjMo0xmbGxsVY/veonF680s+dMNYsnDNObpLxEgMOPR6cn8V5uBYXV9XbHUQPIStHvBzJEZISIBAJLgbVWPrmIhIlI+EcfAwuBY70NqwbOtpMVOF2GO8br+Lw3WTojGQH+uE/P6n1Jt0VvjGkFVgCbgBzgDWPMcRF5VkSeBRCRYSJSBHwD+CcRKRKRCCAO2CEiR4B9wHpjzMb++mZU39lyooxhEcE6t42XGR4VwvyxcbyRVUhzq8vuOGqA+FvZyRizAdjQ6bmVHT4uo21Ip7NaYPL1BFQDr7HFyfZTVTx4Q6IO23ihZbNS2HKinI3Hy7h3sk5r4Qv0zlj1GTtOV9HQ4uT2cXF2R1H94KaRMaREh7Jqt94p6yu06NVnbDlRTniQv64k5aX8/ITHZiSz71w1uWV1dsdRA0CLXn2K02XYdrKcW8YMJdBf/3l4q4cykwj09+P1vXpW7wv0J1l9yqGCi1RdbtZhGy83JCyQuybGs+ZgMVeaWu2Oo/qZFr36lC0nyglwCLeM1nsZvN2yWSlcbmrl7cOWb4tRHkqLXn3MGMPmE+XMSosmIjjA7jiqn01LjmJsfIQuSuIDtOjVx/IrL3O26goLddjGJ4gIy2Ylc6K0lkOFNXbHUf1Ii159bPOJcgAWaNH7jPumJDAoyJ9VuiiJV9OiVx/bcqKciQmRxEeG2B1FDZCwIH/un5rAuuxSLl5ptjuO6ida9AqAitpGDhXU6LCND1o2K4XmVhd/PlBkdxTVT7ToFQBbc9rmKL99vBa9rxk9LJzpqYNZtfc8Ll2UxCtp0SugbRKzpCEhjI4LtzuKssGyWSmcv1DPjrwqu6OofqBFr7jc1MrOvAvcPlbnnvdViyYMIzosUN+U9VJa9Ir3cytodrq4Q4dtfFaQv4OHpyexNaec0ksNdsdRfUyLXrH5eDlDwgLJTB1idxRlo8dmJGOAP+4r7HZf5Vm06H1cc6uL905WsGDsUBx+Omzjy5KGhHLLqFhW7yugxamLkngTLXoft/vMBeqaWnXJQAW0vSlbUdfE1vab55R30KL3cZuOlxEa6GDuyBi7oyg3cMvooSREhbBKpy/2Klr0PszlMmw5Uc4to2MJDnDYHUe5AYef8NjMZHbmXeBM5WW746g+okXvww4X1VBZ18TCcTpsoz7xcGYS/n7CH/cV2B1F9REteh+26XgZ/n7CrWOG2h1FuZHY8CBuHxfHmoPFNLfqm7LeQIveRxlj2Hy8nNnp0USG6Nzz6tMenp7EhSvNbMvRN2W9gaWiF5FFIpIrInki8kIX28eIyG4RaRKRb/bkWGWPvIr2uef1ahvVhXkZscRHBrN6v15T7w26LXoRcQAvAYuBccCjIjKu027VwFeBn/fiWGWDj+aev32s3g2rPsvhJzx0QyLbT1dSUqN3yno6K2f0M4A8Y8wZY0wzsBpY0nEHY0yFMWY/0NLTY5U9Nh0vY0pSFMMig+2OotzUQ5lJALyZpdMXezorRZ8AdPz9raj9OSssHysiy0UkS0SyKisrLX561RslNQ1kF11ioc5to64haUgoc9NjeCOrUKcv9nBWir6r++Kt/le3fKwx5hVjTKYxJjM2Ntbip1e9saV92EbvhlXdeWR6EsU1DezM1+mLPZmVoi8Ckjo8TgRKLH7+6zlW9ZONx8pIjw0jPXaQ3VGUm1s4Po6o0AD+pG/KejQrRb8fyBCRESISCCwF1lr8/NdzrOoHFXWN7D17gbsmDbc7ivIAQf4O7p+awObj5bqmrAfrtuiNMa3ACmATkAO8YYw5LiLPisizACIyTESKgG8A/yQiRSIScbVj++ubUd1792gZLgP3TIq3O4ryEI9MT6LZ6eKtQ8V2R1G95G9lJ2PMBmBDp+dWdvi4jLZhGUvHKvusyy5hdFw4GbpkoLJozLAIJidF8af9hXxhbqquQuaB9M5YH1J6qYH95y5yt57Nqx56JDOJ3PI6jhRdsjuK6gUteh+yPrsUgLu06FUP3TM5npAAB3/arxOdeSIteh+yLruUcfERpOnVNqqHwoMDuGtSPGsPl3ClqdXuOKqHtOh9RMGFeg4X1nD3ZD2bV73zyPQkrjQ7WX+01O4oqoe06H3EmkNFiMB9U6ze1KzUp2WmDCYtNow39Jp6j6NF7wOMMaw5WMyc9GiGR4XYHUd5KBHhkcwkss5fJK+izu44qge06H3A/nMXKaiu54FpXV4Bq5Rln5uWiL+fsHqfntV7Ei16H/CXA0WEBjp0bht13T5afeovB4tobHHaHUdZpEXv5Rpb2t48WzwhnrAgS/fHKXVNy2alcLG+hQ36pqzH0KL3cpuOl3G5qZUHbtA3YVXfmJMeTVpMGKv2nLc7irJIi97Lvb63gKQhIcwaEW13FOUlRITHZiZzsKCGEyW1dsdRFmjRe7FT5XXsO1vNspkp+Pnp/CSq7zx4QyJB/n6s2qtn9Z5Ai96LrdpznkB/v4+XhFOqr0SFBnLP5OH89VAxtY2dVxBV7kaL3ktdaWplzcFi7p4Yz5CwQLvjKC/05OxU6pudegOVB9Ci91J/PVzM5aZWHp+VYncU5aUmJkYyPXUwv919DqeuKevWtOi9kDGG3+8+z7j4CKYlR9kdR3mxp+eOoLC6ga055XZHUdegRe+F3s+t5GRZHU/fOEIXiVD96vZxcSREhfDajrN2R1HXoEXvhV5+P5/hkcEsmaLrwqr+5e/w48k5Kew9W83xEl2UxF1p0XuZrHPV7DtXzZfmpRHg0P+8qv89Mj2ZsEAHr2w/Y3cUdRXaBF5m5Qf5DAkLZOn0ZLujKB8RGRLA47NSeOdICQUX6u2Oo7qgRe9Fckpr2ZpTwVNzUgkJdNgdR/mQZ24cgb+fHyu359sdRXXBUtGLyCIRyRWRPBF5oYvtIiK/aN+eLSLTOmw7JyJHReSwiGT1ZXj1aT/ZeJLIkACenJ1qdxTlY4ZGBPNgZiJ/ziqivLbR7jiqk26LXkQcwEvAYmAc8KiIjOu022Igo/3PcuDlTttvNcZMMcZkXn9k1ZVdeVW8n1vJ87emExkaYHcc5YO+PC+NVpeLVz/UsXp3Y+WMfgaQZ4w5Y4xpBlYDSzrtswT4nWmzB4gSEV2cdIC4XIYfvXuShKgQntCzeWWTlOgw7p08nFV7Cqio07N6d2Kl6BOAjvc4F7U/Z3UfA2wWkQMisvxqX0RElotIls/vUtcAAAuYSURBVIhkVVZWWoilPvJOdglHiy/x9wtHERygY/PKPl9fMIoWp4tf/i3P7iiqAytF39UdN53vd77WPnONMdNoG955XkTmdfVFjDGvGGMyjTGZsbGxFmIpgNrGFn64IYdx8RG68LeyXWpMGA9PT+IPewv0Chw3YqXoi4CO0x8mAiVW9zHGfPR3BfAWbUNBqo/8bGMulXVN/OhzE3UqYuUWvjY/A4ef8OLWU3ZHUe2sFP1+IENERohIILAUWNtpn7XAE+1X38wCLhljSkUkTETCAUQkDFgIHOvD/D7twPlqVu09z1NzRjA5See0Ue4hLiKYp+am8tfDxRwr1rtl3UG3RW+MaQVWAJuAHOANY8xxEXlWRJ5t320DcAbIA34DPNf+fBywQ0SOAPuA9caYjX38PfikxhYnL/zlKMMjQ/j7haPsjqPUpzx380gGhwbyf9Yex6UzW9rO0mrRxpgNtJV5x+dWdvjYAM93cdwZYPJ1ZlRd+Lf1JzhdcZnfPj1DF/1WbicyNIAXFo/hH/6czZ8PFvGwLn5jK70z1gOtyy5h1Z4CvjwvjZtH6RvXyj09OC2RG1IG8+N3T3KpXlehspMWvYc5V3WF7/zlKFOTo/jmHaPtjqPUVfn5Cf+6ZDw19c386N0cu+P4NC16D1JT38zTv92PwyH856NTdXZK5fbGD49k+bx0Vu8vZPPxMrvj+CxtCg/R1Opk+e8PUFTdwCufzyRxcKjdkZSy5Bu3j2JcfAQvrDmqd8zaRIveAzhdhr9/4wj7zlbz84cnM2PEELsjKWVZoL8fv3h0CleaWvnmm9l6FY4NtOjdXKvTxTfeOMy67FL+8c4x3DtZV41Snmfk0HD++e5xbD9VyU835dodx+fodXlurMXp4htvHOGdIyV8e9EYls9LtzuSUr32+MxkckprWflBPmmxYXrJ5QDSondTdY0tPPf6QT48XcULi8fw7M1a8sqziQj/cu94zl+o5x/XHGVYRDDz9PLgAaFDN26opKaBh1buZnf+BX7ywEQteeU1Ahx+vPT4NDLiwnnmd1lsyym3O5JP0KJ3M9tPVXL3f+6g+GID//OFGTyia78qLxMZEsAfvzSTMcPC+fLvD7Auu/MciaqvadG7iRani3/fnMuT/72P2EFBvPX8XG7MiLE7llL9Iio0kFXPzGRKUhQr/nCIH797klany+5YXkvH6N3AybJavvnmEY4V1/LgDYl8f8kEXdxbeb2I4ABWPTOTf113gpUf5HO48CI/e3AySUP0HpG+Jm3zkbmXzMxMk5Xl/euINzQ7+dX7eaz8IJ+I4AB+cP8EFk3QFRiV7/nzgSK++/YxXMbwd7dl8KWb0gj01wGHnhCRA1dbl1uL3gbGGDYcLeOHG3IormlgyZThfPfucUQPCrI7mlK2Kalp4F/fOcHG42UkRIXwlVvSeSgzkSB//e3WCi16N2GMabuSZuNJjhRdYsywcL5373hmpkXbHU0pt7H9VCUvbj3FoYIahoYHsXR6Eo/MSCYhKsTuaG5Ni95mLpdha045Kz/I52BBDcMjg/n67aN4YFoiDl3+T6nPMMawM+8C/7XjDO+fqkSAmSOiuXtyPIvGD9PffrugRW+T5lYXbx8u5tfbz5BXcZnEwSF86aY0HpmeRHCA/jqqlBWF1fW8eaCIddklnKm8gsNPmJ0WzV2T4lkwNo7YcC190KIfcKWXGli9r5DV+wsor21izLBwvnJLOndNjMdfpxZWqleMMZwsq2N9dinrsks4d6EeEZicGMXt4+KYP3Yoo+PCEfHN35K16AeAy2XYfrqS1/cWsC2nHAPMy4jlqbmp3DIq1mf/8SnVHz4q/W055WzJqeBIYQ0AiYNDWDA2jgVj45gxYohPXbmjRd+Pcsvq+OvhYtYeLqG4poHosEAenp7Eo9OTSY7W64GVGggVtY1sO1nBtpxyPjxdRVOri/Agf+aNjuW20UO5MSOGuIhgu2P2Ky36PlZYXc+Go6W8daiYk2V1OPyEmzJieGBaIneMH+ZTZxFKuZuGZic78qrYeqKcbScrqLrcBEDG0EHcmBHDjSNjyEwdQmRIgM1J+5YW/XVqanWy/+xF3sut4P3cCvIrrwAwNTmK+6YkcNekeGL0KgCl3I7LZThRWsvOvCp25FWx72w1Ta1tUy2kxYYxJSmKqUlRTEkaTEbcII++SOK6i15EFgH/ATiAV40xP+60Xdq33wnUA08ZYw5aObYrdha9y2U4X13P0eJLHC2qIbvoEkeLL1Hf7CTQ349ZadHcOjqW+WPidGhGKQ/T2OLkYMFFDhXUcKighsOFNR+f8Yu0jfGPjB3EyKGDSI8dRMLgEIZHhTA8MsTtpyW5VtF3O9eNiDiAl4DbgSJgv4isNcac6LDbYiCj/c9M4GVgpsVj+5zTZWh1uXC6DC1O0/bY6aKhxUldYyt1ja1cbmqlrrGF8tomimvqKalppKSmgaKLDVxuagXalkAbFx/BQzckcvPoWGalRRMaqNMDKeWpggMczEmPYU5624SBxhiKaxrILrrEqfI68iouk1dxmZ35F2hu/fQka0PCAomLCGZwaACDQwOJ6vB3VGggUSEBhAQ6CPL3IzjAQXCAH0H+DoIDHAQF+OEQweEniICfSPsfBuRCDSutNQPIM8acARCR1cASoGNZLwF+Z9p+PdgjIlEiEg+kWji2z0z8l01cbmqlp6NRkSEBJESFkDg4lJkjhjA2PoKJiZGMigsnQC+HVMpriQiJg0NJHBzKnRM/mWfK6TKU1DS0/bnU8PGJYHltEzX1zZwsq6WmvoWahhacfbAGrl97+Q8ND2LXd+Zf9+frzErRJwCFHR4X0XbW3t0+CRaPBUBElgPL2x9eFpEBXVgyeyC/WM/FAFV2h/AQ+lpZo6+TNQP6OuUD8o+9PjzlahusFH1Xv1d0/l/Y1faxcmzbk8a8ArxiIY/PEZGsq429qU/T18oafZ2s8ZbXyUrRFwEdV/FNBDovCXO1fQItHKuUUqofWRmA3g9kiMgIEQkElgJrO+2zFnhC2swCLhljSi0eq5RSqh91e0ZvjGkVkRXAJtoukXzNGHNcRJ5t374S2EDbpZV5tF1e+YVrHdsv34l30yEt6/S1skZfJ2u84nVyyxumlFJK9R29dlAppbycFr1SSnk5LXo3JiIPichxEXGJyGcu8RKRZBG5LCLftCOfu7ja6yQit4vIARE52v73bXbmtNu1/j2JyHdEJE9EckXkDrsyuiMRmSIie0TksIhkicgMuzP1lBa9ezsGfA7YfpXtLwLvDlwct3W116kKuMcYMxF4Evj9QAdzM12+TiIyjrYr4sYDi4BftU9fotr8FPieMWYK8N32xx5FJ25xY8aYHOh6LgwRuQ84A1wZ4Fhu52qvkzHmUIeHx4FgEQkyxjQNYDy3cY1/T0uA1e2vy1kRyaNt6pPdA5vQbRkgov3jSDzwXiAteg8kImHAt2mbLM6nh2164AHgkK+WfDcSgD0dHn80hYlq83Vgk4j8nLZRkDk25+kxLXqbichWYFgXm/63Mebtqxz2PeBFY8xlX1misJev00fHjgd+Aizsj2zupJevk+WpSrzVtV43YD7wv4wxfxGRh4H/AhYMZL7rpUVvM2NMb/7BzAQeFJGfAlGAS0QajTG/7Nt07qOXrxMikgi8BTxhjMnv21Tup5evk5VpTrzatV43Efkd8LX2h28Crw5IqD6kb8Z6IGPMTcaYVGNMKvD/gB96c8n3lohEAeuB7xhjdtqdx42tBZaKSJCIjKBtXYl9NmdyJyXAze0f3wactjFLr2jRuzERuV9EioDZwHoR2WR3Jnd0jddpBTAS+Of2S+MOi8hQ24La7GqvU/u0JG/Qtk7ERuB5Y4zTvqRu50vAv4vIEeCHfDKdusfQKRCUUsrL6Rm9Ukp5OS16pZTyclr0Sinl5bTolVLKy2nRK6WUl9OiV0opL6dFr5RSXu7/A2YLCY3ZJfJAAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#it's definitely got lost:\n",
"kdeplot(trace['center'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Comparison to other measures:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(-5.771172458391055, 0.34871679186296234)\n"
]
}
],
"source": [
"#Now we can compare to another technique, based on the autocorrelation function. This is from \n",
"#Use of the Weighted Histogram Analysis Method for the Analysis of Simulated and Parallel Tempering Simulations\n",
"#First it calculates the autocorrelation time constant, then calculates SEM based on 'effective number of samples'\n",
"#then just apply the 95 CI to that SEM. \n",
"z95 = 1.959963984540054\n",
"\n",
"#this is the Chodera lab 'statistical inefficiency' measure.\n",
"def statistical_inefficiency(corr, mintime=3):\n",
" N = corr.size\n",
" C_t = sm.tsa.stattools.acf(corr, fft=True, unbiased=True, nlags=N)\n",
" t_grid = np.arange(N).astype('float')\n",
" g_t = 2.0 * C_t * (1.0 - t_grid / float(N))\n",
" ind = np.where((C_t <= 0) & (t_grid > mintime))[0][0]\n",
" g = 1.0 + g_t[1:ind].sum()\n",
" return max(1.0, g)\n",
"\n",
" #and these actually return the values for SEM:\n",
"def ci_from_chodera(timeseries):\n",
" autocorrelation_time = statistical_inefficiency(timeseries)\n",
" #original n\n",
" n = len(timeseries)\n",
" #calculate SEM and return CI using effective n:\n",
" sem = np.std(timeseries) / np.sqrt(n/autocorrelation_time)\n",
" #return CI:\n",
" return timeseries.mean()-sem*z95, timeseries.mean()+sem*z95\n",
"\n",
"print(sem_from_chodera(timeseries))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment