Skip to content

Instantly share code, notes, and snippets.

@narrowlyapplicable

narrowlyapplicable/WAIC&WBIC with pystan-lowiter.ipynb

Created August 27, 2018 12:51
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save narrowlyapplicable/580f2d3610190b7bc04693466e409224 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save narrowlyapplicable/580f2d3610190b7bc04693466e409224 to your computer and use it in GitHub Desktop.
Download ZIP
[補足]WAIC&WBIC with pystan - iterを少なく取った場合
Display the source blob
Display the rendered blob
Raw
WAIC&WBIC with pystan-lowiter.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# WAIC&WBIC with pystan - 補足：iterが少ない場合"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- \"WAIC&WBIC with pystan2.ipynb\"の結果、重畳したLorenz関数型ピークの本数選択にはWBICが適しており、WAICでは選択を誤る可能性が示唆された。\n",
"- 参考として、MCMCのiterが少ない場合の結果も示しておく。\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=1000, warmup=400, 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=1000, warmup=400, 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.660633</td>\n",
" <td>-90.208390</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>3</td>\n",
" <td>-3.230468</td>\n",
" <td>-107.446927</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>4</td>\n",
" <td>-3.229345</td>\n",
" <td>-101.171823</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>5</td>\n",
" <td>-3.228514</td>\n",
" <td>-99.363180</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" n_peak WAIC WBIC\n",
"0 2 -2.660633 -90.208390\n",
"1 3 -3.230468 -107.446927\n",
"2 4 -3.229345 -101.171823\n",
"3 5 -3.228514 -99.363180"
]
},
"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": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAHwCAYAAAC7apkrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmYVNWBv/H3dhU7iCggtvsCChLFDbeIa+IS9+gZNRM1Zn7EmUySySQZF7JndDRmJsuYxDiJMZvLiYoYNbgrasQlLmgE2dxbRRRZZK3q+/ujClNAN90N3X2rq9/P83S6695T934bK/SX2+eeStI0RZIkSVJJXdYBJEmSpGpiQZYkSZIqWJAlSZKkChZkSZIkqYIFWZIkSapgQZYkSZIqWJAlSZKkChZkSapySZJ8NkmSVUmSDFhr+7T1bP/NWtueS5KkmCTJ7k0c/9tJksxea1s+SZIvJEnyeJIki5MkWZgkydNJkkxIkmRQe35/klRtLMiSVP3uAfLAIas3JEkyGNgNeLOJ7aOBeyu2HQgMBX4FjG/pZEmS9ABuBy4GInA4sAcwAdgfOHtjvyFJqmb5rANIktYvTdNXkiSZAxwB3FbefDjwPDClie0JFQUZ+BzwB+B64K4kSf4jTdOl6znlF4GPAQelafpoxfaXgTu8giyp1nkFWZK6hnspFeHVDgfuK3+svX1GmqZvAJTL7GnAb9I0fRx4AwgtnOvTwH1rleMPpWm6YIO+A0nqIizIktQ13AuMTpJkaPnx4cD9wIPAqLW231PxvLOAF9M0fbb8+De0PM1iBPBCu6SWpC7IgixJXcN95c9HJEmyNbAj8GCapu8B0yq2D2fNgjyeUile7XfA2CRJRq/nXAmQtltySepiLMiS1AWkaTofeJbSdIojgKfSNF1Y3n1/xfYi8ABAkiQfBUYB/50kSSFJkgLwGpBj/VeRX6R0A6AkdUsWZEnqOlbPQ149/3i1+yu2P1lRnD8H3E1pBYoxFR9fAj6dJEmfZs7ze+DwJEkOaGqnN+lJqnUWZEnqOu4FtgdOZs2CPAXYprz9XoAkSTYDTgV+l6bp85UflJZ760Xp5r2m/Lh8nDuTJPlqkiT7JEmyXZIkRydJcgulec2SVLMsyJLUdUwBVgG9gYdXb0zTdBHwV2AAf59/fDalucST1j5ImqYfUFrnuMlpFmmargKOAb4BnE7pRsDngP8CHmfNOc2SVHOSNPU+DEmSJGk1ryBLkiRJFarqnfRCCJcDxwMrgTnAZ2KM7zcx7mVgMaW7tQsxxn06M6ckSZJqV7VdQb4bGB1j3B2YCVy4nrGHxRjHWI4lSZLUnqrqCnKM8a6Kh1Mp3YEtSZIkdZqqKshrORe4oZl9KXBXCCEFfhFjvKq5g4QQxlO+UzvGuHe7p5QkSVJXkrQ4oLNXsQgh3AMMa2LXhBjjpPKYCcA+wCkxxnUChhDqY4wNIYShlKZlfCHGOKUVp08bGho2Ir0kSZK6qvr6eqjGgtySEMLZwHnAETHGpa0Y/21gSYzxB604vAVZkiSpm2ptQa6qm/RCCEcD5wMnNFeOQwj9QggDVn8NfBx4vvNSSpIkqZZVVUEGrqD0TlB3hxCeCSFcCaUpFSGEO8pjtgAeDiE8S+kdnW6PMU7OJq4kSZJqTdVNsehgTrGQJEnqprrkFAtJkiQpaxZkSZIkqYIFWZIkSapgQZYkSZIqWJAlSZKkChZkSZIkqYIFWZIkSapgQZYkSZIqWJAlSZKkChZkSZIkqYIFWZIkSapgQZYkSZIqWJAlSZKkChZkSZIkqUI+6wCVQgiXA8cDK4E5wGdijO83MW5T4JfAaCAFzo0xPtqZWSVJklSbqu0K8t3A6Bjj7sBM4MJmxv0YmBxj3BXYA5jeSfkkSZJU46rqCnKM8a6Kh1OBU9ceE0LYBBgHnFN+zkpKV5wlSZKkjVZVBXkt5wI3NLF9R+Ad4NchhD2AvwJfijF+0JnhJEmSVJs6vSCHEO4BhjWxa0KMcVJ5zASgAPyhiXF5YC/gCzHGx0IIPwYuAL7RzPnGA+MBYowMHjx4478JSZIk1axOL8gxxiPXtz+EcDZwHHBEjDFtYsjrwOsxxsfKj2+kVJCbO99VwFXlh+n8+fPbHlqSJEldXn19favGVdVNeiGEo4HzgRNijEubGhNjfAt4LYSwS3nTEcALnRRRkiRJNa6qCjJwBTAAuDuE8EwI4UqAEEJ9COGOinFfAP4QQpgGjAEu6fyokiRJqkVJmjY1i6FmpQ0NDVlnkCRJUgbKUyySlsZV2xVkSZIkKVMWZEmSJKmCBVmSJEmqYEGWJEmSKliQJUmSpAoWZEmSJKlCp7+TnmpbPp8nzfWmsTGhri4lKS6nUChkHUuSJKnVLMhqN/l8nsUre3DJNVOZt2AZQwf14aJz9mVATyzJkiSpy3CKhdpNmuvNJdc8wbwFywCYt2AZl1zzBGmud8bJJEmSWs8ryB3s3HHnZB2h01w+8coPy/Fq8xYsY17DfL528nkZpeo8V0+5JusIkiSpHXgFWe1m8XsLGTqozxrbhg7qw+KFSzJKJEmS1HYWZLWb+L+/5gsn7fphSR46qA9fDGO46bF5JJvXZ5xOkiSpdZxioXYz45kZXPPtH/DZL3yGAZsNZPF7C/n1d3/I7L4jyB14AsVn7id97cWsY0qSJK2XBVntasYzM/juZ89fc2OP58ntezT5vY6k2G8gjTMezyacJElSK1RVQQ4hXA4cD6wE5gCfiTG+v9aYXYAbKjbtCHwzxvijTguqtlm1guKjf4I9DiG3y74k/QZSfPo+aCxmnUySJGkd1TYH+W5gdIxxd2AmcOHaA2KML8YYx8QYxwB7A0uBiZ0bU22WNlJ85n6KLzxK3dYjyB14IvR0+TdJklR9quoKcozxroqHU4FTW3jKEcCcGOMrHZdK7alx1lOkHywkt9eR5MedSmHqbbDk/ZafKEmS1Emq7QpypXOBP7cw5nTguk7IonaUNsyh+MgtkOtB/uBPkgzeKutIkiRJH+r0K8ghhHuAYU3smhBjnFQeMwEoAH9Yz3F6AifQxDSMtcaNB8YDxBgZPHjwBiZXe0oXvE3hoRvJ73ccuQOOp/jsA6Svzsg61kbxtSVJUm1I0jTNOsMaQghnA+cBR8QYl65n3InA52OMH2/D4dOGhoaNjdgm3emd9DZIvie5fY+mbug2FGf+lcbpU7NOtMF8Jz1JkqpbfX09QNLSuKqaYhFCOBo4HzhhfeW47AycXtH1FVZSnHobjS//jdyIvcntcxTU5bJOJUmSurGqKsjAFcAA4O4QwjMhhCsBQgj1IYQ7Vg8KIfQFPgbcnE1Mtau0keKzD1B8/hGS+p3IHXQS9OzT8vMkSZI6QNVNsehgTrGocsmWO5Lb60hYsYzCY7fB4gVZR2o1p1hIklTduuQUCyl9c255hYtcaYWLIVtnHUmSJHUzFmRVnfT9eRSm3AjLlpDb/ziSbUdmHUmSJHUjFmRVp2VLKDx0E+k7r5Pf83DqRh2QdSJJktRNWJBVvQqrKD52O8WXniM3fC9y+x4Nuap680dJklSDLMiqbmlK47QpFJ97uHQD30EnQa++WaeSJEk1zIKsLqFx7rMUH/8zyYDNyI87FQZslnUkSZJUoyzI6jLSt16i8PBEqKsrr3CxTdaRJElSDbIgq2tZ+A6FB2+EpYvI7X8cddvvlnUiSZJUYyzI6nqWL6Hw8M2k814lt8eh1O12EK1Y81uSJKlVLMjqmgqrKD5+B8W508jtPIbcWFe4kCRJ7cOCrK4rTWl87iGKzz1EMmx7ch89GXq7woUkSdo4FmR1eY1zp1F87A6S/oNKK1xssnnWkSRJUhdmQVZNSN9+hcJDNwMJ+YNPIRm6XdaRJElSF2VBVu1YNJ/ClBthyUJy+x9L3Q6js04kSZK6IAuyasvyDyg8PJH0rVfI7X4IdaM/iitcSJKktqi62/5DCJcDxwMrgTnAZ2KM7zcx7svAPwEp8Fx53PLOzKoqVVxF8fE/k44+kNxOY0j6bULxybuhuCrrZJIkqQuoxivIdwOjY4y7AzOBC9ceEELYCvgisE+McTSQA07v1JSqcimNzz9C8dkHSbbYjvxHT4be/bIOJUmSuoCqu4IcY7yr4uFU4NRmhuaBPiGEVUBfoKGjs6nraXz5edKli8jtexT5cadSeOx2WDg/61iSJKmKVV1BXsu5wA1rb4wxvhFC+AHwKrAMuGutYv2hEMJ4YHz5eQwePLgD46oapfNepfDQzeT3+wT5j55C8cm7SN9+ud3P42tLkqTakKRp2uknDSHcAwxrYteEGOOk8pgJwD7AKTHGdK3nDwJuAv4BeB/4I3BjjPH3LZw6bWjo3AvN5447p1PPp/Xo1Zfc/p8gGTiYxucfoXHutHY9/NVTrmnX40mSpPZVX18Prbh7P5MryDHGI9e3P4RwNnAccMTa5bjsSOClGOM75fE3AwcCLRVkdWcrllJ8eCK5vT9G7iMHQ7+BND7/MGTwj0RJklS9qm6KRQjhaOB84JAY49Jmhr0K7B9C6EtpisURwJOdFFFdWbFQWuFi1IHkhu9J0m8gxSfvhIIrXEiSpJJqXMXiCmAAcHcI4ZkQwpUAIYT6EMIdADHGx4AbgacoLfFWB1yVUV51QY0v/IXiMw+QDNmG/EdPgd79s44kSZKqRCZzkDPkHGStIRmyDbl9j4JigcLU22HhOxt8LOcgS5JU3Vo7B7karyBLnSZ95zUKD90MjUXyHz2ZZNgOWUeSJEkZsyBLi9+jMOUm0sXvkRt7DHU7jck6kSRJypAFWYLSCheP3EL65hxyow+ibvdDIGnxNzCSJKkGVd0qFlJmigWKT9xJOmohueF7k/TbhOITd0JhZdbJJElSJ2rTFeQQwmblZdia2nd0+Q08pC6t8YWpFJ6+j2TwVuQPPgX6DMg6kiRJ6kRtnWLxdWDvZvbtCUzYuDhSdUhfnU7x0dugT3/y404l2XRo1pEkSVInaWtBPg74RTP7rgJO3Lg4UvVI579O4aGboFggd9BJJFvumHUkSZLUCdpakIfFGOc3s+89YIuNzCNVl8ULKEy5kXTRu+THHkPdzntmnUiSJHWwthbkBSGEXZrZNwJ4fyPzSNVn5TKKj9xC4+uzyO12ILk9DoXEBWAkSapVbf0pPxH4SQihT+XG8uMfUnr7Z6n2NBYp/vUuijOfpG773cjtfxzke2adSpIkdYC2LvP2DeA+YG4IYTLwJrAlcBTwGvCt9o0nVZfG6Y+RLllIbsyh5Md9ksLU22Dp4qxjSZKkdtSmK8gxxsXAgZSKcm9gn/LnbwAHl/dLNS19bQbFv9wKvfqWVrgY5NR7SZJqSZKmadYZOlPa0NDQqSc8d9w5nXo+daL+m5Lf/zjo3Y/iU/fyq+u/nnUiSZK0HvX19QAtvlVum6ZYhBDObWlMjPHqthxT6rKWvE9hyo3kxh5Dft+juOOxBo4ZuyWJb1EtSVKX1tY5yJ9uYX8KbHBBDiFcDhwPrATmAJ+JMa6zMkYI4UvA/6P0L4D/izH+aEPPKW2UlctL0y3GHM6kh+GdBSv41Me2I59zlQtJkrqqNhXkGONhHRWk7G7gwhhjIYRwGXAhcH7lgBDCaErleCylIj05hHB7jHFWB2eTmtZYpPjU3Rx/2qHcPrWBdxet4HMn7Ey/3m3996ckSaoG7fITPISwGXAGcHaMceyGHifGeFfFw6nAqU0MGwlMjTEuLZ/7QeBk4Psbel6pPZxw0FYMHdSL3975MpddO50vnDKcIZv2zjqWJElqow0uyCGEPPAJ4GzgWOAN4Mp2ygVwLnBDE9ufBy4OIWwOLCuf+8n15BwPjAeIMTJ48OB2jCj93eDBgzlu3GB22GYo3//dU3z/+hmc/497s+v2g7KOJkmS2qDNq1iEEPamVIrPAHKU3jzkk8CIGOO8Vjz/HmBYE7smxBgnlcdMoLSE3CkxxnUChhA+C3weWAK8ACyLMX65FfFdxUId5uop13z49dsLlnPFxJm8t2glZx+1A2NHbp5dMEmSBHTcKhbPAzsCdwCfA26LMa4MIRzb2mPEGI9s4RxnA8cBRzRVjsvH+BXwq/L4S4DXW3t+qTNsMag3558xiitvnc2v7pjLO++v4Nj9XeFCkqSuoK232vcFipSmNiwFVrVnmBDC0ZRuyjth9RzjZsYNLX/eFjgFuK49c0jtoX+fPF/65Aj2G7k5t/7lDa6Z/BKrCo1Zx5IkSS1o6zvp7Uhp3vFKSvOD3woh/C+ld9Nrj3ccuQIYANwdQngmhHAlQAihPoRwR8W4m0IILwB/Aj4fY1zQDueW2l2PfB2fOWYHTjhwK6a+8C4/vmkmS5YVso4lSZLWY4PfSS+E0IfS3OOzgMOBF4Gfxhh/1n7x2p1zkNVhKucgN+Xx6e/ymztfYrNNevKvJ49gi0GucCFJUmdq7RzkNl1BDiXDAGKMy2KMv48xfhzYDvg98IUNyCp1C2NHbs6XT9uFD5YVuezaF5j1+uKsI0mSpCa06QpyCGEmsBOld7mbAjwIPBhjfLVj4rU7ryCrw7R0BXm1d95fzv/ePIv5C1dw1lHbs/8olx6UJKkztPYK8oYs87YFMK78cTAwmtIayFMoleVftjVsJ7Igq8O0tiADfLC8wJW3zmbma4s57oB6jjug3hUuJEnqYB2yzBtAjPFt4I/lD0IIm1J6I45/B84EqrkgS1WhX+/SChe/v/sVbnu0gXkLlnPWUTvQI9/WhWUkSVJ7a3NBDiEkwBj+fhX5QKABiMBD7ZpOqmH5XB1nH7U9Qwf1YtLDb/DuopX8y4k7079vj6yjSZLUrbX1jUJuA/aitGLFw8BVwDkxRu82kjZAkiQcu189Qwb24prJL3HpddP515OHM2yzPllHkySp22rr73N3AVYAL1G6UW+25VjaePvuujn/HnZl+Yoil103nRdfW5R1JEmSuq32uElvMPAIpekVD8cYn2nvkO3Im/TUYdpyk15z3nl/OT+dOIt576/g0x/fngN2c4ULSZLaSxY36X0dGALk2npMSSVDNu3Nf5wxkl/8aQ7XTH6Jtxcs54SDtqLOFS4kSeo07XGT3keBTYEngavbNZ3UDfXtneeLpwznD/e8wp8fe5N33l/B2UftQM8ernAhSVJnaOtNerdTWrWiJ/AYpTcKuQJ4NMa4vP3jSd1TLlfHpz++PVsM6s3ND73Oe4tW8C8nDWeAK1xIktTh2noF+SHgYuCJGOOqDsgjqSxJEo4auyVDNu3F1X+ey6XXlla42HJzV7iQJKkjtfkmvS7Om/TUYdrjJr3mvPTmEn52yyxWFVM+d/zOjNxukw47lyRJtaq1N+k5qVHqAnbYsj8XfGoUg/r35Cc3z+Th597JOpIkSTWrzTfpdbQQwveAE4FGYB6lNyJZ57JvCOFsSitnAPxnjPE3nZdS6nybb9KL/zhjV6760xx+d9fLzFuwnJMO3toVLiRJamfVeAX58hjj7jHGMcBtwDfXHhBC2Az4FrAfMBb4VghhUOfGlDpfn155/vXk4YzbfQh3PvEW/3fbHFauasw6liRJNaXqCnKMsfItxPoBTU2SPgq4O8b4XoxxAXA3cHRn5JOylsvVceaR23HqIdvw9MwF/HecwaIPvGdWkqT2UnVTLABCCBcDZwELgcOaGLIV8FrF49fL25o61nhKb2RCjJHBg31nMnWMzn5tnXH0EHbcZgg/uuEZvn/9i1x0zt5su8WATs0gSVItyqQghxDuAYY1sWtCjHFSjHECMCGEcCHwr5SmU1RqatJlk8txxBivAq5aPWb+/PkbmFpavyxeWzttkeMrYVd+esssLvzZo3zu+J0Ytf3ATs8hSVJXUF7FokWZFOQY45GtHHotcDvrFuTXgUMrHm8NPLDRwaQuaPth/bjwzJFcMXEW/3vzTM44cjvG7T4061iSJHVZVTcHOYQwvOLhCcCMJobdCXw8hDCofHPex8vbpG5ps0168bXTRzJyu4H84e5XuPHB12jsXmucS5LUbqquIAOXhhCeDyFMo1R8vwQQQtgnhPBLgBjje8D3gCfKH98tb5O6rT69cnz+5OEcssdQ7n7yLX5x62xWripmHUuSpC7Hd9LrYL6TXvfRke+k1xZpmnLf02/zx/tfY9st+vL5k4YzsH/PrGNJkpQ530lP6qaSJOGIvYbxzyfuzJvvLue/rp3OG+8szTqWJEldhgVZqlF77DyIr52+K2ma8v3rp/P8SwuzjiRJUpdgQZZq2LZb9OOCM0cxZGBvrpg4kweemZd1JEmSqp4FWapxgwb05Kun78roHQZy3b2vEO9/lcbGbnXvgSRJbWJBlrqB3j1z/MuJwzl8z6Hc+9TbXHnrbJavdIULSZKaYkGWuom6uoR/OHw7Tj98W6bNfZ8f3DCDBYtXZh1LkqSqY0GWupnD9tyCz580nHkLlnPptS/w2jxXuJAkqZIFWeqGPrLjpnzt9JEkwOXXT+e5ue9nHUmSpKphQZa6qW2G9uWCT41ii0G9+ekts7jvqbezjiRJUlWwIEvd2Kb9Sytc7L7jptxw/6tcf98rrnAhSer2LMhSN9erR47zTtiZI/fegvufnsfPJs1yhQtJUrdmQZZEXV3CaYduy5lHbMffXlrI5ddPd4ULSVK3ZUGW9KFDxgzl8yePYP7CFfzXH17g1bc/yDqSJEmdzoIsaQ2jdxjIf5w+klxdwuXXz+DZ2QuyjiRJUqfKZx1gbSGE7wEnAo3APOCcGGNDE+MmA/sDD8cYj+vclFJt22pIaYWLn06cyc8nzebUQ7fhiL22IEmSrKNJktThqvEK8uUxxt1jjGOA24BvNjcO+HTnxZK6l4H9evDVf9iVMcMH8ccHXuO6e1+h6AoXkqRuoOoKcoxxUcXDfkCTP5FjjPcCizsllNRN9eyRY/zxO/HxfYbx4LPv8NOJs1i2whUuJEm1reqmWACEEC4GzgIWAodt5LHGA+MBYowMHjx44wNKTajl19bnPjmEHbd5jasm/Y3/uXEWF529N0M27ZN1LEmSOkSSpp3/K9MQwj3AsCZ2TYgxTqoYdyHQO8b4rWaOcyjw1TbMQU4bGtaZztyhzh13TqeeT9m5eso1WUfocNNfWciVt86hZ486Pn/ScLYf1i/rSJIktVp9fT1AizfUZHIFOcZ4ZCuHXgvcDjRZkCV1rpHbDeT8M0ZyxcSZ/OCGGXz22B3Zc/igrGNJktSuqm4OcghheMXDE4AZWWWRtK76wX244MxRbD2kD7+4dTZ3PfEmWfwmSpKkjlKNc5AvDSHsQmmZt1eA8wBCCPsA58UY/6n8+CFgV6B/COF14LMxxjszyix1K5v068G/n7Yr10yey01TXmfe+ys44/BtyeWq7t/ckiS1WSZzkDPkHGR1mO4wB3ltjWnKpIffYPLjbzJyu0343PE70adXNf67W5Kk1s9B9nKPpA1WlyScfPDWnPXx7XnxtcVcdt105i9ckXUsSZI2igVZ0kY76CND+NInR7BwySouvfYFXnpzSdaRJEnaYBZkSe1i12034fwzRtKrR47/jjP468z3so4kSdIGsSBLajfDNu/DBWeOZNuh/bjqT3OY/JgrXEiSuh4LsqR2NaBvD7582i7su8tmTHz4dX5318sUi41Zx5IkqdW83VxSu+uRr+PcT+zIkEG9uGPqm8xftILPHb8z/Xr7V44kqfp5BVlSh6hLEk48aGvOOXoHZr++hO9fN5133l+edSxJklpkQZbUoQ7YbTD/duoIFi1dxWXXTmdOgytcSJKqmwVZUocbsc0mXHDGKPr0yvE/cQZPzHg360iSJDXLgiypU2yxWW/OP2Mk2w/rxy9vn8sdUxtc4UKSVJUsyJI6Tf++Pfi3U3dh7MjNmPTIG/zmzpcouMKFJKnKeEu5pE7VI1/HucfsyNBNe3Pbow28u3Al552wM/36+NeRJKk6eAVZUqdLkoTjD9yKc4/dkblvLuGy66Yzb4ErXEiSqoMFWVJm9hu5Of926i4sWVbg0munM/v1xVlHkiSp+qZYhBC+B5wINALzgHNijA1rjRkD/BzYBCgCF8cYb+jsrJI23vCtB3DBmSO5YuIsfnjji5x11A7sN3LzrGNJkrqxaryCfHmMcfcY4xjgNuCbTYxZCpwVY9wNOBr4UQhh084MKan9DB1UWuFixy37c/Udc/nTX95whQtJUmaqriDHGBdVPOwHrPNTMsY4M8Y4q/x1A6UrzUM6J6GkjtCvT54vnTqC/Udtzm2PNvDrP7/EqoIrXEiSOl/VTbEACCFcDJwFLAQOa2HsWKAnMKeZ/eOB8QAxRgYPHty+YaUyX1vt46v/OISbHpjDdXfNYtGyRs7/x70Y0K9n1rEkSd1IksWvMUMI9wDDmtg1IcY4qWLchUDvGOO3mjnOlsADwNkxxqmtOHXa0NDQ8qh2dO64czr1fMrO1VOuyTpCTXlixrtcM/klNhvQk389eQRbbNY760iSpC6uvr4eIGlpXCZXkGOMR7Zy6LXA7cA6BTmEsEl539dbWY4ldSH77ro5mw3oyc8mzeay617gvBN2ZsQ2m2QdS5LUDVTdHOQQwvCKhycAM5oY0xOYCPw2xvjHzsomqXPttFVphYsBfXvwoxtn8ujf5mcdSZLUDVTjHORLQwi7UFrm7RXgPIAQwj7AeTHGfwICMA7YPIRwTvl558QYn8kgr6QONGTT3vzHGSP5xZ9mc83kl0iSHAeN2YbGNKGuLiUpLqdQKGQdUzUmn8+T5nrT2OjrTB3H11n1ymQOcoacg6wO4xzkjlUoNnL/M+8yZtd6fnLD08xbsIyhg/pw0Tn7MqDnKn+oqN3k83kWr+zBJdc84eusylV2mDXaTNrkl6uf1OS+5upQ2tyBN/L5PXr0YGXam//6ja+zztTaOcgW5A5mQe4+LMgdL9ezPxddOZV5C5Z9uG3ooD58+fQ9uXXKTFb/8Gnqr7XV29ImfnKmzT2GddZjbupnXdrM2DZlaWZ8ukbcNU+43iyt+D7bmr3ZMtLE2Mrx68tJU99fE8dc979R899fk1ma+ENrbvy/nDqGK2+ets7r7LxTdueK+PTah2m2cK2TobnS1oqSt/bD5kpXa0tis89vRclbz/9RhB8pAAAgAElEQVQj2v789f3Zra/k1oCLzhnLLyc9t87r7JJ/3p/iiiUZJqttVX2TnqT29cePn591hE7xiev/c40fJgDzFiwjAZ5+8pU1B5d/oibr+9Fa8cN5nb8tm/hJ//cxzf/kbvJ8TbTBpDX7mjjH5iO3aSYTJMmaz0ya+BGweszqXUn5f5Ik+fuxKp63emsCrHrxxfXmbOp7b+rP7MNzr/d7T5s8R7PnWetYrX5eE+cb2GNMk6+zgT1SdnhzesXW9fw3a+G8a2Ze37HWfG5z52n2HGscqvXH6n3YYeuOqXhBrbm9+XMmzQxs6rX54faXnl0jS9MZm36QNDdmfX+ObXz+OtFb9fzK7aVjbTGgrsnXWePSJfDA1c08u8YccW7WCZplQZbUZaxcsJihg/qsc8UlefNtPjrlzxkm6zynfe+yzM793pd+ltm5O1O/hUc0+Trrt/AdTn3v0QyTdZ7NDj07u5Mvvyu7c3eiXGPTf5/VrVpKMcNcKqm6VSwkqTnTfn4jXzl5FEMH9QFKP0y+cvIopv38xoyTqZYsu/l6Ljh15BqvswtOHcmym6/POJlqSTLjIS468yNrvM4uOvMjJDMeyjiZwCvIkrqQt6fNYdrFv+Qr/3wqPQcNYOWCxUy7+Je8Pa3JN9KUNsjyWbPo/duf891TTocBm8DiRSz77c9ZPmtW1tFUQwrvvsGAGZO55MyDaezRl7pVS0lmTKbw7htZRxMWZEldzNvT5nD3P2c3zUDdw/JZs1h+2feyjqEaV3j3DXik9JsJp1VUF6dYSJIkSRUsyJIkSVIFC7IkSZJUwYIsSZIkVbAgS5IkSRUsyJIkSVIFC7IkSZJUwYIsSZIkVai6NwoJIXwPOBFoBOYB58QYG9Yasx1wM5ADegD/G2O8srOzSpIkqfZU4xXky2OMu8cYxwC3Ad9sYsybwIHlMfsBF4QQ6jszpCRJkmpT1V1BjjEuqnjYD0ibGLOy4mEvqrPoS5IkqQtK0nSd/pm5EMLFwFnAQuCwGOM7TYzZBrgd2Bn4Wozxp80cazwwHiDGuHeHhZYkSVJXkLQ4IIuCHEK4BxjWxK4JMcZJFeMuBHrHGL+1nmPVA7cAx8cY3273sNogIYQnY4z7ZJ1DtcvXmDqDrzN1Bl9n1SeTKRYxxiNbOfRaSleJmy3IMcaGEMLfgIOBG9shniRJkrqxqpu7G0IYXvHwBGBGE2O2DiH0KX89CDgIeLFzEkqSJKmWVd1NesClIYRdKC3z9gpwHkAIYR/gvBjjPwEjgf8OIaSU5pH8IMb4XFaB1aSrsg6gmudrTJ3B15k6g6+zKlOVN+lJkiRJWam6KRaSJElSlizIkiRJUoVqnIOsLqq8NvVvKS3h1whcFWP8cbapVGtCCL2BKZTeJCgP3Li+pSCljRFCyAFPAm/EGI/LOo9qTwjhZWAxUAQKLvdWHbyCrPZUAL4SYxwJ7A98PoQwKuNMqj0rgMNjjHsAY4CjQwj7Z5xJtetLwPSsQ6jmHRZjHGM5rh5eQVa7iTG+CbxZ/npxCGE6sBXwQqbBVFNijCmwpPywR/nDu43V7kIIWwOfAC4G/j3jOJI6kQVZHSKEsD2wJ/BYxlFUg8q/9v4rpbea/2mM0deZOsKPgP8ABmQdRDUtBe4qL137ixijS75VAadYqN2FEPoDNwH/FmNclHUe1Z4YYzHGOAbYGhgbQhiddSbVlhDCccC8GONfs86imndQjHEv4BhKUxPHZR1IFmS1sxBCD0rl+A8xxpuzzqPaFmN8H3gAODrjKKo9BwEnlG+guh44PITw+2wjqRbFGBvKn+cBE4Gx2SYSOMVC7SiEkAC/AqbHGP8n6zyqTSGEIcCqGOP75becPxK4LONYqjExxguBCwFCCIcCX40x/mOmoVRzQgj9gLryfTv9gI8D3804lrAgq30dBHwaeC6E8Ex520UxxjsyzKTasyXwm/I85DogxhhvyziTJG2ILYCJIQQodbJrY4yTs40k8K2mJUmSpDU4B1mSJEmqYEGWJEmSKliQJUmSpAoWZEmSJKmCBVmSJEmqYEGWJLUohPBACOGfss4hSZ3BgixJkiRVsCBLkiRJFXwnPUmqUiGEl4ErgLOA7YDJwNkxxuXNjD8U+D3wM+DfgSXAhBjjH8r7ewEXAwHoBUwEvhxjXBZCGAT8DtiP0s+GR4DzYoyvN3GeLYE7gd/GGH/QXt+vJFULryBLUnULwNHADsDuwDktjB8GDAa2As4Grgoh7FLedxkwAhgD7Fwe883yvjrg15SK+LbAMkrlfM0wIWwPPAhcYTmWVKu8gixJ1e0nMcYGgBDCnyiV25Z8I8a4AngwhHB76anhP4H/B+weY3yvfLxLgGuBC2OM7wI3rT5ACOFi4P61jjsK+Hp5/HUb+X1JUtWyIEtSdXur4uulQH0L4xfEGD+oePxK+TlDgL7AX0MIq/clQA4ghNAX+CGlq9WDyvsHhBByMcZi+fGngNnAjRv2rUhS1+AUC0mqLYNCCP0qHm8LNADzKU2b2C3GuGn5Y2CMsX953FeAXYD9YoybAOPK25OKY327fJxrQwi5jvwmJClLFmRJqj3fCSH0DCEcDBwH/DHG2Aj8H/DDEMJQgBDCViGEo8rPGUCpQL8fQtgM+FYTx10FnAb0A34XQvBniKSa5F9uklRb3gIWULpq/AdKK1HMKO87n9IUiakhhEXAPZSuGgP8COhD6QrxVEorZqwjxrgSOAUYClxtSZZUi5I0TbPOIElqB6uXeYsxbp11FknqyvyXvyRJklTBVSwkqQsJIVwEXNTErocorXMsSdpITrGQJEmSKjjFQpIkSapgQZYkSZIqWJAlSZKkChZkSZIkqYIFWZIkSapgQZYkSZIqWJAlSZKkChZkSZIkqYIFWZIkSapgQZYkSZIqWJAlSZKkChZkSapySZJ8NkmSVUmSDFhr+7T1bP9NkiTnJEmSVnwsS5JkRpIkX1lr/LeTJJm91rZ8kiRfSJLk8SRJFidJsjBJkqeTJJmQJMmgjvtuJSl7FmRJqn73AHngkNUbkiQZDOwGvNnE9tHAveVNRWDL8sco4IfApUmSfLq5kyVJ0gO4HbgYiMDhwB7ABGB/4Ox2+r4kqSpZkCWpyqVp+gowBziiYvPhwPPApCa2J/y9IJOm6Vvlj5fSNP0FMA3YZz2n/CLwMeCoNE1/kKbpE2mavpym6R1pmh4P/KZdvjFJqlIWZEnqGu5l3SJ8X/lj7e0z0jR9Y+0DJCWHAiOBv6znXJ8G7kvT9NGmdqZpuqBt0SWpa7EgS1LXcC8wOkmSoeXHhwP3Aw8Co9bafk/F83JJkixJkmQJsLL8nJ+laXrDes41AnihXdNLUhdiQZakruG+8ucjkiTZGtgReDBN0/coTZlYvX04axbkIjCm/LEHcCZwVpIk317PuRIgbd/4ktR15LMOIElqWZqm85MkeZbSdIqewFNpmi4s776/YnsReGCt51auUPFCkiTbA99JkuS/0jRd0cTpXqR0A6AkdUteQZakrmP1POTV849Xu79i+5MVxbk5BUoXSHo1s//3wOFJkhzQ1E6XeZNU6yzIktR13AtsD5zMmgV5CrBNefu9az8pSZJh5Y9tkyQ5Hvg34N40TRc1c54fl49zZ5IkX02SZJ8kSbZLkuToJEluAc5qv29JkqqPUywkqeuYAqwCegMPr96YpumiJEn+CoxlzfnHADlKayVD6crx68BE4FvNnSRN01VJkhwDfJ7SihbfoTR1Yw7wR1zmTVKNS9LU+zAkSZKk1ZxiIUmSJFWo6ikWIYQxwJWUfp1YAP4lxvh4COFQSu8e9VJ56M0xxu9mk1KSJEm1pKoLMvB94Dsxxj+HEI4tPz60vO+hGONxmSWTJElSTar2KRYpsEn564FAQ4ZZJEmS1A1U+xXkfwPuDCH8gFKZP7Bi3wEhhGcpleavxhj/1tQBQgjjgfEAMca9OzivJEmSqlvS0oDMC3II4R5gWBO7JlBa+P7LMcabQggB+BVwJPAUsF2McUl56sUtlN5edR0xxquAq8oP04YGL0JLkiR1R/X19a0aV9XLvIUQFgKbxhjTEEICLIwxbtLEuJeBfWKM81s4pAVZkiSpmyoX5BavIFf7HOQG4JDy14cDswBCCMPKhZkQwlhK38e7mSSUJElSTcl8ikUL/h/w4xBCHlhOeS4xcCrwzyGEArAMOD3GWL2XwiVJktRlVPUUiw7gFAtJkqRuqlamWEiSJEmdyoIsSZIkVbAgS5IkSRUsyJIkSVIFC7IkSZJUwYIsSZIkVbAgS5IkSRUsyJIkSVIFC7IkSZJUwYIsSZIkVbAgS5IkSRUsyJIkSVIFC7IkSZJUwYIsSZIkVbAgS5IkSRUsyJIkSVIFC7IkSZJUwYIsSZIkVbAgS5IkSRUsyJIkSVIFC7IkSZJUwYIsSZIkVchnHWB9Qgh7AFcC/YGXgU/FGBeV910IfBYoAl+MMd6ZVU5JkiTVjmq/gvxL4IIY40eAicDXAEIIo4DTgd2Ao4GfhRBymaXUh/L5PLle/Ul6DCDXqz/5fFX/G0ySJGkd1d5edgGmlL++G7gT+AZwInB9jHEF8FIIYTYwFng0k5QCSuV48coeXHLNVOYtWMbQQX246Jx9GdATCoVC1vEkSZJapdoL8vPACcAk4DRgm/L2rYCpFeNeL29bRwhhPDAeIMbI4MGDOyxsd7dkWZFLrnqIeQuWATBvwTIuueYJLv/CwWy6qRf4JUlS15B5QQ4h3AMMa2LXBOBc4CchhG8CtwIry/uSJsanTR0/xngVcNXqMfPnz9+4wG107rhzOvV8Wbp84pUfluPV5i1YxhuvvsnXTj4vo1Sd5+op12QdQZIkrUd9fX2rxmVekGOMR7Yw5OMAIYQRwCfK217n71eTAbYGGto/ndpi8XsLGTqozxoleeigPixeuDjDVJIkSW1T1TfphRCGlj/XAV+ntKIFlK4mnx5C6BVC2AEYDjyeTUqtFv/313zhpF0ZOqgPUCrHXwxjuOnx+STbjsw4nSRJUutUdUEGzgghzARmULpC/GuAGOPfgAi8AEwGPh9jLGaWUgDMeGYG13z7B3z24CFccs5H+OzBQ/j1d/6HGc/MIL/n4eT2OxZ69c06piRJ0noladrk1N1alTY0dO5MjO40B3l96nbcnbpRB0BhFcVnHyB9c27Wkdqdc5AlSapu5TnITd3LtoZqv4KsGtE4dxqFByLpssXkxx5Dbq8jIN8z61iSJEnrsCCr8yxZQHHKTRRnPEGy1Qjyh51OMnjrrFNJkiStwYKszpU20vji4xQfugkaC+QPOpG60R+FOtdJliRJ1cGCrEyk78+j8ECkOHcauZ32IH9oINl0SNaxJEmSLMjKULFA43MPUfjLJMj3JHfwJ6kbsQ8kviwlSVJ2bCLKXPrO6xTuv570jdnkRu5H7uBToP+mWceSJEndlAVZ1WHVCopP3UPhickk/QaSPyRQt8NHsk4lSZK6IQuyqkraMIfCfdeRvttAbvdx5A44AXr3yzqWJEnqRizIqj4rllKcehvFZx4g2WwY+cPPINl6RNapJElSN2FBVtVqfOVvFB64gXTRe+T3/hi5fY6CHr2yjiVJkmqcBVnV7YOFFB+eSPGFR0m23KF0NXnotlmnkiRJNcyCrC4gpXHWUxQevBFWLid/wPHU7X4I5HpkHUySJNUgC7K6jkXzKTz4R4qznqZu+93IH/oPJIOGZZ1KkiTVGAuyupbGIo0v/IXiIxMhScgdfDJ1I/f3zUUkSVK7sVWoS0rffZPCA9eTvjqD3Ii9yR9yGgzYLOtYkiSpBliQ1XUVVlF85n4Kj90OvfqW3lxkpzFAknUySZLUhVmQ1eWlb71M4f7rSOe9Qm70QeQOOhH6Dsg6liRJ6qIsyKoNK5dTfPzPFJ66h2TgYPKHnU6y7cisU0mSpC7Igqyakr72IoX7byBdMI/8noeTG3ss9OqTdSxJktSFWJBVe5YtpviXSRSfe4hk6DbkDzuDZMsdsk4lSZK6CAuyalbj3GkUHoywbDH5sceS2/NwyPfMOpYkSapy+awDrE8IYQ/gSqA/8DLwqRjjohDC9sB04MXy0KkxxvMyCanqtngBhSk3UbfLPtQN35v84K0pPn0v6fw3sk4mSZKqVFUXZOCXwFdjjA+GEM4FvgZ8o7xvToxxTHbR1GWkjTTOeJz07VfI7XUk+YNOojj7GRqnT4XGYtbpJElSlan2KRa7AFPKX98NfDLDLOri0gVvU3jgBopzp5HbeQz5QwIMHJx1LEmSVGWqvSA/D5xQ/vo0YJuKfTuEEJ4OITwYQji486OpSyoWaHzuIQp/uRV69CQ/7lTqRuwNiW8uIkmSSjKfYhFCuAcY1sSuCcC5wE9CCN8EbgVWlve9CWwbY3w3hLA3cEsIYbcY46Imjj8eGA8QY2TwYK8YCtJ3XqNw//Xkdh9HbuT+JFtsT/Gpe+CDhRt8TF9bkiTVhiRN06wztEoIYQTw+xjj2Cb2PUBprvKTLRwmbWho6Ih4zTp33Dmdej61XVK/M7k9DoG6HI1/+wuNLz+/Qce5eso17RtMkiS1q/r6eoAWf21c1VMsQghDy5/rgK9TWtGCEMKQEEKu/PWOwHBgblY51bWlDbNLb1X9bgO5PQ4hd8Dx0Ltf1rEkSVJGqrogA2eEEGYCM4AG4Nfl7eOAaSGEZ4EbgfNijO9llFG1YPlSilNvo/jsAySbbVl6q+qtds46lSRJykCXmWLRTpxioZb1G0huryOp22wYja/PojjtQVi1osWnOcVCkqTqVhNTLKRMfLCQ4sM3U5w+laR+R/KHn0EydNusU0mSpE5iQZaakqY0zvwrhSk3wsoV5A84nrrdx0Eu84VfJElSB7MgS+uzcD6FByPF2U9Tt/1o8of+A8mgLbJOJUmSOpAFWWpJY5HGv/2F4iO3QF0duYNPoW7kfpD4fx9JkmqRP+GlVkrfbaBw//Wkr75IbsQ+5MedCgM2yzqWJElqZxZkqS0Kqyg+cx+Fx+6APv3IH3IadTuNyTqVJElqR95xJG2A9K2XKNz3Frkxh5IbfRDJsO2Zv3AFgwf2yjqaJEnaSF5BljbUymUUH/8zhafuJRk4hO/99nkeef4dutna4pIk1RwLsrSR0tdmULj/erYd2o/f3vkyP580m0VLV2UdS5IkbSALstQeli3my2EXTj1kG/728kK+e83zPDN7QdapJEnSBrAgS+2kLkn42D7DuOgfR7Fp/578fNJsfjP5JZatKGYdTZIktYEFWWpnWw3uywWfGskx+23Joy/M53u/fZ6Zry3KOpYkSWqlFgtyCGGzEMLRzew7OoQwqP1jSV1bPlfHSR/dmq+dPpJcXcL/xBe58YFXWVVozDqaJElqQWuuIH8d2LuZfXsCE9ovjlRbdqrvz9fP2o1xewzh7r++zSW/f4FX3/4g61iSJGk9WlOQjwN+0cy+q4AT2y+OVHt69chx5pHb84VTRvDB8gKXXjudOx5roNjocnCSJFWj1hTkYTHG+c3sew/Yoh3zSDVr9A4D+ebZoxmz86ZMevgNfnD9dOYtWJ51LEmStJbWFOQFIYRdmtk3Ani/HfNINa1/nzzjj9+Zzx67I2+9t5zv/fZvPPjsPN9cRJKkKtKagjwR+EkIoU/lxvLjHwI3dkQwqZaNHbk53zx7NDtv1Z9r73mFn9w8k/eXrMw6liRJAvKtGPMN4D5gbghhMvAmsCVwFPAa8K2OiyfVrkEDevLFT47gwWff4cYHX+M7v3meM4/Yjn133TzraJKkTpDP50lzvWlsTKirS0mKyykUClnHEq24ghxjXAwcSKko9wb2KX/+BnBweb+kDZAkCYeOGco3Pr0bWwzqzS9vn8svb5/DB8v8C1KSalk+n2fxyh5c9POpjL/0Pi76+VQWr+xBPt+aa5fqaEk3m/uYNjQ0dOoJzx13TqeeT9m5eso1G/X8YmPK5Mff5LZHG9ikb56zj9qBUdsPbJ9wkqSq0diYkuvVn6//4jHmLVj24fahg/pwyT/vT3HFkgzT1bb6+nqApKVxLf4zJYRwbktjYoxXty6WpObk6hI+sX89o3cYyK/vmMuPb5rJoWOGcsq4renVI5d1PEnqHPd2nUqRprCsWMcHxRwfFPN8UMjxQTHHkkK+tK38+IO1Hi8r5rjkXz66RjkGmLdgGY1Ll8ADXefPYKMc0WLFzExrruN/uoX9KbBR/yVDCKcB3wZGAmNjjE9W7LsQ+CxQBL4YY7yzvP1o4MdADvhljPHSjckgVYvttujHRf+4G7c88jr3/vVtXnhlIecesyM7bNk/62iSVJPSFFY01q1RYpcU8iwt5lhSWLP8VpbepcUcjeu5GNmnrki/fIF+uSL98kWG9Fz54eOehQ8YOqjPOleQ61YtpdgZ37TWq8WCHGM8rBNyPA+cwlpvSBJCGAWcDuwG1AP3hBBGlHf/FPgY8DrwRAjh1hjjC52QVepwPXvUEQ7dlj123JRrJr/EZddN55j9tuS4/evJ5Vqz+IwkdU8rG5Nyia28optjaTFf/lxReisKbzFt/u/WXnXFD0tuv1yRQT2Xr/G4f75A31yR/vki/XIF+uWL9M0Vya3nF/n5Ofdx0ZlHc8m1zzFvwTKGDurDRWd+hGTG5A74U1FbbdBM8BBCT2A0MDfGuNHrIMcYp5ePu/auE4HrY4wrgJdCCLOBseV9s2OMc8vPu7481oKsmrLLtpvwzbN344b7X+WOqW/y/NyFfObYHanfvE/LT5akLqzQmDQ5TWFJMcfSQo4lxXz5c26NQrxqPUW3R9JYLrWlEjus1wr69fv74365Iv1za1717Zsr0qOu/e/XKrz7BgNmTOaSMw+msUdf6lYtJZkxmcK7b7T7udR2rZmDvAml6Q+jgEcpvb30w8AOwNIQwknx/7d379FV1ne+x9872blfCCHkHkjCJYJ4qVWw3qZaBbxMvSz91Vbxdjqu02nX2FNdHVtrW1ttx7NmOfW0HUfHo+CF0a9axlYFkdZ71UrnOBRQUAiSEAkggeyE7ECSff54noQdmkiAZD/J5vNaa6/kuexnf7frYe+Pv/wuZitGqL4K4O247UZ/H3hTzMXvnzPQBZxzNwI3ApgZRUVFI1CmCCN6b91ydSnvrNnKvy1Zzc8eW8tV86Zz4WnVpKQcdJyBiEigurt7aOvYR1vHPiJ79tG2Zx+RPXv/ajuyxztn9/Y62rtT6ewZfOxFaqinXwtuUfpeJoc7/NbcbrJTu/zW3O5+gTh9BILukej6dAu8+QTAUdmtYjRnsqG0IN8HFALPApcAX8Hr+/sgcANwF3DQgOycWwGUDnDoNjN7dpCnDfTtH2Pg6ekGvOvN7AG8UA8Q27FjsFWzRY7MSN9bU0rC3L7gWB5dvomFz3/AW6u2cN38GgrzM0b0dUVEAHpiMTo6u2nv6KI92vvwttv67evyz+mmPdpFR+fg0S8UgpzMsP9IJTcrTElue1+ozfVbcHPjWnRzUrvJSOkhpPaBMS+ITObPYnFQQwnI5wG1Ztbmd2XYBvzazLqdc/8K3DmUFzKzc4dUUX+NQFXcdiXQO0/bYPtFklZ+Thp/f8lU/rh6B0++vJk7Fq3hynMmcerMCYT0bSEiQxCLxYju7TkgzHohd0+0m7ZoF3uiB4Tejm72dHbxWTPDZmekkpPlhd3crDRKCjPjwm+Y3Kww2ZlhcrNS+/ZlZqSScuBn1+/fGNn/ACJDMJSAnGlmbQBm1uKcazOzbn+7xzk3kiOGfgssds7dgzdIbxrwJ7yW5WnOuRpgC95Avq+NYB0io0YoFOL04yZSNymfh5duZOGyet77qIWrz6smLzst6PJEJEFisRh7u3r6tdb2hd1oF3s6umiLdveF3T3+/vZoNz09gyfdzPSU/cE2K8yE/Az/99S/Drt+y292ZlhdviSpDCUgh/wgGhps+0iLcM5dCvwSmAg875x7z8zmmdka55zhDb7rAr7ZG86dc98CXsSb5u0hM1tzpHWIjCVF4zK42R3Dij9v5dk3t3DHotUsmFvDCVMKgi5NRA7Rvq4DW3S7+4fdQVp4u7oHD7rp4RRyMve36pYVZfXrztDX2psZJjtrf9jVTDkiQ1hJzznXc5BrxMxsrKxioJX0ZMQc6Up6R2LL9j08tHQjjds7OH1WEe7sSWSmj5V/liKjTzgcJpaaSU9PiJSUGKHuKF1dB18Cvru7h/b4frrxrbsHtPDGB+G9XYN/1YZTQ/1adHMyUw/Y3t+6m+vvy84Ik542RoPuGFooRI5QAAuFDNtKesBXgVfNbOuRFiUiI6NiYja3fm0mz73VxIvvfsIHm1u57vxaplfmBV2ayJgTDoeJ7E3jZwvf7puf9rsLTqapeTeN2yK0d3Qf0GVhf9iN7h18QFpK74A0P9QW5mVQVdw/7OZm9vbT3d/Kmx5O0RgDkQQbSkD+KTDFObcBeA14BXjNzDaPZGEicmjSwilcemYlx9WOY+Gyeu558gPOPbmUi0+vIC08RluSRBKoJxZjy/YOxo8fz92Pvtu3wtm2lg7+96Mr+frFx2EvryUEZPe24maFyc9Oo6wwq1+3hf6D0rywm5meqqArMkYMZSW96c65EuAs/3ELsNA5twUvML9qZg+ObJkiMlRTK/L4wYJjeebVBl5auZU1m7ylqquKs4MuTca4nTfdFHQJw6oHaE4rYGNmCRszSqnPKKYjNYOf//3p/Zb/BS8kV6Tv4/ZGI7NnHykDzyw6qA7/MVYU3ntv0CWIBGpIK+mZWTPwlP/AOVeAt/jGd/Bmj1BAFhlFMtNTueq8ak6YWsCiFzfx88fX8renlTPvlDKNNJejVl8gzihhY2YJ9RkldKR684gX7otwbEcDtdFm8ndNoXh8Vr+QXDw+i9TILrJ79gZUvYgk0pACsnMuBJzI/lbk0/DmHTbg9RGrTkSOyKyaAn507SwWr9jEf76xhUJ6mh8AABS/SURBVL9s3M1182soHp8ZdGkiI64H2NYvEBezJ9W798d3RZjZ0UBtZzO10WYKuvf0PW/vM4u59Zpv8E9Pv9/XB/nWy2fQ8ch9Ab0TEUm0oSw1/RxwErAOb4npB4DrzCwywrWJyDDIzQrzdxdN4cQPdvIfv/+YOx9dw+V/U8WZx09Uf0hJKjGgOW0cGzNKqPdbiNv7AnEbMzq2UOMH4vHd7YNeJ/rhh2Q+ch8/uexKyMuHSCsdj9xH9MMPE/RORCRoQ2lBrgM6gXpgA/CRwrHI2BIKhZg9YwLTKvNYtKyex1d8zH9v2MU1c6sZl5sedHkihyUGbAuP8/oQHxCIC7raqOvY0tdC/FmBeCDRDz8kevdPR6BqERkLhjJIb9oBg/S+7ZwrAt7E617xhpm9N7JlishwGJ+Xzj9cPp1X39vGM681csei1Vx1bjWfrysMujSRg4oB28P5fYPqNmbuD8TjutqZfkAg1t9HRORwHekgvR/grX6nFQlExoiUUIizP1fCjMn5PLy0ngee28DsDS1cec5kcjKH9JEgkhC9gbg+s4QNfreJttQswA/E0SZqo72BuE2BWESGzeEO0jsDKABWAlryRmQMKi3M4rtfncHSd5p4/q0m1jdEuG5+DTMmjwu6NDlKxYAdfS3EXiCO+IE4v6udqdGt1Ea3UtvZTGGXArGIjJyhDNJ7Hm/WinTgHeBV4FfAW2YWHdnyRGQkpaaEuOgLFcyqKeDhpRv5xdPrOftzxVx2ZiXpafrDkIwsLxDnUe/PMrExo4RI2JuvO69rD7XRrUyJNlPTuZUJCsQikkBDaUF+HbgLeNfM9o1wPSISgOrSHG67+liWvNHIH/6rmbWbWrnhghqqS3ODLk2SSAz4NJzXbx7i1vhA3NlM7e5majubmdAVUSAWkcAMZZDePyWiEBEJVnpaCl85exLH1xawaNlG7l78PhecWs4Fc8pITdVS1XLoYsDOcG7fgLqNGcW0hnMAyOvuoCbaTG1rM7XRrRQpEIvIKKIROSLSz4zJ+fzw2lk8+fJmnnurib9s3MUN59dSOiEr6NJklIsBLam5bMwsYYPfQrzbD8S53R3egLrW1dREm5nY1apALCKjlgKyiPyV7Mww159fywlTCnjspY+587E1XHpmJWd/roQULS4iPi8Q5/jzEJey8YBAXOMH4loFYhEZYxSQRWRQJ00vZEpFHo8ur8debmDVht1cO6+awvyMoEuTgPQFYr8f8a6w1089pztKbbSZmtY11EabKe7arUAsImOWArKIfKZxOWl885JpvPmXHdgrm/nJI2u48pxJzJkxQUtVHwU+be1kfUOEdQ0RPii/hJa4QFzT2cxZrWup7WymeJ8CsYgkDwVkETmoUCjEGcdPpG5SHguX1fPw0nre+2gXV587mdzstKDLk2G0s7WTdQ0R1jdEWN8YYcfuTgByMsNU793JGXGBWEM3RSRZKSCLyJBNLMjkZncML63cym//uIU7FkVYMLeG46cUBF2aHKaWyF7WNbT2tRLvD8SpTKvM40snlVBXlUdZURa7vr0w2GJFRBJEAVlEDklKSoh5s8s4tmYcDy/dyK//80POOK6IK744icx0LS4y2rVE9rK+McK6za2sb4ywfZcXiLMzU5lemcc5JxVTV5VPeVGWBmSKyFFLAVlEDkvlxGxu/dpMfvfHLSx/dysfbI5w/fwaplbmBV2axNnVtrevdXh9QyvbegNxhtdC/MUTvUBcMVGBWESklwKyiBy2tHAKl51VxfG1BTy8rJ5/fvIDzjullC+fVkFaWD1Ug7C7ba8XhhsjrGtoZVuLF4izMrwW4r85sZjplXlUTswmJUWBWERkIIEHZOfcFcCPgRnAbDNbGXfse8D/ALqBfzCzF/39m4CIv7/LzE5OcNkiEmdqZR63X3MsT73SwPJ3t7J2026uP7+WyonZQZeW9Ha372N9XB/i5pYoAJnpqUyrzOWs44upq1IgFhE5FIEHZGA1cBlwf/xO59xM4ErgWKAcWOGcm25m3f4pZ5vZjoRWKiKDykxPZcHcak6cWsAjL9bzs8fW8uXTK5h7cqmC2TBqbd/nDapr9Gaa2LqzfyA+8/iJTK/Ko0qBWETksAUekM3sfQDn3IGHLgaeMLNOoN459xEwG3grsRWKyKE4rraAH107i8dXfMyS1xtZtWEX159fw8SCzKBLG5Na9+zzplzzW4k/6QvEKUytyOP0WUVMr8qnqjibVAViEZFhEXhA/gwVwNtx243+PvBWOF3unIsB95vZA4NdxDl3I3AjgJlRVFQ0QuXK0U731n5FwPevL+W195p48LdrufPRtVx34QzOPaVSi4scxO62TtbW72T1Ru/RuK0N8FqIZ1QXcu7sycyqLaSmPJ/U1MT2896Z0FeTIAX5eaY/DR89RvP3ZkICsnNuBVA6wKHbzOzZQZ420LdozP95upk1OeeKgZeccx+Y2WsDXcQPz70BOrZjh/7pycjQvfXXjq3K4PYFM1n4Yj3/tmQ1b/53Awvm1jAuR4uL9Irs2ceHjZG+xTmaPu0AICPNayE+5cxK6qrymFSSE9dC3EVLi+KqjBx9nkkiBHGflZeXD+m8hARkMzv3MJ7WCFTFbVcCTf71en9uc84twet6MWBAFpFgFeZn8O3L63jl/23jN683cMfC1Vx93mROml4YdGmBaNuzj/VbIn2D6pp27A/EUypymTNjAtOr8phckp3wFmIREfGM5i4WvwUWO+fuwRukNw34k3MuB0gxs4j/+1zgJwHWKSIHkRIKcc5JJcyYnM9DSzdy/+82cOrMXXzl7ElkZ47mj6Ej197R5Q+o8wbWNW73AnF6OIWpFbnMPqaQ6VX5VCsQi4iMGoF/MznnLgV+CUwEnnfOvWdm88xsjXPOgLVAF/BNM+t2zpUAS/xBfWFgsZktC6p+ERm6sglZ3PrVGbzwzie88HYT6xoiXDuvhhmT84Mubdi0d3R5XSb8ULxlewcxvDmjp1bkcvEZhdRV5jG5NIewArGIyKgUisViBz8recSampoS+oI3nHVdQl9PgvPQawuDLmFMqf+kjYeX1tPcEuWck0q49IxK0tPGXmBsj3qBeL3fh7hx+56+QDylPJe6qjymV+VRnQSBeOdNNwVdgiRI4b33Bvfiv38ouNeWxPrSDQl/Sb8P8kFHiwfegiwiR6easlx+sGAmv3m9kT/8V3Pf4iLVpTlBl/aZ9kS7+HBLmzcXcUOExm29gThEbXkuf3taBXVVXguxVhMUERmbFJBFJDDpaalcec5kTphSwKJl9dy9eC0XnFrOBXPKRk1/3I7OLj5sbGN9QyvrGiI0+IE4nBpiSnkuF51WTl1VPtUKxCIiSUMBWUQCN2PyOH547Sye+MNmnnuridX1u7l+fg2lE7ISXktHZxcfbWljXUOEdQ2tXiCOeYG4tjyXC79QTt2kfGoUiEVEkpYCsoiMCtmZYW64oJYTphbw+EubuPOxNVx2ZhVf/FwxKSO4uEhHZzcf+dOurW+M8HFz+/5AXJbLhaeWM70qj9qyXAViEZGjhAKyiIwqn59eyNTyXB5ZvoknX97Mqg27uGZeNYX5GcNy/eheLxD3LsyxubmdHj8Q15TlcsGccuqq8qgpyx2TgwZFROTIKSCLyKgzLjedb106jdf/sp2nX2ngJ4+s4atfmszsYwoPealqLxC39c1D/PFWLxCnpoSoKcvh/Dm9LcQ5pKeljtA7EhGRsUQBWURGpVAoxFnHFzNjUj4PL63noRc28t5HLVx/wTSysnPo6QmRkhIj1B2lq6ur73nRvd1sbPJmmVjX0D8QV5fmMH9OGdMr85lSrkAsIiIDU0AWkVFtYkEmt3zlGJa/u5X3G9ppjYb56aK32dbSQfH4LL537Sns3NnOynXbWbe5lU3Ne+jpiZGSEqKmNId5s8uoq8qjtjyXDAViEREZAgVkERn1UlJCzJ9TxnlfyObHD/6JbS3ecs3bWjr4+aJ3+frFx/Hiu1upLslm7sml1FXlMaVCgVhERA6PArJIEnhq7j8GXUJCXPjEnX3huNe2lg4mpHZz1vLfEO7upgtY4z+S0RXL7w66BBGRpKch2iIyZuxtiVA8vv/cyMXjs+j6dBfh7u6AqhIRkWSjgCwiY8aq+57m5ktn9oXk4vFZ3HzpTFbd93TAlYmISDJRFwsRGTOaV21g1V0PcvM3Lid9fB57WyKsuutBmldtCLo0ERFJIgrIIjKmNK/awEvfUD9cEREZOepiISIiIiISRwFZRERERCSOArKIiIiISBwFZBERERGROArIIiIiIiJxFJBFREREROIoIIuIiIiIxFFAFhERERGJE/hCIc65K4AfAzOA2Wa20t8/AXgaOAVYaGbfinvO54GFQBbwAnCTmcUSW7mIiIiIJKPR0IK8GrgMeO2A/VHgduCWAZ5zH3AjMM1/zB/JAkVERETk6BF4QDaz981s3QD7283sDbyg3Mc5Vwbkm9lbfqvxI8AlialWRERERJJd4F0sDkMF0Bi33ejvG5Bz7ka81mbMjKKiopGtTo5aurckEYK8z3YG9sqSaEHeZzsCe2VJtNH8vZmQgOycWwGUDnDoNjN79hAvFxpg36D9j83sAeCB3vN27NA/PRkZurckEXSfSSLoPpNECOI+Ky8vH9J5CQnIZnbuMF6uEaiM264Emobx+iIiIiJyFAu8D/KhMrNPgIhz7lTnXAi4BjjUVmgRERERkQEFHpCdc5c65xqBLwDPO+dejDu2CbgHuM451+icm+kf+gbwIPARsAFYmtiqRURERCRZBT5Iz8yWAEsGOVY9yP6VwKwRLEtEREREjlKBtyCLiIiIiIwmCsgiIiIiInEUkEVERERE4iggi4iIiIjEUUAWEREREYmjgCwiIiIiEkcBWUREREQkjgKyiIiIiEgcBWQRERERkTgKyCIiIiIicRSQRURERETiKCCLiIiIiMRRQBYRERERiaOALCIiIiISRwFZRERERCSOArKIiIiISBwFZBERERGROArIIiIiIiJxFJBFREREROIoIIuIiIiIxFFAFhERERGJo4AsIiIiIhInHHQBzrkrgB8DM4DZZrbS3z8BeBo4BVhoZt+Ke84rQBnQ4e+aa2bbEli2iIiIiCSpwAMysBq4DLj/gP1R4HZglv840FW9YVpEREREZLgEHpDN7H0A59yB+9uBN5xzU4OoS0RERESOToEH5CPwsHOuG3gGuNPMYgOd5Jy7EbgRwMwoLy9PYImw7KPlCX09OTrdtPrRoEuQo0D5U08FXYIcDRb8IOgKRBITkJ1zK4DSAQ7dZmbPHsYlrzKzLc65PLyAvAB4ZKATzewB4IHDeA05As65lWZ2ctB1SPLSPSaJoPtMEkH32eiTkIBsZucO8/W2+D8jzrnFwGwGCcgiIiIiIodizE3z5pwLO+eK/N/TgIvwBvqJiIiIiByxwPsgO+cuBX4JTASed869Z2bz/GObgHwg3Tl3CTAX+Bh40Q/HqcAK4N+DqF0+k7q1yEjTPSaJoPtMEkH32SgTisUGHNsmIiIiInJUGnNdLERERERERpICsoiIiIhInMD7IEvycM5V4c0mUgr0AA+Y2b3BViXJxjmXCbwGZOB9hj1tZj8KtipJVs65VGAlsMXMLgq6Hkk+/nirCNANdGm6t9FBLcgynLqAm81sBnAq8E3n3MyAa5Lk0wmcY2YnACcC851zpwZckySvm4D3gy5Ckt7ZZnaiwvHooRZkGTZm9gnwif97xDn3PlABrA20MEkq/qqZbf5mmv/QaGMZds65SuBC4C7gOwGXIyIJpIAsI8I5Vw18Dngn4FIkCfl/9v4zMBX4tZnpPpOR8Avgu0Be0IVIUosBy51zMeB+fwVgCZi6WMiwc87l4i0B/m0zaw26Hkk+ZtZtZicClcBs59ysoGuS5OKcuwjYZmZ/DroWSXqnm9lJwPl4XRPPCrogUUCWYeYv4PIM8LiZ/SboeiS5mdku4BVgfsClSPI5HfiyP4DqCeAc59xjwZYkycjMmvyf24AlwOxgKxJQFwsZRs65EPB/gffN7J6g65Hk5JybCOwzs13OuSzgXODugMuSJGNm3wO+B+Cc+yJwi5ldHWhRknScczlAij9uJwdvxeCfBFyWoIAsw+t0YAHwF+fce/6+75vZCwHWJMmnDFjk90NOAczMngu4JhGRw1ECLHHOgZfJFpvZsmBLEtBS0yIiIiIi/agPsoiIiIhIHAVkEREREZE4CsgiIiIiInEUkEVERERE4iggi4iIiIjEUUAWEZGDcs694pz7etB1iIgkggKyiIiIiEgcBWQRERERkThaSU9EZJRyzm0CfgVcA0wGlgHXmll0kPO/CDwG/CvwHaANuM3MHvePZwB3AQ7IAJYA/8vMOpxz44FHgTl43w1vAv/TzBoHeJ0y4EXgETP75+F6vyIio4VakEVERjcHzAdqgOOB6w5yfilQBFQA1wIPOOfq/GN3A9OBE4Gp/jk/9I+lAA/jBfFJQAdeOO9fjHPVwKvArxSORSRZqQVZRGR0+z9m1gTgnPsdXrg9mNvNrBN41Tn3vPdUdyfwd8DxZrbTv97PgMXA98zsU+CZ3gs45+4CXj7gujOBH/jn/8cRvi8RkVFLAVlEZHTbGvf7HqD8IOe3mFl73PbH/nMmAtnAn51zvcdCQCqAcy4b+Be81urx/vE851yqmXX721cBHwFPH95bEREZG9TFQkQkuYx3zuXEbU8CmoAdeN0mjjWzAv8xzsxy/fNuBuqAOWaWD5zl7w/FXevH/nUWO+dSR/JNiIgESQFZRCT53OGcS3fOnQlcBDxlZj3AvwP/4pwrBnDOVTjn5vnPycML0Lucc4XAjwa47j7gCiAHeNQ5p+8QEUlK+nATEUkuW4EWvFbjx/FmovjAP/aPeF0k3nbOtQIr8FqNAX4BZOG1EL+NN2PGXzGzvcBlQDHwkEKyiCSjUCwWC7oGEREZBr3TvJlZZdC1iIiMZfo/fxERERGROJrFQkRkDHHOfR/4/gCHXseb51hERI6QuliIiIiIiMRRFwsRERERkTgKyCIiIiIicRSQRURERETiKCCLiIiIiMRRQBYRERERifP/Ae+wFWdJUAv9AAAAAElFTkSuQmCC\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,-2.5])\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([-115, -85])\n",
"\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*[図]上段：候補となるモデル（ピーク本数）に対するWAIC値, 下段：〃に対するWBIC値*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 6. まとめ\n",
"- iterを多めに取った以前の実験とは異なり、WAIC・WBIC共にピーク本数＝３で最小値を取った。\n",
"- WAIC最小となるモデルはn_peak=5からn_peak=3に変わったが、その差はかなり小さく、MCMCのサンプリングの影響の範囲と考えられる。\n",
"- WBICは変わらず真の本数を選択できており、この問題に対しさほどのサンプリングは要さないと考えられる。\n",
" - 真のモデルを選択する類の問題には、一致性を持つWBICを用いるのが適当である。"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment