aboutsummaryrefslogtreecommitdiff
path: root/statistic/try_HMC.ipynb
blob: bb11871c53fea34e12e4a317c5a7299f5e709389 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import deque\n",
    "from scipy import stats\n",
    "import sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 11\n",
    "_lambda = 13\n",
    "B = 1\n",
    "N = 100000 + B\n",
    "epsilon = 0.01\n",
    "L = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "def h(x):\n",
    "    return (_lambda * x) - ((alpha - 1) * np.log(x))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "def h_prime(x):\n",
    "    return _lambda - ((alpha - 1) / x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "def leapfrog(theta, p):\n",
    "    half_next_p = p - (epsilon / 2) * h_prime(theta)\n",
    "    next_theta = theta + epsilon * half_next_p\n",
    "    next_p = half_next_p - (epsilon / 2) * h_prime(next_theta)\n",
    "    return (next_p, next_theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hamiltonian(theta, p):\n",
    "    return h(theta) + ((p ** 2) / 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "def r(theta1, p1, theta2, p2):\n",
    "    return (\n",
    "        np.exp(hamiltonian(theta1, p1) - hamiltonian(theta2, p2))\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100000 / 100001 100.0 %  \n",
      "acceptance ratio: 0.999990000099999 \n"
     ]
    }
   ],
   "source": [
    "last_theta = 2.5\n",
    "data = deque([last_theta])\n",
    "accept_count = 1\n",
    "\n",
    "for i in range(2, N):\n",
    "    if not i % 1000: # 進捗用\n",
    "        sys.stdout.write(\"%s / %s %s %%  \\r\" % (i, N, np.round(100 * (i / N), decimals=2)))\n",
    "\n",
    "    start_p = p = np.random.normal()\n",
    "    start_theta = theta = data[len(data) - 1]\n",
    "    \n",
    "    for j in range(L):\n",
    "        p, theta = leapfrog(theta, p)\n",
    "        \n",
    "    if np.random.rand() < r(start_theta, start_p, theta, p):\n",
    "        data.append(theta)\n",
    "        accept_count = accept_count + 1\n",
    "    else:\n",
    "        data.append(start_theta)\n",
    "\n",
    "        accept_count = accept_count + 1\n",
    "        \n",
    "print(\"\\nacceptance ratio: %s \" % str(accept_count / N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fb932d0d510>"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3wUZf4H8M83jQChJzQBAwhSFASpggoiAorl7rz7iXdiOU899dSzHfZ+1vM8KwoqNrBgQ+lSFKSGFnoHE2poAQIh7fn9sbPJltmd2d3ZMpPP+/XKK7szszPP7Mx+55mnjSilQEREzpIU7wQQEZH1GNyJiByIwZ2IyIEY3ImIHIjBnYjIgVLiteHMzEyVnZ0dr80TEdnSsmXLDiilsoyWi1twz87ORk5OTrw2T0RkSyKy08xyLJYhInIgBnciIgdicCciciAGdyIiB2JwJyJyIMPgLiItRWSOiKwXkbUicrfOMgNEpFBEVmp/j0cnuUREZIaZppBlAO5TSi0XkToAlonITKXUOp/l5imlhlufRCIiCpVhzl0ptUcptVx7fQzAegCnRTthgWzcewyvztiIA8dPxSsJREQJL6QydxHJBtANwGKd2X1FZJWITBWRzgE+f4uI5IhITkFBQciJBYAt+4/j9dlbcKioJKzPExFVB6aDu4hkAPgawD1KqaM+s5cDOF0p1RXAGwC+01uHUuo9pVQPpVSPrCzD3rNERBQmU8FdRFLhCuyfKaW+8Z2vlDqqlDquvZ4CIFVEMi1NKRERmWamtYwAeB/AeqXUqwGWaaotBxHppa33oJUJJSIi88y0lukH4DoAq0VkpTbtYQCtAEApNRrA1QD+LiJlAE4CuEbx4axERHFjGNyVUvMBiMEybwJ406pEERFRZNhDlYjIgRjciYgciMGdiMiBGNyJiByIwZ2IyIEY3ImIHIjBnYjIgRjciYgciME9Cqas3oMfVu2OdzKIqBozM/wAhej2z5YDAC7v2jzOKSGi6oo5dyIiB2JwJyJyIAZ3IiIHYnAnInIgBnciIgeydXA/cqIEXy/Lj3cyyAIVFQq/HTwR72QQOYatg/tdn6/EfV+twraC4/FOCkXonZ+34oKX52DzvmPxTgqRI9g6uO8rLAYAlJRXxDklFKnF2w8BAHYdORnnlMSGUgrfrshHcWl5vJNCDmXr4E5kV3M3FeCfX6zCy9M3xjsplnv8+zXIHjU5KusuLi3HVzl54COajTG4U8J6adoGZI+ajIoK5/2Qj54sBQDsP3Yqzimx3scLd0Zt3S9N24gHJuZizsb9UduGUzC4a7bsP47sUZOxeNvBeCeFNKN/3goAcF5op3DtP+Yqij1+isVZRhwR3HccKIp4HQu3HgAA/JDLAb+IyP4cEdxv+3Q5jp8qi3cyyALh5tL3FhZzJE4iD7YN7r71Kac8Wh0s2HIAH8zfHuMUUSQkws9fO2YR/jFhBVufWGRP4UnkH2a/A6tt2HsUny2OXp2EJ9sFdzERBa4duxhP/7gu+omhhLG70NWE0gmNKApPlGLnwciLGiPR9/nZ6P/iHMvW91VOXmV5eXU29LV5eOTbNTHZlu2COyUepRRW5R1h8zSLDHntF1z48tx4JwMAcKioJOJ17D9ajAcm5uLmj3IsSJE9fbJwBzo8NjWm22Rwp4hNWrUbV771K37M3RPvpDjC3qOJk8Md9r9fIl5HqdaU9YAFzT7tmn14fNJaFJfGtrOlY4J7OAddKf/xTKau3ovsUZNReKLUmoRVA1v3u4Z/2FYQ36IEAFBBzoTyCoUXpm7AgePOa1seLfuOJuZ3FWkdTXVg6+Ae7IdsxqeLduKCl+dgZd6RymkHtdvQbQecNV5NeYXCqTIbVDaGeUhF+7lv3hf4uM3bXIDRP2/FI9+uDm8jRDZi6+AeqWU7DwMAth84bq6m1saue38xznx0WlS3EcnF1qqv/8q3fsUqj4u1pwqtTqCkjGMRhSvv0Al8u4IjsdpBtQ7u1cmCrVHseZtgF8a8EJrw/XXcUtz44ZIopiZ88RgCubxCYfravQHn/+7tX/HPL1bFMEUULscEdzbUoHDM2rAfczYWxDsZui542bqmiGZ9+Ot23PrJsoDzDxyPvPUMxYZjgjvF3sKtBzF7w76QP7dgywF8sfQ3S9OSYDcPcWFFB649hdFvqXO0uBR3TVjBRgtRZhjcRaSliMwRkfUislZE7tZZRkTkdRHZIiK5ItI9Osm1p39NzI1pOeUr0zei979/ivp2RoxZhJvGVbVdNnv3dO3YxfjX1/qVmpFWkofCSWPHf7M8Hx0em4atNnhwzUe/7sCkVbsxZt620D8cwelx6yc5aPvwlPBXYDNmcu5lAO5TSnUE0AfAHSLSyWeZYQDaaX+3AHjH0lSapJRC3qHE6zL9RU5eTMsp35yzJaZN2KzINAdbRzQ6R+07Wox+L8y2fL2hsuqGY8Za1x3Upr1VT7LafeQkvszJs2gLicF98Q/nTm362n0od+Dw0YEYBnel1B6l1HLt9TEA6wGc5rPYlQA+Vi6LANQXkWaWp9bA+CW/4fyXYl9OGY5xv25H9qjJ7NUZhFhc1uL5VVeHtu5/HrsYD07MxbHi+BZ/ROMcF7Z0NxRSmbuIZAPoBmCxz6zTAHhmEfLhfwGAiNwiIjkiklNQYG0lloLCsh2HLV1nNLnHvnFSRsIuu2L1RSNRuXuEhnJcovnN+H7vRSUcyTWaTAd3EckA8DWAe5RSR31n63zE75xSSr2nlOqhlOqRlZUVWkqjrHr83KMjEWJlKEmYvWE/nv1xfdTSkgiUUjiW4MNgf/jrjngnwdFSzCwkIqlwBfbPlFLf6CySD6Clx/sWAKI6uPZVb/2KkxG0Drhp3FLM3sBHdTmRmVKAhTZ84tax4lLUSU81teznS51V1k6hM9NaRgC8D2C9UurVAItNAjBSazXTB0ChUiqqo0hFEtgBRC2wV+vejxaUrQZaRXFpOUrLo/fdjp23DSdLEnd4hvmbD+DsJ2dg/uYDQZdzf30bPSpW3e4YvxyfL7G2CWqieu+XrViwNfh35XRmimX6AbgOwEUislL7u1REbhOR27RlpgDYBmALgDEAbo9OcoNIkALfV2duitq6/ztzE16YugElZRV4fur6uFeUuVlRuWVUDt7hsWn407sLI95OIM9OXo9XZ26M2vojtWTHIQBAzs5DuvPNFI1Nzt2DUd8EH1cnEYrYgjGbf/j3lA24doxv1WD1Ylgso5SaD4MiTeWqDr/DqkRFW5nFOcCvcvKgFPCnni2xO8y20/uPFaOkrAItGtQKuMz/Zm0GALTJrI13f96G4pJyPHXlWWFtz45W/KY/ZgxgTSWpUx/VWFxajjo1TJXAxl15hUJpeQXSU5O9ppeVV6C0vCqymz3c4xf/hmt7tzK9faUUVuYdQbdWDUx/xgxB7POfjumherCoBOv2eNfznghQG//C1A2WbvuBibl48OvcsD7rbibW67lZpp98U1rhujiVlCfI7QrFzGs/bTbVtHDcgh2Vr3s9Nwtvz91qaTqi1YT3xnFL0eEx/wHubvt0GTo+HvrAdw/rjAC6MsDAcgAwcVk+fvf2Akxdbf9nEzgmuA/73zxs8ClnPBKge/PPm+I/lkgkOU13h5VEk+iXGqd0KQinFcyklfZ4ePgvAX6bP623ro7sqrd+DThvq/ZMgu1xfsyhFRwT3O3MTJdxz5xSJBcnpRQ+mL8dR05YNwBUIpTTJkASYiZ4b97Q1lVWXhH1jnTu1e86chLtHzX3qLkiC4vI9h8rtvxu3Q4Y3DVWBahw1uMeVx4A3p+/3e+hGst2HkLrh6wZE2P5b4fx9I/r8K+vc7H8t8OYtzl2dzF5h04YtnhxB4JVeUeqVVfxSIVz3pWUVeCMR6bihWlVgS+SO8rtB4qCFnm4t2lG5yemh50OXw9OzMXon60tlgIQdv1arDC4w9rb9b1BRtVzt26ZvSHwGBfP/LgOb8/xPhH/8I5+K5FwfoentB9X4clS/P7tBbjuff+xzH/eVIDznp8V8iiDwb7HIydKcP5Lc/D492t15/s2Tb0yyK2zruqUddf0f3E27vl8Rdif/27FLgDAuz+HMYCXjoGvzA1a5BEv0WiePDl3D857YbZh09R4qpbBPZr5wcXb9ZuqAcBXOa6RIT1HUtQTSauN7FGTw/pccWk5HvpmNQ4eP4Wnf1iL3YXFyDf50At3XC1XCt+v3IUKnQvXsWLXPkV6p3Dvlyvx4rTgt9gRPA8q7E9aJZRHIeYfPonvfMrSp6/dqxtw9EbbnLspPp34Qj0+nuMAJUq9yYrfXHfb6/f4dtZPHNUyuCeSBDlX8cOq3Ziw5LeIyibfmbsVd3++El8vt3Z4Y8/v6Jvlu/COxS0/3MzeCVVUKN0LWDgmrdrt1arrg/k7Il7fX963tn333I3+F4GuT83QbYkSDT2e9R++Ov6X4cRXLYP7lv2JP+Z1KBZuPYhjxaURlQEqn/+ROKxTWRvLHNeExdHthXn2k9Mx4JW5Ea3D8/s4VFT1fZnpmHb2kzNwNEod2PSC5trd/rnTo8VlGO/xPUcrPVZWrFY39ujZECbPeHKoqATpqUmolRbaLh8qKkHD2mmV70vKKnDg+Ck0r1/TolSGTwAcLirBiDGL0L5JBpJi2GzFzI/5tk+WYWXeESx6eBBmrHM9lzP/cPALULlS2LTPv+u8Ec89DzRujFXXl6KSchQFeW5Abv4R1K+ZhlaN/DuklVcoLAlSdGe2PfrmfVUZlBHvLcKaXYWmPufJymFzuzw5w7J1efKtWC2xuAPitoLjqFczFY0yanhNT5Tin0hUm5x792dm4pL//hLSZ2au24fuz8zEQu3h0uUVCu0fnYrzXpgdsINUrLkrSDftO+7Xzt8qG/ceq3zow6myclRUKKzOrwom/wkw5MK0tXux92gxCk+W6ub+9Px35qaQj5NZr0yP7vACny3eiexRk3HFm78GfP7p6J+3YsSYRZi/xbqKuIXbDib8CJBWcTdEsCofc9F/fsYFHs+A8Fxv4clSZI+ajBlBHhieyKpNcAeMc42+lmrjeeTmu5p3TV1T1WutuNS6HESilx8Oee0XPDgxF6fKynHmo9O8ms4F4r4gAq7y2W+1lhlGonWB2rL/GDaGcEewce+xkMvV35q9xXCZrVqR4L6j0X9WaXUV6l1MUYAB49zFt2/P3YqXpm0I6+4onqpVcA/G9xZ1xrp9fp07Qm1SFe7TfjxTEu+mVp5NO4tLXPtvNLKgQDBizKKopsvTj7m7cbQ4eM510irj7uTu7z03/wiGvPYL3olC22i3yEe4DO3Co1vMoJerSPSchgnD35hvarm35hhfjAHX3cLbc7fid28nXjPPYGwX3GN17r0zdyu+iHBM7PEWVOwdLAp8gYjk1nTRNv1y31NlFZVdsPcWnkKf52fpLpdIFV13jg+/rbeeXdodnmfRk9UCff+xFMrp86+JuZiQAMMFe16jFmw5gB0HjIcJ0OtTUlZegZeDFNO9PH1j5edi+dB2K9kuuEdD4clSlOtkbYxyg1bQy1EdKy5D3qET2G7ixA3mie/X+E1zB25PJ0rKKssVf1hV1W56+W/6jy08WlwWtPdoIgxHEI5I012q853kHz6B7FGTA46ZEg/zNhcge9Rk/JjrfTfz3YpdATs0fZGTh4cMhguOLcG1YxebarW0R6djodHQxwCwZLv9HujiicEdwFM/rMO/J4f32LWPF+6ofG1V2+cvcvJw/ktzMNDEiRssIH20cKfftMe+8w/4j367xtKBmRJNcWk5vl9prsw/EgXH/O+ylmvDFH+Zk4dDRSX4xmTdQzR9tMD/vACA17UhpX0t3xn42cSLLXqilVIqZu3mAZjqi+EeimHNrsCNAXJ3FSJ71GTkBWk95bmuWHJ0cA9lQKRgT3Z6fuoG3a74h4tKvLrTP/mDftf6SMTiVjjPZE9UT3eMXx7xdvcdLcbrs8wNYRuOsvIKTFuzFy9P34idB0PfR6vdHeZQASdLyv0Cbyyb6s0K8tSy/3svcN3K9LV7sdPk6IrHT5VZUoypx0yxypETJUGHDgGA0nKFw0XefTjcd7qJMNKsL0e3c7fSjoNFXrlkpRR8q8SiEYgToWxWjxU3Kf8YvwJLdhzCRR0ah72OQD+qv45biib10jF+8W9ISYpurun1WZvxsc5dki+9nL0ZZiv+golHqfGtnyxDcpJg7v0DLFtntI5kn+dnobi0An8f0DbocgP/MxdPXdE54PyRHyxBnfQUPDG8ExrXTbc6mSFxfHA/WVJe+XALK63fcwwdmtXxmlZq8PCMn9bvw3V9T/ebHu4t28HjJej9b/0KT0+Rlt2Hyuz+nCh11WmEkwvdfqAIzeqlY2mADkGeuc0yk1ci3xZT09buxfYDRWidWVt3+YoKhaQkMf1oxXD2s6JChTTejBnRfBatL726mZ/WR/48glBGeTTzvbubNnueAXp3lEdOlOLuz1cGXI+7bmVy7h58cUsf9G7TyHQ6reboYpldh0+i4+PTvHrPhTuwli8FFfKPdcHWg/j7p5EXZ7iZaRP+57GLDMvuA+1HuLf+sShdHPjKXNz31SpL1/nJop2Ytsa7kjHYdzdxmfkxdDbvOx5SO3u3QK2VIsmJW12EYFSc4SvXglZInkMLK6WwYMuBgHVe4Z7HoXzu2cnr/KYZDXAXbY4O7rd9uixq677/q1UhD4kL+OdaImlmVWbijuTXLeFXeP33J+8caVGceuUWntQf6mBy7h68P3+7pdu6Tefi+/bcLSg8WYqSsgqvCulA6fLkHu/HKLBPW7MXk3P92+LvD7Mox9fBMPtc+NK7KYuk/iV71OSAnYP0MmJ6v5aZ6/bh2rGLMW7BDtz3pbUXfLP0OjUuD/LM31hwdLGMlb1Ife07esp0r8tgTgboHWfGcYuaauYEaQ3h6ZNFxuXK4Xhuin+ux5O7h7CeYBXhVnlp2kZ8lZOPey5u5/UdPDdlPfoEue2uqFCmKxSjmREBgMMBHjnpFslIm5GcwwAwZfUew7LuYNwX0J0Hi3RbwegNZFcdODrnbiWlrB1oya37MzOj/pgzq1idTHcTM6NKY6Px72Nh+4EiHDzuHyRu+ND/YSdubR6egglLIusIFysvTtuAbTGum3H7IXd3VIdjOKFz8bHJTy4ijs65RztXZ8UJEkmrE6PcmNXM1vvatROTkWidT6vygpdBj5nnX/QUr+Ckt12jCl+jtOYdOomLX43OYHHx9o8JK7w6BsaSo4O7lYb9bx7q1UwN6TPhNn2j+NLrrRxNh4IMMWEHer2eoyWW+YbXAnTqCkTv9x6vwA6wWCYkvhVoC7YGH9Rrp0GvNbdPo1SWHS92z7hbPVZNdbcySJ1JqGJ52Q31Ae3TE2xoYAb3CAR7XqqbmZEkV0VxgCorVYdyymDW7vY/TgeLqmdlXSjummDdxVKvfsrKMaAiKVJ8VGdoj3hicI+yCgdFRLPNNuMxjkYsTFkdpZxZGN+XXSrhY8FsJzIzotFoIl5sF9ztlFPauv943Nrdkn2s32PuKVUUnFMzFeGyXYVqOB2H4mXXkZPYFcFDq6l6iNaAWU4U7CHwVtzNxKujXjTYLudOiY8ZqOgrsKjHqd0Ea5kzY13kY9Z8+OuOiNeRKBjcyTyTGSPG9uhjix5/eg/lqM4Y3Mlyj31v/bj2RBQaBncybYfJcVKIKP4Mg7uIfCAi+0VEtxGniAwQkUIRWan9PW59MikRTF8beZkmEcWGmdYy4wC8CeDjIMvMU0oNtyRFREQUMcOcu1LqFwCJ+aw3IiLSZVWZe18RWSUiU0Uk4AMGReQWEckRkZyCgsR7oCwRkVNYEdyXAzhdKdUVwBsAvgu0oFLqPaVUD6VUj6ysLAs2TUREeiIO7kqpo0qp49rrKQBSRSQz4pQREVHYIg7uItJUtEEdRKSXts7wH9xJREQRM2wtIyITAAwAkCki+QCeAJAKAEqp0QCuBvB3ESkDcBLANYpD1hERxZVhcFdKjTCY/yZcTSVjgiO/EREZs10PVYZ2IiJjtgvuRERkjMGdiCjGjpyI/kOHGNyJiGLsucnro74NBnciohgrq4h+g0LbBfeiU855DBYRVU/frtgV9W3YLrivyDsS7yQQESU82wV39o8iIjJmu+BORETGbBfcmW8nIjJmu+BORETGGNyJiByIwZ2IyIFsF9zZWIaIyJjtgjsRERmzXXBnxp2IyJj9gjvLZYiIDNkuuBMRkTEGdyIiB2JwJyJyINsFdxa5ExEZs11wJyIiY7YL7oqNIYmIDNkuuBMRkTEGdyIiB2JwJyJyINsFd7aWISIyZrvgTkRExmwX3JlzJyIyZrvgTkRExhjciYgciMGdiMiBbBfc2UOViMiYYXAXkQ9EZL+IrAkwX0TkdRHZIiK5ItLd+mRWYYUqEZExMzn3cQCGBpk/DEA77e8WAO9EniwiIoqEYXBXSv0C4FCQRa4E8LFyWQSgvog0syqBfumJ1oqJiBzEijL30wDkebzP16b5EZFbRCRHRHIKCgos2DQREemxIriLzjTdDLZS6j2lVA+lVI+srCwLNk1ERHqsCO75AFp6vG8BYLcF6yUiojBZEdwnARiptZrpA6BQKbXHgvXqY6E7EZGhFKMFRGQCgAEAMkUkH8ATAFIBQCk1GsAUAJcC2ALgBIAbo5VYIiIyxzC4K6VGGMxXAO6wLEUG2ImJiMiY/XqoMrYTERmyX3CPdwKIiGzAdsE9N/9IvJNARJTwbBfcS8uZdyciMmK74E5ERMYY3ImIHIjBnYjIgRjciYgciMGdiMiBGNyJiByIwZ2IyIEY3ImIHIjBnYjIgRjciYgciMGdiMiBGNyJiByIwZ2IyIEY3ImIHIjBnYjIgRjciYgciMGdiMiBGNyJiByIwZ2IyIEY3ImIHIjBnYjIgRjciYgciMGdiMiBGNyJiByIwZ2IyIEY3ImIHIjBnYjIgRjciYgciMGdiMiBTAV3ERkqIhtFZIuIjNKZf4OIFIjISu3vZuuTSkREZqUYLSAiyQDeAjAYQD6ApSIySSm1zmfRL5RSd0YhjUREFCIzOfdeALYopbYppUoAfA7gyugmi4iIImEmuJ8GIM/jfb42zdcfRCRXRCaKSEtLUkdERGExE9xFZ5ryef8DgGylVBcAPwH4SHdFIreISI6I5BQUFISWUiIiMs1McM8H4JkTbwFgt+cCSqmDSqlT2tsxAM7VW5FS6j2lVA+lVI+srKxw0ktERCaYCe5LAbQTkdYikgbgGgCTPBcQkWYeb68AsN66JBIROcv57TKjvg3D1jJKqTIRuRPAdADJAD5QSq0VkacB5CilJgG4S0SuAFAG4BCAG6KV4DZZtbGtoChaqyciirr6tdKivg3D4A4ASqkpAKb4THvc4/VDAB6yNmn6MjNqMLgTka3pVWRajT1UiYgciMGdiMiB7BfcfRthEhGRH/sFdyIim5MYFLrbLrjXrZka7yQQESU82wX3B4acGe8kEFVbHZrWiXcSyCTbBfeaqcnxTgIRUUTYFDLOBndqgmt7t4p3Mhyvdlr4F+znf3+2hSmJnivPaW7p+u4a1M7S9ZH1OjevG9ftOzq4ZzeqVfn64o5NQvrs6yO6YczIHjj7tHoAgNPq17Q0beE6o3GGpetLSYo8D/GPi86I6PP9Dbpib3/+0ojWH6tjN7RzU93pTeum43/XdMOOFy4zXEev1g1NbauLdl6Sy5z7B+hOv7C9/hhWbbNqW56Gz27u7fX+kk7650Os2C64h1LL/MxVZwEARvY9HWOv7xHSdsINCFd0tTaH5ist2fuQ1auZig9uCG3fACAtxbWeey6OTQ7wWe1YAMBr/3eO17xBIV54zWjZsOr4Na+f7jUvM8PV9fvewe0rp6WnGv8UburXOuj8Rhn+XcpH/6U7Ztx7geG63fq0aWRquYs6NDa9zmjqf0b0xkgJliFrWjfd6+6ldaZ+sG5cp4bu9KFnWR94Wzao5fW+Q7PA9RPR2L4v2wX3UHRr1QCf/rU3HrmsY8if7RjkwFht07PDMPu+C3Fjv+yQP7vkkUFo1zj0tP7t/KpAdcfAtgCACwLkcsLhe3HMzKj6kXXyuV1N8rliP6pzvL66rS8WjLqo8uJ2cUdXcMuokYKxI3vgiq7N0e+MqsA478GLgqTOtb1relYNdtq+ifd3+OENPVEjxfvnEej4BMoIzHtwIIae1Qx1061v4ZVkcMcVyg1ZKMs+cXlnDDiz6jy5XTt3zBp3Y0/Ty2bV0R9/pWnddCx6eJDXxRkAnvvdWXji8k5e0wacqX8RvG+wfsOMLI+LwUc39dJdJlBGQATo06bqzmtI56aVFdApSYKXru5SOW/oWc38Pm81Rwb32we0xfd39ENGjRT0b5eJGimhlek+cmlH1ErzHnanR3YDdG9V38pkVkpLSUKbrAw8cXnnymkPDDkT2/59Kb66ra/Xsr5lt6HsW6BcjFuv7AYB5zWpWwNzA9z66vHNWQ7pbC53/tjwTvhrf/8ccs/shmhevyayM125o7sHtcfzvz8bl53dDBd3aoLXR3TDZzf3QW+PYg13Dl18qq8ev7wT6qanoEHtquDxh+4tvJYZ2KExlj56cdC03ty/Ne4ceEbVdnyCZMuGtXQ+FRv927kC8O0DAgfftlm18elfe2Pd00Mx45/m7i46Na+LcTdWBT2BYP6/BppOV6BgCwDj/+ZdrNGwtn5w/+6OfrrT/9z7dNzoc3cVKIcc6OLoObVuuv7QWxueGaY7Hai6u3Nf8L+7ox8uPbspfn5wIP7UI7bPMLJdcE/Vcm5dW1YFWt+c0439WnvNdzNTnjnhb328gsugDo2RUSMFt13YFt/crn9SeWpQK3gurZ2JMvMmdWvgjoFnIClJ0DO7Id68thsAoE1mbdxyQRvDzwdyx0D9snH3hSw9SEukibedh+wAt76+Pr6pl1cOavWTl0CClKe18Sj/PKNxBkQkYFmpW1pKEkb0auX3I/3opl5YpgXlS8/Wzx1d0bU5cp8cUnkuAUArnUBcNz3Vr5z8TI8c/qPDO+H+IWfiqSvPQteW9dGhaegVaH7npAq9C/Y7f+7uN230X7pj1n0Xep0vnZr5p69/u37wWOwAAA3QSURBVEykpyZXbtaznsozaE++qz8m39Uf9Xz6mYgALTyKI36690LT6e7ok57z2mbi3sHt8eDQM5Hz6MW6x69Li3poWi/db7qeVU9cguQkwZ0Bzvtg6gQI7EYGdmiMizs2wZiRrqLS9NRkvP3ncytjVKB6mWiwXXBvWi8dr/yxK8aOrCpn9i2nzAqQQx3vU+Hh6/KuzdG3bSOvgNG4bjrWPDXE70T09Z8/dsXc+wdg1LDgRUBGJ+a0e87H1Lu9c1F9tf07v12mV5B0nyhN6po72Ztp23YXabjdfH5r3H9Je4zsmx3ws4FyOp63oYDrjuOC9llI8QicdXyKJXzjl2fxx2la+XjPAHcRv+vmymEHOsbpqclopBUBubdzrraul6/ugveu032ODAYGKcN2F1dl1EjBtHvO95t/Tsv6+P6Ofn7FOHqaasdqcCfXnUyDWqlYMCpYEZK/h4Z18Hp/UUf/tNdKS0HbrAzUr5WGlY8Pxh0D2+J9n7oZz6IytzSPfWjRoBYuaJ+FC9tnoXPzeujc3L8S1/esCFbh7173Q8M64OObemF4F//gfdegdrh9wBnIzKhRecfl2bbet9ISgNfdmiffC5Gbbws43zuPD2/siWn36N/JvKC1zlr88CDdStzU5CSMvb6HX9Gj2+jrzjVVsW6F8C5PcXb1ud630MrkgDMpyUmY8Lc+yDt8Ag9OzK2c3rt1QyzefsjUOhrVTsPBohK/aX/wSZOeLi3q4Y0R3XDO0zMBAE9c3skvoOrl/hpl1MD8fw2sDAxu7jifluLarxFjFnnNX/zwIDwwMRfX9z0dd01YgcGdmmDRQ4NQv1YqjhWXIWfHYYzo1Qo1UpJx50XeFauPD++Ep39cZ7hPntpk1g54d2CkTnoKGtVOwxkG9Qe3XdgGN/XPDqk4qlm9dFM/qOv7no71e47hLJ+WKGNH9sDewuLKYpwvbumDejp3aO5WP+e1bYSruuk9Ztg175sVu7wucM2DVN6/9IcuWJF3GBOWVD3G+C99Tq983bJhTaQkBb+o1K+VhgeGuC4Im58bhldnbsI7c7eiZ7Z/EdalZzdDq4aF2Lz/OADXXVgwwe7I3F79U1fc++UqNNTGML/1QldRkVHu2PN3veSRQUhPTfbLKMy+70LTmZv7BrfHP3SakLbwqQgdqBUd7T9aXDnt1gvb4MSpcvyuu+u4mt1mPNkyuPt6cEgHfLN8F4DAlSBufds2Ql808gruz151Fgb/9xe/HK2exnXT/YJ7IM3rpeNYcRmOnSoDAIzsm436tdLw0tVd8ODEXJzfLgvJJmuzfE9AwBUo3PRyuk3qplf+ONc+PRRA1Z1Demoyvri1r99n3HzLi82k0jegfX5LH6zfc9RvuRYNvIOZUsDqJ4eY2IIrmJgN7O6ctFHwc3vqyrN0p6elJKGVR3FF7wAtWprVq2l4EbmsSzN8s2IXurWqj5/W7/OrD/D1p54tsXn/Ma9ptWu4frarn7wEqclJSE4SrH1qCNJTk9H24Sl6q6mUmpyEWjrFb40yamD1k5egdlqKYWWtJ3dsn66T0535zwtQIyUZqSn66+vWqgF2vHAZbv5oKUb0CtyfRETQuI5+MG2TZU3T4CWPDMJfxi7Gs1fp95vo1zZTt8HBFV2bY9Kq3TqfiD9bB/eMGik4fqrMq6jDqKxWT7smdbDhmaFBy5z1nN8uEw1rp+H687J15/+q3W5PWJKHhdsO4vda8PtTj5YRVa5Mvft8fDB/u1cOLiU5CW9e2w13jl8R9noBV6uU/UdPYXCnJrjhvGyMW7ADAJARRhlknzaNvIrMfvxHfzTKSEPtGinY8cJluPR/87BOJ/hb5Z7B7ZGWkuR3pxdPgzo2wY4XLsPU1Xu8prfOrI2Ssgrdz7TVAlj7Jhl4Y0RV+bpnLtYd8CPhmys2w90S6EydYQnaafUTu4+cDLqOsdfrt6BxZzBuuSB4E1QrNK6Tjhn/1K8vaFynRsCWZK+P6IZzT2+AJyatDVgBHC+2Du6T7uyHpTtcxSn3Dm6Pc3QqUY24y63NBvZX/tgF/525CS/+oQvqpKd6lVH6ct+yXtu7laU9XTs2q4uX/9jVb/rwLs0xvEtzZI+aHPa6PW/Ve2Y3xLgFO9C7dUO/5nzj/9Yb145ZHNK6fYs7gumRba4zTzAZNVLw4NAOxgvGgW9Borv89t2ftwIArjqneWVA+b+eLXFG4wxLvhMAuLhTE/xn5iYMsaByTy+oA/rNakMdCVGvQtuMs06rizW7qjINV5/bAh8t3BGwqCwQd4V7oH10u/687IAZvHiydXBvk5VReVsWbnfst3RaGgTTuXm9gDkNt0s6NcGMdfvCSk8iGdSxMYZ2bopHh/tXEjer5ypa6d26EUb0aoVZ6/fj/3qGdjdy9mn1sG7PUd2KL7OdeSJ1edfmKK/QzzHHgm/Aczflu6l/68rgIiKWBXbAlTmIZqXeglEXxTUXO/5vfbCvsKq8PDuztuliP0+dm9fFo5d1DPmikChsHdwjMeu+CyGA6TLvULx5bXcc18rZ7Sw9NRmjfVqXPHF5J7TJykDrzNr45YGBaNGgJpKSBFPu9m9FYuTpqzpjRO9WcW0P/saIbnHbtp60lKTKCsdwvDGiGzIsKKKJhG8Fcayfr1M3PdVUx7GvbuuL+ZsPBJwvIrj5/PCbHsdbtQ3ubS2qiNGTlpKEhimJVf5mFc9OIp6VjOGokZIcVlGaE4TRnN2Uy6M8/IXbe9ed69ez2EgsRkIMRc/shl7FkE5TbYO7k519Wj2s3lUY72SQCbF4Ik80XBLDzjgUHgZ3B/ry1r6OKBYa2fd0fLxwZ7yTQWRLDO4OVDMtGTUjGCM9UTx1RWc8NryTqY4yduOucIxnfUOsuI9eugPOSTthcKeEJSJITXZeYAdcnek+vKGn4Vj2TtCsXjruv6Q9ruhqz1YndsXgThQnwcazcRIR8RvegqLPdgOHERGRMQZ3IiIHYnAnInIgBnciIgdicCciciAGdyIiB2JwJyJyIAZ3IiIHEhWt4emMNixSACDcgUMyAQQeq9OZuM/VA/e5eohkn09XShk+ci5uwT0SIpKjlOphvKRzcJ+rB+5z9RCLfWaxDBGRAzG4ExE5kF2D+3vxTkAccJ+rB+5z9RD1fbZlmTsREQVn15w7EREFweBORORAtgvuIjJURDaKyBYRGRXv9IRCRFqKyBwRWS8ia0Xkbm16QxGZKSKbtf8NtOkiIq9r+5orIt091nW9tvxmEbneY/q5IrJa+8zrkiDPqBORZBFZISI/au9bi8hiLf1fiEiaNr2G9n6LNj/bYx0PadM3isgQj+kJd06ISH0RmSgiG7Tj3dfpx1lE/qmd12tEZIKIpDvtOIvIByKyX0TWeEyL+nENtI2glFK2+QOQDGArgDYA0gCsAtAp3ukKIf3NAHTXXtcBsAlAJwAvARilTR8F4EXt9aUApsL1GMo+ABZr0xsC2Kb9b6C9bqDNWwKgr/aZqQCGxXu/tXTdC2A8gB+1918CuEZ7PRrA37XXtwMYrb2+BsAX2utO2vGuAaC1dh4kJ+o5AeAjADdrr9MA1HfycQZwGoDtAGp6HN8bnHacAVwAoDuANR7Ton5cA20jaFrj/SMI8YvtC2C6x/uHADwU73RFsD/fAxgMYCOAZtq0ZgA2aq/fBTDCY/mN2vwRAN71mP6uNq0ZgA0e072Wi+N+tgAwC8BFAH7UTtwDAFJ8jyuA6QD6aq9TtOXE91i7l0vEcwJAXS3Qic90xx5nuIJ7nhawUrTjPMSJxxlANryDe9SPa6BtBPuzW7GM+wRyy9em2Y52G9oNwGIATZRSewBA++9+uGag/Q02PV9nery9BuBBABXa+0YAjiilyrT3nums3DdtfqG2fKjfRTy1AVAA4EOtKGqsiNSGg4+zUmoXgFcA/AZgD1zHbRmcfZzdYnFcA20jILsFd71yRdu15RSRDABfA7hHKXU02KI601QY0+NGRIYD2K+UWuY5WWdRZTDPNvsMV060O4B3lFLdABTBdSsdiO33WSsDvhKuopTmAGoDGKazqJOOs5G47qPdgns+gJYe71sA2B2ntIRFRFLhCuyfKaW+0SbvE5Fm2vxmAPZr0wPtb7DpLXSmx1M/AFeIyA4An8NVNPMagPoikqIt45nOyn3T5tcDcAihfxfxlA8gXym1WHs/Ea5g7+TjfDGA7UqpAqVUKYBvAJwHZx9nt1gc10DbCMhuwX0pgHZaDXwaXBUxk+KcJtO0mu/3AaxXSr3qMWsSAHeN+fVwlcW7p4/Uat37ACjUbsmmA7hERBpoOaZL4CqP3APgmIj00bY10mNdcaGUekgp1UIplQ3X8ZqtlPozgDkArtYW891n93dxtba80qZfo7WyaA2gHVyVTwl3Tiil9gLIE5EztUmDAKyDg48zXMUxfUSklpYm9z479jh7iMVxDbSNwOJZCRNmZcalcLUy2QrgkXinJ8S094frNisXwErt71K4yhpnAdis/W+oLS8A3tL2dTWAHh7rugnAFu3vRo/pPQCs0T7zJnwq9eK8/wNQ1VqmDVw/2i0AvgJQQ5uerr3fos1v4/H5R7T92giP1iGJeE4AOAdAjnasv4OrVYSjjzOApwBs0NL1CVwtXhx1nAFMgKtOoRSunPZfY3FcA20j2B+HHyAiciC7FcsQEZEJDO5ERA7E4E5E5EAM7kREDsTgTkTkQAzuREQOxOBORORA/w8w01sju49Y3gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pd.Series(data).plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fb933aad190>"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxV9Z3/8dfnrtmTm40t7DsqokRQcMGqFa1KbbVitbb+tA7j1JlWx5+2zqi1tlPHx29q3YYyal2q4griVpVNlD1AWMMSQkgCZCf7dpfv748EBjGQG7jJuffm83w88vAuJ+e8T258c3KW7xFjDEoppSKfzeoASimlQkMLXSmlooQWulJKRQktdKWUihJa6EopFSUcVi04PT3dDBs2zKrFK6VURNqwYUOlMSajs/csK/Rhw4aRk5Nj1eKVUioiicj+E72nu1yUUipKaKErpVSU0EJXSqkoYdk+dKVU5PJ6vZSUlNDS0mJ1lKgVExNDVlYWTqcz6O/RQldKdVtJSQmJiYkMGzYMEbE6TtQxxlBVVUVJSQnDhw8P+vt0l4tSqttaWlpIS0vTMu8hIkJaWlq3/wLSQldKnRIt8551Kj9fLXSllIoSWuhKqYhTU1PD888/D8Dy5cu55pprenX5y5cvZ9WqVUefz507l1dfffWU5pWQkBCqWHpQVJ2YP2D4/cd5lNW1kOGsZYC7luT0kcyefpbV0VQfd6TQ77777h5bhs/nw+HovCKXL19OQkIC06ZNA2DOnDk9lqM7tNDVt7y+Zj+rC6r4ak8l5zhW8S+Z8zknfhf4gEYoOTgFz9THiB9ypdVRVR/14IMPsnfvXiZNmoTT6SQ+Pp4bbriBbdu2MXnyZP72t78hImzYsIF7772XhoYG0tPTefnllxkwYAC5ubnMmTOHpqYmRo4cyUsvvYTH42HGjBlMmzaNlStXct1113HbbbcxZ84cioqKAHjqqacYNGgQc+fOxW6387e//Y1nnnmGJUuWkJCQwL/+67+Sn5/PnDlzqKiowG63884779CvXz9mzZrF4cOH8Xq9PP7448yaNSvkP5cuC11EXgKuAcqNMWeeYJoZwFOAE6g0xlwSypCqdy3dWc7SnaX8vxF/5QcJC6hxDGOl+z6K2zJprc7jMvfHxH89k6LMuxjynefBZrc6srLShl/C4dzQztMzCSY/dcK3//jHP7Jt2zZyc3NZvnw5s2bNYvv27QwcOJDp06ezcuVKpk6dyj333MMHH3xARkYGb731Fg899BAvvfQSt912G8888wyXXHIJDz/8ML/97W956qn25dXU1PDll18C8OMf/5hf/epXXHjhhRQVFXHllVeSl5fHnDlzjhY4wJIlS45mu+WWW3jwwQe5/vrraWlpIRAI4HK5WLBgAUlJSVRWVnL++edz3XXXhfzAcjBb6C8DzwKd7iASkRTgeWCmMaZIRDJDF0/1tueW5bNkZxkvjJ3H5e6P2Jn8f8hNe5CAuAFwDoG3q/+FrP2P8iPmsfK1fRSOmcstF4ywOLnqy6ZMmUJWVhYAkyZNorCwkJSUFLZt28YVV1wBgN/vZ8CAAdTW1lJTU8Mll7Rvd/70pz/lxhtvPDqvm2666ejjxYsXs2PHjqPP6+rqqK+vP2GO+vp6Dhw4wPXXXw+0XxwE7Rdi/eY3v2HFihXYbDYOHDhAWVkZ/fv3D9FPoF2XhW6MWSEiw04yyY+B940xRR3Tl4cmmuptucU1PPnZLh4Z9SmXuz9iq+ef2Zp637em65+aSkvKU8zfNZDZrv+mOu83+Ke+id2mp7H1SSfZku4tbrf76GO73Y7P58MYwxlnnMHq1au/MW1tbe1J5xUfH3/0cSAQYPXq1cTGxgaVwxjT6euvv/46FRUVbNiwAafTybBhw3rkKttQnOUyBvCIyHIR2SAit51oQhG5S0RyRCSnoqIiBItWofSnL3ZzkWcvP437C0XxV7HVc+8Jp3XYbPjHPcBi/01cG/MW89/+wwl/mZUKtcTExJNuKQOMHTuWioqKo4Xu9XrZvn07ycnJeDwevvrqKwBee+21o1vrx/vud7/Ls88+e/R5bm7uSZeflJREVlYWCxcuBKC1tZWmpiZqa2vJzMzE6XSybNky9u8/4Qi4pyUUhe4AJgPfA64E/l1ExnQ2oTFmnjEm2xiTnZHR6fjsyiJ//CSPlbsP8R9Zz9Bsz2Rtxn9CF/v3RISK0X9gd+BsrvX+gYVfrzrp9EqFSlpaGtOnT+fMM8/k/vvv73Qal8vFu+++ywMPPMDZZ5/NpEmTjp5q+Morr3D//fczceJEcnNzefjhhzudx9NPP01OTg4TJ05kwoQJzJ07F4Brr72WBQsWMGnSpKP/MBzx2muv8fTTTzNx4kSmTZtGaWkpt9xyCzk5OWRnZ/P6668zbty4EP40/pcEs1XVscvlo84OiorIg0CMMebRjucvAn83xrxzsnlmZ2cbvcFF+Lj8v77kCnmDB/rNY0X/eZTEB38GS2xbETOLriSncRzp1yxl/MDkHkyqwkFeXh7jx4+3OkbU6+znLCIbjDHZnU0fii30D4CLRMQhInHAVCAvBPNVveCNtUU88elOSivL+afMtzgUe3G3yhyg2TUE/1m/46KEjcxf8GdavP4eSquUOpkuC11E3gRWA2NFpERE7hCROSIyB8AYkwf8HdgCrANeMMZs68nQKrRyS2r4ecYHJEgtuWn/95Tm8WHT9ym1jeHn8c/xyte7QpxQKRWMYM5yuTmIaZ4EngxJItWrAsZQcLCMeSMXURx/JYfdp3YVqBEH2/s9wmWHbqFq05/5H/sD/PxiPZVRqd6kY7n0cSXVTVwZ9ymJtnryUv7htOZVFnch+53TmJM+n1W7CkMTUCkVNC30Pm77wRruyPiAMtdkKmMmn/b8dmbeT6qjjhG189lf1RiChEqpYGmh92HGGBIOL2Ooq5Q9nttDMs+qmHM56DqPO9IX8uIK3ZeuVG/Swbn6sC0ltcxK+Ih6PN0+s+VkdqfezYy222nd+ya1zWeRHBv8PRFVZHpjbVFI5/fjqUO6nCYhIYGGhoajz19++WVycnJ49tlnefTRR/ntb3/Lnj17GDVqFAB/+tOfuPfee1m/fj3Z2dk0NDRw3333sXjxYmJiYkhLS+PJJ59k6tSpIV2X3qRb6H3YV1u3cnnSWvYl3kBAXCGb78G4S6mwj+Z2zzu8sz60/6MrFayzzjqL+fPnH33+7rvvMmHChKPP77zzTlJTU9mzZw/bt2/n5ZdfprKy0oqoIaOF3oe5DryHQwIUpdzU9cTdIcKetLsZF7uf3blv4w/okACq933/+9/ngw8+AKCgoIDk5GSOXKG+d+9e1q5dy+OPP47N1l6DI0aM4Hvf+55leUNBC72Pqmlq41w+pzgwijrX6JDPf3/CtTQ7+nOV6z2W79Lx2lToNTc3M2nSpKNfx1++n5SUxODBg9m2bRtvvvnmN0ZR3L59O5MmTcJuj66hn7XQ+6hNeVvIjs+jIO7qHpm/ESe7E27iksSNvPrF8h5ZhurbYmNjyc3NPfr12GOPfWua2bNnM3/+fBYuXHh0SNtopoXeR7UUtO9bPJz+gx5bxr6kmwDhvMACqhvbemw5Sp3Itddey2uvvcaQIUNISko6+voZZ5zB5s2bCQQCFqYLPS30PmpY48cUBsbR5B7eY8tocg5in+sibvR8wSdb9OCo6n2xsbE88cQTPPTQQ994feTIkWRnZ/PII48cHfZ5z549R/e5Ryo9bbEPKjuYx3jXTtYln9q4Ld1RnPoTZnjv5OC2d+CCX/f48pQ1gjnN0CqzZ8/u9PUXXniB++67j1GjRhEXF3f0tMVIFtTwuT1Bh8+1zpbPf83Eyj+Sf/5m1lWk9OiyxPiYufcCtjRkMeLGrxicGtejy1O9Q4fP7R1WDJ+rIkxS5cdsbxnDiGGnNhBXdxhxkJ9wIxclbOKZj5aH/AIUpdT/0kLvY0xLJUPMNvbHXIKtl+4BWpp6IzYxDKqN7P2TSoU7LfQ+pnLPh9jEsFMu7LWt5XrXcPaZM7ki7guqGlp7ZZmq5+k9ZHvWqfx8tdD7mIbCj6j2JWFPO69Xl1uU9AMmxO6jrnRjry5X9YyYmBiqqqq01HuIMYaqqipiYmK69X1dnuUiIi8B1wDlnd1T9JjpzgPWADcZY97tVgrVO0yA9PplfNk0mfSk3j04WZn2fXx1jzOi8QMgxEMNqF6XlZVFSUkJFRUVVkeJWjExMWRlZXXre4I5bfFl4Fng1RNNICJ24Angs24tXfUqU5VDIofZKdMZKL2z//yIVnsa2wIXcInrC2qbWkmOc/fq8lVoOZ1Ohg/vuWsY1KnpcpeLMWYFUN3FZPcA7wE6aEcYq937IQEjVCVeasnyCxOvZ6Crkm25H1qyfKWi3WnvQxeRQcD1wNwgpr1LRHJEJEf/VOt9vgOfsKV5FJkZ3fszLlSaM79HUyAGU/iWJctXKtqF4qDoU8ADxhh/VxMaY+YZY7KNMdlHhrFUvaS1mtSWXNa1TCUz0ZrdHcYexybfNMb7FuP1ei3JoFQ0C0WhZwPzRaQQuAF4XkS+H4L5qlAqX4GNALXJM5Be3n9+rP3xV5PmqCF/+yeWZVAqWp12oRtjhhtjhhljhgHvAncbYxaedjIVUvVFX9AccNNv+HRLc7RlXkVTwE3rXt3tolSodVnoIvImsBoYKyIlInKHiMwRkTk9H0+Fir90ORsax3HeyIGW5nDFJJDrm87Q5s8g0OVeOqVUN3R52qIx5uZgZ2aM+dlppVE9o7WapJY8NjT/hLT9h7FZuMsFYLPtcqbZllJXvJSkoVdYmkWpaKJXivYF5SuwiaHYOcXyMgdoTr+S5oCbqh2vWx1Fqaiihd4HNBQvpjngxpsy2eooAGSmpvFlw3mkHf4YTHTdMUYpK2mh9wG+Q8vY0DiOwekeq6MAYBNhm1xKEpWYynVWx1EqamihR7vWapJa88hpnkj/5O4N9NOTDidfjs/YqN6tw/4oFSpa6NGufAU2wmf/+RGD+g1kbcOZ2A4ssjqKUlFDCz3Kbdn4AS0BF23J4bH//IiUOBebzAw8vj1Qt8fqOEpFBS30KJfSlMPmpjFhs//8WP4B1wLgLVpgcRKlooMWejTzNTMwkEdu8zgGJMdaneZbJo6byPbmETQVvG91FKWighZ6NKvegEP8FNkmYe+l+4d2x/nD01hafz6JDeuhRUdeVup0aaFHsbayrwGoi8+2OEnnYl12DiVeiY0AHPjI6jhKRTwt9CjWUPIV+1oH4EkbZHWUTr2xtohq1wRK2jJpKdT96EqdLi30aGUMMbVr2dg4niGe3r1/aHeM7pfIF3VTcZQvBl+T1XGUimha6NGqoYC4QBW7/GcR5w7m1rHW6J8Uw8rmaThMCxz63Oo4SkU0LfQoZSpWAVARE17nnx9PRGhInkatPwF/iV5kpNTp0EKPUnUlK6j3x2JLnmB1lC6N6Z/Kivpz8JXoYF1KnQ4t9CgVqFjFpqZxDE5LtDpKl0ZmJvB141Tc3nI4vMnqOEpFrGDuWPSSiJSLyLYTvH+LiGzp+FolImeHPqbqFm8Dya072dE2gXSLbgjdHU67DW+/7xIwgtHTF5U6ZcFsob8MzDzJ+/uAS4wxE4HfAfNCkEudjsMbsRGgIf7csBqQ62TOnzCO3KYxNBd+aHUUpSJWl4VujFkBVJ/k/VXGmMMdT9cAWSHKpk6Rt3wtADH9p1qcJHg1TV6W1Z9HXP0GaC6zOo5SESnU+9DvAD490ZsicpeI5IhITkVFRYgXrY6oP7iakrZMRg0ZaXWUoCW4HeTJRe1PDp3wV0gpdRIhK3QRuZT2Qn/gRNMYY+YZY7KNMdkZGRmhWrQ6jrMmh81Nozl7cLLVUbrFkX4upd5Umvbp6YtKnYqQFLqITAReAGYZY6pCMU91iloqSfQVs8s7jqV5kTXg1RkDk1lWl42j/AsIeK2Oo1TEOe1CF5EhwPvAT4wxu08/kjot1esBKHWchUTIAdEj0hLcbPRPw2UaoOJrq+MoFXG6vCZcRN4EZgDpIlICPAI4AYwxc4GHgTTg+Y4C8RljwnN4vz6gtXwtTiM0JkTm2aP1yTNoDTgIFC4itt+lVsdRKqJ0WejGmJu7eP9O4M6QJVKnpfHgaqpas0hPi8xjFKMG9mdt4VlMKvqQ2Kl/sjqOUhFFrxSNJsYQU7eBLU2jyfKE3x2KgjEgOYYN3mkkefdCQ4HVcZSKKFro0aSpmLhAFbv944kP4xEWT0ZEcAy+BoDW/Xq2i1LdoYUeTarWAVDhmmhxkNNT5xjK3pYsavcutDqKUhFFCz2KNJetw2vstCWcaXWU0zI0LY6vmqbgaVgF3gar4ygVMbTQo0hzWQ57WobQ35NidZTTYhOh0DUDJ168B7+wOo5SEUMLPVoYQ0z9ZrY1j2RgSmQeED2WLfNC6v2xlO/Se40qFSwt9GjRfIA4U02pfQIxTrvVaU7bsEwPq5smE1/5ORhjdRylIoIWepQw1RvbH3jOsTZIiDjsNsoTLyOFMvzVm62Oo1RE0EKPEjmblhAwQmFglNVRQiZzzPcBOLjjXYuTKBUZtNCjRELjFgpaB5HpSbc6SsiUtCSztWkUzfv1phdKBUMLPUpk+vPY0TKSAckxVkcJmRinnY3+aYyUbdCqg3gq1RUt9GjQUkmalLHfjMVpj66P9GDcd7BLgNoC3UpXqivR9X9/H3XkgGi1K7IvKOqMK3MKlb5kavJ1GACluqKFHgVqDrTfQ7Q1IbIv+e9M/+Q4VjWeR1rdUgj4rY6jVFjTQo8CLeXrKW7rR4qnn9VRQk5E2OO8mASpxVex2uo4SoU1LfQo4K7fwo7mEfRLip4Dosdq8HwHn7HpVaNKdaHLQheRl0SkXES2neB9EZGnRSRfRLaIyLmhj6lOyFtHamA/hYHoOyB6RFbmQDY0TcB+6FOroygV1oJpgJeBmSd5/ypgdMfXXcB/n34sFbTD7VdRljvOsDhIz4l12dnkm0Y/fx4LVq6zOo5SYavLQjfGrACqTzLJLOBV024NkCIiA0IVUJ1c/aH2A6JNCWdZnKRnlca33180tXaxxUmUCl+h+Bt9EFB8zPOSjtdUL2g4tJ5yr4eElCFWR+lR8ZmTKGnL0EJX6iRCUejSyWudDo8nIneJSI6I5FRUVIRg0cpZ1z5kbjRdIdqZASmxfN04lVGBNeBvtTqOUmEpFIVeAgw+5nkWcLCzCY0x84wx2caY7IyMyLwrfVjxNePx5lPgHx0VQ+aejE2Evc6LiJVm/GVfWh1HqbAUikJfBNzWcbbL+UCtMeZQCOarulKzFbv4KYviA6LHak69hJaAi8rdevqiUp3p8tbwIvImMANIF5ES4BHACWCMmQt8AlwN5ANNwO09FVZ9U3N5DrFAQ2x0HxA9Yli/DFZvP4tJZX+3OopSYanLQjfG3NzF+wb4p5AlUkGrO7iWVl8C7uQRVkfpFXEuB7n+6Vzqfxrq9kDSaKsjKRVWovNKlD5CanLZ3jyCAVFwD9FglSa0n77YWKiDdSl1PC30SBXw4mnbyR7vaBJjnFan6TVp/caxp2Uwjfu00JU6nhZ6pKrNw0kbpfbxVifpVQNTYvmqcQqehtXgbbA6jlJhRQs9QrVVbgCgLjb6xkA/GZsIBa5LcIqXwCG9yEipY2mhR6iaA+tpDriRxLFWR+l1vtRp1PnjqM5faHUUpcKKFnqEMtUb2dk8jAGeBKuj9LqR/Tx8XX8O7vK/g+n0omSl+iQt9EhkDEktO9jjHYUnru8cED0i3u1gc2AaiYEyqNlidRylwoYWeiRqLCSWeupizkCks6F0ol95x+mLTXr6olJHaaFHIH9V+02hbWnnWJzEOv37DWdz02iatdCVOkoLPQIdPrAev7GRnjXZ6iiWGeSJZU3LVDzNG6BFR+5UCrTQI5K3cgN7W7MYn9Xf6iiWsYlQm3Y1NgyBkg+tjqNUWNBCj0DOuq3sbBnBmoKT3Ugq+o0eeyElbRnU733P6ihKhQUt9EjTWkW6rYwSxmC39c0DokdcPCaTJXVTiataCr4mq+MoZTkt9AgTqNoEQLWrb4yBfjJpCW72ui7DSQuU6lWjSmmhR5iaQ+sBaInvG2OgdyV1xBXU+eNpLdSbXiilhR5hmstyONiWTpJngNVRwsLF4wayrC4bDnwIAb/VcZSylBZ6hHHVbyGvZTj9kqL7ptDBmpSVwpq2C3H7q6BytdVxlLJUUIUuIjNFZJeI5IvIg528P0RElonIJhHZIiJXhz6qwtdMqm8f+wNjcNr13+I31hYxf30xB2Jn4DUOfEW620X1bV22gojYgeeAq4AJwM0iMuG4yf4NeNsYcw4wG3g+1EEVmI6bQlc4j//x923DBwxgdcNZtO1fqIN1qT4tmM28KUC+MabAGNMGzAdmHTeNAZI6HicDB0MXUR1R33FAtCm+b42B3pWRGfEsa7iAuNYCPvp6qdVxlLJMMIU+CCg+5nlJx2vHehS4VURKgE+AezqbkYjcJSI5IpJTUaGXa3dX3aH11PnjcCWNsjpKWHHYbeyLuQyAQY2fW5xGKesEU+idXb1y/N+1NwMvG2OygKuB10TkW/M2xswzxmQbY7IzMjK6n7aPs9duIa95BAP70E2hg9Wv/yi2NI0is/Yzq6MoZZlgCr0EGHzM8yy+vUvlDuBtAGPMaiAGSA9FQNUh4Ce1LY9ixuJ22q1OE3bG9ktkSf35DPBtgeZSq+MoZYlgCn09MFpEhouIi/aDnsePWVoEXAYgIuNpL3TdpxJKDfm4pYXmBN1/3pkYp51dju9gE4Mp0SF1Vd/UZaEbY3zAL4DPgDzaz2bZLiKPich1HZPdB/xcRDYDbwI/M0ZPNwilxtL2A6LujL47ZG5XEvqdQ3FbP+ry37U6ilKWcAQzkTHmE9oPdh772sPHPN4BTA9tNHWswwfW4Qw4GDD4XIpr9YrIzkwYmMznGy/gp66Poa0WXMlWR1KqV+nVKRHCVOeyp3UIE7L00MSJxDjtbLddgQMvAd3tovogLfRIYAwprdvZ5x9NWoLb6jRhzdnvAg61pVG7+02royjV67TQI0FLKYlUUx+rB0S7Mm5ACp/XTyexegl466yOo1Sv0kKPAC3lOQDsD4zmjbVFFqcJb26nnfKU63DQhq9Yb02n+hYt9AhQVbwOAG/iRIuTRIZJ515NmTeV6rw3rI6iVK/SQo8A3qqN7G/tT5pHr64NRmldG4sbLiSlZgl4G6yOo1Sv0UKPAAmN29jZOpLkWKfVUSKC3Sbsi52JS1qpL1hodRyleo0Werjz1pNOEQdlLCJ9+6bQ3RGfNYMKbwqVO163OopSvUYLPcx5K9tvCl2lN4Xuln7J8az3z2BA0zLwNVodR6leoYUe5iqK1gDQlqg3he4ux7AfESOt7N+qQwGovkELPcwV7/2aCm8KSZ5hVkeJOOdfcD0VPg/1u16zOopSvUILPcz1929nR+toPPEuq6NEnI+2lLPOfzljfF9SX6eDf6rop4UeznxNZNkKKWa8HhA9RVUZN+Cy+di25q9WR1Gqx2mhh7HW8g3YJUC1Wy/5P1X2tCkUewcRe/AtdERnFe200MNY2f6VADQnnmtxksglNhvbnN9jonMT2/J3WB1HqR6lhR7GWstzqPCmkOwZanWUiNY04EfYxFC44SWroyjVo4IqdBGZKSK7RCRfRB48wTQ/EpEdIrJdRHQQjRBIaMxlR+sYkuL0gOjpaI0bTb5/HCMaP+DFr/ZZHUepHtNloYuIHXgOuAqYANwsIhOOm2Y08GtgujHmDOCXPZC1b/E1kmn2USLjrU4SFYoSr+eM2L0cLN5odRSlekwwW+hTgHxjTIExpg2YD8w6bpqfA88ZYw4DGGPKQxuz72ksy+k4IKoXFIVCTcYP8Bsbg2re14OjKmoFU+iDgOJjnpd0vHasMcAYEVkpImtEZGZnMxKRu0QkR0RyKir0vOCTKStsPyDalqQHREOhxZFJvkzh8tgl5BRWWx1HqR4RTKF3dgL08Zs4DmA0MAO4GXhBRFK+9U3GzDPGZBtjsjMydCjYk/GW51Du9ZCUqgdEQ6U87QaGuMtYs1bvN6qiUzCFXgIMPuZ5FnCwk2k+MMZ4jTH7gF20F7w6RUlNm8n3jSHe7bA6StQ4lHQ1zSaOgdVvUt/itTqOUiEXTKGvB0aLyHARcQGzgeM3cRYClwKISDrtu2AKQhm0T/E1kkkhh2PPtjpJVPHZ4tntvpqZiV/x+wXrrI6jVMh1WejGGB/wC+AzIA942xizXUQeE5HrOib7DKgSkR3AMuB+Y0xVT4WOdrUH1mGXAPbUyVZHiTql6TcTb2/BU7XA6ihKhVxQf88bYz4BPjnutYePeWyAezu+1Gkq37+SZCB96DSq9cb1IVUVM5lDZjiXuz8mv/wRRmUmWh1JqZDRK0XDkK8yhzJvKmOHj7U6SvQRoTD5JibH72Tp2i+tTqNUSGmhh6Hk5i0U+MeSGKP3EO0JBz034jN2YotfwesPWB1HqZDRQg833gb6U0ihGc8ba4usThOVWh3p7HbM4KqEL/gy74DVcZQKGS30MFNdshqbGGpiJ1odJaqVpt9CuqOW3ZvmWx1FqZDRQg8zlYXt+3Xbks6zOEl0K4u/hHrJYHzT25TXt1gdR6mQ0EIPN1XrKGwdQErqQKuTRDUjDrzDbueShA0sWb/a6jhKhYQWejgxhoyWTezyTcBp14+mp6WefQ8GG/a9f9EBu1RU0NYII/6GYjy2Sg7Y9QrRXhGXxcHEmXzX/TFr9xR3Pb1SYU4LPYyU7lsOQH28jrDYG95YW8S2hNtIcTSwa81frI6j1GnTQg8j9SVf0RpwYk+bZHWUPuNw3FSKA6OY4n2DgvJ6q+ModVq00MOIuzaHna2j8CQkWB2l7xChIPV2xscWsvSr96xOo9Rp0UIPFwEvAwM7KHVOQqSzIehVTyn1/JAmksiqeJGapjar4yh1yrTQw0RD2Ubc0kYgbYrVUfocvy2W7bE3cnniKh5/53Or4yh1yrTQw0TZ3uUApA292NogfdTBjNsRYEztazS2+qyOo9Qp0UIPE5WFyyj3etha/a0796le0OgczB7XFcz2fMzbq7dbHZ6wBC4AABEmSURBVEepU6KFHiaGm01sbj2LWL3lnGUKM/+JJHsjtVufp8XrtzqOUt0WVKGLyEwR2SUi+SLy4Emmu0FEjIhkhy5i9PPX7yfTXk6xQ88/t1J1zNnss0/hpsT3eG/9XqvjKNVtXRa6iNiB54CrgAnAzSIyoZPpEoF/BtaGOmS0O7in/UBcXcL5FidR+zLuZoCril1r5+lY6SriBLOFPgXIN8YUGGPagPnArE6m+x3wn4AOXddNDcXLqPfHEpNxjtVR+rzSuBmUyhh+kvw2CzfqcAAqsgRT6IOAY3+zSzpeO0pEzgEGG2M+OtmMROQuEckRkZyKiopuh41WSQ1r2dp6BklxMVZHUSLkZ/yC0THF7Fj7V/wBHbRLRY5gCr2zq1yO/paLiA34E3BfVzMyxswzxmQbY7IzMjKCTxnFAs2VDJICCm26/zxcFCdcQ7kMZ3bcX/l0q97RSEWOYAq9BBh8zPMs4OAxzxOBM4HlIlIInA8s0gOjwSnZsxiAmvipFidRRxixszvjl4yNKWLzqhcJ6Fa6ihDBFPp6YLSIDBcRFzAbWHTkTWNMrTEm3RgzzBgzDFgDXGeMyemRxFGmtmgZbQEHtgy9QjScFCdcS71rFD90/ZW/bzvY9TcoFQa6LHRjjA/4BfAZkAe8bYzZLiKPich1PR0w2sXXrmandyxJ8UlWR1HHMGJnc8o9jIvdz9dL5+q+dBURgjoP3RjziTFmjDFmpDHm9x2vPWyMWdTJtDN06zw4gdZ6BpNHZYxunYej4oRrKZWR3Jn8Ep9sKbI6jlJd0itFLVS48+84xcdem94QOhwZsZOX+RtGuA+yd/VT+PS8dBXmtNAtdLjg77QGHNj76YBc4epQ/GXsk3O5Je5lPt60x+o4Sp2UFrqFPPUr2NZ2JjGxiVZHUSciwq4BD5HhrKEq54+6la7Cmha6RaqqihnhyGe/c5rVUVQXqmOz2WG/jJvi5/P7dxdbHUepE9JCt8jere3Hkxs9M6wNooKyZ+AjOGx+Jh/+Tx3jRYUtLXSL+A4uod4fjy1tstVRVBAaXUNZ4/oZ1yQtYclXC6yOo1SntNAt4PP5GepdxX7HFMSm459HikODfkWFP52h+35NY0ur1XGU+hYtdAts272JQc4yZMDlVkdR3RCwx7Mq8UHGu/ew/ov/sDqOUt+ihW6Bkrz2/efDz9ALbSNNY/8b2BGYSnbNkxyuLLA6jlLfoIXeywIBQ8rhz6kwA4lLP8PqOKq7RFiT/nvs4qd08R1gdEgAFT600HvZ5sISzovZyE77pbyxTm+gEIliPKN53/tzxgeWs2/TK1bHUeooLfRetiv3A9w2Lw3pV1kdRZ0G76h/ZmfraDw7foWv8ZDVcZQCtNB7VSBgiK/8mIZAAjV6/9CI5nK5qTxrHrE0cuCzW3XXiwoLWui9KLe4mvPda9hjvxAjTqvjqNM0ffIlLDC/YGjLUorW/9nqOEppofemTZs+J8NZQ41nptVRVAi8ua6YhqF3s7ZpEum7f83h0q1WR1J9nBZ6L/EHDLaDi/AbO5VJ37E6jgqROLeLjQOfojXgomHx9/G2NlgdSfVhQRW6iMwUkV0iki8iD3by/r0iskNEtojIEhEZGvqoke3xj3cwPWYVhfZz8dqTrY6jQig5bThbhj7HINnHlgU/xgR0rBdljS4LXUTswHPAVcAE4GYRmXDcZJuAbGPMROBd4D9DHTTSHT6Qw5iYIipSrrU6iuoBB2IuZpH/DiYHPmT5p3+wOo7qo4LZQp8C5BtjCowxbcB8YNaxExhjlhljmjqergGyQhszsjW3+ZnQ9hE+Y+dAkhZ6tGoa/Ws2eqdwYc2j/PdbL1gdR/VBwRT6IODYK2BKOl47kTuAT08nVLT5YschrklaToHzQlrtqVbHUT3F5mDXyP/hoG8Qs1vvZeP2tVYnUn1MMIUunbzW6Um3InIrkA08eYL37xKRHBHJqaioCD5lhMvb/BEDXZWUpf7Q6iiqhxlnCmuHvIKI4Fn/Q/KLC62OpPqQYAq9BBh8zPMs4ODxE4nI5cBDwHXGmE7HFjXGzDPGZBtjsjMyMk4lb8QprW1haP0Cmk0cB+O/a3Uc1Qt8cSNYmjmXAY4KvItnUlHddzZelLWCKfT1wGgRGS4iLmA2sOjYCUTkHOAvtJd5eehjRq4FOflcnbySwtjv4rfFWh1H9ZKWlOmUnfUKo117OLBoJo1Nejqj6nldFroxxgf8AvgMyAPeNsZsF5HHROTI+K9PAgnAOyKSKyKLTjC7PiUQMBTkzifJ3sghj+5u6WuGnj2b3cOfZpJrI/nvXYmvranrb1LqNIixaAyK7Oxsk5OTY8mye8vK/EpYchnjEyv5bMTXGLFbHUlZwFbwAj8KPM5e+1RG/XAx4oy3OpKKYCKywRiT3dl7eqVoD1q67iumJ25hX8otWuZ9WGDEnbzo+zdG+tey++2Loa3G6kgqSmmh95CyuhaGVr2MzzjYn3yT1XGUxWLH3sGfGh9muNlCxfvnEagvsjqSikJa6D3kjRWbuMHzGXtir6XF0TfO6FEnZhMh86zbeaLlv3B7D1K3KJum0nVWx1JRRgu9B9Q2e7HtnUecrZWC9DlWx1Fhwm4TRpx5HV8MfJdGL9gXX0T+muesjqWiiBZ6D3hr9Q5u9bxPvedyat3jrI6jwoiI0JpwBu+kL2BzywRGFfyCz/56PYeqq62OpqKAFnqItXj9tO14mjRHHYnnPW51HBWmMjMGkzfuHT7y3sqV7oU0fjCJlz5cQG2z1+poKoJpoYfYq8s3cWvS2xxOuQLSp1odR4Uxp9NF3bjf84Hnr2S4m7m19kfMf/F2/roij1af3+p4KgJpoYdQcXUT7p2Pk2hvYnXCvbyxVs9kUF1rTP0OXwxbTH7MTP4h9XUu3XsZv3rqCRZuOkAgoPcqVcHTQg+hvyz6kB+nfsTO+JuocR8/ZLxSJ9Zm95A7+DmWDHidxBgXzw98COfqm/iH59/k6z2VVsdTEUILPUQ+2VLCD72P4LcnsTPzfqvjqAhVFnchXwz7nM2eX3FFygb+2/MTDnz2Y3750kK2H6y1Op4Kc1roIbCrtJ4dS/6dc+J3sT7tUVrtaVZHUhEsYIthe+ov+Wjo1+xO+ik/8CznSfeNbHn/Jp6Y/x4lh3VMGNU5LfTTdLixjT/Pf4lfZr5KfsxVFCfO6vqblApCiyOD3MxH+XjYCvITfsQPPct5IHADRe9M5433n2NfuW6xq2/SwblOQ22zl7uef4fnMv4JmyuJxUM/wWdLtDqWilJufzWDqv7G6NpXSbNVUO71kMNM7CNu5bxzLiM1wW11RNULTjY4lxb6KapqaOVXL3/Mb+PvoZ+7gSWDF1LvGml1LNUHiPGSWvMZGZVvc5asxCk+9rQMZgsX05Qxk6xRl5E9PJ3EGKfVUVUP0EIPsZ2ldfzhzQ/4necB+rvq+HLQq1TGdPrzVapHOX3VpFQuZFD9x4yRTTjET7Uvia8bJlHinIK9/wxGj5pM9vA0krTgo4IWeogYY5i/vpivl/0P/zHwKWLdbpb2e4GqmHOtjqYUTn8tGQ1fklb7GVlta/BI++mOFd4UtjSPpsxxJt7kc4jrdx6DB41i3IAkUuJcFqdW3aWFHgIb9h/mxU8/53uBZ/heykoqnGewasBcGp1DrI6m1LcZQ6K3kLSmVSTUrWFgYAepvgJsEgDaS76gNYtSMwRv/BgS0sfTf/BERg+bQHx8ksXh1cmcdqGLyEzgz4AdeMEY88fj3ncDrwKTgSrgJmNM4cnmGQmF3tjqY/HWfezKXcB5vre5JHEjXnGz03M3eZ5/JCC6daMihz3QRErrDuIaNhPftJUkbwEZppBk2zfPlqkNJFEnmTQ5BtLqHEAgZhD2uExccRntX/GZxCZmEpfQD7c7FhGxaI36ppMVuiOIb7YDzwFXACXAehFZZIzZccxkdwCHjTGjRGQ28AQQPnd1MAZMgEAgQMAE8AcCGGNo9flpavPS3FRPW3MVrY3VHD58iJaafLw1u+jn28ZVsbuYFeejxqSzOeUe9qXcSoujn9VrpFS3+W1xVMVmUxX7zS5w+WtwNO7GX7sHmkqI8R0iMXCIdNsBBjlySWupgxPcZKnBH0ujiafFxNJKHF6Jw2eLw2dLIGCPxzgSCNjjwObu+HIhdjfY3YjNjTjciN2NzR6DzeHCZndjtzuw2x2IzQHYEJsDEXv7c5sNEQc2mw3Egc1mB5sdW8f7YrNhExsicvTLZrMhCLajr9kQmyCAiA2bTRDkmO+xIdI+hn37NETMP1pdFjowBcg3xhQAiMh8YBZwbKHPAh7tePwu8KyIiOmJ/TnFC2D1T9pLmo4vY4DA0dcCpv01m3xz8baOryMrHQMkn2AxLU43B5wjyYn5KS2pF1EWdyFG9KCSij5t9hTakqZA0hQAWoFa2rfe/AFDS2sjgeYqaKvC4a3G6a/G5T9MjP8wMaYGV6ABJ024TBNumnD7K0iW/cRKM3G2ZuJsLdg7dvVEmoAR2ttFOtqmvdiN+d+CN0jH60KMs7NLezr5x2D8fTDxsZDnDabQBwHFxzwvAY4fRvDoNMYYn4jUAmnANwahEJG7gLs6njaIyK5TCQ2kHz/v0Gul/d+sHcCLPbuo4PXCeocdXee+IUzX2Rz331D5HfC7U13noSd6I5hC7+xvjePXLphpMMbMA+YFscyTBxLJOdE+pGjWF9db17lv0HUOjWAu/S8BBh/zPAs4eKJpRMRB+54MvQWLUkr1omAKfT0wWkSGi4gLmA0sOm6aRcBPOx7fACztkf3nSimlTqjLXS4d+8R/AXxG+2mLLxljtovIY0COMWYR7TuZXxORfNq3zGf3ZGhCsNsmQvXF9dZ17ht0nUPAsguLlFJKhZYOn6uUUlFCC10ppaJEWBe6iMwUkV0iki8iD3byvltE3up4f62IDOv9lKEVxDr/TEQqRCS34+tOK3KGkoi8JCLlIrLtBO+LiDzd8TPZIiIRPxpaEOs8Q0Rqj/mcH+7tjKEmIoNFZJmI5InIdhH5l06miarPOsh1Dt1nbYwJyy/aD8DuBUYALmAzMOG4ae4G5nY8ng28ZXXuXljnnwHPWp01xOt9MXAusO0E718NfEr79Q7nA2utztwL6zwD+MjqnCFe5wHAuR2PE4Hdnfx+R9VnHeQ6h+yzDuct9KNDDhhj2oAjQw4caxbwSsfjd4HLJFIGXehcMOscdYwxKzj5dQuzgFdNuzVAiogM6J10PSOIdY46xphDxpiNHY/rgTzarzI/VlR91kGuc8iEc6F3NuTA8T+Ibww5QPsQFJF8h+Zg1hnghx1/jr4rIoM7eT/aBPtziTYXiMhmEflURM6wOkwodewePQdYe9xbUftZn2SdIUSfdTgXesiGHIggwazPh8AwY8xEYDH/+xdKNIu2zzkYG4GhxpizgWeAhRbnCRkRSQDeA35pjKk7/u1OviXiP+su1jlkn3U4F3pfHHKgy3U2xlQZY1o7nv4P7WPQR7tgfheiijGmzhjT0PH4E8ApIukWxzptIuKkvdheN8a838kkUfdZd7XOofysw7nQ++KQA12u83H7E6+jfZ9ctFsE3NZxBsT5QK0x5pDVoXqSiPQ/cjxIRKbQ/v9qlbWpTk/H+rwI5Blj/usEk0XVZx3MOofysw5mtEVLmPAccqBHBbnO/ywi1wE+2tf5Z5YFDhEReZP2I/3pIlICPAI4AYwxc4FPaD/7IR9oAm63JmnoBLHONwD/KCI+oBmYHeEbKwDTgZ8AW0Ukt+O13wBDIGo/62DWOWSftV76r5RSUSKcd7kopZTqBi10pZSKElroSikVJbTQlVIqSmihK6VUlNBCV0qpKKGFrpRSUeL/A5Ax7hCnZC8TAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "sns.distplot(list(data)[B:],bins=100, label='HMC', ax=ax)\n",
    "arr = np.arange(0, 2.5, 0.01)\n",
    "sns.lineplot(\n",
    "    data=pd.DataFrame(\n",
    "        data=stats.gamma.pdf(arr, alpha, scale=1/_lambda),\n",
    "        index=arr,\n",
    "        columns=['theoretical']\n",
    "    ),\n",
    "    ax=ax,\n",
    "    palette=['orange']\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}