narrowlyapplicable/WAIC&WBIC with pystan2.ipynb

## WAIC&WBIC with pystan2.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# WAIC&WBIC with pystan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 人工データに対し、pystanを用いてWAICおよびWBICを算出した。\n",
    "  - データは <https://github.com/narrowlyapplicable/peak_separation_whth_WBIC> で作成した信号ピーク模擬データを使用した。  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. 必要なライブラリのインポート"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- pythonのバージョンは3.6.6\n",
    "- 使用ライブラリ群は以下\n",
    "  - numpy 1.15.1\n",
    "  - pandas 0.23.4\n",
    "  - pystan 2.17.1.0\n",
    "  - matplotlib 2.2.2\n",
    "  - seaborn 0.9.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as sct\n",
    "from pystan import StanModel\n",
    "import pickle\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "plt.style.use('ggplot')\n",
    "sns.set_palette('deep')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. データの作成"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Lorenz関数(Cauchy分布の確率密度関数)型のピーク信号が、複数個重畳して測定される状況を考える。\n",
    "- scipy.statsを用いて3本のLorenz関数を重ね、さらにガウスノイズを付加したデータを作成した。\n",
    "  - 2つのピークが重畳 & 1本のショルダーピークが存在する、という状況を想定している。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "x = np.arange(-10, 10.5, 0.5)\n",
    "k_true = 3\n",
    "peak1 = sct.cauchy.pdf(x = x, loc = 0, scale = 1.0)\n",
    "# peak2 = 0.16 * sct.cauchy.pdf(x = x, loc = 1, scale = 0.4)\n",
    "peak3 = 0.09 * sct.cauchy.pdf(x = x, loc = -0.7, scale = 0.3)\n",
    "peak4 = 0.36 * sct.cauchy.pdf(x = x, loc = 3.0, scale = 0.6)\n",
    "\n",
    "data = pd.DataFrame(data={\"x\":x, \"y\":peak1+peak3+peak4+np.random.normal(size=x.shape[0], scale=0.01)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4HNWZ+PtvL9p3qW3J2mxZlizb2BgwNluIAwbMEjvkFx+MkwkzmBByw8wdyMydZCaZzJDcXEISCDMhi0MyQxZiDpCAQwCzh0AwNjYY77Jky5bUsqzN2rfurvtHd4u2rKUl9VItvZ/n4UHddarqdXV1v1WnzmIxDAMhhBDCbKzRDkAIIYQYiSQoIYQQpiQJSgghhClJghJCCGFKkqCEEEKYkiQoIYQQpmQPppBSai3wMGADHtVa3z9s+V3AlwE30AXcqbU+qJSaBxwCjviK7tBa3xWi2IUQQkxj4yYopZQNeAS4BqgDdimltmmtDwYUe1xr/VNf+XXAg8Ba37JqrfXy0IYthBBiugumim8lUKW1Pqa1HgC2AusDC2itOwJepgDS+1cIIcSUBFPFVwDUBryuA1YNL6SU+jJwLxAPXBWwqEQp9T7QAXxda/2XEda9E7jT93KL1npLcOELIYSYroJJUJYR3jvnDklr/QjwiFJqE/B14DagASjWWrcopS4CnlFKLRl2x4UvIfmTkgH8bAL/BiGEELFnpNxylmASVB1QFPC6EHCOUX4r8BMArXU/0O/7e7dSqhooB94ba4dO51ibD47D4aC5uXnK24mEWIk1VuIEiTVcJNbwmGmx5ufnB1UumGdQu4AypVSJUioe2AhsCyyglCoLeHkjcNT3/ixfIwuUUvOBMuBYUJEJIYSY0ca9g9Jau5RSdwPb8TYz/6XW+oBS6j7gPa31NuBupdQaYBBow1u9B3AlcJ9SyoW3CfpdWuvWcPxDhBBCTC8WE063YUgVnznFSpwgsYaLxBoeMy1WXxVfSJ5BCSGEMBnDMOjr68Pj8WCxjPtbHzKNjY309/ePW84wDKxWK4mJiZOOTxKUEELEoL6+PuLi4rDbI/szbrfbsdlsQZV1uVz09fWRlJQ0qX3JWHxCCBGDPB5PxJPTRNntdjwez6TXlwQlRBhVO7t44V0n1c6uaIcipplIVutNxVTiNHf6FSKGVTu7eOjJwwy6DGxWJ5+9Zi4XL8whPk6uC4UIhiQoIcKksrYDl9vbStbtMfjV9hp+/VINszMT+fSVhSxfkMXhk+2c2N3CgjkJlOanRjliIcxFEpQQYVJelI7d5sTtNrBaLdx4ST4ut0F9cw8piXaqnV3899NHcXkM4uwW7tlQIUlKhFW1s4vK2g7Ki9Jj4lyTBCVEmHT1usjLTmLJvAyWlWae84PwwrtO3B7vHZbLZVBZ2xETPxrCnH7wxOFz3rtoYTarl89mYNDNA787TF1zD4YBFks9hY5krrowl8vOc9DVM8jP/lh91rpfuaVi3H0+8MADZGdnc8cddwBw//33M2vWLDZv3hySf5NUhgsRJvXNPdSe7uHGS+aMmHjKi9Kx270PkC0W72shwqW334V/XAbD8L6eqltvvZUnn3wS8LYq3LZtGzfffPOUt+snd1BChElrxwBpSXbi40buM1Kan8o9Gyp48s/1OJu7mJeXEuEIxXQy1h1PfJyN228s5aEnD+N2G9hsFm6/sXTowik1OS6oO6bhioqKyMrKYv/+/TQ1NbFkyRKys7Mn/W8YThKUEGHS2jFAVnr8mGVK81P5zCdK+d5v36e6vlPuokTY+C+IQv0M6tZbb0VrzenTp9m4cWNItuknVXxChElrRz856Qnjlju/zMGG1UXkZk+ut70QwSrNT+X6VfkhfdZ5/fXX8/rrr7N3715Wr14dsu2C3EEJETa52YlBVdslJdhZc1FeBCISIvTi4+O57LLLyMjICHoIpGBJghIiTL60vmz8Qj59A24+qGpjQUEajozx77qEMAuPx8OePXv42c9CPxG6VPEJYQK9/W7+54XjvHdEpksTsaOyspLLL7+cK664gvnz54d8+3IHJUQYfFh9hideP8k/fLqc3OzEcctnpcVTnJvM3qo21q6cE4EIhZi68vJy3nnnnbBtX+6ghAiDpvY+mtv7SUkMvk7+/NJMjjd009E9GMbIxHRhwslmRzSVOCVBCREGLR0DxNutpCQFX0lxfmkWBrDv2JnwBSamDavViss19c624eRyubBaJ59mpIpPiDBo7RggOz1+QlMNFM5KIistnpOne7g8jLGJ6SExMZG+vj76+/sjOvVGQkLChGfUnSxJUEKEQWtHP9njdNIdzmKx8I3PLyElUb6WYnwWi2XSM9VOhcPhoLm5OSL7km+CEGFQVphGVtrEEhQgyUmIAEF9G5RSa4GHARvwqNb6/mHL7wK+DLiBLuBOrfVB37KvAZt9y/5Ba709dOELYU4bVhdPet1fbT9OUoJtStsQYjoY9+mVUsoGPAJcDywGblVKLR5W7HGt9VKt9XLgAeBB37qLgY3AEmAt8GPf9oSYtjweY0otl/oG3Ow83IonRlppCREuwTSvWAlUaa2Paa0HgK3A+sACWuuOgJcpgP+btR7YqrXu11ofB6p82xNi2jp8soO7H97N8YauSa2/rDSTju5BTpzqDnFkQsSWYKr4CoDagNd1wKrhhZRSXwbuBeKBqwLW3TFs3YIR1r0TuBNAa43D4Qgm9jHZ7faQbCcSYiXWWIkTohtr/7FeXG6DuYW5OLLGf4g9PNaPr0jnse01HG0Y4OKl5jrecg6Eh8Q6yr6CKDNS+8Vz6h601o8AjyilNgFfB26bwLpbgC3+5aFoIRLJliZTFSuxxkqcEN1YTza0YLWAMdhFc/P4d0EjxbogP5Ud+5xce2Ho5tYJBTkHwmOmxZqfnx9UuWASVB1QFPC6EHCOUX4r8JNJritEzGvtGCAzNR6bdfJ9Uy5f6sDZ3IvbY0xpO0LEsmAS1C6gTClVAtTjbfSwKbCAUqpMa33U9/JGwP/3NuBxpdSDQD5QBuwMReBCmJW/k+5UXLI4Nqp7hAincRtJaK1dwN3AduCQ9y19QCl1n1Jqna/Y3UqpA0qpD/A+h7rNt+4BQAMHgReBL2ut3WH4dwhhGheUZbFqcc6Ut+NyezjRKA0lxMxlMeGAg4bTOfVawJlWpxsJsRInTI9Yn3mrju27TvGDLy0n2SQdeKfDcTWjmRar7xnUuHXXMlisECE06PLQ0TMYkpGml5Zk4PEYHKhpD0FkQsQeSVBChNCJxm7++ScfcKCmY/zC4yiZk0pakp291TK6uZiZJEEJEUKtHQMA5EyxkQSA1Wph6fxM9h9vx+32THl7QsQaSVBChJA/QU21FZ/f+aWZ9Pa7+e0rJ6h2Tm5kCiFilSQoIUKopbOflEQ7CXGhGXIyKcGK3WbhrweaeejJw5KkxIwiCUqIEGrtGAhJ9Z7fsYZu3B4DwwC326CydurPtoSIFeZouyrENHHF0lkMukL3vKi8KB2b1YnLbWCzWSgvSg/ZtoUwO7mDEiKELijLYuWiqXfS9SvNT+USX6ffv7+5jNL81JBtWwizkwQlRIgMurwjP/QNhHawlEVzMwBITowL6XaFMDtJUEKEyKnWPr7zm4McDHHH2tysRABOt/WFdLtCmJ0kKCFCpKWjH4Ds9ISQbnd2lnd7pyRBiRlGEpQQIRLKTrqBEuJsZKXF09LeH9LtCmF20opPiBBp7egnzm4lNSn0X6tv3raExPjQ9K0SIlZIghIiRFo6B8hOi8diCf0Eg0kJ8lUVM4+c9UKEyLUr8ujqdYVl2zWnunhtz2k2rC4iLVla84mZQZ5BCREiJXNSWTo/Myzb7u5z8+6hFhpapaGEmDkkQQkRAi63hw+OtnGmayAs2/c3NW+UBCVmEElQQoRAS/sAP9lWxaET4RkrLzstHrvNQmNrb1i2L4QZSYISIgRaO71NwEPdxNzParUwOytR+kKJGUUSlBAh0DI0D1RoO+kGKpyVFLZtC2FG0opPiBBo7RjAYoGs1PC1sNt8Q2nYti2EGQWVoJRSa4GHARvwqNb6/mHL7wXuAFxAE3C71vqEb5kb2OcrelJrvS5EsQthGq2d/WSmxGGzSaWEEKEyboJSStmAR4BrgDpgl1Jqm9b6YECx94EVWusepdSXgAeAW3zLerXWy0MctxCmctOlBVy5bHZY99Hc3s9jLx7n+lVzWDwvI6z7EsIMgrmDWglUaa2PASiltgLrgaEEpbV+PaD8DuBzoQxSCLNzZCTgyAjf8yeApAQblXWdnDc/QxKUmBGCSVAFQG3A6zpg1RjlNwMvBLxOVEq9h7f6736t9TPDV1BK3QncCaC1xuFwBBHW2Ox2e0i2EwmxEmusxAmRjdXjMXhxx0mWlmZTlJs24fWDjdUBpKfE09FL1D4HOQfCQ2IdZV9BlBlpYDFjpIJKqc8BK4CPB7xdrLV2KqXmA68ppfZprasD19NabwG2+Lfd3NwcRFhjczgchGI7kRArscZKnBDZWNs6B/jFHw+yac1ckmwTr+abSKyzM+OpcZ6J2ucg50B4zLRY8/PzgyoXzBPdOqAo4HUh4BxeSCm1Bvg3YJ3WemheAK210/f/Y8AbwAVBRSZEjGj1zwOVFp4+UIFysxJplL5QYoYI5g5qF1CmlCoB6oGNwKbAAkqpC4CfAWu11qcD3s8CerTW/UopB3A53gYUQkwbrZ3+eaDC+wwKYH5+Kme6BnG5PdilxaCY5sY9w7XWLuBuYDtwyPuWPqCUuk8p5W8y/j0gFXhSKfWBUmqb7/1FwHtKqb3A63ifQR1EiGnko0664b+DumLpLP7h/5RLchIzQlD9oLTWzwPPD3vv3wP+XjPKen8Flk4lQCHMrrWjn+REm0woKESIyUgSQkzRp68s4poVeRHZl8cw+M//3c9F5dmsu7wgIvsUIlokQQkxRYnxkbt7sloseAyDBhnVXMwAUpEtxBQ981YdR+s6I7a/3KxETktLPjEDSIISYgp6+ly88G4DNae6I7ZPb1PzfjzGiN0RhZg2JEEJMQWRbMHnl5udyKDLQ1tneGbvFcIsJEEJMQX+PlCR6KTrV5KXyhVLY2NYHCGmQhpJCDEF/lEkItFJ169odjJ/c21JxPYnRLTIHZQQU9DePYjdZiEtObLXeh7DoLffFdF9ChFpcgclxBR86opC1q6cg8Uy0pjK4fP9rYeJj7Pyj59ZGNH9jqfa2UVlbQflRemU5qdGOxwR4yRBCTFF0RhBIicjnur6rojvdyzVzi4e1IdxewzsNif3bKiQJCWmRKr4hJiC375cw57K1ojvNzcrkdaOAQYGPRHf92gO1bTjchsYBrjdBpW1HdEOScQ4SVBCTNKgy8ObHzbhbIn8qA65WYkYwOkz5umwOycncehvm81CeVF6FKMR04EkKCEmqS2C02wMl5edBEBjq3kS1EULc1hYlEZKok2q90RISIISYpKi0UnXb3ZWAjddms+cnKSI73sspQWp9PS7KZ6dHO1QxDQgjSSEmKTWzsjNpDtcQpyNT15mrtHMv//EYRLirGxaMxcZhUmEgtxBCTFJDS19xNksURtyqLvXRe3pnqjse7j+QTdV9Z3MzU3hymWziY+TnxYxdXIWCTEJ1c4u3vigEZfH4L9+X0m1M/JNvp99u54f6MMYJrhdqTvdg2HA3Nxk6pt6OCXTgYgQkAQlxCRU1nZEvUl1bnYCvf1uOnuiP6LEiUbvndzcvBQefrqS7TtPRTkiMR1IghJiEvxNqC1Er0l1bpa3gcQpE8wNdaKxm/SUODJT48nNTpQ7KBES0khCiEkonp2MYUBFcRrrLi+MSpPqvGxvv6PG1j7KC9Mivv9AxbnJQ60Z87IS2R2Fzsti+gkqQSml1gIPAzbgUa31/cOW3wvcAbiAJuB2rfUJ37LbgK/7in5ba/1YiGIXImpa2r0t+C5d4ohaf5/stHjsNguNJriDuvrCvKG/c7MT6e5z09UzSGpyXBSjErFu3Co+pZQNeAS4HlgM3KqUWjys2PvACq31MuAp4AHfutnAN4FVwErgm0qprNCFL0R0NPkS1KzMxHFKho/VamHzDfO5bEl054bqG3AzMOgeep2b5T0mZqh6FLEtmGdQK4EqrfUxrfUAsBVYH1hAa/261trf3nUHUOj7+zrgZa11q9a6DXgZWBua0IWInqYzvgSVEflRJAJdWJ5NviO6nXXf2tfE//3fe+jqGQRg/pxU7r65jHyTdSIWsSeYBFUA1Aa8rvO9N5rNwAuTXFeImBBnt1KcmxzxeaCGa+scYNfhFtzu6A0ae7Kxh/SUuKHqvJQkO0vnZ5KcKI+4xdQEcwaNNNHNiB0vlFKfA1YAH5/IukqpO4E7AbTWOBxTr7Kw2+0h2U4kxEqssRInhD/Wm69ycPNVi0KyranEuu9EHY/+6Rj/de/HyHWE/1nYSLHWNR+krCjrrPcP1bTS2TPIysW5YY9pNHK+hkckYw0mQdUBRQGvCwHn8EJKqTXAvwEf11r3B6y7eti6bwxfV2u9Bdjie2k0NzcHEdbYHA4HodhOJMRKrLESJ8ycWFPivH2gDh87RYIlM5RhjWh4rH0DbpxN3VxUlnnW+0+/VoWzuZf5syM/V5bfTDkHIi0Usebn5wdVLpgEtQsoU0qVAPXARmBTYAGl1AXAz4C1WuvTAYu2A98JaBhxLfC1oCITwqQMw+Cb/7Ofqy7MZfXy2VGNxd8gobG1F0rDn6CGO3m6BwOYm3f24LC5WYnsrT6D2+3BZpPulmJyxj1ztNYu4G68yeaQ9y19QCl1n1Jqna/Y94BU4Eml1AdKqW2+dVuBb+FNcruA+3zvCRGzOroHTdG0G7zPe1KT7FGLx5Eez4bVRczLO7t6MS87EY/HoLk9OuMUiukhqKeYWuvngeeHvffvAX+vGWPdXwK/nGyAQpiNvwWfI8ot+PxysxI5FaV5obLTE1hzUd457wc2Nc/Njl5TfBHb5N5biAka6gNlkgT18eWzKc1PjcqAtQdq2mnvHjznfX9SkiGPxFRIghJigprO9GMBcjIiPw/UcNXOLn790nFeeu8UDz15OKJJqrffzX89Xcnb+5rOWZaSaOebt53HJ5ZHrxWfiH2SoISYIEdGAhcvysZugof/gaOquyI8qvrJ090AFOemjLg835Ek80KJKZGzR4gJuuw8B5tvKI12GIB3VHW7zdvd0GIhoqOqnzjlTVBzc0ee3r2yrpNn3qqLWDxi+pEEJcQEeTzRnyDQrzQ/lXs2VJCREkdBTnJEB6492dhDTno8aaMMCFvT0MUL7zbQ3Rv9+apEbJIEJcQE9A24ufvh3bz+fmO0QxlSmp/KykU5NLT2njVoa7idaOwetXoPIDfbPPNVidgkCUqICWhp78ftMUhNMtc4cxXFabjcBlX1kWskcffN5ay7fPShNfOyPpqvSojJMNe3TAiTM8M0GyNZUJBGvN061EcrEsbr3+TIiMdqtdDYJk3NxeRIghJiAswyzcZwifE2HvzyBcTZI1Mpsv94Oy0d/Vy5bBYWy0hjQoPNZmVWRgJnus7tJyVEMCRBCTEBze39JCfYSDFZFR8QseQE8Nf9TdSc6ubj5489FuHX/2Yx8XHRGzBWxDbzfcuEMLGywrSozwE1mtaOfn7x/DHWrpzD0vnhHTj2RGMPc/NGbyDhJ8lJTIU0khBiAlYszOamS80552ZachwnGrs5eCK8nXW7e100t/czd4wWfH4nG7vZ8scqWjoi92xMTB+SoIQIkscwaOscwGOYpx9UoDi7lQX5aRw5Gd4E5R9BYrQOuoEGXB52V7bhbJaGEmLiJEEJEaQznQN8dcte3hph7DmzWFicRn1zLx094WuY0NjmHYtwrD5Qfv6m5tEabV3ENnNWpgthQmZtwReoojgdqOfIyQ4ursgJyz5WL5/NqkU5JCWM/3wpNTmOlMTozVclYpvcQQkRJLP2gQpUnJvC8gWZpCSG99ozmOTkl5cdvfmqRGyTBCVEkJrO9GO1WshKi/40G6OxWS18aX0Zi+dlhGX7Hd0D/OgPlVTVdwa9TnFuMnH2kftKCTEWqeITIkjN7f3kpMdjs5r/x7ar14XdZiExPrTNvKvr29l3rH3EWXRHs/GquSGNQcwccgclRJAuXeLgpkvzox3GuBpb+/jKj99nT2VbyLd9rN7bQrB49vgt+ISYKklQQgTpvJIMLlnsiHYY45qVlUBakp0jYZi8sLq+ndmZCSRP4BlXR/cg3338IHsqW0Mej5jeJEEJEYRBl4dqZxd9A5GbzmKyrBYL5UVpHD7ZgRHiPlvVde1BddANlJxoo+ZUN7VNPSGNRUx/QV0GKaXWAg8DNuBRrfX9w5ZfCfwQWAZs1Fo/FbDMDezzvTyptV4XisCFiKSGll4e+N0hvriulAvLsqMdzrgqitPZXdlGY1sfeb55maZq0OVhVlYSpQUTmxTRbrPiyEyQlnxiwsZNUEopG/AIcA1QB+xSSm3TWh8MKHYS+Fvgn0bYRK/WenkIYhUiaoaamGeYt4l5IG9/KDhysjNkCerk6R4uXDiLguyJN7zIzUqUeaHEhAVzB7USqNJaHwNQSm0F1gNDCUprXeNb5glDjEJEXbOvk67DxJ10A83KTODz181jkS9RTVW1s4uHnjyM221gs1m4Z0PFhKaXz8tK5PDJDjyGgXWU6TmEGC6YBFUA1Aa8rgNWTWAfiUqp9wAXcL/W+pnhBZRSdwJ3AmitcTim/iDabreHZDuRECuxxkqcEPpYO/tPkZ4SR1FBbsi26Reu47p+9ayQbeuV91sZdPmeZ7kN6lvdrFoWfMznV7g402OQkpZJSmJcyOIay0w+X8MpkrEGk6BGutyZyJPXYq21Uyk1H3hNKbVPa10dWEBrvQXY4t92c3PzBDY/MofDQSi2EwmxEmusxAmhj7XuVDvZafFh+feH67j29rt5/2gbZYWpUx79ovbUGQAsFrDZLBRk2yYUc1menbIb5tLb1U5vhGaln8nnaziFItb8/OC6awTTiq8OKAp4XQg4gw1Ea+30/f8Y8AZwQbDrCmEW668o4NNXFkY7jAnpG3Dz2PbjfFB1ZkrbaescYHdlG4vmprHp2vIJV+8FMutI8MKcgrmD2gWUKaVKgHpgI7ApmI0rpbKAHq11v1LKAVwOPDDZYIWIlpI5k/tBjqastHhysxI5UtvBNSuCH/lhuG1v12MYBp+7poSK0oJJXz1/61f7Kc1PY9MaGVlCBGfcOyittQu4G9gOHPK+pQ8ope5TSq0DUEpdrJSqAzYAP1NKHfCtvgh4Tym1F3gd7zOog+fuRQjz6up18d6RVjrDOIVFuCwsTqOythO3e3Ltl7r7XHxQ1cYnls+ecgOROLuVU20yL5QIXlD9oLTWzwPPD3vv3wP+3oW36m/4en8Flk4xRiGi6mRjNz9/rpqvqIWkJUfmAX+oVBSn8+beJmoaeyZVLZeSaOc/b1+KPQTjD+b6WvIJESwZSUKIccTCNBujWViYhgVvkp2ojp5BDMMgPTluQkMbjSY3O5EzXYMxMRqHMAdJUEKMo/lMP3abhYzU2Lp7Au+EgQ98aTmfuGBizeM9HoOHnzrCz5+rHr9wkPyz68rkhSJYkqCEGEdzez+OjISY7WCaPolqyXcPtVDX1MvyBVkhi6NodjIfP382CXHysyOCI2eKEONoOtMfMyNIjKS1o5+fbqsKepLBgUEPz75Vx9zcZFZUhG7cwVmZiWxaMzdkQy+J6U8mLBRiHHetX4Brkq3gzCApwc4HR9vo7nPxqSsKx20s8eqeU7R1DXL7DfNDftfodnvo7nOTnhJ71aUi8uQOSohxODISYvqq39nibdpdWdvJ97ce4v2jo8/L5DEM3jvSyrL5mZQXhWYcv0Df14f59q8PUO2M0HASIqZJghJiDM3t/bz83inau2OvD5RfZW3H0IBlHgO2/LGajlH+PVaLha9uWszfXDsv5HFUO7uoaeimvXuQh548LElKjEsSlBBjOObs4qk/19Ld64p2KJNWXpSO3WbBaoE4m4VrVuQNVbFt39lAZV0nhmHQ2TPIwKCHOLs1LFVwlbUdeHwjHbnchjdxCjEGeQYlxBia2/3TbMRHOZLJK81P5Z4NFVTWdlBelD70DKpvwM0rexrp+EsdJXNSGBj00NE9yF3rFrCgMC3kcXgTpROX2/DN+hv6KkQxvcgdlBBjaDrTT0ZKHPFxE5+kz0xK81O5flX+WQ0kEuNt/L+bl7Hp6rm0dgxQ39xLZ6+LHz59JCzVb95EuZA4m4Ul8zImPeCsmDkkQQkxBn8fqOkqPs7Kx5fPZvXy2UPz6rjDWP22oCCNz107j7UrJz94rZg5JEEJMYbm9n5mZU7fBOW3sDgdu937nMpmC2/12yWLHZQWhL4KUUw/8gxKiDHcd/tSBgZjtw9UsEZ7ThUObreHI7WdZKbFk58Tu833RfjJHZQQY4izW0lJmhnXcSM9pwoHjwGPPHOUt/c1hXU/IvZJghJiFDWnutGvn4zpPlBmFGe3Mi8vhap66QclxiYJSohRHHN28eqexmiHMS2VFaZxsrFbpt4QY5IEJcQomtv7SYizkp48M6r4ImlBQRoeA443THyeKjFzSIISYhT+JuaWGJ1mw8xK81OwANXO4EZYFzOTXBoKMYqmMzOjiXk0JCXY+cZtS2J6EF4RfpKghBiBYRgMuDySoMKowJEc7RCEyUkVnxAjsFgs3H7DfFKS7DLqdpi0dQ6w9bUT1Df3RDsUYVJB3UEppdYCDwM24FGt9f3Dll8J/BBYBmzUWj8VsOw24Ou+l9/WWj8WisCFCKfq+k4eeuoILreB3ebkng0VMnZciFkt8Pr7p8lJT5C7KTGice+glFI24BHgemAxcKtSavGwYieBvwUeH7ZuNvBNYBWwEvimUipr6mELEV6/eaWGQZeBYYR3bLqZLCM1nlmZCUFPRS9mnmCq+FYCVVrrY1rrAWArsD6wgNa6Rmv9ITB8TJjrgJe11q1a6zbgZWBtCOIWImya2/txNvdhtRKRselmsgUFaVTVd2EYRrRDESYUTBVfAVAb8LoO7x1RMEZat2B4IaXUncCdAFprHA5HkJsfnd0Dh6BnAAAgAElEQVRuD8l2IiFWYo2VOGFqsW7ffQSrBf5p0wXUN3WzpCSbhXPDd+M/U47rSC6o6OWdA830G0kUzgptFepMPq7hFMlYg0lQI3UCCfZyJ6h1tdZbgC3+5c3NzUFufnQOh4NQbCcSYiXWWIkTJh/roMvDy++eZOn8TEpz7ZTmZgDusP67Z8JxHU1euoWMlDiO1zaSaO0L2XZhZh/XcApFrPn5+UGVCyZB1QFFAa8LAWeQcdQBq4et+0aQ6woRce8fbaOz18Xq5bOjHcqMMDsrge9+8XzpDC1GFEyC2gWUKaVKgHpgI7ApyO1vB74T0DDiWuBrE45SiAgpmZPCTZfmUzFXnjlFgiQmMZZxG0lorV3A3XiTzSHvW/qAUuo+pdQ6AKXUxUqpOmAD8DOl1AHfuq3At/AmuV3Afb73hDClWZmJfPKyAqzywxkx+4+3868/30tb50C0QxEmE1Q/KK3188Dzw97794C/d+Gtvhtp3V8Cv5xCjEJExF8+bCI3K0Fa7EVYapKdlo4Bquo7ubgiJ9rhCBORkSSEAHr7XejXT7LjYEu0Q5lximYnkxBn5WidjNghziYJSgjgnQMtDLg80jgiCmxWC/PzU6XDrjiHJCgx4xmGwZ/3nqZkTgrFuSnRDmdGKitIw9ncS3efK9qhCBOR0czFjHektpNTrX387dqSaIcyYy2el86ZrgEGXcMHoxEzmSQoMeN19brIz0lixcLsaIcyY5XMSaVkjgzGK84mCUrMeCsWZnNReZb0yYkyj2HQ0t7PrMzEaIciTEKeQYkZrbG1D7fHkORkAtveruc//ne/VPOJIZKgxIzldnv4gT7Mr7Yfj3YoAijJS8HlNqg51R3tUIRJSIISM9YHVWdo7x6UZ08mUVqQBsBRaW4ufCRBiRnrjb2nyUmPZ8m8jGiHIvCOKJGfk0SVdNgVPpKgxIz07sEWKms7WVKSgdUqz5/MYkFBKtXOLjwemcBQSIISM1C1s4v/3X4MgHcONFPtlCt2s/jYslncceP8oCecE9ObJCgx41TWduCfYdztNqis7YhuQGJIcW4KyYl2XtrVIBcOQvpBiZmnvDANu82C221gs1lk9HITqXZ28aA+jMttEGd3cs+GCkrzpQPvTCUJSswogy4Pv3nlBNeuyCPObqW8KF1+AE2ksrYDl9t7ezvo8t7dyuczc0mCEjPKzsMtOJt72fDxIhZL6z3TKS9KJ87uZNDlTVIZKXFRjkhEkzyDEjOGYRi8/N4pCmclsUimdDel0vxU7tlQwQ2r5pCaZOPFnafoH3RHOywRJZKgxIxxoKadhpY+rlmRJ0MbmVhpfirrryjkCzct4HRbH8/vcEY7JBElUsUnZoyXdp0iMzWOi2XkiJhQUZzOF24qlarYGUwSlJgx1l9RSFfPIDabVBzEiot8FxODLg+9A27Sk+WZ1EwiCUrMGNIaLDYZhsFDTx7BarVw74aFMvLHDBJUglJKrQUeBmzAo1rr+4ctTwB+BVwEtAC3aK1rlFLzgEPAEV/RHVrru0IUuxBBaeno54V3G7jxknyy0uKjHY6YIIvFwseWzeJ/XzzOCzu9n6OYGcat61BK2YBHgOuBxcCtSqnFw4ptBtq01guAh4DvBiyr1lov9/0nyUlE3Ku7G3l7fzOGIQPoxKpLFuewsiKb5/5aLyNMzCDBVMavBKq01se01gPAVmD9sDLrgcd8fz8FXK2UkvtwEXXdfS7e2tfEyopsstMToh2OmCSLxcKmNXPJSovnF3+qprffFe2QRAQEU8VXANQGvK4DVo1WRmvtUkq1Azm+ZSVKqfeBDuDrWuu/DN+BUupO4E7f+jgcjgn9I0Zit9tDsp1IiJVYzRLnkRNtHDjeypKSbBbOzRqxjD/WN9+opn/Qw4Y1FTgc5uz7ZJbjGoxox/qVTQn84o+HqG8zcDa3B3UOxAKJdZR9BVFmpDuh4XUlo5VpAIq11i1KqYuAZ5RSS7TWZ43OqbXeAmzxr9fc3BxEWGNzOByEYjuRECuxmiHOQyfaefipSgwgzmbhHjXyWG0Oh4OGU6d57q3jLJ6bTmrcQNRjH40Zjmuwoh1rTgr8n4/N4aEn9uJyG9htllHH64t2rBMx02LNzw/uOWIwVXx1QFHA60JgeM+5oTJKKTuQAbRqrfu11i0AWuvdQDVQHlRkQozgaF3X0NXRoNvgxZ0No84dNODysHxBJtetnBO5AEXYVdZ14nIbGIaMRj/dBZOgdgFlSqkSpVQ8sBHYNqzMNuA239+fAV7TWhtKqVm+RhYopeYDZcCx0IQuZqIlJRnE2S1YLGCxwIfVZ/ju7w6NOBxOSqKdTWvmUVFszqo9MTnlRenYfE3NLRYZjX46GzdBaa1dwN3AdrxNxrXW+oBS6j6l1DpfsV8AOUqpKuBe4Ku+968EPlRK7cXbeOIurXVrqP8RYvozDIPHX6kB4J4NFay/vIB/vqWCzTfMp2ROCglxNgDcvrupqtozVNV3RitcEUal+ancqyqYnZmAzWYhN1Mav0xXFhM2vTWczqmPvTXT6nQjIZpxvnuohV8+f4zPXjOXK5fNHrFMfXMPj/zhKFcsncU7B1vp6RvggS8uN/3IEbHy+YO5YnW29PKtx/Zz5fmzufXquecsN1Os45lpsfqeQY3b0tvc31whgJ4+F0+9cZJ5eSlcsXTWqOU8HrBa4dm36znd1ktvv4eaxp4IRioiKT8niY8tm8Wbe09zqqU32uGIMJAEJUzv2bfr6exxsWnNXKxjjEJeNDuZS5d8lMAMQx6gT3efvKyASxY7SIi3RTsUEQaSoISp1TX18Oe9p/n48tnMzU0Zt3xFcTpxdgtWCzKd+wyQlhzHbWtLZAiraUoGixWmNic7kQ2ri7h0cXAdA/0T3tW3uinItskAsTOEs6WXN94/zcarimUw2WlEEpQwLcMwsNmsXH1h3oTWK81PZdWy2HnoLKauoaWXP+89TdHsZD62bPTnlCK2SBWfCEq1s4snX6uK2ECdXT2D3PerAxw60R6R/YnYdmFZFgsKUnn27Tr6BmSK+OlCEpQYV7Wzix88cZitLx/lQX04Iknq93+p41RLL+kpMkGdGJ/FYmHD6iI6e1y88G5DtMMRISIJyuSqnV288K4zqlMMvHekZagDrCsCQ8tUO7t4e38zV1+UR4EjOaz7EtPHvLxUVi3K4ZXdp2jp6I92OCIE5BmUiVU7u3joycO+QTGdow6KGW5zc1OwWMDfpzsvOzFs+zpa18nPn6smNcnGTZfKxHRiYj51RQHZ6fGkJMpP23DVzi4qazsoL0qPmcZD8ima2JGTHUODYrpc3juXaJxYlyx2MCszkcr6Po7UNFM4Kzx3NdXOLn741BFcbgOb1UJ9c2/MfJGEOWSnJ/CpKwqpdnaxs7JaWnL6VNV18uCTR/B4DOz26F3sTtS0TFDVzi7e3N8e8ydnQrx16K7FAIpmj98PKJR2HGzmaF0XG68q9rWMmxfSlnGGYVDX1Mu7B5txG5CebB+qSvR3so3lz09ER7Wziwe1t+bBaoXPrpnLZefNGrOT93R3+GTH0Hdr0GVwqKY9Jr5b0y5B+U9Ot9vAbh99rphYUN/US7zdwqVLHLy5t4m9VW2cV5IRkX2faunl8VdOUJybck6/kvrmHtq7Blk8b+KxVDu72FvVRu+Am6r6LpzNvVitFi4sy2JFeRZ2mwW325BOtmLSKms7cLu9P8YeD/z6pRP88a9Ozi/NYt1l+aQmz5yGN/2DbjweWDQvgxd3NTDo8h6Xt/c3c978DOblmfu3cdolqMpab7UYfPRAPxYT1KDLw56jbVxYns2mNfOw2628truRy5c6wn5SDQx62PJcNXF2K5tvmD80tYHf1ldP0tLRz7c2Lztn2Vj8z9T8X5L8nERuvXouK8qzhn407tlQEXP15MJcyovSsdudQxc6167Iw9nSx/tHW9mw2ju13Z7KVhpa+zAMg0VzM6bludbd5+JHv68kzm7lng0Lh75b8XYrL713igf1Ef6/L5xPSpJ504B5I5uk8qJ07DbnUJIqL0yLWixTeSh54Hg7vf1uVi7KAeCTlxbw3uFWHn/lBF/dtDisveX1Gyepb+7l7z9dPuIQMmsuyuXHz1axp7KViytygt7ugeNnhj4XiwVWLsph9fKzRyYvzU+dlj8WInJGG03E7TGGLqie39FAbVPP0N/3bFhIWRR/K0KtvWuAh5+upLGtjztuLMVisZz13bp0iYNjDd1Dyam1o5/sdPNNWzLtmpn754q5cOEsDANaOgaiEoe/79Czb9Xz0JMT7zu083ALaUn2ocn2khJsqE8Uc15JJp4JTpEykabqbZ0D7DrcwnUX541anbi0NJPZWQm8sruRYKdr2XfsDK/sacRm9Y6TZ5cqPBFGpfmpfHp16VkXO4F3+xeUZw3N9eD2GPz02SoO1kyPTuHN7f1874nDNLf3c/fN5VxQlnVOmeRE+9D3e3dlK9/45T4ef7WG53dEt0vLcNPuDgq8J+eK8+byLz/6C1tfO8GiuemkRbjeubL2o4eS7klUNV6zIo+LyrPP+lKtWJgNCycWR3V9J9/X3tY7Vks9n7ysgBsuGb35dlZaPN/4/HlkpY5+vKwWC2suzOPxV09QXd/FgnGuPLv7XPz6pRpy0hPYeFUxx5xdUoUnoqqiOJ0X3vVWA1otFhLjrQHfVw9WqwVLDDaqMAyDnz9XTXefi3s2LKRkzvjfsbLCNEryUvjzB00AxO1wco8yx7P7aZmgwHu1dNt1JRyoaY94n4i+ATflhWnYbRZcbmNS01KXzEmlZM7Iy/YdO8Pxhm7WXV4w5jYGXR4ef/UEHt8Xz2PAn/eeHkpQb+1rYlZGAvPzUznW0MWbH5zmqgtzKS0Yv6rj0iU5vLizgVNtfeMmKP36STp7Brn75jKKc1NYKHdOIsr81YD+KviSPG9fP4Btf62nqr6LlRU59PS7wn4xFYpWx4GPE/52bQkewwi6k3t6chxLSjKoqu/CAAbdBq+8d4rSdQsmFUsoTdsEBTAnJ4k5OUkAeAwjIs1MPR6DR545SlZqPPeqCv7n+WP09Lsonh1836HX329kXl7KqFc/lbWdvPTeKc4ryWD+GCf0qdY+Glr7sFktvoFXLWzyzTzqdnvQr5+kf9AzlEgBPqg+w71BXD3Fx9n49h3jN5LYW9XGjoMt3HhJPsVBTJchRKSM9rxzdmYif/mwicdfPQGAzerk7pvLJtVqdTwHatr5yTNHvc/HbJNrdRzYrD5ukn2c/A1L/L8DhRP4vQqnafcMaiT7j7fzrccO0NXrCvu+XtzVQGVtJxXF3quuz14zl+4+N+8ebAlq/a6eQfQbtew52jZqmRsvzSczNe6su6NAHT2DgHcCv+/csYyv3FLBussLuGdDBecv8NZH22xWvvvF5fxf6xdQFHAyejzBD2XkT3xNZ/pGLXOkrpPCWUnccMkot4NCmMzlS2dx1QW5Q6/dHoNn3qoHvFVo3VP8Hentd7PjYDM/+kMl//37SgbdBh7jo0cBE/XyroahxDLZbfjvKNdfXsA/31LBjb5aljc/PM2uw8H9doXDtL6D8stIieNUWx9PvnGSv7t+ftj2c8zZxR/frufihdlcusTbuq2iOJ15eSmcHuNHPNDuo214PAYrx2gdlxhvY8PqYn7+XDV/3nuaTwR8mfYdO8Ojf6rms2vmsXJRDpmp8WSmxo94RZWUYOP8Bd4m3g89eXhS/Y+ee8fJS7tOcf+dIzdXVauL6RtwY7fNiGshMU34+w253QZWq4U1F3m/YydP93D/44eoKErjwvJsMlPjqGvqGbUacHhL3t++XMNfDzTjchtkpcWzojyL96vOeC80LVBV34nL7Qnq+2IYBs++Xc/7VWewWMDC1CbpHH5HaRgGuw63UlnbyQdVZ7j16rmkRrhJelB7U0qtBR4GbMCjWuv7hy1PAH4FXAS0ALdorWt8y74GbAbcwD9orbeHLPogFc1O5vqVc/jTDicrFmazdH5myPfR2+/iF88fIystnk1r5g49YLVYLPw/GyuwBfkDvetwK3OyEymclTRmuYvKs/hLcTrP+hJiSpKd198/jX7jJIWzkifUZHZ4ffxEqgcuLMviuXecvPnhaa5f9VHji0MnOshIiSPfkUSiTMctYsxo34nUJDvXrshjT2Urv3m5Zqh8nM3bsKClvZ+/7Gsi3m5lYNBDVX0nHoOhqreM1DiuPH82KxZmUzInBavFQrWzi/pWNzX1rby9v5mHn67krk8uGLd/0m9eruGtfc1csdTBJYtyqApx4yOLxcI/fmYh23c28Md3nFTWdXLNRbkkJkVulJ5xfzWVUjbgEeB6YDFwq1Jq8bBim4E2rfUC4CHgu751FwMbgSXAWuDHvu1F3A2XzCHfkcRvXq6hpy/0VX1NZ/pxuT1svmE+ycMaZfiT06nW3jGbZbd29HO0rpOLK3LGbUFksVi49api1OpinC293P/4IZ54/STL5mfyT7dUTHgK7NL8VK5flT/hk65gVjKL56bz+vuncbk9AHT2DPKLP1Xz65eOB90MXQizGek7kZOewM0fK+S+25fyiQs+6sPnCqga93gMOnoGaWjtxV8D7x804KZLC7jlE96hw/zPxP1N4j9/XQm33zCfY84uvvu7Q2NWnQMsK83iU1cU8Llr5lFWlD6p7+94bFYLN1ySz9c2LSLebuHpN+t4/KXKSXWdmYxgLutXAlVa62Na6wFgK7B+WJn1wGO+v58CrlZKWXzvb9Va92utjwNVvu1FnN1m5bbr5tHePciuI63jlp/oNBfFuSl8e/OyUVvA7Tt2hm/+z36O1HaOuo1TrX2kJNq5uCI7qH3m5SSRm53Iw09XUnOqG6sFrr04L+J3LNesyPMe18Pe4/q7V0/Q0+/mc9fOi8mmukKMx2KxcHFFDnH2s/v1rVyUwz9vXMS/fW4JX1pfds7y8axalMM/fmYhXb0uvv/EYQYGz5588XRbHzsPeZ8JnV+ayfWr8iPyHSvOTeHSJQ4seGc1mOyzrokKpoqvAKgNeF0HrBqtjNbapZRqB3J87+8Ytu7YbaPDaF5eKt/4/BL6Bjz86Z16CmYlk50WT0ePi86eQXr73Vx1YS7Vzi6+/8Rh78i/Nif3blg4auI53dbHe0daWbtyDnH20fN9RXE66cl2tu9sGOp8O9zieRl870vLJzR8UGB/K/BOV7EgiGbiobRobjr5jiR2HGwhzm5ld2Ub668okLmcxLQ2XtX4ZKvOywrT+OqmRdQ19xAf99HFZnV9Jz9+tgqLBZaVZkb8QnTR3Axe3NkQ0bEyg0lQI/1aDq+3Ga1MMOuilLoTuBNAa43D4QgirLHZ7fYRt9PS3cb9j+9kYNAz4nqfvnoR9fvbh1rHudwGP362mmtXFXH5sjnMm/PRh+Jye/jeEztoaO7mxo+Vk5Mx9jxJn/zYfH67vZKO/jjmF3zUZNVut5OekUWc3Trhq6GV59l4/t0GXC4PdruVlecV4XCc23M8FEY7pgBfu+1iahs7+eETeylwpPDZtecF/dwtHMaK1Wwk1vCIRKwOh4NVyya/3G94rA4HLC7z/v32hw28truOD6tayEyN5z+/sIp8R+S7bDgcDv4zI4NDJ86waG4mC+eG53cmUDAJqg4oCnhdCDhHKVOnlLIDGUBrkOuitd4CbPG9NEIxpYPD4Rhxaoid+50MurzJyQJcvCib1efPJj0ljrTkONpaWyjIthFn9/YNslos5GTE8cyfj1F5ooW//3Q5AO8fbeXVPY1U1XXxxU+WYgx20dw8dnXgigUpPP26lSdePswXbio9K9bfvrCfnYda+MbnzyM+Lvgf9pwU+MfPLBy6SstJcYd0SoxAox1TgNqAgWAb23p478DJqPZEHytWs5FYw2O6xPrqzhr2Vp8BoLNngFpnE/H0RjK8ITkp8KkrS2hubp7Ssc3PD24y0mAS1C6gTClVAtTjbfSwaViZbcBtwDvAZ4DXtNaGUmob8LhS6kEgHygDdgYVWZj4B5P136auXp4b1K15R8/gUOOKD6ra+Om2agCsFshIDa5BQnKinSvPn82be5vo7XeRlOA9/P7mnFlp8RNKToHxRntYksBR5P19qaIdkxDTQcmcFD6sPoPB5IZNi2Xj/hpqrV3A3cB24JD3LX1AKXWfUmqdr9gvgBylVBVwL/BV37oHAA0cBF4Evqy1dg/fRyT5k4+/4+poH/TwFjzpyXHkZXubfp9s7D6r7EQeFq69eA7fuWPZUHICOO7soLGtb0Ijg5uNN/F7HwjLXE5ChI53lIeZ+d2ymLAZsOF0nlMLOGHhvL33z2vkvwubzNAihmHg9hjYbVb+tLOZP71dw/fuWm7quVnGO6ZTmV4k1KZL9Y7ZSKzhMdO+W74qvnEfuJv319DEptKxFXyNK7YeZtHcdNZdXsDbHzawZF6GqZNTMMxQ1SjEdDRTv1ux/YsYRVM5Yew2K9np8bzxwWmuWZHH5k8uwhiMzkNPIYQwKxkgLUquu3gOvf1unn27jvqmbixhnCFXCCFikdxBRcm8vBQKZyXz5w+asFiasE/yWZYQQkxXcgcVRf5pLiI5dIgQQsQKSVBR9LFls4iTptlCCDEiqeKLotL8VO5RFdS3uiM2fL0QQsQKSVBRVpqfyqplsdNfQwghIkWq+IQQQpiSJCghhBCmJAlKCCGEKUmCEkIIYUqSoIQQQpiSJCghhBCmJAlKCCGEKZlyPqhoByCEECLsxh0h24x3UJZQ/KeU2h2qbYX7v1iJNVbilFglVok1JmIdlxkTlBBCCCEJSgghhDlN5wS1JdoBTECsxBorcYLEGi4Sa3hIrCMwYyMJIYQQYlrfQQkhhIhhkqCEEEKYUszOB6WU2gD8B7AIWKm1fi9g2deAzYAb+Aet9fYR1i8BtgLZwB7gb7TWAxGI+wlgoe9lJnBGa718hHI1QCfef4NLa70i3LGNEMN/AF8Amnxv/avW+vkRyq0FHgZswKNa6/sjFuRHMXwP+CQwAFQDf6e1PjNCuRqidFzHO05KqQTgV8BFQAtwi9a6JlLxBcRR5IsjD/AAW7TWDw8rsxp4Fjjue+v3Wuv7IhlnQCw1jPGZKqUseI/7DUAP8Lda6z1RiHMh8ETAW/OBf9da/zCgzGqidFyVUr8EbgJOa63P872X7Yt5HlADKK112wjr3gZ83ffy21rrx0IRU8wmKGA/8GngZ4FvKqUWAxuBJUA+8IpSqlxr7R62/neBh7TWW5VSP8Wb0H4S7qC11rcExPoDoH2M4p/QWkd7JsOHtNbfH22hUsoGPAJcA9QBu5RS27TWByMVoM/LwNe01i6l1HeBrwH/MkrZiB/XII/TZqBNa71AKbUR7zl6y7lbCzsX8BWt9R6lVBqwWyn18gif6V+01jdFIb6RjPWZXg+U+f5bhfd7vipSgflprY8Ay2HofKgH/jBC0Wgd1/8FfoT34sTvq8CrWuv7lVJf9b0+63vlS2LfBFbgHWhht+/cPieRTVTMVvFprQ/5PvDh1gNbtdb9WuvjQBWwMrCA74rqKuAp31uPAZ8KZ7zD+WJQwO8iud8wWAlUaa2P+e5At+L9DCJKa/2S1trle7kDKIx0DOMI5jitx3sugvfcvNp3nkSU1rrBf4ehte4EDgEFkY4jhNYDv9JaG1rrHUCmUmpOlGO6GqjWWp+IchxDtNZvAq3D3g48J0f7nbwOeFlr3epLSi8Da0MRU8wmqDEUALUBr+s498uVg7dqzTVGmXD7GNCotT46ynIDeEkptVspdWcE4xrubqXUh0qpXyqlskZYHszxjrTbgRdGWRat4xrMcRoq4zs32/Geq1GjlJoHXAC8O8LiS5VSe5VSLyillkQ2srOM95ma8RzdyOgXp2Y5rgC5WusG8F64ALNHKBO242vqKj6l1Ct468GH+zet9bOjrDbSFefwtvTBlJm0IOO+lbHvni7XWjuVUrOBl5VSh31XOCE1Vqx4q0K+hffYfAv4Ad4f/0BhPZaBgjmuSql/w1tF9dtRNhOR4zqCqJ+XE6WUSgWeBv5Ra90xbPEeYK7WukspdQPwDN4qtGgY7zM123GNB9bhrYYezkzHNVhhO76mTlBa6zWTWK0OKAp4XQg4h5Vpxnubb/ddqY5UZtLGi1spZcf7/OyiMbbh9P3/tFLqD3iriEL+QxrsMVZK/Rx4boRFwRzvkAjiuN6G9yHv1VrrEb8gkTquIwjmOPnL1PnOkQzOrXKJCKVUHN7k9Fut9e+HLw9MWFrr55VSP1ZKOaLxzDSIzzRi52iQrgf2aK0bhy8w03H1aVRKzdFaN/iqRU+PUKYOWB3wuhB4IxQ7n45VfNuAjUqpBF9LvTJgZ2AB34/X68BnfG/dhrflTKSsAQ5rretGWqiUSvE9nEYplQJci7dRSEQNq6e/eZQYdgFlSqkS35XhRryfQUT5Wsj9C7BOa90zSploHtdgjtM2vOcieM/N10ZLtOHke+71C+CQ1vrBUcrk+Z+PKaVW4v0taYlclENxBPOZbgM+r5SyKKUuAdr91VZRMmrtiVmOa4DAc3K038ntwLVKqSzfY4Brfe9NmanvoMailLoZ+G9gFvAnpdQHWuvrtNYHlFIaOIi3qufL/hZ8SqnngTt8V1z/AmxVSn0beB/vFzJSzql/Vkrl4216fAOQC/xBKQXez+hxrfWLEYzP7wGl1HK8t+s1wBeHx+prNXc33hPSBvxSa30gCrH+CEjAW8UDsENrfZdZjutox0kpdR/wntZ6G95z8NdKqSq8d04bIxHbCC4H/gbYp5T6wPfevwLFAFrrn+JNoF9SSrmAXmBjNJIpo3ymSqm7AmJ9Hm8T8yq8zcz/LgpxAqCUSsbbkvOLAe8Fxhq146qU+h3eOyGHUqoOb8u8+wGtlNoMnAQ2+MquAO7SWt+htW5VSn0L70UYwH1a65Dc+ctQR0IIIUxpOlbxCSGEmD3nDc4AAAEaSURBVAYkQQkhhDAlSVBCCCFMSRKUEEIIU5IEJYQQwpQkQQkhhDAlSVBCCCFMSRKUEEIIU4rZkSSEiGVKqVK8Pe/X+OZdygc+BD6jtX4jqsEJYRIykoQQUaKU+gJwL95Bg/8A7NNa/1N0oxLCPKSKT4go0Vr/HDiKd66lOXinOBFC+EiCEiK6fg6cB/y31ro/2sEIYSZSxSdElPgmBNyLd+qX64GloRoFWojpQO6ghIieh4HdWus7gD8BP41yPEKYiiQoIaJAKbUeWAvc5XvrXuBCpdRnoxeVEOYiVXxCCCFMSe6ghBBCmJIkKCGEEKYkCUoIIYQpSYISQghhSpKghBBCmJIkKCGEEKYkCUoIIYQpSYISQghhSv8/OsNBGvn4N6QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data.plot(x=\"x\", marker=\".\", linestyle=\"--\") #pandasのplot使用\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. モデルの定義\n",
    "- モデルには以下の事前知識を仮定する。\n",
    "  - ピークの幅には上限があり、Cauchy分布確率密度関数のscaleは高々2程度である。\n",
    "  - 観測ノイズは正規分布に従い、そのsdは高々0.05である。\n",
    "- これらの条件を反映し、弱情報事前分布として半t分布を用いる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. WAIC用のstanコード"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "code_waic = \"\"\"\n",
    "data {\n",
    "  int N;\n",
    "  int K;\n",
    "  vector[N] X;\n",
    "  vector[N] Y;\n",
    "}\n",
    "\n",
    "parameters {\n",
    "  vector[K] mu;\n",
    "  vector<lower=0>[K] sigma;\n",
    "  real<lower=0> s_mu;\n",
    "  real<lower=0> s_noise;\n",
    "}\n",
    "\n",
    "model {\n",
    "  mu ~ normal(0, s_mu);\n",
    "  sigma ~ student_t(4, 0, 2);\n",
    "  s_noise ~ student_t(4, 0, 0.05);\n",
    "  for (n in 1:N) {\n",
    "    vector[K] line;\n",
    "    for (k in 1:K){\n",
    "      line[k] = log(sigma[k]^2) + cauchy_lpdf(X[n] | mu[k], sigma[k]);\n",
    "    }\n",
    "    target += normal_lpdf(Y[n] | exp(log_sum_exp(line)), s_noise);\n",
    "  }\n",
    "}\n",
    "\n",
    "generated quantities {\n",
    "  vector[N] log_likelihood;\n",
    "  for(n in 1:N){\n",
    "    vector[K] line;\n",
    "    for (k in 1:K){\n",
    "      line[k] = log(sigma[k]^2) + cauchy_lpdf(X[n] | mu[k], sigma[k]);\n",
    "    }\n",
    "    log_likelihood[n] = normal_lpdf(Y[n] | exp(log_sum_exp(line)), s_noise);\n",
    "  }\n",
    "}\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# コンパイル済みのモデルがあれば読み込み。なければコンパイルし保存。\n",
    "try:\n",
    "    with open('./model_peakwaic.pkl', \"rb\") as f:\n",
    "        stanmodel_waic = pickle.load(f)\n",
    "except FileNotFoundError:\n",
    "    stanmodel_waic = StanModel(model_code=code_waic, model_name=\"model_peakwaic\")\n",
    "    with open('./model_peakwaic.pkl', \"wb\") as f:\n",
    "        pickle.dump(stanmodel_waic, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2. WBIC用のstanコード\n",
    "- modelブロックでは、target記法により逆温度 $1/log(N)$ [Nはデータ数] の事後分布を定義している。\n",
    "- generated quantitiesブロックにおいて、上記事後分布の下での対数尤度を算出している。\n",
    "- stanコードは以前の実装例 <https://github.com/narrowlyapplicable/peak_separation_whth_WBIC/blob/master/wbic-mix-cauchy_lpdf.stan> そのまま。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "code_wbic = \"\"\"\n",
    "data {\n",
    "  int N;\n",
    "  int K;\n",
    "  vector[N] X;\n",
    "  vector[N] Y;\n",
    "}\n",
    "\n",
    "parameters {\n",
    "  vector[K] mu;\n",
    "  vector<lower=0>[K] sigma;\n",
    "  real<lower=0> s_mu;\n",
    "  real<lower=0> s_noise;\n",
    "}\n",
    "\n",
    "model {\n",
    "  mu ~ normal(0, s_mu);\n",
    "  sigma ~ student_t(4, 0, 2);\n",
    "  s_noise ~ student_t(4, 0, 0.05);\n",
    "  for (n in 1:N) {\n",
    "    vector[K] line;\n",
    "    for (k in 1:K){\n",
    "      line[k] = log(sigma[k]^2) + cauchy_lpdf(X[n] | mu[k], sigma[k]);\n",
    "    }\n",
    "    target += 1/log(N) * normal_lpdf(Y[n] | exp(log_sum_exp(line)), s_noise);\n",
    "  }\n",
    "}\n",
    "\n",
    "generated quantities {\n",
    "  vector[N] log_likelihood;\n",
    "  for(n in 1:N){\n",
    "    vector[K] line;\n",
    "    for (k in 1:K){\n",
    "      line[k] = log(sigma[k]^2) + cauchy_lpdf(X[n] | mu[k], sigma[k]);\n",
    "    }\n",
    "    log_likelihood[n] = normal_lpdf(Y[n] | exp(log_sum_exp(line)), s_noise);\n",
    "  }\n",
    "}\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# コンパイル済みのモデルがあれば読み込み。なければコンパイルし保存。\n",
    "try:\n",
    "    with open('./model_peakwbic.pkl', \"rb\") as f:\n",
    "        stanmodel_wbic = pickle.load(f)\n",
    "except FileNotFoundError:\n",
    "    stanmodel_wbic = StanModel(model_code=code_wbic, model_name=\"model_peakwbic\")\n",
    "    with open('./model_peakwbic.pkl', \"wb\") as f:\n",
    "        pickle.dump(stanmodel_wbic, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. MCMC実行 & IC算出\n",
    "- ピーク本数は少なくとも2本以上であることから、ピーク本数候補は2,3,4,5の４通りとした。\n",
    "- この４通りのについて、WAICとWBICを算出した。\n",
    "  - こちらもWBICについては以前の実装例そのまま。\n",
    "  - WAICについては導出がやや煩雑となるため、waic()関数を作成している。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1. WAIC算出"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def waic(samples):\n",
    "    tE = - np.mean(np.log(np.mean(np.exp(samples), axis=0)))\n",
    "    fVar = np.sum(np.mean(samples**2, axis=0) - np.mean(samples, axis=0)**2)\n",
    "    waic = tE + fVar/samples.shape[0]\n",
    "    return waic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/lib/python3.6/site-packages/pystan/misc.py:399: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  elif np.issubdtype(np.asarray(v).dtype, float):\n"
     ]
    }
   ],
   "source": [
    "dict_waic = {}\n",
    "for k_cand in [2,3,4,5]:\n",
    "    #pystan用にデータを辞書型にまとめる\n",
    "    standata = {\"N\":data.shape[0], \"K\":k_cand,  \"X\":data[\"x\"], \"Y\":data[\"y\"]}\n",
    "    fit = stanmodel_waic.sampling(data = standata, iter=10000, warmup=4000, seed=1234)\n",
    "    ms = fit.extract()\n",
    "    dict_waic[k_cand] = waic(ms[\"log_likelihood\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_waic = pd.DataFrame(data={\"n_peak\":list(dict_waic.keys()), \"WAIC\":list(dict_waic.values())})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2. WBIC算出"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/lib/python3.6/site-packages/pystan/misc.py:399: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  elif np.issubdtype(np.asarray(v).dtype, float):\n"
     ]
    }
   ],
   "source": [
    "dict_wbic = {}\n",
    "for k_cand in [2,3,4,5]:\n",
    "    #pystan用にデータを辞書型にまとめる\n",
    "    standata = {\"N\":data.shape[0], \"K\":k_cand,  \"X\":data[\"x\"], \"Y\":data[\"y\"]}\n",
    "    fit = stanmodel_wbic.sampling(data = standata, iter=10000, warmup=4000, seed=1234)\n",
    "    ms = fit.extract()\n",
    "    dict_wbic[k_cand] = - np.mean(np.sum(ms[\"log_likelihood\"], axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_wbic = pd.DataFrame(data={\"n_peak\":list(dict_wbic.keys()), \"WBIC\":list(dict_wbic.values())})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. 結果"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 算出されたWAIC, WBIC値は以下の通り。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>n_peak</th>\n",
       "      <th>WAIC</th>\n",
       "      <th>WBIC</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2</td>\n",
       "      <td>-2.684415</td>\n",
       "      <td>-95.823023</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3</td>\n",
       "      <td>-3.228095</td>\n",
       "      <td>-109.945743</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4</td>\n",
       "      <td>-3.236460</td>\n",
       "      <td>-100.927044</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>5</td>\n",
       "      <td>-3.239078</td>\n",
       "      <td>-98.677328</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   n_peak      WAIC        WBIC\n",
       "0       2 -2.684415  -95.823023\n",
       "1       3 -3.228095 -109.945743\n",
       "2       4 -3.236460 -100.927044\n",
       "3       5 -3.239078  -98.677328"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.merge(df_waic, df_wbic, on=\"n_peak\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*[表]候補となるモデル（ピーク本数）に対するWAIC値およびWBIC値*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- これらの値を棒グラフとし下図に示した。\n",
    "  - ただしWAIC値はピーク本数3~5での変化が小さいため、この領域を拡大した形で示した。そのためピーク本数＝２における値が上方に見切れている。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAHwCAYAAAC7apkrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmcZHV97//Xqarep3sWhtlYZJBFFlkUAUEIIFHiRty+PzVXmIv8CN4kehOTXGFwSQyEPMyNxtwrCT9jEJDgFxXJFS9GEERRZBFkHfZtGJhhmL33qjq/P+r0TE1P90wP092nq+f1fDzq0afO+Z6qT82jhn7znc/5niRNUyRJkiTVFPIuQJIkSZpKDMiSJElSHQOyJEmSVMeALEmSJNUxIEuSJEl1DMiSJElSHQOyJEmSVMeALEkNIEmSTyRJMpgkSeew/Q9sZ/+3hu17MEmSSpIkR4zw+l9MkuTJYftKSZL8SZIkdyVJsjFJkvVJktyXJMnSJElmj+fnk6SpxIAsSY3hZqAE/M7QjiRJ5gKHAS+NsP9w4Ja6fScA84B/Bc7b0ZslSdIE3AhcDETgNOBIYClwPHD2rn4gSZqqSnkXIEnasTRNn0uS5Cng7cAPs92nAQ8Bt4+wP6EuIAN/CHwbuBb4zyRJ/jJN057tvOWngN8FTkzT9Fd1+58FfuQMsqTpzBlkSWoct1ALwkNOA36aPYbvX5am6YsAWZj9MPCtNE3vAl4Ewg7e6+PAT4eF483SNF37mj6BJDUAA7IkNY5bgMOTJJmXPT8NuBX4GXDosP031513FvBYmqa/zZ5/ix23WRwEPDIuVUtSgzEgS1Lj+Gn28+1JkuwN7A/8LE3TNcADdfsPZOuAfB61UDzkKuDYJEkO3857JUA6bpVLUgMxIEtSg0jTdDXwW2rtFG8HfpOm6frs8K11+yvAbQBJkrwNOBT4n0mSlJMkKQMvAEW2P4v8GLULACVpt2NAlqTGMtSHPNR/POTWuv331AXnPwR+Qm0FiqPqHp8GPp4kSdso73M1cFqSJG8d6aAX6UmazgzIktRYbgH2A97P1gH5dmCfbP8tAEmSzAE+BFyVpulD9Q9qy721ULt4byT/mL3Oj5Mk+fMkSY5JkuR1SZKckSTJD6j1NUvStGRAlqTGcjswCLQCvxjamabpBuBeoJMt/cdnU+slvmH4i6Rp2k1tneMR2yzSNB0Efg/4HPARahcCPgj8LXAXW/c0S9K0kqSp12BIkiRJQ5xBliRJkurkfie9EMKXgDOBKrAKWBJjXDFszFHAZUAXtauzL44xfic7tpjanaHmAL8BPh5jHJi8TyBJkqTpZCrMIH85xnhEjPEoardJ/fwIY3qAs2KMhwFnAF8NIczKjv0d8JUY44HAWuATk1G0JEmSpqfcA3KMcUPd0w5GWJg+xvh4jPGJbHsFtZnmPUMICbUljb6bDf0W8PsTW7EkSZKms9xbLABCCBdTWzJoPXDqDsYeCzQDTwF7AOtijOXs8HJgr+2cex7ZFdsxxjfveuWSJElqMMmOBkxKQA4h3AwsGOHQ0hjjDTHGpcDSEMIFwB8DXxjldRZSu0Xq2THGajaDPNyoy3LEGC8HLh8at2LFitGGTohzTl4yqe+nfHzz9ivyLmHaeGVdHxf964N8/B378bY37pl3OZKkBrdo0aIxjZuUgBxjPH2MQ6+hti7nNgE5hNCVHbsoxnhntns1MCuEUMpmkfcGJjf1Spowe3S1UComrFzTl3cpkqTdSO49yCGEA+uevg9YNsKYZuB64MoY43VD+2OMKbXbq34o23U2IyyIL6kxFQoJ82a1snKtAVmSNHmmQg/ypSGEg6kt8/YccD5ACOEY4PwY47lAAE4G9gghLMnOWxJjvB/4H8C1IYS/Ae6jdvtUSdPE/DmtrFjdm3cZkqTdyO58Jz17kDUh7EEeX9f/fDn/ec/L/K9PvYliMfd/9JIkNbCsB3mHF+n520bSlLZgTivVasrq9d7/R5I0OQzIkqa0+bNbAXjZPmRJ0iQxIEua0ubPqQXklWvsQ5YkTQ4DsqQpraO1RGdbyZUsJEmTxoAsacqbP6eVl10LWZI0SQzIkqa8+bNdC1mSNHkMyJKmvPlzWtnYU6anr5x3KZKk3YABWdKUtyBbycJZZEnSZDAgS5ryhlaysA9ZkjQZDMiSprw9Z7ZQKCTOIEuSJoUBWdKUVywW2HNmCyudQZYkTQIDsqSGMH92q3fTkyRNCgOypIYwf04rq9b2Ua2meZciSZrmDMiSGsKC2a2UKylrNg7kXYokaZozIEtqCK5kIUmaLAZkSQ1hKCCvXNubcyWSpOnOgCypIXS2lWhvKbqShSRpwhmQJTWEJEmYP8eVLCRJE8+ALKlhzJ/d6gyyJGnCGZAlNYz5c1pZt2mQvoFK3qVIkqYxA7KkhrFgdu1CvVW2WUiSJpABWVLDcKk3SdJkMCBLahjzZrWSACudQZYkTSADsqSG0VQqsMfMFi/UkyRNKAOypIYyf7ZLvUmSJpYBWVJDmT+nlVVr+0jTNO9SJEnTlAFZUkNZMLuV/sEq6zYN5l2KJGmaMiBLaiiuZCFJmmilvAsIIXwJOBOoAquAJTHGFcPGHAVcBnQBFeDiGON3smPfBo4BBoG7gD+MMTq1JE1T87O1kFeu7eOQ13XlXI0kaTqaCjPIX44xHhFjPAr4IfD5Ecb0AGfFGA8DzgC+GkKYlR37NvAG4I1AG3DuJNQsKSezZjTR0lRg5ZrevEuRJE1Tuc8gxxg31D3tALa58ibG+Hjd9ooQwipgT2BdjPFHQ8dCCHcBe09guZJyliSJK1lIkiZU7gEZIIRwMXAWsB44dQdjjwWagaeG7W8CPg58ejvnngecBxBjZO7cubtWuDQCv1cTb9+FM3n8+XX+WUuSJsSkBOQQws3AghEOLY0x3hBjXAosDSFcAPwx8IVRXmchcBVwdoyxOuzw14HbY4w/H62OGOPlwOXZ03T16tU7+UmkHfN7NfFmtSe8sraXFS+torlpKnSKSZIawaJFi8Y0blICcozx9DEOvQa4kRECcgihKzt2UYzxzmHHvkCt5eIPd7FUSQ1gwZxWUuCVdX3stWd73uVIkqaZ3KdeQggH1j19H7BshDHNwPXAlTHG64YdOxd4J/DREWaVJU1DQytZuNSbJGkiTIUe5EtDCAdTW+btOeB8gBDCMcD5McZzgQCcDOwRQliSnbckxng/8M/Zeb8KIQB8P8b415P7ESRNps1rIXuhniRpAiS78e1a0xUrVux41Dg65+Qlk/p+ysc3b78i7xJ2C5/9l/s5aJ8uznnX/nmXIklqEFkPcrKjcbm3WEjSazF/ThsrnUGWJE0AA7KkhjR/disvr+ljN/5XMEnSBDEgS2pIC+a00jdQYUNPOe9SJEnTjAFZUkMaulDPW05LksabAVlSQ9q81Jt9yJKkcWZAltSQ5nQ101RKWOlayJKkcTYV1kGWpJ3W3NTE0iXHMqO9mWJLiaTSR7lsP7IkadcZkCU1nFKpxMaBJr7+vbtZtbaXebPbuHDJW+hsxpAsSdpltlhIajhpsZVLrqiFY4BVa3u55Iq7GUxbKFe847wkadc4gyxNI9e943/kXcKkePe1f7M5HA9ZtbaXtZsGueDr99I00E9Lfy8tfX209PfS3F/72TL0s6+P5oE+Cg28hvKH//Pv8i5BkqYtA7KkhjOwdiPzZrdtFZLnzW6jvHot+z/xMP0trfS3tjHQ0sqmzi4GmltJC8P+wSxNaR7op3lYcN4SpGvbTQP9DR2kJUk7z4AsqeE8cNl3+czSc/mf1z+yuQf5M+8/lAcu/gb7P/3UNuNTYKC5hf6WNvpb2+hvaWWgpbX2PPu5sWsWA82tkCTDTk6HzUDXbfdtmZ1uHugn2eadJUmNyIAsqeGsfOApHrj4G3zmkx+ieXYnA2s38sDF32DlA9uGY4AEaBnop2WgHzauG/V1q0nCQHMLA3XBub91S5Dua21jw8w5DLS0bvse1SrNA32bg3NLf98IrR29NA0OGKQlaYozIEtqSCsfeIqffHJ8+3ALaUprfx+t/dtfW7maJMNmoLdsD7S20dvewbrZcxlsbtnm3KRaHbEnuv55c3+fQVqScmRAlqSdVEhTWvt6ae3b/m2uq0lhS4De3NqxZXa6p30Ga+fsSbmpeZtzk2plm+Bcm5WuBeoVr/Yyq6OJtpYiyfC2EEnSLjEgS9IEKaRV2vp6aOvrgfWjj6sUCvS3tGWz0q1b9Ur3t7TRPaOLNXvM2ypI33fFQwA0lRJmdjQza0YTMzuamDWjmZl127M6mpg5o5nW5oJBWpLGyIAsSTkrVqu093bT3tu93XGVYnFzaD76r85hffcg6zYNsn7TAOu6B1n+Si8PPbOe/sFt14JuaSpkwXkoNGdhuqMpC9e1YN3aXJyojylJDcOALEkNolip0N7TTXtPN295wx6jjusbqLB+0yDrugdqP7MAvWHTIOu6B3luZTdrnxpksLxtkG5tLmyZkR4K1PUz0tl2c5NBWtL0ZUCWpGmmtblI65wi8+dsu9rGkDRN6Ruo1Gagu2shen0WoNdvGmDdpkGeXtHN+u61DJa3XQe6raW4uX1jtPaOmR1NNJW8YaukxmNAlqTdUJIktLWUaGspsXCPtlHHpWlKT39l80z0+u7BbWann1jex7pNg1Sq2wbpjtZiLUR3DJuRntHErGy7q6OJUtEgLWnqMCBLkkaVJAkdrSU6Wkssmrv9IN3dV9ncD715RnrTIOu7a8H6pef7WN89SHWEID2jrbS5F3qr9o66fumu9hJFg7SkSWBAliTtsiRJmNFWYkZbib32HH1cNU3Z1Fve7oz0i6t72NA9yPAcnQCd7aXNFxRuHai3XHDY2d5EseCKHZJeOwOyJGnSFJKErvYmutqb2Gde+6jjqtWUjb3lzf3QQwG6vl/6+VU9bOweZPh8dJJAV/uWNo6ZI1xkOGtGMzPaSxR2sPRdqVQiLbZSrSYUCilJpY9yuTwOfxKSpjIDsiRpyikUEmZ21MLsvvNHH1eppmzoHsxmogeyiwy3zE6v3TjAMy9vYmPPtqG2kEBXx7Dl7ur6pRfN7aS1vY2//dadrFrby7zZbVy45C10NmNIlqY5A7IkqWEVCwmzO5uZ3dkMdIw6rlypbg7S267cMcDq9f08+eImuvu2BN8LlxzLV667m1Vra3dMXLW2l0uuuJtPfvAI/vWGhygVCzSVCjQVE5pKBUrFhOamAk3FAqW6/bXnW7abSsnmMc2loXOz/XVjSsXaa3qDF2nyGZAlSWO25tOfzruEXTIze4ymTIGNxTY2FNtY0PLmzeF4yKq1vXSWoOu5xygnRcpJkb7sZzkpMli3PfTYVaVqmRIVmtIqpbRMU1qhlFYppZXNj6a67VJ2vGmUY01phRI7HjvvH/7eXm7ttgzIkiRlSlSZXelmdqWb5g1rmDe7bauQPG92GzPWr+Ljq28f0+tVgQqFrQLz4Ha3xzCWoee1sT2FlhHCeYHBpEQ12YVVP75yD4WEEWe4S9kseHNp29nxzTPr2Zih2fWmUt344bPvw2bcm0oJpVJhhz3i0kQxIEuSNILe71/LZ8/6JJd+99HNPcif/dAh9F552ZhfowAUqNKUViEdnLhiR1EhoZIUtp3ZZntBvRa8m97zPgbLVcrlKoOVlMFylcFKlXI5ZSDbHihX6e6rbQ+WU8qVam1cuTZ+2wX9dk6pmGwJ53XBe6T2li0BfGztLU2lpDY+5/YWLwSdmnIPyCGELwFnUvsf7VXAkhjjimFjjgIuA7qACnBxjPE7w8b8E/BfY4wzJqVwSdK01vfEE7ReeRl//YGPQGcXbNxA75WX0ffEE3mXNmZFUopphea0stPnzjn+/F167zRNqVZTBsppFqyz4DwUoitptm/rAF4fsMuV+nOy8XX7egcqbOgZ3DqcV4bO3dV4ztaBuphsG8Drwnn9zHfTNqG7sM1rNZUKzOpqpb29lUuv9ELQqSb3gAx8Ocb4OYAQwqeAzwPD/1b2AGfFGJ8IISwC7g0h/DjGuC477xhg1mQWLUma/vqeeIK+v/tS3mU0pCRJKBYT2orQxq73Yu+sappS2Ry8s5/lYeG8siWMD82ODw/pI86OZ9v9veUtYwbrg3464g1xhrtwybF87cp7trkQ9JJPHg/lTRP9R6TtyD0gxxg31D3tgG3/RSbG+Hjd9ooQwipgT2BdCKEIfBn4GPD+CS5XkiRNpFu+OS4vU8geTePxYsXs0Tz2UyoplKsJg2mBwWpCOfs5mCYMVgsMpgmLuo4b8ULQas8muG18/hymtLefk3cFo8o9IAOEEC4GzgLWA6fuYOyx1L6iT2W7/hj4jxjjSyGEHb3PecB5ADFG5s6du4uVS9vye6XJkNf3bE0u76o85PUdW53Lu46/YgLFYkoLo7e3FCubRrwQtDDYs52zpo+p/PtyUgJyCOFmYMEIh5bGGG+IMS4FloYQLqAWeL8wyussBK4Czo4xVrN2iw8Dp4yljhjj5cDl2dN09erp8tdQU4nfK00Gv2eaaH7HJl6y7Odc+LEzuOSaB7f0IH/sjSTLbsq7tEmRx3ds0aJFYxo3KQE5xnj6GIdeA9zICAE5hNCVHbsoxnhntvto4ADgyWz2uD2E8GSM8YBdr1qSJGnilF99kc5lN3HJx06i2tROYbCHZNlNlF99Me/Sdnu5t1iEEA6MMQ5dEvw+YNkIY5qB64ErY4zXDe2PMd5I3cx0CGGT4ViSJDWK8qsvwh3XAuwWbRWNIveADFwaQjiY2jJvz5GtYJGtTHF+jPFcIAAnA3uEEJZk5y2JMd6fQ72SJEmaxnIPyDHGD46y/x7g3Gz7auDqMbyWayBLkiRpl+zwHpQhhDkhhDNGOXZGCGH2+JclSZIk5WMsN2m/CHjzKMeOBpaOXzmSJElSvsYSkN8D/Msoxy6ndptoSZIkaVoYS0BeEGMcbaG6NcD8caxHkiRJytVYAvLabJWJkRwErBvHeiRJkqRcjSUgXw98LYTQVr8ze/4V4LsTUZgkSZKUh7Es8/Y54KfA0yGEm4CXgIXAO4EXGOW20JIkSVIj2uEMcoxxI3ACtaDcChyT/fwccFJ2XJIkSZoWxnSjkBjjIPCN7CFJkiRNWzsMyCGEc3Y0Jsb4zfEpR5IkScrXWGaQP76D4ylgQJYkSdK0sMOAHGM8dTIKkSRJkqaCMfUgjySEMAf4KHB2jPHY8StJkiRJys9OBeQQQgl4N3A28C7gReCfJ6AuSZIkKRdjCsghhDdTC8UfBYrUbh7SB7w1xrhq4sqTJEmSJtcO10EOITwE/BxYBPwhsCDG+Amgd4JrkyRJkibdWG413Q5UqAXiHmBwQiuSJEmScjSWO+ntT63veAD4DvByCOGfqN1NL53Y8iRJkqTJNZYZZGKMt2dtFQuAzwAHA53AbSGE/zaB9UmSJEmTaix30gvA7THGl2OMvcDVwNUhhL2As4A/Ab4+sWVKkiRJk2Msq1j8DfD6EMJTwO3Az4CfxRifB/42e0iSJEnTwlh6kA+itoLFUmoX6n0GeDqE8FwI4aoQwrkTXKMkSZI0aca0DnKMcSVwXfYghDALOA/4M+BjwDcmqkBJkiRpMo31RiEJcBRwcvY4AVgBRGprJEuSJEnTwlgu0vsh8CbgMeAXwOXAkhjjxgmuTZIkSZp0Y1nm7WCgH3gGeAp40nAsSZKk6WqHM8gxxgNDCPPZ0l7x30MIc4E7qLVX/CLGeP/ElilJkiRNjl29SO8iYE+gOFEFSpIkSZPptV6k9zZgFnAP8M1dKSCE8CXgTKAKrKLW37xi2JijgMuALqACXBxj/E5dbX8DfDg7dlmM8Wu7UpMkSZJ2XzvsQQ4h3AisoXaB3pnAA8BHgJkxxhNjjBfuYg1fjjEeEWM8Cvgh8PkRxvQAZ8UYDwPOAL6azWIDLAH2Ad4QYzwEuHYX65EkSdJubCwzyD8HLgbujjEOjncBMcYNdU87gHSEMY/Xba8IIayi1tqxDvgk8LEYYzU7vmq8a5QkSdLuYywX6V060UWEEC4GzgLWA6fuYOyxQDO1FTUAXg/8PyGE9wOvAJ+KMT4xyrnnUeudJsbI3Llzx+cDSHX8Xmky5PU9W5PLuyoPeX3HVufyrsrDVP59OaYe5F0VQrgZWDDCoaUxxhtijEuBpSGEC4A/Br4wyussBK4Czh6aMQZagL4Y4zEhhA9Q64k+aaTzY4yXU1vHGSBdvdq/hhp/fq80GfyeaaL5HdNEy+M7tmjRojGNm5SAHGM8fYxDrwFuZISAHELoyo5dFGO8s+7QcuB72fb1wL/tQqmSJEnazY3lRiETKoRwYN3T9wHLRhjTTC38XhljvG7Y4R8Ap2XbvwM8jiRJkvQaTcoM8g5cGkI4mNoyb88B5wOEEI4Bzo8xngsEasvL7RFCWJKdtyS7QcmlwLdDCH8KbALOneT6JUmSNI3kHpBjjB8cZf89ZGE3xng1cPUo49YB756wAiVJkrRbyb3FQpIkSZpKDMiSJElSHQOyJEmSVMeALEmSJNUxIEuSJEl1DMiSJElSHQOyJEmSVMeALEmSJNUxIEuSJEl1DMiSJElSHQOyJEmSVMeALEmSJNUxIEuSJEl1DMiSJElSHQOyJEmSVMeALEmSJNUxIEuSJEl1DMiSJElSHQOyJEmSVMeALEmSJNUxIEuSJEl1DMiSJElSHQOyJEmSVMeALEmSJNUxIEuSJEl1DMiSJElSHQOyJEmSVMeALEmSJNUp5V1ACOFLwJlAFVgFLIkxrhg25ijgMqALqAAXxxi/kx17O/BlamF/U3b+k5P3CSRJkjSdTIUZ5C/HGI+IMR4F/BD4/AhjeoCzYoyHAWcAXw0hzMqOXQb8QXb+NcBFk1G0JEmSpqfcZ5BjjBvqnnYA6QhjHq/bXhFCWAXsCazLxndlh2cCK4afL0mSJI1V7gEZIIRwMXAWsB44dQdjjwWagaeyXecCPwoh9AIbgOO3c+55wHkAMUbmzp2768VLw/i90mTI63u2Jpd3VR7y+o6tzuVdlYep/PtyUgJyCOFmYMEIh5bGGG+IMS4FloYQLgD+GPjCKK+zELgKODvGWM12/ynwrhjjr0MIfwH8A7XQvI0Y4+XA5dnTdPVq/xpq/Pm90mTwe6aJ5ndMEy2P79iiRYvGNG5SAnKM8fQxDr0GuJERAnIIoSs7dlGM8c5s357AkTHGX2fDvgPctOsVS5IkaXeV+0V6IYQD656+D1g2wphm4HrgyhjjdXWH1gIzQwgHZc9/F3h0omqVJEnS9DcVepAvDSEcTG2Zt+eA8wFCCMcA58cYzwUCcDKwRwhhSXbekhjj/SGE/xf4XgihSi0wnzPZH0CSJEnTR+4BOcb4wVH230PWSxxjvBq4epRx11ObXZYkSZJ2We4tFpIkSdJUYkCWJEmS6hiQJUmSpDoGZEmSJKmOAVmSJEmqY0CWJEmS6hiQJUmSpDoGZEmSJKmOAVmSJEmqY0CWJEmS6hiQJUmSpDoGZEmSJKmOAVmSJEmqY0CWJEmS6hiQJUmSpDoGZEmSJKmOAVmSJEmqY0CWJEmS6hiQJUmSpDoGZEmSJKmOAVmSJEmqY0CWJEmS6hiQJUmSpDoGZEmSJKmOAVmSJEmqY0CWJEmS6hiQJUmSpDoGZEmSJKlOKe8CAEIIXwLOBKrAKmBJjHHFsDGvA74PFIEm4J9ijP+cHXszcAXQBvwI+HSMMZ20DyBJkqRpY6rMIH85xnhEjPEo4IfA50cY8xJwQjbmOOCzIYRF2bHLgPOAA7PHGZNQsyRJkqahKTGDHGPcUPe0A9hm9jfGOFD3tIUs3IcQFgJdMcZfZc+vBH4f+L8TVrAkSZKmrSkRkAFCCBcDZwHrgVNHGbMPcCNwAPAXMcYVIYRjgOV1w5YDe41y/nnUZpqJMbJo0aKRhk2Ym578z0l9P+1+Pv3QVXmXoGlu0XXX5V2CpruPX5R3BRJJmk5Oq24I4WZgwQiHlsYYb6gbdwHQGmP8wnZeaxHwA+C9wL7A38YYT8+OnQT8ZYzxveNZv16bEMI9McZj8q5D05ffMU0Gv2eaaH7HppZJm0EeCrBjcA21WeJRA3I2c/wwcBJwB7B33eG9gRUjnihJkiTtwJS4SC+EcGDd0/cBy0YYs3cIoS3bng2cCDwWY3wJ2BhCOD6EkFBr07hh+PmSJEnSWEyVHuRLQwgHU1vm7TngfICsv/j8GOO5wCHA/wwhpEAC/H2M8cHs/E+yZZm3/4sX6E0ll+ddgKY9v2OaDH7PNNH8jk0hk9aDLEmSJDWCKdFiIUmSJE0VBmRJkiSpzlTpQdY0k61ZfSW1pf2qwOUxxn/MtypNJyGEVuB2ajcOKgHf3d7ykNJrFUIoAvcAL8YY35N3PZpeQgjPAhuBClB2qbepwRlkTZQy8JkY4yHA8cAfhRAOzbkmTS/9wGkxxiOBo4AzQgjH51yTpqdPA4/mXYSmtVNjjEcZjqcOZ5A1IbLl917KtjeGEB6ldofDR3ItTNNGjDEFNmVPm7KHVx1rXIUQ9gbeDVwM/FnO5UiaJAZkTbgQwn7A0cCvcy5F00z2T9/3Urv9/P+OMfod03j7KvCXQGfehWjaSoH/zJax/ZcYo8u9TQG2WGhChRBmAN8D/nuMcUPe9Wh6iTFWYoxHUbuD5rEhhMPzrknTRwjhPcCqGOO9edeiae3EGOObgN+j1o54ct4FyYCsCRRCaKIWjr8dY/x+3vVo+ooxrgNuA87IuRRNLycC78suoroWOC2EcHW+JWm6iTGuyH6uAq4Hjs23IoEtFpog2W2//xV4NMb4D3nXo+knhLAnMBhjXJfdhv504O9yLkvTSIzxAuACgBDCKcCfxxj/S65FaVoJIXQAhexanQ7gHcBf51yWMCBr4pwIfBx4MIRwf7bvwhjjj3KsSdPLQuBbWR9yAYgxxh/mXJMk7Yz5wPUhBKhlsmtijDflW5LAW01LkiRJW7EHWZIkSapjQJYkSZLqGJAlSZKkOgZkSZIkqY4BWZIkSapjQJYkjVkI4bYQwrl51yFJE8mALEmSJNUxIEuSJEl1vJOeJE1xIYRngf8FnAW8DrgJODvG2DfK+FOAq4FrHSBrAAAgAElEQVSvA38GbAKWxhi/nR1vAS4GAtACXA/8aYyxN4QwG7gKOI7a74g7gPNjjMtHeJ+FwI+BK2OMfz9en1eS8uYMsiQ1hgCcASwGjgCW7GD8AmAusBdwNnB5COHg7NjfAQcBRwEHZGM+nx0rAP9GLYjvC/RSC+dbFxPCfsDPgP9lOJY03TiDLEmN4WsxxhUAIYT/Qy3c7sjnYoz9wM9CCDfWTg1/A/y/wBExxjXZ610CXANcEGN8Ffje0AuEEC4Gbh32uocCF2Xj/30XP5ckTTkGZElqDC/XbfcAi3Ywfm2Msbvu+XPZOXsC7cC9IYShYwlQBAghtANfoTZbPTs73hlCKMYYK9nzPwCeBL772j6KJE1ttlhI0vQ0O4TQUfd8X2AFsJpa28RhMcZZ2WNmjHFGNu4zwMHAcTHGLuDkbH9S91pfzF7nmhBCcSI/hCTlwYAsSdPXX4UQmkMIJwHvAa6LMVaB/w/4SghhHkAIYa8QwjuzczqpBeh1IYQ5wBdGeN1B4MNAB3BVCMHfJZKmFf+jJknT08vAWmqzxt+mthLFsuzY/6DWInFnCGEDcDO1WWOArwJt1GaI76S2YsY2YowDwAeAecA3DcmSppMkTdO8a5AkjaOhZd5ijHvnXYskNSL/j1+SJEmq4yoWktSAQggXAheOcOjn1NY5liS9RrZYSJIkSXVssZAkSZLqGJAlSZKkOgZkSZIkqY4BWZIkSapjQJYkSZLqGJAlSZKkOgZkSZIkqY4BWZIkSapjQJYkSZLqGJAlSZKkOgZkSZIkqY4BWZIaQJIkn0iSZDBJks5h+x/Yzv5vJUmyJEmStO7RmyTJsiRJPjNs/BeTJHly2L5SkiR/kiTJXUmSbEySZH2SJPclSbI0SZLZE/dpJSlfBmRJagw3AyXgd4Z2JEkyFzgMeGmE/YcDt2S7KsDC7HEo8BXg0iRJPj7amyVJ0gTcCFwMROA04EhgKXA8cPY4fS5JmnIMyJLUANI0fQ54Cnh73e7TgIeAG0bYn7AlIJOm6cvZ45k0Tf8FeAA4Zjtv+Sngd4F3pmn692ma3p2m6bNpmv4oTdP3At8alw8mSVOQAVmSGsctbBuEf5o9hu9flqbpi8NfIKk5BTgE+OV23uvjwE/TNP3VSAfTNF27c6VLUuMwIEtS47gFODxJknnZ89OAW4GfAYcO239z3XnFJEk2JUmyCRjIzvl6mqbf2c57HQQ8Mq7VS1KDMCBLUuP4afbz7UmS7A3sD/wsTdM11FomhvYfyNYBuQIclT2OBD4GnJUkyRe3814JkI5v+ZLUGEp5FyBJGps0TVcnSfJbau0UzcBv0jRdnx2+tW5/Bbht2Ln1K1Q8kiTJfsBfJUnyt2ma9o/wdo9RuwBQknY7ziBLUmMZ6kMe6j8ecmvd/nvqgvNoytQmSVpGOX41cFqSJG8d6aDLvEmazgzIktRYbgH2A97P1gH5dmCfbP8tw09KkmRB9tg3SZL3Av8duCVN0w2jvM8/Zq/z4yRJ/jxJkmOSJHldkiRnJEnyA+Cs8ftIkjS12GIhSY3ldmAQaAV+MbQzTdMNSZLcCxzL1v3HAEVqayVDbeZ4OXA98IXR3iRN08EkSX4P+CNqK1r8FbXWjaeA63CZN0nTWJKmXoMhSZIkDbHFQpIkSarTMC0WIYQjgX8GZgDPAn8QY9yQHbsA+AS1f/77VIzxx3nVKUmSpMbWSDPI3wA+G2N8I7Xeub8ACCEcCnyE2nJEZwBfDyEUc6tSkiRJDa2RAvLB1C5OAfgJ8MFs+0zg2hhjf4zxGeBJahepSJIkSTutYVosgIeA9wE3AB+mtpwRwF7AnXXjlmf7thFCOA84DyDG+OYJq1SSJElTVbKjAVMqIIcQbgYWjHBoKXAO8LUQwueB/wAGsmMjfcgRl+aIMV4OXD40ZsWKFbtWsCRJkhrGokWLxjRuSgXkGOPpOxjyDoAQwkHAu7N9y9kymwywN2DylSRJ0mvSMD3IIYR52c8CcBG1FS2gNpv8kRBCSwhhMXAgcFc+VUqSJKnRNUxABj4aQngcWEZthvjfAGKMDwMReAS4CfijGGMltyolSZLU0HbnO+nZgyxJkrQbyXqQd3iRXiPNIEuSJEkTzoAsSZIk1TEgS5IkSXUMyJIkSVIdA7IkSZJUx4AsSZIk1TEgS5IkSXUMyJIkSVIdA7IkSZJUx4AsSZIk1TEgS5IkSXUMyJIkSVIdA7IkSZJUx4AsSZIk1TEgS5IkSXUMyJIkSVIdA7IkSZJUx4AsSZIk1TEgS5IkSXUMyJIkSVIdA7IkSZJUx4AsSZIk1SnlXcBYhRCOBP4ZmAE8C/xBjHFDCGE/4FHgsWzonTHG83MpUpIkSQ2vYQIy8A3gz2OMPwshnAP8BfC57NhTMcaj8itNkiRJ00UjtVgcDNyebf8E+GCOtUiSJGmaaqQZ5IeA9wE3AB8G9qk7tjiEcB+wAbgoxvjzkV4ghHAecB5AjJG5c+dObMWSJElqOEmapnnXsFkI4WZgwQiHllLrMf4asAfwH8CnYox7hBBagBkxxldDCG8GfgAcFmPcsIO3S1esWDGO1UuSJGkqW7RoEUCyo3FTKiCPVQjhIODqGOOxIxy7jVqv8j07eBkD8gQqlUqkxVaq1YRCISWp9FEul/MuS5Ik7cbGGpAbpgc5hDAv+1kALqK2ogUhhD1DCMVse3/gQODpvOpULRxvHGjiwsvu5LxLf8qFl93JxoEmSqVG6uiRJEm7q4YJyMBHQwiPA8uAFcC/ZftPBh4IIfwW+C5wfoxxTU41CkiLrVxyxd2sWtsLwKq1vVxyxd2kxdacK5MkSdqxhmyxGCeT3mJxzslLJvX98vLl6/+ZC694cJv9l/zXI/mL3z8vh4om1zdvvyLvEiRJ0gimXYuFGsfGNeuZN7ttq33zZrexqdpM4bATodScU2WSJEk7ZkDWuIv/9G/8ye+/YXNInje7jT95/yFc9+8/ofD6Iymd/l8o7Hc4JDv8HzhJkqRJ51VTGnfL7l/GFV/8ez7xJ/+Vzjkz2bhmPVd84cssu38ZdM2lePiJFI/8HQqL30jl4TtIVz2fd8mSJEmb2YM8iXaXHuSxSBYspnjYiSQzZlJd+RyVh+6ATWvzLmtc2IMsSdLUNNYeZGeQlYv05Wcor3qOwuIjKBx8DKVTP0L12YeoLrsLBvvzLk+SJO3GDMjKT7VK9an7qb6wjMIbjqWw+HAKex9E9bF7qD7zIKTVvCuUJEm7IS/SU/4G+qg+cDvlW79Dum4VxTe+jdJpHyVZsF/elUmSpN2QAVlTx8Y1VH71fyj/6oeQVikd926KJ5wJXXvkXZkkSdqNGJA15aSrnqN863eoPHA7ycy5lE4JFI88BVradniuJEnSrrIHWVNTWqX6zINUlz9O4eBjKCx+I6W9DqT6+D1Un34AqpW8K5QkSdOUAVlT22A/1YfuoPrswxQPPYHiYSdQ2O8wKg//ivSlp/KuTpIkTUO2WKgxbFpH5a4fUf7lDVApUzr2DIonvh9m7pl3ZZIkaZoxIKuhpK8sp3zbd6jcfxtJ5yyaTgkUjz4NWtvzLk2SJE0Ttlio8aQp1ecepvriExQOejOF/Y+ktOgAqk/8hupT90OlnHeFkiSNSalUIi22Uq0mFAopSaWPctnfY3kzIKtxlQeoPvKrWn/yYW+leMhxFPY7lMojd5Iufzzv6iRJ2q5SqcTGgSYuueJOVq3tZd7sNi5c8hY6mzEk58wWCzW+ng1U7v4x5V98H/r7KL35dyme9EGS2fPzrkySpBFVKlUGaeGSK+5m1dpeAFat7eWSK+4mLbbmXJ2cQda0kb76EuWfRZJ93kDx0OMpnfwhqssfp/LIr6B3U97lSZLG4pZv5l3Ba1auJmwsF9lQLrGh3MSGwVJte7DExnK2nT3vrpT42/924uZwPGTV2l6qPZvgtsb9cxizt5+TdwWjMiBr2klfWEZ5xVMUDjyawgFHU1q4P9Un76P6xH1QGcy7PElSAxmsJpvD7cbBEuvrtjcM299TGTlWtRQqdJXKdDWVmd/Sz4Ed3XQ1lWkpb2Le7LatQvK82W0UBntwtf98GZA1PVUGqS67i+pzj1A89K0UD34LhdcdSuXRO0mfX5Z3dZKkHNWH3vVDs7sjBN6N5SZ6KsURX6M1C72dTWUWtvZzUFN3LQRnj86mMjNLg3Q1lWkupCO+RumpW7nwY2dwyTUPbulB/tgbSZbdNJEfX2NgQNb01ruJyr0/ofr0AxQOfxulo99OuviNVB76BemrL+VdnSRpnAxUkxFbGWo/m2rBN9vXWx099M5sKtNZKrNXaz9dWejtzILuzFLt2PZC784ov/oinctu4pKPnUS1qZ3CYA/Jspsov/riLr+2do0BWbuFdO1KKj//HtW9DqR46Fspve0DVFc8SeXhX0HPhrzLkySNYCj01gfejeWm2uzusP19o4Te9mKlFnBLZfZu66WzszbDOxSEh1ofukplmsYh9O6s8qsvwh3XAthWMYVMuYAcQvgw8EXgEODYGOM9dccuAD5B7Tv0qRjjj7P9ZwD/CBSBb8QYL53sutUY0hefoPzyMxRefxSFA99Eaf5iqk//lurj90J5IO/yJGna6x+ssKF7kA09ZTZ2D7K+Z5AN3YNs7Knt27Bi/83hd/TQu6WVYZ+2XmZ21loauoYF3s6cQq8a35QLyMBDwAeAf6nfGUI4FPgIcBiwCLg5hHBQdvh/A78LLAfuDiH8R4zxkckrWQ2lUqb6+D1Un3+E4iHHUzzwTRT2fQPVR++i+vwjkPofU0naGX0DlVrAzQLvxu4yG3oGa4/uQTb2lFmfheD+weqIr9HRWqSrvYmuJOV17b2bWxmG9/R2liqUDL2aYFMuIMcYHwUIIQw/dCZwbYyxH3gmhPAkcGx27MkY49PZeddmYw3I2r6+Hir3/ZTK0w9SPPxtFI86hcL+WX/yK8vzrk6ScpOmKf2D1dpMb3cWdHvKdTO9g1ngre0bKI8ceme0lehqb6Kzo8TihR21AJw9n9mRbbc30dleolTMbs3QwMu8afqYcgF5O/YC7qx7vjzbB/DCsP3HTVZRmgbWv0LljuupLnw9xcPeSumEM6m+/AyVh38Jm9blXZ0kjYs0TekbqLBhO7O7Q/s29JQZHCH0JtRCb2cWbvdc2EJnRxMz25vo6qiF3a72Um27rUSx6P3I1JhyCcghhJuBBSMcWhpjvGGU05IR9qWMfDfAEf/tJYRwHnAeQIyRuXPnjqFa7S7Sl56ivPJZCvsfQeGgYyid+hGqzzxE9bG7YbB/zK/j90rSZEnTlJ6+Mus29bNu0wDrN/WzbuMA6zb117Y3DbB+45ZjI830Jgl0dTQzc0YLs2a0ss/8mczqbGHmjGZmzdjyc9aMZro6mic89K6e0FfXVDKVf1/mEpBjjKe/htOWA/vUPd8bWJFtj7Z/+PteDlyePU1Xr/avoYapVmo3FXlhGYU3HEdh/zdS2Ofg2prKzz4M6cj/jFjP75WkXZGmKT39FTYOa23YMru7pbVhQ88g5cq2c0JJAp1ttZncrvYmXr+onc72mZtnd7uyGd+u9iZmtJUoFEaag6o3QGVggLVey6xxlMfvy0WLFo1pXCO1WPwHcE0I4R+oXaR3IHAXtZnlA0MIi4EXqV3I97HcqtT00N9L9be3UX3mQYqHn0jxiJMpLD6cysO/JF35XN7VSZokpVKJtNhKtZpQKKQklT7K5fJOv05tpreyVcitbZe37untHmRjb3nE0FtIyFoYaj28C+a0bgm6HVlrQ3sTnR1NzGgdS+iVNJopF5BDCO8H/gnYE7gxhHB/jPGdMcaHQwiR2sV3ZeCPYoyV7Jw/Bn5MbZm3b8YYH86pfE03G16l8sv/oLpgP4qHnUjp+PdQXfU8lYfugI1r8q5O0gQqlUpsHGjikivu3HKXsyVvobMZyuUy1aHQu83s7uCWPt+6Gd9KdYTQW0joai9tDr57zW3bana3s72Jro5a8O1oK1FIDL3SZEjS3XdJq3TFihE7MSbMOScvmdT30zhLChQWH07h4LdAUzPVZx+muuwuGOjbatg3b78in/okjatiywwuvKwWjofMm93Gf/vQkXz13+9lQ0+Z6giht1hI6Mxmc7e0M2wJwfUtDu2tRUPvcK5isft4+zmT/pZZi8UO/9JNuRlkacpKq1SffoDqC49ReMOxFPY7nMLeB1F97B6qzzwA1R33J0uNbs2nP513CROuL2niwfZ9eetFf7pVOAZYtbaXjlLKAS8/yoxKHzOqvcyo9NFZ6WNGpZcZ1T7aqgM7/u0LDGSPqWrOP/5j3iVIuTEgSztrsJ/qgz+n+sxDtf7kw0+ksN/hVB65g/SlZ/KuTtJrUCHh8dZF3NexmEfb9qZcKHF4zwDzZrdtM4PcuW4VH1xz53ZeTVKjMyBLr9WmtVTu/CHVefvW+pOPfRfVV5bzwqoe9pnXnnd1knYgBZY378F9HYt5oH0/uoutdFT6eEv3kxzd/Qyzrv01nz3rk1z63Uc39yB/9kOH0HvlZXmXLmmCGZClXZSuep7yKy9QeN1hFN5wLBdf9TAnHD6XM9+2NzM7mvIuT9Iwa4sd3NexmPs7FvNK00xKaYVDepZzdM/THNS7gmK2lH7/E6/SeuVl/PUHPgKdXbBxA71XXkbfE0/k/AkkTTQDsjQe0pTqsw9RXf44v/f5i7j1vlXc89gazjhuIae/aQHNTd5NSspTb9LMg+37cl/HYp5tnQ/A4r6VnLThEd7Y8zyt6eCI5/U98QR9f/elySxV0hRgQJbGU3mAD5+yLycfOY/v3/4CN/ziRX7xwCu8/6S9OebgOSRerS5NmjIFHm/b0ldcSYrsObied6y7j6O6n2V2pTvvEiVNUQZkaQLMn93KJ888kGXPb+C6257nGzc+za33reLDp+zD4oUz8i5PmrZS4IXmuVlf8evoyfqKj9v0BEd3P81eA2vGtMKEpN2bAVmaQG/Yt4ul/+UwfvnQam64YzmXXvMoxx2yB+8/aW9mdzbnXZ40bbxamsF97bW+4lebuihVyxzau5yju5/mwL6XNvcVS9JYGJClCVYoJLztiD1588FzuOmuFdx870p+88Ra3vmWBbzjLQtoaSrmXaLUkHoKzTzY/jru61jMcy3zSNKUxf0rOWXDQxze88KofcWStCMGZGmStLUUef9J+3DSEbX+5B/+agW/eLDWn3zsIXt4Ny1pDMoUeKxtL+7rWMyytr2oJEXmD6zjjLW/4cieZ5lV6cm7REnTgAFZmmRzZ7Zw3nsP4MnlG4m3Pc+//d9n+OlvVhJO3ZcD9urMuzxpykmB55r35L6OxTzY/jp6iy3MqPRy/MbHeVP30ywcXGtfsaRxZUCWcnLA3p189g8O5dePvMoPfrGcL1+7jDcfNJsPnLwPc2e25F2elLvVpc7aesXti1nT1ElTtcxhvS9wdPfTvL7vZfuKJU2YnQrIIYQ5wLExxptGOHYG8OsY49rxKk6a7gpJwlsPm8ubDprNf979Mj+++2V++9Q6Tn/zAn7vuIW0NtufrN1Ld6GFB7K+4hda9iRJU17f9zJv3/AAh/W8QEtazrtESbuBnZ1Bvgh4FdgmIANHA6cDf76rRUm7m5amIu89YS9OPHwuP/jFi9x010v88qFXOPPEvTnh8LkUCv4DsqavQQosa9ub+zoW81jbXlSTAgsG1vJ7a3/DkT3PMLPSm3eJknYzOxuQ3wOcMMqxy4E7MSBLr9mcrhbOedf+nHr0POKtz3PVT57l1vtXEk7Zl4P37cq7PGncVIHnWvbkvo79ebD9dfQVmuks93DixmUc3f00CwfX5V2ipN3YzgbkBTHG1aMcWwPM38V6JAGLF87gLz96CPc8tobv376cf7juMY46YBYfPHkf5s1uzbs86TV7pdRV6yvuWMza0gyaq4N1fcUrKdhXLGkK2NmAvDaEcHCM8bERjh0E+L/80jhJkoS3vGEPjnz9bG7+zcvc9OuX+OIVD3Ha0fN41/GLaG/1Gls1hk2FFh5o34/7OhazvGUuSVrlgL6X+d1193NY7ws0p5W8S5Skrezsb9jrga+FEH4/xri5KSyE0AZ8BfjueBYnCZqbCrzruEWceNhcbrjjRW6+dyW/euRV3nvCXpx0xJ4U7U/WFDSYFHm0bW9+07GYJ1oXUU0KLBxYw7vW3suR3c/SVbWvWNLUtbMB+XPAT4GnQwg3AS8BC4F3Ai8AXxjf8iQNmTmjmbPeuZhTjprHdbe9wL/f8hy33b+SD5+yL4ftNzPv8iSqwLMt87P1ivelv9BMV7mbt218lKO7n2bB4Pq8S5SkMdmpgBxj3BhCOAE4G3g7cAy1VS0+B1wVYxwY/xIl1dt3fgd/Fg7m/ifX8b3bX+Br33ucwxfP5EO/sw8L92jLuzzthlaWZm7uK15f6qC5Osgbe57n6O6nWdy/yr5iSQ1np5sYY4yDwDeyh6QcJEnC0QfO5vDFM7n1vpXceOdL/PW3HuLkI+fx3hP2Ykab/cmaWBsLrfy2o9ZXvKJ5DwpplQP6XuKMdfdxqH3Fkhrczt4o5JwdjYkxfvO1lyNpZzSVCrzjLQs5/rC5/J9fvsjPfruKux59lXe/dRGnHDWPUrGQd4maRgYGK9yfXWz3ZOtCqkmBvfpf5T1r7+aI7uforPblXaIkjYudnWb6+A6Op8AuBeQQwoeBLwKHULtr3z11xy4APgFUgE/FGH+c7X8W2JjtL8cYj9mVGqRG09XexB+cvt/m/uTrbnuB23/7Ch/8nX04Yv+ZJIkX8um1qVZTHnthI79+dDW/eXwt/XPfxqzyJk7e8AhHdz/NvPKGvEuUpHG3sz3Ip05UIXUeAj4A/Ev9zhDCocBHgMOARcDNIYSDYoxD/4536nbWaJZ2C3vNbefTHzyIh55Zz3W3vcDX///27j26yvrO9/g7CeFOSMgFCAkX5WIiIopCKorYWrSjtajl19Xa2zrT4xrPdJ1Oa8+Z1VpnPDMyHdeaM+0cO+OStjNWq/X8xKHWeq21gHIMeEdM8IIiNwFzQe7X7PPH3oy7CIORJM/O5v1aay94rvsb1i87H558n9/z6zc4bXQJ82bXUlM5MOny1ItsbNlNY1MrK5pb2bbzAP37FnHOpGHU//5XjN23BX83ISmfnXCjYgihLzAZeCvGeMLzIMcYmzPnPXLT54B7Y4z7gLdDCG8C04FnTvQ9pXxSUFDAGaeUUj+mhCUvv8eDz2zk5rte5fwzKrli5ihKBhYnXaJy1Ps797NidRuNTS1seG8PhYUFnD62hHmzK5hySil9iwtpe3hL0mVKUrfrbA9yCen2h3rSwXQB8DQwDtidmR/5ia4uMmMU6UdZH7Yhsw7SrR2PhxBSwO0xxgVHO0EI4VrgWoAYIxUVFd1Uqk5muTSuwpwqPnP+BOLv3+SxxnU891obn79oPJfNHENxn6Kky1MO2Lv/IMtf3cKSFzfxypstdKRgfM1Q/vSzY5k5ZSRDB/f7o/3bEqpTPS+pzzJ/FXzyyKWfl0fq7BXk24BhwAPAXOALwD+RntHivwDzgeMG5BDCE8CIo2y6Icb4wDEOO1oT5eG5g2bGGDeFEKqA34UQVscYlx65cyY4Hw7PqZYWvw3V9XJxXH3uE1XMmDSEhUvWc9ejr/HIM2u5elYNZ00osz/5JNTRkWL1uu0sb27lxTfa2Xegg/KSvlw6YyQz6soZMSw9XeCBvTto2bsj4WqVlFz8LFN+SWKMVVdXf6T9OhuQPw2cEmPcGUK4F9gK/HOM8VAI4V+Amz/KSWKMF3fyfSF9xbg2a7kG2JQ53+E/t4YQFpFuvfhQQJZOZiOGDeCbV06kae373LdkPbc/uIYJNUOYN7uWMcMHJV2eesD6rbtZ3tTCitVtvL/rAAP6FXHuaeU01Jdz6qjBFPqfJUkCOh+Q+8cYdwLEGNtDCDsP3yQXY+wIIXTnfRu/Ae4JIfwj6Zv0JgArQgiDgMLMQ0wGAXOAv+nGOqRerX7sUH4wuoRlr7zHA8s28sNfNtFwegVzzx9F6eC+SZenLta+Yz8rVreyvKmVjS17KCosYPK4oTTUl3PGKaUU9/F2O0k6UmcDckEIYRwftDt8aPlECwohXAncClQCD4UQXooxXhJjfDWEEIEm4CDw55kr18OBRZmb+voA98QYHz3ROqR8VlRYwKwzqzj3tGE83Pguv39hCy+83sYl00fy6Wkj6FtsaOrN9u4/xItvtNPY1Mpr67aTAsaNHMQXPzWGcyYN80EyknQcnf2UHASsOWJd9vIJP080xrgIWHSMbfNJ9zlnr3sLOPNE31c6GQ3o14erL6xl1pmV3L90A79ZtpGnVr7HVRfUcO5pw+xP7kUOdaRofmc7y5taePHNbRw42EHF0H78SUM1M+rLGV7WP+kSJanX6GxA/iKwJMa4uTuKkZSMytL+/NkV43l9/Xbi4vX8/OG3ePLFLYTZozmlenDS5ekYUqkU67fuprG5lWebW9m++yAD+xfxidPLaagr55Tqwf4nR5I+hs4G5L8FTg0hrCF9E9xiYGmMcV1XFyap502sLeH719TzTFMLv356I7f8qpnppw3jygtqGFbS7/gnUI9o276PFc1tNDa38G7rXooKCzjjlKE01FcwedxQ+4ol6QR19kl6EzM9v7Myr+8Cd4QQNpIOzEtijD/r+jIl9ZTCwgJmTq5k2sRhPLriXZ54fjMvvtnOnHNGMufcEfTv6/zJSdiz7xAvvNHG8qZWXl+/gxRwavVgrrl4DNMmDmOQfcWS1GU6/YkaY9wC3Jd5EUIoJf3wje8AXyI9J7KkXq5/3yLmnl/DBVMqWbR0Aw81buLpV95j7vk1NJxe7pRgPeDQoQ6a3tlOY1MrL69p58DBFFWl/bj8vGpm1JVTWWpfsSR1h04H5BBCATCVD64in0d6PuIIPNWl1UlKXHlJP75x+V0NuNEAABCvSURBVKlcdFYVcfF6fvHY2/zhpXR/8oSaIUmXl3dSqRTvbEnPV/zs6jZ27DnIoP59mDm5khl15YwbOci+YknqZp191PRvgbOB10g/YnoB8PUYo49akvLcqaOG8JdfquPZ1W0sWrqef/i/qzl7QhlXzarxSmYXaN2+j+VNrSxvbmVz2176FBUw5dRSGurKOX3cUPoU2VcsST2ls1eQJwH7gLdJT+/2puFYOnkUFhQwo66cs8aX8vhzm3lsxWZWvrWNT549nD+ZUc2AfvYnd8buvQd54fV2GptbeWND+qN0/KjBfPnTY5k2sYyB/e0rlqQkdPYmvQlH3KT3FyGECmAZ6faKp2OML3V9mZJySd/iIi7/xChmTq7kgWUbePzZzTzzagtXzBzF+ZMrKSy0BeBYDh7q4NW16fmKX16zjYOHUgwv68/nZo5iel05FUOdLUSSktaVN+n9gPTT77yEJJ0kyob05euXnsLsqcO5b/E67v7dOyx+cSvzZo+mbkxJ0uXljFQqxdrNu2hsauXZ1W3s2nuQIQP6cMGUShrqyxkz3L5iScolXXGT3vlAKfAc8K9dWp2kXmHsiEF89wun8cLr7dy/dD0/XvgaU04p5fMX1jJ82Mnbn9zy/gd9xVva91Lcp4AzTy2job6c+jElFNlXLEk5qbM36T1EetaKvsByYAnwE+CZGOPeri9PUm9RUFDAtEnDmHJqKb9/fguPrNjETb9YxUVTq7jsE9UMOkn6aXftPcjzr7XR2NTKmk07AZhYO4RLzh3B2RPLGNDv5Ph3kKTerLOf1E8B84FnY4wHuqEeSb1ccZ9CLp0xkvMmV/DAso08+cIWGpta+Ox5o5g1pTIvr5oePNTBqrfep7G5lVfeSvcVjxzWn7nnj2JGXblPIZSkXqazN+n9fXcVIim/lAwq5itzxjJ7ahX3LV7HvU+uY/FLW5k3u5bJ40qTLu+EpVIp3np3F8ubWnjutTZ27T3EkIF9uPDMKhrqy6mtGmhfsST1Uv6uT1K3qq0ayLfnTeLlNdu4f8l6bv33N6gfW8K8C0dTXTEg6fI6bWv7XpY3p/uK39u2j+I+hUwdX0pDfTl1Y4ZS5AwektTrGZAldbuCggKmji9j8rihLH5pK799ZhN/e+cqLphSxRXnVTN4YHHSJf6ndu453Ffcwlvv7qIAmDS6hMsaqjlrQhn9+zp5jyTlEwOypB7Tp6iQi6eNoKGunAef2cTSl7eyYnUrlzVUc9FZVTn1tLgDBzt45a1tNDa1surt9znUkaK6YgBXXVDD9Lpyyob0TbpESVI3MSBL6nGDBxbzxU+N4cKpVSxcvJ6FS9az9OWtXD2rljPHlybWu5tKpVizcSeNza08/1obu/cdomRQMRedVUVDfQU1lQPsK5akk4ABWVJiqssH8N+vnsiqt7excPF6bvvNm0yqHcK82aOprRrYY3Vsafugr7jl/X307VPIWRPS8xWfNrrEJwNK0knGgCwpcZPHlVI3uoSlK9/jwf+3ifl3vcrMMyq4YmYNQwd1T3/yzt0HeDYzX/HazbsoKIDTRpfw2fOqmTrevmJJOpkZkCXlhKKiQi46azjT68p5qHETf3hxK8+ubuMzM0Zy8bQRFPc58f7kAwc7WLkm01e89n06OlLUVA7g8xfWcu5pwygdbF+xJMmALCnHDOrfhzB7NBdOqeL+pev59dMbeWrle1w1q5ZpE8s63QPckUrx5oYdNDa18vzr7ezdf4jSwcVcfPZwGurLGVXZc60ckqTewYAsKScNH9af/zZ3As3vbGfhknX89LdreLJ6MOGiWsaOGHzc4ze37qGxuZUVza20bt9Pv+JCzp5Yxoy6cibV2lcsSTq2nAvIIYR5wE1AHTA9xvhcZn05sBA4F7gjxvjNrGOmAXcAA4CHgW/FGFM9W7mk7lA3poQbvnw6y1a18MCyDfzw7mYa6sv5wqdOYciQwXR0FFBYmKLg0F7atu/h2dWtNDa1sm7LbgoKoH7MUOaeX8OZ40vpV2xfsSTp+HIuIAOrgKuA249Yvxe4EZiceWW7DbgWaCQdkC8FHuneMiX1lMLCAi6YUsk5k4bxyPJNrHtvH7sOFPPD2xrZ2r6HqrIBfPuLZ/OLh95h9bp2aqsGMm92LeeeVt5tN/lJkvJXzgXkGGMzQAjhyPW7gKdDCOOz14cQRgIlMcZnMst3AnMxIEt5Z0C/Iq6aVUuqaBB/9dPlbG3fA8DW9j386Fcv8N1rzmbXjvd75SOsJUm5I+cC8scwCtiQtbwhs+5DQgjXkr7STIyRioqK7q9OJ50kx9VtZ/9pYu/dky679+b/CMeHbW3fQ/GuXSz70k3JFNXDrnvh54m8b1si76okJPVZ1pLIuyoJuZzDEgnIIYQngBFH2XRDjPGBTp7uaHfaHLX/OMa4AFhweJ+WFr8N1fUcV91vf/sOqsoG/FFIriobwP72HQlW1bMcZ+pujjF1tyTGWHV19UfaL5GAHGO8uAtPtwGoyVquATZ14fkl5ZiVty3k+hu+wf9e1PQfPcjXX1nPyvk/S7o0SVIe6PUtFjHGd0MIO0IIDcBy4KvArQmXJakbbVm5hpXzf8b1132evmVD2N++g5Xzf8aWlWuSLk2SlAdyLiCHEK4kHXArgYdCCC/FGC/JbFsLlAB9QwhzgTkxxibgOj6Y5u0RvEFPyntbVq7hd9fdknQZkqQ8lHMBOca4CFh0jG1jj7H+OT489ZskSZLUaYVJFyBJkiTlEgOyJEmSlMWALEmSJGUxIEuSJElZDMiSJElSFgOyJEmSlMWALEmSJGUxIEuSJElZDMiSJElSFgOyJEmSlMWALEmSJGUxIEuSJElZDMiSJElSFgOyJEmSlMWALEmSJGUxIEuSJElZDMiSJElSFgOyJEmSlMWALEmSJGUxIEuSJElZDMiSJElSFgOyJEmSlKVP0gVkCyHMA24C6oDpMcbnMuvLgYXAucAdMcZvZh2zGBgJ7MmsmhNj3NqDZUuSJCmP5FRABlYBVwG3H7F+L3AjMDnzOtI1h8O0JEmSdCJyKiDHGJsBQghHrt8FPB1CGJ9EXZIkSTp55FRAPgH/FkI4BNwP3BxjTB1tpxDCtcC1ADFGKioqerBEnSwcV+oJSY2ztkTeVUlIaoy1JPKuSkIu/7zs8YAcQngCGHGUTTfEGB/4GKe8Jsa4MYQwhHRA/gpw59F2jDEuABZkFlMtLX4bqus5rtQTHGfqbo4xdbckxlh1dfVH2q/HA3KM8eIuPt/GzJ87Qgj3ANM5RkCWJEmSjqdXT/MWQugTQqjI/L0YuJz0jX6SJEnSx5JTPcghhCuBW4FK4KEQwksxxksy29YCJUDfEMJcYA7wDvBYJhwXAU8AP02idkmSJOWHnArIMcZFwKJjbBt7jMOmdVtBkiRJOun06hYLSZIkqasZkCVJkqQsBmRJkiQpiwFZkiRJymJAliRJkrIYkCVJkqQsBmRJkiQpiwFZkiRJymJAliRJkrIYkCVJkqQsBmRJkiQpiwFZkiRJymJAliRJkrIYkCVJkqQsBmRJkiQpiwFZkiRJymJAliRJkrIYkCVJkqQsBmRJkiQpiwFZkiRJymJAliRJkrIYkCVJkqQsfZIu4EghhHnATUAdMD3G+Fxm/aeBvwf6AvuB/xFjfDKzbRpwBzAAeBj4Vowx1ePFS5IkqdfLxSvIq4CrgKVHrG8BPhtjPAP4GnBX1rbbgGuBCZnXpT1QpyRJkvJQzl1BjjE2A4QQjlz/Ytbiq0D/EEI/YBhQEmN8JnPcncBc4JEeKViSJEl5JecC8kd0NfBijHFfCGEUsCFr2wZg1NEOCiFcS/pKMzFGqquru73QbI+++XiPvp9OPt9addfxd5JOQPV99yVdgvLdV36QdAVSMgE5hPAEMOIom26IMT5wnGNPB24B5mRWFRxlt6P2H8cYFwALOlGqTlAI4bkY4zlJ16H85RhTT3Ccqbs5xnJLIgE5xnjxxzkuhFADLAK+GmNck1m9AajJ2q0G2HRiFUqSJOlklYs36R1VCKEUeAj4Xoxx2eH1McZ3gR0hhIYQQgHwVeA/vQotSZIkHUvO9SCHEK4EbgUqgYdCCC/FGC8BvgmMB24MIdyY2X1OjHErcB0fTPP2CN6gl0tsaVF3c4ypJzjO1N0cYzmkIJVyumBJkiTpsF7TYiFJkiT1BAOyJEmSlCXnepCVH0IItcCdpKfz6wAWxBj/KdmqlE9CCP1JP3GzH+nPsoUxxr9OtirloxBCEfAcsDHGeHnS9Si/hBDWAjuAQ8BBp3rLDV5BVnc5CFwfY6wDGoA/DyHUJ1yT8ss+4JMxxjOBqcClIYSGhGtSfvoW0Jx0EcprF8UYpxqOc4dXkNUtMtPvvZv5+44QQjPpJxw2JVqY8kaMMQXszCwWZ17edawulZl//zJgPvCdhMuR1EMMyOp2IYSxwFnA8oRLUZ7J/Or7edJTQP5zjNExpq72Y+B/AkOSLkR5KwU8HkJIAbdnnvqrhNlioW4VQhgM3A/8RYxxe9L1KL/EGA/FGKeSfoLm9BDC5KRrUv4IIVwObI0xPp90LcprM2OMZwOfId2OOCvpgmRAVjcKIRSTDsd3xxj/Pel6lL9ijNuAxcClCZei/DITuCJzE9W9wCdDCL9MtiTlmxjjpsyfW4FFwPRkKxLYYqFuknns98+B5hjjPyZdj/JPCKESOBBj3BZCGABcDNyScFnKIzHG7wHfAwghzAa+G2P8cqJFKa+EEAYBhZl7dQYBc4C/SbgsYUBW95kJfAV4JYTwUmbd92OMDydYk/LLSOAXmT7kQiDGGH+bcE2S1BnDgUUhBEhnsntijI8mW5LAR01LkiRJf8QeZEmSJCmLAVmSJEnKYkCWJEmSshiQJUmSpCwGZEmSJCmLAVmS9JGFEBaHEL6RdB2S1J0MyJIkSVIWA7IkSZKUxSfpSVKOCyGsBX4CfBUYAzwKfC3GuPcY+88Gfgn8C/AdYCdwQ4zx7sz2fsB8IAD9gEXAt2OMe0IIZcBdwAzSPyOWAX8WY9xwlPcZCTwG3Blj/Ieu+nolKWleQZak3iEAlwLjgCnA14+z/wigAhgFfA1YEEKYlNl2CzARmAqMz+zzV5lthcC/kQ7io4E9pMP5HxcTwlhgCfATw7GkfOMVZEnqHf5PjHETQAjhQdLh9nhujDHuA5aEEB5KHxpuBv4rMCXG2JY5398B9wDfizG2AvcfPkEIYT7whyPOWw/8ILP/r07w65KknGNAlqTeYXPW33cD1cfZvz3GuCtr+Z3MMZXAQOD5EMLhbQVAEUAIYSDwI9JXq8sy24eEEIpijIcyy9cAbwILP96XIkm5zRYLScpPZSGEQVnLo4FNQAvptonTY4ylmdfQGOPgzH7XA5OAGTHGEmBWZn1B1rluypznnhBCUXd+EZKUBAOyJOWv/xVC6BtCuAC4HLgvxtgB/BT4UQihCiCEMCqEcEnmmCGkA/S2EMIw4K+Pct4DwDxgEHBXCMGfJZLyih9qkpSfNgPtpK8a3016JorVmW1/SbpFojGEsB14gvRVY4AfAwNIXyFuJD1jxofEGPcDVwFVwL8akiXlk4JUKpV0DZKkLnR4mrcYY03StUhSb+T/+CVJkqQszmIhSb1QCOH7wPePsukp0vMcS5I+JlssJEmSpCy2WEiSJElZDMiSJElSFgOyJEmSlMWALEmSJGUxIEuSJElZ/j8hgAatM4rymQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### matplotlib figureの準備\n",
    "f, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7))#, sharex=True)\n",
    "\n",
    "### WAIC\n",
    "# ベースライン調整\n",
    "baseline_waic = 4\n",
    "tempDf_waic = df_waic.copy()\n",
    "tempDf_waic[\"WAIC\"] += baseline_waic\n",
    "# plot WAIC-values\n",
    "ax1.set_title('WAIC')\n",
    "sns.barplot(x = \"n_peak\", y = \"WAIC\", palette=\"rocket\",data=tempDf_waic, ax=ax1, bottom=-baseline_waic)\n",
    "sns.lineplot(x=[0,1,2,3], y = \"WAIC\", palette=\"rocket\",data=df_waic, ax=ax1, marker=\"o\")\n",
    "ax1.set_xlabel('n_peak')\n",
    "ax1.set_ylim([-3.3,-3.2])\n",
    "\n",
    "### WBIC\n",
    "# ベースライン調整\n",
    "baseline_wbic = 150\n",
    "tempDf_wbic = df_wbic.copy()\n",
    "tempDf_wbic[\"WBIC\"] += baseline_wbic\n",
    "# plot WBIC-values\n",
    "ax2.set_title(\"WBIC\")\n",
    "sns.barplot(x = \"n_peak\", y = \"WBIC\", palette=\"rocket\",data=tempDf_wbic, ax=ax2, bottom=-baseline_wbic)\n",
    "sns.lineplot(x=[0,1,2,3], y = \"WBIC\", palette=\"rocket\",data=df_wbic, ax=ax2, marker=\"o\")\n",
    "ax2.set_xlabel('n_peak')\n",
    "ax2.set_ylim([-120, -90])\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*[図]上段：候補となるモデル（ピーク本数）に対するWAIC値, 下段：〃に対するWBIC値*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. まとめ\n",
    "- WBICはピーク本数＝３で最小値を取った。\n",
    "- 一方で、WAICは本数が増えるごとに単調減少傾向にあった。\n",
    "- 真のピーク本数は３本であったから、WBICいずれも真の本数を選択できていたことがわかる。対してWAICは、複雑すぎるモデルを選択してしまっている。\n",
    "  - WAICは予測性能が最良のモデルを選ぶ指標であり、「真のピーク本数」を推定するものではない。\n",
    "  - WBICは一致性をもつ規準であり、実際に真のピーク本数を当てることができている。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:anaconda]",
   "language": "python",
   "name": "conda-env-anaconda-py"
  },
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}