aboutsummaryrefslogtreecommitdiff
path: root/statistic/try_HMC.ipynb
blob: 2609f9d34e5f859fa9c621611b459d242fa484ec (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
226
227
228
229
230
231
232
233
234
235
236
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 12,
   "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": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 11\n",
    "_lambda = 13\n",
    "B = 1\n",
    "N = 100000 + B\n",
    "epsilon = 0.01\n",
    "L = 100\n",
    "initial = np.full(1, 2.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def leapfrog(theta, p, h_prime):\n",
    "    half_next_p = p - (epsilon * h_prime(theta)) / 2\n",
    "    next_theta = theta + epsilon * half_next_p\n",
    "    next_p = half_next_p - (epsilon * h_prime(next_theta)) / 2\n",
    "    return (next_p, next_theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hamiltonian(theta, p, h):\n",
    "    return (\n",
    "        h(theta)\n",
    "        + (\n",
    "            np.sum(p ** 2) # 各次元がもつ運動エネルギーの足し合わせ\n",
    "             / 2\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def r(theta1, p1, theta2, p2, h):\n",
    "    return (\n",
    "        np.exp(hamiltonian(theta1, p1, h) - hamiltonian(theta2, p2, h))\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hmc(h, h_prime):\n",
    "    data = np.empty((N, len(initial)))\n",
    "    data[0] = initial\n",
    "    accept_count = 1 # data[0] は受容とみなす\n",
    "    for i in range(1, 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(0, 1, len(initial))\n",
    "        start_theta = theta = data[i - 1]\n",
    "\n",
    "        for _ in range(L):\n",
    "            p, theta = leapfrog(theta, p, h_prime)\n",
    "\n",
    "        if np.random.rand() < r(start_theta, start_p, theta, p, h):\n",
    "            data[i] = theta\n",
    "            accept_count = accept_count + 1\n",
    "        else:\n",
    "            data[i] = start_theta\n",
    "            \n",
    "    print(\"\\nacceptance ratio: %s \" % str(accept_count / N))\n",
    "    \n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100000 / 100001 100.0 %  \n",
      "acceptance ratio: 0.999810001899981 \n"
     ]
    }
   ],
   "source": [
    "def h(x):\n",
    "    \"\"\"\n",
    "    目標分布(事後確率)関数の対数関数のマイナス。必ずスカラー値を返却させる\n",
    "    \"\"\"\n",
    "    return (_lambda * x) - ((alpha - 1) * np.log(x))\n",
    "\n",
    "def h_prime(x):\n",
    "    \"\"\"\n",
    "    目標分布関数の対数関数の微分のマイナス\n",
    "    \"\"\"\n",
    "    return _lambda - ((alpha - 1) / x)\n",
    "\n",
    "data = hmc(\n",
    "    h,\n",
    "    h_prime\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fbac21509d0>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3wUZf4H8M83jVBDSagBQkdAQTqiiCAicOrd2dATz3ao2MvvDvt5nifWs55g74JdlCYiVTGQUKRLgAChBgKhhJD2/P7Y2bB9Zndny0w+79eLF9nZycyzmdnvPP0RpRSIiMheEmKdACIiMh+DOxGRDTG4ExHZEIM7EZENMbgTEdlQUqxOnJ6errKysmJ1eiIiS8rNzT2glMrQ2y9mwT0rKws5OTmxOj0RkSWJyHYj+7FahojIhhjciYhsiMGdiMiGYlbnTkQUK+Xl5SgoKEBpaWmsk+JXamoqMjMzkZycHNLvM7gTUY1TUFCA+vXrIysrCyIS6+R4UUrh4MGDKCgoQLt27UI6hm61jIi0FpH5IrJBRNaJyF0+9hkqIsUiskr792hIqSEiioLS0lI0adIkLgM7AIgImjRpElbJwkjOvQLAfUqpFSJSH0CuiMxVSq332G+xUuoPIaeEiCiK4jWwO4WbPt2cu1Jqj1JqhfbzUQAbALQK66xh2LT3KF74YRMOHDsZqyQQEcW9oHrLiEgWgDMBZPt4e5CIrBaRWSLS3c/vjxeRHBHJKSwsDDqxAJC3/xhe/ikPRcfLQvp9IqJ4MHv2bHTp0gUdO3bEpEmTTD++4eAuIvUAfAngbqXUEY+3VwBoq5TqCeAVAN/4OoZS6g2lVF+lVN+MDN3Rs0REtlRZWYnbbrsNs2bNwvr16/Hpp59i/XrPmu7wGAruIpIMR2D/WCn1lef7SqkjSqlj2s8zASSLSLqpKSUisolly5ahY8eOaN++PVJSUjB27Fh8++23pp5Dt0FVHLX6bwPYoJR6wc8+zQHsU0opEekPx0PjoKkpJSKKgMe/W4f1uz0rI8LTrWUDPHaRz9ppAMCuXbvQunXr6teZmZnIzvZV2x06I71lBgMYB2CNiKzStj0IoA0AKKUmA7gMwK0iUgHgBICxiouz1hhVVQoKQGJCfPc+IIoXvsKj2b13dIO7UmoJgIBnVUq9CuBVsxJF1jLivwuRf7AEW/4zOtZJIQpaoBx2pGRmZmLnzp3VrwsKCtCyZUtTz2HZuWVYLogfWwqPo7KKF4TIqH79+mHz5s3Ytm0bysrKMHXqVFx88cWmnsNy0w/E+bgDIiJdSUlJePXVVzFy5EhUVlbihhtuQPfu5pYgLBfciYjsYPTo0Rg9OnJVmZatliEiIv8Y3ImIbIjBnYhqpHjvrR1u+iwb3B09q4mIgpeamoqDBw/GbYB3zueempoa8jEs16DKzjJEFK7MzEwUFBQg1AkMo8G5ElOoLBfciYjClZycHPIKR1Zh2WoZIiLyj8GdKIoKj57EnHV7Y50MqgEY3IlCcKS0HPuPBL++5TVvZePmD3NRWl4ZgVQRncLgThSCIc/MR///zAv693YUlQAAquK0lwbZh2WDO78bFEuHS8pjnQSigCwX3DlxGBGRPssFdyIi0sfgTkRkQwzuREQmKDx6EtNX7451MqoxuMehqiqFQU/Nwzcrd8U6KURk0HXvLsOdn67E4ZKyWCcFgIWDu517y5RWVGJPcSke+GpNrJNCFJbCoydxtLRm9CzaU+wY9xAvS05aMLizuwyRVfR78kec99zCWCejRrJgcCeieLVsWxEOHjvptu2Ax2sr2VtcGjfVLMFicCci01wxZSkuff2XWCfDNAOfmodBT/0U62SEhMGdiEyVf7AkIsetqlIxWVzjhEXnAWJwJyJL6PDQTIx6aXGsk2EZDO4+ZE2cgf/7fHWsk0FELpQCNu49GutkWIZlg3uk11D9PLcgoscnayg6XoayiqpYJyPqNu09WiM/t51YLrgHM3HY3VNX4vp3l0UuMeTloleW4OuV9nkw9n5iLm79KDfWyQhKZZUKq2/5nuITGPniIvzzu3UmpoqizXLBPRjfrNqN+ZvidwFcf6w8QGvNrmLcM81eVVrzNu6PynmUUigpqwj7OP+esR6n//MHnCgLrSHQOZ3xiu2Hwk4LxY6tg7vVRWt642MnK1BeySJ4rL29ZBu6PToHe4uDX+HJlXPaCqv28iBzMLgTejw2B7d8aK2qBzuasWYPAGDX4RMxTgnZAYN7DCmlsP3g8VgnA0D0qh5qukh3BCBysmxwj1a99Mw1e5A1cQa2Fh4z/difLtuJc59dgJz8ItOPbRez1uzBTC1HaycScI6k2D4A4mHa2tzth5A1cQb+OT02jbpHS8uRt9/a3S51g7uItBaR+SKyQUTWichdPvYREXlZRPJE5DcR6R2Z5EZ/2jBnUXnd7iOmH3vlDkeD1dZCc3LvudsPodRm9ay3frwCEz5eEetkREW8TIn3+oItsU4Clmw+AAB475f8mJz/mreX4fwXFsXk3GYxknOvAHCfUuo0AAMB3CYi3Tz2GQWgk/ZvPIDXTU0l6So4VIJLX/8FD35t7jTBP6zbi/8tyDP1mNFUUVmFqjiZgpXiR05+EY4E6C66eufhKKYmMnSDu1Jqj1JqhfbzUQAbALTy2O0SAB8oh18BNBSRFqantoYJpurpaKmjC916k0sY4z/MxTOzN5l6zGjq+NAsjP8wJ6rnzN1+CD0em4NDx91nE1y/+whKy/V7JRm57hWVVXhyxnoUHjVvxsWqKoXdJjXmRis4fplbgG9X7ULe/qMBqzdHvLAQT83cAAAoLa/EZZOX4qb3zLkvcvKL4qbtzFVSMDuLSBaAMwFke7zVCsBOl9cF2ja3ylIRGQ9Hzh5t2rQJLqVBeGvxVrRtUjdix480s7tATvg4FzPX7EX+pDFhHeehr9fg4+wdYR8n2n7cENnG4l/yDqC0ohLDujYDALw2Pw/HTlZgxY5DGH5as+r9Xpr3e8DjSBAXftHmQry5eBt2FJVgyri+oSXcw+RFW0x7kP+y5SB6tm5oyrECuc9jmhB/9+bm/cewef8xPDD6NFRoJbl1u4tNScNlk5cCABrXTTHleGYx3KAqIvUAfAngbqWUZ/bQ113plf9QSr2hlOqrlOqbkZERXEqD8O8ZG/C3D6KbWwuHswdF/oHjuOiVJdh2wNxcwMw1e005zsfZOwDAdvX64br6rWzcEGQucOhz87GmwD245GqDhowU2Kq0AkBFpXlVTku3HAzr961c/fXtql24OcolvEgzFNxFJBmOwP6xUuorH7sUAGjt8joTQOyb3OOca0bt65UFGPrcAqzZVYzX5ud5vR9Pbv+kZjRwRtK+IyfxwtzIVHf5C7EFh0rwRQTnTPooe3vEjh1pd01dhTnr9sU6GaYy0ltGALwNYINS6gU/u00HcK3Wa2YggGKllK36r0U6T/JzXni5Jk9KKRw/6T6U/WSFOTnuRVpPhnjV54m5+McXv2FPcc0eDOSZN7hyyq+4//PVhu+DYKqJAGD34fBG1prh2MmKiHRbtiIjOffBAMYBGCYiq7R/o0XkFhG5RdtnJoCtAPIAvAlgQmSSG/wNF4gz+A19dj6GPb/A9/lMO1tgZvfbf3vJNnR/bI5bgHtnSb65J4lTB4+XYVrOzpitoBPughKRGsNRqC13V3yiPGBPESv7y5u/Ytjz8bVm68mKSgx9dj4W/h7dea50G1SVUkugE+OU426+zaxERcOi3wtx7TvL8OnfBkZs5RijlAK+XGFucXn2Wkc9e8GhU8HdtLpy61atRpRSyi3zEa/Vav2fnAcAeP7ynri0T2aMU+NbqCN5VxeY00jqKWviDKx8ZAQahdBoWnDoBPIPluDx6evw9YTBOFZWgVYNa0cgle4sO0I1XNnbHNUguduDGx16oqwSa3eZcwMFHqUYHGeaik/YM0dmBVOX7wz4vpnX2wxTFoU3WGnMy4uruxcC1ptaIdjUrjHhez/8hQUYPCk6JcoaG9xDde9nq/CHV5aYuiL6xK/CH3j08k+bAQB7wpxR0CxlFVVxu9jDibJK3aqT1xdswR9f+zmo4+4ocpQAF2ttElv2h9bryV/alucXYeqyHUEfb/baPej88CzTr8e63UcwZdFWU48ZDbF8xB44Zl7c0GOb4L778AnsP+o7sBUePYkKk6a0XaFNGWBkMIpdGcmh9Xz8B5zx+JwopCY4uw6fwGmPzsaHvwbu2fH07I1YFeJAHGc/6pwgS4V6Lp+81CsjYCT3+fTsTQEDu1IKX+QWeHWrjK9yxinllVW45NUl+DkvcMN+VZXCphq8LJ9lg7tn5uasST9V1yV66vfkj/j3jA0+34uVG99b7ncmRisv1uF0orwy7h6AWwqPYbs2hsDqk5EZrc+fu36f7riJGWv24P7PV2PpVnN7bEXK3uJSrC4oxj++/C3gfq8v3IKRLy7yGk/gVOJnMZNwG8TjheWCu/OenvjVbzh4zPjQ6x/WnRrIc9vHK/Da/ODqG82+4PM27scBP+mftdacQUdmqrTwABXAMUBo+PML8f7SfFOP6znFgKtQ69jN/EsbGcznXHnJbpwlr90BusRuO3Dcq8vwbD/fv1CvS6yeFZYL7k7rdh/BU7M2em3P9bM02G6XuugZQeTazOx6aWVVLndoODfr3uJS/LIl+v3knXN/rN5pbm+KYyfDXxYvVow+fMorq/DB0nz9qk0LPv/Pe26BVxfFiF3TKIcSywZ3f37aaJ9RZqHeC5/l+O61EQ/fvVEvLcLVb3pOTeSw/0gp5vupqiqvrLJc/aln24RV8wnv/5KPR79dp9tOQfHFdsE9Uu6augrTlhvrqXCyohL7j0Sv10rR8TK3dTcn6tRFminY6qpDWhVA0fGy6n73Sik8/M0a9P/PPFz/3nKfvzdp1kaMfHER8k2edycehVuMN/r7eg8b5/vO7rXOmUfNsqag2O900tGsyrBJFbuXGhfcnb1dQvGPL411Wbz9k5Xo/x/fjbuR0PuJuSj3M4GUc7HkSNm0Tz833etfP+Dp2e5VaL2fmIu/vrMMAFClgI9+DfzgdF63gwHquM3wSfYOXPXGryH9rmddrVmLsAChDUCLVklheX4RdhYFPxDwoleXhDULpesAvWB4/l1enrfZ8O8eOl5WvZBIvLN8cP8ytwAfLs03vH+kBvm49p2eu96cqqFwMhTOG9isAVeuXNNlZFbCwyXlPlf3yd5mvKtgtHJXD369JmCvkZ/zDvgtwW07cNzt7715v3lznITaLTMaLp+8FOc8Mz/gPuPezo5qL5RFQQz1N5JBcbru3WW45u1sS8yMavngft/nq/HIt7FZZ9FpZ1EJTnt0Nj4JMMCkuKQct32yIiojSKsUsDw/9BJKPAs1N/rtKsckpaGOonR2J/zLW9luJTjP9Hj2vAjXzDV7AgYSpRSUUsjeejCqwbPgkHdOPdDZF4eR291ZVIKvXKbnMDJ+4LsIrAOrlMLv+xwP7Cqdv/XXKwuQV/1wd+xbcjK6DwTLBXfXL1Ow05deMWUpfgwiV210atut2hffXxcqAHhj8RbM+G1PUKUMOsXIDI/lAZbU8+wRIRCUV1bhv3N/R0mZfkB2ViGFI9gHU/bWg5jw8Qq3If6+jvXQN2tx5Ru/4vvfItN33zWOlVVU4fTH5uDspwPn1H1575f8kFYK+9P/fsa9n51alGNvsX4X6M9zC0wrQRvl+nC9Z9pqnP+C+wRme6PYDgcEuRKT1S3bVoRlQVQFfP/bHlzcs6XP96Ixj8aJKBf9KiqrkJSo/7yPxSCPfUccX+jKKuU1QRfg6Op47rML8OczPVeA9G/a8p14ad5mnKyowsRRXQPua9YIZyOc99YRrQFz1+ETfqulNu09ivmbHA+unT5y04GUVVbpLg/37+/X460l2wA4em8dPH4SR0MsnTz+3XoA/ldL8ieaQ/YB/yWQ0vIq3Rw7ED/dpy2Xc9cT6bgTzcmflHJfHef3fUeRNXEGsoMYSRjMnyMSjZVZE2fovv97EHWel09eiqe1RjjXHkLnPrsAAPCVTgOy6/U7qQ3JN2ue+2D9vOVgUJkNX3aHMZfQ1kLHA9FfLxgRVAd2p1hNoxwPbvkot/qeCSReRrjaLrjbzVVvnuq54ZxL4wETJhqLJ8GOyJ26fAe+yC3AwKfmIXf7IVMncYumsooqXDFlaayTYVoJ8Q2Dk4hd81Y2fvExL0xZRRWmLNyC8iiWksLlOqtsvOTYnSwX3MPJqYRircciup71Zkop0xba9efCFxe5vd4aal9vlxzFnuITGPb8gpBWu4+HfMlyLce7bncxpulMtRtJnl9oU/422kGcVUE/bthvqBrQM8NoxhS1kbAk7wDu+HSl1/Z3f96Gp2ZtxPu/5Ef0/IFK35NmbXRpCNUXTkNxpFmuzn3hpsiuZO9Jr6/y9NW7dfvqbi08hnUhNCQ5bYzAyMzPlhdga+Fxw1PIrtzh3hVPKYWvVuxC2yZ1qrfpLTDsbwKryqrQc2qPTV8XlW6SIoKb3vceYBWJvNrU5TvRrWWDkOu2nSb5mJ4jXvi6ZM6eRsdPVgZ8lG07ELll9CYvNDbnVKB77rkffje8byRZLrjrfZ0qojTB1cnyKsxZt9fQUz4Sy36VllciNTlRd78ik+rRPasPFm8+gPs+X+22TW+B4fOeW+Bze7CTuAGnGh2D/eI4f69KKa/BK/562jj9uEE/Y+E5SC53+yH0adsIQOAco+vEbMu2FaHXv+bqnssfs+p8PddEnbPevAntwrkvP8uJ3CLfZvg0hDn3I8Fy1TJ6jNb7+fLd6t1+Z2r09OTMDbj5w1zkxKg/+Yj/Gntg+BvFF6ivcFlFFSbN2uh3AiWlzB+KHozDJeUoOh7aeAFn3MveVuQ15sCMuVM8S3GrtcFHG/YcQVmAumTXXhi+uswVhzBzo9F8jr8Rmp4BeO2u0EufgVRUVuHeaauQF6GFrWNdE/5rmI3mobJgzj1y7vh0JXq1bmhoX+dw62AXGg51cNH8Tfuru5I5zh98XXmZy2hSZ12hr+//F7kFmBzjhq1py3cgs1Edv+//uMH8Psy/Rmg+851FJRj10uKwjuHZSyNQKcDZhXLDnsgEY7Mt/L1Qt5dTvHnw6zV4aeyZfueKd/XIN2ujkCJvlgvukW6QNjJYBgi9Hnzh74UhLXf2pAmLjeh1OXRmHp1BPVBw/zg7sjMEGp3Hx0yRmkdfb1TyvA37MKRzhuHjKaUwP8ptT8FasGk/bh7S3tC+N77v3lbz8k+b4379gG9X7cbt53X0u+BOPLBdtUy0eTaUuo6K81f32fnhWRFNUzT8ssUaq/Z48hVoBaJbAtvlp1eRkczGPp2RiZ7BTc/01bvxwdL4nn43e1tRyAujRDqw5+sM3DLq4leDW2M32hjcTWZk5RszvbXYeBvDTxv3404fXdD8icY8OL6c8LP8mRn8DUIJdXoBvQXJD5eUBR28fXHtCqn3sIgXu0KctdFsrtOCHC4pC7h6VjA8xweY1XnBLJYL7pFuHHEOcw+FWTmCYAS7Nuz01bvx37mnumr5GpXqLHGEMg+IGV788Xf9nUy081CJV1dPo+7/bHXA91/+yfd85eGIkwGQlrHSZUZNM+YIsgrrBfdYN30HEEojp1HBDKzQ85JL7wjXXPLAp+ZhustsemZOWRuMYOdICVc4E0yFPKDMQzC3tZn3Qk3gmqNebaAB1C4sF9zJXF979FJ4JYiFCyJl5pr4WyA81ly7WH7uZzbUIzGqRvMnXjJinh0Dolnw0ZtbKZIs11vGSpSKnxucrE2vbh8ApizaGrPSFsUfy+XcozkroxkiOUOcv8Wkw3GopAx7LNJgR95+iuOueRRdzLlH0NGTFaavzOMqmKlyjTpwrAxTFoY+ypdC0/Eh63ePjVclHr2v4r0PvVksl3M3W6SGVAPAOU9Hdu7rp+J4Yiii/UdD73kWSQ9+ba8ps/2xXHA3uw579rrINd4dKa2I+mpKRPHC37xGsRav6TJbjQ/ukVZRWTOKgESe4mVFoprKcsGdiIj0MbhH2EiPVZSIaooVIY76JXPoBncReUdE9ouIz3krRWSoiBSLyCrt36PmJ9PlfBbrCklEFAtGukK+B+BVAB8E2GexUuoPpqSIiIjCpptzV0otAhCbpUSIiCgkZtW5DxKR1SIyS0S6+9tJRMaLSI6I5BQWFoZ2JtbKEBHpMiO4rwDQVinVE8ArAL7xt6NS6g2lVF+lVN+MDOMrz7hibCci0hd2cFdKHVFKHdN+ngkgWUTSw04ZERGFLOzgLiLNRRxDi0Skv3ZMa67BRkRkE7q9ZUTkUwBDAaSLSAGAxwAkA4BSajKAywDcKiIVAE4AGKs4NI2IKKZ0g7tS6iqd91+Fo6skERHFCY5QJSKyIQZ3IqIoq4rCnPKWC+5itWkhiYg8vBiFtYotF9yJiKzu+9W7I34OywV3dsQhItJnueBORGR10ciiWi64s86diEif9YJ7rBNARBSmouNlET+H5YI7EZHVFZ8oj/g5GNyJiGyIwZ2IyIYY3ImIbIjBnYjIhiwX3NkTkohIn+WCOxER6WNwJyKyIQZ3IiIbYnAnIrIhBnciIhticCcisiHLBXf2hCQi0me54E5ERPoY3ImIbIjBnYjIhhjciYhsiMGdiMiGLBfcuYYqEZE+ywV3IiLSx+BORGRDlgvurJQhItJnueBORET6GNyJiGyIwZ2IyIZ0g7uIvCMi+0VkrZ/3RUReFpE8EflNRHqbn0zXE0b06EREtmAk5/4egAsDvD8KQCft33gAr4efLCIiCoducFdKLQJQFGCXSwB8oBx+BdBQRFqYlUAiIgqeGXXurQDsdHldoG3zIiLjRSRHRHIKCwtNODUREfliRnD3VQuufO2olHpDKdVXKdU3IyPDhFMTEZEvZgT3AgCtXV5nAthtwnF9EraoEhHpMiO4TwdwrdZrZiCAYqXUHhOOS0REIUrS20FEPgUwFEC6iBQAeAxAMgAopSYDmAlgNIA8ACUAro9UYomIyBjd4K6UukrnfQXgNtNSpONIaXm0TkVEZFmWG6E6d/2+WCeBiCjuWS64ExGRPgZ3IiIbYnAnIrIhBnciIhticCcisiEGdyIiG2JwJyKyIQZ3IiIbYnAnIrIhBnciIhticCcisiEGdyIiG2JwJyKyIQZ3IiIbYnAnIrIhBnciIhticCcisiEGdyIiG2JwJyKyIQZ3IiIbYnAnIrIhBnciIhticCcisiEGdyIiG2JwJyKyIQZ3IiIbYnAnIrIhBnciIhticCcisiEGdyIiG2JwJyKyIQZ3IiIbMhTcReRCEdkkInkiMtHH+9eJSKGIrNL+3WR+UomIyKgkvR1EJBHAawBGACgAsFxEpiul1nvsOk0pdXsE0khEREEyknPvDyBPKbVVKVUGYCqASyKbLCIiCoeR4N4KwE6X1wXaNk+XishvIvKFiLT2dSARGS8iOSKSU1hYGEJyiYjICCPBXXxsUx6vvwOQpZQ6A8CPAN73dSCl1BtKqb5Kqb4ZGRnBpZSIiAwzEtwLALjmxDMB7HbdQSl1UCl1Unv5JoA+5iSPiIhCYSS4LwfQSUTaiUgKgLEAprvuICItXF5eDGCDeUkkIqJg6faWUUpViMjtAOYASATwjlJqnYj8C0COUmo6gDtF5GIAFQCKAFwXwTQTEZEO3eAOAEqpmQBmemx71OXnBwA8YG7SfEuvVwsHjp3U35GIKE5l1K8V8XNYboRqv6xGsU4CEVFY6qYkRvwclgvur13dO9ZJICKKe5YL7gkJvnpmEhGRK8sFdyIiq0uMQiaVwZ2IyIYY3InINro2rx/rJMQNBnciiroXr+ylu0/juilBH7duLUO9u2sESwb3J//UI9ZJILKtWkmRDwuiU+U8ZVwf5Dx0ftDHVcpz2quay5LBPb1eeAMAnvhj+A+Hb24brLtPer3gcx7+eH4ZLuzeHJ/8bQDm3Xdu0MeaMLQDXriiZ1C/0zMzzef2O4d3MvT7H904IKjzhePnicP8vlePObuoCKbB0Fef75Hdm4fUM+7czk2D/p1g1U+1xj1kyeAeysO5WYNTD4TayeEPIOjVuqHuPjed0z7s8wDAE5d0d3v98U0DMHlcH5zVIR0dMuoFfbx6qUn4c+9MmNFgf++Izrr7XNE3E2d3SkejOsnhn9CAVg1r+33PWerrn9U4KmmJlOcv74kxZ7TAdWdlRfQ8j/yhW/XP9U18MJ6R6fj+PDCqK/5+YVfTjuv6PXfq6fJdbZde1+fvBfMwyggzcwkAF3RvHvYx9FgyuIdiUPsmPrd/eesgw8fIbOQ/aOjpaeBhAAC3Du3gtW3coCx8PcFRUqiTkojBHdNDTgcApNV2BNnOzUJrfPIstjsD5tAuGfjhniFe+zu/OOKjLH7fiM648ex2IaXDqDGnO+a1+/cfe6B5g9Tq7Y+6BC6jPrvZ+P3iz2V9MsM+BuAY0PfPi089+J+7vCd+++cFYR/3Upf0pUSoiqZdel2sfGQExg8xJwPklJTonV7XDEhTP8P+3/prX7/HvO28DljxyAjDaXhw9KmHled35ZWrzsSKR0bg/gu6GD5eqGpMcPenT1tjObhOTethyT/ci/urH7sAOQ/7rxc8WwvCDVKT8K2BahwAaNO4js/tvVo3xIL7h3qlwZd3r+vnc/vwrk3xzKVnYGy/NobSclV/n2uueGmRdipgdm5WvzqYOqVoX7iJo7xzaHcM74SHRp9m6DyhalDbkeP0fLY4C4BX9W+DDf+60NCxzCj1maG5y998ZPdmAICzOjRBg9RTpSN/1QdPX3q62+vrB2e5vb7QJVc5oN2p70eWn1yvr8Dn/L1AVaiN6qb4fOCH4vzTHH8D1wy40Wu17MHhOK+L/+qcW87tgMQQ07nqUfeH7UU9W6Jx3RT2czdTsssT3bPR5Yd7huAXj3paIzmWtNrJXjdvp6aOapIHR3dF60a+AzUAXN4nE1lN/L/v5Fr9k5Ve11APgi7N66NvW+85eBISBFf0a119Y13Yw71o6NkWUSspuEAmXj843D/SkUu5om9r5E8ag+wHh/s9xmMXOXLT08YPRP6kMX73u8tHXf/VA9pgyT/OM91P/6oAAA9ASURBVJze5ERHQhvUTkJtP3N9DO7oXuJTXuvUhO9/f+mt2yZx7aC2aJ9eFy+N7YXPbh7kVnqbMq4v8ieNQUuP6qivJwz2+Xe6sl8bt+11U3w/BFqkpaJzs/rV7SrDT2uK7+8422s/X/fkm9f2xey7zzHUJjT69BZo07hOwOo0z2N7undEZzSum4IhnU8tAtRXm4dKL4w21Upzw7t6B/jEBEH91GSk1Un2eggCwMNjvDMmfVy+e6731fz7h+qkxFy2C+6edcCTr+mN8UPa4+Exp4rgrrkewJHbbNmwNp657Aw00HI7U65xX2/EaN321PED8fzlPXHDYO+qBmfgAhwPhn4+6n09b8QmIXQHM+rOYZ2w2iVnMW5gW7eG4ntc/paBclh6bSD1U93r2ps1cP/7izga1Z74Yw9cP7gd1j0+EgP8VKOdOuapgFQ7ORF5T47Cf/50OjK1B6qvulfAfQmxK/u1xi3ndsCdw/w3Cjeue+o4L191JlIjkHMffXqL6kDkzy3ndsBP9w/FJb1aoX+7wKVN5z2TVju5+hpm1K+FMae3wJRxjvv6nhGd3aqoPC24fyhm3+VdxdajVRp+vNc7YL/+F/c5n+rWSkLX5g2QVju5Ojc9/XbHw8bzgZNRvxYW/f08/OlMX6t3eju3cwY6ZNTFe9efKqF2a9kAKx4ZgfR6tfDJ3wb4LQ0+PKYbFv/9PJ/Vh9cMbOv2ekjnDKx4+FSp5NZzHVWmZ2Sm4cMb++Od6/ripnPau2Wkfp44zKs2YGB7x2t/9f2RYo1m3xA5c34X9nBUE9RPTcLR0gr0aJmG809ripKySrf9r+jbGjPX7MGCTd7ru17Us6Xu+a4e0AZN6tWqrrOsrcXlW4d2BABcP7gdHv9uvc/fvahnS3y32rHA1Tmd0rF48wEDn9A31y/tG+P6YPyHuT73S0gQpHk0cvZq3RCrH70Ax8sqquvmAeCFK3piysKtmJaz0/Mw1ZwPgAdHn4byiir8sH6fofSKCNa5VIv46qvcqmFt7Dp8Ajee3Q73X9AFK3Yccnvfs641+8HzkTVxRoCTOkomvqqKXDkD5WMXdcPF2j1wcc+WGH5aU2Q2qo3Wjeqg/3/mAQA6Nq2HF6/shbTayTjnmfnVx7h2UFt8sHR79WvX6/NnLaC5PjvvGNYR367ajbZN6lTfB6HUCjh/56sJZ6F1ozpe08w6SyH1fFTf+KuCARyf85E/dMOctXsxUiv9jfKoivOle8u06obUcKQkJWDefUP9vn9Wh3Sc1SEdizZ7f49P99PrC4BXziq9Xorb96Npg1RMv30wOjer7/aQf+XqM/HRr9tx/wVdfGaC3r2uf0ymKbd8zr2FRy7c+Xqcx1MYOJVjTBDBW3/th0/+NtBrH2euqFWj2sifNAajtJvXyJfLs4EuJSkB+ZPG+Gwk9cx91U5OqD7Ph0F2Gzzbo4HVtQtZo7op1e93bOq/9OHa0JRWJ7m6iH/P+Y6cX9smdfG41mvHNegDwID2TXB6qzT8/UJH9UurhrXxho+iczicpa1RPZqjttao/OO9jtyXfi+c4CPjncMcD+TGdVOQP2kMrncpib181Zm4pFcr9GnbuLpI79SjVRpau7SbPDi6K67s52i7OL1VGr6ecFZ1Vdvwrk3xgo/BPHcO74RFfz/PLScpIXwGp95tGvmcP9xZ4rq4Z0s8cUl33RKBqxvPbofPbhnk1hjeKcD9ZYTzOza0i//1lYOZ8rutVu3peb8aTo+Pv/kZmQ29Sm8t0mrj/0Z29Vu6rZ2S6HZPRIulc+4XdGuG+y7ogpEvLgIAnNclA5f1yUR6vVo+b5CPbhyARZsLvXKrrm4Z0gEXndHS62I4vwidmtbzqqIZ0a0Z5q7fZ6jInj9pDA6XlKFhnRTMDZCz/fcfe+Dhb9biBgM9ST64oT86PjQTVS51Ds560JTEBPRu2whL8g74HZzizNn5ctf5nXDX+Y5idGJCIp64pDuGdM7Ate8sQ4J2M9erlYTvfNTFfn/H2ablWK47Kwu52w+hvcvfvmPT+njyTz0wpFPgxdb/b2QXVFUpXNo7E6t2HtY914w7z0anpvWRlJgQdm+O8UM6QCmFCUM74OoBbZDZqA7mb9oPwHe7Tkpiglv7kFMwYyaMtgrUS03C/qMnkZQoGDcoy3BJy59Zd52Djg/N8tqekpSA0vIqra1L/yHVu82pAD51/ECMfeNXAI7v96tBTPn98JhuGN61meGeah21e6t/VmMsyy8yfB5fBrRrjFoxbny3ZHAf1rUpxvZrjXtGdEazBqlY9egIpNVOrn5ynuejYQRw5P6u6Bu4B0hCgrgFds+H8Vwf9Y2vXd0bJWUVhtPfsI73F/WWcztgxY7DGNHNUVK4ZmBbrzrAQGnu1bohVuw4XN0o9exlPXHOb7txRmYa5m10BBN/deOuXyY94wZlAQAW/p9+w2WPVgGKwEG6qGdLn1Vjfxmg/zdqVCcZT192huFzdW/pSLfRAVpOQzv7fsiIiFtf7iGdMjBhaAe3XG9SgiOg3+JSynPeeuef1sxnFz89emH0gxv6Y/bavWha31H6uP28jsjJP4SeHlUn55/WFC/P24xhfr5XTv7S+M1tgzF33T7dz+BMr+t9OtCl7SUrvW5Q0wukJidWx4LOzerh933HAu7funEdbH5yFL5ZuSvs4D7NhC6z4bJkcE9JSsCkS099WX0Fy2hKSUpASlLwaXAtKrfPqOezocqoVo3qYMWOw3hIa71Pq5NcXTVl5RnwXxrbK6RBa4Bj8NRnOQWG9//x3iFB9xACgO9uPxvHTla4VRl8fssgvw2WiQniNXAnMUEC9g6KhMxGddwG2g1o3wQbnvDuEnpGZsOw0ta1eQN0bd5Ad7+/DWmPPcWluPEc99JqVpM6yD9Yggla25WrXyYOw5HSct1jf37zWdhdfMJtW/1aSTh60j1T5qvUZFWWDO52cff5nfG/BVtMOdbjF3dHi7RUXNCtmSnHiwYjXd8u6WWsB4UvT/35DDx6UXe3ulBnwB3go365Y9PQBnX5aqTz1RMqWL21Xhg3nJ0V1O/Fcn6V+qlJuGOYdxA29rvJePZy72kxvrltMLYeOO6z3aBlw9poCf37KK1Osld17Jx7hmBLYeDcvJUxuMdQSlICNj5xISqqwv8yNq6bggcjPBjITDPvPMerMdxsiQniNZdMVnpdLLh/aEwauIKVXq9WWDlmswYIBWPNP0eafsyGdVLQu435pfOWDWt7jQ0AgE7ayG1fGQArYXCPsUj0m7aCbi31i+mREqibH1Gv1g3x6wPD/Y6VsAr7VDBFiLMO1spVcVcPaIOuzevjqv7Gph0ga3P2KEpKtHJrS2w1T0uNScnHTMy563jsom5onpZaPXeFFTVrkIrZd3uPyCN7euvavlhVcNhtnhmqeRjcdTSsk4J/mDglKVGkNaqbEnAiLKoZLFzZQERE/jC4ExHZEIM7EZENMbgTEdkQgzsRkQ0xuBMR2RCDOxGRDTG4ExHZkMRqBjkRKQSwXXdH39IBhL4OnTXxM9cM/Mw1Qzifua1SKvAKNYhhcA+HiOQopcxdxy3O8TPXDPzMNUM0PjOrZYiIbIjBnYjIhqwa3N+IdQJigJ+5ZuBnrhki/pktWedORESBWTXnTkREATC4ExHZkOWCu4hcKCKbRCRPRCbGOj3BEJHWIjJfRDaIyDoRuUvb3lhE5orIZu3/Rtp2EZGXtc/6m4j0djnWX7X9N4vIX1229xGRNdrvvCxxslaYiCSKyEoR+V573U5EsrX0TxORFG17Le11nvZ+lssxHtC2bxKRkS7b4+6eEJGGIvKFiGzUrvcgu19nEblHu6/XisinIpJqt+ssIu+IyH4RWeuyLeLX1d85AlJKWeYfgEQAWwC0B5ACYDWAbrFOVxDpbwGgt/ZzfQC/A+gG4BkAE7XtEwE8rf08GsAsAAJgIIBsbXtjAFu1/xtpPzfS3lsGYJD2O7MAjIr159bSdS+ATwB8r73+DMBY7efJAG7Vfp4AYLL281gA07Sfu2nXuxaAdtp9kBiv9wSA9wHcpP2cAqChna8zgFYAtgGo7XJ9r7PbdQYwBEBvAGtdtkX8uvo7R8C0xvpLEOQfdhCAOS6vHwDwQKzTFcbn+RbACACbALTQtrUAsEn7eQqAq1z236S9fxWAKS7bp2jbWgDY6LLdbb8Yfs5MAPMADAPwvXbjHgCQ5HldAcwBMEj7OUnbTzyvtXO/eLwnADTQAp14bLftdYYjuO/UAlaSdp1H2vE6A8iCe3CP+HX1d45A/6xWLeO8gZwKtG2WoxVDzwSQDaCZUmoPAGj/OxfA9Pd5A20v8LE91l4E8HcAVdrrJgAOK6UqtNeu6az+bNr7xdr+wf4tYqk9gEIA72pVUW+JSF3Y+DorpXYBeA7ADgB74LhuubD3dXaKxnX1dw6/rBbcfdUrWq4vp4jUA/AlgLuVUkcC7epjmwphe8yIyB8A7FdK5bpu9rGr0nnPMp8ZjpxobwCvK6XOBHAcjqK0P5b/zFod8CVwVKW0BFAXwCgfu9rpOuuJ6We0WnAvANDa5XUmgN0xSktIRCQZjsD+sVLqK23zPhFpob3fAsB+bbu/zxtoe6aP7bE0GMDFIpIPYCocVTMvAmgoIknaPq7prP5s2vtpAIoQ/N8ilgoAFCilsrXXX8AR7O18nc8HsE0pVaiUKgfwFYCzYO/r7BSN6+rvHH5ZLbgvB9BJa4FPgaMhZnqM02SY1vL9NoANSqkXXN6aDsDZYv5XOOrinduv1VrdBwIo1opkcwBcICKNtBzTBXDUR+4BcFREBmrnutblWDGhlHpAKZWplMqC43r9pJT6C4D5AC7TdvP8zM6/xWXa/krbPlbrZdEOQCc4Gp/i7p5QSu0FsFNEumibhgNYDxtfZziqYwaKSB0tTc7PbNvr7CIa19XfOfyLZSNMiI0Zo+HoZbIFwEOxTk+QaT8bjmLWbwBWaf9Gw1HXOA/AZu3/xtr+AuA17bOuAdDX5Vg3AMjT/l3vsr0vgLXa77wKj0a9GH/+oTjVW6Y9HF/aPACfA6ilbU/VXudp77d3+f2HtM+1CS69Q+LxngDQC0COdq2/gaNXhK2vM4DHAWzU0vUhHD1ebHWdAXwKR5tCORw57RujcV39nSPQP04/QERkQ1arliEiIgMY3ImIbIjBnYjIhhjciYhsiMGdiMiGGNyJiGyIwZ2IyIb+H7nfIb6HftMSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pd.DataFrame(data).plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fbac1e6be50>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU9f3v8dfnzJJ9I4Q1wQCybxGjKKgouKBW0dZWrFXbqz9/1tvdeutt77XWtr8u9v6s1lprq8UFXIqiaFEEFRc2CauECAlLFhIgCdn3mfnePxJoxJAMMJMzy+f5eOTBzJxvznmfTPLhzPme7/eIMQallFLhz7I7gFJKqcDQgq6UUhFCC7pSSkUILehKKRUhtKArpVSEcNq14YEDB5rs7Gy7Nq+UUmFp06ZNVcaYjJ6W2VbQs7OzycvLs2vzSikVlkSk+ETL9JSLUkpFCC3oSikVIbSgK6VUhLDtHLpSKnx1dHRQVlZGa2ur3VEiVmxsLJmZmbhcLr+/Rwu6UuqklZWVkZSURHZ2NiJid5yIY4yhurqasrIyRo4c6ff36SkXpdRJa21tJT09XYt5kIgI6enpJ/0JSAu6UuqUaDEPrlP5+WpBV0qpCKEFXSkVdmpra3n88ccBWL16NV/60pf6dfurV69m7dq1x54/8cQTPPvss6e0rsTExEDF0k5R1bvFG0qOPf76jBE2JlHq344W9Lvvvjto2/B4PDidPZfI1atXk5iYyMyZMwG46667gpbjZOgRulIq7Nx3333s2bOHnJwc7r33XhobG7nhhhsYP348N998M0fvxLZp0yZmz57N2WefzRVXXEFFRQUAW7du5bzzzmPq1Klcf/311NTUAHDxxRfz05/+lNmzZ/PII49QWVnJV77yFc455xzOOecc1qxZw/79+3niiSd4+OGHycnJ4aOPPuKBBx7gD3/4AwBFRUVceumlTJs2jenTp7Nnzx4aGxuZO3cu06dPZ8qUKbz++utB+bnoEbpS6vRs+gHUbA3sOtNy4Ow/nnDxb3/7W3bs2MHWrVtZvXo18+fPJz8/n2HDhjFr1izWrFnDjBkz+O53v8vrr79ORkYGL730Ej/72c94+umnufXWW/nTn/7E7Nmzuf/++/nFL37BH//Yub3a2lo++OADAL7+9a/zwx/+kAsuuICSkhKuuOIKCgoKuOuuu0hMTOTHP/4xAO++++6xbDfffDP33Xcf119/Pa2trfh8PtxuN0uXLiU5OZmqqirOO+88rr322oB3LGtBV0qFvXPPPZfMzEwAcnJy2L9/P6mpqezYsYPLLrsMAK/Xy9ChQ6mrq6O2tpbZs2cDcNttt/HVr3712LpuvPHGY49XrVrFzp07jz2vr6+noaHhhDkaGho4cOAA119/PdA5OAg6B2L99Kc/5cMPP8SyLA4cOMChQ4cYMmRIgH4CnbSgK6VOTy9H0v0lJibm2GOHw4HH48EYw6RJk1i3bt3n2tbV1fW6roSEhGOPfT4f69atIy4uzq8cR0/1HG/RokVUVlayadMmXC4X2dnZQRllq+fQVa8SO4oZV/s0OdW/gV2PQeNeuyMpRVJSUq9HygDjxo2jsrLyWEHv6OggPz+flJQU0tLS+OijjwB47rnnjh2tH+/yyy/nscceO/Z869atvW4/OTmZzMxMXnvtNQDa2tpobm6mrq6OQYMG4XK5eP/99ykuPuEMuKelz4IuIk+LyGER2dFLm4tFZKuI5IvIB4GNqGzRXgtrb+Xakos4u/oXjK15EjZ9F7PsTHxrb4X2GrsTqiiWnp7OrFmzmDx5Mvfee2+PbdxuN0uWLOEnP/kJ06ZNIycn59ilhs888wz33nsvU6dOZevWrdx///09ruPRRx8lLy+PqVOnMnHiRJ544gkArrnmGpYuXXqsU7S75557jkcffZSpU6cyc+ZMDh48yM0330xeXh65ubksWrSI8ePHB/Cn8W9yoo8IxxqIXAQ0As8aYyb3sDwVWAvMM8aUiMggY8zhvjacm5tr9AYXIaq5DN6dg2nYy9M1N/DUwSuQhBGk+g5wTfwy7sh4DSshE2vOSkgeY3daZYOCggImTJhgd4yI19PPWUQ2GWNye2rf5zl0Y8yHIpLdS5OvA68aY0q62vdZzFVoWryhBLe3lisOzCfOW8Wt+37DbnMW15wzlDMHJeL1jWZL6VRu2DWTZ0b9kqSVF2JdsRYSR9kdXSlFYM6hjwXSRGS1iGwSkVtP1FBE7hSRPBHJq6ysDMCmVUAZH+cf/gFxHWUsKLyftgGzuGv2KMYMTkJEcDoszskewITJl/HVov+iubUZ3+ovQXvvnUxKqf4RiILuBM4GrgauAP6viIztqaEx5kljTK4xJjcjo8d7nCobjal/nuHN7/PLA3cQP3w2z98+g3j3Fz/ETRmewrfnX8Mde/83pr4QszF4o/WUUv4LREEvA942xjQZY6qAD4FpAViv6k/N5Uyt/h1rm3JY0fEVHr95OnFuxwmbt7T7iM+6lIcPLkCKF8P+F/oxrFKqJ4Eo6K8DF4qIU0TigRlAQQDWq/qRZ+vPEF8bvzz0PW45P5vEmL6HKMwZP4gPnLezpXkc7Z98n3+u/ZTFG0o+N/+LUqr/+HPZ4gvAOmCciJSJyO0icpeI3AVgjCkA3ga2A58AfzfGnPASRxWC6gqw9j/Ls1VXcf60GaTFu/36NkuEr5ydzZ8a7sHRUc3U6oeCHFQp1Rt/rnK5yY82DwH61xym2rb8HI8vhlX8D64d1DmVp79H2W6nxbDs83m+9Eq+IYspTP0m9e4zgxlXhaBAfyrzZ2bPxMREGhsbjz1fuHAheXl5PPbYYzzwwAP84he/oLCwkDPP7Px9fPjhh/nRj37Exo0byc3NpbGxkXvuuYdVq1YRGxtLeno6Dz30EDNmzAjovvQnHSka7Rr34ip/heerr+K8Sac22GHC0CTe8N5Bq9fN5MrfBTigUqdmypQpvPjii8eeL1myhIkTJx57fscddzBgwAAKCwvJz89n4cKFVFVV2RE1YLSgR7m6Lb/H47Oozfo2g5JiT2kdIsLMyRP5W9X1ZLe+Q2qbdqEo+1133XXHpqndu3cvKSkpHL26bs+ePWzYsIFf/epXWFZnGRw1ahRXX321bXkDQQt6NOuoJ6bkWd6sm01qevZprWpIciybYr9BozeOcUceD0w+pXrR0tJCTk7Osa/jh+8nJyeTlZXFjh07eOGFFz43i2J+fj45OTk4HCe+kiscaUGPYu1FzxMrLayzvkpSrOu01zcpeySLjlzJyOY3dRIvFXRxcXFs3br12NeDDz74hTYLFizgxRdf5LXXXjs2pW0k04IerYyheeefyW8ZRUrWBQFZZVZaHG+0fg2PcUDB/wvIOpU6Hddccw3PPfccI0aMIDk5+djrkyZNYtu2bfh8PhvTBZ4W9Ci0eEMJyz9aQWrbTt5oupqRAwNzk1oRYdSIsbxy5BJ8e56GlkMBWa9SpyouLo7f/e53/OxnP/vc66NHjyY3N5ef//znx+YwLywsDNqt4fqL3uAiSg2sXorHWFQPmM+IAN4GKycrlad238CN6Suh8C8w9YGArVuFrlC+gfiCBQt6fP3vf/8799xzD2eeeSbx8fHHLlsMZ31OnxssOn2ufRavL+bSovMpaBrG7gkvk+DHqNCT8crmMn7g+D6zBlZgXbcfLD1uiDQ6fW7/ONnpc/WUSxRKa93CIKuCTVwe8GIOcG72AJ6tmofVegDK/xXw9SuleqYFPQqlVy+lzeekdkBwrrnNTIujImEOVd4MTOETQdmGUuqLtKBHG+PjzJblfNiYS/bQYUHZhIiwYMYonq+6FCpW6CWMEcqu07XR4lR+vlrQo4w5/BEDrCo2+i4n1hW8QRXzc4bxesOVGASK/ha07Sh7xMbGUl1drUU9SIwxVFdXExt7cqO3tbcqylQXvECiz039gCuCup2kWBfnTZ7G+5Xncsmef2BN/aV2jkaQzMxMysrK0DuPBU9sbCyZmZkn9T36FxZlHAffYl3jVM6cOCTo2/r6uWfw2LNzmJv0X3BwFQybF/Rtqv7hcrkYOXKk3THUcfSUSzSpLyTNV8KnZibxQbi65XhTMlM4lDSHOm8S+/Ke0BtfKBVkWtCjSHXh0s5/U+b22za/du4Y3qi9gMzGFTh9jX1/g1LqlPlzx6KnReSwiPR6FyIROUdEvCJyQ+DiqUBqLX6TotZMBg7uvwEh1+YM4436ubhoJavxrX7brlLRyJ8j9IVAryc/RcQB/A5YEYBMKhg8TQxqWc8mz0zSE2P6bbOJMU5M+vkUtw3ljPpX+m27SkWjPgu6MeZD4Egfzb4LvAIcDkQoFXjtB1bikg5aBwb36paenDsynVdrLmFI23poLuv37SsVLU77HLqIDAeuB/ocEigid4pInojk6eVO/atq11IavXFkT7i8X7a3eEPJsa9hqXGsN/OwMJjil/tl+0pFo0B0iv4R+IkxxttXQ2PMk8aYXGNM7tFbQal+YAzx1StZ25TD/hqPLRGyRkzl0+bRNOx+wZbtKxUNAlHQc4EXRWQ/cAPwuIhcF4D1qkBp3EuqqWAXM3A57LmwacrwFN5pvIjkpjxoKrYlg1KR7rT/uo0xI40x2caYbGAJcLcx5rXTTqYCpnpvZ191bXJg7kx0KlwOi4qUawBY855O2KVUMPhz2eILwDpgnIiUicjtInKXiNwV/HgqEOqL36GyI5WUjCm25hiVPYVtzWMYUveGrTmUilR9Dhc0xtzk78qMMd88rTQq4BavL+aK+jVsbJ1GelL/Xa7Yk9R4N3m+S7nd+gvtdXtxp4yyNY9SkUZHika4xPa9pFtV7LPOQQJ4q7lTVZ/R2b2yb/NCe4MoFYG0oEe4uNqPAHvPn3c3cMg4drSMwV3+qt1RlIo4WtAjXFrjWg52DCA5Y6LdUQBwWEK+Yw4jJZ+66n12x1EqomhBj2TGMMq3kU2tOaQl2Hv+vLumjGsB+GzTIpuTKBVZtKBHMF9tAQOsIxQ7zrE7yufEDJhIqScLd4Ve7aJUIGlBj2AVRe8A0JAyy+YknyeWRXXqlUx2bqK44oDdcZSKGFrQI1jzgQ86rz8f2H/T5fora+pNuMTLZ5tesjuKUhFDC3oES2nK49O2SSTHu+2O8gXpZ1xEjS+duMN62kWpQNGCHqFM80EGSRmljhy7o/RMLA4lX8501wb2VOjMm0oFghb0CHVo33sA1CSEVofoUYs3lFDkvpRERws7NuvUP0oFghb0CFVb8iFtPheOAWfbHeWEGlMvpMkXj7vidbujKBURtKBHqLja9exsG0tacpLdUU7IJzEUWBeS61pD4cFau+MoFfa0oEcibytDfTupcJ2FFQLzt/SmJu0qMly1/HPFEhZvKLE7jlJhTQt6BGooX49bOiAjtK4/70ltyhw8xsHwpnftjqJU2NOCHoEO7X0fgIyRF9sbxA8djmSKmM6MmHVUN7bZHUepsKYFPQKZw2vY1zaMyaPG2B3FL1XJlzE+rpiD5Z/ZHUWpsObPHYueFpHDIrLjBMtvFpHtXV9rRWRa4GMqvxlDeutmdnkmsnRLeAyrr0m9HICB9atsTqJUePPnCH0hMK+X5fuA2caYqcAvgScDkEudIm9jKQMcNVQ4p9odxW8N7pEcNFmc5VjD4fpWu+MoFbb6LOjGmA+BI70sX2uMqel6uh7IDFA2dQoO7f8QgIb4s2xOcnIOxM/l/MTtvLtD50hX6lQF+hz67cBbJ1ooIneKSJ6I5FVW6nDvYGg4sJZ2nxNJC9Eh/ydQk3o5MVYHFbvetDuKUmErYAVdRC6hs6D/5ERtjDFPGmNyjTG5GRkZgdq06sZdl8dnraNCekBRT6rizqHZJDCs6V3qmjvsjqNUWApIQReRqcDfgfnGmOpArFOdAp+XIZ4dFPkmhPyAouP5xE2J+0IuSdrIuwUH7Y6jVFg67YIuIiOAV4FbjDG7Tz+SOlWe2p3ESQuHXOHTIdpddeplDHYdYefOD+yOolRY8ueyxReAdcA4ESkTkdtF5C4Ruauryf1AOvC4iGwVkbwg5lW9qNz/MQANCdNtTnJqDiZcgg8h+cgKWtq9dsdRKuw4+2pgjLmpj+V3AHcELJE6Zc0Va6n3xuMeON7uKKekzZFOY8LZXNS8gQ92VzJv8hC7IykVVnSkaIRYvKEEObKRHS1jGZgcZ3ecU5Ywaj458YWs3ZFvdxSlwo4W9Ajh8LUywipir5kYdh2i3TkyrwHAlC+nw+uzOY1S4UULeoRIbs3HKV4qXVPsjnJaFu9KpU4GMzNuPev36gVTSp0MLegRwt2wFYDmxPAaUPQFIhxMvJQLE7ewaofOj67UydCCHiESmz7liCeZhLSRdkc5bRUJc0l0tFC9dxU+n7E7jlJhQwt6hBjk2cHO1tEMTIq1O8ppOxQ3kw5iONu1ht+/rVPqKuUvLeiRwNtGlrWXYjM+rDtEj/JacVTEzmRO0kbyK+rsjqNU2NCCHgE8R7bjFC9Vrkl2RwmYg4mXckbMQZoq8zFGT7so5Q8t6BGgqmQ9AC1JkXNvkfKEOQDkutaw61CDzWmUCg9a0CNA08EN1HsTiE0Nj1vO+aPZOYxq53jmJG1kxY5DdsdRKixoQY8AMQ3byI+QDtHuKhLnkpuwk493FtodRamwoAU93Pk6GOTdTYlvXER0iHZXnjAXp/gY3PQBpUea7Y6jVMjTgh7mPDX5uKU9ojpEj6qOyaFF0piTvJEV+TpHulJ90YIe5ipL1gHQHEEdokcZcXAw4RIuTdnEOzsO2B1HqZCnBT3M7fnsQxq9cbhSw3PK3L4cSJhLslWPr2o9lQ1tdsdRKqT5c4OLp0XksIjsOMFyEZFHRaRIRLaLSHjeXSFMZXjy+axtFAOTwnfK3N5UxF2IESeXJG1k5U692kWp3vhzhL4QmNfL8iuBMV1fdwJ/Of1Yyi8+LyOkkP0R2CF6VIcjhUMxuVyeksc/1uyzO45SIa3Pgm6M+RA40kuT+cCzptN6IFVEhgYqoDoxT91nxFmtVLkm2h0lqMoT5jImZh+tdftoavPYHUepkBWIc+jDgdJuz8u6XvsCEblTRPJEJK+ysjIAm45uh4vXAtCUEOZT5vbhQPxcAGYnfsK6PTpHulInEoiC3tNn/R4n3zDGPGmMyTXG5GZkZARg09GtsfwTWn1unGkT7I4SVA2uUdQ7s7kseSMf7NYDAaVOJBAFvQzI6vY8EygPwHpVH1z1W/msdSTpSQl2RwkuEcoT5nB+4nbWFxbrZF1KnUAgCvoy4Nauq13OA+qMMRUBWK/qjfExyFsQ0R2i3R2In4tb2hnRvoF9VU12x1EqJDn7aiAiLwAXAwNFpAz4OeACMMY8ASwHrgKKgGbgW8EKq/7NU1dEgjRx2DWJRLvD9IPKuHNplwTmdp12GZURDXut1Mnps6AbY27qY7kB/mfAEim/HCpew3CgOWFaVBR0n7g5GH8Rl6Vu5Me7DvOtWeF/qz2lAk1HioaphvJPaPM5caRNtjtKvymPn0OGo4raAxtp7fDaHUepkKMFPUy56rZS1D6StKRoOD7vVB5/CQAXJmzgt2/pvUaVOp4W9HBkDIM8O6mwJkRFh+hRrc4MKt3TmJu8kT2HG+2Oo1TI0YIehjwN+0my6mlPjrwZFvtSkTCXaXG7qa4u7buxUlFGC3oYqtjfOWVu4pBcm5P0vwMJc7DEMIG1VDXq7ItKdacFPQzVl3+C11hkjjrP7ij9rsY9mQYZxJykT1hTVGV3HKVCihb0MOSo28b+9kyyB0Xh9AkiHEycy+ykzawv1PFrSnWnBT0MDezYSYU1HsuKng7R7soT5pDoaKGx7H2dBkCpbrSghxlP82EGWocpMeNYvKHE7ji2OBh3AR24ybHW6DQASnWjBT3MlO/rnDK3MX6KzUns47XiOeA+n8uS1/Nxoc6+qNRRWtDDTM2BTwDwpUb2HOh9OZg8jxExh9hXtN7uKEqFDC3oYUZqtnKgPYPE5EF2R7FVecKl+BAG1r6Fx+uzO45SIUELephJb9/BPu/YqBoh2pNW5yBKZBqzE9byhxW77I6jVEjQgh5GPK31DLXKOOiI7HuI+utg8jwmx+3hSGWh3VGUCgla0MNI2b71WGJojI+eGRZ7U5k8D4DMppU2J1EqNGhBDyNHyjYA4Es5y+YkoaHBPZIy3yhynR/S1OaxO45StvOroIvIPBHZJSJFInJfD8tHiMj7IrJFRLaLyFWBj6rMkS0c8aTgTs7qu3GU2Bd7Kecm7GBToZ52UarPgi4iDuDPwJXAROAmETn+JO7/AV42xpwFLAAeD3RQBalt+ZQyFsvSD1ZH1aVfjVN8HNm91O4oStnOn8pwLlBkjNlrjGkHXgTmH9fGAMldj1OA8sBFVAAd7a1kOfbSFMUDinpSHzeNSu9ABta+ZXcUpWznT0EfDnSffLqs67XuHgC+0XUT6eXAd3takYjcKSJ5IpJXWakj/E7GC6vexi0e9nvH2h0ltIiQb83hbPdGyip19kUV3fwp6D1d8Hz8jEg3AQuNMZnAVcBzIvKFdRtjnjTG5BpjcjMyonCmwNPgqtsGgCclukeI9qQm7UrirDZ2bX3V7ihK2cqfgl4GdO+Fy+SLp1RuB14GMMasA2KBgYEIqDqltu2gyReLlaJH6MdrHXAhDb5EXBXL7I6ilK38KegbgTEiMlJE3HR2eh7/l1MCzAUQkQl0FnQ9pxJAw8xn7POMxrIcdkcJOUZc7I+5hMl8RENzi91xlLJNnwXdGOMBvgOsAArovJolX0QeFJFru5rdA/yHiGwDXgC+aXSi6oDp8HgY7dxDuUywO0rI2uWaywBnPTu3aeeoil5OfxoZY5bT2dnZ/bX7uz3eCcwKbDR1VPH+TznT0UK9e5LdUUJWR8ZltO5z075/CZz/ZbvjKGULvaA5DBwu6bwpdEeydoieiHEmssVzHuM7VuD1eu2Oo5QttKCHgY6qzXQYByZVJ+Xqzb74q8hwHqEof4XdUZSyhRb0MJDU8inFnmywYu2OEtI6Bl9Nq89NY+Fiu6MoZQst6CGuvcPLCNlNuYy3O0rIc8Ums803k+zmt8Gnp11U9NGCHuL2lhYy0FlLbYx2iPqjafD1pDuqObz3XbujKNXvtKCHuEPFXR2iOkLUL6Om3Uirz011wSK7oyjV7/y6bFHZp61yMwCeZJ2Uyx9ri9sY3DyDHPNm52kXHYiloogeoYe4hKbtHPRl4nUk991YAVDgupx0xxGayz+wO4pS/UoLeghr7fAyQnZRG6u3nDsZTRlX0OKLoTL/ObujKNWvtKCHsL+u3ESW+xDFPp2Q62QMSx/Eh425pFW/oVe7qKiiBT2EOes6z5+3Jus9RE+GwxK2O+aRTDWeCr3aRUUPLeghLKV5OwDtSVrQT1bzwHnUe+Op2vG03VGU6jda0EPYMF8+5Z5hdDjT7I4SdrIHD2RlwwWdp108zXbHUapfaEEPUa0dXsa4dlH2hftxK3+4HBalqV8mhmZ8pa/bHUepfqEFPUTtKtlHlvsQR2L0+vNTNXri1ZS3D6T+s4V2R1GqX2hBD1GH968BoC1lus1JwtclE4bwZv3FJNe8C616Ay0V+fwq6CIyT0R2iUiRiNx3gjZfE5GdIpIvIjrd3WnqqNwIQFvSNJuThK/EGCdlKV/Gwotv/4t2x1Eq6Pos6CLiAP4MXAlMBG4S+fyJXREZA/xvYJYxZhLwgyBkjSopzdso9WThcaTYHSVsLd5QQmvCRApasmne/YzdcZQKOn+O0M8Fiowxe40x7cCLwPzj2vwH8GdjTA2AMeZwYGNGl+Z2DyOtz7RDNAAmDE1mWd0cEhs3QX2h3XGUCip/CvpwoLTb87Ku17obC4wVkTUisl5E5vW0IhG5U0TyRCSvslLPaZ7I7v1FDHNXURs71e4oYS/G6eCzmKvxGcG7V6cCUJHNn4IuPbxmjnvuBMYAFwM3AX8XkdQvfJMxTxpjco0xuRkZGSebNWocLl4LaIdooGQOH8vaxqm0Fy0E47M7jlJB409BLwOyuj3PBMp7aPO6MabDGLMP2EVngVenoKpkLT4jtGuHaECMGZzIG41XEddeCofeszuOUkHjT0HfCIwRkZEi4gYWAMuOa/MacAmAiAyk8xTM3kAGjSaZvnxKvWfgsRLsjhIRnJZFacLl1HiS2Lv+EbvjKBU0fRZ0Y4wH+A6wAigAXjbG5IvIgyJybVezFUC1iOwE3gfuNcZUByt0JGto7WCcezcHLO0QDaTJWYNZWnMJI5pXQGuV3XGUCgq/rkM3xiw3xow1xow2xvy667X7jTHLuh4bY8yPjDETjTFTjDF60e8p2r1/F4NdR6iN1dMtgXRGejzvtH4JJx2w/3m74ygVFDpSNMRUdXWItmuHaECJCGnDp7OlaRxtu54Ec3y/vlLhTwt6iDHVeXiMRVuizuESaGeNSOPlmiuIaSqA6g12x1Eq4LSgh5i01m2UerPxWnF2R4k4iTFO9sRdTbMvFk/h3+yOo1TAaUEPIXXNbYxz7qTU0gFFwTI5O5NltRfB/heho8HuOEoFlBb0ELJ3z2ZSnY3UxJ9td5SINSojgY9883GaZsw+7RxVkUULegipLf0IgI7Uc21OErksES6YcTXbm8+kJf8R7RxVEUULeghx1nxCgy+B9vixdkeJaNdPz+SVxi8T37JLR46qiKIFPYQM6dhKqUwB0bclmF7dfIDKAfOp7EilYft/2x1HqYDRyhEiamqrGeXaR3PyOXZHiQpnjxzKP+uuJKHqLWjYY3ccpQJCC3qIKN3zIQ7xET90lt1RokKc20HHyDvxGouarQ/bHUepgNCCHiKayjvvIZo55mJ7g0SRW+bM5J2GC4ktWQgdjXbHUeq0aUEPEXF1n1DiySI5ZbDdUaLGgAQ3DWd8mzhpomzzE3bHUeq0aUEPAcbnY4RvO4fcZ9kdJaos3lBCS2Iu21vGYhU+pje/UGFPC3oIKC39lAHOOvaZKSzeUGJ3nKgS43aSF3sbw6xiCvIW2R1HqdOiBT0EHCpaBUBr6vk2J4lOzuyvcaBjCFbBbzE+PUpX4cuvgi4i80Rkl4gUich9vbS7QUSMiOQGLmLkk6qPqfEkY6XqTfYAfAAAABCmSURBVC3s4HS6Wev+FuOcO3nq1WfsjqPUKeuzoIuIA/gzcCUwEbhJRL5QeUQkCfgeoPOSnqRhbRsp8EzBshx2R4la7SNuo8qTxsTav+D16XQAKjz5c4R+LlBkjNlrjGkHXgTm99Dul8DvgdYA5ot4DTUlDHMeoNSlH2rsJM44Nsbcwsz4jaxev8LuOEqdEn8K+nCgtNvzsq7XjhGRs4AsY8ybva1IRO4UkTwRyausrDzpsJGodNc7ADQk6flzu9Vn3UGDLwFXwa9o9+i5dBV+/Cno0sNrxz6TiogFPAzc09eKjDFPGmNyjTG5GRkZ/qeMYO3lq2nyxuJM11vO2c3rSCEv9ptcFLeGdz5ebnccpU6aPwW9DMjq9jwTKO/2PAmYDKwWkf3AecAy7Rj1z4DmDezyTsbtjrE7igIqh/0njSaJ5MJf09zusTuOUifFn4K+ERgjIiNFxA0sAJYdXWiMqTPGDDTGZBtjsoH1wLXGmLygJI4g7c3VZMoeahJm2B1FdfE4UtgU9y0uSljPH19cbHccpU5KnwXdGOMBvgOsAAqAl40x+SLyoIhcG+yAkaz4s3ewxJCYNcfuKKqbQ0PvpN6XzKzWP1HX3GF3HKX85td16MaY5caYscaY0caYX3e9dr8xZlkPbS/Wo3P/1Je8S5vPybhJc+2OorrxWElsTbyT2Yl5vPWuHqWr8KEjRW2UXPcBO9onsXxnnd1R1HEqBv8Hld5BTDr8aw7XN9sdRym/aEG3SUt9BWNcReyx9HLFUOS1Ytmadg9T4gr5YOWf7Y6jlF+0oNukeOcbANQkz7Y5iTqRqoFfpULGcl7D/2PvoWq74yjVJy3oNmkvW0GtJ5GYQXp1Z6gy4mDzgPvIch/ik7d/jjE6JYAKbVrQ7WAMQ1s/Zmv7dNwut91pVC/qUueyndlcYz3Fu5u32B1HqV5pQbdBY1U+GdZhil0z7Y6i/LAn85e4LB/eTT+msU0HG6nQpQXdBsX5nVd7NqZdbG8Q5ZeWmDP4JPZ2rkh8n/9+/km74yh1QlrQbeCrWElZxxCSM8bbHUX56cCw73PYN4Sb5Lds2X/Q7jhK9UgLej/zdrSR7fmEEtdMHFZP856pUOS14tgy+DeMiS0l/52f0Obx2h1JqS/Qgt7PXl7+PEmOZgqderliuKlOnsNWx5e4MWERL6/8l91xlPoCLej9LKPubVp8MZjBl9odRZ2CPVm/pE2SmHrgXnaV19gdR6nP0YLen4xhMh+wqS0Xd2yi3WnUKWhzDGDDgF8wLX43H7/+Q71dnQopWtD70cHidQxxVrI3Vo/Ow9mhtOvYbF3FbcnP8eZ7r9odR6ljtKD3o4r8l/AZoX3QPLujqNNUNOK3VPsGkVP2XfZVVNgdRylAC3q/Sq5+m0/bJhKfPMzuKOo0eRwprBv0MJmuQxS/fRsdetWLCgFa0PtJWdlnjHbupsitc59HiubUmRQN/jEXx6zkw+UP2h1HKf8KuojME5FdIlIkIvf1sPxHIrJTRLaLyLsickbgo4a3vZufA6Bt8NU2J1GBNG7ub9jJBVxY/18UfPqO3XFUlOuzoIuIA/gzcCUwEbhJRCYe12wLkGuMmQosAX4f6KDhzBjD4CNL2OsdB8k6OjSSLP6kjM3DH6XSm076lpupqtxrdyQVxfw5Qj8XKDLG7DXGtAMvAvO7NzDGvG+MOXpbl/VAZmBjhreiwo2Mc+/myKCv2R1FBYEVm867GX8lXpqo/teVdLQ12B1JRSl/CvpwoLTb87Ku107kduCtnhaIyJ0ikicieZWVlf6nDHOHP30Kr7E4M/cOu6OoIHGmn8VL7t9xpquIoteuw3h1VkbV//wp6D1NONLjaAoR+QaQCzzU03JjzJPGmFxjTG5GRob/KcNYW0cHIxtfZ3N7LsuL7E6jgikuez6L2r/PBO975L9xG+gNMVQ/86eglwFZ3Z5nAuXHNxKRS4GfAdcaY9oCEy/8rVv/OsNchyhN+YrdUVQ/sMZ/j1XcxuTmxRS8/R2746go409B3wiMEZGRIuIGFgDLujcQkbOAv9JZzA8HPmZ4MsbQsvtZWnyxtA+5xu44qh9YlsVFX32K9zzXMaHmcT5d8b/sjqSiSJ8F3RjjAb4DrAAKgJeNMfki8qCIXNvV7CEgEfiniGwVkWUnWF1U2VRYzIXuVVSkfAmfI8HuOKqfLNl8gNJRf+DdlkuZUv0Qm9/8np5+Uf3C6U8jY8xyYPlxr93f7bFOTtKDT1Y/Sm5iC6sTbrY7iupnLpeL8glPsG7fDzm//k+sfaWRGdf/HYdDx/Kp4NHfriDZV9nIXOcr7PONpz7+LLvjKBs4HC72jHqEd9q/zMz2f7Bh8ZdobGnu+xuVOkVa0IPkrZULGRdbwv4B/wNE70wUrSzLQdX4/+ZffJuZjrcoeukCyiu/cE2BUgGhBT0IdhyoI6fxbxwxGRwacL3dcZTdRKgbfR+vxvyGia5P8fwrl4KCj+1OpSKQFvQgeHXlK8xM3E5R2h34xG13HBUiWjO/ztK0Z4mz2sjedClrVz2M0c5SFUBa0ANs3Z5qLm57jBYrjf1pt9gdR4UYb/osVmQtZz+TmHn4R2x64Woa6qNn1LQKLi3oAdTm8bJw6UIuStrCjpS78Vh6qaL6Iit+KJtHv8xr3v8kx6yg6bVJbNvwgt2xVATQgh5Af3lvN99Le5w6Gcze1FvtjqNCmFgumsf+lEXJL9Lii2Xanq+zZdEciks/szuaCmNa0AOk6HADVdsfZ1LcXnZk/B+8VqzdkVQYcA2awZoxK3nFezfjzVoGrj6Lt17+AfsO6mkYdfK0oAdAm8fLr5e8y08G/4Ny9wyKE3WYv/Kf5YyjbexPaLpsGyXu87jS8wjxK8bz+F9/wO/e3KYdp8pvWtAD4MFlO7jV+hXxTi95Q36v152rU/LO/ji2ZD/DsoGLqHGcwd1Jj/CNqst48ql7eCPvM1ra9b6lqnd+Df1XJ/ZyXimJ+/7IJUM3wdmP0ViTbXckFeYaUy5gc/Isihs/Yuzhh/hP98PUF/yVJesv5+CQ25kxdQbnj07HpdMIqOOIXR/ncnNzTV5eni3bDpQPdley6NXH+UvWg5QmXsnawX/Wo3MVcGktm8msfIrx7W/jEg8bGifzdtNcPMO/wiVTxzLrzIHEOB12x1T9REQ2GWNye1ymBf3UvL/rME/88xkWZt9PU8xo3sv8p16mqIIq1nOIkXUvkln3KhlmP+0+Jx80nM3HLedhhl7FzClTuXhcBrEuLe6RTAt6gL2yqYzl7/yDx7L+izb3UN4b/jKtzkF2x1LRwhjS2ndwRv1ShjUsJ9VUALCteQyftJxFx4CZpJ1xCTmjsxk3OAnL0k+NkUQLeoA0tHZw/9JtDC3/Ez8e8jw17rF8MOw5LebKPsaQ0r6LYU2rGFi3kqHeT3GKF58RdrWewda2KbQk5pA4ZDrDzjibyVmDSY3X6SjCmRb009TS7mXRhmLWrV/Gd9P+Qk78bnwjbmSJ4+d6mkWFFIevhQGtW0hqWE9a0wZG+LYRJy0AeIzF3rZMSn2jaYs/k7i0sWQMm8ywzMmkpg7D0k7WsHDaBV1E5gGPAA7g78aY3x63PAZ4FjgbqAZuNMbs722doV7QvT7D1tIaVm4poG3/K1yX+CbT4gtpkAy2Z/xfihOv1Q5QFfqMj6SOYuKbd+Bq+JSU1p0MNoUMsg5hyb//9hu9cVR602mUQbQ4B9PuHgpxw3AnDSMhcRDJKYNJThlEUvJgLHey/u7b6LQKuog4gN3AZXTeMHojcJMxZme3NncDU40xd4nIAuB6Y8yNva233wq6MYD5979dj9s9Xpra2mlq66C5uY7W5mqO1Byi5kgFLUc+w2ooZErsTibE7sMSw2HrTEoG3MKepK/hteKDn1upILJMG9JUTEftbmLa9pHsKSXRd5hkKhkgVWQ4qnBbnh6/12MsGk0yTSaJdomnnXg6JA6PxOOx4vBYCXiteDwSh0fiMJYbI2584sJYbsRyY8QFlhssF+Jw43TF4nbH4nbFEONy4na5iHW5iXE5EMuJZTkQsRDLgSUWxrKwxIFlOTtfFwdYDqyu18WysMRCBAQLSzrv9yoCIgJ0/YcUhv8x9VbQ/bkO/VygyBizt2tlLwLzgZ3d2swHHuh6vAR4TETEBON8TulSWHfLFwo0+D73mjEG4cSbd3d9pfW00AnNqQmUWxPZkjyfqsSLqI45KyzffKV64pMYSByLlTiWDjo/Vld3b2AMtFfja67A11qNaT+Co6MGp7cWt7eWGF8d8dQRSwsx0kqsVJEsrcRJC3HSSpzVSozVYc/OnSaf6fw7Nxz999/+/Zpw/MLjlx1r09XswYPfYVl9590677hgJD+6fFzAs/tzhH4DMM8Yc0fX81uAGcaY73Rrs6OrTVnX8z1dbaqOW9edwJ1dT8cBu04x90Cgqs9WkSXa9jna9heib5+jbX8hMPt8hjEmo6cF/hyh93RYevz/Av60wRjzJPCkH9vsPZBI3ok+ckSqaNvnaNtfiL59jrb9heDvsz/d2mVAVrfnmcDxN0U81kZEnEAKcCQQAZVSSvnHn4K+ERgjIiNFxA0sAJYd12YZcFvX4xuA94Jy/lwppdQJ9XnKxRjjEZHvACvovGzxaWNMvog8COQZY5YBTwHPiUgRnUfmC4IZmgCctglD0bbP0ba/EH37HG37C0HeZ9sGFimllAosHRqmlFIRQgu6UkpFiJAu6CIyT0R2iUiRiNzXw/IYEXmpa/kGEcnu/5SB5cc+f1NEKkVka9fXHXbkDBQReVpEDneNZehpuYjIo10/j+0iMr2/MwaSH/t7sYjUdXt/7+/vjIEkIlki8r6IFIhIvoh8v4c2kfYe+7PPwXmfjTEh+UVnB+weYBSdgzq3AROPa3M38ETX4wXAS3bn7od9/ibwmN1ZA7jPFwHTgR0nWH4V8BadYx3OAzbYnTnI+3sx8KbdOQO4v0OB6V2Pk+icRuT43+lIe4/92eegvM+hfIR+bMoBY0w7cHTKge7mA890PV4CzBUJ6/H5/uxzRDHGfEjvYxbmA8+aTuuBVBEZ2j/pAs+P/Y0oxpgKY8zmrscNQAEw/LhmkfYe+7PPQRHKBX04UNrteRlf/KEca2OM8QB1QHq/pAsOf/YZ4CtdH02XiEhWD8sjib8/k0hyvohsE5G3RGSS3WECpeuU6FnAhuMWRex73Ms+QxDe51Au6AGbciCM+LM/bwDZxpipwCr+/QklUkXae9yXzXTO1TEN+BPwms15AkJEEoFXgB8YY+qPX9zDt4T9e9zHPgflfQ7lgh6NUw70uc/GmGpjTFvX07/ROQd9JPPn9yBiGGPqjTGNXY+XAy4RGWhzrNMiIi46C9siY8yrPTSJuPe4r30O1vscygU9Gqcc6HOfjzu3eC2d5+ci2TLg1q4rIc4D6ozpuolmBBKRIUf7gUTkXDr/Rqt7/67Q1bUvTwEFxpj/PkGziHqP/dnnYL3P/sy2aAsTmlMOBJWf+/w9EbkW8NC5z9+0LXAAiMgLdPb4DxSRMuDngAvAGPMEsJzOqyCKgGbgW/YkDQw/9vcG4Nsi4gFagAVhfpAyC7gF+FREtna99lNgBETme4x/+xyU91mH/iulVIQI5VMuSimlToIWdKWUihBa0JVSKkJoQVdKqQihBV0ppSKEFnSllIoQWtCVUipC/H8nCs/Lhw95dwAAAABJRU5ErkJggg==\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(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
}