narrowlyapplicable/structural-changes-withCRP.ipynb

## structural-changes-withCRP.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ディリクレ過程に基づく構造変化推定（MLPノンパラベイズ本より）\n",
    "- MLPシリーズ「ノンパラメトリックベイズ」の§6の内容をPythonで実装した。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## モデル\n",
    "- クラスタ毎に線形回帰モデルを仮定し、この回帰係数および観測ノイズの生成過程にCRPを用いる。\n",
    "  - Gistだと下の数式が正しく表示されないかもしれません。\n",
    "  - 線形回帰モデル\n",
    "  $$y_t \\sim {\\cal{N}} (\\theta_t^{\\mathrm{T}}x_t, \\sigma_t^2) $$\n",
    "  - クラスタ毎の係数の生成過程\n",
    "    - ステップtにおけるクラスタの総数を$K_t^+$と置き、第kクラスタの出現回数を$n_{t,k}$と表記する。\n",
    "    - $k\\in\\{1,...,K_{t-1}^+\\}$、すなわち既存のクラスタに属する場合\n",
    "  $$P(\\theta^{(k)}, \\sigma_t^{(k)2}) = \\frac{n_{t-1,k}}{t-1+\\alpha}$$\n",
    "    - $k=K_{t-1}^+ + 1$、すなわち新規クラスタの場合\n",
    "  $$(\\theta^{(K_{t-1}^+ +1)}, \\sigma_t^{(K_{t-1}^+ +1)2}) \\sim {\\cal{N}}(\\mu,\\Sigma)IG(\\frac{n_0}{2},\\frac{\\tau}{2})$$\n",
    "  $$P(\\theta^{(K_{t-1}^+ +1)}, \\sigma_t^{(K_{t-1}^+ +1)2}) = \\frac{\\alpha}{t-1+\\alpha}$$\n",
    "    - $\\mu, \\Sigma, \\tau$については下記の通り事前分布を設定する。\n",
    "  $$\\mu \\sim {\\cal{N}}(\\mu_0, V_0), \\qquad \\Sigma^{-1} \\sim W(\\tau_0, \\Sigma_0), \\qquad \\tau \\sim Ga(\\frac{m_0}{2}, \\frac{\\tau_0}{2})$$\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. 準備\n",
    "### 1.1. ライブラリ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as st\n",
    "import seaborn as sns\n",
    "\n",
    "plt.style.use('ggplot')\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2. シミュレーションデータ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 「§6.4 実験例」と同じ条件で作成する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ y_t = 1 + 0.5t + u_t, \\quad u_t \\sim  {\\cal{N}} (0, 0.3)　\\qquad (1 \\leqq t \\leqq30) $$  \n",
    "\n",
    "$$ y_t = 25 - 0.3t + u_t, \\quad u_t \\sim  {\\cal{N}} (0, 0.1)　\\qquad (31 \\leqq t \\leqq 60) $$  \n",
    "\n",
    "$$ y_t = 1 + 0.1t + u_t, \\quad u_t \\sim  {\\cal{N}} (0, 0.2)　\\qquad (61 \\leqq t \\leqq90) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(1,91,1)\n",
    "y = np.zeros(90)\n",
    "y_tru = np.zeros(90)\n",
    "\n",
    "y[:30] = 1 + 0.5*x[:30] + st.norm.rvs(loc=0, scale=np.sqrt(0.3), size=30)\n",
    "y[30:60] = 25 - 0.3*x[30:60] + st.norm.rvs(loc=0, scale=np.sqrt(0.1), size=30)\n",
    "y[60:] = 1 +0.1*x[60:] + st.norm.rvs(loc=0, scale=np.sqrt(0.2), size=30)\n",
    "y_tru[:30] = 1 + 0.5*x[:30]\n",
    "y_tru[30:60] = 25 - 0.3*x[30:60]\n",
    "y_tru[60:] = 1 +0.1*x[60:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAFRCAYAAADgqHO9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8U1X6x/HPTZuWbqxVAVFRURERFINSRMZdQNxGvBXFBRdEFEEZ/LmM4z7juLEKyiAquNAj7vsyjuJS0ahYHUHHcUFBhU4pS1uatLm/P5JiWwrdc5Pm+369+tIkN8nT0xuePOece47lOA4iIiKJwON2ACIiItGipCciIglDSU9ERBKGkp6IiCQMJT0REUkYSnoiIpIwlPQkoVmW5ViWNSZK7/WDZVl/boHXOd+yrIqWiKmJ7/+wZVlvuvX+Is2hpCeJrhuwxO0g6mJZVo9IUj6y1kN5wK4uhNRklmVVWJZ1vttxiCS7HYCImxzH+dXtGBrLcZwyoMztOETikSo9adMsyxpiWdb7lmVtivx8blnWCdUer9G9Gbk90bKsPMuySizLWmVZ1ijLsjpYlvVY5DW+syzr9GrP6Rl53pBa7/2tZVk37SC2syzLWmZZ1gbLsgoty3rJsqx9qx3yU+S//4q8/g+R523TvWlZ1gjLsj6xLKvcsqy1lmXNsSwro9rjD1uW9aZlWeMsy/rRsqyNlmU9Z1nWTvW0X6dqbfGbZVm3AVatY46zLOtty7KKIr/LO5ZlHVrt8R+AJOChyO/hVHvtRyNtXGZZ1teWZU2xLKvG64u0JCU9abMsy0oCngeWAQMiPzcBpfU89XrgZaA/8CKwEFgMvAEcDLwELLQsq0szQ0wFbo3EdRxQCbxkWVZK5PEBkf+eTrgbdmBdL2JZVj/Cv+dS4CDgPGAkcH+tQwcCRwEnAsMix95dT4wLgEOAk4CjgZ7AabWOyQTuAwYBg4H/AK9Wa5+Bkd9tcuT36Fbt9/8COBXoQ7gtbgbOrycmkaZzHEc/+mmTP0AnwAGO3MExDjCm1u3p1W7vFLlvVh2vOzJyu2fk9pBar/0tcFO12z8Af95BLJ0jr3N45HaPuuInnBQqqt1eBHxU65hTgBCwR+T2w8A6ILXaMdcAv+wgnl6R9z+u2n0pwGrgzR08zwOsB86udl8FcH4D/mYzgDfcPnf003Z/VOlJm+U4znpgPvCaZVmvWJZ1jWVZ+zXgqZ9Xe411hKuUglqvGwB2bk58lmUdZFnWM5ZlfW9Z1iZgVeShPRr5UgcQrvKqe4dwN2SfavetcBynvNrt1cAuO3jdqud+UHWH4zgB4OPqB1mWtadlWYsi3bkbgY1AB+r5PSzL8kT+Jssj3bubgfH1PU+kOZT0pE1zHOdiwt1zbwB/AL60LOuSep4WbMB9Dr9/fkKR/9Yei/Ju7w0sy0oHXo+8zgXAoYS7AR3C1VRjbW+7lOr3B+p4bEfjZw0dW3sR2B24jHAX50HAWur/PaYA1wKzCHfvHkT4S0pTfn+RBlHSkzbPcZwvHce513Gc4cCDwLgWfot1kf92r7rDsqyd2fFlBfsT7jq93nGcfzmOs4Jwt2n1RFOVpJLqef9/E07o1f2BcFL7qp7n1ve6EB6nAyAy3jiw2u0uhCvCOxzHec1xnK+ALWxbBQfY9vcYCrzqOM6DjuN85jjOt8A+zYhXpF5KetJmWZbVy7Ksv0dmcO5hWVYOcATNSwTbcMKXELwPXG1ZVn/Lsg4hPPmlfAdP+zHy+ETLsva2LOsYwuNZ1SuzQmAzcLxlWV0ty+q0nde6CxhgWda9lmX1tixrGOHq6THHcVZt5zkN+b2+JTxB5j7Lso6yLKsP4Uosq9ph6wkn/Ysty9o30sZPsO0lFd8DR1mW1d2yrOzIfV8DR0Zee9/IzNDDmhqvSEMo6UlbVkK4clgMfAM8RXh86vJWeK8LCCeoDyLvNw/4ZXsHO45TCIwh3K33b8KzKP/E712lOI4TItxlaBO+fOGz7bxWAXAy4eruc8ITW14iPD7WXBcAywl3Yb5DeBzwmVoxngHsTXjc82FgOtv+7lMIdzN/z++V8a2R13wOyCdc6c5sgZhFtstyHO2cLiIiiUGVnoiIJAwlPRERSRhKeiIikjCU9EREJGEo6YmISMKIx62FNN1URETqUu8qQvGY9FizZk2jn5OdnU1hYWErRJMY1H5Np7ZrHrVf8yRK+3Xv3r3+g1D3poiIJBAlPRERSRhKeiIikjCU9EREJGEo6YmISMJQ0hMRkYShpCciIglDSU9ERBKGkp7ENb/fz6xZs/D7/W6HIiJxIC5XZBGBcMLLzc0lGAzi9XrJy8vD5/O5HZaIxDBVehK38vPzCQaDVFZWEgwGyc/PdzskEYlxqvQkbuXk5OD1egHwer3k5OS4HJGIxDolPYlbPp+PvLw88vPzycnJUdemiNRLSU/ims/nU7ITkQbTmJ6IiCQMJT0REUkYUenetG17ATASWGuM6Vvt/onA5UAF8JIx5upoxCMiIokpWpXew8Cw6nfYtn0UcArQzxhzAHB3lGIREZEEFZWkZ4xZChTVuvtS4A5jTHnkmLXRiEVERBKXm2N6+wJH2La9zLbtd2zbHuhiLBKn3nsvhSFDdubzz71uhyIiccDNSxaSgU7AIGAgYGzb3ssY49Q+0LbtccA4AGMM2dnZjX+z5OQmPU/CYrH9HAfuvDOZ77/3cPHF2XzwQZCuXd2Oalux2HbxRO3XPGq/mtxMej8DT0eS3Ee2bYeAbGBd7QONMfOAeZGbTmFhYaPfLDs7m6Y8T8Jisf3eeiuVTz7pwoUXbubxx9M57TRYsqSQ1FS3I6spFtsunqj9midR2q979+4NOs7N7s1ngaMBbNveF0gB2v5fRlqE48C992bRo0cFf/7zRmbMKObTT1O45pqOONv0FYiIhEXrkoUngCOBbNu2fwZuBBYAC2zb/hIIAOfV1bUpUpd3303ls89SuOOOYlJS4MQTt3DVVZu4994s9t8/yLhxJW6HKCIxKCpJzxgzejsPjYnG+0vbEq7yMunWrRLbLt16/5VXbmLFimRuvbU9++5bwZFHlrsYpYjEIq3IInHngw9S+PjjVC67bFON8TuPB8477y26dPmVceOy+O9/k9wLUkRikpKexJ3p07PYZZdKRo8urXG/3+/n/PNH8b//DaGkZAOjR6fzzjvLtbO6iGylXRYkrixblsIHH6Ry000baNeu5mNVm8qGQt/h8eSyZs3rjBnjAe4hJSWJvLy8rcdpKyKRxKSkJ3HllltCpKdvok+f94CDazxWc1PZfPr3X8SyZWOB2wgGr2PJkiU8+eSTBINBvF4veXl5SnwiCUbdmxI3Hn30vyxfvjNlZbdx7rmjtumyrNpUdurUqeTl5XHddV1ISrofuBqP5zwAgsEglZWVBINB8vPzXfgtRMRNqvQkbjzwwM7AOhznvq1Jq3alVntTWWM+4cor/8svv8zjwAPfx+t9EgCv10tOTk40wxeRGKCkJ3Fh+XIv3323H8nJN+A4WxqctAYNOoSXXvJw4okO99xzOHPmPMc33/xLY3oiCUpJT+LC9OlZdOwYYu7cQXz++dRGJa3OnUM89FARJ5+czcyZR7JkSV/S0sKP+f1+TWwRSSBKehLzvvwymTfeaMfUqRsZOvRghg49uP4n1dK7dwWzZhVzwQWdufrqjsycWcwnn/jJzc3VxBaRBKKkJzFv+vQs2rcPccEFzVta7IQTtjB16kbuuqs9ffoEqajI3zqxBahzjFBE2hbN3pSY9tVXybzyShoXXlhC+/bNX5p10qTNjBxZxu23t8frPQWv10tSUpImtogkCFV6ElNqj7HNmJFFZmaICy/c3CKvb1kwbVoxP/zQhenTD+Xuu1/i55/f0JieSIJQ0pOY4ffXHGO7666XeemlI7n88s106tRyG3CkpzssWLCeESOyufvuobz00v507KgNPkQSgbo3JWZULSNWdfH43LmdSEtzWmWboF13rWT+/CJWr07i0ks7UVERvt/v92utTpE2TJWexIzqy4glJfVh5cr+XHJJCZ07h1rl/QYODPL3vxdz1VWduPXW9px00luazSnSxinpScyoWkYsPz+fZcsmkJ/vMH58y4zlbU9ubhlffeVl/vxMVq8OaDanSBunpCcxxefz0aXLYdx1185ceGEJ2dmtU+VVd8MNG/nmm2TeeOOPJCUdAbyr2ZwibZTG9CTmzJ6diddLq1d5VZKTYc6c9fToESI9/RUuueQ2dW2KtFFKehJTfvopiSVL0jn77BJ22aX1q7wqnTo5PPxwERUVKbz77lUccMDAqL23iESPkp7ElFmzMvF44NJLo1PlVbfPPhXcd996vvzSy1VXdcTRVQwibY6SnsSM1auTMCad0aNL6dYtelVedcceW861127i+efTmDUr05UYRKT1aCKLxIz77gsnmcsui36VV92ECZtZsSKZv/+9Pb17Bzn++HJX4xGRlqNKT6Jqexd///KLhyeeSMe2S9l110qXoguzLLjrrmL69w9w+eWd+PprfTcUaSui8mm2bXsBMBJYa4zpW+uxPwF3ATsZYwqjEY+4o/YyY9VnSM6dm0koBJdf7m6VVyUtDR58sIjjjuvIH/+YxKxZn3D00f3dDktEmilald7DwLDad9q2vRtwHLAqSnGIi2ovM5afnw/A2rUeHnssg1GjStl9d3ervOpWr/6IkpLjKS7O4Lzz0vjww0/cDklEmikqSc8YsxQoquOhacDVgObJJYCqZcZqb+Uzd24mwSBMnLhtlefmWpj5+flUVn4AjCMUOpLbb+8S9RhEpGW5Nlhh2/bJwGpjzOe2bbsVhkRR9WXGqrbyKSz0sHBhOqeeWkbPnjWrvB11h0bD72uBPg4czKefXsntt39O+/aLtRWRSJxyJenZtp0OXA8c38DjxwHjAIwxZGdnN/o9k5OTm/Q8CWup9hs2bBjDhv3e033vvUmUl1vceOO2r19QUFBjLcyCgoIaz21tw4YN47XXXmPp0qUcfvghXHNNEXPm9MHj+ZDU1Bm8+uqrDBo0qN7X0bnXPGq/5lH71eRWpbc3sCdQVeX1AD61bftQY8yvtQ82xswD5kVuOoWFjZ/vkp2dTVOeJ2Gt0X5FRR7mzt2ZU04po0uXYmq/fL9+/bbuuuD1eunXr1/U/4a9evWiV69eAAwePIuPPsolFDKUlw/ilVde2frYjujcax61X/MkSvt17969Qce5kvSMMV8AO1fdtm37B8Cn2ZuJZd68DMrKLLp0eQC/f+9tugvr6g5101FHHcycOWcQCLyL4zzDQQf9VO9z/H4/BQUF9OvXz/X4RQQsJwprLdm2/QRwJJAN/AbcaIx5sNrjP9DwpOesWbOm0TEkyred1tLS7VdcbOHzZVNe/gKWZcfN/nV+v5/HHitkyZKxDBu2hQceWI9nO9PB3B6TbCv02W2eRGm/SKVn1XdcVCo9Y8zoeh7vGY04JHY8+GAmZWVePJ5b42r/Op/Ph88HvXtv5JZbOjB9epCrrqr72sLql2hU3Y7130+krdNSExJ1GzdazJ+fwaBBv7B8+UqCwaS4279u3LgSVqzwcs897endu4IRI7YA4equqju2+k7w8fb7ibRVSnoSdQsWZLBxo4ebb05iy5bYGbNrDMuCO+4opqBgCxMmZHLXXQXsuefGbboz8/LyNKYnEkOU9CSqNm+2+Mc/MjnuuC307VsB+OI2GXz5pZ8ffphEMPgukyfvyRln3LlNd+bEiRMZNmxYQoypiMQDLTgtUfXwwxkUF3uYPHmT26E0W35+PhUVPwGnAjvzzjsTSU5O32bFGRGJHar0JGpKSiweeCCDo47awkEHBd0Op9l+H7NbjscznrVrH+GEEwo46KAHGDw4vrprRRKFkp60quoTO/z+P1BUlMTkyXUtwxp/al9H+MYbm5g9uydDh/4fPl+p2+GJSB2U9KTVVL9OLTm5Pe3a/cIRR5Tj88V/lVclfAlDuKIbMGATK1d6+ctfOrDPPhUcfnjA5ehEpDaN6UmrqX6dWiBwPhs2pHLllfE/lrc9Hg/Mnr2evfaq4JJLOvHjj0luhyQitSjpSaupGvPyeNJxnKn07VvIYYe17eonK8vhoYeKcByLCy7ozObN9S4QISJRpKQnraZqzOuYY54AuvGXv7gdUXTsuWcl999fxH/+k8wVV3QkFHI7IhGpoqQnrerAA318+eWJHHpoOYMHt+0qr7ojjghw000bee21NG6+Wd2cIrFCE1mkVRmTzi+/JHHPPcVYCdbTN3ZsCStWJHPHHRnsvns7Tjlli9shiSQ8VXrSagIBmDUrkwEDAgwdWu52OFFnWXD77Rs4/PAQV13VkS++8LodkkjCU9KTVvPUU+msXp3MlVduSrgqr0pKCixeXEHnziHGju3MunX6yIm4SZ9AaRXBIMycmUn//gGOOirxqrzqdt4ZHnqoiPXrLc48M4Vp0+bg9/vdDkskISnpSat45pk0Vq1KZvLkxK3yquvbt4IrrviElSs7c889e2HbuUp8Ii7QRBZpUX6/n/ffX8ajj17LAQcEOe64xK7yqvN4nsayXsNx/kwgUKBNZUVcoKQnLaZq2bHy8lE4TibnnvsRltXD7bBiRk5ODikpZ1Je3hfHuYe0tI/cDkkk4ah7U1pMfn4+gUAljnMdUIDjPOt2SDHF5/NhzGKuuuoz9tijhHvvPYzvvtM1fCLRpKQnLSYnJ4ekpFxgf7zeOxk8eJDbIcUcn8/HlCmXsHhxGR6Pw9ixndm4UYOeItGipCctZsAAH9263Ud29m8YY2u8agd2372SefPW88MPyVx2WScim62LSCtT0pMW88or7Vi1qj0335zCoYcq4dVn8OAAt966gbfeascdd2S5HY5IQlDSkxbhODB9ehZ77VXBSSeVuR1O3OjTZykHH/whc+Zk8dRTaW6HI9LmafamtIg33kjlq6+8zJixniTNzWiQqtmugQBY1mtMmXI4e+9dwUEHtZ1NdkViTVSSnm3bC4CRwFpjTN/IfXcBJwEB4L/AWGNMcTTikZblODBtWhY9e1Zw6qmq8hqqapPdUKgSj8cmPf0rzjknk9Gjp3H88X01JirSCqLVvfkwMKzWfW8AfY0x/YBvgGujFIu0sLfeSqWgIIUrrthEsvoOGqxqk92kpCRSUjYwduxzFBVVct99x2Lb52rFFpFWEJWkZ4xZChTVuu91Y0xF5OaHgK5ijkNVVd5uu1Xwxz+qymuMqk12p06dSl5eHu3afYPHcx5wKIHALD74IN/tEEXanFj5Xn4BkOd2ENJ477yTymefpfD3vxfj1c45jebz+Wp0Y6akTKe8/EYc52Z+++3fLkYm0ja5nvRs274eqAAe28Ex44BxAMYYsrOzG/0+ycnJTXqehNXVfo4Ds2cns9tuDpdemk5qarpL0cW2hp57w4YN47XXXuOdd5by1lvrWLiwD6edVsGwYU4Uooxd+uw2j9qvJstxovOBsm27J/Bi1USWyH3nAeOBY4wxpQ18KWfNmjWNfv/s7GwKCwsb/TwJq6v93nsvhdzcbG6/vZjzz2/ony/xNOXcKy21OPXUbFatSuLFFwvp1aui/ie1UfrsNk+itF/37t0B6l3eyLXr9GzbHgb8H3ByIxKexJDp07Po2rWSM8/Un6+lpac7LFhQREqKw/nnd6a4WEuVibSEqCQ927afAPKB/Wzb/tm27QuB2UAW8IZt28tt274/GrFIy8jPTyE/P5WRI7/iH/+YpZmGraBHj0rmz1/Pzz8ncdllnahI3GJPpMVErXuzBal70wW12y83twtffulQVtaNiopNeL1e8vLydG1ZHZp77j3+eDpTp3Zk3LjNnHjiW+Tn55OTk5Mwba3PbvMkSvs1tHvT9YksEn8+/jiF995L5eijX+KddzZRGVktWZuito6zziplxYpk5s3LZMGCF3Cch/QlQ6SJtPamNNr06Zl06VLJJZew9eJqr9dLTk6O26G1WTfeuJE99viWiorZVFYOJBgMkp+v6/hEGkuVnjTKZ595efvtdlx//UaGDDmYvLy8hOtuc0NyMtx2238491wPjvMMyclD9CVDpAmU9KRRpk/PIisrQGnpPfj9h2xzcbW0nqOP7s/MmV8zZUoOPXr4OeCAMiDuxuRFXKXuTWmwL77w8uab7Sgru42ZM28nNzdXszaj7I9/3I8HHijhu+86MGVKB+JvHpqIu5T0ZBt+v59Zs7a9DGH69ExSU8sIhWZRWVmpcSWXHH98Of/3f5t47rl0Zs/OdDsckbii7k2poWqPt2AwWGOGYEGBxauvpjF69EqeeWYLwaAmr7jp8ss3s2JFMn//exb77Rfk+OPL3Q5JJC6o0pMaqvZ4q13J/e1vSWRlhbjhhg41dgbQeF7r2l7VbVlwzz0bOPDAIJdf3omvv9b3V5GG0CdFaqja4w3YWsmtXJnM0097mDx5Ex06OJq8EiXbq7qrpKU5PPhgESNG7MTYsZ158cV1dO6sQT6RHVGlJzXU3uPN5/Mxc2YmmZkOF1202e3wEsr2qu7quncPMX9+Eb/8ksT48Z0JBl0IVCSOKOnJNnw+HxMnTsTn8/Htt8k8/3waEyaE6NRJVUQ0Vd9ZfXvjp36/n/z8e7n00uW8/34qt9zS3oVIReKHujdlh2bMyCQtzWHSpEq3Q0k4VVX39i7+r9n9OZ1TTvmSBQv2Zs2a17n00mR1QYvUQUlPtuu775J49tk0xo0rITs7hQRYszbm7Gj8tHr3J0Bm5q14PGfx6qsn8dZbw3nyySuU+ERqUfembNesWVmkpMD48RrLi0W1uz89nhAwGviOQOBxXn11hdshisQcJT2p048/JvHUU2mMGVPCTjuF3A5H6lB70tGoUaNISSnF4zkNSOHVV8dTWqrNZ0WqU/em1Gn27EySk+HSS1XlxbLa3Z9VY4ApKSu59dbDOO+8Eo44YhaDB2tBcIk94YlY0V2wXklPtvHzz0kYk84555TQtauqvHhSPQmuWvUVDz98APn5qaSm5pKXlwegXTEkJtR3HWprUdKTbcyenYnHAxMmqMqLZ7vs8iiW1QfHuYlA4AuWLFnCk08+GfV/ZETqUnsiVrQ2odaYntSwZo2HvLx0cnNL6d5dVV48Gzw4h5SUicBHhEKPsH79bvVe7C4SLQ25DrU1qNKTGubMySQUCi9oLPHN5/NhzEJef/1NHn+8Hx99dB3JyY8Av2qxcHFdfdehthYlPdnqt988PP54BrZdSo8euhi9Laga4xs+fBOnn57NPvssZ/jwaQwZcpi6NsV1bqzjq+5N2Wru3EwqKlTltUUHHxzkrruK+fLLbH799TolPImK7e0S4iZVegLAunUeFi1K5/TTy9hjD1V5bdHpp5excmUyc+Zksf/+Qc47r9TtkKQNc2t2Zn2ikvRs214AjATWGmP6Ru7rDOQBPYEfANsYsz4a8UhNfr+fv/61A4HALkycuKnOxwsKCujXr19MnLTSdNdcs4mVK7385S8d2GefCgYPDrgdkrRRbs3OrE+0ujcfBobVuu8a4J/GmH2Af0ZuS5T5/X5s+zKWLfMBT1BUtGybx3Nzc7npppvIzc2NqW4KabykJLjvvvX07FnBuHGdeOmlr2Ku+0niV/XuTLdmZ9YnKknPGLMUKKp19ynAI5H/fwQ4NRqxSE35+fkEApcBacDt20xjb8iebhJf2rd3eOihIgKBSi65pCt33jlXX2ik2aq+IN91113k5uYCbLM3ZyyM8bk5preLMeYXAGPML7Zt7+xiLAmrb9+hOM4QLOtJUlK+2+bbWF07qUv822uvSkaOXERe3vk4zkMEAvbWLzRasUWaoq7uzKp9OWHbMb6//vVVPvnEx5Ah5Zx88paoxRkXE1ls2x4HjAMwxpCdnd3o10hOTm7S89q6r746HkjissuKOeOM1xg0aFCNx4cNG8Zrr73Ge++9x5AhQ7Z5XOoXq+fehAn78PTT/0cweDcezy3stlsHzjzzTAKBACkpKbz66qsx8feO1faLF41pvw8//JClS5cydOjQRv/thw8fzowZM7aeP8OHD6/xvgUFBQQCyYRCZ1JZOYGrrhpIWppDnz6pZGdnNuq9msPNpPebbdvdIlVeN2Dt9g40xswD5kVuOoVN2NgtOzubpjyvLduwwWLWrF0YMaKMa689CaDONurVqxeDBg2isLBQbdgEsXru9erViyefHMpf/vIxBQXX8sYbjxEIBKisrCQQCPDKK6/Qq1cvt8OM2faLFw1tv+bOtuzVqxeLFy/e2lPQq1evre/7zTfJLFt2FqHQ5UBHLOsrLrroC668MpsOHZwW2auze/fuDTrOzaT3PHAecEfkv8+5GEtCWrAgg02bPEyatO2MTUkMAwf6ePZZsO0AL7+cS1LSLMCvruwE1BKzLatfbF5eDq+8ksaiRel8+GEqXu9ODB26mh498rDtbgwc6AOclv416hWtSxaeAI4Esm3b/hm4kXCyM7ZtXwisAs6IRiwS/kb39tufMm/edZxwQhl9+1a4HZK4KDUV5s8vYvjwnQgG32L06Gkce+yBGtNLMC01fv/jj0k8+mg6ixenU1SURM+eFVx//UZsu5TsbA9wUgtG3XhRSXrGmNHbeeiYaLy//K6qC6O8fAqOk8Lxx78D7ON2WOKynXYK8dBDRZx6ahc+/PBPXHnl/9wOSaKsOWthVlTAm2+2Y9GidN5+ux1JSQ7HH7+Fc84p5YgjyvHE0NpfcTGRRVpO+BKFVBznSuBl1q17FyU9ATjwwCD33lvMhAmduf76Dtx11wYsbbyeUOpbC7P2pq9r1nh44okMHn88nV9/TaJbt0qmTNnI6NGldOsWm7u0KOklmJycHDweD6FQF1JS/k5OzhS3Q5IYcsopW1i5chMzZ2bRp08FF1xQ4nZICc2NncV3FEtubi6BQAVJScMZMGAgH3+8C44DRx1Vzt/+VszRR5eTHONZJcbDk5bWp89AMjJOoHPnb5g+fYrrHySJPVOnbuLrr5O56ab29OoVZOhQLVXmhlhYu7J60n3zzS8oL78Sx7mIUGgvvvhiMxMmbOass0rjar1eJb0Es2hROhs2pPLww52V8KROHg/MnFnMySe9ArmEAAAgAElEQVRnc+mlnXnxxXXsuWf8/KPWVri9dmV4icJcgsHBwAAs6xocJwl4B6/3RhYutPF6HZ5/PjYq0YZS0ksgZWXh7YMOP7ycQw/Vt3fZvszM8FJlI0bsxNixnXnhhUKysqI/vTyRtcZqSA3tLi0utpg500N5+afA/sB6fL58xo4N8NNPb5CTcwbguF6JNoWSXgJ54okM1q1LYu5cbWYh9dtjj0rmzSti9OguXH55JxYsKCIpye2oEkdL7yxeX3ep48Bnn3lZtCiD559PY8uWkVjWMuACUlKe4YYbHokc3xuAWbNmxeQuCvVR0ksQW7bAffdlMmhQOTk5qvKkYQ4/PMAtt2zg+us7cuedWVx7rRYyiKbG7iy+o0pue92lmzdbPPNMGosWZfDvf3vJyAhxxhmljBlTwpYtv5Gfvws5OY9s83rxui6vkl6CyMsLTymePl1VnjTOeeeVsmKFl9mzs9h//wpOPbXM7ZCkDvVVcrWTVNeuJ3DNNR14+uk0Sko89OkT5I47ijnttDIyM6u6srefdFu6Eo0WJb0EEAjA7NmZ+HwBhgxRlSeNY1lw660b+M9/kpkypSN77llB//5Bt8OSWuqb+OLz+Vi48EkeeaSMb789hsmTO9OuncPIkWWce24JAwYEG31dZmMr0VigpJcA7r57LWvWdOfii/OxrD3cDkfiUEoK/OMf6xkxIpsLLujMyy+vY5ddYvPi40RVV3ej3++noKCAzp0H89lnh7JkyQiKiz3svXeQm27awKhRpXTqlFgTlCzHadgvbNv2vcBCY8zy1g2pXs6aNWsa/aREXan9ww8/YdSoA3GctaSm/gFjmjbDKlHbryXEc9vVHiP697+TOeWUbHr3rmDJkkLatWv9GOK5/aKt+t8rGLQYPdoQDF4IHElycogRI8oZM6aEwYMDbW61ncguC/X+Vo2p9LzAa7ZtrwMWAY8ZY35uWngSLQ8+uAXH6QlcTkVFMG5mWIn7tjdGNGNGMePGdeaaazoybVpxm/vHM575fD523vkwHnssnQULIBg8Cfgey7qOSy/N4JprxrodousavAyoMWYi0B24BjgIWGHb9pu2bZ9r23b0dgCUBquogE8/HYFlfYrH82pczbAS91UfIwoGg1t3Vj/xxC1MmbKRJ59MZ968DJejFAh/1l9/PZVzzunM4ME7M2dOJgceuBmv92Q8nn1JTQ3vnCGNHNMzxlQCLwIv2rZ9APA48DAwx7btxcCNxpjVLR6lNMmzz6bx668ZXHddBaHQ1LiaYSXu29GU9MmTN7NihZfbbmvPypVPcfbZ2Tq3oqB2d/Ovv3p44ol0Hnssg19+SaJr10omT97M6NEl7Lqrg99/EQUFh9KvXz/9fSIaPKYHYNt2e8L73o0B+gFPAY8Q3g9vCnC0MaZfK8RZncb0GqCyEo48cmfatXN4/fV1ze6CSrT2a0nx3HY7uu7r3Xc/Y/To3XCc3UlJGcqTT97WKv+wxnP7taSaCz6fwCGH/IOPP+5KZaXFH/4Q3sbnuOO2bLPgc6K0X4uP6dm2vQQ4AVgK3A88a4wpr/b4VcCGRkcqreKFF9L47rtk5s0r0piLNFntKenVk+Dy5flYlsFxPiQQMPzrX/GxDFUsacwuCv/8ZwHl5VfgOBcTCvXi8883c8kl4QWftTZqwzWme/ND4HJjzK91PWiMCdm2vUvLhCXNEQrB9OmZ7LdfkOHDt7gdjrQRtSe23HzzzaSkrCEQsAmFXmfp0vFcdVVlXC9V1tpb+VR/faDetSsdBz76KIVFi9J54YVrIws+v4vXeyuLFp1BTs6AFo+xrWtw0jPG3N2AY0qbF460hJdeasd//uNlzpyimNqxWOJb7Yuf169fv3VFjo0bv2LOnP789a+bueGGjS5H2jStvZVP7dc/44wztnsx+YYNFk89lc6iRel8842X9u1DnHdeKQcd9BGrV79BTs7p+HxKeE2hi9PbmFAIZszIolevICNHqsqTllPXxJbq3Z9lZZu5//5MevcOcsYZ8bdUWWtv5VP79YFt2vPzz70sXJjOs8+msWWLh4MOCnDPPes55ZQtpKU5hBd77t1iMSUiJb025rXX2rFihZeZM9fHdTeTxJ761lq88caNfPONl6uv7shee1VwyCHxtVRZa2/lU/v1R40axahRo3jnnU8oLz+NG27oQ0FBCunpIU4/vYxzzinlwAPjqw3jQaNmb8YIzd7cDseBE07YiaKics45528cfvhhLfZNNRHar7UkUtu99dbnTJx4KJaVwRtvFNOtW/OXKotm+zV2TG9Hx9fVXQpsPT4jYxCPPprBU0+lsWmTh969g5xzTgmnn17WonsXJsr51xorskiMe+ONVP79by/JyRO4556HmDkzfjZ2lPjn9/u5+OJcAoF9CIXe58wzU3n11TLS0tyOrOEas4ByfWOAdXWXXnzxRH788Qhuuy2djz9OJTXV4cQTwws++3yNX/C5tSfetEVKem2E44TH8jp0KGLTpkWEQvG1saPEv6p/5EOhL/B4zuHbb5/m6qtTmTmzbS5VVt8YYPXuzKSk3qxceSGHHNKV4mIPe+5ZwQ03bMC2S+ncuWlVXWtPvGmrNLevjXj77VSWL09hzJhVpKRYJCUladkxiaqqf+STkpJISXmNMWO+5umn0xkz5nP8fr/b4bW46r9vXZ+1/v19TJq0lB49VhIIfMmLL+7NkCHlLF5cyNKlaxk/vqTJCQ+2v0yc7JjrlZ5t21cCFwEO8AUw1hijaYeN4Dhw771Z7LprBX/6084cf3z8bewo8a/2RBfH+Q+PP/4lb799Ou+/fzpLltCmzsftTexZvTqJxx5L54kn0lm7tju77lrB1VdvZPToUnbeecdjnI3prozXncvd5mrSs217V+AKoI8xpsy2bQOcSXg9T9mO2h+Md99N4dNPU/jb34pJSYnPjR2lbah+7s2aNQuYDexFMPgIL7wwH6BNfSGr+n0rK+HNN1NZtCiDt95KxXHgmGPKOeecYo46qrxBM6kb210ZrzuXu831So9wDGm2bQeBdKDxUzMTSO0PxuLFeUybdgJdu1aSm6u1ASR25OTkkJIynUDgdEKhZTzzzFgWLepDRcU618agmjvxo/bz166tWvA5ndWrk9lpp0ouu2wzY8aU0qNH45YGa8p1gvqC23iuJj1jzGrbtu8mvGB1GfC6MeZ1N2OKdbU/GHl5v/LRR6nccssGUlNdDk6kmuqVSMeO33DddTmEQguB4UD093Zs7sSPmgs+H8vAgQ/y0Ue7UFFhMWRIOX/5y0ZOOGELkR7HRlN3ZXS43b3ZCTgF2BMoBp60bXuMMebRWseNA8YBGGPIzs5u9HslJyc36XmxZvjw4cyYMYNAIEBKSgorV9p07epwxRVppLXi3PC20n5uSOS2GzZsGMOGDQOgqOi/3HnncVjWPaSkXMvw4cMb1C4t1X4FBQU1vjAWFBRsja0hPvzwP5SXXx5Z8HlfPvuslMsvD3HhhZXsu68FZEZ+mmbYsGG89tprLF26lKFDhzJo0KAmv1Z1iXz+1cXt7s1jge+NMesAbNt+GhgM1Eh6xph5wLzITacpF1q2lQs0e/XqxeLFi8nPz6dDhxO59tqO3HjjBkpKSigpab33bSvt5wa1XdikSVmsWPFfXnhhEhddNJRevXZpULu0VPv169evRiXVr1+/el/XccDv97JwYQYvvDA5suDzB3i9f2PhwtMZPDi8/mVL/Xl79epFr169Iq/ZMi+aKOdf5OL0ermd9FYBg2zbTifcvXkM0PbmNrewqn78s87qTHZ2Jeeco7E8iQ+zZ6exYcMW7r//II455n8MHBiI2ns3ZuLHpk0WTz2VxqOPZrBihZesrBBnn13GwQd/xC+/vE5Ozml1Lvisi8Vjn9tjessi+/R9ClQAn/F7RSc78MknXt55px1//vOGyEK0IrEvORnmzl3PiSfuxEUXdeLllwvZddfo7QVX38SPL77wsmhROs88k0ZpqYcDDwxw553FnHpqGRkZDrBf5Gdbulg8Prhd6WGMuRG40e044s306Vl06lTJueeqypP40rGjw8MPF3HSSdmceWYqp556L3/4g3uzEEtLLZ5/vh2LFmWwfHkK7dqFOPXUMs49t5T+/Ru+4HNr79LQGhKxMnU96Un9ap+Yn3/u5a232nHNNRsj3z5F4ss++1QwefJH3HrroUybdiD33ZeLMdGtjL7+OplHH01nyZJ0Nm70sO++QW69dQOnn15Khw6N/1zF2+zLRK1MlfRiXF0n5n33HUfHjiHOP78VZ66ItLJg8Dks61kc5w4CgS+2LqNV/Que3++noKCAfv36tcg/yOXl8PLLaSxalM6yZamkpDiMGBGu6g49NNCsNULj7WLxeKxMW4KSXoyrfWI+++wPvP56Gn/608YW3X5EJNrCF6/nEggciOPcwpo1T9f4gnfzzTdz4403NqoS2V533Q8/hJcGW7w4naKiJHr2rODPf96AbZfRpUvztz+qEk8Xi8dbZdpSlPRiXO0Tc+XKUbRvH+KCC1TlSXzz+XwYk8fSpR/z/PPrWbz4RCoq9iMU+hyAl19+uVGVSO1ekccfNxQVHc6iRem88047kpIcTjhhC2PGlHLEEeV4Eny5/XirTFuKkl6Mq35idut2PJMmdefKKzc1acxBJNZUVUajR5dz3HFpFBc/g8czCK93IyNGjGDZsmVAwyqR33tFuhEKjePcc4eyeXMHsrKKOeusH5gypSNduzauqmvrEz3iqTJtKUp6caDqxBw/vhOZmSEuvHCz2yGJtKhu3UIsXLiJP/5xd7p1+4Bp075i0KBD6N27d4PG9EIhSE09GccZDIzAcSx23fV7vvvuQkpKXuDppz3k5ubRtWvjlx1LtIkebV2CF/jx45tvknnxxXaMHVtCp06q8qTtGTAgyN13b2TVqr154YWjgfAXvquvvnq7yebNN7/grLMKOOSQDtx8cw5ZWcczaNC7zJv3FqedNp9Q6FlCoUCT9pvTfnVtkyq9GFRXl8rMmZmkpTmMG6exPGm7Ro0qY8UKL/ffn0nv3sE6VxtyHMjPT2HGjHLee+8oIAWP522mTs1kwoTupKTsC4DfX9qsiRqJOtGjrVPSizF1dal07HgYzz2Xxvjxm+ncueVmmonEouuu28g33yTz5z93YJ99Khg5Mnx/cbHFkiXpLFqUzrffemnXrhTLmoPj3I9lfUtS0lRSUiZufZ3mTtRI1IkebZ2SXoyp69qZb789lpQUh0suUZUnbV9SEtx333pGjszm4os7MW9eCGM68vzzaWzZYjFgQIDp09fTvfv7nHvutVu/INZViTV3okYiTvRo65T0YkztLpWePY/hrrvSuPDCErKz667y2voMM0k87ds7PPRQESedtBOjRnnJyEjijDNKOeecEg44oCJy1ABVYtJoSnoxpnaXyhNP5OD1wvjxdc/Y1Awzaav23ruSJ574Hz/91Ikjj1xHZua2E7hUiUljKenFoKoP8qpVSSxZks5555Wwyy51V3mJupSQJIb+/YMcc0yIwkLNWJaWoaQXw2bPzsTjgUsv3f51eZphJiLScEp6MWr16iSMSeess0rp1m37MzY1w0wSicavpbmU9GLUnDmZAFx2Wf2rr2hcQxKBxq+lJWhFlhj0yy8eHn88Hdsujequ0iKxrDVWSPH7/cyaNQu/398CEUo8UKUXg+bOzaSyEi6/XGtsilRp6fFrVY6JSUkvxqxd6+GxxzIYNaqM3XdXlSdSpaXHrzXzOTEp6cWY++/PJBCAiRM3uR2KSMxpyfFrzXxOTEp6MaSw0MPChemcdloZe+6pKk+kNWnmc2JS0osh8+ZlsGWLxRVXqMoTiQbNfE48mr0ZI4qKLB56KINTTimjVy9VeSIircH1Ss+27Y7AfKAv4AAXGGMSbrfG+fMzKSuzuOIKzdgUEWktsVDpzQBeNcb0BvoDK1yOJ+qKiy0WLMjgxBO3sN9+FfU/QUREmsTVSs+27fbAUOB8AGNMAAi4GZMbFizIYNMmD5MmaSxPRKQ1ud29uRewDnjItu3+wCfAJGNMwuyWunGjxfz5mQwfXkafPqryRERak9tJLxkYAEw0xiyzbXsGcA1wQ/WDbNseB4wDMMaQnZ3d+DdKTm7S81rb/PkeNmzwcOONsRlflVhtv3igtmsetV/zqP1qcjvp/Qz8bIxZFrm9hHDSq8EYMw+YF7npFBYWNvqNsrOzacrzWtPmzRbTpu3CscduYbfdioix8GqIxfaLF2q75lH7NU+itF/37t0bdJyrE1mMMb8CP9m2vV/krmOAr1wMKaoeeSSD4mIPkydrLE9EJBrcrvQAJgKP2badAnwHjHU5nqgoLbW4//4MjjpqCwcfHHQ7HBGRhOB60jPGLAcSbkmEhQvTKSpKYvLkIrdDERFJGLFwnV7CKSuzuP/+TI44ohyfT1WeiEi0KOm54I47Clm3Lonhwz92OxQRkYSipBdlH3zwKfPnZwNvc8stx2nHZhGRKFLSi7J58yqBbsAtBINB8vMTbplRERHXuD6RJZGUl8Onnx6PZb2PZS3VxpUiIlGmpBdFxqTzv/+lcfPNHsrKpmrjShGRKFPSi5JAAGbNymTAgAAXXrgHljXR7ZBERBKOxvSi5Kmn0lm9Opkrr9yEZbkdjYhIYlLSi4KKinCV179/gKOOKnc7HBGRhKWkFwXPPJPGjz8mM3myqjwRETcp6bUCv9/PrFmz8Pv9VFbCjBlZHHBAkOOOU5UnIuImTWRpYX6/n9zcXILBIF6vlwkT3uP777vzj38UqcoTEXGZkl4Ly8/PJxgMUllZieN4eOSRHvTuHWTYsC1uhyYikvDUvdnCcnJy8Hq9JCUlkZRk87//7cKkSZvw7KClq3eHiohI61Gl18J8Ph95eXl88MGHLF58DSkpQU48cftVXu3u0Ly8PF2wLiLSSlTptQKfz0evXn/ixx/bM2nSZj77bPuVXPXuUK3FKSLSulTptQLHgWnTsthrrwq6d393h5VcVXcooLU4RURamZJeK3jjjVS++srL9Onr+eij3ys5CFd21ZNeVXdofn6+1uIUEWllSnotwO/3b01ahxziY9q0LPbYo4LTTitj+fL6Kzmfz6dkJyISBUp6zVR7IsrUqf+ioKA7d99dTHKyKjkRkViipNdMNa/Lgwcf7E6PHhWMGlW69RhVciIisUGzN5up5nV5w1izZncmTtxMpEdTRERiiCq9Zvr9urx8nntuKhs2VHLGGaX1P1FERKJOlV4L8Pl8DBgwhZUrO3P55ZtITXU7IhERqUtMVHq2bScBfmC1MWak2/E0xfTpWeyySyVnnqkqT0QkVsVKpTcJWOF2EE314Ycp5OenMmHCZtq1czsaERHZHteTnm3bPYATgflux9JU06dnsdNOlZx9tqo8EZFYFgvdm9OBq4Gs7R1g2/Y4YByAMYbs7OxGv0lycnKTnlefDz+0ePddL3/7WwW77dalxV8/VrRW+yUCtV3zqP2aR+1Xk6tJz7btkcBaY8wntm0fub3jjDHzgHmRm05hYWGj3ys7O5umPK8+N93Umc6dPYwatY7CQqfFXz9WtFb7JQK1XfOo/ZonUdqve/fuDTrO7e7Nw4GTbdv+AVgMHG3b9qPuhtRwn33m5V//asf48SWkp7fdhCci0la4WukZY64FrgWIVHp/MsaMcTOmxpg+PYuOHUOcd16J26GIiEgDuF3pxa0vvvDy5pvtGDduM5mZqvJEROJBLExkAcAY8zbwtsthNNj06Zl06BBi7FhVeSIi8UKVXhP8+9/JvPpqGhddtJn27VXliYjECyW9JpgxI4usrBAXXKAqT0QknijpNdLXXyfz8svtGDu2hI4dVeWJiMQTJb1Gmjkzk/R0h4sv3ux2KCIi0khKenXw+/3MmjULv99f4/5vv03muefSGDu2hM6dVeWJiMSbmJm9GSv8fj+5ubkEg0G8Xi95eXlbdz2fOTOTdu0cxo3TWJ6ISDxSpVdLfn4+wWCQyspKgsEg+fn5AHz/fRLPPJPGueeW0qVLyOUoRUSkKZT0asnJycHr9ZKUlITX6yUnJweAWbOySEmB8eM3b7f7U0REYpu6N2vx+Xzk5eWRn59PTk4OPp+PVauSWLIkjfPPL2HVqo+22/0pIiKxTUmvDj6fr0Yimz07k+RkmDBhM08++Xv3J4S7Q5X0RETig5JePVavTsKYdM46q5SuXUNbuz+BGt2fIiIS+5T06jF7diaO45CRMQu/v0+d3Z8iIhIflPR2YM0aD48/3g7HeZAHHvgzCxb8PoanZCciEn80e3MH5s7NJBQCx7ljm0sYREQk/ijpbcdvv3l47LEMjj76Z1JS1mxzCYOIiMQfdW9ux9y5mVRUwM03t6OwUGN4IiJtgZJeHQoLPSxalM4f/1hGz56V9OypMTwRkbZA3Zt1eOCBDAIBi4kTN7kdioiItCAlvVqKijw8/HAGp55axt57V7odjoiItCAlvVrmzcugrMziiiu0X56ISFujpFfN+vUWDz2UwciRW9hnnwq3wxERkRampFfNgw9msnmzh0mTNJYnItIWKelFbNhg8eCDGYwYUcb++6vKExFpi1y9ZMG27d2AhUBXIATMM8bMcCOWBQsy2LhRVZ6ISFvmdqVXAUwxxuwPDAIus227T7SD2LTJYv78TI4/voy+fVXliYi0Va4mPWPML8aYTyP/vwlYAewa7TgeeiiD4mIPkydrxqaISFvmdqW3lW3bPYGDgWXRfN+SEot58zI4+ugt9O8fjOZbi4hIlMXEMmS2bWcCTwGTjTEb63h8HDAOwBhDdnZ2o98jOTm5zuc9/LCH9euTuOmmUJNeN1Fsr/2kfmq75lH7NY/arybLcRxXA7Bt2wu8CLxmjLm3AU9x1qxZ0+j3yc7OprCwsMZ9ZWUWhx22M337Bnn88aJGv2Yiqav9pGHUds2j9mueRGm/7t27A1j1Hef27E0LeBBY0cCE1yL8fj/5+fmsXXs2//tfN668cn203lpERFzkdvfm4cA5wBe2bS+P3HedMebl1npDv99Pbm4ugUASodAVHHjgOgYO1FieiEgicDXpGWPeowHlaEvKz88nGAwSCl0EdKNv33nAyGiGICIiLomZ2ZvRkpOTQ3JyJnANlvUeubld3Q5JRESiJOGSns/nY+zYpcCu3HSTh4EDtTmsiEiiSLikFwjA88/34ZBDAlx4YU+3wxERkShyeyJL1BUXe+jdu4ILLijBiowmVs3mzMnJwedT5Sci0lYlXNLbeecQixb9fk1e1WzOYDCI1+slLy9PiU9EpI1KuO7N2qpmc1ZWVhIMBsnPz3c7JBERaSUJV+nVlpOTg9frBcDr9ZKTk+NyRCIi0loSPun5fD7y8vI0picikgASPulBOPEp2YmItH0JP6YnIiKJQ0lPREQShpKeiIgkDCU9ERFJGEp6IiKSMJT0REQkYSjpiYhIwlDSExGRhKGkJyIiCcNyHMftGBor7gIWEZGosOo7IB4rPaspP7Ztf9LU5+pH7ae2U/vF60+CtV+94jHpiYiINImSnoiIJIxESnrz3A4gzqn9mk5t1zxqv+ZR+1UTjxNZREREmiSRKj0REUlwbX4TWdu2hwEzgCRgvjHmDpdDimm2be8GLAS6AiFgnjFmhm3bnYE8oCfwA2AbY9a7FWcss207CfADq40xI23b3hNYDHQGPgXOMcYE3Iwxltm23RGYD/QlfInSBcDX6Pyrl23bVwIXEW63L4CxQDd0/m3Vpiu9yD8+9wHDgT7AaNu2+7gbVcyrAKYYY/YHBgGXRdrsGuCfxph9gH9GbkvdJgErqt3+OzAt0nbrgQtdiSp+zABeNcb0BvoTbkudf/WwbXtX4ArAZ4zpS/iL/pno/KuhTSc94FDgW2PMd5FvNouBU1yOKaYZY34xxnwa+f9NhP/B2ZVwuz0SOewR4FR3Ioxttm33AE4kXKlg27YFHA0siRyittsB27bbA0OBBwGMMQFjTDE6/xoqGUizbTsZSAd+QedfDW096e0K/FTt9s+R+6QBbNvuCRwMLAN2Mcb8AuHECOzsYmixbDpwNeGuYYAuQLExpiJyW+fgju0FrAMesm37M9u259u2nYHOv3oZY1YDdwOrCCe7DcAn6Pyroa0nvbqu0Nd01QawbTsTeAqYbIzZ6HY88cC27ZHAWmPMJ9Xu1jnYOMnAAGCuMeZgoAR1ZTaIbdudCFfEewLdgQzCQzu1JfT519aT3s/AbtVu9wDWuBRL3LBt20s44T1mjHk6cvdvtm13izzeDVjrVnwx7HDgZNu2fyDclX404cqvY6S7CXQO1udn4GdjzLLI7SWEk6DOv/odC3xvjFlnjAkCTwOD0flXQ1tPeh8D+9i2vadt2ymEB3WfdzmmmBYZg3oQWGGMubfaQ88D50X+/zzguWjHFuuMMdcaY3oYY3oSPtfeMsacDfwLGBU5TG23A8aYX4GfbNveL3LXMcBX6PxriFXAINu20yOf46q20/lXTZu/ON227RGEv20nAQuMMbe7HFJMs217CPAu4enOVeNS1xEe1zPA7oQ/XGcYY4pcCTIO2LZ9JPCnyCULe/H7lPHPgDHGmHI344tltm0fRHgiUArwHeFp9x50/tXLtu2bgVzCs7A/I3z5wq7o/NuqzSc9ERGRKm29e1NERGQrJT0REUkYSnoiIpIwlPRERCRhKOmJiEjCUNITEZGEoaQnIiIJQ0lPREQSRpvfRFakLbFte2/Cy+sda4z51Lbt7kABMMoY87arwYnEAa3IIhJnbNu+GLgKOAR4BvjCGPMnd6MSiQ/q3hSJM8aYfwD/IbweajfgencjEokfSnoi8ekfQF9gViIvHizSWOreFIkzkQ1+Pye8Zcxw4EDtOCDSMKr0ROLPDOATY8xFwEvA/S7HIxI3lPRE4oht26cAw4DxkbuuAgbYtn22e1GJxA91b5yNQrQAAABCSURBVIqISMJQpSciIglDSU9ERBKGkp6IiCQMJT0REUkYSnoiIpIwlPRERCRhKOmJiEjCUNITEZGEoaQnIiIJ4/8BBB9L3qOsk28AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize = (7,5))\n",
    "ax.plot(x, y, marker=\".\", color=\"k\", linewidth=0)\n",
    "ax.plot(x, y_tru, color=\"b\")\n",
    "ax.set_xlabel('x'); ax.set_ylabel('y')\n",
    "ax.set_title('simulation data')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. サンプリング用関数群\n",
    "- ギブスサンプリングを行うための関数群を定義する\n",
    "  - 本来クラスにまとめた方が簡潔になるが、試行用なので略。\n",
    "  \n",
    "まず潜在変数zのサンプリングから"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1. 潜在変数z\n",
    "- 既存のクラスタ数を$K^+$として,  \n",
    "$k \\in \\{1,...,K^+\\}$のとき\n",
    "$$p(z_k | y_t,x_t,\\theta^{(k)},z_{1:T}^{\\t}) \\sim {\\cal{N}} (y_t | \\theta^{(k)\\mathrm{T}} , \\sigma^{(k)2}) \\times \\frac{n_k^\\t}{T-1+\\alpha}$$\n",
    "$k \\notin \\{1,...,K^+\\}$のとき\n",
    "$$p(z_k | y_t,x_t,\\theta^{(k)},z_{1:T}^{\\t}) \\sim {\\cal{N}} (y_t | \\theta^{(new)\\mathrm{T}} , \\sigma^{(new)2}) \\times \\frac{\\alpha}{T-1+\\alpha}$$\n",
    "$$(\\theta^{(new)}, \\sigma^{(new)}) \\sim {\\cal{N}} (\\mu, \\Sigma) IG(\\frac{n_0}{2}, \\frac{\\tau}{2})$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sampler_z(t, yt, xt, theta, sigma_y, z, alpha, T, mu, sigma, n0, tau):\n",
    "    k_plus, n_t = np.unique(z, return_counts=True)\n",
    "    prob_z = np.empty(k_plus.shape[0]+1)#P(z_t = k) not notrmalize\n",
    "    for k in k_plus:\n",
    "        ## culculate n_k^¥t\n",
    "        if(z[t]==k):\n",
    "            n_kt = n_t[k_plus==k] - 1\n",
    "        else:\n",
    "            n_kt = n_t[k_plus==k]\n",
    "        ## culculate probability s.t. z_t = k\n",
    "        prob_z[k] = st.norm.pdf(x=yt, loc = np.dot(theta[k_plus==k], np.r_[xt,1]),\\\n",
    "                                    scale=sigma_y[k_plus==k])*(n_kt/(T-1+alpha))\n",
    "\n",
    "    ## create new cluster\n",
    "    theta_new = st.multivariate_normal.rvs(mean=mu, cov=sigma, size=1)\n",
    "    sigma_new = np.sqrt(st.invgamma.rvs(a=n0/2, scale=tau/2, size=1))\n",
    "    prob_z[-1] = st.norm.pdf(x=yt, loc = np.dot(theta_new, np.r_[xt,1]), scale=sigma_new)*(alpha/(T-1+alpha))\n",
    "    \n",
    "    ## sampling\n",
    "    prob_z /= np.sum(prob_z)\n",
    "    z_sample = np.random.choice(np.r_[k_plus, k_plus[-1]+1], size=1, p=prob_z)[0]#random.choice?\n",
    "    return z_sample, theta_new, sigma_new\n",
    "    ### return new z[t]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- クラスタが追加された場合に$\\theta^{new}$と$\\sigma^{new}$も返す必要がある\n",
    "  - 当座は毎回$theta^{new}$を返す仕様とし, mainの側でサンプリング値$z_k$が新クラスタならば$theta$に$theta^{new}$を追加する処理を加えておく。\n",
    "  - クラス化してしまえば内部で処理できるが、今回は手をつけていない。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2. パラメータ$\\theta$\n",
    "- これはクラスタ毎に定まる。$\\cal{T}_k = \\{t|z_t=k\\}$に対して、\n",
    "$$p(\\theta^{(k)} | y,x,z,\\sigma^{(k)}, \\mu, \\Sigma) \\sim {\\cal{N}}(\\mu_k, \\Sigma_k).$$\n",
    "ただし\n",
    "$$\\Sigma_k^{-1} = \\sum_{t \\in \\cal{T}_k}{\\frac{x_t x_t^{\\mathrm{T}}}{\\sigma^{(k)2}}} + \\Sigma^{-1}$$\n",
    "$$\\mu_k = \\Sigma_k(\\sum_{t \\in \\cal{T}_k}{\\frac{y_t x_t}{\\sigma^{(k)2}}} + \\Sigma^{-1}\\mu)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sampler_theta(k, y, x, z, sigma_y, mu, sigma_inv):\n",
    "    #print(mu, sigma_inv)\n",
    "    sigma_yk = sigma_y[k]\n",
    "    t_k = np.where(z == k)[0]\n",
    "    x = np.c_[x.copy(), np.ones(x.shape[0])]#X_t = [X_t, 1]\n",
    "    sigma_k_inv = sigma_inv.copy()\n",
    "    mu_k_tmp = np.dot(sigma_inv.copy(), mu) #Sigma^-1 * mu\n",
    "    for tt in t_k: # Sigma_(t in T_k)\n",
    "        x_tt = x[tt][:,np.newaxis]\n",
    "        sigma_k_inv += np.dot(x_tt,x_tt.T) / (sigma_yk**2)\n",
    "        mu_k_tmp += (y[tt]*x[tt]) / (sigma_yk**2)\n",
    "    sigma_k = np.linalg.inv(sigma_k_inv)\n",
    "    mu_k = np.dot(sigma_k, mu_k_tmp)\n",
    "    del mu_k_tmp\n",
    "    return st.multivariate_normal.rvs(mean=mu_k, cov=sigma_k, size=1)\n",
    "    ### return new theta[k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3. 観測時ノイズの分散$\\sigma^2$\n",
    "- これもクラスタ毎に定まる。\n",
    "$$p(\\sigma^{(k)2} | y,x,z,\\theta^{(k)},n_0, \\tau) = IG(\\sigma^{(k)2} | \\frac{n_0+n_k}{2}, \\frac{\\tau_k}{2})$$\n",
    "ただし\n",
    "$$\\tau_k = \\tau + \\sum_{t \\in \\cal{T}_k}{||y_t - \\theta^{(k)\\mathrm{T}}x_t||^2}$$\n",
    "- 下記の関数ではsqrtを取り$\\sigma^{(k)}$をサンプリングしている。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sampler_sigma_y(k, y, x, z, theta, n0, tau):\n",
    "    t_k = np.where(z == k)[0]\n",
    "    n_k = t_k.shape[0]\n",
    "    tau_k = tau\n",
    "    for tt in t_k:\n",
    "        resid = y[tt] - np.dot(theta[k],np.r_[x[tt],1])\n",
    "        tau_k += np.dot(resid, resid)\n",
    "    return np.sqrt(st.invgamma.rvs(a=(n0+n_k)/2, scale=tau_k/2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4. $\\mu$ ($\\theta^{(new)}$の平均)\n",
    "- サンプリング分布は\n",
    "$$p(\\mu | \\theta^{(1:K^+)}, \\Sigma, \\mu_0, V_0) = {\\cal{N}}(\\mu | \\mu_+, V_+)$$  \n",
    "  ここで$\\mu_0, V_0$は\n",
    "$$V_0^{-1} = K^+\\Sigma^{-1}+V_0^{-1}$$\n",
    "$$\\mu_+ = V_+(\\Sigma^{-1}\\sum_{k=1}^{K^+}{\\theta^{(k)}}+V_0^{-1}\\mu_0)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sampler_mu(theta, sigma_inv, mu0, v0_inv):\n",
    "    vp = np.linalg.inv(theta.shape[0] * sigma_inv + v0_inv )\n",
    "    mup = np.dot(sigma_inv, np.sum(theta, axis=0)) + np.dot(v0_inv, mu0)\n",
    "    mup = np.dot(vp, mup)\n",
    "    return st.multivariate_normal.rvs(mean=mup, cov=vp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5. $\\Sigma$ ($\\theta^{(new)}$の共分散行列)\n",
    "- サンプリング分布は\n",
    "$$p(\\Sigma^{-1} | \\theta^{(1:K^+)}, \\mu, \\nu_0, \\Sigma_0) = W(\\Sigma^{-1} | \\nu_+, \\Sigma_+)$$\n",
    "ただし\n",
    "$$\\nu_+ = \\nu_0 + K^+$$\n",
    "$$\\Sigma_+^{-1} = \\Sigma_0^{-1} + \\sum_{k=1}^{K^+}{(\\theta^{(k)}-\\mu)(\\theta^{(k)}-\\mu)}^\\mathrm{T}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sampler_Sigma_inv(theta, mu, nu0, sigma0_inv):\n",
    "    nup = nu0 + theta.shape[0]\n",
    "    sigmap_inv = sigma0_inv\n",
    "    for ii in range(theta.shape[0]):\n",
    "        tmp = (theta[ii] - mu)[:,np.newaxis]\n",
    "        sigmap_inv += np.dot(tmp, tmp.T)\n",
    "    return st.wishart.rvs(df=nup, scale=np.linalg.inv(sigmap_inv))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.6. $\\tau$ ($\\sigma_y$の生成過程のパラメータ)\n",
    "- サンプリング分布は\n",
    "$$p(\\tau | \\sigma^{(1:K^+)}, n_0, m_0, \\tau_0) = Ga(\\tau | \\frac{m_+}{2}, \\frac{\\tau_+}{2})$$\n",
    "ただし\n",
    "$$m_+ = m_0 + n_0K^+$$\n",
    "$$\\tau_0 = \\tau_0 + \\sum_{k=1}^{K^+}{\\frac{1}{\\sigma^{(k)2}}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sampler_tau(sigma_y, n0, m0, tau0):\n",
    "    mp = m0 + n0*sigma_y.shape[0]\n",
    "    taup = tau0 + np.sum(1/sigma_y**2)\n",
    "    return st.gamma.rvs(a=mp/2, scale=taup/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. サンプリング"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. 初期値設定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ハイパーパラメータ群"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 1\n",
    "n0 = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "事前分布のパラメータ群"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu0 = np.array([0, 0])\n",
    "v0_tmp = np.random.uniform(-10, 10, (2, 2))#np.array([[0.5, 0], [0,4]])\n",
    "v0 = np.dot(v0_tmp, v0_tmp.T)\n",
    "v0_inv = np.linalg.inv(v0)\n",
    "\n",
    "nu0 = 2\n",
    "sigma0_tmp = np.random.uniform(-10, 10, (2, 2))#np.array([[0.5, 0], [0,10]])#np.random.uniform(-1, 1, (2, 2))\n",
    "sigma0 = np.dot(sigma0_tmp, sigma0_tmp.T) # positive definite\n",
    "sigma0_inv = np.linalg.inv(sigma0)\n",
    "\n",
    "m0 = 0.5\n",
    "tau0 = 2 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "パラメータの初期値"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
      " 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]\n",
      "mu :  [-9.56246008  7.33209689]\n",
      "Sigma :  [[120.50695055 -49.34876924]\n",
      " [-49.34876924  23.14239448]]\n",
      "tau : 0.24918271440471546\n"
     ]
    }
   ],
   "source": [
    "T = x.shape[0]\n",
    "\n",
    "theta = np.array([[0.5,1], [-0.3, 25]])#, [0.1, 1]])\n",
    "sigma_y = np.sqrt(np.array([0.2, 0.2]))#, 0.2]))\n",
    "z = np.repeat(np.array([0,1]), 45)#np.repeat(np.array([0,1,2]), 30) #np.zeros(T, dtype=\"int\")\n",
    "print(z)\n",
    "mu = st.multivariate_normal.rvs(mean=mu0, cov=v0)\n",
    "sigma = st.wishart.rvs(df=nu0, scale=sigma0)\n",
    "sigma_inv = np.linalg.inv(sigma)\n",
    "tau = st.gamma.rvs(a=m0/2, scale=tau0/2)\n",
    "print(\"mu : \", mu)\n",
    "print(\"Sigma : \", sigma)\n",
    "print(\"tau :\", tau)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. ギブスサンプリング\n",
    "- 本では少なくとも12000ステップのサンプリングを行っていたが、ここでは2000ステップ（うちバーンイン1000ステップ）と設定した。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4min 27s, sys: 1.26 s, total: 4min 29s\n",
      "Wall time: 4min 31s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "n_step = 2000\n",
    "burnin = 1000\n",
    "z_sample = np.zeros((n_step, T))\n",
    "theta_sample = np.zeros((n_step, T, theta.shape[1]))\n",
    "\n",
    "n_cluster = theta.shape[0]\n",
    "z_unique = np.unique(z)\n",
    "for step in range(n_step):\n",
    "    #print(step)\n",
    "    for tt in range(T):\n",
    "        z_new, theta_new, sigma_new = sampler_z(tt, y[tt], x[tt], theta, sigma_y, z, alpha, T, mu, sigma, n0, tau)\n",
    "        z[tt] = z_new\n",
    "        if(z_unique is not np.unique(z)):\n",
    "            z_unique = np.unique(z)\n",
    "            # クラスタ数が増えた場合の処理\n",
    "            if(n_cluster < z_unique.shape[0]):\n",
    "                n_cluster += 1\n",
    "                theta = np.r_[theta, [theta_new]]\n",
    "                sigma_y = np.r_[sigma_y, sigma_new]\n",
    "                #print(z_new, theta, sigma_y)\n",
    "            # クラスタ数が減った場合の処理\n",
    "            elif(n_cluster > z_unique.shape[0]):\n",
    "                n_cluster -= 1\n",
    "                theta = theta[z_unique]\n",
    "                sigma_y = sigma_y[z_unique]\n",
    "            #Zを0,1,2,...に置き換える処理\n",
    "            for ii, z_val in zip(range(z_unique.shape[0]), z_unique):\n",
    "                z[z==z_val] = ii\n",
    "\n",
    "    for kk in range(n_cluster):\n",
    "        theta[kk] = sampler_theta(kk, y, x, z, sigma_y, mu, sigma_inv.copy())\n",
    "        sigma_y[kk] = sampler_sigma_y(kk, y, x, z, theta, n0, tau)\n",
    "    mu = sampler_mu(theta, sigma_inv, mu0, v0_inv)\n",
    "    sigma_inv = sampler_Sigma_inv(theta, mu, nu0, sigma0_inv.copy())\n",
    "    tau = sampler_tau(sigma_y, n0, m0, tau0)\n",
    "\n",
    "    z_sample[step] = z\n",
    "    theta_sample[step] = np.r_[theta, np.repeat([np.repeat(np.nan,theta.shape[1])], T-theta.shape[0], axis=0)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 旧型のMacBookPro(Mid2012)でも、4分ほどで実行できている。\n",
    "- sampler_thetaに渡すsigma_invは、copy()を取らないと内部で値が変化してしまう。\n",
    "  - 内部でも代入はしていませんが、他のベクトルと内積を取った時に値が変わった？\n",
    "  - これで丸一日を溶かしました。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.unique(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. 結果\n",
    "### 4.1. クラスタ数の分布"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "サンプリングにより得られたクラスタ数の配分は以下。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([2, 3, 4]), array([ 12, 957,  31]))\n"
     ]
    }
   ],
   "source": [
    "print(np.unique(np.max(z_sample[burnin:].astype('int')+1, axis=1), return_counts=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "これを図示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFopJREFUeJzt3XuQXnV9x/F3ZAGlqAEWgU2iYEm9MV4ZpDL1hqOClOBUvuIFA0YztXjFKmgdsV46WDsi3o1gm9QLfEUsUamI4LUVFFKsVVqHAiWBmLAkoIjKBLZ/nN/qM5tN9mGf3bO/s7xfM8/sc37n95zzzY/z7Idz2XMWjI2NIUlSbR4w1wVIkjQZA0qSVCUDSpJUJQNKklQlA0qSVCUDSpJUJQNKklSloak6RMRngGOAzZl5SGnbGzgfOBC4EYjM3BoRC4CzgaOBu4CTMnNd+cxy4B1lse/NzNUz+0+RJM0n/exB/RPw/AltpwOXZeZS4LIyDXAUsLS8VgKfgN8H2hnAU4HDgDMiYq9Bi5ckzV9TBlRmfhfYMqF5GTC+B7QaOK6nfU1mjmXmFcDCiDgAeB5waWZuycytwKVsH3qTGfPly5cvX/PyNaUpD/HtwH6ZuREgMzdGxMNK+yJgfU+/DaVtR+3biYiVNHtfZCZ33333NEtsDA0NsW3btoGW0RZrnXldqRO6U2tX6oTu1NqVOmFmat1tt936W9dAa9negknaxnbSvp3MXAWsGu8zOjo6UEHDw8MMuoy2WOvM60qd0J1au1IndKfWrtQJM1PryMhIX/2mexXfpnLojvJzc2nfACzp6bcYuGUn7ZIkTWq6AbUWWF7eLwcu6ml/RUQsiIjDgTvKocBLgOdGxF7l4ojnljZJkibVz2XmXwCeCQxHxAaaq/HOBDIiVgA3AceX7hfTXGJ+Hc1l5icDZOaWiHgP8KPS792ZOfHCC0mSfm/KgMrMl+xg1pGT9B0DTtnBcj4DfOY+VSdJut/yThKSpCoZUJKkKhlQkqQqGVCSpCrN9B/qSpqmTS982lyXAMAun1471yVIgHtQkqRKGVCSpCoZUJKkKhlQkqQqGVCSpCoZUJKkKhlQkqQqGVCSpCoZUJKkKhlQkqQqGVCSpCoZUJKkKhlQkqQqGVCSpCoZUJKkKhlQkqQqGVCSpCoZUJKkKhlQkqQqGVCSpCoZUJKkKhlQkqQqGVCSpCoZUJKkKhlQkqQqGVCSpCoZUJKkKhlQkqQqGVCSpCoZUJKkKhlQkqQqDQ3y4Yh4E/AqYAz4CXAycABwHrA3sA44MTPvjojdgTXAU4DbgBdn5o2DrF+SNH9New8qIhYBrwcOzcxDgF2AE4D3A2dl5lJgK7CifGQFsDUzDwbOKv0kSZrUoIf4hoAHRcQQsAewEXg2cEGZvxo4rrxfVqYp84+MiAUDrl+SNE9NO6Ay82bgH4CbaILpDuBq4PbM3Fa6bQAWlfeLgPXls9tK/32mu35J0vw27XNQEbEXzV7RQcDtwBeBoybpOlZ+Tra3NDaxISJWAisBMpPh4eHplgjA0NDQwMtoi7XOvK7UCbBprgsophqvLo1pV2rtSp3Qbq2DXCTxHOCGzLwVICIuBJ4GLIyIobKXtBi4pfTfACwBNpRDgg8FtkxcaGauAlaVybHR0dEBSmy+bIMuoy3WOvO6UmdNphqvLo1pV2rtSp0wM7WOjIz01W+QgLoJODwi9gB+AxwJXAV8C3gRzZV8y4GLSv+1ZfoHZf7lmbndHpQkSTDYOagraS52WEdzifkDaPZ8TgNOjYjraM4xnVs+ci6wT2k/FTh9gLolSfPcQH8HlZlnAGdMaL4eOGySvr8Fjh9kfZKk+w/vJCFJqpIBJUmqkgElSaqSASVJqpIBJUmqkgElSaqSASVJqpIBJUmqkgElSaqSASVJqpIBJUmqkgElSaqSASVJqpIBJUmqkgElSaqSASVJqpIBJUmqkgElSaqSASVJqpIBJUmqkgElSaqSASVJqpIBJUmqkgElSaqSASVJqpIBJUmqkgElSaqSASVJqpIBJUmqkgElSaqSASVJqpIBJUmqkgElSaqSASVJqpIBJUmqkgElSaqSASVJqtLQIB+OiIXAOcAhwBjwSuB/gPOBA4EbgcjMrRGxADgbOBq4CzgpM9cNsn5J0vw16B7U2cDXM/PRwBOAa4HTgcsycylwWZkGOApYWl4rgU8MuG5J0jw27YCKiIcATwfOBcjMuzPzdmAZsLp0Ww0cV94vA9Zk5lhmXgEsjIgDpl25JGleG+QQ3yOBW4F/jIgnAFcDbwD2y8yNAJm5MSIeVvovAtb3fH5DadvYu9CIWEmzh0VmMjw8PECJMDQ0NPAy2mKtM68rdQJsmusCiqnGq0tj2pVau1IntFvrIAE1BDwZeF1mXhkRZ/OHw3mTWTBJ29jEhsxcBawanz86OjpAic2XbdBltMVaZ15X6qzJVOPVpTHtSq1dqRNmptaRkZG++g1yDmoDsCEzryzTF9AE1qbxQ3fl5+ae/kt6Pr8YuGWA9UuS5rFpB1Rm/gJYHxGPKk1HAj8D1gLLS9ty4KLyfi3wiohYEBGHA3eMHwqUJGmigS4zB14HfC4idgOuB06mCb2MiBXATcDxpe/FNJeYX0dzmfnJA65bkjSPDRRQmXkNcOgks46cpO8YcMog65Mk3X94JwlJUpUMKElSlQwoSVKVDChJUpUMKElSlQwoSVKVDChJUpUMKElSlQwoSVKVDChJUpUMKElSlQwoSVKVDChJUpUMKElSlQwoSVKVDChJUpUMKElSlQwoSVKVDChJUpUMKElSlQwoSVKVDChJUpUMKElSlQwoSVKVDChJUpUMKElSlQwoSVKVDChJUpUMKElSlQwoSVKVDChJUpUMKElSlQwoSVKVDChJUpUMKElSlQwoSVKVhgZdQETsAlwF3JyZx0TEQcB5wN7AOuDEzLw7InYH1gBPAW4DXpyZNw66fknS/DQTe1BvAK7tmX4/cFZmLgW2AitK+wpga2YeDJxV+kmSNKmBAioiFgMvAM4p0wuAZwMXlC6rgePK+2VlmjL/yNJfkqTtDLoH9SHgrcC9ZXof4PbM3FamNwCLyvtFwHqAMv+O0l+SpO1M+xxURBwDbM7MqyPimaV5sj2isT7m9S53JbASIDMZHh6ebokADA0NDbyMtljrzOtKnQCb5rqAYqrx6tKYdqXWrtQJ7dY6yEUSRwDHRsTRwAOBh9DsUS2MiKGyl7QYuKX03wAsATZExBDwUGDLxIVm5ipgVZkcGx0dHaDE5ss26DLaYq0zryt11mSq8erSmHal1q7UCTNT68jISF/9pn2ILzPflpmLM/NA4ATg8sx8GfAt4EWl23LgovJ+bZmmzL88M7fbg5IkCWbn76BOA06NiOtozjGdW9rPBfYp7acCp8/CuiVJ88TAfwcFkJnfBr5d3l8PHDZJn98Cx8/E+iRJ8593kpAkVcmAkiRVyYCSJFXJgJIkVcmAkiRVyYCSJFXJgJIkVcmAkiRVyYCSJFXJgJIkVcmAkiRVyYCSJFXJgJIkVcmAkiRVyYCSJFXJgJIkVcmAkiRVyYCSJFXJgJIkVcmAkiRVyYCSJFXJgJIkVcmAkiRVyYCSJFXJgJIkVcmAkiRVyYCSJFXJgJIkVcmAkiRVyYCSJFXJgJIkVcmAkiRVyYCSJFXJgJIkVcmAkiRVyYCSJFXJgJIkVWlouh+MiCXAGmB/4F5gVWaeHRF7A+cDBwI3ApGZWyNiAXA2cDRwF3BSZq4brHxJ0nw1yB7UNuDNmfkY4HDglIh4LHA6cFlmLgUuK9MARwFLy2sl8IkB1i1JmuemHVCZuXF8DygzfwVcCywClgGrS7fVwHHl/TJgTWaOZeYVwMKIOGDalUuS5rVpH+LrFREHAk8CrgT2y8yN0IRYRDysdFsErO/52IbStnHCslbS7GGRmQwPDw9U29DQ0MDLaIu1zryu1Amwaa4LKKYary6NaVdq7Uqd0G6tAwdUROwJfAl4Y2b+MiJ21HXBJG1jExsycxWwanz+6OjoQPUNDw8z6DLaYq0zryt11mSq8erSmHal1q7UCTNT68jISF/9BrqKLyJ2pQmnz2XmhaV50/ihu/Jzc2nfACzp+fhi4JZB1i9Jmr8GuYpvAXAucG1mfrBn1lpgOXBm+XlRT/trI+I84KnAHeOHAiVJmmiQQ3xHACcCP4mIa0rb22mCKSNiBXATcHyZdzHNJebX0VxmfvIA65YkzXPTDqjM/D6Tn1cCOHKS/mPAKdNdnyTp/sU7SUiSqmRASZKqZEBJkqpkQEmSqmRASZKqZEBJkqpkQEmSqmRASZKqZEBJkqpkQEmSqmRASZKqZEBJkqpkQEmSqmRASZKqZEBJkqpkQEmSqmRASZKqZEBJkqpkQEmSqmRASZKqZEBJkqpkQEmSqmRASZKqZEBJkqpkQEmSqmRASZKqZEBJkqpkQEmSqmRASZKqZEBJkqpkQEmSqmRASZKqZEBJkqo0NNcFSJKmds+rj53rEhpf/vfWVuUelCSpSgaUJKlKBpQkqUqtn4OKiOcDZwO7AOdk5pmzub5NL3zabC6+b7t8eu1clyBJndLqHlRE7AJ8DDgKeCzwkoh4bJs1SJK6oe1DfIcB12Xm9Zl5N3AesKzlGiRJHdD2Ib5FwPqe6Q3AU3s7RMRKYCVAZjIyMjLYGr921WCfb9nA/94WdaXWrtTZpW21M2NKd2qdss6Kto+2xrTtPagFk7SN9U5k5qrMPDQzDy39B3pFxNUzsZw2XtZ6/62zS7V2pc4u1dqVOme41im1HVAbgCU904uBW1quQZLUAW0f4vsRsDQiDgJuBk4AXtpyDZKkDmh1DyoztwGvBS4Brm2a8qezvNpVs7z8mWStM68rdUJ3au1KndCdWrtSJ7RY64KxsbGpe0mS1DLvJCFJqpIBJUmqUmcftxERS4A1wP7AvcCqzDx7Qp8FNLdVOhq4CzgpM9eVecuBd5Su783M1XNc68uA08rkncBrMvPHZd6NwK+Ae4Bt5RL8uarzmcBFwA2l6cLMfHeZ19ptrPqs9S3Ay8rkEPAYYN/M3NLimD4Q+C6we6nhgsw8Y0Kf3cu/5SnAbcCLM/PGMu9twIpS5+sz85LZqPM+1Hoq8CpgG3Ar8MrM/L8y7x7gJ6XrTZk5a8+H6LPWk4AP0FyQBfDRzDynzGvl+99nnWcBzyqTewAPy8yFZV5rY9pTzy7AVcDNmXnMhHmtbqtd3oPaBrw5Mx8DHA6cMsltk44ClpbXSuATABGxN3AGzR8JHwacERF7zXGtNwDPyMzHA+9h+xORz8rMJ87WL9L7UCfA90otT+wJp7ZvYzVlrZn5gfE6gbcB38nMLT1d2hjT3wHPzswnAE8Enh8Rh0/oswLYmpkHA2cB7wco/54TgMcBzwc+XsZ5Lmv9D+DQsp1eAPx9z7zf9GwXs/2LtJ9aAc7vqWk8nNr8/k9ZZ2a+qWc7/QhwYc/sNsd03BtoLmKbTKvbamf3oDJzI7CxvP9VRFxLc6eKn/V0Wwasycwx4IqIWBgRBwDPBC4d/2UVEZfSDOoX5qrWzOx9CtgVNH8j1qo+x3RHfn8bK4CIGL+NVT+fbaPWlzBL/313pmx7d5bJXctr4pVJy4B3lfcXAB8te//LgPMy83fADRFxHc04/2Cuas3Mb/VMXgG8fDZqmUqf47ojz6Ol7/806nwJTXjOiYhYDLwAeB9w6iRdWt1WOxtQvSLiQOBJwJUTZk12a6VFO2mfdTuptdcK4F97pseAb0TEGPCpzJz1yzynqPNPI+LHNH9k/dflTwWmvI3VbJlqTCNiD5pfQK/taW5tTMv/SV4NHAx8LDN3uJ1m5raIuAPYp7Rf0dNv1rfTPmrtNXE7fWBEXEWzd3tmZv7L7FXad61/ERFPB34OvCkz19Py97/fMY2IRwAHAZf3NLc6psCHgLcCD97B/Fa31S4f4gMgIvYEvgS8MTN/OWH2ZLfTGNtJ+6yaotbxPs+i+eKf1tN8RGY+mebw2SnlCzdXda4DHlEOWXwEGP/CVDumwJ8D/zbh8F5rY5qZ95TDN4uBwyLikAldqtlO+6gVgIh4OXAozTmecQ8vh0tfCnwoIv54jmv9CnBgORz5TWD8PFOr49rvmNIcIrsgM+/paWttTCPiGGBzZl69k26tbqudDqiI2JXml9PnMvPCSbrs6NZKrd9yqY9aiYjHA+cAyzLztvH2zLyl/NwMfJlm13lO6szMX2bmneX9xcCuETFMpWNanMCEwzdtjmnPOm8Hvk2zN9fr92MXEUPAQ4EtzOGtwXZSKxHxHOBvgGPLIZ3xz4yP6fXls0+ay1oz87ae+j5Nc2If5mhcdzamxc620zbG9Ajg2HIB0XnAsyPisxP6tLqtdjagynHPc4FrM/ODO+i2FnhFRCwoJybvKOcuLgGeGxF7lZOjzy1tc1ZrRDyc5uToiZn58572P4qIB4+/L7X+1xzWuX/pR0QcRrMN3UbPbawiYjeaL9usPaWxz//+RMRDgWfQXHk43tbmmO4bEeNXZD0IeA7w3xO6rQWWl/cvAi4v5y7WAidExO7R3B5sKfDD2aiz31oj4knAp2jCaXNP+17lCi/K/7AcwSydf7wPtR7QM3ksfzjx39r3v8///kTEo4C96Dln0/aYZubbMnNxZh5I8/29PDMnnmNsdVvt8jmoI4ATgZ9ExDWl7e3AwwEy85PAxTSXmF9Hc5n5yWXeloh4D80vVYB3Tzj8Mxe1vpPmWO7HIwL+cOnzfsCXS9sQ8PnM/Poc1vki4DURsQ34DXBC2UC3RcT4bax2AT6Ts3sbq35qBXgh8I3M/HXPZ9sc0wOA1eU8xAOa0vKrEfFu4KrMXEsTtP9cTixvofnlQGb+NCKS5pfSNuCUCYd/5qLWDwB7Al8s4zd+6fNjgE9FxL3ls2dm5qz9Mu2z1tdHxLE0Y7cFOAla//73Uyc0F0ecV75L49oe00nN5bbqrY4kSVXq7CE+SdL8ZkBJkqpkQEmSqmRASZKqZEBJkqpkQEmSqmRASS2KiG9HxKvmug6pCwwoqWMi4qSI+P5c1yHNNgNKup8p91CTqueGKk2h3Dzzo8ArgEcAXweWZ+Zvd/KZZcDfAo+kefLsKRNvpxQR7wIOHr/fWXlsyA3AruVRBifR3AJrX2CU5gmw64BP0tyk906aW2ItLPdsex8QNE9v/TLN4yV+E81TkD9Lc/f5NwGX0twmSqqae1BSf4LmLtQHAY+n3Ndt0o7NTXTXAG8BFgJPB268TytrbmL7YeCozHww8DTgmsy8FvhL4AeZuWeWR4PTPNn0T2ie2nowzbN43tmzyP2BvWkCduV9qUWaK+5BSf358PijDyLiKzRBsCMraG6We2mZvnma67wXOCQibup9gvBE5c7urwYe3/OU2L8DPk/zqPvxZZ3R+3gMqXbuQUn9+UXP+7to7ui9I0uA/x1kZeXu6y+m2VvaGBFfi4hH76D7vsAewNURcXtE3E5zGHLfnj637uyQpFQj96Ckmbce6OfJp7+mCZZx+/fOzMxLgEvKc4TeS/PQvT9j+yeVjtI8+uRxmbmjvTUfW6DOMaCkmXcu8I2I+CrwLZpnAj04Myc+qO4a4LTysMo7+MPhOCJiP+CpwGU04XMnMP58nU3A4ojYLTPvzsx7I+LTwFkR8drM3BwRi4BDSshJneQhPmmGZeYPaR6OeRZN8HyH5uKEif0uBc4H/hO4Gvhqz+wHAG+meWz2FpqnAv9VmXc58FPgFxExWtpOo3kw5xUR8Uvgm8CjZvQfJrXMBxZKkqrkHpQkqUqeg5KmISLeDrx9klnfy8yj2q5Hmo88xCdJqpKH+CRJVTKgJElVMqAkSVUyoCRJVfp/lsmmY/b1tg0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.hist(np.max(z_sample[burnin:].astype('int'), axis=1)+1)\n",
    "ax.set_xlabel('n_cluster')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 構造変化の数と同じクラスタ数=3が多数を占めており、わずかに2個・4個の場合も存在した。\n",
    "  - 初期値はクラスタ数2であったから、正しくクラスタ数を推定できたことがわかる。\n",
    "- MLP本の図（P.96, 図6.4）とは少し異なり、クラスタ数5以上は出現しなかった。\n",
    "  - サンプリングのステップ数を増やせば出現する？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2. 傾き(勾配)の事後分布"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_after_burnin = z_sample[burnin:,:].astype('int')\n",
    "theta_after_burnin = theta_sample[burnin:,:,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "時刻0における傾きの事後分布は、下記のようにして取得できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000,)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad = [theta_after_burnin[ii,z_after_burnin[ii,0]] for ii in range(z_after_burnin.shape[0])]\n",
    "grad = np.array(grad)\n",
    "grad.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.7/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1a1c917630>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmUnFd95vHvrW4trd4XqTdtrdUyxrvxgu0YbBMPMJg5JHecBAIJQZMDCUnIMiQzE5hwJsskE3ASkhMDGSAOAxfHGLPYBhM7sg2WZWNJXiRr37vVi6Repd7qzh9vtaRu9VL7+3bV8zmnT3V1V731u6ern7593/vea7z3iIjI/BcLuwAREckOBbqISIFQoIuIFAgFuohIgVCgi4gUCAW6iEiBUKCLiBQIBbqISIFQoIuIFIjSPL+eLksVEUmPmesB+Q50Tpw4kdLjGxoa6O7uzlE1+aE2RIPaEA3zvQ1h1N/S0pLU4zTkIiJSIBToIiIFQoEuIlIgFOgiIgVCgS4iUiAU6CIiBUKBLiJSIBToIiIFQoEuIlIg8n6lqIgUjviWxy/5Wuz2e0KoREA9dBGRgqFAFxEpEAp0EZECoUAXESkQCnQRkQKhQBcRKRAKdBGRAqFAFxEpEHNeWGSt/Sfg3UCnc+6KxNfqgG8Aq4FDgHXOnc5dmSIiMpdkeuhfBqZe+vVJ4EfOufXAjxL3RUQkRHMGunNuC3BqypfvBb6S+PwrwHuzXJeIiKQo3bVcGp1z7QDOuXZr7bKZHmit3QxsTjyWhoaG1AosLU35OVGjNkSD2pB9QxUVl3xtyRz1Ra0NqYpy/TlfnMs59wDwQOKu7+7uTun5DQ0NpPqcqFEbokFtyL74wMAlXxuao76otSFVYdTf0tKS1OPSneVy0lrbDJC47UzzOCIikiXpBvqjwAcTn38Q+HZ2yhERkXQlM23x/wF3AA3W2mPAp4A/B5y19sPAEeDnc1mkiIRrunXPJXrmDHTn3C/M8K07s1yLiMwT3nvYvxvGRjGXXRl2OZKgHYtEJCX+zCl4/mno6gjuL1iIWXtZuEUJoEv/RSQF/uwQPPYQ9J2BW94OTa3w/L/je7rCLk1QoItIKnbtgPFx+Nn/FPTKb70bFi+Gf38cPzIcdnVFT4EuIknxI8Ow5zVYuQZTXQuAKVsShPpgPxzeH3KFokAXkeTseQ1GR+BN107++rJmKK+EowfDqUvOU6CLyJz8+Bjs3gnNyzH1Syd9zxgDK9qg4xh+dDSkCgUU6CKSjEP74OzQpb3zCctXB2Pr7UfzWpZMpkAXkbkd3h8MqzS1Tv/9xmZYuAiOHcprWTKZAl1EZuWHz0HHMVixOhhemYaJlUDLSjh2CB8fz3OFMkGBLiKz27U9GE5Z3jb741ashuFzsP+NvJQll1Kgi8is/PYXYMHCYFhlNi0rIRbD79ian8LkEgp0EZmRj4/jd26D1pXBsMoszMJF0NCI37crT9XJVAp0EZnZgT3Q3zv3cMuEuqVw9AB+XOPoYVCgi8iM/I4XoKQEWlcm94T6ZTAyEpxElbxToIvIjPwrL8L6NwXDKclIXHTkD+3LYVUyEwW6iEzLDw7A8cOYjW9O/klVNbCoDA4r0MOgQBeR6R3YDYBZf3nSTzHGwKo1eAV6KBToIjItv/f1YPx89YaUnmdWroOjB3ViNATasUhEzrt471D/0o+hph6/9anUDrJ6HTw5Au1Hkp8dI1mhHrqIXMKPj0PPSVjalPJzzaq1wTG0PnreKdBF5FKnuoPL/ZfNcXXodJa1wOKyYIVGySsNuYjIpbrag9s0euj+2R9AdS3+1ZcmDeHEbr8nW9XJDNRDF5FLdXVARRVmSXl6z69bCqe6tfJininQRWQS7z10tsOy1Hvn59XWQ3wc+vuyV5jMSYEuIpMN9MG5s7A0jfHzCTV1wW3v6ezUJElRoIvIZKe6gtv6Zekfo6o2uO09lXk9kjQFuohMdroHjIGa2rQPYRYsgPIK9dDzTIEuIpOd7oGqGkxJhpPgqmsV6HmmQBeRyU73BCc1M1VVC71ngpOskhcZ/Qm21v4O8GuAB14BfsU5dy4bhYlI/vmRYRjshxQW5JpRTR2Mj8FAP1RWZX48mVPaPXRrbSvwceB659wVQAlwX7YKE5EQnEmcxKxryPxY1Ykx+D4Nu+RLpkMupUCZtbYUWAKcyLwkEQnN6Z7gtiYLQy4Tga5x9LxJO9Cdc8eBvwKOAO1Ar3PuB9kqTERCcLobFi6CdK8QvYhZtDhY0+WMpi7mS9pj6NbaWuBeoA04A3zTWvt+59yDUx63GdgM4JyjoSG1f+VKS0tTfk7UqA3RoDbMravvDKZhGUsqK7NyvKG6BvxAH+UVFSxJ1D3ffw5Rrj+Tk6J3AQedc10A1tqHgVuASYHunHsAeCBx13d3d6f0Ig0NDaT6nKhRG6JBbZidj8eJ93TB2ssYGBjIzjErquDgXvr7+xlK1D3ffw5h1N/S0pLU4zIJ9CPATdbaJcBZ4E7gxQyOJyJh6umEsdHsTFmcUF0LoyNwdih7x5QZZTKGvhV4CPgpwZTFGBd64iIy3xw7FNxmO9BBJ0bzJKN56M65TwGfylItIhIiPxHo1XXZO2i1FunKJ10pKiKB9qPBGugLFmTvmGVLoHQB9J3J3jFlRgp0EQHAd7ZDVXVWj2mMCY7Z35vV48r0FOgicmFTi8rsBjoQHLNfPfR8UKCLSLDeytnB3AR6VQ0M9OPHRrN/bJlEgS4iFzaFrshRD9176D6Z/WPLJAp0EQnGzwGqcrAqYlVNcHtSSz3lmgJdRILxc2OgIneB7hXoOadAFxHo6oDa+sx3KZqGWbQ4WPBLgZ5zCnQRwXe1w9Lm3L1AVQ2+U4Geawp0EYHOdsyyHAZ6ZTUo0HMu+/9fiUjkxbc8fv5zPzIC/b34gT5Mrl6wqgYO7sEPD+fqFQT10EVkIHEVZy7moE+YOPbE9EjJCQW6SLHrz0Oga+piXijQRYpdf19wW5mDKYsTEmvE6MRobinQRYpdfy8sLsMsWJizlzALFgZro6uHnlMKdJFi19+b2+GWCcuadXFRjinQRYpdngLdLGvR1MUcU6CLFDE/NgZDOVplcarGVug7Q3xoMPevVaQU6CLFbGDihGgeeuiNwYVL4+1Hc/5axUqBLlLMzk9ZzOEMlwmNrQCMn1Cg54oCXaSY5WMO+oSlTQCMKdBzRoEuUsz6e2HhomBFxBwzCxdB3VINueSQAl2kmPX35We4ZUJji3roOaRAFylm+ZqDnmCWNauHnkMKdJEi5cfHYbA/r4FOYyt+oB8/MbtGskqBLlKsBvuDzZvz2kNvCT7RFaM5oUAXKVb9+ZuDfl5jEOhaAiA3FOgixSqfUxYnNDRCrEQ99BxRoIsUq/5eKC2FxWV5e0lTWkpJY7PWdMkRBbpIsUrMcDEmZxvPTaukeYXWRc+RjPYUtdbWAF8ErgA88KvOuZ9kozARybH+Xqipy/vLljQvh9dexnuf9z8mhS7THvr9wOPOucuAq4BdmZckIrnm4/FgYa6KPF5UlFDashKGz0Hv6by/dqFLu4dura0Cbgc+BOCcGwFGslOWiOTU2UGIx/N7QjShpGV58MnJE6H8h1DIMhlyWQN0Af/XWnsV8BLwW845LXYsEnV9IcxwSShtWQkE+4uajVfk/fULWSaBXgpcC/ymc26rtfZ+4JPA/7j4QdbazcBmAOccDQ0Nqb1IaWnKz4katSEa1IYLzoyeYxgob2omVlGReWFJWJKou8QYKF1AWf9pKufhzyPK76NMAv0YcMw5tzVx/yGCQJ/EOfcA8EDiru/u7k7pRRoaGkj1OVGjNkSD2nDBeHcXxGIMeoMZGMhCZXMbePhBACoqKqC8gqEXf8zZimpit9+Tl9fPljDeRy0tLUk9Lu2Tos65DuCotXZj4kt3Aq+nezwRyaP+XqiowsRCmrlcWXNh2EeyJqNpi8BvAv9irV0IHAB+JfOSRCTn8rzK4iWqquHEEbz34dVQgDIKdOfcduD6LNUiInngvQ8CPbElXCiqaiA+DoP5Ge4pFrpSVKTY9J2BsbFwe+gTr91/JrwaCpACXaTYdLYHt/ncqWiqqprgVuPoWaVAFykyvmsi0EPsoZctCRYG61MPPZsU6CLFprMdjIHyytBKMMYEM1005JJVCnSRYtPZDuWVmJKScOuoqtaQS5Yp0EWKjO/qCHf8fEJlNQz048fGwq6kYCjQRYpNZ3u44+cTqmrAx6GnM+xKCoYCXaSI+MF+GBqIRqBP1KDNLrJGgS5STDo7gtsoBHpi6qI/eTzkQgqHAl2kiJzf+i0Kgb5oMSxcpA2js0iBLlJMJuagh7BT0VTGGKiqwSvQs0aBLlJMOjugph5Tmum6fFlSVQPtx8KuomAo0EWKiO9qh2VNYZdxQXUtnOnBnxsKu5KCoEAXKRLee2g/hglzlcWpJtZ06dCJ0WxQoIsUi/5eGOyHlhVhV3JBdS0AvkPDLtmgQBcpFomxatMUoUCvrIKSEmhXDz0bFOgiRcK3Hwk+iVAP3cRKYGmTeuhZokAXKRbtx2BRGdRGbMf6puXQfjTsKgqCAl2kSPj2o9C8PJj/HSGmaTl0tuPHx8MuZd5ToIsUi/ajmOblYVdxqablMD4G3SfDrmTeU6CLFAE/NAhnTkHzyrBLuYRpSkyj1Dh6xhToIsUgMUYd2R46mrqYDQp0kSJwPiwjNMNlgimv0BIAWaJAFykGJ45C6QJoaAy7kuk1LVcPPQsiskJPNMW3PD7t12O335PnSkQy49uPQmNLMO87gkzTcvyLz+K9j9wsnPlEPXSRYtB+FNMSvROi5zUvD3ZS6j8TdiXzmgJdpMD54eFg386mCJ4QTTATY/vHj4RbyDynIReRQtdxDLyH5uidEIVgaNMPDQafP/MEJrEJh4Y2U6ceukiB84f3AWBWrQm5klmULQm2oztzKuxK5jUFukihO7wPlpTD0uawK5mRMQZq6qD3dNilzGsZD7lYa0uAF4Hjzrl3Z16SiGSTP7QXVq2L/uyR6lo4vF8zXTKQjR76bwG7snAcEckyPzoCxw9jVq8Lu5S51dTByDCc1XZ06coo0K21y4F3AV/MTjkiklVHD8L4OGb1+rArmVt1XXDbq3H0dGXaQ/8c8AdAPAu1iEiWTZwQZdU8CPSaRKDrxGja0h5Dt9a+G+h0zr1krb1jlsdtBjYDOOdoaEhtcf3S0tKUn5MtQxUV0359yTxqQ7aoDdGQaht6O44yUl1Lw4bLJo1Lz/TezoeSWAkV07y+Ly9ncHEZpYP9LK6oSPn3LF+i/D7K5KToW4H3WGvfCSwGqqy1Dzrn3n/xg5xzDwAPJO767u7ulF6koaGBVJ+TLfGBgWm/PjSP2pAtakM0pNqG8d2vwsq19PT0TPr6TO/tfKioqGBghtf3VbWMdncyNjCQ8u9ZvoTxPmppaUnqcWkHunPuD4E/BEj00H9vapiLSHj8ubPQfgxz3S1hl5K8mlo4uBfvfdiVzEuahy5SqI4cAB/HzIfx8wk1dTA6AmcHw65kXsrKpf/OuaeBp7NxLBHJjvMnROfDlMUJ1Toxmgmt5ZKG6ZbV1boTEjn7d0NdA6a6NuxKkqeZLhnRkItIAfLxOH7Pq5iNbw67lJSYxWVQVg6ne+Z+sFxCPXSRAnHxf47+dA/09+Ln4yX0tfVwOpozXKJOPXSRQjSxnVtja7h1pKOuAc6cxo+Ohl3JvKNAT5L3PpgGJjIfdJyAiipMRWXYlaSutgF8HNqPhl3JvKMhlzn4vjOwc1vwC3J2EP/WOzFrNoZdlsiMfDwOJ4/DqrXAzHvjRlZtPQD+6EHMygiv4R5B6qHPZdszwQJHy5qhqgZeeSn4hRGJqtPdwVzupnk43AJQWQ0lpXDsYNiVzDsK9Fn4gX44cRQuvxpz+zvgqrdA35kg4EWiquN4cDsfx88BE4tBbT1ev2cpU6DPZn9imfe1m4LblWuC3sOrL+nSZImujuNQXYtZUh52JemrrYejB/R7liIF+gx8fDy4MKN5xfkTSyYWgyuuhVPdQc9dJGL8+Dh0npi3vfPz6hpgaDD4XZOkKdBn8voOGByAdZsmf71tQ7A/4+6d4dQlMpvOdhgbg5YVYVeSmdrE8rRHD4RbxzyjQJ+Bf/aHsGgxrGib9HVTUhIMvZw8EfSGRKLk+GGIxaBpediVZKamHozROHqKFOjT8KMj+B0vwOr1QYBP1dgK42PQfTL/xYnM5vhhaGzFLFgQdiUZMQsWwLIWvGa6pESBPp3D+2BsdOZpX42JxeZPnshfTSJz8P29wSys1lVhl5IVZkUbHN4fdhnzigJ9Gn5vYnbLsuZpv28WLQ7G+Camh4lEwfHDwW2BBDpt66GnE993OuxK5g0F+jT8/l3Bv62Ly2Z+UFMrdHXgx8fyV5jIbI4fhqoaTFV12JVkxfkrsg/sCbeQeUSBPoX3Hvbvwqy7bPYHNrZAfFzj6BIJfvhcsDxFy8qwS8melWuhpAR/UIGeLAX6VB3HYaD/wsVEM2lsAWM07CLRsHtn0MFYXiDDLYBZuAhaVyvQU6BAn8Lvex0As+7yWR9nFi5KjKPrxKiEz+94ARIzQwqJWbMBDu4JLvSTOSnQp9q3Cyoqk1vYqKkFujvwYxpHl/D4+Dh++1ZoWTX9NNv5rG0jnDsL7fpPOBkK9Cn8vl2wdhMmmZ1elrVAPA6nunJfmMhMDrwB/b2XXARXCMyaDQD4g2+EXMn8oEC/iO87A50nMHONn0+oXxrc9ijQJTz+5a3BcrOtBXRCdMKylmCpDY2jJ0WBfrEDuwEwU9dvmUlZOZQtUQ9dQuO9x29/Hja+OTivU2BMLAZtG/AH1ENPhnYs4sKOLn77VjCG+NEDmPYjcz7PGIOvWwo9nbkuUWR67Uehsx1z971hV5Izpm0j/nsOf+7s7NeGiHrok/R0BetIl6bwd65+KfSdCeYBi+SZf/l5AMxVN4ZcSe6YNRuDPUY17DInBXqC9z4YOqlbmtoT65aC91rmU0Lhd7wAbRswiX04C9L6TcEFRrt2hF1J5CnQJ5wdDKZHpRro9csA8FpESPLMn+6Bg3swVxdu7xzALF4CbRsV6ElQoE/oSeyMUp9aoJsliROjh/floCiRmfkdWwEKPtABzKar4PA+/GB/2KVEmk6KTjiVOLE5sVNKKuqW4g8p0CW//Mtbg2l9zfN8d6IZTExWAPAj54Khzd2vwHW3hFhVtKmHPuFUd3BCNJ2NAeqXQsdxnRiVvPFDg/DGK5hrbkzuIrj5rmEZLFiA37U97EoiLe0eurV2BfBVoAmIAw845+7PVmF5d6rrwsYVqapbFpyFP3oA5lgDRiQb/KsvwfgY5uqbwi4lL0ysBN/Yin9dgT6bTHroY8DvOuc2ATcBH7PWzss082eHgh3G65ald4DEuLtOjErebN8KldWQuDS+KDQvD/Yg6OoIu5LISjvQnXPtzrmfJj7vB3YBSaxoFUETV3rWpzF+TuLEaHUtaBxd8sCPjuJfeRFz9Y2YWIEtxjWb5mDja797Z8iFRFdWxtCttauBa4Ct2The3k2sxVKb4pTFi61ci9dMF8mHXdvh3NmimN0ySVUt1Dbgd24Lu5LIyniWi7W2AvhX4Ledc33TfH8zsBnAOUdDQ2q94NLS0pSfk6qevtOMV9dSUVeX9jHil1/J4De/TH1F+SWXJ+ejDbmmNkRDaWkpC1/ZxnB5JQ233TXpJP5QRUWIlSWvJFZCRZq1jt96J0OPPUxd2SJi5ZVZriw5UX4fZRTo1toFBGH+L865h6d7jHPuAeCBxF3f3d2d0ms0NDSQ6nNSNd7ZDkubGBgYSPsYZmmwlG739m2XbI6RjzbkmtoQDfWVlZx77kewah1d3/lG2OWkpaKiIu3fNfPmG+A736D7ye8Re+tdWa4sOWG8j1pakpuwkfaQi7XWAF8Cdjnn/jrd44TN9/fB4EDqV4hOtWptcDydGJUcGv7pT2B0FFavC7uUcKxeD0ub8C9sCbuSSMqkh/5W4APAK9baiblEf+Sc+37mZeXRkUQAZxjopqZeJ0Yl5849+yQsLoPG+Tn/IFPGGMwNt+MfewjfexpTXRt2SZGSdqA7554F5v0VDf58oGdhTGzVOp0YlZzx54YYfvE5aFsfrBNepMxbbsd/3+FffA5z57vDLidSivddkeAP74OKKsyixRkfy6xaqytGJWf8jm0wMhwMOxQx07oSWlfht2nYZaqiD3SOHMh8/DzBrFp34YpRkSzzzz9FrKERljaFXUrozI13wP7d+GMHwy4lUoo60P3gAHR1pLzC4ox0YlRyxPd0wWsvU/a2dxbH2i1zMLe/AxYuwj/xSNilREpRBzrZHD9n4sRonU6MStb5554EoOwujRkDmPJKzG3vwG/bgteevucVdaD7LM1wmWSVrhiV7PLx8SDQL7uSkmXNYZcTGebue8F7/JOPhl1KZBR1oHN4P9QtzerGs6ZtPXQcww+lf5GSyCS7dsKpLsxt7wi7kkgx9cswN9yG3/KDYPhUinuDC394//lx72wxazcF+5MeeAOuuC6rx5bi5J/9IZRXFs1SubO5eNMLAN/QCMNn8T98BPPe94dUVXQUbQ/dnx2CzhOYldkNdNo2QCyG37cru8eVouTPnMK//DzmpjvS23ylwJm6huD6jx8+orF0ijjQOfAGAKYtu+tJm8VlsLxNgS5Z4f/tuxCPY96uk6EzuvZmiHv8w18Nu5LQFW2g+z2vQSwGay/L+rHNuk1wcA9+bCzrx5bi4c+dxf/oO7CiDb97B/EtjzP0A03Tm8pUVGLuvhe/9d/xB/eEXU6oijfQ974Kq9Zl9YToees2BVf06aIHyYB/9ofB++jyq8MuJfLMO38OqmqIf/0L+Hg87HJCU5SB7kdH4OAezPo35eT4Zu2m4HU07CJp8uPjwXS8Zc2YpY1hlxN5/oUtcMW1cOAN4v/0WeJbHr/kBGoxKMpA58AeGBvDbMhRoNc1BHPbFeiSJr/tGejpVO88FWs2BqtQvvx8MOmhCBVloPu9r4IxsC53e1qbdZvw+3cFUxhFUuDHRvGPfg2Wt8Hy1WGXM28YY+DG22FsFF58LuxyQlGcgb7nNWhdjSnP4ZZd6zbBmVNBL0skBX7LE9DVQex9v6x1W1JkqmuD6z8O7cWfOBJ2OXlXdIHux8Zg/+6cDbdMMOuD3r9/45Wcvo4UFn9uCP/db8DGN8Obrg27nPnpimuhqga2bsGPDIddTV4VXaBzZD+MDOc80GldDXUN+O0v5PZ1pCBMnMSLf+H/QH8vrNmAf+aJsMual0xJCdz4MzDQh/+eC7ucvCq6QPdvvBp8sj534+eQ2CrrqrfA6y/jh4urlyDp8X1n4LWfwqq1mAbNbMmEaWqFNRvxT3wLf7x4hl6KL9Bf/gmsXIOpyv1ehObqG2FkmJGd23L+WjK/ee/h+achVgLX3xp2OYXhultgcRnxr/4tPj4edjV5UVSLc/mujmD++fs+mPVjTzfn1Y+PQ9kSzr3wDLRl/4pUKSD7dsHJE3DTHZgl5WFXUxDM4jK479fwX/os/snvYN7x3rBLyrmi6qH7xFQmk6cekCkpwVxxHSPbni2aHoKkzvd0wks/hsaWYHaUZI258Q648gb8Iw/iT54Iu5ycK65A37YF1mzM7/jk1TcS7z0dXMwkMoUfGSb+938a3Ln5bZqmmGXGGGLv/yiULiD+5b8p+I5V0QS67zgGRw9ibsjv+KS54jooLcVvfz6vryvR573HP/j3wUblt96FqawOu6SCZGrrMfd9BPa9HkwJLWDFE+jbngVjMNflOdCXlLPwqhvwP3kqWENGJME/8TD+J09h3vOLGF0RmlOxW96Oufnt+O9+A79rR9jl5ExRBLqPx4PFe9Zfjqmtz/vrl9/7i9B3Bv/cj/L+2hJN8e85/L9+BXPDbZh32bDLKQrml34dmpYT/8Jf4c/0hF1OThRHoG97BjqOYW772VBef8EV10LbBvwPvhXMfJGi5eNx4t/6Z/wjD2Ju/BnMhz+BiRXFr2HozKLFxP7Lf4WRYeKf+zR+sD/skrKu4Kct+tFR/Lf+OdgkYPhsKEtqGmOI3fM+4v/wZ/iXnsO85fa81yDh86e6iP/T5+CNVzC33o35wEcxsZKwyypo0/6+33Y3PP0Y8fv/J7FPfCY3eyKEpPAD/anvQU8nsV/+E3x3R3iFXH0jNC3HP/YQ/vpb1SsrcBcHiR8bg72vwc4Xg/M4v/wbQaBrRksoTPMKzEd+n/g//gXx+z9N7Nc/GSzqVQAKOlX8YH+wlsObrsGEvK60icUw7/p5OHYI//3iWl+iWPmzQ/jXXoZHHgyWc61rIPbHnyN22zsU5iEz195M7CO/B0f2E//M7+D3vh52SVlRsD10f3aI+N9+BobPEnvfh8IuB0hc5PDay/hvfw3fuhpzzU1hlyRZ5s8O4Xduwz/9GBw7DD4eXDB0692Yplb87p343TvDLlMILjCMNS0n/g9/Rvyv/gjz1rsw774v2KBmnsoo0K219wD3AyXAF51zf56VqjLkzw0Rv//TcGgvsc1/gFnRFnZJQGIB/g98DN9xnPiXPkvsdz+DadsQdlmSAe89HD+Mf/3loDe+51UYG4OyJbDpSlh7GaamLuwyZYpJY+tveydsfwH/3JP4H/8I2jYEH40twfmv2+8Jr9AUpR3o1toS4PPA3cAxYJu19lHnXGj/u/j4OH7bs/jvfB262oMwv/bmsMqZllm4iNhH/4j4n/0+8T//A8yd/zGYh1xAJ2YKmR8ZDgL88H448Ab+9e3Qeyr4ZvMKzNvehbn2ZuLHD2tYZZ4wCxfBW27DX34VvPISHNobrK2zuAzf1Eo87jEr10DLCsziJWGXO6tMeuhvAfY55w4AWGu/DtwL5CzQvfcwPp74GIXBAeg9je9shwO7g1+urg5oXUXsN/8Yc0X4GwTEtzzOUEUF8YGBSV+Pfepv8A9/Ff/Db+OfexKz6epgrL+xFWrrYUkFlC6A0lKIxRQOszi/zd/57f78pJuJ7/mxMfzY6KXfm7gzPg7nhmBoMPg4O4jv64Wek9B1Er97B/SevvA6ixZDUyvmvb+EufxqTN3S80c0Rbhbzny1mNs0AAAFeklEQVRnKqrg5rfhr78Vjh0KPk6ewP/LP1x4u9Qt5fSqtcTrG6GmDqpqMFU1UFkNFZWwYCEsWBD87paU5v33NpNAbwWOXnT/GHBjZuVMb/zzfwqvvAjjYzM/qGwJtG0k9r4PwjU3R34WiSmvwHzgo/hb3o5/5ong3/WXnmPaHUiNgZJSiKX75kjjeVPeiJ3GJL8/6qTHTQnXqaE76dOpwTzH91KU9maAxkB1XfBHdnkb1Cc2AS+vDP4lv/XutGuS6DELFkDbemhbj/ee2OXXwPFDwbrqJ44S7zqBf+2nMBJc+T3rO7KkNHj/GEPsjz+HaVqe09ozCfTpUuKStllrNwObAZxztLS0pPxCK/7X36X8nGnd96vZOU4aamb6RksL3Pb2fJYi+ZDD99qM76V5ZN614arw/9tPRibd2GPAiovuLwcuWZ/SOfeAc+5659z1BH8EUvqw1r6UzvOi9KE2RONDbYjGx3xvQ4j1zymTHvo2YL21tg04DtwH/GIGxxMRkQyk3UN3zo0BvwE8AewKvuRey1ZhIiKSmozmoTvnvg98P0u1zOSBHB8/H9SGaFAbomG+tyGy9ZukZy6IiEikRXtun4iIJC0ya7nMtYyAtXYR8FXgOqAH+M/OuUP5rnM2SbThduBzwJXAfc65h/Jf5eySaMMngF8DxoAu4Fedc4fzXugskmjDrwMfA8aBAWBzmFc4T5XskhrW2p8Dvgnc4Jx7MY8lzimJn8GHgL8kmFAB8HfOuS/mtcg5JPNzsNZa4NMEU7Z3OOdCnRgSiR76RcsI/AfgcuAXrLWXT3nYh4HTzrl1wGeBv8hvlbNLsg1HgA8BX8tvdclJsg0vA9c7564EHgL+d36rnF2Sbfiac+7NzrmrCer/6zyXOaMk68daWwl8HNia3wrnlmwbgG84565OfEQtzOdsg7V2PfCHwFudc28CfjvvhU4RiUDnomUEnHMjwMQyAhe7F/hK4vOHgDuttUnNzcyTOdvgnDvknNsJxMMoMAnJtOEp59xQ4u7zBNcfREkybei76G45c1zsl2fJ/C4AfIbgj9G5fBaXpGTbEGXJtOEjwOedc6cBnHNpX4ycLVEJ9OmWEWid6TGJKZO9QP43CJ1ZMm2IulTb8GHgsZxWlLqk2mCt/Zi1dj9BKH48T7UlY876rbXXACucc9/NZ2EpSPZ99D5r7U5r7UPW2hXTfD9MybRhA7DBWvuctfb5xBBNqKIS6NP1tKf2mpJ5TJiiXl8ykm6Dtfb9wPUE46BRklQbnHOfd86tBf4r8N9zXlXyZq3fWhsjGHL83bxVlLpkfgbfAVYnhu6e5MJ/31GRTBtKgfXAHcAvAF+01oa6qkFUAj2ZZQTOP8ZaWwpUA6fyUl1ykloKIeKSaoO19i7gvwHvcc4N56m2ZKX6c/g68N6cVpSaueqvBK4AnrbWHgJuAh611l6ftwrnNufPwDnXc9F75wsEkx2iJNlM+rZzbtQ5dxB4gyDgQxOVWS7JLCPwKPBB4CfAzwH/5pyLUg+4EJZCmLMNiX/3/xG4JwpjhtNIpg3rnXN7E3ffBewlOmat3znXC5zfUsda+zTwexGb5ZLMz6DZOdeeuPsegqvNoySZ3+dHCHrmX7bWNhAMwRzIa5VTRKKHPtMyAtbaP7HWvifxsC8B9dbafcAngE+GU+30kmmDtfYGa+0x4OeBf7TWRmqphCR/Dn8JVADftNZut9Y+GlK500qyDb9hrX3NWrud4L30wZDKvUSS9Udakm34eOJnsIPgHMaHwql2ekm24Qmgx1r7OvAU8PvOuZ5wKg7oSlERkQIRiR66iIhkToEuIlIgFOgiIgVCgS4iUiAU6CIiBUKBLiJSIBToIiIFQoEuIlIg/j8Cu4pJU6wpoAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(grad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "この結果を全時刻で統合する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "grad_all = np.empty((T,n_step - burnin))\n",
    "for tt in range(T):\n",
    "    grad = [theta_after_burnin[ii,z_after_burnin[ii,tt]] for ii in range(z_after_burnin.shape[0])]\n",
    "    grad_all[tt] = np.array(grad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XmUZNld4Pfv22PPfa29uqp6X9XqRd1Sd6NGaiGQNDA8CaEZMcDIwzEwGDOeGY9HZmR7jgbbAxh0bDQakGHA8IRBtIWkBkm9qdfqtaqrqqu69srKzMo99oi3+o+XGZlZGVm5RsTLrPs5p09XRryMd/Nl5PvFvfd3f1cKggBBEARBiBq51Q0QBEEQhHpEgBIEQRAiSQQoQRAEIZJEgBIEQRAiSQQoQRAEIZJEgBIEQRAiSQQoQRAEIZJEgBIEQRAiSQQoQRAEIZLUVjegwUSZDEEQhGiSVjpguwcohoeH13R8d3c3ExMTDWrN9iCu0eqI67QycY1Wth2v0eDg4KqOE0N8giAIQiSJACUIgiBEUmSG+EzTfAL4XUABvmZZ1pfrHGMCv0k4t/S2ZVmfbWojBUEQhKaJRA/KNE0F+ArwMeAW4GdM07zlqmMOAv8aeMiyrFuBX2t6QwVBEISmiUSAAu4DTluWddayLBv4c+CTVx3zT4GvWJY1DWBZ1liT2ygIgiA0UVSG+HYAlxZ8PQTcf9UxhwBM03yBcBjwNy3L+m5zmicIgiA0W1QCVL18+KvXMKnAQeBRYCfwvGmat1mWNbPwINM0vwB8AcCyLLq7u9fUEFVV1/w91xtxjVZHXKeViWu0suv5GkUlQA0BuxZ8vRO4egHTEPCyZVkOcM40zZOEAevwwoMsy/oq8NXZL4O1rh/YjmsONpu4RqsjrtPKxDVa2Xa8RqtdBxWVAHUYOGia5j7gMvAZ4OoMvW8CPwN83TTNbsIhv7NNbaWwKkEQ4FU8nKJDZaZCajCFltBa3SxBELaYSAQoy7Jc0zR/GXiKcH7pDy3LOmaa5peA1yzLenL2uY+Ypnkc8IB/YVnWZOtaLcyp5quURkuUJ8o4BQe35BL4AQEBiqow+c4k7Qfa6bqlC0lesbqJIAgCAFIQbOtydYEodbQxlekKU+9O4ZZc2va3kdmToae3h4mJCcqTZSaOTFCZqaDoCoqmLPs6ru2iqAq99/SS6E008SdoHfFeWpm4RivbjtdodohP1OITVs/3/HBoruRQnihTGCpgF2y0uIYkS0wcnWDy2CT2zTZjF8eozlRREyp6Ul/xtVU9fKuNvDLC/h/fjySJnpQgCNcmAtR1rpqtkj2bpXSlhFtxw9zJAGRDRtGURcFHjYVvl8JoAc/20JJrn1cK3IDKVIV4V3yzfgRBELYpEaCuQ9V8lezpLMXRIm7FRY2pyIq86kQGWZHX3QNS4yrZc1kRoARBWJEIUNeRmbMzZM9ksQt2LSitZnhuM0myRGW8Uve5wA9EEoUgCDUiQK1g/O1xIPzkr8ZVkoNJZCUqFaJWL/ADpo5NoRhK04PS1ZyKQzVXxcgYix4feXmEnjt71jV0KAjC9rP17rRNljuXozhSZPq9aa68doWJI9HMpgn8gMkTk0ydnKr7fHmyjO/6TW5VfWpMJXsmu+gxp+SQu5CjNFZqUasEQYgaEaBWQZIlFE1BS2rkLuSwi3arm7RI/nKeC09dYOb0DNmz2brH5M7nUOPR6DDLikxpfHEgmjgygdFuUJ4ot6hVgiBEjQhQa6TGVMZfH1/yuFN2GHtzjOJoEd9pTk/FKTlceuYSY6+NISkSqqHiVTwqU0vneMoT5UjN77gFF6fkAOHPURwtIisydj5awV8QhNaJxkfqLUSSJcpTZYqjRZL9SWA2UDx9CUmSyJ3PgRwGsv77+4m1xxrSjunT00wem0Q11EU9IzWuMnN6hv77+muPVXNVvIqHnIzO5xElppA7l6Pr1i4mjk6gGOEiX6foEASBWCclCILoQa2HltAYf2ucIAjC4PSDS8iKjKzKaEkNLR5O8l959QqbXanDKTtc/MFFpo5P1RbQLiTJEqWxEoE/f97suWwtAESFrMoUR4q4ZZfiSLGWeOK7Pm7ZbXHrBEGIAhGg1smzPSaOTITBSZWXBgpJwq24TL1bP2lhvabfncarerVFs/X4rk9pdH6OpzxeRlaj96u28zZXXr+yKHjKskxlon4a+rUsDMiCIGwP0btrbRFzmWj1gtPCY2ZOzdTmWjZDdaa6Ypq7GleZORtuk+VWXJz85p1/M8m6TOlKadHPo8SUdWXyXXrm0soHCYKwpYgAtQFaaukQ29UUQ2H01dFNOV/gB6tKIpAkicpkBd/1yV/MI2vR/DUrmoKeXrwmS5Ik7MLaEiWq2SqFoYIYGhSEbSaad65tRJIl7KxN9lz99O+1qGarBN7qh7Lyl/IUhgsoerTmn1biFNfW48uezWK0GeSH8g1qkSAIrSACVBOocZWJtyeoZqsbep3C5cKqkx3UuEr2bHbD52wF316aKOG7PpXp+nNTpYkSalxdNO8mCMLWJwJUkyhxhUvPXNrQOp/KZGVNyQ521kZaecuV6JFYspZr+tR03aFSt+zi5sNgVs1WNz1rUhCE1hEBqkkkSUKNqVx65tI1kyaCIKhbkigIVjf/tJCW0iJTPWIt1JhK6cp8byjwA7LnszglZ0ngyp7PIhuzKeq2vyV7jIIg1Lf17l5bmCSFJZMu/eASAx8YwM7ZtW3SvaqH7/oETtgD2Pfj+xb1lpyCg+d4W24+aT0kWVoUjHPnc/iOj5bQmDw+yY6Hd9SeK42Wajv5KnGF3IVcwxZHC4LQXKIH1WSSLCGrMkPPDIXbpU9W8G2/FrzUhIqkSEy/N73o+7ZissNGLEyUmH5vOlyULEmUx8u1+SnfWdxjkhWZyuTa11AJghBNIkC1gCRL6CkdNabWLemjGGEZoIXzKeXxcq2ncD3wqh5uxQ03VSzNJ0wohsLk8UkgLJJ79fWz83ZkqrYLgrAxIkBFlO/45C7kal/bueusiKoE1ekqUyemFu0PJasyhcsFfNenMFRYOscWQHlMVEQXhO1ABKiIUuMqM++F1SCcUjhHdT1RY2HR28rM0iE7SZKYeW+mbtV2Na6K9VCCsE2IJIkIc4oOpfESds5GUrdguvgGSLJE9kKWRHdiyXOKoTB1YqruxytJkpZdLyUIwtYiAlSEqXGVyeOTYfKEcf39qhI9iWW33VBiCpJS/zmn6OBW3GsW1BUEIfrEEF+ESZJEdapadyjrenCtPaFkVV4+eGlKuC+XIAhbmghQEafGVZGVtkaKoTD17hR28TpLLBGEbUYEqIiTZAktoa18oLCIGlMZ+eGI2CdKELYwEaCEbUmSJVzHZfzIeKubIgjCOokAJWxbqq6SO5dbVNdPEIStQwQoYVvTkhqjh0fxPTGPJwhbjQhQwrbnu77YK0oQtiARoIRtT42r5C6KtHNB2GpEgFolbeo8nS9+DckV+w1tNZIkUZ0WvzdB2GpEgFoFNTtMz3P/J4nLb6NPnmt1c4R1cCuu2MxQELYYEaBWoJYm6HnuKwRyeKm07EiLWySshxpXyZ7LtroZgiCsgQhQ1xBkp+l75Q/Adxl/5Ffx9BRqTgSorUhWZMoTYhsOQdhKRIBaRlAq4P/O/4hSzTH58C/htg3gtvWLHtQW5uSvv21LBGEri0y5Z9M0nwB+F1CAr1mW9eVljvuHwDeA91uW9VrDGqRq0DvI2OATuF17AXAyAyQuHIYggGsUMhWiSVIl8kN52m9ob3VTBEFYhUj0oEzTVICvAB8DbgF+xjTNW+oclwZ+FXil0W2SdAPll/4VlZ4ba485bQPIbgWlPN3o0wsNoBoqxeFiq5shCMIqRSJAAfcBpy3LOmtZlg38OfDJOsf9T8BvAS3Zf8LJDADXWaKE52CMnUTNj4U9xy2uMl0RBWQFYYuIyhDfDuDSgq+HgPsXHmCa5t3ALsuyvmWa5m8s90KmaX4B+AKAZVl0d3evqSGqqi76nivpK+hJPfzCOARAqjqNnsms6XWjSspdQT/+PdSRd/HbBvC69+J170WqFtDOvYY69DaSE34e8JOduIO3EOy4hfauffiZfpCj8hlndWzZJu7FSfWmGn6uq99LwlLiGq3ser5GUQlQ9SZ0ah9zTdOUgd8Gfm6lF7Is66vAV+deY2JiYk0N6e7uZuH3FPIFVG/+MiVjGbyxc+RyW7sygXHlXVKnniY+epxAUqh270e98h6xc/Ojp56RorjrfVQGbkMpz2CMncK48CbKez/EAHxFx+nYRaX3EKV9D+IlOlr3A61S4Aecf/M8/Vp/w8919XtJWEpco5Vtx2s0ODi4quOiEqCGgF0Lvt4JDC/4Og3cBjxjmiZAP/CkaZqfaGiiRB1OZgBtK6ea+x5tR/+G9Kmn8WIZsrd8jOL+h/DjbQDI1QLazBCBrGJ37wdpvodUvOFhCHza/QL2pRNo05fQpy6QOf5dMse/S2XwNgo3fJBq302RTSKRZAk7LzYyFIStICoB6jBw0DTNfcBl4DPAZ+eetCwrC9T6uKZpPgP8RrODE4SJEsmzL0DgL7p5bwWSXaTr5a8Tu/IuhQMfYuaOT4GyeDNE30iFAWbZF5HxO3ZSUjKwNxyFVQoTJM+9SPLcS/QMH2Xmzn9A4dCPNPJH2RCRai4IW0Mk7rCWZbnALwNPASfCh6xjpml+yTTNT7S2dYs5mQFkz0EpTra6KWui5kbo+97/hjF+mql7P8vM3T+9JDitl5fqJnf7Jxj5+Jeo9N5I+sTfR7pmoW/7IlFCELaAqPSgsCzr28C3r3rsi8sc+2gz2lSP2zabyZcbxUv1tKoZa9b+5l8iORXGH/1V7K59jTmJopG79cfoffq3SZ75IYUbP9yY82xQEAR4VQ81Hpm3vyAIdUSiB7WVOJlwcn1LpZoHPvrUBcq77m5ccJpld++n0nuI9MnvgxfRuZ4AnLLT6lYIgrACEaDWKNDiuImOpYkSQet3bJXsUt21Smp+DNmtYnfsbko7crd8DKWaJ3X2xaacb60kRcLJiwAlCFEnAtQ6OJkB1AU9KLmcZeD/+x/o+cFvY4y915I2xYaPMvjkvyZ55rklz+nTFwGwO5sToOyeA1S7D5A6+X3wohcIZF3GzkW0dycIQo0IUOvgtg2g5a+AH2aDtb/9V8hOGaU0Rc+z/wfdz/4e+uT5prXHGD1B10t/iBT4xIePLXlem7qIr+i46b6mtSl3yxOo5RmS5xtelWrNZEXGLbutboYgCCsQAWodnMwAku+iFieIjRwjcekNcjd/lNGPfZGZO38SLTtMz9O/jTZzueFtMcZO0f3Cf8LJ9FPcfS/6xJla4JyjT1/Ead8JstLw9syp9h6i2rWP9Lt/t6Q9UeBWRYAShKgTAWodnNlMPn3qAu1vWDjpPvI3fhgUjcKhx7jy0X8DkkLy7A8b2g594gxdP/wD3FQ3Ex/6r6kM3oHs2ejTC6pG+R7azFDThvdqJIn8TR9BLU0THz7S3HOvgl9t/ZyhIAjXJgLUOrjpfgIk2t7+a9TSFNPv+8yiNUW+kaK0624SFw7X6thttvjQW/Q8+xW8RDvjj/xyuMC25wAAxvj8PJiWG0X2nKYlSCxUGbgFN9FJ8kxjA/V6eHb0enWCICwmAtQ6BKqOl+xCqRYo7nsQezYwLFS84WFkt0ri0uubfPKA1Mnv0/nSH2J37GT8sV/Dj4WFa/1YGifdjz5+una4Npsg4TS7BwUgyRT3f4DY2KmwGnqEeI4nFusKQsSJALVOdscuPCPFzB31dgUBu3Mvdttg2HvYrG0qfI/2Nyzaj3yT8s67GH/kV/CN9KJDqr0HMBbMQ+lTF/HVGG6LFhUX9z1IIMlheago8UXJI0GIOhGg1mn6fZ9m7PH/jkBP1j9Akije8DD6zBDa9IVNOWfi0uukzv6Q/I0fZuqBn6tbqqjacwDZraLNDAFhgoTdsatldQP9WIbyjjtJnH85UinnQRCIxbqCEHEiQK1ToCdX3F6itPtefEUndWZzeg/6+Bl8LU729k8sG3CqPQeB2Xkoz0Gbudya4b0Fijc8hGKXSAy92dJ2LCRrMnZWrIUShCgTAaqBAi1Oac/7iV96PazysEH61Hnszr3X7A35sQxOuhdj/DRadgQp8LA79mz43BtR7TmEk+qNVLKEoilisa4gRJwIUA1W3P8QsueQuHB4Q68juVW07Ah218rBptpzEGP8DPrUeQDszl3X/oZGkySKNzyEMXmuKWvDVkOSJbyKmIMShCgTAarBnI5d2B27SZ/6AZJTXvfr6FMXkQjCHtQKwnmoCslzL+HpSbxE17rPu1lKe+4nkNVIJUuIxbqCEG0iQDXBzF0/hVKapuMNa90ZffrUOYBVBqhwHkqfGcLp2B2J3W19I0ml/2aMsVOtbkqNyOIThGgTAWoNfHd9G93Z3fvJ3fpjJC6+RuLC+mrT6ZMXcFI9+MYyWYML+PE2nNm08qZXkLgGp20QtTAOfjR6LmKxriBEmwhQa+CUnHUXGc3f/BEqPQdpf+MbqPkri5/0vWv3rIJgPkFileZ6Ua2oILEcJ9OPFPiRWbQbuAG+J0oeCUJUiS1FV8l3fZL9Scrj65xHkmSm7v88fX/3ZTpf+iPyN/0oxvhpjPH30PJXCJAIVINA1SkP3sHM+z5d+1alPI1SyWF37V316co77yQx9GbDNyhcC3dus8fcCG7bYItbA77v41U95IT4nCYIUST+MlfJczx67txYNQY/3sb0fT+Lnr1M1ytfJ3HxMG6yi9zNT5C/+SMU9z2Ak+4jdfaHi3pZ+mS40HdNPaj+Wxj+1G/hx9IrH9wkTrqPAAktN9rqptSIbTcEIbpED2oVAj8g0ZtAS2lI8sYSDioDtzH+yK/iq3rdLTDkSo6Bb32R1OnnmLn7p4EwQSKQNZz21vc6NkTRcFPdqNloBChZlanOVIl3xVvdFEEQ6hA9qFVwSy5dt3UhSRKytvFLVu09iNO5p+7+TH4sQ2n3PSTOv1JLS9cnL2B37AR563+ecDMDaLmRlQ9sAkVTsPNisa4gRJUIUKtgdBoYaQMARW/8pn+FA48gu9VwN1rfQ5++tKb5pyhz2gYik8knyRJeWWTyCUJUiQC1Artg03lzZ+1rWW/8JXM691Dt2kfqvWfRZoaQfGdN809RFrVMPpFqLgjRJQLUCuLdcZJ982uPmtGDAigcfAS1OEHm2LcBtk0Paj6TLxrzUCJACUJ0iQC1gsGHFycmKLrSlI3uyjvuwou1ER89jhfL4MWvXTl9q5jP5Fs8DxW7fIT+b30Rya02tT2imoQgRJcIUCvQU/qir7WMhu82YXGnrFA48EFgNr08AuWKNsUymXzJcy+ilqebPvTnu35zfp+CIKyZCFBrpCf1pt3Qivs/gK/FqfTd2JTzNUuYyTcfoCSnTOzKSQDU4kRzGxOIXpQgRNXWz1tuMi2pQZM+cPtGmpGPf4lA1Vc+eAtxMv3ERt4JM/lkldjIMaTZrD6lMNncxgRhCSstuXR3YkEQWkv0oNZIMRRo4mhboMVatl17ozhtA4sy+eJDb4fzbHqi6T0oWZfFxoWCEFHb687XBLIuI22X+aAWWZjJJ7k2sdHjlHfcgZfsRi02twclqzJOwWnqOQVBWB0xxLdGkiQ1ZS3UdrYwky+QFWTPprzjLmS7hD51saltkWRJpJoLQkSJO+06NGst1La1IJMvPvQWnp6g2nMAN9mFUpoKtx9pIt8RWXyCEEWiB7UOsi7jVcSn7o1wM/3oM0PI1QLlnXeBrOAmu5ECH6U8g5ds3jb1niN+l4IQRaIHtQ6KJnpQG+VkBlCLE8huhfLOOwHwUmFQanaiROA0fuG1IAhrJwLUOiiGQnCtHXCFFTltAwD4aoxKb7jOy012A81PNRc9KEGIJhGg1kFLN6maRJN4todX9XBKTtO2QJ/L5KsM3ApKuAbJi7cTSHLTe1BiDkoQoikyc1CmaT4B/C6gAF+zLOvLVz3/68AvAi4wDvy8ZVkXmt5QwvJHvuNv+aE+3/HxPZ+2/W103tSJZ3tMHJmgOFJEjav4jo/neqiGihpTa4HMlm1810dW1//5xsn0Ux68g8LBR+cflBW8RGfTU82DINjwzyMIwuaLRIAyTVMBvgL8KDAEHDZN80nLso4vOOxN4F7Lskqmaf4S8FvAp5vf2uZWk2gUp+yQ3pGm586e2iaMsioz8MAAdtFm6sQUsY4YyYEkWmJxlYX2VDtHnjyCW3bXH6RllcmH/umSh91U89dC4YfBWgQoQYiWSAQo4D7gtGVZZwFM0/xz4JNALUBZlvX0guNfBj7X1BYuIOtyU6tJbDbP8Uj2Jum7t6/u83pSp//e/mW/X42p7Hp0F8MvDlOZrqAam/c2cpNdxIfe3rTXW40gCPBsDzUelT8HQRAgOnNQO4BLC74emn1sOb8AfKehLboGxVC2bDWJwA+QVZm+++oHp9WSZInBhwZJDiRxK5u3O66b7EaxC7Xt7ufok+fBa0zFB0mWcMut3+FXEITFovKRsd7dvm6anGmanwPuBR5Z5vkvAF8AsCyL7u7uNTVEVdVVfc9U99SWHBJySg4HPn4APbn+ArQLr1HPEz0Mvz7MzNkZtPjGC66qPbsAaKOKnwmDqJwdJfWD/53K+03sOz6+4XNczTVcUnqKru7NXXu12vfS9Uxco5Vdz9coKgFqCNi14OudwPDVB5mm+Tjwb4BHLMuqu7OdZVlfBb46+2UwMbG2jLDu7m5W8z2FcmGZEBpdTtGh/4F+cuUclFc+fjlXXyN9jw7TMHlhcsNBSpMTJIDK2AUqWjsAqVMvhk+ef5Pc3g8u+R41G75V3LbBJc+thu/5SJclgs7N/YWu9r10PRPXaGXb8RoNDq7ubzUqAeowcNA0zX3AZeAzwGcXHmCa5t3AHwBPWJbV3F3t6lB0ZcvtI2R0GKQGUg157d67esGH/FAeNbb+t5U7W0FCLcz/QcYvHwXAmDyH5NqLtx8JArpe/BqBojH2kX+9rnPKioxbFUN8ghA1kRijsizLBX4ZeAo4ET5kHTNN80umaX5i9rD/FUgB3zBN8y3TNJ9sUXOBrVmP7+psvM3We08vqR2pDQXuQE/ga4laJp9czaNPnqPatQ/Jd9EnTi86Xs2NoBXG0bPDKKXpdZ/Xt7d4WqYgbENR6UFhWda3gW9f9dgXF/z78aY36hoUQyHIBbVkicAPcIoOejqamwt6jofe1vi29d7Ty7mRcxt6DTfVVVusGxs+hkRA9o5P0fPs7xEbfZdq/y21Y+OXj9T+bYyeoLT/A+s6p1isKwjRc80AZZrmn7CKmRbLsv7xprVoi9BSGv7I/GJdt+zSe08vE0cnUHQFSY5Wlp9ne8S74g0/jyRJ6Gl9Q70oN9mFNhPOK8WHj+LGO7C79lHtvoHYlXfJLjg2fvltql37UErTxEaPrz9AbaPKIIKwXaw0xHcaODP7Xxb4FGGlh6HZ7/0kMNPIBkaVntZrn7qDIMDoNGjb18bux3e3bI8h3/VxSvVTsSUktFRztjWPdcY2dMP3kt2opSkkt4px5V0qg7eBJFHpuwktN4JcDkOUUpxEnxmivOMOKv03E7tyct1bdYghPkGInmv2oCzL+ndz/zZN8yng45ZlPb/gsYeBf9u45kWXmlAJ/LBzOdd7gnCeZ/fjuxl5ZYTKZKWpc1W+5y+7maKkShtKXliL1M4U0+9No6fWN6ToJruQfJfEhVfDzQwHbweg2n8THP0bYlfepbT3/trwXnnHnWgzl0mdewl98hx2z4E1n9Nzt1bCiyBcD9aSJPEAYQWHhV4BHty85mwdCxfrakmNRE+i9pwkS/Te3dvUeQ237NL3vr5l1zdtxhql1TLajA2tEXNT4ZqP1HvP4Ksxqr0HAXDaBvGMFMaVk0A4/+RkBvBSPVT7biSQZGKjx5d93WsJvKBphXIFQVidtdxF3gT+vWmacYDZ//8vwFuNaFjUzQUot+zSfrB9yfNqTEVLNyco+K5PojdBajCFntZrPbtF7Uk0Lx9Gkjc2nDiXaq7lx6j03wyyOvfCVPtuInblXeRKHn3iDOUd4V5SgRbH7tpPbPTE+k4aiGE+QYiatQSonwMeArKmaV4hnJN6GLjuEiQgTAaQdRlZlcnsztQ9Jt4db8qn8sAPanX1En2JJaWHAj9oep25WEesbqBcDS/RSTBbXKS8445Fz1X6bkSp5kmfeAqJYNHzlYGbw116y1nWLBCZfIIQNasOUJZlnbcs6wPAAeATwAHLsj5gWdb5RjUu6oIgILM/s2zGXmZvpuE13pySQ+89vbWK5LHO2JK8S6/qhY83UXJHEqe8ztp5soKX6CCQZCoLUsoBKn03AZA68zxuohOnfef8c7PHrqsXJbGpNQUFQdi4NU8UWJZ1EXgVGDJNUzZNMxKLfVsh0ZOg42DHss8bbUZDExPciktmT4bU4Hx1CDWu1oLVHN/ziXU0N0DFO+PI8vrfGnbnHioDtxLoiUWP+/F2nEw/UuCHvacFRXudth14scy65qEkVcIpNKYYrSAI67Pqu6dpmoOEezZ9CLh60mXrlVXYBP33Lb8lBYTDgEaHgZ2zN736uVt1SfYn6b27d8lzakIlcOe7UbIioyabO8Qnq7PnXOeo2dQDPwdB/SHCMN18tDb/VCNJVPpvIX757TDdXF7921JWZFHRXBAiZi0fcf8AsIEPAwXgHuBJ4J81oF3bRmb35g/zzQ3Z9b1/mf2cUjrBgpu7rMvISvM7ukbGWNSONZHkZQNM4cCHyN38Uezu/Uueq/TfjOyU0afOr+l0siqvf0hSEISGWMtd6wOE26y/BQSWZb1NuC/Tf9uQlm0Tib7Epvae5koWDT44uOzrxrvjixYKt2ojvuRgsiG9Ei/VQ+62Hw+D2FUqfTfhqzHa3v4m+Gs7t8jiE4RoWUuA8oC5v/gZ0zR7gCLX3ljwuierMka7sWmv57t+GJyuUUop3hNflJHWzDVQCyV6Ek3fkiTQE0zf+1mMqfO0HfmbNX2vKHckCNGylgD1CvBjs/9+CvgL4K+A1za7UdtNoj+B52xOpQLVWJoEcTUtqdWG9HzXR8u0JkAphtKS3lsnE9ezAAAgAElEQVR5193kDzxC+r1niA2tfpnewnk7QRBaby0B6h8Bz87++9eAHwDvcNW+TcJSmd2ZTRs+Ws0CWEmWakkRXsUj3t34IrHLaUYF9Xqyd3wSu2M3nYf/FKUwvqrvaUX9REEQlreqj7emaSrA7zK/lXoZ+J8b2K5tRY2rm5JF57v+qocLtaSGk3dAYt018TZDojfB5MQkqtHknpSiMfngz9P3979F10t/yNjj/6LunNVCYohPEKJlVT0oy7I84COsO2lYaL+hfcMJA17VI9mfXNWxsY4YvuMjKRJKrHWrAJIDyZbd+L1kFzN3/xT6zBDG2KkVj/c9f93VLwRB2HxrGeL7beDfmabZmgmNLa5tfxtqXF1/2jWAFC7+XY14bxy36qLG1U1fg7UWWkJr6e7DpZ134asxEhdfX/lgXwzzCUKUrGXc5VeAfuDXTdMcJ8zPkghTznc3onHbiSRJ9N7by+VnL6Ml1xfj1djKCRJzjIwRzkU1sUjscvS0jldp0Y1f0SnvvJP40FtM32OCsvy1D4IgzH5sbtENQRCWsZa71+ca1orrRLwzTnIwSWWysq7tKNZSIVxWZRRDaVmK+ULxrji5C7kNbcGxEaXd95I8/wqxkWNUdt617HESEm7FRU+3bs5OEIR5awlQH17m8appmnuB71qWdWXjTdreeu/u5fx3zi+6Wc/NTSnG8lvF+56/6uG9OYqhYHRu3hqs9UoOJpk+Nd2yAFXtOYhnpElcfP3aAUqVcIoO9DSxcYIgLGstd4xDwL8EHiOsaP7Y7Nd3A78EnDVN84lNb+E2o+gKHTd1YOdsvKqHmlDpvrOb3vf1kuhLoMbVulW13YpLoj9R5xWXl+hLNL1IbD1Gm4GktG4eDFmhtOse4iPvIDnl5Q9TZdySqMe3mdyKSzVXbXUzhC1qLT0oGfiMZVl/PfeAaZqfBD5rWdYDpml+Hvgy8N1NbuO203Gog3hXnFhnbFGPKb0zDUDuQo6xN8fQEvPDc5IkEWtfW7DpvLGztYFh1twGhq1cCFvefS/p088Sv/w2pb0P1D1GUiSckqjHtxLf9bFzNkaHsWwCTnmyzNTxKcqTZQI/QItpGJ0GmT0ZkMEpODgFB6/HQ+6Xr1kZZTPa67t+3Z0FfM9n+uQ0gR8gqzKSKqFndGLtsUXJPYEfUM1VqUxW8GwP3/EJ3KB2vJ7R0eIaSkype00CP5zfDPxw5+bAC8+n6MqieWXf82tVYGRVbtmoQ1SsJUB9FPiZqx77FvAns//+L8Dvb0ajtjtJkq65eDa9K83EOxOLHltNBYmrRenNHWuPURorNfRGdC125x7cZDeJi68vH6AkicARaeYrGX11lOJoEVmTibXHiPeG72U7b+PbPnbBxi26qAl10YesarbKyMsjAMja7M13BkpHSww8NLCh9XqBHyz73po5M8PUiSl2fGgH8c75vzvP9hh6bgiv4oUf5IL5QIIU/s2pCRXP8XCLbhhUNBlJkcJzSbPHnw+DiiRJ4bIOQ0GJKyiaglf18CoeXtWrZfAGBBBQa68kSUiqROAG4TFzb0EpnBcd7xgnX8jXvpY0iURvuNXP3PV1Ky75S3mKw8VaIAy88PUkOWzXXHWZuSCJD8jUArMkS/PBdcGlnHts7ueTFInM7gzxnsYXAFhLgDpDOJS3MAj9s9nHAboJa/MJGyTJEqnBFMWRYi3INGv7+EZJDibJXcitO4NxwySJ0u73kT7xd8iVHH6s/i7InivSzK8lez5LebxcCyZu2WXm9Ex441TDG5yEVPf3LElLH1cNlaAccPF7F+m+vZv2G67eySfc6dizvbAnltIW9VDmekDTJ6fZ+8TeuqW1ipeLaAmNy89dpveeXjK7Mzglh6FnhgCWLINY+LVXDd8Py5Xsqt34rxrc8Ks+ftWvHbOqbNpl/jRUQ0WpLmhjEP5MuTM5tEx4PeycXUuMWtg2aS7SBIsXokuSVNskKfDCYLYWgR9ELkD9IvBXpmn+S+AyYZFYD/jJ2edvBP7t5jbv+tV5Uye582Hmm+/5Wz6zLN4VX/SprBVKu+8lc+IpEpfeoHDw0brHiB5UyCk5qDF1Ua/EKTlMvDWx5GaraBtb5ybJElpCY/KdSabenar1FJDCm2rghb2KgADVUDE6DNr2teEUnfB4P0wImjk9Q/ft3Yte26t6VHNVtISGltAYe2OMykSF/OU8irZ8UlLUyaqMnJLBD69Lyz74NdiqA5RlWW+YpnkQeAAYBEaAlyzLcmaffw54riGtvA6p8fAP0bd9vIpHYmBtCRJRI2uzGxi28P7vZvqx2waJDR9dNkBtVlHfrSoIAqaOTzF1cgrVUOm5q4fUjhRBEDDy8khDq5LU66XU28fMztmMvDyCrCzuMRQuF5YEqNzF3KLX0BIahcsFFL3+XJEQLWtaxTkbjJ5vUFuEq3Qe6mTk1XDMfq0JElFkpMPdhVvJ6dhNbOTYss9fz/X4nJLDyEsjuCW3NoQ3eniU2OkYRruBW3QXBYRWkSRp0dzWHLfsUpmqEOuc/1spXCosaXMUfgZhdaIziy4skRhIhFtWGGpLywVtlkR/om4KfTM5mX6Uah7JLtV9fm6C+Xrjll3OP3Ue3/EX3cC1hIZbdsmey0b+xq7GVaZPTde+nhveE7YuEaAiTJIk0jvTyPHt8WtK9iXXPBm72dxMHwBabrT+AcH8xPj1pDRWQpbrp3vPzRFFnSRLlMZLtQ8YVw/vCVuP+O1FXMehDtr3Lc1s2orUuNryT+FOuj9sS75+0ZOAIKwmcZ0pT5RbWvV+swR+QGGoAEBhaOnwnrC1iAAVcYqhkNlbPyV6K2p1tpGX7CSQNbTcSN3nFV2hMlVpcqtaz87b2yJpQItrZM9l8WyPalYM7211rS91LVxXFENp7RCaJOOke1Fz9XtQiqZclzc2t+RGamH3RlSmK+HarAhUUdlOJKeCmh9Dy48i2YOEm1s0lghQQlOpcZXqTLWl60/cTD/65Lnln9/gxpJbjVcNKx1slwAlazITRyYWZfNdNwK/7s7Rkl1CqWQJFA1fixOoMSTfRc2Phf8VxpDtEpJrI3k2smcjuQ6SV0VyHeRqHrU8U3u98i2PAg81/McRAUpoKj2jk7uYa/4W8As4mX7il95AcqsE6tJq7y3bu6pFqrlqWH5nm1A0BSmzTXtPQQDLDMW2HfkmqZM/wNcT+LE0npFGdsqoxUnkaxRJBgiQCLQYgaIRKHoYyFQj/H88jtM2QDHdh5Ppw033oezcQaoRP99VRIASmkpv08NaZy3cBcTN9CERoObHcDp2LX2+7IY1zLbBnMxqlK6UWvqBoRG2bG/Q94iNnkDy5tcLarJH++gZ9Jkh1OwIhYOPkrv9JxZ9m1KaJvXes9g9N+Ck+5GreZRqHi+Wwe7ah5vswou3I/kOslNBciogSbjpXpx0H26qG5TVV6tRlOa8X7bXu1KIPDXW2i3oYUEmX260boAK/ACv4i1bf227sbP21r2hbwFytUj63b8HfHw9ia8ncdN9VHtuWDQcp2ZH6Hz1T9BnLi15DV9LYLfvwGnfQfrk9yjteT9uZn4OKH3yexD4TL3/c3jJrmb8WE1xffwFCpGhxtSW1+Rz0z0EkoyWv0K9gY+AALtgXz8BqmjPFxUVNlfg0/nyH2GMv0cgq8gLekZuspvCDQ9T2nMfiQuv0PbO3+KrMSbv/zxO+47accn2brKeCpKEXMnT/50v0fb2XzP5wV8CQC5nSZ59keLe+7dVcIIIBajZzQ5/l7DG7tcsy/ryVc8bwB8D7wMmgU9blnW+2e0UNkaSJWS9xZ/WZRU31Y26zGLduVTzRM/Wrn+4GoEf4JWvn97ipggCYsNHSJ57GclzwjkhScJJ9ZK77eME2nyV78w73yI2dpKpez9Lad+D4NnIdglj/DSpM8/TfuSbtB/5JgDlwTuYft9n8GPpxadLZSCXA8CPpcnd8gTtR75JbOQYlYFbw95Z4JO/+SPNuwZNEol3pWmaCvAV4EeBIeCwaZpPWpZ1fMFhvwBMW5Z1wDTNzwD/Afh081srbJSqqy0vJ+Sm+5etJqFoSstrBjaLU3DCvYGElQUBsZF3yBz7NvrMEG6iAy/eDkGAFPikrjxPfOQYkw/+E5yO3cSH3iLz7t9T2P9QGJwAFB0/rlPefS/l3feiZodJXDiM076T8q57lk2AWKhw8BGSZ1+g7e2/wmkbIHX2RUp77sNLdq/4vVtNJAIUcB9w2rKsswCmaf458ElgYYD6JPCbs//+S+D3TdOULMvaPulH1wlZl1ueKedk+oiNvAO+B/LSagPXS6p5aby04e0yrgueQ/cP/y9iY6dwk91Mvf9zlHbfu+i9o0+cofPlr9P7/f9I/qbHSb33DNXOvczc9VPLvqzbNkjujk+urS2ySvbOn6T7hT+g55nfg8Ajd/NH1/uTRVpUAtQOYOHM4BBw/3LHWJblmqaZBbqARVvPmqb5BeALs8fR3b22TxWqqq75e643G71G5d4ylZnWVmvQ+vYhvevTThk/M7jk+YBgw++DrfBeKp8s097d3rLEFUVRyGSiUSlFKs2gvfs0zk2PESQWlxeLvfQn6GOnqNz/s9i3/AiqrLKk1Zm7KQ8eJPbcfyJz4in8eBv2R/45mWTnhtpV9xqlH8Q9/yLq5aPYBx8mObh/Q+dYq3hbvCnv7agEqHp/HVf3jFZzDJZlfRX46tzzExMTVx9yTd3d3az1e643G71GRadIfjrf0kKemtZGHKgMn6YiL13R4VZcxsfHN3Tj3grvpfHL4wRu6wYhMpkMudn5lUbRZoZw2nZcc/hMH3+Prpf+CKWaR3nvBSYe+RW8RAcA8UtvkDn+PfIHHyO7+wEo1K+EPyd7/88TH3gNp30HrqfW5o/Wa7lrVLr9U7QFMHPwcbwGX8OrVZUqsYn1L4QeHFz6obCeqOSWDgEL8313AsPLHWOapgq0AVNNaZ2wqfTM7FqoFnLT165qHrhBy4chm8Etbe+hTG3qAn1//x+IX3q9/gFBQOrd79Hz7O/ja3Gm3v85lEqenqd/B6UwgZofo+O1/4dq1z6yqx2KkyTKe96P27a6m/B6ueleJh/6QtMz93zXJwia86EmKj2ow8BB0zT3EW4n/xngs1cd8yTweeAl4B8CPxDzT1uTltZa+qkdIFAN3EQH2jJVzYFtn2rull18198We40tR5+6AEDy3MuUd9+75PmOw/+F5IVXKe28i+l7P0ughVUTup/7Cr1P/w6+HieQFaYe+Cd15yqvF4Ef4JQdtIRGelea9oPN2WEhEj0oy7Jc4JeBp4AT4UPWMdM0v2Sa5idmD/vPQJdpmqeBXwf+VWtaK2yUFtdavhYKwpp8y6aaG9u/qnl1ZvsXxdWy4UCMMXYKpbR4wEWbvkTywqvkD32YqQd+vpYe7nTsZvzRfw6Bj5q7wtT9n68N921XXtXDKTq4JRfPCUcOZFVGTajoGZ1Ef4Ldj+1m70f30n17d7iesQki8/HQsqxvA9++6rEvLvh3BfjpZrdL2HyKobS0WOwcJ91Pcvx03QKbsiZjZ7d3qnnxSrFpN5pW0bLDuMlu1OIEiQuHyS/Idku99zS+apC7+SNL5qfctkHGHv8N1MIE1d5DzW72urkVFySWlK5yyuEeZ5IsEbgBAQGKquC7Pnpap/1QO5k9mfBvM0Ilvrb3u1OIJEmWkLXWd97dTB+y56CUpuquIWn19vSN5uSdSHxQaJjAR8sOU9z7AHr2Monzr5C/KQxGcjlL4uLrFG74IIFef0G2l+jES2wsA6/Z1LhK2/42SmMl7IKNb/vE2mN03tJJsj+JrMq4FRen4FDNVkn0J9CTq6/B12wiQAktoRjKku3fK9kKekJvWvByZmuZabkr9QPUNl0LFfgB0+9NU52pbp8dZ+tU+VaKU8huFadtEKdjJ52H/xR98hx2935Sp5+DIKBw8NHWtHcDAj+o+8HC93ySg0nab2in/Ybl54jUmIoaU4l3x5c9Jipa/zFWuC7VuzHKslwb/26GuQC13DyUV/Galq3UDIEfMHVyinPfOcf0yemGBye5nCVx/hXkcrah58kc+zZ9T/37cKh2AS17GQCnfZDyzrvxFZ3E+VeQXJvk2R9S3nEHXira69Tm+I6PU3JQYsqylT+8ikdmVzTWlG0W0YMSWkLRlUUpzkEQYHQY2DPNm/cJ9CSekV42k2+7VTW//MPLVLPVcH5CW99ryOUsHa/9GXbPAQr7H1o6PBYE6BNnSJ15nvjQW0iBj69oFA88Qu6mxwn0JBDuzqrNDCGp+9nI52Rt+iLp499FIkDLDuO075x/LjtMgISbGSBQDco77yJx6Q3cTB+KXaJw6LF1n7eRvKoHEkiKhKzIyKpMZm+G9gPtKIbC+NvjFC4XllSglzUZo6OF+9g0wPb4yxO2HC2pUZ4o1/7IvIpH+4F2JnLNXdjqZvqW7UERbJ9U8+zZLNXp6sZ+liCg442/IDZ6gvjocdLHv0tx34NUBm9Dy46gT11AnzqPWpzE1+IUDjxCeccdJM++SOrk90mefYHKwK2o2WG07AgSAV6mj9xjv77sPNA1+R4dh/8MX0+g2EWMsVOLA9TMMG6qu7YpZWnvAyQvvErbkSexO3ZjdzWv+kLgBwResOLwte/6GG0Ggx8crCUrXL3gO7M/w8zZGXR18dyRntEjleCwGbb+X56wJc0t1p0LUIEfkNqRYvLoZFPb4aT7SVx6vf4chqFQmdz6Vc29qsfE0YkNB9r4pTeIDx9l5o5PUe27kdSpH5A68zzp088C4MbbcTp3k7vpI5R3v68WGOyeA+Rvepy2d76FceUkTvtO8jvuxItlaH/zL+l85f9m8uH/qu5W5XO0mSF8NbZoSC598vvo2ctMfOAXaTv6JMaVkxQO/cj892QvhxUkZlV7bsBNdKKWpsgfemxVhVk3i1t2aT/Ujp21cYoOTtFB0ZVFvaAgCAj8gIEHB64ZaIy0ES7VWMBzPNJ70st8x9YlApTQEnpaXzSWrsbDidtmb8XhZvqQnTJyNY8fWzx+L2vytqhqPvrK6JquqzF6nPjIMXI3f7R2TeRKnvY3v0G1c284NCbJTN/3j8ne/gm07AhO+44l128ht22QyYe+sOTxWMwg/uIfkzn2bXK3/fiS55XSFG1vf5PE0JsEkkxx/0Pkbv0Ysl0ic/w7lHbcRWXHncSunCRx4dVa8V/JraIWJijtef/8i0ky+UM/QvLCq5R33r3q67FRQRAQ64rRfet8cPVsj9GXR6nMVGqp/m7ZZecjO1eVJBTviVMeL9eSJTzbI71LBChB2BRKbPF6C70tHK5QdKWpW3HMZ/KNUq1zg3VKTtPa0gj5S3kq05VV956SZ1+g/fW/QCIgfvF1Zu75aco776H9zW8gu1Wm3//ZRT0dP95ONb7+qgLOTT+CP/IemRNPYbfvpLLzLgh8lHKWxIVXSZ94CoDsLR9DqeZJnn2BxIXDeLEMgaIxc0+4NLLae4jUmefRpy5gd+9HzYVDiM5V5YaKBx+hePCRdbd3PbyyR89dPYseU3SFHR/awdS7U0y9O4UkS3Tf3k2sY3X17dr2t5G/mEdLhj0pzdDQEuucWIwwEaCElli4WNeturTd0AY0fyuOhZl89RZkVmequBV3Sy5oreaqjL85vrrgFARkjn2bzInvUum/meytH6fjzW/Q9fLXqXY9izF5juxtP4GbGdjcRkoS03f/9Ox253+M//Zfo5RnkGYz8ko77yJ7xz/Am60IXjjwCG1HnyQ+fJSp9/9srddW7TlIgIQxdgq7ez/aTFhBYuEQXyMEQYDv+viOT+AH6Kmla4qUmEKit/4wcedNnST6EhSGCtdMDb9arCNWW4wbBAF6R3TXMm3E1vurE7YFSZJqNeACJyA1GFYUVzSlqQHKj7Xhq7FlM/kUXWHq5BS9d/Y2rU0b4Xs+ufM5cudyVPPVJXMVdQUBHa/9GcnzL1Pc9yDT93waZIWxx/4b0qe+T+bYd7Dbd5G/8cONabSiMfmBX6D97b8mkORwgWyyE7t9J07nnkWHupl+Jh/6AnK1iG8ka4/7RhKnfQfG2CnytzyBlh3GV41aYGsEp+SQ3pVGT+toKY3AD7jy2pVFPRnP9mi7oe2ac0qxjtiqe06Lvq87RmWqglf1SN+8/Yb3QAQooYUUIyy1osSU2qd8xVAIckHzspEk6ZqZfLIqU7hcoOeOnshnSJWnygy/MAx+OKe32goB+uQ5kudfJnfj4+Ru/8R88oCskL/pIxT33BcmPDSwWKofbw8Lsq72+AXBaU6191C4ANez0bOXw+G9ayRerFcQBLhllx0f2kG8c/FiVztrkz2bra0xC7yAjgONqeOX2Z+hcLkAEiT6t3Yiz3LEQl2hZWQjfPsZbfNrN7SMhu82dysOJ9237LYbAL7tUxguNLFFa1caK3H5ucuohrrmbD1j/DQAhRs/XDezzY+31wqpRlm19xCS72JMnEWbGV4y/7QZ5oLTzkd2LglOAJ23dKKlNXwvHPKL98YbVhkl3hVH1mW0lLZtK9KLACW0jGIoeLZHrHt+eENPNn+vKDfTj1LJITnlus9rCY2ZUzNNbdNyCsMFLn7/Itlz2VoWZGG4wPCLw+ueJNcnzuCk+/GNpRs3biXVngMEkkziwqvITqkhAcqreOx6bNeyQ3KSJDH44CC+6+OWXbpubtxeTZIkEetc3/DgViGG+ISW0ZIabtkltSO16DGavJfhwkw+u2tf3WOqM1Xsot3Swppe1ePKa1dQYyoTRyeYODpBrCNGeaJcy+Zas8DHmDxHqYlp140SqAZ2514SF98AwGnf3AQJt+LSfUf3oh5/PYqhMHDfAONHxjHaG1vZoeNgx7buZmzjH02IOiNjQMCim6tiKE1/V85n8i2/eaESU5h+d7rhbQmCgNHDo+QuLd3Ce+SVkdpw0VzBT7fsrj84AVp2BNkpY/fcsO7XiJJq7yGkIEyy2ewelCRLZPasrtZdoi/B7g/v3tTz1xPvjtcdatwuRIASWkZLaRgdxqLkA1mXm56M4CW7CGT1mvNQsiJTHC42dI1WEASMvjpKcbTI+OvjTJ2Y32Bvrvq4rGzun6w+cQaAanfzyv40UqUvXCrgJjrXNW/mlJ3a3kkL+Z5PciC5pu1JtvVWJk0iApTQMlpCo/PGxWnAktSCvaIkGSfdi5pfPkBBGEAmTzSuFNOV169QHiuHiQ4JlelT04y+Nko1V2Xy2GRD1mIZE2dw4+14icbNlTST3bmXQNbW1XvyPZ9Ed4J4V3xJxXC37NJ589baG2o7EHNQQsvImkx699L1G82uJgHgpvvQpy9d8xg1ppI9ncXJOfTd17epvZmxt8YoDi/e4VaNqxRHiosqBmyqIECfOIvdvb+pdekaStGYuv8f4SbnA27gBwRBsOLvy3d8+u7tIwgCzn/3PHJi/vhYZ2xbVmqIOhGghMhpdjUJCOeh4kNvgWeDsnwihBpXKU+VufS9Swx+cHBdN63KdIXc+RxO0cEtu3hVj8AP6vaQVEOFBs2zK6Vp1PIM+W0yvDfn6jp7c+Wq6lV5qB1Tdui5vae2fqn9hnZy53MohoJTchi4bZMraAirIgKUEDnNriYBYaq5RICWH1u0ZUM9ihb28C5+7yJGu4Gsysi6jKyEgdWtuLgVl4nUBGW/jJ7SiXXFqExWKI2VwsSGhFabo2jVGhZjdv7J7t4eCRLLkSRpyd5JCwV+gJ7WadvfVnus69Yu8hfzQJi4kxjYngtho07MQQmRo8SaP8S3mky+hSRZQo2peBUPp+BQnapSulLCztvhNiKKjKIp+FWf8kSZiSMTlMfLyIqMntIjMYGuT5zBV2MNWS8UJZIsXXMjP7fq0n9//5Lv6b6zGztnk96VjnwVke1KBCghcuZW4l+tkUHLTfUQIF0zk28lywUdSQqDWRSC0kLGxJkwe68B5YCiRDEUYp2xuhVK5vYhq7e+Lb0zTbwvHq41Elpie78zhS2pXjUJp+TgFBq49YWi4aa6V8zk2y7kajFcmLzN5p/qUWIKyYFkuJX6Vdyye821TTs/uLM2LyU0nwhQQuTUrSYRgJps7JSpm+lHW+UQ31anT54FoLrN558gTDQx2oy6PVhJkYh3bd+FrludCFBC5CiGAlfdS9S4umKJmY1y0n2o+THw3YaeJwqMiTMEsop91XYW203gByjxcO8xLbU043K5wCVEgwhQQuTIurzkpqFn9HAeoYGFZO2eA0iBR/zykYado6UCH23qPJmjT5I49wp2xy5QtvfaHs/xah9s9DadIJifx/Qdn3iv6D1FmUgzFyLn6moSnu2R6c0Q74ozeWISXWtMwdZK/824yS5Sp5+nvOuehpyjVeIXDtN29EnU8gyBJFPtOUDutp9odbMaLnAC9Lbw/ZIaSFEYKtTWrrm2S2b36mrrCa0hApQQSQurSfiOT7I/iRpXkeUGdvolmcIND9N+5G9Qs8O42yH9OghIn/webUefpNq5l9ztP0F54FYCfemGf9uSRG1X4VhXDBYkgmqG1pgKHcKmEUN8QiTJ+vxbU9ZktKSGrMhr3oxvrUp7HySQNVKnn2/oeZoi8Gl76/+l7eiTlHbdw/ijv0ppz33XT3AiTP1XYmEW3lwFeAjrKuodrds6RVgdEaCESFK0+dRePa3XFkou94m3mq1uynl9I0lp9/tIXHh12Q0Mo06yixijx+l66T+TPv0s+YOPMXX/57f1fFNpolR3nZyiK4sW2WqZ2eG9skt619I6kEK0iCE+IZKUmEKQDQj8gOTA/Cd+o8Ogmq0uKl3jOR6yJuO7/jVL2qxW4cAHSZ5/meT5VygcfHTDr9cMkl0kc+IpYsPvoBXGAQgkmZk7PhVu5b6N+a5PvDuOW3KX1Eac6z3NiXXGyGVzIEGy7/rpSW5VIkAJkaSlNfxhH8/2SA7O30gSfQmm351GTi0IUFWPwbyJg5UAAA/aSURBVIcHGX5hGF3d+LCN07Gbaudekqefp3DgQ9GptOC7xIfewo9lsDv3EKgG+B7Jsy+QOfa3yHaZysCtlPY+gN21B7tj97r2RNpqvKrHzgd3cunppdXor15kmxxIMv3uNHpGb/62LsKaiQAlRNJcNQlZlhdtm21kjCUD00bGINmbRDM2bwireOBDdL76xxhjp6j23bRpr7sSbeo8+uQFyrvvwTfmh6DU/Bidr3y9tiVIIMk47TuRPBstN0ql9xAzd/1UUxI7Aj8g8IIVb/C+54dFf+X5RIVGMNoNYh2xJdXg61WINzLheynWFWtYe4TNIwKUEElz1SS0Nm3RPj6yJodbUMzyXb9WqsboDIf/NqOwZ2nnXbS9/Ve0v/mXTN33OZzOvRt+zWuRy1najv4NyQuHAWg7+iTF/Q+Rv/HDxMZO0v6GBZLC5AM/R6DG0CfOYkychcBn4gO/SGXwjqbt6eRVPFI7U+Qv5VET6pLr7ZZdkMO07o6bO7CzNqOHR1GNpfUIF65LWg+n5NB3Tx8QvmcWljPyHb+WYj5HksMlDCK9fGsQAUqIpLlqElffYCAc/nNLYbUHr+qR2RfebDJ7Moy8PLI5qcOKxtT9n6fj8J/S+/3/SPHAB8ne9uMbHjKT3CpqbhSlkkfyqkiug1qcIHXqaaTAJ3fTj1LecRep08+QOv0sqdPPIgU+1e79TN3/ebxEuKtrZeDWjf+M6+BVPTpu7qDzxk7aD7Zz5fAVnKJDEARh1fB2g7Yb2sjszdQ+WOhJnT0/uofhF4Zxqy4SEr7rY2QMAi+cZ1xvNQdZlUntTIXnyeiUrpRqr+U5Xt0q5j1391yzurkQHSJACZEk6zJu2SU1kFrynJ7RsfN2mHaeUmsT4/HeOJKyeb2Iat9NXPnof0/mnb8ldfo54kNvUe05iBdvx4tn8GIZAkUnUDQCVYfAR3YqSG4V2amgywGZwjSyU0Yp51BzI6iFCSSW9hrK/397dxsj13XXcfx7587Tzs7OPngfvE/2Oq5xYpnGSYMVBWidpJXitmoAhT+gQs1DqRBIPJQKAn3Fi4ogISBCqKhNQSlCag6lVSNVFKmhESCkChwnDSUqaYtjOzbYTuxsbK935s69vLjjfbbXa+/OvTP395FWnjt7Z/b47Ln7m3vuuedMvpOL7/xxmtVhAC4c/Aiz+95P9bv/TLOnxqU9D0KufZOWRlG06swoiiL8ks/gD8Sze5f6S0w/PM3sa7Pky3kqo5XrBk2hUmDHwzs4//J5/LJP/0w/fsmnWqjyontx3Q8VURQRXA6WTVcUhRGV7Ys/szJaYfa12YX24Hneml2LfZMavdcpEg8oMxsCngFmgOOAOecurNjnAPBpoAY0gU85555pb0mlna7NJtEzvPqMpbK9wsXvXSRfytM7uTiAIufnKPWX1py1+lZFhR7euucxruz8Ifq//VWKb76GP/ctvPDmZlYveTnCQoWwVKUxMMmVHffR6J+gWRkkyheJ/CJhobzmvUnN6jBvHfiJTfu/3KwwCPFyHo2rjWXdcsFcwNR7ppYFl+d59M/0X++tlvFyHiN3jyx7rtxfpneyl6vnr153BGYURTTnm4wdHOPcsXN4eY+cn6Mx12B87+JKtytvxPU8TzORd7jEAwp4HHjOOfeEmT3e2v7dFftcAT7inHvVzCaAo2b2j865i+0urLRP/67+Nf/AlAfKeJ5HMB8wsHtg2fd6J+JRWpv9h6kxtJPzP/qr8UYU4TXm8OffxgvqeM0GXrMOeISFHqJCiTBfprptlNnLV9t2bWgzhM14qP70Q9M0601O/9tpgrmAXC7uSisPbv7ggtEDoxz/h+PXDajmXJPJQ5OUB8pURiq8/q+v05hrUOorUepb7KrLl/PLfu+50uo5HaWzpCGgHgUOtR4/DTzPioByzv33ksenzewsMAIooLrY6L2jaz7vF33yxTzkoVhdfo2qb7qPN/7zja395Ox5RMUKQXGdZcDzJfA25wbidrh2o+vUoal4Gft8jh0P7+DcS+e4dPISowfW/n3cLr/oM7BnYOGseKnGlQYTD0xQHoiD0S/5TD84zdljZ9cMy6UDJVa+l3SeNNwIMOacOwPQ+veGR4GZHQSKwPfaUDZJ0I0+/foVn8rI6oDIl/MU+rp3xoStEoURYRAy/eA0fnEx3D3PY/TAKLs+sGtTboK+nqE7hxbmWYyiiMaVBvXLdbbft53K6PLfs5fzGHvXGP13rO5aLPYVF4J25U260nna8hHDzL4ObF/jW5/c4PuMA38DHHHOrbnugpl9DPgYgHOO4eHhDZU1n89v+DVZk4Y6au5p0j/TT7l/9afo+d3zXDp9adnw9CT4vk+t1hnDmeffnmf34d1tX7xvaVsqPFjg7Etn6RnqobazRnWsumzKq5tR2Ffg1IVTFCoFekd7E2+nmyENx1tSvNu9D+F2mdl3gEPOuTOtAHreObd3jf1qxN1/f+ic+7ubfPvo9OnTGyrP8PAw58+f39BrsibtdXT1wlVOfuPkqu6/dqvVaszOziZahpsVNSNmDs+0/edudlsK5gKOf+04uUKOobuGVl2j7ERpP95uxcTEBKxalnS1NHTxPQscaT0+Anxl5Q5mVgS+DHx+A+EkGVUeLFMdr9JsbN5ovm4WhRE9I90xJVK+J0+umKNZb275Csyy9dIQUE8A7zOzV4H3tbYxs/vM7KnWPga8G/h5M3ux9XUgmeJKJxg7OIaf99ec4VqWC+aCNa/ndKpitQge5CsaJNHpEu/i22Lq4tsCnVJHjSsNTnz9xKr52NqlU7r4wiBk5vDMpkwRtVFb0ZbOvnCW2ROz7P7Q7q4YZt4px9tGdFIXn8iWKFQKjN8/TuPKzd1Um0XXuveSCKetUhmrLNxsLJ1NASVdrTJaYXDvIMF8kHRRUqnbuvcgXvNp5bpQ0pkUUNL1BvcMJl2E1PJLPuWh7lp6It+TZ+AdnT96TxRQkgG5fI7qRJWwueatc5kVhRHloXJXde9dc21CW+lsCijJhG37tsWL58mCYC6gf3d3de9Jd1FASSbke/KUt5Vve4G8bpIrrj1bvEhaKKAkM4buHIpXe82YKIxoXG4QzAVEYUSumKPYV2R4/3BXdu9J99CdbJIZldEK+Z5km3wURQtrFq21/HkYhISNcGE9rJudoDUMwni1Ws/DwwO/tbx5MUd5W5m+6T56hnu2dMJXkc2mgJJMGXjHAG98+4223LwbhRH1y3X8QjxSzi/58RlLKyPCICQKooXZLq6d2RT7ikRBRP3tOo0rjYXA8vzF10aNiGajSRTEq9yWBkv0jvdSrBUTnyRXZLMooCRTajM13nzlzdt+n2AuXn682L9iQtownng1iiIqoxV69vZ07Ug5ka2mgJJMyflxl1d9tn5LodGsN8GDkQMj1HbeeCmNbpyiRqSd1BcgmdM33XdLgyWCuYDqZJWZR2bWDScRuX06g5LM6R3r3fBrgvmA6lSV0Xu2ZtlzEVlNZ1CSOblCPBjhZoVBSLFWZPRehZNIOymgJJN6hnsIg/WnPorCCDyYfGBSAx1E2kwBJZlUm6kRXF3/OlRzvsnUe6bIFXSoiLSbjjrJpGKtuO69UM35JoN7B7V0g0hCFFCSSZ7nURoo3XCfiEjLNogkSAElmVWdrF53uHkYhPRN9WlqIJEE6eiTzOod773u7ObNepOhu4baXCIRWUoBJZnlF/01h5tHYURlrNKW+fpE5PoUUJJp5eHyqpV2g7mAbfu3JVQiEblGASWZVttZoz5bX7gnKooiSoMlSn03HkAhIltPASWZVh4oM/PIDNXJKp7vUZ+ts22fzp5E0kCd7JJ5pf4SI3ePAPHoPY3cE0kHHYkiSyicRNJDR6OIiKSSAkpERFJJASUiIqmkgBIRkVRSQImISCopoEREJJUUUCIikkoKKBERSSXvessNdImu/s+JiHQwb70duv0Mytvol5kdvZXXZelLdaR6Uh2pjjbha13dHlAiItKhFFAiIpJKCqjVPpN0ATqA6ujmqJ7WpzpaX2brqNsHSYiISIfSGZSIiKSSAkpERFJJK+ouYWaPAE8CPvCUc+6JhIuUODObBj4PbAdC4DPOuSfNbAh4BpgBjgPmnLuQVDnTwMx84D+A151zHzSzXcAXgCHgBeDnnHP1JMuYJDMbAJ4C9hPfo/iLwHdQO1rGzH4L+ChxHb0M/AIwTgbbks6gWlp/XP4COAzsA37GzPYlW6pUCIDfds7dBdwP/FqrXh4HnnPO7QGea21n3W8AryzZ/iPgT1t1dAH4pURKlR5PAl9zzt0J3E1cV2pHS5jZJPDrwH3Ouf3EH5Z/moy2JQXUooPAd51z3299MvkC8GjCZUqcc+6Mc+6F1uO3if+oTBLXzdOt3Z4GfiyZEqaDmU0BHyA+Q8DMPOAh4IutXTJdR2ZWA94NfA7AOVd3zl1E7WgteaDHzPJABThDRtuSAmrRJHByyfap1nPSYmYzwD3AN4Ex59wZiEMMGE2waGnwZ8DvEHeDAmwDLjrngtZ21tvTHcA54K/N7JiZPWVmvagdLeOcex34Y+AEcTC9BRwlo21JAbVorak3NAa/xcyqwN8Dv+mcm026PGliZh8Ezjrnji55Wu1puTxwL/Bp59w9wGUy3p23FjMbJD6r3AVMAL3Elx1WykRbUkAtOgVML9meAk4nVJZUMbMCcTj9rXPuS62n/8/MxlvfHwfOJlW+FPhh4ENmdpy4a/gh4jOqgVY3Dag9nQJOOee+2dr+InFgqR0t917gf5xz55xzDeBLwANktC0poBb9O7DHzHaZWZH4wuSzCZcpca1rKZ8DXnHO/cmSbz0LHGk9PgJ8pd1lSwvn3O8556acczPE7eafnHMfBr4BPNbaLet19L/ASTPb23rqYeC/UDta6QRwv5lVWsfetXrKZFvSTBJLmNn7iT/5+sBfOec+lXCREmdmPwL8C/Fw12vXV36f+DqUA3YQH1Q/6Zx7M5FCpoiZHQI+0RpmfgeLQ4OPAT/rnJtPsnxJMrMDxINIisD3iYdP51A7WsbM/gD4KeIRtMeIh5xPksG2pIASEZFUUhefiIikkgJKRERSSQElIiKppIASEZFUUkCJiEgqKaBERCSVFFAiIpJKCigREUklLVgokmJmtpt4Gq73OudeMLMJ4FvAY8655xMtnMgW00wSIilnZr8MfBx4F/Bl4GXn3CeSLZXI1lMXn0jKOec+C7xKPP/hOPDJZEsk0h4KKJHO8FlgP/DnWZgkVATUxSeSeq3FIl8iXnLhMPCDWZ/xW7JBZ1Ai6fckcNQ591Hgq8BfJlwekbZQQImkmJk9CjwC/ErrqY8D95rZh5MrlUh7qItPRERSSWdQIiKSSgooERFJJQWUiIikkgJKRERSSQElIiKppIASEZFUUkCJiEgqKaBERCSV/h96sz1c2cpQ9gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, np.mean(grad_all, axis=1))\n",
    "ax.fill_between(x, np.mean(grad_all, axis=1) +np.std(grad_all, axis=1), np.mean(grad_all, axis=1) -np.std(grad_all, axis=1),\n",
    "               alpha=0.3, color=\"purple\")\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('grad')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 元にしたMLP本の図（P.96, 図6.3）をほぼ再現できた。\n",
    "  - x=30,60辺りで変動が起こっているらしいことは分かる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}