Changepoint Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Standard Data Science Helpers\n",
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "import pandas as pd\n",
    "import featuretools as ft\n",
    "\n",
    "# Graphic libraries\n",
    "import matplotlib as plt\n",
    "%matplotlib inline\n",
    "plt.style.use('fivethirtyeight');plt.rcParams['font.size']=18\n",
    "import seaborn as sns\n",
    "\n",
    "# Extra options\n",
    "pd.options.display.max_rows = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Closing Value</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1986-01-02</th>\n",
       "      <td>25.56</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1986-01-03</th>\n",
       "      <td>26.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1986-01-06</th>\n",
       "      <td>26.53</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1986-01-07</th>\n",
       "      <td>25.85</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1986-01-08</th>\n",
       "      <td>25.87</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            Closing Value\n",
       "Date                     \n",
       "1986-01-02          25.56\n",
       "1986-01-03          26.00\n",
       "1986-01-06          26.53\n",
       "1986-01-07          25.85\n",
       "1986-01-08          25.87"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv('data/crude_oil_prices.csv',\n",
    " parse_dates = ['Date']).set_index('Date')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAD8CAYAAABDwhLXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8lfXZ+PHPlYQskhBC2BuNggNxodaBC7difdwLqa1tta3VWm1rrX1a7dP++jxV21q12iJ1j6qodS8UB2KRioAaBGTPkED2un5/3Hfg5Jz77J1c79eLFzn3WVdOkvu6v+v6iqpijDHGJFtOugMwxhjTO1jCMcYYkxKWcIwxxqSEJRxjjDEpYQnHGGNMSuSl403r6upsapwxxvRw/fr1E9/b1sIxxhiTEpZwjDHGpIQlHD/V1dXpDiEiFmdiWZyJZXEmVk+J0xKOMcaYlLCEY4wxJiUs4RhjjEkJSzjGGGNSIi3rcIwxvcfSbW0c9swmAEr7CCsuHEpejoR5lumJrIVjjEmqM17asvPrHW3KrC8a0hiNSSdLOMaYpNrc3Nnt9o/er0tTJCbdLOEYY4xJCUs4xhhjUiJswhGRv4vIJhH51OO+60RERaTSvS0i8kcRWSYin4jIAckI2hhjTPaJpIVzP3CS/0ERGQlMBVb5HD4ZqHL/XQHcFX+IxhhjeoKwCUdV3wZqPO66Dbge8N1qYBrwD3V8AJSLyNCERGqMMSarxTSGIyJnAGtV9T9+dw0HVvvcXuMeM8YY08tFvfBTRIqBG4ETvO72OBZys7VMrIKaiTF5sTgTy+JMrOrqauraAIo978sUmRRLKNkSZyixVBrYDRgL/EdEAEYAC0RkMk6LZqTPY0cA60K9WFVVVQwhJE91dXXGxeTF4kwsizOxuuK8Y9EOYHvA/ZnyPWTb55npEr49gaouUtVBqjpGVcfgJJkDVHUD8CxwqTtb7VCgTlXXxxC3MaYHuPmjwGRjeq9IpkU/ArwP7Ckia0Tk8hAPfwFYDiwD7gWuTEiUxhhjsl7YLjVVvSDM/WN8vlbgqvjDMsYY09NYpQFjjDEpYQnHGJNSZ40tSncIJk0s4RhjUqpvnu2F01tZwjHGpFTIhXmmR7OEY4xJipYO79QS7Ljp+SzhGGOSQoPklSeWN6U2EJMxLOEYY5Ki0zrPjB9LOMaYpOi0fGP8WMIxxiSFJRzjzxKOMSYpLOEYf5ZwjDFJYfnG+LOEY4xJio5g09RMr2UJxxiTFCu2d6Q7BJNhLOEYY5Liua+Cr7dRa/30SpZwjDFJ8ea6lqD3NbZbwumNLOEYY5KiIDf4fdtaOlMXiMkYlnCMMUlxxujg2xCc8dIWapptjKe3iWSL6b+LyCYR+dTn2O9F5DMR+UREnhaRcp/7fioiy0TkcxE5MVmBG2MyW5+c4NsQLN/RwX2fNaQwGpMJImnh3A+c5HfsVWAfVZ0IfAH8FEBE9gLOB/Z2n/MXEQnRsDbG9FTtYSYG/ObjHSmKxGSKsAlHVd8GavyOvaKq7e7ND4AR7tfTgEdVtUVVVwDLgMkJjNcYkyU6bJjG+EnEGM43gBfdr4cDq33uW+MeM8b0Mi+vaU53CCbD5MXzZBG5EWgHHuo65PGwkO3q6urqeEJIikyMyYvFmVgWZ2K9v7E17GMy4XvJhBgikS1xhhJzwhGR6cBpwHG6axXXGmCkz8NGAOtCvU5VVVWsISRFdXV1xsXkxeJMLIszsSI9Oab7e8mmzzNb4gwlpi41ETkJuAE4Q1Ubfe56FjhfRApEZCxQBXwYy3sYY4zpWcK2cETkEeBooFJE1gA348xKKwBeFRGAD1T1O6q6WEQeB5bgdLVdpao22d4YY0z4hKOqF3gc/luIx98K3BpPUMYYY3oeqzRgjDEmJSzhGGMSzopBGy+WcIwxCdeWAQnnn8sb2f/JDRz//CY+q21LdzgGSzjGmCRoTXOVgYa2Ti6fs40VOzr4aHMbN31Yl96ADGAJxxiTBOlu4dzol2BeXRt8bx6TOpZwjDEJ19YZvFJ0Ktz/RWP4ByVAY3sn8za2sLnJVn9EIq7SNsYY42VZY3oTTirsaOtk5IPrd95+d9og9q7ok8aIMp+1cIwxCXf14sJ0hxCRN7bk8v2525i9sinq5579ytZuty9+Y2uQR5oulnCMMQnV3pkBU9Qi8MHGFm74rIAHqhuZ/mYNr6yOrrr1vE3di5Ou2GHdauFYwjHGJFRzR3YknCve3tbt9rmvRd5C2RLh9tj/2drK4c9s5OLXt7Kx0RKSJRxjTEK1eCSc/Sszb2xjVX1gAtgaYSK58LWasI/p6FSmPLuZxdvaeX5VM8c8tynqGHsaSzjGmITyOmd/vCU7Fl5e8kb4RALw4ebwe/089mX3mXLrGm0LVEs4xpiE8mrhZIv3NrayPs6ur9YOZc66ZmZ+3hBw37K67Ei8yWIJxxiTUOkew5mzLr5Fnnd+Wh/y/tqW4C2VTlVOfXEz017eyvzNgcnloKd6d7eaJRxjTEK1pjnhXPZWZIP/xXnea4VeXRPdbLXuz23xTDTGYQnHGJNQ6W7hbGsJ//5tnUpju/fjCoMkoi6hXv2ORTvCvne2TBtPBks4xpiECnIeT7sN7tjMG2ubGThrXdDHbQ9TebQzxN4L720MP5mgvRfPHbCEY4xJqGV17ekOwdP4xzbw1Y52znoldJdbuAWcqz2mU0cjQ/NxSoRNOCLydxHZJCKf+hyrEJFXRaTa/b+/e1xE5I8iskxEPhGRA5IZvDEm87wU5Yr9VPrtwvBdXlOHF4S8P94ew1AtpJ4ukhbO/cBJfsd+AryuqlXA6+5tgJOBKvffFcBdiQnTGJMtMjnhLN0WfkB/vwH5Ie9vi3MMphf3qIVPOKr6NuC/GmoaMMv9ehZwps/xf6jjA6BcRIYmKlhjTOZSVV4OkmwOHRT6JJ4qC7eGTzjhEkpbnBmjF88ZiHkMZ7Cqrgdw/x/kHh8OrPZ53Br3mDGmB1NVLnmjhvOC1CP738PKUxxR7La2dHLF2zUc+vRG7lwcuCYn3llmvbhHLeH74XjNJwz58VZXVyc4hPhlYkxeLM7Esjhj9+mOHJ5fFXxLgoKtK4HigOPJ+V4C3ycaD1bvKklz44d17Na2gUU7cljWkMPpg9t5c2seEHttuOovv6Q8hqdn4s89WrEmnI0iMlRV17tdZl3LZ9cAI30eNwIIPv8QqKqqijGE5Kiurs64mLxYnIllccbnwfl1QPAV+lVVVTB3rffxRPN4n3ic/3HRzq8fWx9/EdJrq/vxzrRB4R/oI1N/7v7CJcVYu9SeBaa7X08HZvscv9SdrXYoUNfV9WaM6bnyQpxJ7jqyf+oCyQKLanpvJYJIpkU/ArwP7Ckia0TkcuC3wFQRqQamurcBXgCWA8uAe4ErkxK1MSajzPy8Meh9w4pzUxhJZrnziOwZu0qFsF1qqnpBkLuO83isAlfFG5QxJrvUhChomd978w0XVfXlqrm1AcdfXNXEyaOKPJ7Rs1mlAWNMUvXJCV2bLNGCFeXMJBe8Htm+Oz2NJRxjTFL1SfFZJlOKY+5e5nQgjSzpxU08P5ZwjDFJlSepbXFkSvHQw4c4i10zZdFrJrCEY4xJquoUFvNcVd8e9Ur+83ZLzlhK3z5Oor1y75KkvH42soRjjEmqiQPiX7sSiYVbWpn4xMaon3f1vqVxve+3xvf1PP7hJmergkkp+v6zgSUcY0xSVRSk5jRz9HObo37OK4c00jeOSQYlecLvDu3ned9H7s6fkuIuxUxmCccYk1Rdi0Kv2Tfzupb694FRcQzq331Uf3JE+OPhgettnj2pcufXZX0s6YAlHGNMnMLt79I1aeCyPbt3PSWy5RPLHjNFuU5c8bRAThvtjP9cukdfnpw6YOfx/gXCUUN37avz0qkDY36PniTRxTuNMb3Mi6tC73/T1cIp9bvK70hg2eS31rVE/ZymeHdS83P8iEJqZ3gXxx/t0YrqVCWnl3W3WQvHGBOXj8PsMeM2JAJOrolcLhNq2+g/eXR3pVpfj8VI8e6rk40s4Rhj4hMmcXR1WeX6nW1iaWCs3NHO7xZu5+kVjajbQuoIk7laI8hsqVgr418BId6dQ7ORdakZY+ITYa+Q/wLQaLvUdrR1MunJXdOeZ7CNtRcPpTHMSs/57vTkUPqmYFDfv6J2u7VwjDEmOv1DDP6fOmrXpmy5fuf0jihPuH//rCHg2Murm6lvC51wwt0PsZ0Ivxlk/U0wfcRaOJZwjDFxGVIUeBq5anQrtxxcxt+mVOw8FpBwojzfvrM+cGLArxZspyFMC+ewIQUh7weYvzl8Kwhgj35Op9CoklyumRjdglH/YZzWXtjCsS41Y0xcZq9s6nb7nHFFXDaskaqq7idk/6LRCqhqxNOSt3lsgaAKTWESTiSvXtsaWfZ764yBLN/ewejSXEqjrEraJ9daOJZwjDEx61Tl2a+6T4vO92/KuESEHOk+O61DIdxC/znrWpj28hbP+76q72Dqv0JXGEjU9OtZx1RQnJfDPhWxdQxtaOzodruuFzZxrEvNGBOzOo+WQZB843lfJN1qwZJNpGJtSPiWvLnryP5MGxNfkU//adAPVgffJbWniivhiMg1IrJYRD4VkUdEpFBExorIPBGpFpHHRMRqcxvTQ3nllrc9xlq6BCac4NlAVVm4JbKxlWC+s1fwgf2iUJkRpzrA1unD2Dp9GBfsXhxXHF7uXRo4CaKniznhiMhw4AfAQaq6D5ALnA/8DrhNVauAbcDliQjUGJN52j0SxsodHR6PdOQGTI0O/to/fK82poKcXU4fXciPJpYGPcn5jinNPrEy4P48gdwcITfFO5b2ZPF2qeUBRSKSBxQD64FjgSfd+2cBZ8b5HsaYDBVuwN5fpFOjNzV1MOuL2LucRvTN5YFjBzCwKDforAHfZDm6NLD0TJgGkIlBzAlHVdcC/wuswkk0dcC/gVpV7dpxaQ3gXVzIGJP1lm8P3prx4t9YCFZ0c3V9dK/r78YDynze0ztz+I6peDViwi0ojVa+jZjHPktNRPoD04CxQC3wBHCyx0ND/tSqq6tjDSFpMjEmLxZnYlmc0fvGvCKCNSE84+zs/vjqL5dT4THKu2BrLhB+/UwwB3aupevta7fkAYFvUpqrO+Nc0yxA90kBn65YTd9tiZtJ9v0xefzf8u5xRPOzzKSfe6zimRZ9PLBCVTcDiMhTwNeAchHJc1s5I4B1oV6kqqoqjhASr7q6OuNi8mJxJpbFGb2FW1rZ2hY4xvKNPfsCjZ5x5v97fbeaLqPHjmNIcWB31jsdDTjXsdH7/LwhDPZ5zatGdXL7yvUBs8T27F8ANFFVVUX9llag+/ey26gRVA2OPen5OyC/CZbXdDsW6c8y2T/3TlWeWN5EXUsnF1UVexYbjUS4pBhPI28VcKiIFIuzcus4YAnwJnC2+5jpwOw43sMYk0Kqyo0f1lE+cy3lM9fyp0U7gj728jk1nsevnRh8o7VIp0Xv3T+2a+GyfOmWbADKC3L4yxH9Qz6v0GPApiDBkwU2N2fuuptffrSdb7+9jevn1XH2q8Erb8crnjGceTiTAxYAi9zX+itwA3CtiCwDBgB/S0CcxpgUeGdDK3curt95+6aPtrO12Xs85UuP8ZspQwsYURI8WQTOUvPOOLGuidwepGLAObsFTmtWn97+8eWBMe9d0Se2IIIYWBh4uq3PkD0K/vjprp/5+xtbWbG9PcSjYxfXMJaq3qyq41V1H1W9RFVbVHW5qk5W1d1V9RxVjX5nJGNMWtz4YV3AsVfXRP4nvHu/0C2TwEkD3o+LZEuBePm+g/htE/3rg8soSPA0teNHFAYcezhDF3/WeJQRSgQrbWOMoblduWXBdhbVBG6mFk13flmYMv+RTotuTfBunF78G1eX7tGXS/eIrgJ0NLwS2K0fb+eKvYJ3QaZLspYe2UQ9Ywx/+7yBP/t0pfnqE8XZpyRMdvKf7tzm0aXW3qncvSTxq/B/uG/3E/uPJ0VX7TkZvEoDZYJkLUGyhGOM8exK6+LVs6RBxl7CrTXxX9ryxtru3XVf1LYx/c0a5oQojxOr7+xVwv6VzrjMWWOLOH54YBdXqh01NHGz4BJpdUN866CCsS41Y0xIy3cEDiBvD7Kp2Rd10Q0237Okniv3dloe9y2t58cf1IXbsTpmQ4pzeeO0gXSq02UU6bYIyXTaqPQnvc1NgcnlniX1nD46vmKlXqyFY4wJ6ab52wOOvbK62eORsE+UM7t8z/nXJSDZXDEh9BiMiFMbLV3J5ut+FaczoUPtTI9q3L89pNzjkfGzhGOMiVptkHnL48ujTDiJCMbHtydk3gC8r0F+u6OmYG5ESA1tnSzeFtgq3b0sOZ1flnCM6SVUlYerG/jxB7V87FP2P9h4TCiDigKrA0D0m52tCFFZOha7hZmWnW6R1pJLlZs/Cmy9AuQmKTNk9k/HGJMwj37ZxJVznXIx93/ewKfnOCVgWmI452/NwFXz/zO5X7pDCMt/4Wua8w33feY9GzBZlbIt4RjTS/zcZyZaWyf8duF2Lqnqy0VvRF/K5Nr3veucpaOk/8VVxRw/vJBpY9I/AB9OQAsnPWGEFazCdtyvm5RXNcZknK1+q8efWdnEsc9vZn1j+NNepF0/+w2IbYPftRFOw33o2IqAYzdMKuXMsUUZMessHP+Ek+4xnFSzhGNML7WtJfKzXSR7w/z3QWWUF8R2Srl7ifeiU183TCrlWI+1M5mfZnbxz9ufbQus7NCTWcIxxoTV7F6Kb27qoHzm2oD7TxtVyNX7xr5yP1wD6oFjK/jp/mUUeM9VyBqvrOk+nfzx5U1piWNRTRunvxj79t2xsjEcY3qB6rrorqTL8qVb5eWuraQPfXqT5+Ov2y++MjEtYfqWump5eo0tJHhjzqRaWpucKszRqG/r5MjZ3j/HZLMWjjG9QLRbNvuX0m9qV2pbOgPGgbr0DVO0M5yxYdZ9VHqU9u/SloLK0olSmAEttHRWqLaEY0wvEM3g9KGD8gM2JNvepju71bwMiHHspsvgouDPH16cy+RBwScjZMiWMhH56f5l3W7H+7nF4v2NreEflCSWcIzpBUIlC38XVRUHrD4//vnNjH9sQ9DnVMR56R4qvvuPqehWsXpMaff3GleaPSMDB1R2T5zjY9zZNB5Pr0zPuBFYwjGmV2iKYqAj0RuPRaI1RI/fwX6tmzu+Vk5FQQ59cuD2r5VTmJc989T8Jz2kYt+faM37+qCkvXb2XBoYY2IWzTjHjjT0UQVr4XhVD5gyrJBlFwwBkrdAMVny/RbidJWkU1X+sqQBVeWKCSXkp2MFLVBz2bCkfqZxJRwRKQfuA/bBKXz6DeBz4DFgDLASOFdVt8UVpTEmLlfN9a4M4CUnypUt9xzVP9pwAgSbpXaCx7bMkH2Jpov/ZnZtHcq6hg72enxXd+XP529ny/Rh5CVr280Qkv25xtuldgfwkqqOB/YDlgI/AV5X1Srgdfe2MSZNPvXYNjqUaAtwnjQy8pIy39nLe/uAYC2coizqLotEvn+XWifdkk2Xh5elfibZhbsXJ/09Yk44IlIGHAX8DUBVW1W1FpgGzHIfNgs4M94gjTGxe+CL6LZrjmZdy8kjC+kXbptPH+ft5n1SC9bCyYRpxInk36XWEqSr8wfvRt4ijUaortVfHlQW9L5EiadLbRywGZgpIvsB/wauBgar6noAVV0vIiFHoKqrq+MIITkyMSYvFmdi9dQ471ka3ZVrZdNGILJWy82jaqiurvG8zyvO/E6A4oDHvfVVARCYXVavXM7WJCedVP7cN7cIsGsTtu3N7QQrzuMfVyLiXN/c/f191a1ZTvCNxhMjnoSTBxwAfF9V54nIHcTQfVZVVRVHCIlXXV2dcTF5sTgTq0fHOTewFE0o/3XAOBZ1bueOT8PXNxu/h3csweJs61R4b93O23ninAM+DhLjxD13T2pRzlT/3CuaO2D+ri60uvbg35tvXImIs1OV/31nGxA4LfrMMUVUVQ2P6/UhfFKMZwxnDbBGVee5t5/ESUAbRWQogPt/emooGNPLvbK62bPuWTgiwnf3Ts7Omf4nnFDdd6+eOjArKkBHw3/SQDBnjE78Vgur6jt47EvvNTj3TYl/4kckYk44qroBWC0ie7qHjgOWAM8C091j04HZcUVojAmrU5W/LK7njJe2cM+Sera3dnLua9Hvc/Ot8c6g/pDi5PRjeZ1vvXYcHVmSG7D+pieIdI1TMqr1/L+FOzyPf37ekJTNiIt3Hc73gYdEJB9YDszASWKPi8jlwCrgnDjfwxgTxoPVjfzM3WDt7fUt3baQjsYtSd4106vF8pVHnbf9B/RJahzp0ifCS/xItoOIVrCZb4OTdHHhJa6Eo6oLgYM87jountc1xkTHf1bTo0G6ToK5cf9SZozv2+0K/O4j+/Odd4IvoQs2xTlam5sCF5pePyn5M6bSIdJ1LsFmr0Vi5Y52Ojpht37hT+97lad27b+VtjGml/FfqHn9pFJ+PKmMSr85yOfu5j2bqcu3JyRmnOf2RYFdPRNSfCLMNG3RFfemvVP5xls1lM9cy6QnN3LgUxuZ+nz44fMXThkYY4SxsYRjTC8zdXgByy4YwqxjKljwX4P52f7erYkcEb44f4jnfd8c3zfslgKR+teq7puS5QrkpmGVfSaJptgqwF2L63lqRfdW7fzNbZz9yhbPMbIuse7QGitLOMZkuVAnFC+FeUJlYS7TxhQxLkzSGFSUy6xjKgKOj+ibvH5//8WRvdEnUVaHuOmj7Z7HX1vbEpCIuvw9RTPTfPXudqsxPUBTlFfDRVEWhnx/Y0vAscsnJGb8xgRX29IZUQvkhVWhx+u+OWeb53TsAwamfhagJRxjslyo0v5eol3b8vb6wIRTksQaZ9Em0Gxz9T4lES2q3d4WmHDaOpUfvlfLcyubOGxIAYMKc3ggzA6eClz6ZmA1iPIoShIlinWpGZPlWpO8xfJ39wqcHNDTFmSmUqRrnLxaJa+vbeah6ka2tykvr24Om2xC6ZefhmrUKX9HY0xC7WiNPOEcWBn9+pZJlT1vAWY6RVoBe8JjG7jmvW3dNs87/zXvunWxSMdFgyUcY7Lc7xZ6Dxh7iaUtNKE8jyFFu04VxwwriOFVTJdodlSd+XkjTyxP/VYFyWIJx5gstqimjceXR77Ic1GUs58A8nKEv06p4IDKPhw5JJ/fHZLcagQ9XbSTNpKxVUEahm8AmzRgTFb76bzoTkax7h591NAC3jg9eXvd9yaN7dH/EBo7oK41cVt/v5mmn6W1cIzJYh9sjK5m2r4V6a9RdkgPLMoZjZYoZxUC/HN9Hu9tCJwtGKu90/R7YC0cY7JYtDUeDxuc/pN9cQ/bNjpasWyb/XlDDvVrEpdw0sUSjjG9REmecMOk0nSH4blFQW9SEEORhpWNOby8ObqtwoMZXJS+ji3rUjOmF/jFgWW8M20QAwpTV4o+mHCbkGVCUkymqSOi31ztkPIY+uGC8C/emkrWwjGmh3vmxAEcPSzxO0jGKtwMqW+O79llc0oi3RTHx8fbI3/O9/cpYcm2Nl5fG9gFd93EUo4ckr5p7dbCMSZLrW+M7Ko3k5INhF+HMrAo/a2wZIt2HGvRjsg/kxNGFPLE1AGe9/38wLK0VuK2hGNMFvrL4nomPLYh3WHEpLeP4QDc51Gpec8INkyLxOjS3Ig3eku1uBOOiOSKyMci8rx7e6yIzBORahF5zN1+2nioae7glgXb+cMnO7qVrzAmFFXduZ10NopyN4Ue6ZRRRcw5Y9fmZ+X5wq8Pjn9B7SmjChlV4iSuaWO6t2wzoasyESn1amAp0LWL0++A21T1URG5G7gcuCsB79PjnPfaVuZvdlZ+f17bxj1HBe47Yoy/aKoFjM/AnTNDLV/80cTE7CKaDfYbkM+mS4exqKaNUSW5VBbGd/3/zxMGcPTQXeMzP5pYyqtrWmhsV/rlC9dnwGSMuL5DERkBnArc594W4FjgSfchs4Az43mPnmp1ffvOZAPwWJR70Jve650NkS/2vG9K5l3EhFoLNLo08xJkMuXnCgcOzGdgUW5cxTTvPao/xw0v7DY+M3FAPu+cMYj7pvTnw68PZlAGjI3F26V2O3A9uy5aBgC1qtru3l4DDI/zPXqkZ4LswmdMOA98Efl6jH0yoLKAvwt3D961E+9Vfm8VbKO23frlcfa4YgZHuCVCssV8OSEipwGbVPXfInJ012GPh4bssa2uro41hKRJRUy/+KgI/49rwdJqornAy8TPzovFGVpbJ9y6LJ8PtuXyj0nNDCoI/iezpkn4rLbd875cUW7fq4VZa/owuED54djWtH72od+72PPomKY1pDrkzPr99P5cwtq6huqmzB8ci6f9ejhwhoicAhTijOHcDpSLSJ7byhkBrAv1IlVVVXGEkHjV1dUpiUnnrg04VjJ0DFXlkV2RpirOeGV7nHM3tDBvYysnjSxMWv2p8pm7fhdOnV/EiguH0j/IFevBMwN/b7r88qB+XLJPKZdMTniIUQv7c/f4/QeYsMfuKd2nJeN+P4N8LqFMHV7A1P0yoyMpXPKOuf2qqj9V1RGqOgY4H3hDVS8C3gTOdh82HZgd63v0NoW9vMZUpnlnfQunvbiFXy/YzuGzNzF7ZRPNCZ5NqB5Tti5+Y6vnY59eEXpflJNGZtZ6m1j09p1Ez9utKOLHvnxKJY8cV8Ejx3uvuclEyegwvQG4VkSW4Yzp/C0J79Ej1TQnrvy4id+Vc7d1uz39zRqGPLCOhlhr/Hv4o8fe9u96TApYXd/OjLe2BRz3VRrDCnaTWa6ZWMq40tyI1iodMriAk0cVkZdFC5sS8huqqm+p6mnu18tVdbKq7q6q56hq9pc4TZFf/TvynRtN8qgqCza3srreeyX/v1Y1J+y9Pt4SOMV5VEngAO+9S8NPFMii844JYnx5HxacPYTNlw5LdyhJYZdEGeSNdZab02lTUweN7Z1cObeWY5/fHPRx//ufHQl7z2dWBs5WvGzPwFlcXi0hf+GKYprsEa78zMAsnc3Xuya9GxPELdX5zJ4bWamYc8ZF3s8eSlun93iQ/6kkkh0ijxpaEHSigclOY0tzWbHDu5V95xGz48gFAAAcvUlEQVTpq/gcD/sNTZNgFXNbOzJ/amNPs3BLK7M3Rn7t1TdBYyVvBWnR/nlxPR2dyqc1bWxp7mDYA+tDvs6tk/vx8HGZt8DTxOeBY4NPBjghSyeIWAsnTfYo78OnHiVKBv1jHZ+eM5gRJfajSZVnv4puEe7PPqzjgMo+HDo4vjLv57zqPRttS3MnA2aFXE0AOBt5vXHaoLRtF5xovzsk/lpiPUmwRbtfy4BdW2NlLZw0qGnu8Ew2Xf68OHx/vUmcHW3RtypPemEL9THOVlNVaprj31Br8blDekyyOWpoAVdMSH9xyWxQlaCq0ulgCScN7gkz4+juJYnZStZE5tU1sc06i+XntHBLK/3vX8e4R+LbWmDv/nlUZsDunYlyzriiXr8GJ1J9+2Tv52QJJw1+tzBxs5xMfNY3drAyyMBsOLcs2M4ej67nhg9q2dLcwd8/a+DNtcGT11c72jn6ueCz36LxwikDwz/I9EiLa7xLG2WD7G2bGZMAU0NMf47EpqZO7lna0K3V+qfDy7lkj13dQ39dUs9dS+qDzjiK1soLh9Iv3D7NWWZ4357TWku2Oeuzd/lEz/qtNSYKtS2drGlITBLw9f13a3d+vXJHO9fPq0tYsgHol5+9XSpdbtx/194sY0pzOXpYfBMwTHawFo7plVo7lCnPbkrqe9S3dTLpyY0Jfc0JJR09Yqzjuv1KGd43l/WNnVy6R3HGbomcbkJguf0r987eyRVZmXBaO5RnVjZRli+cOKKwR/wB9lbtncqVc7fx0qpmnj2pkkmVqZny+dSKJr4KUromUbyqCMTrvonZ253iS0S4sCp7T5ypUtJHAmZRHjkke1uDWdmlduHrW7ni7W2c/1oNv/nYBuCzUacqz65sonLWOh7/sontbcrRz23m4y2R72YZj7dT0A9+44d1CX29104bGHTBsOmZvC6lE1ywPKWy7td3TX07r63ddbL4fQLrWqVCp0c5+p5u6bY2fvZhLXcvqd9ZzqXi/nVc+mZNwGO/+07oisiJEs0amp8fUBbTe9S1JvZnfWBlz1hzYyJ3/IjsrCgQTNZ1qXmdJ9o7NS0luhvaOrnugzo+2NjCmWOKuOnAsm590arKzM8beXNdMyeOLOSi3Ys9449EW6fyn61tjCnNzar1F5/XCxfP3TVWsmBzK/uH6DYLtptlIk18YgOrIuhOu7iqmJ/tX8awvrncsiC5lby9+urPHVfE48udbrl7jupvXce90E0HlPGU33b0xw3P3i61rEs4Y8vy6Jcv3a4ea1s703ISfuzLJh5Z5myKdduieo4bUcgRPv2r18+r21lW/rmvmhlbmsekAdFfpbZ0KIP/savUyS0Hl3FilvzOXbywe6HLx5c37TyJBvNlXTu7JWk19dqGjoiSDcCBlfkMc6frDi/OZW1j5GM+HUEKcwYz54yBHPVs9ynaf51SwV+nRPUypocZW5bHk1MHcLZbBmnWMRUU52Vdx9ROWRl5hV9V3Fg2LmuP8oTg5dr3a7vdvuGDXbeb2zVgD5NTX9wSU//rY1923+nx5/O30+Rz7rtvaT3lM9dy8etb2djY4bmLZDb5++e7PrdOVWZ93sBN8+v4si7+1s8HGyMfuznLpyr0p+cO5rr9SjltVCEfnTUo7HNfWB18AeiUoQWMK911gXTV3iVMHJDPwQN3XYw8erwV4zSO40cUMv+IRmpnDGfamMRUKk+XrGvhgJNwfNc11LZGl3Aeqm7gqrlOcvjW+L78/rDyhMTV5JNN3lrvfcLxv/LtXyD8/tByvjmn+9jFt+bUcOjgfH70vvfA85ytuUwEHviiges+cB7z/Kpmnl+1q2TKI8dVcPKo9P2CxjpedefievLE2Rfm+VVN3DTf6c7606f1rLtkaFxXeFsivDjply/dFleKSLexnK3Th7G6oYO1DR388L1aqv2S4SVvBI5PdbnpwDIOGpjP9tZOOhXK3QuoV04dyNqGDsrycyiz2QGmB4r5t1pERorImyKyVEQWi8jV7vEKEXlVRKrd/xO+cUNRXve+7KYoLnxbOnRnsgG497MGFoUopBmN5Ts6du5dcsci7wKc/le+eSJ4nZefWN4UNNkA3PSF06fmu8jQ3wWv19DQ1smGxg5W7kh9OYx74qgJd8en9ez/z407k02XC18PfiKPxA3zvD/Tfx7YxG8m92N4cS5Thhbw3pmDQ75Obo4wpjSPw4cU8OHXw7d4AMrzhX+dXMlBA50xrLL8nJ3JBpykNqIkz5KN6bHiaeG0Az9S1QUiUgr8W0ReBS4DXlfV34rIT4CfADfEH+ou/kOn0bRwvKo0Hzl7E1umD0vIxINwe5d8b273BJGXE/s0x5vnh592O/zBwHhWXDg07GZdvq3AC3cv5jeT+3U7OQK8sKqJC1+vYdqYQu4/uqLboHZrh/LTBE8LBmcPmbrWzphKu9yzJHgV7gF9lCsnlHDl3iVRv24kg/n9C4QVF/bMbYONiVTMl1Kqul5VF7hf7wCWAsOBacAs92GzgDPjDdLfOxu6r9V43G+MI5i61k6OC1I7q3LWOt6JYm3G7AQt6lvf2Bl058dw7ohg22EvYx9ez+r64C2eRTVt3VqBDy9r5Ntvd29ZvLO+ZWdrY/bKZvrfv46fzKvlu+9sY1ldW9DP2Uu0031/HmMiC9a6Aeib5M7lcmu1GJOYSQMiMgbYH5gHDFbV9eAkJSCy/oY4/GtV8AHadze0MOPNGr7+8hZGPxS69XH6S1uIpOfp6RWNTPdYQxKrxjSs5Nr3iY3M3dDCgs2tAWMtT3ok8JfXtPDXJfXM39TKQ9UNnP7SloDH3L2kgUeWNXLQU5si7qacfWIlr5waXeXjB6oju8Do0tDWya0hpjV/a3z8K97/dHjoccDdy7JyuNSYhJJ4ZzSJSAkwB7hVVZ8SkVpVLfe5f5uqdhvHqaur2/mm1dXVUb/nwXOLA47NPyLwJLShRTh9fvSD5scNaOf1rc4J4uoxrVw8YlcWqm2D0+YX0dKZuDURL05u5OQPA7+nVHrz0EZK8kAVJr+bmlhyUOYd4bQUvX6mobx/eCN5EfwIVOGo94toDvLzOnlgOzfv0UpuAn6cob6HO/dpZnJ5jIuwjMkiVVVVO7/u169ft7+suC67RKQP8E/gIVV9yj28UUSGqup6ERkKhKyQ6BtcpA75YjPzNu3qVhtYmOP5Os8s3A5EX4mgK9kA3LEyn/1GD+ascc7J5Mfv19LSmdgN0g7bu4qj12wJusc9wLbLhtH//tDbDi88ezBjSvMon7k26hiO+aCY2hnD3andqdkA7t0zB1PV3+lOe6+yja89E3kxzcPeLd75/Yby/sYWmjsDW2NAtxlv1dXVMf0udjM3+Od+4UHjErJwMyFxpoDFmVjZFGco8cxSE+BvwFJV/YPPXc8C092vpwOzY32PYM7drXurZXNzJ+Uz11I+cy3vbWhhW0snLR3KKzHu5OjvG3O20djeiapy72eJPRl3de1fs2/gYHWuwIkjCthwyTBEhPHlwU+uW6cP23nyff202DbnamjrDLsbaTzOGF3IlKEFfHH+EGpnDGdC/11jN3v178OJI6JbzTrpyY28tDr0WNrVIWbxpWoBXY5ENrHAmJ4unr+4w4FLgGNFZKH77xTgt8BUEakGprq3E+qjzcHHB055cQtjH17PuIfXMz/E45ZdMIRnT6rkmAj34Rj2wHoeXhbd2MG2y4ZRO2M4M/YM3tVy29ec3schxYGVErZeNpzHplZS6PYd/W2K92LA83YrItdnht2BA/P5nt9sq/N2K+KwwaErMXvNaEukWcdUMPukSgYVeVeFeGxqJQvPHsyfj+g+HrLx0uCzu85/rWbnbL3rP6ilfOZaDvznhp3T078Islj0j2HGXGIxqsT7+3p3WtKHMY3JCjF3qanqXLyLmQIcF+vrRiKSvvuGEAPx704bRGVhLkcNzeWooQURd0FdNdf7avlHE0v49l4l7PHorkWXW6YP23lVe+3EUmZ+7p2sxrqtkj3L+3DIoPydXYXX7Vca8Ni9K/rwj2MquhW9nHPGQPYbEJhIbpncj+/uXcLvF27niKEFnO12CcbS3ZYIU4cXRHSVP6Y0jzGleVzsV7p+zcVDGREkId7xaT3/qWnb2SX55fYOTntxC7NPqgz6Phftnvhxqn+dXMm+T3Tf/+bh4yq6teSM6c2ycurM4CBXyJE4aGAf9q5I7Ang5weUISLUzhjuef/Ikjxm7FnsmXQm+UwJfuqEATz2ZRMVBTlMG+NdJfaMMUXUzhgeUZ/u8L653H5493W3r502kOPj3Fa5S9f3GyyJPTF1AD+eu5mzqsq4KcaKy11K+oRujPuPfy3Y0sZ7G7y3OjhxREG3FmGijCzJY/G5Q3h7fQuHDc4PO75kTG+TlX8RM8b35X8/iX4ywMkjC/nVwYEnvrGluTFvARwsyfi77Wv9PROO7zhC3z45fCMBU3RDOWhgPrUzhlPT3MG179eF3SRs/SXDeOCLBnIExvfvw5XvbOOC3Yv5sU8LrHbGcFQVEWFLcwfrGjrYp6IPOSI8cWAzVVUjExL7/zukH9eHWEvj77zXtgYcG1USmIQTaXjfXC5IQuvJmJ4gK1ejDe8bfQvn1weV8cjxA6jqF9i6+cuRsZ2AfnFgdFftW6cP21l4tH+BMD+CIpDJUlGYy/3HhC4QeUlVMUV5whV7lfDNCSUcMaSAT84Zwk/3LwuoytDVXVZZmMvEAflJ2TI43mQ8eWA+n5wzhKEe42XGmOTLyhZOLL4bomTJYYMLeOHkSt7b2MrY9k1o/6EcM6yAvy5t4HcLg7ekoh0HyM0Rll84lG0tnfTNE/ITsfgjibzGkdIpL0dYeeFQxjwc2+SGb02wLY2NSaesbOFEY1BRDssvGBK2TtrXhhRw3X6lTCzr5OxxxQwozOWn+5fx8HHBWwGDimL7+PoX5GRMsqk+f4jn8dsOK2d0Bo5B9MsXyvNj++z2i2EvImNM4mTeGSVCxXkStiTMg8dWcNjgfCri2JztlFFFfHbeEMY/tqHb8dknVvaItRUDi3J3jkPVtTp13TJ5R1ER4eVTB3LI05EvEu0yIsi0ZWNMamRtC+eZEweEfcxpo4sYkICT55DiXO72GeeZf9YgpkS4fieb9MvPyehk02XP8j7UzhjO9ZOi6/LL5p0SjekJsvYvcPKg0Cf8RHefnL+7U/qldsZwz4kHJvWu9xtjeu/MQUGrMdziMTvRGJNaWdulBvDDfUu4PchGZ/dNSd7UV5MZ8nIC1z7dO6WCI2d372574eRKvjak57VIjck2WZ1wvrtXCYtq2liyrY3Lx5dwzb4lrKrvYHBxjnWf9FL7VvRhwyXDmLO+hQn98xhVktW/4sb0KFn91zi4OJd/ntC9fMlY23ek1yvME04c6V2pwRiTPtYMMMYYkxKWcIwxxqSEJRxjjDEpYQnHGGNMSljCMcYYkxKiGro8TDLU1dWl/k2NMcakVL9+/brV/7IWjjHGmJSwhGOMMSYl0tKlZowxpvexFo4xxpiU6LUJR7JkM5tsiFNEsqqeUKZ/piJS7P6f6XFmRdn0TP8cu4jI3iKS8TWZRCTX/T/qz7VXJRwR2VdEzhaRIs3gvkQRmSAihwFkeJyHici9wMHpjiUUETlCRO4SkSshMz9TEckRkQoReQX4MWRmnAAicqiIPAr8XkT2SXc8wYjIIe7v5w0iMjDd8QQjIhNFZC5wCxB+o680EZHDRWQW8HMRqYjl97NXJBwRKXB/8R4ALgF+IyKj0hxWABHp58b5KPBrEblVRHZPd1xeRORbwF+BBcDHXVc9mUZEDgDuAv4NnCIit4nIpDSHFUBVO4F2oB8wTkSOh8y7OheRc3A+z+eBQuBa93jGxCkiuSLyPzi/n+8CBwA3i8jg9EYW1M+BJ1X166q6FjLr8wQQkXHAX4A3gdE456dTo32dXpFwgClAP1WdBHwD2ANoTG9Inn6MM5FjP+DbOFc7Y9IaUXCjgBtV9S5VbVbVjnQHFMRkYL6q3gd8E+fnfoqIVIZ+WlrsBWwA3gFOz9CWeBXwnKo+CNwGTtdahsWZA6wCzlHV+4EfAocCRekMyp/bqt0NqFfV291jU0WkHIi52ypJDgaWup/nj4CFwGkiMjKaF+mxCUdEDhCRPd2brcAx7tdH41xFHisiI9IRmy8RGSsiXX8I9wK/AFDVL4FyYN90xebLjbPA/boC2Af4UESOFZGXReRnInKWe3/a/khE5FwRuVZEvuYeWgCUiMgQVd0AvAFUAoenK0boFuehPoe/AhYDXwCdwEkiMiQtAbp84jzMPfQ5cJaIXA+8DwwD7hSRtHarut18e7g3O4FHVPULESlQ1XXAGpyfe1r5xum2ajcBR4rIqSLyDHAd8EfS3K0qIqeLyPd8fj/nAyNFZKSqbsNpOdYCX4/mdXtcwnFPjP8C7gT+ISLHqepbwCMiMhunO+B+4AzgJ+lKOiIyRkReBO4DHhSRPVX1K1VdJyL57sOagC/TEV8XvzgfFpEJqloDbAUeAs7EaWqvB34hIvul44/E7Ub5BXCDe+geETkdaABW4rRyAeYAdcBI93kpTY4ecd7blaiBSUBfVX0b54/5T8AtIpKXIXGeATwFXA0cBVyqqicBm4H/SkdyFJFy9+/9VeBcESlR1Q5VrQVQ1RYRKQXGAutSHV+IOPu68e0AZgK/Bv6uqifi/K0d6ncxkqo4h4rIczgJrz8wU0ROVNXlOBcY57oP/RxYAgyQKCY69IiE4/fHeB2wUFUPA2bjdKMAXAOsAE5wu1f+BygA9iRFPOKcp6rH4fSL/lpE9nbv6+qeGg6sdp+bsp9ViDjfwDkBjgVuxml9rVPV2ao6E3gBmJaqOH25XXp7Aj9S1T8A/w18H2eTwfXAJBHZS1Xbcf5Yvu4+L6XJ0SPOm4EfuFe964AGEZkJzMBp6Xyiqu0ZEuc1wB6q+jrQjPM5gvN3NhEnuadaX+BlnJ91X+BIj8ccAix2L+ZKRKQqlQG6/OM8yue+53G6zvu7tz8CNgItKYyvy0HAXFU9SlV/DdwBfMu9by6wr4hMdn8/1gKHq2pzpC/eIxIOzuBl14myAWhzj5cBS9wTTQewBTgJQFUX41zlrklDnF3TiBe7sfwZZ6zhQhEZpKod4kwWqFHVj0Xku8BNbt9uOuO8EzgQuALnqvY+4Gyf5w0C3ktRjIjIpSIyxedz2Qj0F5E8VX0Sp3V4PNB1grzFfdxwYL6kaDp3mDifwvl8pwEDgROAHcB+wO+B/UVkTAbE+U83zvPdlsyX7PrZ74/z+aaET5xl7iD7X4HH3RgOEZFh7uO6fr7lwGoRmYHTNZSSSSMRxDkcQFU/wWlRfE+cscWLcbqst6YwzqPdLvPXgX/43L0VqHa//gD4GLhdREqAvYFV4k7jj+i9MmusLzoiMhW4HudK621Vfdxt9p+HMwArwDPAKThXaIXAz4AXccZ0VuB0D9Qm8woySJy/wrn6fsx92K043Si3qupSETkBp/tvFc4v6A9V9fPAV095nL9x47xZVZeJyFM4V+JH41yhX6Wq65MYowBDgIdx+uq/xLli/DbwAzfWP6pqrYiMx5nxd5KqbhCRvwODcRLjBaq6LEPinOA+7gSgRVW3u68xFGhX1c0ZEmfX5zkVp0VzFc4YTj3wPVX9LA1xXq2qW9zHHI7T5TPfndTQ9dwHgIuAWcBt7gk+E+L8SFUf8HnutcA4nIkZ16jqknTFKc5EkDYR+QGwl6p+x+e5fwBG4MxWuzSq85KqZuU/YHdgHs6V4f7uB3ede9+ewFM+j70Z+L379ZHu7bPSFOcjwJVAKXATTnN6Lk5T9mHgB+7zLgJqgOMzNM5r3OeVAeNxuiqTHWOu+/8ewIPu13k4Y0h/w7mSfRmnu6LYvf9xn1j7AAMzOM6r3a9zgJwMjfMJ4Er36xJg3zTG+Sffv3P3+DU4LdkyoMQ9dj5wdobG2Q8o9TneJxPi9HnMc13nIGCQz2NLY3nvbFshngM7Z3ccAvxbVWe7970G/MG9mqnBaUJPUNWlOM3EH4pIjqq+gzPtNJ1x/h/whKr+WkTGqTMgh4i8y65+20dV9aEsiHOHOle2yby6zQN+BeSKyAs4J5MON/Z2EfkeznTiP+Akw/OBoTitsjbcbj5VbcPpCszUOD9wH9uZrBgTEGcrzpomVLUeWJTGOH8ArBORKao6x33avTgn8teBUSIySVUfTVaMCYjzVWC0iOyvquvc39GMiFOcyUubgS9E5FacadBHqzNLbUcsMWTNGI7b/7oGZzYHOL/oF/j0b/cBlrv37wAqcAZjrwbuAV4D1G9APB1x5uE0X29zb69wn3cFcDnONF40yetaEhhnUvtkRWQKzgmuP7DMjbcNOEZEJrsxdOJMEvi9qs4CXgEuFZGP3e8jaSdFizOtcSrOCfSXPk89FadlvhCn9ZW07t0ExfkfN86kzqCLMs7/dp9WCFyGk7xLcVo62+IKJNnNtwQ1AUtwxmKuxjnRjXeP347T9fMu8CDOrKkXcfoiJ+DMCJkFHJqBcf4LGOze/0OcwcyDLc6AWI8ELvG5/Rfgu+4fwr/dYzk4/dFPAiPdY0OAcRZnr4jzcWCMe2wacJTFGXecI3AmMv0DmJSwOFL1DSfgAxvl/v9b4DH361yclswR7u2ROAkmPwvivB8ocG8XW5xB4yzGmb7e1ad8EfA/7tcLge+7Xx+Es9gvXT93i9Pi7ClxPpqsOLKmS01VV7lf3g6MFWcxUgdQp6pz3fu+gzMtOm1lVqKIsxGndhaqmvIyO1kUZ6OqtuiuLsap7BqHmQFMEJHncVpmC1IdXxeLM7FiiTPZ3eW9KM5/Jy3OdGXcOLP1t4E5Prcn4yw+ewEYku74LM6kxJiL0+R/EdjdPbY7zkyqI4Dh6Y7R4rQ4Lc7Q/7JuHY4706xTRJ7EWUHegjMhoFqd+mMZweJMLPdqKx9nsenTOEVYt+J0BWxPZ2y+LM7EsjgTK+1xpjvjxpili4G3cSoH/CDd8VicKYvzUJxFanOBy9Mdj8VpcVqc0f3LuhYOgIhchzOL4gZVTUe9oYhYnIklTqHVS4A/WJzxszgTy+KM4L2zNOHkaJIXxyWCxWmMMbtkZcIxxhiTfbJmWrQxxpjsZgnHGGNMSljCMcYYkxKWcIwxxqSEJRxjjDEpYQnHmCiIyEoRaRKRHSJSKyLvich3uvYWCvPcMSKikqKtrY3JNJZwjIne6apairPF7m+BG3B2yDTGhGAJx5gYqWqdqj4LnAdMF5F9RORUEflYRLaLyGoR+aXPU952/68VkXoROQxARL4hIktFZJuIvCwio1P8rRiTEpZwjImTqn6Is3vqkTjbY1yKU333VOC7InKm+9Cj3P/LVbVEVd937/sZcBYwEGf780dSGb8xqWIJx5jEWAdUqOpbqrpIVTtV9ROc5DElxPO+jbMR1lJVbQd+A0yyVo7piSzhGJMYw4EaETlERN4Ukc0iUoeziV1liOeNBu5wJyDUAjWAuK9nTI9iCceYOInIwTgJYi7wMPAsMFJV+wF34yQQAK/ChauBb6tquc+/IlV9LxWxG5NKlnCMiZGIlInIacCjwIOquggoBWpUtVlEJgMX+jxlM84+JON8jt0N/FRE9nZfs5+InJOa78CY1LL1AMZE7zkRacdJHkuAP+AkDoArgf8TkT8Dc4DHcSYQoKqNInIr8K6I9AFOUtWnRaQEeNQdt6kDXgWeSOl3ZEwK2PYExhhjUsK61IwxxqSEJRxjjDEpYQnHGGNMSljCMcYYkxKWcIwxxqSEJRxjjDEpYQnHGGNMSljCMcYYkxKWcIwxxqTE/wcxYML58DeS9QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df['Closing Value'].plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "d1 = pd.datetime(2000, 1, 1).toordinal()\n",
    "d2 = pd.datetime(2015, 1, 1).toordinal()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['ord'] = [x.toordinal() for x in df.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc3 as pm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Input dimension mis-match. (input[0].shape[0] = 3, input[1].shape[0] = 8223)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-61-9955221579d5>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      4\u001b[0m     \u001b[0mT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mUniform\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'changepoint'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0md1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0md2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mshape\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m3\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m     \u001b[0mu\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mHalfNormal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'u'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mshape\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 6\u001b[1;33m     \u001b[0mgrp\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'ord'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvalues\u001b[0m \u001b[1;33m>\u001b[0m \u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m*\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      7\u001b[0m     \u001b[0my_obs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mNormal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'y_obs'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmu\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mu\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mgrp\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mobserved\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'Closing Value'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvalues\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\theano\\tensor\\var.py\u001b[0m in \u001b[0;36m__lt__\u001b[1;34m(self, other)\u001b[0m\n\u001b[0;32m     52\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     53\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0m__lt__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mother\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 54\u001b[1;33m         \u001b[0mrval\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtheano\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtensor\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbasic\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlt\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mother\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     55\u001b[0m         \u001b[0mrval\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_is_nonzero\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mFalse\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     56\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mrval\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\theano\\gof\\op.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, *inputs, **kwargs)\u001b[0m\n\u001b[0;32m    672\u001b[0m                 \u001b[0mthunk\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0moutputs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mstorage_map\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mv\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mnode\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0moutputs\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    673\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 674\u001b[1;33m                 \u001b[0mrequired\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mthunk\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    675\u001b[0m                 \u001b[1;32massert\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0mrequired\u001b[0m  \u001b[1;31m# We provided all inputs\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    676\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\theano\\gof\\op.py\u001b[0m in \u001b[0;36mrval\u001b[1;34m()\u001b[0m\n\u001b[0;32m    860\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    861\u001b[0m         \u001b[1;32mdef\u001b[0m \u001b[0mrval\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 862\u001b[1;33m             \u001b[0mthunk\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    863\u001b[0m             \u001b[1;32mfor\u001b[0m \u001b[0mo\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mnode\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0moutputs\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    864\u001b[0m                 \u001b[0mcompute_map\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mo\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\theano\\gof\\cc.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m   1733\u001b[0m                 \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0merror_storage\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfile\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msys\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstderr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1734\u001b[0m                 \u001b[1;32mraise\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1735\u001b[1;33m             \u001b[0mreraise\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mexc_type\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mexc_value\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mexc_trace\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1736\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1737\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\six.py\u001b[0m in \u001b[0;36mreraise\u001b[1;34m(tp, value, tb)\u001b[0m\n\u001b[0;32m    691\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mvalue\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__traceback__\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0mtb\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    692\u001b[0m                 \u001b[1;32mraise\u001b[0m \u001b[0mvalue\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwith_traceback\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtb\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 693\u001b[1;33m             \u001b[1;32mraise\u001b[0m \u001b[0mvalue\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    694\u001b[0m         \u001b[1;32mfinally\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    695\u001b[0m             \u001b[0mvalue\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: Input dimension mis-match. (input[0].shape[0] = 3, input[1].shape[0] = 8223)"
     ]
    }
   ],
   "source": [
    "mod = pm.Model()\n",
    "\n",
    "with mod:\n",
    "    T = pm.Uniform('changepoint', d1, d2, shape = 3)    \n",
    "    u = pm.HalfNormal('u', shape=4)\n",
    "    grp = (df['ord'].values > T) * 1\n",
    "    y_obs = pm.Normal('y_obs', mu=u[grp], observed = df['Closing Value'].values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with mod:\n",
    "    step = pm.Slice()\n",
    "    trace = pm.sample(step=step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\willk\\Anaconda3\\lib\\site-packages\\matplotlib\\axes\\_base.py:3604: MatplotlibDeprecationWarning: \n",
      "The `ymin` argument was deprecated in Matplotlib 3.0 and will be removed in 3.2. Use `bottom` instead.\n",
      "  alternative='`bottom`', obj_type='argument')\n",
      "C:\\Users\\willk\\Anaconda3\\lib\\site-packages\\matplotlib\\axes\\_base.py:3604: MatplotlibDeprecationWarning: \n",
      "The `ymin` argument was deprecated in Matplotlib 3.0 and will be removed in 3.2. Use `bottom` instead.\n",
      "  alternative='`bottom`', obj_type='argument')\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[<matplotlib.axes._subplots.AxesSubplot object at 0x0000026BA7303A90>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x0000026BA727F908>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x0000026BA7330AC8>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x0000026BA735A128>]],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAEYCAYAAABBWFftAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8XXWd//HXO0nTNpSUHcpSCkNFFhWBAQRRBGRTARcElAEVRVxBnd8IIovFBXBfGHAGOywji4pIZRksmyhShrIvAimdAoVSVnvbpk1yk8/vj/NNenqbpPe2SW6S+34+HueRc77n+z33c79Neu73nu+iiMDMzMzMzMzWXl21AzAzMzMzMxst3MAyMzMzMzMbIG5gmZmZmZmZDRA3sMzMzMzMzAaIG1hmZmZmZmYDxA0sMzMzMzOzAeIGllkiaZ6kb1Y7jmqQFJKOq3YcZmY2cHxf833NqsMNLDMDmAT8rpICkuZIOmdwwjEzM1srvq9Z1TRUOwAzq76IeKnaMZiZmQ0U39esmvwEy2qOpC9IekJSm6SXJeW/4WqU9FNJr0taKOkHkupzZd8r6c50fpGkP0vao+T6Ienzkq6QtFjS85L+rSTPhpJ+K2lpep1zJV0m6daSfF+S9KSk5ZJaJJ0hqSF3fp6k70i6RFJB0quSzpdUl8uzrqRfSnolXWe2pIN6ifm4kuM+34OkO4F/As5OeUPSlMr+JczMbCD4vub7mg0vbmBZTZH0LeB84N+BtwCHAA/lsnwJWADsCXwZOBU4Pnd+AnAhsBewN9AC/I+kDUte6mzgLmAX4PvA+ZLekzv/X8DbgPcD+wNbAkeWxHoO8K/A6cAOwCnAZ9O1874EvAj8M/AV4Isp7m7TgYOB44C3A3cDN0h68yoVVP57+BAwD/ghWTeMScDzq7memZkNMN/XfF+zYSgivHmriQ1YB1gG/Gsf5+cBM0rS/ge4qp9r1gFvAB/PpQXws5J8TwLfS/tTU54DcufHkP1Hfms6bgJagUNKrnM88I+SmP9Skue7wPy0v116rcNK8jwATC+J+bhy30M6ngOcU+1/V2/evHmr1c33tZXy+L7mbdhsfoJltWQnYBzwp37yPFRy/AKwafeBpG1S94I5kgpAAZgIbF3BdXZMP2d1n4yIDmB2SazjgWslLenegF8CEyVtnMt7T8lr3Q1sIak591p3leS5K71Gf/qtCzMzqzrf11bwfc2GDU9yYbUo+jnX3kve/BcRNwCvAl8g+2auHfgr0FjhdVYXR3feo4Cnezn/ej9l1c+5fJ7+Xh/Kew9mZlZ9vq/5vmbDiH+prJY8ASwn67ddsdQffUfgvIi4JSK6r7fJGsQB8I7ctRuA3XJ5Hk/X3jYi5vSydeby7lVy/XcAL0ZEIV0H4F0lefbNnVtT7UD9anOZmdlg8X1tBd/XbNjwEyyrGRGxRNIPgXMkLQNmknVXOCwivlfGJd4AXgE+I+kZYEPgArL+75XE0SLpj8CFkj6brvk1oJn07VuK9bvAdyWRYm0gG8D89oj4eu6Su6SBw1cCu5MNGj4nXecZSb8F/j291rPA54CdgY9VEncv/g/YR9Jksn71r0dE11pe08zMyuT7mu9rNjz5CZbVmjOBM8hmUnqMrN/6ruUUTP/JHkU2jesjwKXAT8hmZ6rUJ9Pr3wzcSdYPfCbZt3vdr3cu2exJnwYeJuuy8RWyAcB5PyfrKz8b+AVwEfDj3PlPA7cA/52usw/w/oh4cg3izjubrJ/+U2Q308lreT0zM6uc72u+r9kwo4jVdVc1s8GW1iR5kmy2p69VUG4ecElEfHuwYjMzM6uU72tWy9xF0KwKJL2LrI/7g8C6ZN/gTSH79tDMzGxE8X3NbIUh6SIoaXpaWfyxPs5L0s/SFKGPSNo1d+6EtNJ3i6QThiJesyFQD3yTrGvDHcC2wHsi4tGqRmVmZrZmfF8zS4aki2D6VmMJcHlE7NzL+cPIVu0+jGyl8Z9GxJ6SNiDrf7s72SDJ+4HdIuKNQQ/azMxGNUnbA9fkkrYFziKbsWz7lLYe2SKou0iaAvydbHwGwKyIODldazeyb+rHAzcBp0REpPvYNWTf5M8DPhoRbygb5f9TsvteK/CJiHhgUN6omZkNqSHpIhgRd6UbU1+OIGt8BTBL0nqSJgH7ATMj4nUASTOBQ4Cr8oUXLVrkgWRmZiPAxIkTy1nPZkhExFPALtAzXuQF4LqI+El3njRD26JcsWciYpdeLncRcBLZQqs3kd2rbgZOA26LiPMknZaOvw4cCkxN256p/J7dF/N9zcxs5Ci9tw2XWQS3IFvcrtv8lNZXupmZ2UA6gKzx9Gx3QnrK9FFKvtQrlb4QbI6Ie9IXhZcDR6bTRwCXpf3LStIvj8wsoPuLRTMzG+GGSwOrt280o590MzOzgXQMqzak9gUWRkRLLm0bSQ9K+rOkfVPaFmRfAHbLfxm4aUQsAEg/N8mV8ReIZmaj0HCZRXA+sFXueEvgxZS+X0n6nUMWlZmZjXqSGoHDgdNLTh3Lyo2uBcDkiHgtjbn6g6SdWLMvA8su09LS0luymZlV0dSpU/s8N1waWDOAL0q6mqwP+qKIWCDpFrIVv9dP+Q5i1RugmZnZ2jgUeCAiFnYnSGoAPgTs1p0WEW1AW9q/X9IzwJvIvgzcMne97i8JARZKmpTuaZOAl1N6X18srqK/m/jqtLS0rFX50cL14Dro5nrIuB4yg1UPQzVN+1XAPcD2kuZLOlHSyZJOTlluAuYCc4D/BD4PkCa3OBe4L23Tuie8MDMzGyClT6oADgSejIiern+SNk6TYSBpW7IJKuamrn+LJe2Vxm0dD1yfis0AupcYOaEk/fi0TMlepC8WB+G9mZnZEBuqWQSPXc35AL7Qx7npwPTBiMvMzGqbpCbgvcBnS071NibrXcA0SUWgEzg596Xf51gxTfvNaQM4D/iNpBOB54CjUvpNZFO0zyGbpv2TA/SWzMysyoZLF0EzM7MhFxGtwIa9pH+il7RrgWv7uM5sYJV1HiPiNbIZCkvT+/xi0czMRrbhMougmZmZmZnZiOcGlpmZmZmZ2QBxA8vMzMzMzGyAuIFlZmYjlqQxkvaVdHQ6XkfSOtWOy8zMapcbWGZmNiJJegvwNNnyHr9Kye/GM8+amVkVuYFlZmYj1UXAWRHxZqAjpf0ZeGf1QjIzs1rnBpaZmY1UOwH/nfYDICKWkq1FZWZmVhVuYJmZ2Ug1D9gtnyBpD7LFe83MzKrCCw2bmdlIdSZwo6SLgUZJpwMnA5+pblhmZlbL/ATLzMxGpIi4ATgU2Jhs7NXWwIci4k9VDczMzGqan2CZmdmIFREPAJ+vdhxmZmbd3MAyM7MRSdK0vs5FxFlDGYuZmVk3N7DMzGyk2qrkeDOydbCuq0IsZmZmgBtYZmY2QkXEJ0vTJB0CHFuFcMzMzIAKJrmQ9GVJGw1mMGZmZmvpT8CR1Q7CzMxqVyWzCB4IzJN0g6SjJY2t5IUkHSLpKUlzJJ3Wy/kfS3oobU9L+kfuXGfu3IxKXtfMzEYnSduWbDsD3waer3ZsZmZWu8ruIhgRh0vaEDgGOBW4WNK1wOURcVd/ZSXVAxcC7wXmA/dJmhERT+Su/5Vc/i8Bb89dYllE7FJurGZmVhPmAAEoHbcCDwInVC0iMzOreRWtgxURr0XEhRHxDrKBxP8M3CFpnqQzJE3oo+gewJyImBsR7cDVwBH9vNSxwFWVxGZmZrUlIuoioj79rIuICRGxb0TcX+3YzMysdlW80LCkAyT9F3AnsBA4HvgXsidON/dRbAtW7rIxP6X1dv2tgW2A23PJ4yTNljRLkvvWm5mZmZnZsFR2F0FJPyDrHrgIuBz4ZkS8kDs/C3ijr+K9pEUfeY8BfhcRnbm0yRHxoqRtgdslPRoRz5Qbu5mZjQ6Snqfv+0ePiJg8BOGYmZmtopJp2scBH4yI+3o7GREdknbvo+x8Vl6vZEvgxT7yHgN8oeTaL6afcyXdSfa0zA0sM7Pac1y1AzAzM+tPJQ2s75ENIO4haX1gfK4B9GQfZe8DpkraBniBrBH1sdJMkrYH1gfuKXmN1ohoS9PE7wNcUEHcZmY2SkTEn6sdg5mZWX8qGYP1B7InT3lbAtetrmBEFIEvArcAfwd+ExGPS5om6fBc1mOBqyMi3/1jB2C2pIeBO4Dz8rMPmplZ7ZK0i6QvSfpWuqdMkzStzLLb55YAeUhSQdKpks6R9EIu/bBcmdPTciNPSTo4l97rUiSStpF0r6QWSddIakzpY9PxnHR+ysDVipmZVVMlT7C2j4hH8wkR8aikN5dTOCJuAm4qSTur5PicXsr9DXhLBXGamVkNkHQS8GOyxYUPJZto6SDg+nLKR8RTwC7pWvVkPSyuAz4J/DgiflDyejuS9cDYCdgcuFXSm9LpvpYiOT9d62pJFwMnAheln29ExHaSjkn5jl6jijAzs2GlkidYL0vaLp+Qjl8b2JDMzMzK8m/AIRHxQbL1Ej8IfAToWINrHQA8ExHP9pPnCLJeFm0R8X9k63DtQR9LkUgSsD/wu1T+MuDI3LUuS/u/Aw5I+c3MbISr5AnWdOBaSWcAc4F/As4FLhmMwMzMzFZjk4j4S9rvklQXETdL+vUaXOsYVl5/8YuSjgdmA1+LiDfIlheZlcuTX3KkdCmSPYENgX+kbvKl+XuWL4mIoqRFKf+rpYG1tLSswdsZuPKjhevBddDN9ZBxPWTWtB6mTp3a57lKGljnkX0r+AOyGQGfJ2tc/WiNojIzM1s78yVNiYh5wNNkT41eBdoruUgaF3U4cHpKuojsC8RIP38IfIq+lxzprTdI9JOf1ZxbSX838dVpaWlZq/KjhevBddDN9ZBxPWQGqx7KbmBFRBfw/bRZLfnb32jwtxw2ShSnToW99652GDYwLiCbCGkeMI2sq10j8OUKr3Mo8EBELATo/gkg6T+BG9Jhf0uO9Jb+KrCepIb0FCufv/ta8yU1ABOB1yuM28zMhqFKnmB1T6P+NmBCPj0ipg9kUDa8NLS00HTmmdUOw2xAtJ57LkU3sEaFiLg0t39zWtajMSKWVHipY8l1D5Q0KSIWpMMPAo+l/RnAlZJ+RDbJxVTgf8meRq2yFElEhKQ7yMaFXQ2cwIoJOGak43vS+dtLZtA1M7MRquwGlqRvAGcBD7PyelhBNj7LzMxsyEj6CfDriLgPIE0wUWn3wCay2f8+m0u+QNIuZPe3ed3n0vIivwGeAIrAFyKiM12neymSemB6RDyervV14GpJ3wYeBH6V0n8FXCFpDtmTq2MqidvMzIavSp5gnQrsERGPDFYwZmZmFRBwvaSlwJXAlWnq9bJFRCvZ5BL5tH/pJ/93gO/0kr7KUiQpfS7ZLIOl6cuBoyqJ1czMRoZKpmlfBjw5WIGYmZlVIiJOIRvX9Hmy8UyzJN0v6avVjczMzGpZJQ2sM4GfS5okqS6/DVZwZmZm/YmIroiYGRGfAnYmW5vRkzGZmVnVVNI4uhT4DNnMRx1pK7JmCzqamZmtNUkTJB0n6UayqdqLZJNHmJmZVUUlY7C2GbQozMzMKiTpt6Qp1slmATwhIlZZqNfMzGwoVbIO1rMAqUvgprkpbM3MzKphNvC1iHiu2oGYmZl1K7uLoKT1JF0JLAfmpLTD09SzZmZmQyoiznfjyszMhptKxmBdDCwCtmbFOiP3AEcPdFBmZmZmZmYjUSUNrAOAL6eugQEQEa8AmwxGYDbyFLu6aCoUeo6vb22lrlDgjuXLAVhYLLJeocDCYrEnz1sLBabmygDsWSiwcaHA5mn7fWvrSucmFgp0dXUB8GRHBxuk8rcvW8amuf26QoEvLV7cU/aAQoF/S8fFri4+VygwuVBgq0KBLQoFjiqJozenLF7M/ilfa2cnexcK7FMo0NrZycRCgS3TtlWhwCm517592TLG5c5vWSjw/3LnV+eu5cvZbjXxPdDezjn9XPN9hQKbFQoc1M91jiwUuDbV9waFAoXOzrJjPH3xYj7Qy7VfKRbZMdXx1yt4z31p7+riA4UCm6R/v5nLlq31Nc3MzMwGSiUNrEXARvkESZOBssZiSTpE0lOS5kg6rZfzn5D0iqSH0vbp3LkTJLWkzbNDDVMNdXUE9DR+flYs0gy8mo6/u2wZuwGbNmRD/y5bupSlwEsl13kOuGrMGF5sbuYg4Du5BtlzwLrAtanRdnt7O1umc3/t7GSr3P7mwJ8jesrOA/asrwfgsCVLuB+4t6mJ55ubuaepieVlvMcnIthB4oVikV2XLmUz4O7mZl7s6kLA/OZm5jc3c8u4cfwugu8tWdITz6658/Obm/n+uuuW8YqZu4pFtl1Nnmva2ngw937z7l6+nCeAl5qb+VNzc5/X+ENzMx9uamJuRwfjgOZUX+V4JIK3SKukn79sGVOBF5qbOb+M99yefl/6ctKSJSwCXpowgc/U1fGNDk9kamZmZsNHJbMIXgJcK+kMoE7SO4DvknUd7JekeuBC4L1k07zfJ2lGRDxRkvWaiPhiSdkNgLOB3cmenN2fyr5RQew2RBqB17q66Orq4glgD+CN9KH/DxH8YsyYnrzf7uxkWn09p3R28nh7Ozs1NrKks5PXgL0bGwHYtq6OJ9MH7iWdnRSAkySuLRY5Cnigs5Pt0vUe7epiam7/IOAG4I3OTtavr+cl4MCxY7mmtZXZwHPrrMOE1ICY3NDAH/tpeHSbB+wN7NXayockfpoaDLe1t7N5Lt+OjY18pK2Nm7u6OD3Fs0MZ9bdeocDOZA3JInBtYyPvGDeOh7u62Ck1Xn69dCmnd3bSAUwC7lxnHa5bvpxfRDAe2LJQ4LepHGRPzz7c0UFXOvdgUxP7t7ayiGz18GMkfr7uuiwsFnlLaysvNzczs72drVNMT3d0cPyyZbxEtibDefX1/Ms667CwWOSI1lbmAVOBl4EvNqz8X8oZixfz7xGsA+xVKHBtUxNHtbYyh+x35YcNDRzd1MSPlixhelcXY4E3gLl9/Fss7+riBuD+pibq6uo4sKGBi9rbe81rtUHShsBhwKSIuEDS5kBdRMyvcmjDwpOzL2HuI2OrHUbVtbe31Xw9uA4yroeM6wEO/vC0Qbt2JU+wzgd+Q9ZQGgNMB64HflpG2T2AORExNyLagauBI8p83YOBmRHxempUzQQOqSBuG0LjgJe7uvjusmXsBzSTNbCua22lHvjA+PEA/GTJEpqAj6+zDpsDd6WnELe1tbEe0FRfz2vFIn/s6uKM1Ci7ra2NTYFPjh3LrPR6T0HPU5M5wNvq6nr2d62rYyfg0mXL+N+2NiYAE+vruaRY5H3Q07gq9X8dHUzNdUPMewH4YQRfrqvraVwBzO7sXGUdg/FAd4e4OcAM6Oke+J5eutI92dFBAfhIXR3PNTezH/C91HhoAd5RX8+THR18tbOT340dy4LmZtYHzm5t5YR11mEKcMWYMcxvbu5pXAHsP348+wOfkpjf3My6dXXc2NTEc83NPNHUxJWpAXxb7mng7M5O3pz2P7psGV9raGBeczPTx4zh3NRt8OOtrbyN7KnYARL/B7xn7Mr/WX9n3XWZRNZQnNXczIdaW/lniZebm/l+QwPfSE8nH+rq4g3g901NzG1u7vPf4LJly9gM2Do15F5PjTerTZLeTfbfwMeBM1PyVOCiqgVlZmY1r5Jp2gP4SdoqtQXwfO54PrBnL/k+LOldZItFfiUinu+j7BZrEIMNgXFkXQJnRPDHceP4wfLlLIrg58UiR6WGUFdXFz/q6uKi1HDaGngofWif1dnJEmCLQoFXgQ+yolE2q7OTKcBbGxvpXL6cJ9rbeRbYO33Yfh7YN7e/X2MjS9rbubGri2JHB5NTjM8CB9Wt+G7hnwoFFpM9Ir2puZltxoyhJfekrdv9bW2MJRt0OL+kK95TwJ4l3ePmRvR0WXweuGvcOHZMT+Z6c2tbGzsCp06YAMDOdXX8KTUw5gPvHTuWU5cu5SBgj9SQeRPwcorlJeDdfVz/GeDo1KD8S1sbp3V09HTN7G7C3Jt7Gvh34Ni6On7T2koL8JVika8UCgSwIfBYezsPA/+TYt2roYHNOjoYV7fqdzavAe9sbOSx9O91b2qY7tfYSCE1sJ4GvlBX19Nw6uvf4ImuLv4pdzyrWGRSr+/YasRPgKMj4jZJ3b0a7iX7Us/MzKwqKpmmff++tnKK95JWOljkj8CUiHgrcCtwWQVlbZgYD/y+o4MNgJ0bG1lXYkEEDwGnNTUBMG3pUl4G/iVNUHEH2QdsgMciOJxsvM5/NTRwPfRMtPBYBNunRszuwH+2tfEq8O6xY1lYLLIUeEdjIwuLRZan1//kuHE8Ajzc1dXTeChtgjzT3MzewNa9jB/K+3NHB1sB144bxxXpqVy3Z4Hdc0/ElnR2chtwXEMDC4tFlkG/jSuA+3NdHLuPd5B4vlhkDLB+fT1zgF1yjZgnyRpiczs6GE/25K83LwAHpNc/rqODL9fX80JzM2fV1/d0bfx7RM+15wEHjB3LPZ2dHMqKsWMvNDfzSHMzd6W6aEj57ywWmdLL6/5vWxvrA3V1dfyto2Olb0ZubWvrOX4O+OjY1XdV+EcETbnj/4ngqF4adVYzpkTEbWm/+77QTmXd383MzAZUJTehX5Ucb0z2WXU+rHb8/Xzo+TIfYEvgxXyGiHgtd/ifZF0Su8vuV1L2znICtqHXBFwVwTnpg/66wDXAXmQNhOVdXVwcwc2NjbwndWO7qrWVM9KTjLnAJ9IH5o81NXF6ocAVy5bxhQkTmAscmc4dWV/P6Z2dbACMq6vj98uXsxnZB/nbli/vaTRs2tDABsDtwCmp7D7AFV1dnNTZycT6et7o7OQR4MOr+aD+UFcX25I13M5tb+ezxSK7F4tMqqtjIbB/asA8VyxyTGsrbwE+2tTEla2tZU21+SQr1j/48/Ll3A3cN348M9vaesZDbUI2ngvgt62tPAncOH48N6aulb1ZUCxSD2zY0MArxSKtwFHjxrGgWOSCzk52TvnmAe9saGB5VxfLgDePGcNWbW3MiGBJZycT6uu5e/lytm1oYJLEArKZFF/q6uLSiF77/P6lo6On6+R2dXW81NnJ8q4ulnZ1cW5nJ6fW19Pa2UkbsF0vT6xK7VhXxzXp/f/H0qUsAD7X1NR/IRvNnpB0cETckks7EHi0WgENN2/e/dNMnTp19RlHuZaWlpqvB9dBxvWQcT0MrrK/+o2IbfIbMBH4DvCLMorfB0yVtI2kRuAYsiEpPSTle/ocTtZLCeAW4CBJ60taHzgopdkwtA7Z18gnp259EyVeBb6aPjyfvnQp20JP4wpgrzFjeqaizD9pAXgHcGP6QP0CKxoxx44f37MoG8B9qftg935+PNS+ZN3n3pm6n/10wgS2BbZfupStCwV2W7qUPYEjUkxvKxT46/JV5xR8mhXjvb4wYQL7Ake2tnJ3ezsdwNtbW5lcKPDu1lbeJnFr6j53X2cnC2GlKdp/tXTpKtd/FngzMLlQ4Pj2dv69oYGtGhq4r7OT7VOe740fzz3ApEKBbxSL/HrMGJrq63lXYyOLU/oVJde+va2tZ2zVxg0N7ANstXQp729tzZ40pvf0CllXvtva2nq63Z3a1MTWwDZLl7JFocBJ7e2sK/HBpia2BaYsXcqhy5YxFtijl6dnD6WncJCNBdsHmLJkCTu3tvI+ic+vsw53tLev0s2vr3+DU5qaqAc2KxT4QWcnv21s7HmKZjXpa8CvJV0GjJf0S+BS4P9VNSozM6tpij6mdS6rsNQAzI+IzcrIexhZf/l6YHpEfEfSNGB2RMyQ9D2yhlUReB34XEQ8mcp+CvhGutR3IuK/8tdetGiRuwwOoobLLqPpzDNXn9HW2JLOTrZaupQ3ypjJ0NZO67nnUjzBqz1Uy8SJE/vvi1shSVuQTXKxNdlwx/8eDTMIDtR9zd9SZ1wProNuroeM6yEzUPVQem9b237q72XFGPl+RcRNwE0laWfl9k8HTu+j7HSyWQvNRqWZufFIZla+iHgBuKDacZiZmXUru4El6XlWnlyiiWzSuM8PdFBmteaDTU18sNpBmI0Akq6gjImOIuL4IQjHzMxsFZU8wTqu5Hgp8HRErLqgj5mZ2eCYU+0AzMzM+lPJOlh/HsxAzMzMVicivlXtGMzMzPpTSRdBd8swM7NhJa3FeCywOdnyH1fn1sYyMzMbcpXMb/wP4EiyWQDnp7JHpPRncpuZmdmgk/RV4GqymWdvBF4DrpT0tTLLby/podxWkHSqpO9LelLSI5Kuk7Reyj9F0rJc/otz19pN0qOS5kj6mZStTyBpA0kzJbWkn+undKV8c9Lr7DrA1WNmZlVSyRisNwHvi4i/dCdIeidwZkQcPOCRmZmZ9e9rwP4R8Vh3QuptMRP44eoKR8RTwC6pXD3ZcnvXAdsDp0dEUdL5ZDPcfj0VeyYidunlchcBJwGzyGbMPQS4GTgNuC0izpN0Wjr+OnAoMDVte6bye1b07s3MbFiq5AnWXmQ3jrx7ydaCNTMzq4bSSS/mUkZ39l4cQNZ4ejYi/hQRxZQ+C3rW6u6VpElAc0TcE9nikpeT9fiArKfHZWn/spL0yyMzC1gvXcfMzEa4Sp5gPQh8V9JZEbFM0njgW8BDgxOamZlZv84BfiXpHLKu61sBZwJnS+r5AjEiylmv8Rjgql7SPwVckzveRtKDQAH4ZurVsUV6/W7zUxrAphGxIMWxQNImKX0LsoWRS8ssKA2gpaWljPD7trblRwvXg+ugm+sh43rIrGk99LdAcSUNrE8AVwKLJL0BrA/MBj6+RlGZmZmtnV+mn8eSPbVSOv54OqeUXt/fRSQ1AodTsti9pDOAIvDrlLQAmBwRr0naDfiDpJ1yr5u3uqdoZZfp7ya+Oi0tLWtVfrRwPbgOurkeMq6HzGDVQyXTtM8D9pa0FdlsTQsi4rkBj8iGneLUqbSee261wzAbEEXfUEaTbQboOocCD0TEwu4ESScA7wcOSN3+iIg2oC3t3y/pGbLxyfNZuRvhlmQzGgIslDQpPb2aBLyc0rufuPVWxszMRrBKnmAhaUNgP2BSRFwgaXOgLiLm91/SRrS996a4997VjsLMbCUR8ewAXepYct0DJR1CNhHDPYnlAAAgAElEQVTFuyOiNZe+MfB6RHRK2pZsgoq5EfG6pMWS9iIbm3w88PNUbAZwAnBe+nl9Lv2Lkq4mm9xiUXdXQjMzG9kqWQfr3cC1ZN0C9wEuILu5/CvwgUGJzszMrA+SJgJfBt4OTMifi4iDyrxGE/Be4LO55F8AY4GZabb1WRFxMvAuYJqkItAJnBwRr6cynwMuBcaTzR54c0o/D/iNpBOB54CjUvpNwGFkk3S0Ap8s602bmdmwV8kTrJ8AR0fEbWkMFmTf1O0x8GGZmZmt1m/JxlddByxbkwukJ1QblqRt10fea8m+aOzt3Gxg517SXyObobA0PYAvrEHIZmY2zFXSwJoSEbel/e6BuO0VXsPMzGyg7AVsGBEd1Q7EzMysWyXrYD0hqXRB4QOBRwcwHjMzs3L9Fdih2kGYmZnlVfL06WvADZJuBMZL+iXZ2KsjyimcBg3/lKw7xyURcV7J+a8CnyabEvcV4FPdA5gldbKiIfdcRBxeQdxmZjY6fQK4SdK9wML8iYiYVpWIzMys5lUyTfssSW8FjgOmky2QuEc5MwhKqgcuJBtIPB+4T9KMiHgil+1BYPeIaJX0ObJJNI5O55ZFxC7lxmpmZjXhO2RTnc8DmnPpq1uDyszMbNCU1cBKDaTbgIMj4oI1eJ09gDkRMTdd72qyJ189DayIuCOXfxZZQ87MzKwvxwBv8vTmZmY2nJQ1BisiOskWdKxkzFbeFmRPvLrNT2l9OZEVU9wCjJM0W9IsSUeuYQxmZja6zAU8wYWZmQ0rlYzB+hZwkaSzyRpIPV0wIqJrNWXVS1qvXTgkHQfsDrw7lzw5Il5MCzveLunRiHimgtjNzGz0uQKYIennrDoG6/bqhGRmZrWukgbWJenn8axoHCnt16+m7HyyfvLdtgReLM0k6UDgDODdEdHWnR4RL6afcyXdSbaopBtYZma1rXsdqe+WpAew7RDHYmZmBpTRwJK0WUS8RNZFcE3dB0yVtA3wAlm/+Y+VvM7bgV8Ch0TEy7n09YHWiGiTtBGwD9kEGGZmVsMiYm3uS2ZmZoOinCdYTwPNuSnTfx8RH6rkRSKiKOmLwC1kT7umR8TjkqYBsyNiBvB9YALwW0mwYjr2HYBfSuoiGwN2Xsnsg2ZmZmZmZsNCOQ2s0vFT+63JC0XETcBNJWln5fYP7KPc34C3rMlrmpnZ6CWpGTiHbMzuRuTuVxExuUphmZlZjStnVkCvJ2JmZsPRvwO7AtOADYAvAc8BP65mUGZmVtvKeYLVIOk9rPhmsPTYszWZmVk1HATsEBGvSeqMiOslzQb+iBtZZmZWJeU0sF4GpueOXys59mxNZmZWDXXAorS/RNJ6wAJgu+qFZGZmtW61DayImDIEcZiZmVXqYbLxV7cBfwEuBJaQTc5kZmZWFeWMwTIzMxuOPgPMS/tfBpYD65Gt12hmZlYVlSw0bGZmNmxExNzc/ivAiVUMx8zMDPATLDMzG2Ek7SZp59zxxpJ+LelhSRdLmlDN+MzMrLa5gWVmZiPNT4DNcseXAG8C/gPYGbigGkGZmZmBuwiamdnIswPZpBakmQMPBXaOiKclzQD+Bny+ivGZmVkN8xMsMzMbaRqA9rS/F/BSRDwNEBHPk010YWZmVhVuYJmZ2UjzOHBU2j8GuLX7hKQtWLE2Vr8kbS/podxWkHSqpA0kzZTUkn6un/JL0s8kzZH0iKRdc9c6IeVvkXRCLn03SY+mMj+TpJTe62uYmdnI5waWmZmNNF8HfinpdeB9wPm5c0cDd5dzkYh4KiJ2iYhdgN2AVuA64DTgtoiYSrbG1mmpyKHA1LSdBFwEWWMJOBvYE9gDODvXYLoo5e0ud0hK7+s1zMxshHMDy8zMRpSI+CswGXgvsG1EPJU7fSPwlTW47AHAMxHxLHAEcFlKvww4Mu0fAVwemVnAepImAQcDMyPi9Yh4A5gJHJLONUfEPRERwOUl1+rtNczMbITzJBdmZjbiRMRi4P5e0p/qJXs5jgGuSvubRsSCdL0FkjZJ6VsAz+fKzE9p/aXP7yW9v9dYRUtLyxq9oYEqP1q4HlwH3VwPGddDZk3rYerUqX2ecwPLzMxqmqRG4HDg9NVl7SUt1iC9Iv3dxFenpaVlrcqPFq4H10E310PG9ZAZrHpwF0EzM6t1hwIPRMTCdLwwde8j/Xw5pc8HtsqV2xJ4cTXpW/aS3t9rmJnZCDdkDSxJh0h6Ks2ktMpgXkljJV2Tzt8raUru3Okp/SlJBw9VzGZmVhOOZUX3QIAZQPdMgCcA1+fSj0+zCe4FLErd/G4BDpK0fprc4iDglnRusaS90uyBx5dcq7fXMDOzEW5IughKqgcuJBuQPB+4T9KMiHgil+1E4I2I2E7SMWSzQh0taUeyvvE7AZsDt0p6U0R0DkXsZmY2eklqIrs3fTaXfB7wG0knAs+xYkr4m4DDgDlkMw5+EiAiXpd0LnBfyjctIl5P+58DLgXGAzenrb/XMDOzEU7ZxEaD/CLSO4BzIuLgdHw6QER8L5fnlpTnHkkNwEvAxqSpa7vz5vN1l120aNHgvwkzM1trEydO7G1ckpXwfc3MbOQovbcNVRfBvmZY6jVPRBTJForcsMyyZmZmZmZmVTdUDaxyZlIa1FmYzMzMzMzMBttQTdPe1wxLveWZn7oITgReL6esu5yYmdlo4vuamdnINVRPsO4DpkraJq03cgzZDEp5+RmVPgLcHtkAsRnAMWmWwW2AqcD/DlHcZmZmZmZmZRuSBlYaU/VFsqls/w78JiIelzRN0uEp26+ADSXNAb7KisktHgd+AzwB/A8wHXhiTaZ7rwVlTIf/VUlPSHpE0m2Stq5GnNWyuvrJ5fuIpJC0+1DGV23l1I+kj6bfocclXTnUMVZbGX9jkyXdIenB9Hd2WDXirAZJ0yW9LOmxPs5L0s9S3T0iadehjrFWlPt/3WjQ2++dpA0kzZTUkn6un9JH7e+gpK3S/z1/T/8/n5LSa6ouJI2T9L+SHk718K2Uvk36XNiSPic2pvRR+7lRUn26F92QjmuxDuZJelTSQ5Jmp7TB/5uIiBGzAfXAM8C2QCPwMLBjSZ7PAxen/WOAa6od9zCrn/cATWn/c66flesn5VsXuAuYBexe7biHU/2QPUF+EFg/HW9S7biHYR39B/C5tL8jMK/acQ9h/bwL2BV4rI/zh5FNUy5gL+Deasc8Grdy/68bLVtvv3fABcBpaf804Py0P2p/B4FJwK5pf13g6fR/UE3VRXo/E9L+GODe9P5+AxyT0i/O/T89aj83kj2wuBK4IR3XYh3MAzYqSRv0v4khW2h4gOwBzImIuRHRDlwNHFGS5wjgsrT/O+AASbXSl3219RMRd0REazqcRTamrVaU8/sDcC7ZH9/yoQxuGCinfj4DXBgRbwBExMtDHGO1lVNHATSn/YmsOt501IqIu8jGzvblCODyyMwC1pM0aWiiqynl/l83KvTxe5f/LHAZcGQufVT+DkbEgoh4IO0vJusxtAU1Vhfp/SxJh2PSFsD+ZJ8LYdV6GHWfGyVtCbwPuCQdixqrg34M+t/ESGtgrc1077Wg0intT2TFope1YLX1I+ntwFYRccNQBjZMlPP78ybgTZLuljRL0iFDFt3wUE4dnQMcJ2k+2cK0Xxqa0EYEL7sxNFzPsGlELICs4QFsktJrom5SF6+3kz29qbm6SF3jHgJeBmaSPdH9R/pcCCu/19H6ufEnwL8BXel4Q2qvDiBrXP9J0v2STkppg/43MVSzCA6UtZnuvRaU/d4lHQfsDrx7UCMaXvqtH0l1wI+BTwxVQMNMOb8/DWTdBPcje/r5F0k7R8Q/Bjm24aKcOjoWuDQifqhskfUrUh119VK21tTy/89DyfXct1FfN5ImANcCp0ZEoZ8HEaO2LiKiE9hF0nrAdcAOvWVLP0ddPUh6P/ByRNwvab/u5F6yjto6yNknIl6UtAkwU9KT/eQdsHoYaU+wKpnuHa083XstKKd+kHQgcAZweES0DVFsw8Hq6mddYGfgTknzyPrfzlDtTHRR7t/X9RHRERH/BzxF1uCqFeXU0Ylk/dyJiHuAccBGQxLd8FfW/1G21lzPsLC7a0/62d2deVTXjaQxZI2rX0fE71NyTdYFQPry706y+/l66XMhrPxeR+Pnxn2Aw9NnmavJugb+hNqqAwAi4sX082WyxvYeDMHfxEhrYK3NdO+1YLX1k7rA/ZKscVVr42f6rZ+IWBQRG0XElIiYQjZG7fCImF2dcIdcOX9ffyCbKAVJG5F1GZw7pFFWVzl19BxwAICkHcgaWK8MaZTD1wzg+DRT017Aou5uGjagyvk9He3ynwVOAK7PpY/K38E0ZuZXwN8j4ke5UzVVF5I2Tk+ukDQeOJBsPNodZJ8LYdV6GFWfGyPi9IjYMn2WOYbsPX2cGqoDAEnrSFq3ex84CHiMofibWNPZMaq1kc3w8TRZf9ozUto0sg/CkH2Y+S0wh2y9rG2rHfMwq59bgYXAQ2mbUe2Yh1P9lOS9kxqaRbDM3x8BPyJbNuFR0mxEtbSVUUc7AneTzdz2EHBQtWMewrq5ClgAdJB9E3gicDJwcu7358JUd4/W2t/XEP9brPJ7Olq3Pn7vNgRuA1rSzw1S3lH7Owi8k6w70yO5e/xhtVYXwFvJZrt9hOzD9Fkpfdv0uXBO+pw4NqWP6s+NZF36u2cRrKk6SO/34bQ9nrtnD/rfhNIFzczMzMzMbC2NtC6CZmZmZmZmw5YbWGZmZmZmZgPEDSwzMzMzM7MB4gaWmZmZmZnZAHEDy8zMzMzMbIC4gWVmZmZmZjZA3MAyMzMzMzMbIG5gmZmZmZmZDRA3sMzMzMzMzAaIG1hmZmZmZmYDxA0sMzMzMzOzAeIGlpmZmZmZ2QBxA8vMzMzMzGyAuIFlZmZmZmY2QNzAMjMzMzMzGyBuYJkNE5LulHRJSdo3Jc2rUkhmZmZrxfc2q0VuYJmZmZmZmQ0QN7DMzMzMzMwGiBtYZmZmZmZmA8QNLLPhowtQSdqYagRiZmY2QHxvs5rjBpbZ8PEysHlJ2q7VCMTMzGyA+N5mNccNLLPh41bgQEkflbSdpNOAfasdlJmZ2Vrwvc1qjiKi2jGYGSBpDPB94GNAI/Br4DXg+IiYUsXQzMzM1ojvbVaL3MAyMzMzMzMbIO4iaGZmZmZmNkCGpIElabqklyU91sd5SfqZpDmSHpG0a+7cCZJa0nbCUMRrZmZmZma2JobqCdalwCH9nD8UmJq2k4CLACRtAJwN7AnsAZwtaf1BjdTMzMzMzGwNNQzFi0TEXZKm9JPlCODyyAaEzZK0nqRJwH7AzIh4HUDSTLKG2lX5wosWLfJAMjOzEWDixIml6+FYL3xfMzMbOUrvbcNlDNYWwPO54/kpra90MzMzMzOzYWe4NLB6+0Yz+kk3MzOrKkmHSHoqjR8+rZfzYyVdk87f292TQ9IUScskPZS2i4c6djMzGzzDpYE1H9gqd7wl8GI/6WZmZlUjqR64kGwM8Y7AsZJ2LMl2IvBGRGwH/Bg4P3fumYjYJW0nD0nQVdLS0lLtEKrOdZBxPWRcD6O/DoZLA2sGcHyaTXAvYFFELABuAQ6StH6a3OKglGZmZlZNewBzImJuRLQDV5ONJ847Args7f8OOECSx6CZmY1yQzLJhaSryCas2EjSfLKZAccARMTFwE3AYcAcoBX4ZDr3uqRzgfvSpaZ1T3hhZmZWRb2NEd6zrzwRUZS0CNgwndtG0oNAAfhmRPylrxcaDd/0job3sLZcBxnXQ8b1MPLrYOrUqX2eG6pZBI9dzfkAvtDHuenA9MGIy8zMbA2VM0a4rzwLgMkR8Zqk3YA/SNopIgq9vVB/N/GRoKWlZcS/h7XlOsi4HjKuh9FfB8Oli6CZmdlIUs4Y4Z48khqAicDrEdEWEa8BRMT9wDPAmwY9YjMzGxJuYJmZmVXuPmCqpG0kNQLHkI0nzpsBnJD2PwLcHhEhaeM0SQaStgWmAnOHKG4zMxtkQ9JF0MzMbDRJY6q+SDbxUj0wPSIelzQNmB0RM4BfAVdImgO8TtYIA3gXME1SEegETvb4YjOz0cMNLDMzG5UkjQH2AjaPiGskrQMQEUsH4voRcRPZJE35tLNy+8uBo3opdy1w7UDEYGZmw4+7CJqZ2agj6S3A08B/kj1JAng3njTJzMwGmRtYZmY2Gl0EnBURbwY6UtqfgXdWLyQzM6sFbmCZmdlotBPw32k/oKdr4PiqRWRmZjXBDSwzMxuN5gG75RMk7UG2oL2Zmdmg8SQXZmY2Gp0J3CjpYqBR0unAycBnqhuWmZmNdn6CZWZmo05E3AAcCmxMNvZqa+BDEfGnqgZmZmajnp9gmZnZqBQRDwCfr3YcZmZWW9zAMjOzUSct+Nur/FpVZmZmA80NLDMzG422KjnejGwdrOuqEIuZmdUQN7DMzGzUiYhPlqZJOgQ4tgrhmJlZDfEkF2ZmViv+BBxZ7SDMzGx0K/sJlqQvA1dGxKtr8kLpm8OfAvXAJRFxXsn5HwPvSYdNwCYRsV461wk8ms49FxGHr0kMZmZWGyRtW5LUBHwMeL4K4ZiZWQ2ppIvggcB3Jd0JXAH8ISLayikoqR64EHgvMB+4T9KMiHiiO09EfCWX/0vA23OXWBYRu1QQq5mZ1bY5QABKx63Ag8AJVYvIzMxqQtldBNNTo62Bm4FTgZckXSLpXWUU3wOYExFzI6IduBo4op/8xwJXlRubmZlZXkTURUR9+lkXERMiYt+IuL/asZmZ2ehW0RisiHgtIi6MiHeQzcb0z8AdkuZJOkPShD6KbsHK3TLmp7RVSNoa2Aa4PZc8TtJsSbMkuf+8mZmZmZkNSxXPIijpAOA4sidQs4ELgOeAU8iebu3bW7Fe0qKPlzgG+F1EdObSJkfEi6lP/e2SHo2IZyqN3czMRi9Jz9P3vaVHREwegnDMzKxGVTLJxQ/IGj+LgMuBb0bEC7nzs4A3+ig+n5XXJNkSeLGPvMcAX8gnRMSL6efcNAbs7YAbWGZmlndctQMwMzOr5AnWOOCDEXFfbycjokPS7n2UvQ+YKmkb4AWyRtTHSjNJ2h5YH7gnl7Y+0BoRbZI2AvYhe2pmZmbWIyL+XO0YzMzMKmlgfY9sFqYeqfEzPveE6cneCkZEUdIXgVvIpmmfHhGPS5oGzI6IGSnrscDVEZHv4rED8EtJXWRjxs7Lzz5oZmbWG0m7kHVb34hcV/WIOKtqQZmZ2ahXSQPrD8CnWLkb4JbAJcCeqyscETcBN5WknVVyfE4v5f4GvKWCOM3MrMZJOgn4MdniwoeSjRE+CLi+mnGZmdnoV8ksgttHxKP5hHT85oENyczMbK39G3BIRHyQbC3FDwIfATqqG5aZmY12lTSwXpa0XT4hHb82sCGZmZmttU0i4i9pv0tSXUTcDHygmkGZmdnoV0kDazpwraT3S9pR0geA35F1ETQzMxtO5kuakvafBo6QtC/QXrWIzMysJlQyBus8sq4VPyCbcv15ssbVjwYhLjMzs7VxAdkkSfOAaWRfCDYCX65iTGZmVgPKbmBFRBfw/bRZLfnb32hoaal2FGYjUnHqVNh772qHUXMi4tLc/s1p1tvGiFhSvajMzKwWVPIEq3udqrcBE/LpETF9IIOy4aWhpYWmM8+sdhhmI1LruedSdANryEn6CfDr7rUbI6Iddw80M7MhUHYDS9I3gLOAh1l5PawgG59lZmY2XAi4XtJS4Ergyoh4qsoxmZlZDajkCdapwB4R8chgBWNmZjYQIuIUSV8BDiBbxH6WpLlkT7U8dtjMzAZNJbMILgOeHKxAzMzMBlJEdEXEzIj4FLAz2bIiHkdsZmaDqpIG1pnAzyVNklSX3wYrODMzszUlaYKk4yTdSDZVexE4ocphmZnZKFdJF8FL089P59JENgarfqACMjMzW1uSfgscCjwAXAWcEBGvVjcqMzOrBZU0sLYZtCjMzMwG1mzgaxHxXLUDMTOz2lLJOljPAqQugZtGxIJBi8rMzGwtRMT51Y7BzMxqU9njpyStJ+lKYDkwJ6UdLunbgxWcmZmZmZnZSFLJBBUXA4uArVmxWOM9wNEDHZSNTMWuLpoKhZ7j61tbqSsUuGP5cgAWFousVyiwsFjsyfPWQoGpuTIAexYKbFwosHnaft/autK5iYUCXV1dADzZ0cEGqfzty5axaW6/rlDgS4sX95Q9oFDg39JxsauLzxUKTC4U2KpQYItCgaNK4ujNKYsXs3/K19rZyd6FAvsUCrR2djKxUGDLtG1VKHBK7rVvX7aMcbnzWxYK/L/c+dW5a/lytltNfA+0t3NOP9d8X6HAZoUCB/VznSMLBa5N9b1BoUChs7PsGE9fvJgP9HLtV4pFdkx1/PUK3nNf2ru6+EChwCbp32/msmW95juwUGCDQoHmft7vmYsXMznFtnXud7XU9unfbKNCgf0LBYrp98/MzMysVCUNrAOAL6eugQEQEa8Am5RTWNIhkp6SNEfSab2c/4SkVyQ9lLZP586dIKklbZ4BaphqqKsjoKfx87NikWbg1XT83WXL2A3YtCHrmXrZ0qUsBV4quc5zwFVjxvBiczMHAd/JNcieA9YFrk0fhG9vb2fLdO6vnZ1sldvfHPhzRE/ZecCe9dl8LIctWcL9wL1NTTzf3Mw9TU30/tF6ZU9EsIPEC8Uiuy5dymbA3c3NvNjVhYD5zc3Mb27mlnHj+F0E31uypCeeXXPn5zc38/111y3jFTN3FYtsu5o817S18WDu/ebdvXw5TwAvNTfzp+bmPq/xh+ZmPtzUxNyODsYBzfXlz1/zSARvkVZJP3/ZMqYCLzQ3c34Z77l9NY2Xk5YsYRHw0oQJfKaujm90dPSer6GBM+vq2KKfa02pq+O+piZeaG7mAOD09vZV8hS7uri1qYn5zc08t846PAzc1ta22vdho18Z97Wxkq5J5++VNCV37vSU/pSkg4cybjMzG1yVNLAWARvlEyRNBlY7FktSPXAh2YxOOwLHStqxl6zXRMQuabskld0AOBvYE9gDOFvS+hXEbUOoEXitq4uFxSJPkP2DvZE+9P8hglPHjOnJ++3OTqbV1zMWeDx9sF3S2clrwN6NjQBsW1dHd4klnZ0UgKMkrk2Nrgc6O9kunX+0q4upuf2DyBpvb6SnMC8BB44dyzWtrcwGbl9nHSalxt7khgb+2E/Do9s8sj+CvVpbOVji96nMbe3tbJ7Lt2NjIx+RuDk1Fh7t6mKH1V4d1isUeGd6MrN5ocA9qSH5cFcXO6XGy6+XLmVyocCkQoFd01Omy5Yu5RcR3A1smSsH2dOz97e383o690qxyFvSa2yce8q3sFhkk/S0Z2Z7O1un8k93dLBXocD/b+/e4+OqyzyOf55JmiZpm1JApNCWtlIU7ApIt5ablFu5uHJRkLogVaq7Iqsi6C6Iys0iVEXQBQWx2u1aQMBLsLi8SqEIuIXCUq5SUsqtpVyE2mkzuc3Ms3+cX9LpMGlmksnMJPm+X695zTlnzjnznF8mc84zv8uZGGp6FjU3d60/I9SKHRqPswY4tHrbbp0Xbd7M9e78BZgRj7M+meSgUPs0Lh7n1lBbdvWWLUyNxzkgHucDISnNpTWd5o/Aovp6YrEYR1VXs76bdT9VX89z6fR2R+eZO2JEV8JfZ0au+rrqWIzxYZ372tupBaaHz6dUNjPbycw+Y2b/HuZ3M7NxPW2X577zOa/NBTa6+57Aj4Crwrb7ALOBDwLHAteH/YmIyCBQSIJ1E3CHmR0OxMzsQGAhUdPBnkwH1rj7WndvB24BTszzfY8Blrr7O+6+EVhKdEKSClQLvJlOc0VLCzOBBqIE63eJBFXAx+vqALhmyxbqgdNHjGA34M+hFmJZWxs7APVVVbydTHJnOs1FISlb1tbGe4HPDR/OivB+q6Gr1mQNsG8s1jX94ViMDwK/amnhkbY2RgKjq6q4KZnkY8DIbmpnXuzoYEpGM8RM64EfuvOVWIxrM2pjHk2l3nUhXwd0NohbAzRCV/PAw3M0W3uuo4M4cEosxisNDcwEvhcSzybgwKoqnuvo4LxUituHD2dDQwNjgIsTCeaMGMFEYNGwYaxraODA2tqu/R5RV8cRwFlmrGtoYFQsxpL6el5paODZ+noWhwR4WUZt4KOpFB8I059qaeH86mpeamhgwbBhXB4S1tMTCfYlqhU70owXgcOHD9/mmOaNGsVY4I6aGlY0NPCJRIJ/NOPNhga+X13NN0OivCqdZiPw2/p61jY0dPs3WNjSwq7AHiHhecedEe8qya1WQ1diuj2PtbWx2J2zu/lMfHnzZnaOx/l0RweLa2oYU0DNnpSHmR1G9BE4neg+jgBTgJ8W6S3yOa+dSHSeBLgdONLMLCy/xd3b3P1Foq+I6UWKS0REyqyQBOsq4DdEv9gNAxYAfwCuzWPb3YFXM+bXhWXZPmlmT5rZ7WbW2dor322lAtQSNQlsdOei2lpGApvc+UkyyanhQjedTnN1Os0VIXHaA1gVLtpXpFJsAXaPx9ktkeB9bE3KVqRSTAQ+VFNDCni2vZ2XgYPCxfarbK1BeRWYWVPDMbEYS9Jp7u/oYEKI8WXgQ7GtH/33hRqV40PSM2nYMJoaGojFtv33eKytjeFEH751WU3xVgP7ZF3Ir3XvarL4KnB/bW1X88D7ctSW3dPWxj7AuSNHAjA1FqMzDVsHHD18ON9raWEWMD0kMnsBb4ZYXgcO66Zm5QWiBA3ggbY2Tk4k2D0eZ79Egs4U5uGM2sC/AtNiMX6TSNAEfC2ZZFw8zlmh6eDT7e08AVwXYp1RXc2uQG3s3V8pbwOH1NTwdPh7dSamM2tquo7veeCcWKwrcerub/BsOs37MuZXJJOMzXnEkZeB6T0kQw+2tnJcWxtfNeOsEbnTtZ+MGsXfGho4GrgkRzNCqUjXAKe5+7FENxgGeJjiJTL5nJu61ofKHUEAABZLSURBVHH3JFFLkJ3y3FZERAaoQoZpd6IT1jW9eJ9cPyFndxa5E7jZ3dvM7ItEv/odkee2UiHqgN92dLAjMLWmhlFtbWxwZxVwR309AJc1N/Mm8JmODujooB34x7D90+6cANzS0MDiRIK5ySTxVIqGqiqeduf9IYmZBvy8rY2/AYcNH84bySTNwIE1NbyRTNIa3v89sRhXJxLsmk53JQ/ZKcgLDQ2cFI8ztoeajvs7OhgPLK6t5aOtrcxMJDg5HNPLwDkZF/JbUimWAT+truaNZJIWomaD2/NYRhPHzvm9zXg1mWQYMKaqijXAJzKSjueAWbEYazs6qCOq+ctlPXBkeP8zOjqYX1XFnBEjuKG5mR+H5Pav7swM+34JOHL4cG5obeU44PashPD65mbGEzWfA1ieTDIxx/s+0tbGGCAWi/GXjo5triDvaWvrmn8F+FRW7Vcuf3enPmP+f9w5M0dSB1Ei/wZwxHbKfUlLC3M6OrgwFuP8kCxuz76xGL/TABcDxUR3XxamO88Z7RR2/8ftyefc1N06BZ3XmpqaCgirMg2GY+grlUFE5RBROQz8MpgyZUq3r+V9ojGzI7p7zd3v7WHzddD1Yz7AOOC1rH28nTH7c0Jb9bDtzKxtl/fwflIm9cDN7lwSLvRHAbcCM4gShNZ0mp+586eaGg4PzdhuTiS4KDQVWwt8Nlww/3N9PRfG4yxqaeGckSNZC5wUXjupqooLUyl2JKo1+W1rK7sSXcgva23t6g/13upqdgTuBb4atj0YWJRO8y+pFKOrqtiYSvEk8MluLtQ7rUqnmUyUuF3e3s6/JpNMSyYZG4ttcyH/SjLJ7ESCfyDqB7Q4kchrJJjn2Do85/2trTwErKyrY2lbW1d/qF2I+nMB3JZI8BywpK6OJaFpZS4bkkmqgJ2qq3krmSQBnFpby4ZkkvmpFFPDei8Bh1RX05pO0wJ8YNgwxre10ejOllSKkVVVPNTayuTqasaasYFoJMXX02l+5Z6zze8DHR1dTSf3jMV4PZWiNZ2mOZ3m8lSKc6uqSKRStAF7ZvTP684+sRi3huO/sbmZDcDZIcn9+ubNfHLYsK7mkQ+3tzMyHHenh1pb+X1HB98fNYpfNzfz5VSKq6qq+EJWzVXnvupiMXaMxZhQXc3LySQL02k+18PnRCrGs2Z2jLvfnbHsKOCpIu2/x/NaxjrrzKwaGA28k+e2XbZ3Eh8ImpqaBvwx9JXKIKJyiKgcBn8ZFHKl8IusRyPwP0R9s3qyEphiZpPMrIaoc29j5gpmltnS5wSiVkoAdwOzzGxMGNxiVlgmFWgE0c+wXwzN+kab8TfgvHDxfGFzM5OhK7kCmDFsWNdIKZk1LQAHAkvCBfV6tiYxn66r67pnAMDK0HywczqzP9ShRM3nDgkX2teOHMlk4P3NzewRj3NAczMfAU4MMe0bj/NgjuG6n2drf69zRo7kUOCkRIKH2tvpAPZPJJgQj3NYIsG+ZtwTakRWplK8AdsM0f6LMFBEppeBDwAT4nHObG/n+upqxldXszKV4v1hne/V1fG/wNh4nG8mk/x62DDqq6r4aE0Nm8PyRVn7vretratv1XuqqzkYGN/czD8lElFNYzimt4ia8i1ra+tqdndufT17AJOam9k9Hudf2tsZZcbJ9fVMBiY2N3NcSwvDyd0Ub1WohYOoL9jBwMQtW5iaSPAxM740YgT3tbe/q5lfd3+Dr9bXUwXsGo/zg1SK22pqqI7FiKdS/NKd94QYdovHmdXezkaigUMWhjL5eXs7L4QmleemUrQDl6ZSjIvH2S80Ec3c19L2dg4IzSkPSSQ4yYxv5VHTJRXhfODXZrYQqDOzG4BfAd8o0v57PK+F+c6Rb08B7g2tQRqB2WGUwUlEfcMeKVJcIiJSZubdDOvc44bRiEffAja7+9V5rH88UfPCKmCBu88zs8uAR9290cy+R5RYJYl+4Tvb3Z8L254FfDPsap67/zJz35s2bVKTwX5UvXAh9d/+ds8rSq9tSaUY39zMxjxGMpR3+8bmzcRgu8PAJ9Np9t+yhQdHjGD0dvpl5bOvQiQuv5zkHN1dotPo0aN7HnWkSMxsd6JBLvYg6vP03+6+roj77+m8VgssAvYnOq/Ndve1YduLgLOIznnnuvufMvc9mM5rg/2X6nyoDCIqh4jKYfCVQfa5rddt0d09ZWbziJo69JhguftdwF1Zy76TMX0hcGE32y4gGlRDZFBamtEfSQqXzz3FqmMxnsojgS3k/mRS2dx9PTC/H/ff03mtFTi1m23nAfP6KzYRESmfvnb2PRpQj2+RPjq5vp6Tyx2EyABnZovIYxAkdz+zBOGIiMgQVcggF6+y7YmrnmhU7i8VOygREZFeWFPuAERERAqpwToja74ZeN7d333HVBERkRJz90vLHYOIiEgh98G6vz8DERERKaZwe5FPA7sRDYN+S8a9sURERPpFIU0E1bZdREQGBDM7D7gA+CXwODABWGxm8939h2UNTkREBrVCmgj+neh+HncS3bJnAvBxYCHw9na2ExERKbXzgSPc/enOBeGHwqWAEiwREek3hSRYewEfc/cHOheY2SHAt939mKJHJiIi0jfZg16sJY+WGCIiIn0RK2DdGcCKrGUPAwcWLxwREZGiuAT4hZlNMbM6M9sLuBG42MxinY/yhigiIoNRISeXx4ErzKwOIDzPA1b1R2AiIiJ9cAPRABergS3Ac8DpRElWB5AMzyIiIkVVSBPBzwKLgU1mthEYAzxKdMISERGpJJPKHYCIiAxNhQzT/hJwkJmNJxrydoO7v9JfgUnlSE6ZQuLyy8sdhsiAlJwypdwhDEnu/nK5YxARkaGpkBoszGwnYCYw1t3nm9luQMzd1/VHcFIhDjqI5EEHlTsKEZG8mdlo4CvA/sDIzNfcfVZZghIRkSGhkPtgHQbcQdQs8GBgPjAF+DrRcO0iIiKV4jagCvgd0FLmWEREZAgppAbrGuA0d18W+mBBNIrg9OKHJSIi0iczgJ3cXQNZiIhISRUyiuBEd18WpjvvI9JOgc0MRURESuBBYO9yByEiIkNPIcnRs2Z2jLvfnbHsKOCpIsckIiLSV58F7jKzh4E3Ml9w98vKEpGIiAwJhSRY5wN/NLMlQJ2Z3UDU9+rEfDY2s2OBa4naxN/k7ldmvX4e8Hmie5O8BZzVOQqUmaXYmsi94u4nFBC3iIgMPfOA8cBLQEPGcs+5toiISJEUMkz7CjP7EHAGsAB4FZiezwiCZlYFXAccDawDVppZo7s/m7Ha48A0d0+Y2dlEg2icFl5rcff98o1VRESGvNnAXu6+odyBiIjI0JJXghUSpGXAMe4+vxfvMx1Y4+5rw/5uIar56kqw3P2+jPVXECVyIiIivbEW0AAXIiJScnklWO6eMrNJFDYoRqbdiWq8Oq0DPrKd9ecCf8qYrzWzR4maD17p7r/vZRwiIjI0LAIazewnvLsP1r3lCUlERIaCQvpgXQr81MwuJkqQutqxu3u6h20tx7Kc7eDN7AxgGnBYxuIJ7v6amU0G7jWzp9z9hQJiFxGRoeWc8HxF1nIHJpc4FhERGUIKSbBuCs9nsjU5sjBd1cO264g6G3caB7yWvZKZHQVcBBzm7m2dy939tfC81syWA/sDSrBERCQnd59U7hhERGRo6jHBMrNd3f11oC8nq5XAlNDMcD1R5+N/znqf/YEbgGPd/c2M5WOAhLu3mdnOwMFEA2CIiIiIiIhUlHxqsJ4HGjKGTP+tu3+ikDdx96SZ/RtwN1Ft1wJ3f8bMLgMedfdG4PvASOA2M4Otw7HvDdxgZmmiPmBXZo0+KCIisg0zawAuIWpuvjMZTdXdfUKZwhIRkSEgnwQru//UzN68kbvfBdyVtew7GdNHdbPdX4B/6M17iojIkHU9UXP0y4D/JhqZ9hvAHeUMSkREBr98EizdlFFERAaaWcDe7v62maXc/Q9hNNo7gR+VOTYRERnE8kmwqs3scLbWZGXPa8hbERGpNDFgU5jeYmY7ABuAPcsXkoiIDAX5JFhvAgsy5t/OmteQtyIiUmmeIOp/tQx4ALgO2ELUr1hERKTf9JhgufvEEsQhIiJSTF9ga0uLrwBXAjsQ3WpERESk3xRyHywREZEBwd3XZky/BcwtYzgiIjKExModgIiISLGY2QFmNjVj/j1m9msze8LMfmZmI8sZn4iIDH5KsEREZDC5Btg1Y/4mYC/gRmAqulG9iIj0MzURFBGRwWRvokEtCCMHHgdMdffnzawR+AvwpTLGJyIig5xqsEREZDCpBtrD9AzgdXd/HsDdXyUa6EJERKTfKMESEZHB5Bng1DA9G7in8wUz252t98YSERHpF2oiKCIig8l/AHea2c+AFHBIxmunAQ+VJSoRERkyVIMlIiKDhrs/CEwAjgYmu/vqjJeXAF/r63uY2Y5mttTMmsLzmG7WmxPWaTKzORnLl5vZajNbFR679DUmERGpHEqwRERkUHH3ze7+mLtvzlq+2t1fK8JbXAAsc/cpwLIwvw0z2xG4GPgIMB24OCsRO93d9wuPN4sQk4iIVAglWCIiIoU5EVgYphcCJ+VY5xhgqbu/4+4bgaXAsSWKT0REysjcvdwx9NmmTZsG/kGIiAwBo0ePtnLH0Fdm9nd33yFjfqO7j8la5+tArbt/N8x/G2hx9x+Y2XJgJ6I+YncA3/Wsk3Hmea2pqanfjkVERHpnypQpXdPZ57aSDXJhZscC1wJVwE3ufmXW68OB/wIOAN4GTnP3l8JrFwJziU5GX3H3u0sVt4iIDD1mdg/b3rC400X57iLHss6k6XR3X29mo4gSrM8Qnf9yyjyJD0RNTU0D/hj6SmUQUTlEVA6DvwxKkmCZWRVwHVGn43XASjNrdPdnM1abC2x09z3NbDZwFXCame1DNNTuB4HdgHvMbC93T5UidhERGXrc/ajuXjOzN8xsrLtvMLOxQK4+VOuAmRnz44DlYd/rw/NmM1tM1Eer2wRLREQGllL1wZoOrHH3te7eDtxC1IY9U2ab9tuBI83MwvJb3L3N3V8E1oT9iYiIlEMj0Dkq4BzgDznWuRuYZWZjwuAWs4C7zazazHYGMLNhwD8BT5cg5rIZzL9S50tlEFE5RFQOg78MSpVg7Q68mjG/LizLuY67J4luBrlTntuKiIiUypXA0WbWRNQy40oAM5tmZjcBuPs7wOXAyvC4LCwbTpRoPQmsAtYDPy/9IYiISH8pVR+s7bVF72mdfLYVEREpCXd/Gzgyx/JHgc9nzC8AFmSt00zU11hERAapUiVY64DxGfPjgOx7kXSus87MqoHRwDv5bDsYRqUSERHppPOaiMjAVaomgiuBKWY2ycxqiAataMxaJ7NN+ynAvWHY2kZgtpkNN7NJwBTgkRLFLSIiIiIikreS1GC5e9LM/o2o028VsMDdnzGzy4BH3b0R+AWwyMzWENVczQ7bPmNmvwGeBZLAORpBUEREREREKlGparBw97vcfS93f5+7zwvLvhOSK9y91d1Pdfc93X26u6/N2HZe2O79gJvZajNbY2YXZL9PqOm6Nbz+sJlNLNEhVgQzO7aH8jnPzJ41syfNbJmZ7VGOOMulp/LJWO8UM3Mzm1bK+Motn/Ixs0+Fz9AzYYjpISWP/7EJZnafmT0e/s+OL0ec5WBmC8zsTTPLOSqeRX4cyu5JM/twqWOU/mFmO5rZUjNrCs9jullvTlinyczm5Hi9sbvPT6XrSxmYWb2ZLTGz58J365W5tq1keXw3dnt9ZmYXhuWrzeyYUsZdTL0tAzM72sweM7OnwvMRpY69mPryWQivTzCzLRbdsH1gcvcB8yCq/XoBmAzUAE8A+2St8yXgZ2F6NnBrueOusPI5HKgP02erfLYtn7DeKODPwApgWrnjrqTyIWqi+zgwJszvUu64K7CMbgTODtP7AC+VO+4Sls9HgQ8DT3fz+vHAn4gGL5oBPFzumPUo2t9+PnBBmL4AuCrHOjsCa8PzmDA9JuP1TwCLu/v8VPqjL2UA1AOHh3VqgAeA48p9TAUce6+vz8L35BNEI2xOCvupKvcxlbgM9gd2C9NTgfXlPp5ylEPG63cAtwFfL/fx9PZRshqsIunL/bSGgh7Lx93vc/dEmF1BNGjIUJHP5weioZXnA62lDK4C5FM+XwCuc/eNAO6e6warg1k+ZeRAQ5gezbsH9Bm03P3PRE28u3Mi8F8eWQHsYNGNemXgyzz3LgROyrHOMcBSd38nfIcsBY4FMLORwHnAd0sQa3/pdRm4e8Ld7wMI3y3/x8A6P+t+p30oA3d/3N07zxXPALVmNrwkURdfn67Vzewkoh8enilRvP1ioCVYfbmf1lBQ6D3D5hL9mjxU9Fg+ZrY/MN7d/1jKwCpEPp+fvYC9zOwhM1thZseWLLrKkE8ZXQKcYWbrgLuAL5cmtAFB9zUcvN7r7hsAwvMuOdbZ3t//cuCHQCJ7owGkr2UAgJntAHwcWNZPcfYH3e+0eNeonwQed/e2foqzv/W6HMxsBPAfwKUliLNflWqY9mLpy/20hoK8j93MzgCmAYf1a0SVZbvlY2Yx4EfAZ0sVUIXJ5/NTTdRMcCbRr6sPmNlUd/97P8dWKfIpo08Dv3L3H5rZgUSD90x193T/h1fxhvL384BnZvcAu+Z46aJ8d5FjmZvZfsCe7v617L4Ylaa/yiBj/9XAzcCPPaMv+gCg+50W4RrVzD4IXAXMKmJcpdaXcrgU+JG7bxnojc8GWoLVl/tpDQX5lA9mdhTRyeCwAfwLSW/0VD6jiNo+Lw//2LsCjWZ2gkc3EB3s8v3/WuHuHcCLZraaKOFaWZoQyy6fMppLaPbk7v9rZrXAzsBQa06ZS17fUVKZ3P2o7l4zszfMbKy7bwjNPnN93tcR/TjTaRywHDgQOMDMXiK6LtnFzJa7+0wqTD+WQacbgSZ3v6YI4ZZSv97vdIDo0zWqmY0Dfgec6e4v9H+4/aYv5fAR4BQzmw/sAKTNrNXd/7P/wy6ugdZEsC/30xoKeiyf0ATuBuCEIdh/Zrvl4+6b3H1nd5/o7hOJ+qgNleQK8vv/+j3RQCmY2c5ETQYH0q+sfZVPGb0CHAlgZnsDtcBbJY2ycjUCZ1pkBrCps0mVDHiZ5945wB9yrHM3MMvMxoQR9mYBd7v7T919t/C9ewjwfCUmV3nodRkAmNl3iS40zy1BrMWm+532oQxCs9AlwIXu/lDJIu4fvS4Hdz804xrsGuCKgZhcAQNrFMGQJx0PPE80QslFYdllRBfCEF3M3EbUSfIRYHK5Y66w8rkHeANYFR6N5Y65ksona93lDKFRBPP8/BhwNdF96Z4CZpc75goso32Ah4hGTloFzCp3zCUsm5uBDUAH0S+Uc4EvAl/M+PxcF8ruqaH2/zWYH0T9SJYBTeF5x7B8GnBTxnpnhfPzGuBzOfYzkYE7imCvy4DoV34H/ppxfv58uY+pwOPv9fUZUauaF4DVDKDRE4tVBsC3gOaMv/0qBvAovX35LGTs4xIG8CiCFg5CRERERERE+migNREUERERERGpWEqwREREREREikQJloiIiIiISJEowRIRERERESkSJVgiIiIiIiJFogRLRERERESkSJRgiYiIiIiIFMn/Axq+YJq6/QoAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.traceplot(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>mc_error</th>\n",
       "      <th>hpd_2.5</th>\n",
       "      <th>hpd_97.5</th>\n",
       "      <th>n_eff</th>\n",
       "      <th>Rhat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>changepoint</th>\n",
       "      <td>7.355990e+05</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>7.355990e+05</td>\n",
       "      <td>7.355990e+05</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>u__0</th>\n",
       "      <td>inf</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>inf</td>\n",
       "      <td>inf</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>u__1</th>\n",
       "      <td>inf</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>inf</td>\n",
       "      <td>inf</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>u__2</th>\n",
       "      <td>inf</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>inf</td>\n",
       "      <td>inf</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>u__3</th>\n",
       "      <td>inf</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>inf</td>\n",
       "      <td>inf</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     mean   sd  mc_error       hpd_2.5      hpd_97.5  n_eff  \\\n",
       "changepoint  7.355990e+05  0.0       0.0  7.355990e+05  7.355990e+05    NaN   \n",
       "u__0                  inf  NaN       NaN           inf           inf    NaN   \n",
       "u__1                  inf  NaN       NaN           inf           inf    NaN   \n",
       "u__2                  inf  NaN       NaN           inf           inf    NaN   \n",
       "u__3                  inf  NaN       NaN           inf           inf    NaN   \n",
       "\n",
       "             Rhat  \n",
       "changepoint   NaN  \n",
       "u__0          NaN  \n",
       "u__1          NaN  \n",
       "u__2          NaN  \n",
       "u__3          NaN  "
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pm.summary(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "supplied range of [inf, inf] is not finite",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-44-7613a60d8120>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mpm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot_posterior\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtrace\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\pymc3\\plots\\posteriorplot.py\u001b[0m in \u001b[0;36mplot_posterior\u001b[1;34m(trace, varnames, transform, figsize, text_size, alpha_level, round_to, point_estimate, rope, ref_val, kde_plot, plot_transformed, bw, ax, **kwargs)\u001b[0m\n\u001b[0;32m    114\u001b[0m                               \u001b[0mpoint_estimate\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mpoint_estimate\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mround_to\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mround_to\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    115\u001b[0m                               \u001b[0malpha_level\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0malpha_level\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mref_val\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mref_val\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0midx\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 116\u001b[1;33m                               rope=rope[idx], text_size=scale_text(figsize, text_size), **kwargs)\n\u001b[0m\u001b[0;32m    117\u001b[0m             \u001b[0ma\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_title\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfontsize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mscale_text\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtext_size\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    118\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\pymc3\\plots\\artists.py\u001b[0m in \u001b[0;36mplot_posterior_op\u001b[1;34m(trace_values, ax, bw, kde_plot, point_estimate, round_to, alpha_level, ref_val, rope, text_size, **kwargs)\u001b[0m\n\u001b[0;32m    136\u001b[0m         \u001b[0mset_key_if_doesnt_exist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'align'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'right'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    137\u001b[0m         \u001b[0mset_key_if_doesnt_exist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'color'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'#87ceeb'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 138\u001b[1;33m         \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mhist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtrace_values\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    139\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    140\u001b[0m     \u001b[0mplot_height\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_ylim\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\matplotlib\\__init__.py\u001b[0m in \u001b[0;36minner\u001b[1;34m(ax, data, *args, **kwargs)\u001b[0m\n\u001b[0;32m   1808\u001b[0m                         \u001b[1;34m\"the Matplotlib list!)\"\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mlabel_namer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1809\u001b[0m                         RuntimeWarning, stacklevel=2)\n\u001b[1;32m-> 1810\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0max\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1811\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1812\u001b[0m         inner.__doc__ = _add_data_doc(inner.__doc__,\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\matplotlib\\axes\\_axes.py\u001b[0m in \u001b[0;36mhist\u001b[1;34m(self, x, bins, range, density, weights, cumulative, bottom, histtype, align, orientation, rwidth, log, color, label, stacked, normed, **kwargs)\u001b[0m\n\u001b[0;32m   6589\u001b[0m             \u001b[1;31m# this will automatically overwrite bins,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   6590\u001b[0m             \u001b[1;31m# so that each histogram uses the same bins\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 6591\u001b[1;33m             \u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbins\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mhistogram\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbins\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mweights\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mw\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mhist_kwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   6592\u001b[0m             \u001b[0mm\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfloat\u001b[0m\u001b[1;33m)\u001b[0m  \u001b[1;31m# causes problems later if it's an int\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   6593\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mmlast\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\numpy\\lib\\histograms.py\u001b[0m in \u001b[0;36mhistogram\u001b[1;34m(a, bins, range, normed, weights, density)\u001b[0m\n\u001b[0;32m    708\u001b[0m     \u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mweights\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_ravel_and_check_weights\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mweights\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    709\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 710\u001b[1;33m     \u001b[0mbin_edges\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0muniform_bins\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_get_bin_edges\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbins\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mweights\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    711\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    712\u001b[0m     \u001b[1;31m# Histogram is an integer or a float array depending on the weights.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\numpy\\lib\\histograms.py\u001b[0m in \u001b[0;36m_get_bin_edges\u001b[1;34m(a, bins, range, weights)\u001b[0m\n\u001b[0;32m    361\u001b[0m             \u001b[1;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'`bins` must be positive, when an integer'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    362\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 363\u001b[1;33m         \u001b[0mfirst_edge\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlast_edge\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_get_outer_edges\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    364\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    365\u001b[0m     \u001b[1;32melif\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mndim\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbins\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\numpy\\lib\\histograms.py\u001b[0m in \u001b[0;36m_get_outer_edges\u001b[1;34m(a, range)\u001b[0m\n\u001b[0;32m    251\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0misfinite\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfirst_edge\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0misfinite\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlast_edge\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    252\u001b[0m             raise ValueError(\n\u001b[1;32m--> 253\u001b[1;33m                 \"supplied range of [{}, {}] is not finite\".format(first_edge, last_edge))\n\u001b[0m\u001b[0;32m    254\u001b[0m     \u001b[1;32melif\u001b[0m \u001b[0ma\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msize\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    255\u001b[0m         \u001b[1;31m# handle empty arrays. Can't determine range, so use 0-1.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: supplied range of [inf, inf] is not finite"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxkAAAE5CAYAAAD80O3gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XucVmW5+P/PBQOKB8DwrBDWHr+ipOlGhRABrRQzqSCTr4c0s8zUyrR0t78eMw/Uz3Knpak/D78dHvCroo1SuT3hKY+p6LZBVCR0a6BUYCFw/f54HsdhZmAeYM3zzAyf9+vFi2etda+1rnXN4Znrue91r8hMJEmSJKkoPWodgCRJkqTuxSJDkiRJUqEsMiRJkiQVyiJDkiRJUqEsMiRJkiQVyiJDkiRJUqEsMtSlRMTVETGx1nFUKiIaIqJ/O22OjIitqxWTJElSR7PIkDpQZh6Qme+00+xIwCJDkiR1GxYZ6tQi4oiIeCYi/hgR15VX7x0RD0XE7Pd7NSJio4i4OyKejIhnI2J8ef3giHghIn4VETMj4rcR0ae8bffysR+OiMkR8Vx5fc/y8mPl7V8vrx8TEfdHxC0R8XxE/DIiepS3TSqf97mIuKBZ/K9ExKYri6Mc/zDgPyPi6fdjkyRJ6sosMtRpRcROwA+AfTJzF+Bb5U1bAXsBBwLnl9f9A/h8Zu4GjAV+EhFR3lYPXJKZOwHvABPK6/9f4NjMHAEsa3bqo4GFmbk7sDtwTERsV962B/Bd4GPAR4EvlIc6XQDsA3wc2D0iPtfGJbWKIzOnAo8Dh2bmxzPz3dVOlCRJUidTV+sApFXYB5iamX8ByMwF5brh1sxcDjwfEVuU2wbwo4jYG1gObAO8v+3lzHy6/PoJYHD5PomNM/Oh8vpfUypaAD4N7Nzs3o9+lAqEJcAfMnM2QERMoVTsvAfcm5lvldf/J7A3cGuL62kVx5qlRZIkqXOzyFBnFkC2sf6fLdoAHApsBvxrZr4XEa8A67fRfhnQp9l+KzvvCZk5fYWVEWPaiCfbOdbK4n4/DkmSpG7H4VLqzO4GDo6IAQAR8aFVtO0HvFkuMMYCH17VgTPzbeBvETG8vOqQZpunA9+IiF7l824fERuWt+0REduV78X4EjADeBQYXb73oicwCbhvNa7zb8DGq9FekiSpU7MnQ51WZs6MiHOB+yJiGfDUKpr/J3B7RDwOPA38dwWnOBr4VUQsAu4FFpbXX0FpKNOT5fs63gLev8fiYUr3gXwMuB+4JTOXR8RpwD2UejUaMvO2ii8UrgZ+GRHvAiO8L0OSJHV1kdnWaBSp+4uIjTLz7+XXpwJbZea3VtF+DHByZh64sjaSJEmyJ0Prts+UeyDqgFcpPa9CkiRJa8meDElSTUXEVZRmd3szM4e2sT2AnwEHAIuBIzPzyepGKUlaHd74LUmqtauB/VexfRylaaTrga8Bv6hCTJKktWCRIUmqqcy8H1iwiibjgWuz5BGgf0RsVZ3oJElrwnsyJEmd3TbAa82W55bXvd680cKFCx3/K0kdoF+/fpU+E6xJVYuMrvAGMH/+fAYMGFDrMKrO6163eN3dw5r80u+i2rrOTv9+IknrModLtbBgwap67Lsvr3vd4nWri5kLDGy2vC0wr0axSJIqYJEhSerspgFHRMlwYGFmvt7eToLGxsZah9CpmI/WzMmKzEdxvCdDklRTETEFGANsGhFzgTOAXgCZ+UuggdL0tbMoTWF7VG0ilSRVyiJDklRTmTmpne0JfLNK4UiSClDVImP+/PldYkz0utpV5nWX9N1yW5b06F3Rvr2XL+Gvb8ztiLA6nF/vrm/YsGG1DkGSpDZVtcgYMGBAp5/ZpbGxkfr6+lqHUXVe9wfmL+/N/31lUUX7Hzy4H/X1fToitA7l11uSJHUkb/yWJEmSVCiLDEmSJEmFssiQJEmSVCiLDEmSJEmFssiQJEmSVCiLDEmSJEmFssiQJEmSVCiLDEmSJEmFssiQJEmSVCiLDEmSJEmFsshQt3PxxRez1157sddeezFixAg+9KEP8fbbb/OPf/yDffbZh5EjRzJ8+HB+9KMfNe3zjW98g5133pm99tqLA/fek3kvPgvA7Mcf5My9P8LFh4zh4kPGcPflP27a58FfX8b+nxjG8OHDufTSS5vWP/vss3zqU5/iE5/4BF/60pf461//2macv//97xk2bBi77rorF110UdP6V155hX333ZfddtuNo446iiVLlgDwz3/+k6OOOopdd92Vfffdl1dffbXQvEmSJBXFIkPdzoknnsiMGTOYMWMGp59+OiNHjmSTTTZhvfXWY9q0aTz44IM88MAD3H333Tz77LNN+51zzjnMmDGDO+5/lK3/18ea1g/++HBOvP5eTrz+Xvb92skAvDHrBR675f/jlt/fz4wZM5g+fTovvfRS0/nPOOMMHnroIQ488EAuvvjiVjEuW7aMk08+malTp/Loo48ydepU/vu//xuAM888k+OOO44nn3yS/v37c9111wFw3XXX0b9/f5566imOO+44zjzzzI5KoSRJ0lqxyNBae/XVV9l999054YQTGDFiBMcccwz33nsv++23H7vtthtPPPEEixYt4pvf/CZjx45l1KhR/OY3v2nad9y4cey9997svffePProowA88MADfOYzn+GII45g991355hjjiEzVzu2m2++mYkTJwIQEWy00UYAvPfee7z33ntExBpd81sv/4mBH/tX+mywAXV1dYwcOZI77rgDgFmzZjFy5EgAxo4dy+23395q/yeeeIKPfOQjDB48mN69ezNhwgQaGhrITO6//37Gjx8PwKRJk5py1dDQwKRJkwAYP34899133xrlRJIkqaPV1ToAdQ+zZ8/m6quvZsiQIYwdO5abbrqJu+66i4aGBn7yk5+www47sPfee3PJJZfwzjvvsO+++zJmzBg222wzbrnlFtZff31eeukljj76aO69916gNOzo4YcfZquttmK//fbjkUceYcSIEZx22mk88MADrWKYMGEC3/nOd5qWFy9ezO9//3smT57ctG7ZsmWMHj2al19+ma9+9asMHTq0ads555zDBRdcwJ5778NHjjyVut7rATDn2cf52ZfG0HezLTjgO2exxUd3YIuPDmH6JT/i7QXz6bNBT373u9/x8Y9/HIAhQ4bQ0NDAZz7zGW699Vb+/Oc/t4r19ddfZ5tttmla3nrrrXniiSdYsGAB/fr1o66urmn966+/3mqfuro6+vbty4IFCxgwYMAafc0kSZI6ikWGCvHhD3+YnXbaCYAddtiB0aNHExHstNNOzJkzh3nz5nHnnXfyH//xH0Dp/oK5c+ey5ZZbcsopp/Dcc8/Ro0ePpiFHALvttlvTH9Uf+9jHmDNnDiNGjOC8886rKKa77rqLPffck0022aRpXc+ePZkxYwbvvPMOhx12GMOHD6e+vp4zzjiDLbbYgiVLlnDst07ivqv/g32/djJb77Az3//Nk6y3wUb894zfcd1JR3DybX9g849sz+gjT+DLXziQvhtuwNChQ5sKg5///Od8//vf58ILL2TcuHH06tWrVWwr64Foa/37vS2r2iZ1dRGxP/AzoCdwRWae32L7IOAaoH+5zamZ2VD1QCVJFbHIUCHWW2+9ptc9evRoWo4Ili1bRs+ePbn22mupr69fYb/zzjuPzTffnBkzZrB8+XK22GKLNo/Zs2dPli5dClBxT0bzoVIt9e/fn7322ouHH36YcePGseWWWzadc+L/PpyzL/wJAOtvtHHTPjvs9SluO+/7LHp7PhtuMoDdP3cYk7/9dQb0WMLZZ5/N1ltvDcD222/PLbfcApSGTv32t79tdf6tt956hR6OefPmsdVWWzFgwAAWLlzI0qVLqaurY968eU2xvb/PNttsw9KlS/nrX/+6QgEldVUR0RO4BPgUMBd4LCKmZebzzZr9O3BjZv4iInYEGoDBVQ9WklQR78lQVey7775cfvnlTZ/G//GPfwTgr3/9K1tssQU9evTg+uuvZ9myZe0e67zzzmu6sbv5v+YFxsKFC3nwwQc54IADmtb95S9/4Z133gHg3Xff5b777mPw4MEAvPHGG0Cpt+B3v7mdLf9lBwD+9pf/aYr5teeeJHM5G/T/EAB/X/BWaf1rr3H77bc3FTRvvVVav3z5ciZPnsxRRx3V6hp22203XnrpJV555RWWLFnCzTffzLhx44gIRo0axW233QbAlClTmq5h3LhxTJkyBYDbbruNvffe254MdRd7ALMyc3ZmLgGuB8a3aJNA3/LrfsC8KsYnSVpN9mSoKk455RROO+00Ro4cSWYyaNAgbrjhBr761a9y+OGHc9tttzFq1Cg23HDDQs53xx13sM8++6xwvDfeeINvfOMbLFu2jMzkc5/7HKNGjQLgmGOOYf78+WQm2w/dhbHfKg3Jevb3t/Po1Kvp0bOOXuutz6TzLm/6w/4/Tz6KGxYvZL1ePfnxj39M//79AZg6dSpXXHEFAJ/97Gc57LDDgNI9FSeeeCI33XQTdXV1TJ48mQkTJrBs2TIOO+wwhgwZAsBZZ53FV77yFX74wx+y8847c/jhhwNw+OGH8/Wvf51dd92VTTbZhKuuuqqQXEmdwDbAa82W5wJ7tmhzJvDbiDgB2BD45KoO2NjYWGR8XZq5WJH5aM2crMh8lLQcfbK6opqz0yxcuLDTT4XT2Ni41kntirzuD8xf3psbX1lU0f4HD96QAT2WdERoHcqvd/fQr1+/btGVFRFfBPbLzK+Wlw8H9sjME5q1OYnSe9ZPImIEcCUwNDOXv9+mK7zHVFt3+55fW+ajNXOyIvPRtjV5v3G4lCSp1uYCA5stb0vr4VBHAzcCZObDwPrAplWJTpK02iwyJEm19hhQHxHbRURv4BBgWos2c4B9ASJiCKUi462qRilJqphFhiSppjJzKXA8MB14gdIsUjMj4uyIOKjc7LvAMRHxR2AKcGT6NEpJ6rS88VuSVHPlZ140tFh3erPXzwMjqx2XJGnN2JMhSZIkqVAWGZIkSZIKZZEhSZIkqVAWGZIkSZIKZZEhSZIkqVBVnV1q/vz5LFiwoJqnXCPr6uPkve6S9bf+SMX7Ln73XRbMm110SFXh17vrGzZsWK1DkCSpTVUtMgYMGMCAAQOqecrVtq4+Tt7r/sD85b2BRRXtv0GfPgzsgnnz6y1JkjqSw6UkSZIkFcoiQ5IkSVKhLDIkSZIkFcoiQ5IkSVKhLDIkSZIkFcoiQ5IkSVKhLDIkSZIkFcoiQ5IkSVKhLDIkSZIkFcoiQ5JUcxGxf0S8GBGzIuLUlbQ5OCKej4iZEfHrascoSapcXa0DkCSt2yKiJ3AJ8ClgLvBYREzLzOebtakHTgNGZubbEbF5baKVJFXCngxJUq3tAczKzNmZuQS4Hhjfos0xwCWZ+TZAZr5Z5RglSavBngxJUq1tA7zWbHkusGeLNtsDRMSDQE/gzMy8a2UHbGxsLDrGLstcrMh8tGZOVmQ+Surr69dqf4sMSVKtRRvrssVyHVAPjAG2BR6IiKGZ+U5bB1zbN8fuorGx0Vw0Yz5aMycrMh/FcbiUJKnW5gIDmy1vC8xro81tmfleZr4MvEip6JAkdUIWGZKkWnsMqI+I7SKiN3AIMK1Fm1uBsQARsSml4VOzqxqlJKliFhmSpJrKzKXA8cB04AXgxsycGRFnR8RB5WbTgfkR8TxwD3BKZs6vTcSSpPZ4T4YkqeYyswFoaLHu9GavEzip/E+S1MnZkyFJkiSpUBYZkiRJkgplkSFJkiSpUBYZkiRJkgplkSFJkiSpUBYZkiRJkgplkSFJkiSpUBYZkiRJkgplkSFJkiSpUFV94vf8+fNZsGBBNU+5RhobG2sdQk143SXrb/2Rivdd/O67LJg3u+iQqsKvd9c3bNiwWocgSVKbqlpkDBgwgAEDBlTzlKutsbGR+vr6WodRdV73B+Yv7w0sqmj/Dfr0YWAXzJtfb0mS1JEcLiVJkiSpUBYZkiRJkgplkSFJkiSpUBYZkiRJkgplkSFJqrmI2D8iXoyIWRFx6iraTYyIjAin1pKkTswiQ5JUUxHRE7gEGAfsCEyKiB3baLcxcCLwaHUjlCStLosMSVKt7QHMyszZmbkEuB4Y30a7c4ALgX9UMzhJ0uqzyJAk1do2wGvNlueW1zWJiF2BgZl5RzUDkyStmao+jE+SpDZEG+uyaWNED+Ai4MhKD9idnuy+tszFisxHa+ZkReajZG0fXmuRIUmqtbnAwGbL2wLzmi1vDAwF7o0IgC2BaRFxUGY+3tYBfbJ7iU+5X5H5aM2crMh8FMfhUpKkWnsMqI+I7SKiN3AIMO39jZm5MDM3zczBmTkYeARYaYEhSao9iwxJUk1l5lLgeGA68AJwY2bOjIizI+Kg2kYnSVoTDpeSJNVcZjYADS3Wnb6StmOqEZMkac3ZkyFJkiSpUBYZkiRJkgplkSFJkiSpUBYZkiRJkgplkSFJkiSpUBYZkiRJkgplkSFJkiSpUBYZkiRJkgplkSFJkiSpUBYZkiRJkgplkSFJkiSpUBYZkiRJkgplkaEu4Re/+AUjRoxg+PDhXHrppU3rzzvvPIYMGcJee+3FXnvtxW9/+1sAHnnkET7xiU8wduxYZs+eDcA777zDF77wBTKzzXN85jOf4amnnmpafnveHH76xVEAzH78Qc7c+yNcPGks/88XPsHvL5tcOs+M+xk0aBCjRo1i2LBhjBs3jrvuuqtDciBJktRV1NU6AKk9zz//PNdeey133303vXv3ZsKECey333589KMfBeC4447jhBNOWGGfn//851x33XXMmTOHK6+8knPPPZfJkydz0kknERFrFMfgjw/nyIt/zZJ3F3HxIWMZMurTDOoLI0aM4IYbbgDgmWee4dBDD6VPnz6MHj167S5ckiSpi7InQ53en/70J4YNG8YGG2xAXV0dI0eO5I477ljlPr169eLdd99l8eLF9OrVi5dffpnXX3+dvfbaa63j6d1nQ7YZsgvz577catvOO+/M9773PS6//PK1Po8kSVJXVdWejPnz57NgwYJqnnKNNDY21jqEmuis173BBhtw33338dhjj7H++utz++23M2TIEBobG1mwYAHXXnst11xzDUOGDOHb3/42ffv2ZcKECRx77LGst956nHXWWXzve9/j2GOPbfMa31/37rvv8uUvf5k+G/Xlr0uTZe+9R/RoXYcvemcBc559gn2O+S7//OcCFi1atMJxN9lkE2bOnNlp8/m+zh5fR+lO1z1s2LBahyBJUpuqWmQMGDCAAQMGVPOUq62xsZH6+vpah1F1nfm66+vrOfnkk/nud7/LhhtuyLBhw1h//fWpr6/nlFNO4fzzzyciOPfcc7nqqqu45JJLqK+v58ADDwTgwQcfpL6+ng9/+MOce+659OrVix/+8IdsvvnmK1x3nz59uOaaaxi0y57c+Moi3p43h2u+dWhTHK88/QgXTxpLRA/GHHUiW3x0B9ab+wQbbrjhCrlbvHgxvXv37rT5hM799e5I6+p1dwURsT/wM6AncEVmnt9i+0nAV4GlwFvAVzLz1aoHKkmqiMOl1CUcccQR3H///dx5551ssskmTfdjbL755vTs2ZMePXpwxBFH8OSTT66wX2by4x//mO9973tccMEFnHbaaRx88MFcdtllqx3D4I8P58Qp93DCr+9mz4lHrrTdM888w/bbb7/ax5fWVRHRE7gEGAfsCEyKiB1bNHsKGJaZOwNTgQurG6UkaXVYZKhLeOuttwB47bXXuP3225k4cSIAb7zxRlObO+64gyFDhqyw369//Ws+/elP079/fxYvXkyPHj3o0aMHixcv7pA4n3vuOSZPnswxxxzTIceXuqk9gFmZOTszlwDXA+ObN8jMezLz/R/cR4BtqxyjJGk1OLuUuoQjjjiCBQsWUFdXx49//GP69+8PwOmnn85zzz0HwKBBg/jpT3/atM/ixYuZMmUKt9xyCwDf/OY3OeKII+jVqxdXXnllYbE9/PDDjBo1infffZfNNtuMCy64wJmlpNWzDfBas+W5wJ6raH80cOeqDtid7r1ZW+ZiReajNXOyIvNRsrbDi2NlzwzoCAsXLqzeydbQujpm2+v+wPzlvbnxlUUV7X/w4A0Z0GNJR4TWofx6dw/9+vVbs/mYO5mI+CKwX2Z+tbx8OLBHZp7QRtvDgOOB0Zn5z+bbusJ7TLV1t+/5tWU+WjMnKzIfbVuT9xt7MiRJtTYXGNhseVtgXstGEfFJ4Ae0UWBIkjoX78mQJNXaY0B9RGwXEb2BQ4BpzRtExK7AZcBBmflmDWKUJK0GiwxJUk1l5lJKQ6CmAy8AN2bmzIg4OyIOKjebDGwE3BQRT0fEtJUcTpLUCThcSpJUc5nZADS0WHd6s9efrHpQkqQ1Zk+GJEmSpEJZZEiSJEkqlEWGJEmSpEJ5T4bWWmNjI0cddVTT8quvvsppp53GggULaGhooEePHmy22WZceumlbLXVVjzwwAMceuihDBo0CIDPfvazfP/73wfgYx/7GBtvvDE9evSgrq6Oe++9F4Bnn32Wk046iUWLFjFw4EB+9atf0bdvX5YsWcK3v/1tnn76aSKC888/n1GjRrWK8e233+aoo45izpw5DBo0iKuvvrrpgX4t43/llVcZ/fXvsXjh27xw711Ej2DDD23GF8/6D/putiWzH3+Qa086nA9tPYhre/fgc589sKbxS5IkdTb2ZGit1dfXM2PGDGbMmMF9991Hnz59OPDAAznxxBN56KGHmDFjBvvttx8XXnhh0z4jRoxo2uf9P9Dfd/vttzNjxoymP9ABTjzxRM444wweeughDjzwQC6++GIArrnmGgAeeughbr31Vv793/+d5cuXt4rxoosuYvTo0Tz55JOMHj2aiy66aKXxr79BH3Ya+xn2PuJ4vnXjfZx4/b3sMOpT3H35j5v2Gfzx4Zx4/b3ccf+jNY9fkiSps7HIUKHuu+8+tttuOwYNGkTfvn2b1i9evJiINX848axZsxg5ciQAY8eO5fbbbwfgxRdfZPTo0QBsttlm9OvXj6eeeqrV/g0NDUyaNAmASZMm8Zvf/Gal8Q8a/BE22Xog62+0cdP6997tGvFLkiR1BhYZKtTNN9/MhAkTmpbPOeccdtppJ2666Sb+7d/+rWn9H/7wB0aOHMnEiRN54YUXmtZHBJ///OcZPXo0V199ddP6IUOG0NBQmt3y1ltv5c9//jMAQ4cOpaGhgaVLl/LKK6/w9NNPM3fu3FZxvfnmm2y55ZYAbLnllrz11lsrjf+zE77YtDz95+dy/rhdePrOm/nkNz7osZjz7OP87EtjOOqL4ztV/JIkSZ2BRYYKs2TJEu68804+97nPNa37P//n/zBz5ky++MUvcvnllwOwyy678Oyzz/Lggw/yta99jUMPPbSp/fTp07n//vuZOnUqv/rVr3jwwQcB+PnPf84VV1zB6NGj+fvf/06vXr0AOOyww9h6660ZM2YMp512GnvuuSd1dWt2q9H78R8w/gtN6/Y7/geceucf+fi4CTx8/ZUAbL3Dznz/N0/yrRvu5YivfaPTxC9JktRZWGSoML/73e/YZZdd2HzzzVttmzhxYtMQob59+7LRRhsB8OlPf5r33nuP+fPnA7DVVlsBpaFDBx54IE8++SQA22+/Pbfccgv33XcfEydOZLvttgOgrq6O8847jxkzZjBlyhQWLlzIRz/60Vbn33zzzXnjjTcAeOONN9hss81WGv+mm2/Ratsu+09g5n/dAcD6G23MehuU4h/7qf07TfySJEmdhUWGCtNyqNRLL73U9PrOO++kvr4egP/5n/8hMwF44oknyEw+9KEPsWjRIv72t78BsGjRIu655x6GDBkC0DQ8aPny5UyePLlpNqjFixezaNEiAO655x569uzJDjvs0Cq2cePGMWXKFACmTJnCAQcc0G78f5nzQfwv3H8Xmw3+FwD+9pcP4v/jE491mvglSZI6C8dlqBCLFy/mnnvuWWHWozPPPJNZs2YREQwcOLBp22233cZVV11Fz5496dOnD1deeSURwVtvvdU09GjZsmVMnDiRT37ykwBMnTqVK664AihNeXvYYYcBpT/eJ0yYQI8ePdhqq6247LLLms5/wgkn8JWvfIVdd92V73znOxx55JFcd911bLvttk2zOrUV/9LyursuPoe/vPoSET3ov9W2fO4Hpdmlnv397Tw69Wp69Kxjy74bdIr4JUmSOpN4/xPZali4cGH1TraGGhsbmz5xX5d43R+Yv7w3N76yqKL9Dx68IQN6LOmI0DqUX+/uoV+/fms+5Vk31BXeY6qtu33Pry3z0Zo5WZH5aNuavN84XEqSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSVHMRsX9EvBgRsyLi1Da2rxcRN5S3PxoRg6sfpSSpUhYZkqSaioiewCXAOGBHYFJE7Nii2dHA25n5L8BFwAXVjVKStDqqNoXtWWed5dSCktQBzjjjjC49lW1EjADOzMz9ysunAWTmec3aTC+3eTgi6oA3gM2y2ZuYU9hKUsdwCltJUle0DfBas+W55XVttsnMpcBCYEBVopMkrTaLDElSrbX1CVnLXolK2kiSOom6ap2oq3Tnn3XWWdlVYi2S171u8brVycwFBjZb3haYt5I2c8vDpfoBC5o38AnoktR52JPR2lm1DqBGvO51i9etzuQxoD4itouI3sAhwLQWbaYBXy6/ngj8V1brpkJJ0mqr2o3fkiStTEQcAPwU6AlclZnnRsTZwOOZOS0i1geuA3YFlgG9gOXAFZl5fotjrQdcC/wrMB/4Uma+UrWLqYGI2B/4GaX8tZWTk4CvAkuBt4CvZOarVQ+0StrLR7N2E4GbgN0z8/EqhlhVleQjIg4GzqQ0DPGPmfm/qxpklVXwMzMIuAboX25zamY2VD3QKoiIq4ADgTczc2gb24NSrg4AFgNHZuaT7R7XIkOS1FWUp7v9E/ApSkOoHgMmZebzzdocB+ycmcdGxCHA5zPzSzUJuAoqzMlY4NHMXBwR3wDGdNecVJKPcruNgd8AvYHju2uRUeH3Rz1wI7BPZr4dEZtn5ps1CbgKKsxDz6e7AAAT3ElEQVTJ5cBTmfmL8pTaDZk5uBbxdrSI2Bv4O3DtSoqMA4ATKBUZewI/y8w92zuuw6UkSV3JHsCszJydmUuA64HxLdqMp/QJJMBUYN/yJ3HdVbs5ycx7MnNxefERSve9dFeVfI8AnANcCPyjmsHVQCX5OAa4JDPfBujOBUZZJTlJoG/5dT9a3yfWbWTm/bS4x62F8ZQKkMzMR4D+EbFVe8ddZ4uMiPhORMyMiOciYkq5K75lm4Mj4vlyu1/XIs6itXfdETEoIu6JiKci4ply9drlRcS3ytc8MyK+3cb2iIiLy08TfiYidqtFnEWr4LoPLV/vMxHxUETsUos4i9bedTdrt3tELCsPmVDX4HS3rVWSk+aOBu7s0Ihqq918RMSuwMDMvKOagdVIJd8f2wPbR8SDEfFIeShRd1ZJTs4EDouIuUADpU/y11Wr+zsGWEeLjIjYBjgRGFbuFupJ6UbD5m3qgdOAkZm5E7DSP1S6ikquG/h34MbM3LW87dLqRlm8iBhK6VOaPYBdgAPLX9/mxgH15X9fA35R1SA7QIXX/TIwOjN3pvSp3uXVjbJ4FV73+93lFwDTqxuh1pLT3bZW8fVGxGHAMGByh0ZUW6vMR0T0oPTU+O9WLaLaquT7o47S+98YYBJwRUT07+C4aqmSnEwCrs7MbSkNE7qu/L2zLlqj36nrarKg9APVpzwV4ga07gbrrl2H7V13d+weHAI8kpmLy59q3gd8vkWbNeoK7OTave7MfOj973G6zxCKSr7eUPpU6magu/xsrytWZ7pbVjbdbTdTSU6IiE8CPwAOysx/Vim2WmgvHxsDQ4F7I+IVYDgwLSKGVS3C6qr0Z+a2zHwvM18GXqRUdHRXleTkaEr3qZCZDwPrA5tWJbrOp6LfMS2tk0VGZv4Z+DEwB3gdWJiZv23RrNt1HVZ43WfS/boHnwP2jogBEbEBpU8kBrZos0ZdgZ1cJdfdXHcZQtHudZd79T4P/LIG8WntON1ta+3mpDw86DJKBUZ3L6xXmY/MXJiZm2bm4PKNvI9Qyku3vPGbyn5mbgXGAkTEppT+Bppd1Sirq5KczAH2BYiIIZSKjLeqGmXnMQ04ojy0fDilvx9fb2+ndbLIiIhNKH1yvR2wNbBhuQu5uW7XdVjhdXe77sHMfIHSsJjfAXcBf6Q0jWNz3W54RYXXDTTNPHM08P2qBdhBKrzunwLfz8xlVQ5Pa6ncO3U8pWFuL1Aa3jkzIs6OiIPKza4EBkTELOAk4NTaRFsdFeZkMrARcFNEPB0RLf+g6jYqzMc6o8J8TAfmR8TzwD3AKZk5vzYRd7wKc/Jd4JiI+CMwhdK0rV3674KViYgpwMPA/4qIuRFxdEQcGxHHlps0UCo6ZwG/Ao6r6LjdNF+rFBFfBPbPzKPLy0cAwzPzuGZtfklpyMXV5eW7Kc2R/FgNQi5Ehdc9s9zmtfLy7HKbbvPJV0T8CJibmZc2W3cZcG9mTikvv0hpisd2K/Wuoq3rLq/fGbgFGJeZf6pJcB1oJV/vl/mgsNyU0rzfX8vMW2sQoiRJ3U6X/oR6LcwBhkfEBhERlLrDXmjRpjt2HVZy3d2yezAiNi//Pwj4AqVPJZpbo67Azq696y6v/7/A4d2pwGjvujNzu2ZDJaYCx1lgSJJUnLpaB1ALmfloREwFnqQ0jOIp4PJo9nRZSl1ony53HS6jG3QdVnjd3wV+FRHfoTRcqLt0D94cEQOA94Bvlh82dCxAZv6SUlfgAZS6AhcDR9Us0mK1d92nU5ra89JS3cnSzOwONz+2d92SJKkDrZPDpSRJkiR1nHV1uJQkSZKkDmKRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmS1lkRcUP5wWRPR8QrEfF0G23Wj4g/RMQfI2JmRJzVbNvVEfFys2N8vLx+TEQsbLb+9Gb7fCsinisf69sVxHhkRLzV7FhfLer6JUnqKO1OYRsRVwEHAm9m5tA2tgfwM0rTfy6mNOXpk0UHKklrIyLGUPr9dOT76zLzS822/wRY2Mau/wT2ycy/R0QvYEZE3JmZj5S3n5KZU9vY74HMPLBFDEOBY4A9gCXAXRHxm8xsbCf8GzLz+HbaSJLUaVTSk3E1sP8qto8D6sv/vgb8Yu3DkqTqKX9YcjCtH9JIlvy9vNir/G9N5/4eAjySmYszcylwH/D5cgwfjYi7IuKJiHggInZYw3NIklRz7RYZmXk/sGAVTcYD15bfiB8B+kfEVkUFKElVMAr4n5X1KEREz/JQqjeB32Xmo802nxsRz0TERRGxXrP1I8pDrO6MiJ3K654D9o6IARGxAaUe4IHlbZcDJ2TmvwInA5c2O9aE8jmmRsRAJEnq5Cp6GF9EDAbuWMlwqTuA8zNzRnn5buD7mfl483YLFy70qX+S1AH69esXtY5BkqTmirjxu603NwsKSZIkaR1VRJExlw+6+wG2BeYVcFxJkiRJXVARRcY04IgoGQ4szMzXCzhut9fY2N6EMusW89GaOVmR+ZAkqWuoZArbKcAYYNOImAucQWl2FTLzl0ADpZsXZ1GawvaojgpWkiRJUufXbpGRmZPa2Z7ANwuLSJIkSVKX5hO/JUmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSoqMiJi/4h4MSJmRcSpbWwfFBH3RMRTEfFMRBxQfKiSJEmSuoJ2i4yI6AlcAowDdgQmRcSOLZr9O3BjZu4KHAJcWnSgkiRJkrqGSnoy9gBmZebszFwCXA+Mb9Emgb7l1/2AecWFKEmSJKkricxcdYOIicD+mfnV8vLhwJ6ZeXyzNlsBvwU2ATYEPpmZTzQ/zsKFC5tO1NjYWNgFSNK6qL6+vul1v379ooahSJLUSl0Fbdp682pZmUwCrs7Mn0TECOC6iBiamcvbOmDzN8d1WWNjo7loxny0Zk5WZD4kSeoaKhkuNRcY2Gx5W1oPhzoauBEgMx8G1gc2LSJASZIkSV1LJUXGY0B9RGwXEb0p3dg9rUWbOcC+ABExhFKR8VaRgUqSJEnqGtotMjJzKXA8MB14gdIsUjMj4uyIOKjc7LvAMRHxR2AKcGS2d7OHJEmSpG6pknsyyMwGoKHFutObvX4eGFlsaJIkSZK6Ip/4LUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSClVRkRER+0fEixExKyJOXUmbgyPi+YiYGRG/LjZMSZIkSV1FXXsNIqIncAnwKWAu8FhETMvM55u1qQdOA0Zm5tsRsXlHBSxJkiSpc6ukJ2MPYFZmzs7MJcD1wPgWbY4BLsnMtwEy881iw5QkSZLUVbTbkwFsA7zWbHkusGeLNtsDRMSDQE/gzMy8a2UHbGxsXM0wuy9zsSLz0Zo5WZH5KKmvr691CJIkrVQlRUa0sS7bOE49MAbYFnggIoZm5jttHdA3x5LGxkZz0Yz5aM2crMh8SJLUNVQyXGouMLDZ8rbAvDba3JaZ72Xmy8CLlIoOSZIkSeuYSoqMx4D6iNguInoDhwDTWrS5FRgLEBGbUho+NbvIQCVJkiR1De0WGZm5FDgemA68ANyYmTMj4uyIOKjcbDowPyKeB+4BTsnM+R0VtCRJkqTOq5J7MsjMBqChxbrTm71O4KTyP0mSJEnrMJ/4LUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSClVRkRER+0fEixExKyJOXUW7iRGRETGsuBAlSZIkdSXtFhkR0RO4BBgH7AhMiogd22i3MXAi8GjRQUqSJEnqOirpydgDmJWZszNzCXA9ML6NducAFwL/KDA+SZIkSV1MXQVttgFea7Y8F9izeYOI2BUYmJl3RMTJ7R2wsbFxtYLszszFisxHa+ZkReajpL6+vtYhSJK0UpUUGdHGumzaGNEDuAg4stKT+uZY0tjYaC6aMR+tmZMVmQ9JkrqGSoZLzQUGNlveFpjXbHljYChwb0S8AgwHpnnztyRJkrRuqqTIeAyoj4jtIqI3cAgw7f2NmbkwMzfNzMGZORh4BDgoMx/vkIglSZIkdWrtFhmZuRQ4HpgOvADcmJkzI+LsiDioowOUJEmS1LVUck8GmdkANLRYd/pK2o5Z+7AkSZIkdVU+8VuSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSwyJEmSJBXKIkOSJElSoSoqMiJi/4h4MSJmRcSpbWw/KSKej4hnIuLuiPhw8aFKkiRJ6graLTIioidwCTAO2BGYFBE7tmj2FDAsM3cGpgIXFh2oJEmSpK6hkp6MPYBZmTk7M5cA1wPjmzfIzHsyc3F58RFg22LDlCRJktRVVFJkbAO81mx5bnndyhwN3Lk2QUmSJEnquuoqaBNtrMs2G0YcBgwDRq/qgI2NjRWcdt1gLlZkPlozJysyHyX19fW1DkGSpJWqpMiYCwxstrwtMK9lo4j4JPADYHRm/nNVB/TNsaSxsdFcNGM+WjMnKzIfkiR1DZUMl3oMqI+I7SKiN3AIMK15g4jYFbgMOCgz3yw+TEmSJEldRbtFRmYuBY4HpgMvADdm5syIODsiDio3mwxsBNwUEU9HxLSVHE6SJElSN1fJcCkyswFoaLHu9GavP1lwXJIkSZK6KJ/4LUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSCmWRIUmSJKlQFhmSJEmSClVRkRER+0fEixExKyJObWP7ehFxQ3n7oxExuOhAJUmSJHUN7RYZEdETuAQYB+wITIqIHVs0Oxp4OzP/BbgIuKDoQCVJkiR1DZGZq24QMQI4MzP3Ky+fBpCZ5zVrM73c5uGIqAPeADbLZgdfuHDhqk8kSVoj/fr1i1rHIElSc5UMl9oGeK3Z8tzyujbbZOZSYCEwoIgAJUmSJHUtlRQZbX1C1rJXopI2kiRJktYBdRW0mQsMbLa8LTBvJW3mlodL9QMWNG9gd74kSZK0bqikJ+MxoD4itouI3sAhwLQWbaYBXy6/ngj8V7Z3s4ckSZKkbqndIqN8j8XxwHTgBeDGzJwZEWdHxEHlZlcCAyJiFnAWMMzpbj9QwRTAJ0XE8xHxTETcHREfrkWc1dJePpq1mxgRGRHDqhlftVWSj4g4uPw9MjMifl3tGKutgp+ZQRFxT0Q8Vf65OaAWcVZDRFwVEW9GxHMr2R4RcXE5V89ExG7VjlGSpJbanV1qtQ5Wmu72T8CnKA2hegyYlJnPN2tzHLBzZh4bEYcAn8/MLxUWRCdTYU7GAo9m5uKI+AYwprvmpJJ8lNttDPwG6A0cn5mPVzvWaqjw+6MeuBHYJzPfjojNM/PNmgRcBRXm5HLgqcz8RXlK7YbMHFyLeDtaROwN/B24NjOHtrH9AOAE4ABgT+BnmblndaOUJGlFRT/xew9gVmbOzswlwPXA+BZtxgPXlF9PBfaNiO58v0a7OcnMezJzcXnxEUr3vXRXlXyPAJwDXAj8o5rB1UAl+TgGuCQz3wbozgVGWSU5SaBv+XU/Wt8n1m1k5v20uMethfGUCpDMzEeA/hGxVXWikySpbUUXGU5321olOWnuaODODo2ottrNR0TsCgzMzDuqGViNVPL9sT2wfUQ8GBGPRMT+VYuuNirJyZnAYRExF2ig9En+ump1f8dIktThKpldanU43W1rFV9vRBwGDANGd2hEtbXKfERED0pPjT+yWgHVWCXfH3VAPTCGUi/XAxExNDPf6eDYaqWSnEwCrs7Mn5QfGHpdOSfLOz68Tmdd+50qSeoCiu7JWJ3pblnZdLfdTCU5ISI+CfwAOCgz/1ml2GqhvXxsDAwF7o2IV4DhwLRufPN3pT8zt2Xme5n5MvAipaKju6okJ0dTuk+FzHwYWB/YtCrRdT4V/Y6RJKmaii4ynO62tXZzUh4edBmlAqO7j7dfZT4yc2FmbpqZg8s38j5CKS/d8sZvKvuZuRUYCxARm1IaPjW7qlFWVyU5mQPsCxARQygVGW9VNcrOYxpwRHmWqeHAwsx8vdZBSZLWbYUOl8rMpRHx/nS3PYGr3p/uFng8M6dRmu72uvJ0twso/QHRbVWYk8nARsBN5Xvg52TmQSs9aBdWYT7WGRXmYzrw6Yh4HlgGnJKZ82sXdceqMCffBX4VEd+hNDToyO76YUVETKE0VG7T8j0oZwC9ADLzl5TuSTkAmAUsBo6qTaSSJH2g0ClsJUmSJKno4VKSJEmS1nEWGZIkSZIKZZEhSZIkqVAWGZIkSZIKZZEhSZIkqVAWGZIkSZIKZZEhSZIkqVAWGZIkSZIK9f8DiDnwv0D5A7MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.plot_posterior(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['changepoint_interval__', 'u_log__', 'changepoint', 'u']"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace.varnames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "735599.0"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace['changepoint'].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(725008, 736884)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['ord'].min(), df['ord'].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Date\n",
       "1986-01-02    725008\n",
       "1986-01-03    725009\n",
       "1986-01-06    725012\n",
       "1986-01-07    725013\n",
       "1986-01-08    725014\n",
       "               ...  \n",
       "2018-07-03    736878\n",
       "2018-07-04    736879\n",
       "2018-07-05    736880\n",
       "2018-07-06    736881\n",
       "2018-07-09    736884\n",
       "Name: ord, Length: 8223, dtype: int64"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['ord']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Closing Value</th>\n",
       "      <th>ord</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2014-12-31</th>\n",
       "      <td>53.45</td>\n",
       "      <td>735598</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            Closing Value     ord\n",
       "Date                             \n",
       "2014-12-31          53.45  735598"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.loc[df['ord'] == 735598]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Closing Value</th>\n",
       "      <th>ord</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [Closing Value, ord]\n",
       "Index: []"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.loc[df['ord'] == int(trace['changepoint'].mean())]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}