Model Fits and Parameter plots.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: All of the model fits presented here have been pre-run using the script `fit_choice_models.py`. This can take a lot of time to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "sns.set_context('paper')\n",
    "\n",
    "\n",
    "fig_size = (5.0, 2.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Experiment 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chance: 132252.482051\n"
     ]
    },
    {
     "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>LOO</th>\n",
       "      <th>dLOO</th>\n",
       "      <th>LOO_se</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Sticky-Choice</th>\n",
       "      <td>74117.013931</td>\n",
       "      <td>22693.468864</td>\n",
       "      <td>457.742150</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lin-Reg</th>\n",
       "      <td>71422.248629</td>\n",
       "      <td>19998.703562</td>\n",
       "      <td>452.491382</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bayesian-GP</th>\n",
       "      <td>58137.534089</td>\n",
       "      <td>6713.989021</td>\n",
       "      <td>438.439840</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF</th>\n",
       "      <td>55814.695037</td>\n",
       "      <td>4391.149970</td>\n",
       "      <td>434.124566</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kalman</th>\n",
       "      <td>54757.107327</td>\n",
       "      <td>3333.562260</td>\n",
       "      <td>429.144697</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Clustering</th>\n",
       "      <td>54121.496627</td>\n",
       "      <td>2697.951559</td>\n",
       "      <td>424.809565</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Linear/Clustering</th>\n",
       "      <td>53652.000256</td>\n",
       "      <td>2228.455189</td>\n",
       "      <td>422.527828</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF/Kalman</th>\n",
       "      <td>53632.482207</td>\n",
       "      <td>2208.937139</td>\n",
       "      <td>425.894747</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF/Clustering</th>\n",
       "      <td>51423.545067</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>424.499452</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            LOO          dLOO      LOO_se\n",
       "Sticky-Choice      74117.013931  22693.468864  457.742150\n",
       "Lin-Reg            71422.248629  19998.703562  452.491382\n",
       "Bayesian-GP        58137.534089   6713.989021  438.439840\n",
       "GP-RBF             55814.695037   4391.149970  434.124566\n",
       "Kalman             54757.107327   3333.562260  429.144697\n",
       "Clustering         54121.496627   2697.951559  424.809565\n",
       "Linear/Clustering  53652.000256   2228.455189  422.527828\n",
       "GP-RBF/Kalman      53632.482207   2208.937139  425.894747\n",
       "GP-RBF/Clustering  51423.545067      0.000000  424.499452"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loo = pd.read_pickle('Data/model_fits/model_fits_exp_lin.pkl')\n",
    "summary = pd.read_pickle('Data/model_fits/model_params_exp_lin_rbf_cls.pkl')\n",
    "loo.rename(index={'Kalman w/SC':'Simple Kalman'}, inplace=True)\n",
    "\n",
    "#note the \"Simple Kalman\" is modeled without sticky choice and wasn't included in the paper\n",
    "loo.drop('Simple Kalman', inplace=True)\n",
    "loo.drop('Scrambled', inplace=True)\n",
    "\n",
    "# calculate chance\n",
    "n = len(pd.read_csv('Data/exp_linear/rbfpred.csv')) # n_subj * n_trails_per_subj\n",
    "n_arms = 8\n",
    "chance_loo = -2 * n * np.log(1. / n_arms)\n",
    "\n",
    "print \"Chance: {}\".format(chance_loo)\n",
    "\n",
    "loo.sort_values('LOO', ascending=False, inplace=True)\n",
    "loo['dLOO'] = loo.LOO - np.min(loo.LOO)\n",
    "loo['LOO dLOO LOO_se'.split()]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAADKCAYAAACrHYtRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXncnNPd/9+fRLSEIIS6qapdbY2krV3sKqWlSlv0SVuipUVLPfXQoo2l1GMp2gblsdRSSy1JVSOJ2IklsVP7TyxBJRQl8vn9cc7knoy5t2Sumbkn3/frNa9zXec65/p+Z+65v3Ous3yObBMEQRC0Ln0a7UAQBEFQLBHogyAIWpwI9EEQBC1OBPogCIIWJwJ9EARBixOBPgiCoMWJQB8EQdDiRKAPgiBocSLQB0EQtDgR6OvAjjvuaCBevffV0sT3s1e/ukUE+jrw+uuvN9qFIOiQ+H62PhHogyAIWpwI9EEQBC1OBPogCIIWJwJ9EARNz/Dhwxk+fHij3ei1LNRoB4KgKfjXv+Dmm2GbbWCppRrtTVDB2LFjG+1CryZa9EEA8Mwz8I1vpDQIWoxo0QcBwAYbwBtvwIABjfYkCGpOBPogAFhoIRg4sNFeBEEhRNdNEAA89xx897spDYIWIwJ9EAC8/z489lhKmwRJp0q6VdLpFfl/lHS7pNskrZ/z2iSNl3SHpG0b43HQrETXTR2YNm0axx57bKPd6DUcffTR9Te61lpw1131t9sBkjYE+tveXNLvJX3B9r358om2n5W0OnAi8HXg58BRwFTgBmBcQxwPmpJo0QdBc7Ix7cF6HLBR6YLtZ/Phh8BH+Xh94E7b7wBvS1q8s5tLGilpsqTJ06dPr63nQdNRtxa9pAHApUB/4BPAj4GVbF+dr//O9o+r1BsGbGv7qB7YWgY4E/gU6T0eB7zXnftI+jwwxPZ53bUXtABTpsD228NNN6UZOI1nSeDpfDwDWKdKmROAM/JxX9suK78U8HZHN7c9GhgNMHTo0IYpdJ74QM8E1bpb/ueDl5kXd1qWerbovwNcbXsYsDnwH2C30sVqQX4++B1wVra1NfBWdyvafjCC/ALIssvCj36U0ubgLaA013MAFd9hSYcAj9q+LWd9VHb5Y+WDBZt6Bvp3gY0lLWN7FvAtYDtJEyUNknQbgKRN80DTBEl7lipLGiDpOkmfy9f65vyrJC1bVq4vsLztWwFsf2D7znx5PUnX5/svJmkhSZdKmpTThSQNkzQq3+uHku7K9tbMfl6Xz8+uy6e2ADJt2jSmTZtWX6PLLw+/+EVKm4M7gW3y8bbAnAEESdsDmwCjyspPlbSxpP7AANsz6+Zp0PTUM9BfBLwATJA0DrgO+IftYbbLOwlPBL5qeyvgLzlv8Vz/CNuPAjcDW+XuoH62XyurPwjosNPR9s7AWNI/0a6kVtEWwCOkQS0A8o/HN4BNsy9PkQa8Tsjnb0vauCM75X2g7777bpcfTtDO6NGjGT16dH2N/vvfcN99KW0CbN8PvC/pVmA28IKkI/Pl3wGfJf0v/THnnUTqohwHHF9vf4Pmpm599LY/BH4F/ErSt4BDgFkdlH09p7MlQQq459h+JBf5M3AEsAJwjaRBtP8obEMK9tV4OKcvkfpAlwfuz3mTgSHAq/n8s8D9tj8q82Vt4ERJBhYD7unk/c7pA21ra2v5XYpqyciRI+tv9PHHYehQmDwZhgypv/0q2D64Iuu4nL9mlbL/j9RNGQQfo56DsZ8BpuWA/xppoGixKkUtaWnbb0gqPXGcD6wo6Wu2/2r7GUltwB7At23PAIaV2XpZ0ua2b5XUjxTAYe6ttwQ8k6+NAYYC/yy7/gwwWFKfHOT7AE8AF9u+L9uJ6akF0NbWVn+ja62Vgvxaa9XfdhAUTD27bj4P3CZpIqkLZBQwUNKVksrXnh8BXC9pAqklDylA7w/sLWnLnDcW+CAH+Up+DByYbU0gtd6rcQ2wjqRJwHrAVaULuTvpKuCO7MvqpEfio/LClHHAp3vyAQRNTP/+qSXfv3+jPQmCmlPPrptrgWsrsncoO94sl7udNNBUzsSc7l5+S1IXTjVbrwPfrHJpYr5+QVnenlXKlMqdDVQOuu5azWbQy3n5ZTj3XNh332YakA2CmtArux4kHUAKuDs12pfu0NbW1pjVnkH3ee01OPNM2GWXCPRNyJqbharD/NArA30HLe0gmHc22ABefbXrckFDGHHGpY12oVfTKwN9byO0buaPeBoKgvkjtG6CANL0yo02SmkQtBgR6IMA4JOfhLXXTmlBSNpO0n9J6pvXZARBXWipQF8uX5DPf9fNehMl3SLpTkkNWK0TNJyVV4bzz09pAUj6PWmW2YF5Ed4ZXVQJgprRUoG+kh4KpW0DbAqMKMaboKmZNQvefDOlxbCm7cOAksZC36IMBUElLR3oy4TSLpB0et6Rp7ORvU8An8x1lDd8GC9pjKSlJH0ii5rdKOkySSPq8DYWaGbOnFkfgbMpU2DppVNaDLOzfpLzAsHCflGCoJKWDvQVTLS9GR3Pvb8ZmEb7tM2vAC/Y3pqkbf8D4GvAHbZ3BP7VmbEQNasN9913X30EzlZZBf7yl5QWw0+By4B1gSuAQ4syFASVLEjTK0uCZu8BSLqZ9PhcWkG7DbAWSZ7hXGBt4JuSdiB9TneSVuNOzeUf7MxYiJrVhiFDhtRH5GyppWD33bsuN+8MzI0GADpTPg2CWrMgtejnCra2t8kSya+U5T0M9MkzIp4ALsxlNgP+B3iWpIkDaeu2oGAGDBhQH5GzN96Aiy9OaTH8suL8J0UZCoJKWrFFv5ekjbou1iHnAj8EDgbOkDQ+558G/BX4i6S/A++Q9uwMWoHnnoN99kkKlksvXbPbZknubwPrS7ouZ/ermYEg6AYtFehtTyTpyFfmjyg7Hlbl+rCy4/FAKbhX28N2N9uz8nS5Z+bX56BJGDw4bTpS+3n0Y4E7SA2H03PeB8ArHdYIghrTUoG+ToyRtBjwz7ItCjslRM16AX36wKKL1vy2WUZ7BvBTScsBbaS9EMo3vQmCQolA30Ns79B1qaDX8cwzcNRRMGpUITNvJP2KtAPUKsDzpI13tq+5oSCoQgT6OhCiZo2nyyeqDz+EadNSWgw72P5S3gxnKzrYSyEIimBBmnUTBB2z5powcWJKi2FmTj8ibXb/ua4qSDpV0q2STq/IP1LStAq5jwsk3Z3lPL5dU8+DXk8E+iCoD+dJWoQ0e+sWumjRS9oQ6G97c2BhSV8ou3wusFeVanvl6cDxtBDMRcMCfRYgez63QG6vtZpfdwXNOqi7iKSzJE3Isgln5fxjJE3J/h5fO2+DhvPAA7DEEiktANuX2X7P9vW2B9v+TRdVNgbG5eNxwJwpw7ZfpWJdSD6/UNL1kj7TlT/lK7enT5/eg3cS9EYa3aK/KE9t/BlJYqBm9FDQrJJfAnfa3iovlrq87NqhtjcFNpAUm4P3EqZNm9a5Zs7yy8ORRxa2jaCkr+VGzRRJUyVN7aLKkrR398wAluqi/KG2NwF+A5zSlT+2R9seanvooEGDuvQ/6N00OtCXGADMlHRElgu+W9JgSctLugxA0kKlxUuS9s19l7dK2lDSwPxPNEHSGblMSdDsu/naZEnb57yuRM42tX1x6cT2pCplpgIrdPSGQuumuRg9enTnmjmf+hQcfnhKi+FXwN62N7C9vu2uVla/Rfq/IKdvdVbY9ps5vQ0o7E0EvZNGB/p9JE0CzicJPZ1ue0tS/+Nhtl8G+ktanKRFM07SMsAuwBbAV0mt7w1JomVbkRamlHN5fmrYBjisLL8zkTPNOZBukvSwpBXL8voCX6STBVPlLaZFC5ifHfSMkSNHdq6Z8847cPvtKS2Gp4CXe1D+TtJ3FmBb4K7OCksakNM16eJHIVjwaPT0yotsH5UXkpwL3CBpL2A27X2QV5MC+tbAKNI85A2ACWX3uQXYQtKfgb8BF5Vd20HSwaTgvWxZflciZwDY3l7SBbR/VqeQlCsvs/3aPL7voM50qZfzxBOw2WZJAmHIkCJcuBB4StIj+dy2d+mosO37Jb0v6VZgCvCCpCNtHyfp+8ABwEBJS9k+ELhE0lKk/5sfFvEGgt5LowN9ibdJj6cHAIOBVYFz8rWrgEuAfrafkfQ2cK/t3QEk9QP62v5lPn+QuQP9EcCWJK3528vyPyZyVjqWdIekfWyX7lP+OR1qexxBa/G5z8FDD8GqqxZl4XjgO8CL3a1gu/Lp9Licfx5wXkXZnefXwaB1aXSg30fSZqTNPkaRWu6T8gsA2zMlvU+egWB7et4IZBJpTvJ44JY8C6Yf7TMVStyQ73cP3X+kPRb4X0n7knRJniK0SVqbRRaBddct0sITuf88COqO7OaXSs9dMofmPvtex9ChQz158uRGuxF0xksvwdlnwwEHwAofG2NXtSo9ISueLkbStzGA7YPm9761oJHfzxMfeL2Q+/588DKF3LcJ6dZ3s9Et+i6RNBp4rbcG+aCX8OabSY/+m9+sFuhrQay7CBpG0wd623XYXihY4FlvPXj++ZrfVtLqtp8ijUMFQUNo+kDfCoSoWXNRZ8noXYGT+PjeBga+V09HggWXCPRBAPDoo7DXXnDJJWkGTo2wfVI+vND2nCnBkjapmZEg6IJGL5iqOZK2lDQ+r4a9WdImkmaU6dasUVF+hKQn8vW/Slo45z+R73G3pIMqyk6U9LdGvL+gIPr3h402Smkx/KLi/JCiDAVBJS3Vos+rZo8FdsnTMhcHVgMesr2VpC8BhwP7VlQ92fa5ko4ibQZxAzDd9jBJIq1KPKO8bF3eUFA/PvMZ+P3va35bxZ6xQRPQUoGeJGdwse2ZALbfBh5IsRqAB4EVO6gLSUiqkoWJTcBbipK42VyrZT/4AF5/HZZZBhZeuJbmYs/YoOG0WtdNG1lPRNK3c1fNb8uubwE8UaXezyQ9TtLM+XvOG6S0G9CjwM0VZSdKOqEzR0LUrHmpKnD20ENpWuVDD9XUlu0Ztp8H/pLTD4EDga5EzYKgZrRai/5lUrDH9p8l3QEcA6wnaQJpZewBkn4GDAcuyPVOJmmRXAkMBF5l7q6by8s0vrvVdWN7NDAaoK2trflXpS1AVBU3W3VVuO66IiUQTgCGkUT4JgFnApsXZSwIymm1Fv1YYG9JS+Tz0g/ZQ1lbflfbL9s+Oe/Ec0Gpou0PgLOAn5Tf0Gnp8NtU79YJeiFtbW0fFzlbcknYeeeUFkOp/3DRvANUdAcGdaOlAr3t6aTB2Guzdv3ZpJZ6d/kHsI2kT5C7brKmzizbU2rvcdA0vP46nHtuSovhKUn3kr6b/aiBrEIQdJdW67rB9niS0Fk5lefl5S8oO54NlPbm/Ngu0eVlgxbj+edhv/1g8OA0IFtjbO8raSHbsyT1Ab5WcyNB0AEtF+ibkba2tnqvxgx6yoYbwqxZ0Ke2D7mSTi+TG/4ZcILt2ZJOB0bU1FidGT58OABjxoxpsCdBV0SgDwIACfr2LeLO5bNrtiMNygJ0uYF3szN27NhGuxB0k5bqow+Ceebpp2G33VJaWzqacRUzsYK6ES36OhCiZs3JXN1ps2fDu++mtLaskTesV8Xx6rU2FAQdEYE+CABWXx1uvLGIO+9VdnxlB8cdIulUYChwf/nWgpKOJC28+pPto3LeusAfSD8kP7Q9dT59D1qEluq6kTRM0qh8/HVJV+UZDuVlVpZ0cWM8DBY0bN/S0aurupI2BPrb3hxYWNIXyi6fy9w/IgC/Br4F7JGPgwBo0Ra9pE1JrZ3hecpkEHTO/ffDJpvAHXekGTjNwca074E8DtgIuBfA9quS1q4oP9D2iwBliwbniZ5s8dfdsgvQ9n5NR0u16DNrAr8FvgHMylLFk3Lrfq5pFZLulPQ7SVMkfVPS1ZKmStogXz9d0i2SbpW0UlmdMyU9KGnHur+7oCbMnDlzjrgZACuuCCedlNICkLTsPFRbEpiZj2cAS3VRvk8Hx9X8maPFNH369HlwLehNtGKg3x74u+03gFnAV2xvATwGbF1RdmlgFPBl4DSSnOwPgO/n60fY3pK02nb/sjq/Imnl7E8HhKhZc3PffffNLWy27LJw0EEpLYZL5qHOW8CAfDwgn3fG7A6OP4bt0baH2h46aNCgeXAt6E20YqD/PbBJbm33B86TdAuwO1nwrIzptl+1PQ142vb7wDTaW06HS7qV9GPQVlbnNdsv0Yn+Tfk/0qKLLlq7dxfUhCFDhswtbjZzJtx8c0qL4UFJO0saJGmgpIHdqHMnsE0+3pa0L0JnvClpRUltpCeAIABas49+FvBN4CbgGuBJ29+WdBwf1xdxB8eStDQwzPbmkrajfeBrrnK1dT2oFwMGDJhb2Oypp2DbbWHyZBgypAiTQ/OrhPn4E+Zc2L5f0vu5sTEFeEHSkbaPk/R94ABgoKSlbB8IHA1cRvpeHljEmwh6J60Y6LH9pqR9gOuBPpKGklo4T/XgNv8C3sniaDFNrdVZZ50U7Avqo7e9FYCkxWy/04N6B1dkHZfzzwPOqyg7FdhsPl0NWpCWCvS2JwIT8/FjpG0Eq7F3LjPnn6J0bPu50nXgK1VslNcZNt9OB83BJz8Jq3X0dZl/JO1EGutZNA/2j7b9vcIM1oE1N9u20S4E3aSlAn2zEqJmvYAXX4TTToNDDoFPf7oIC0eSWts3ZgXLlYswUk9GnHFpo10IukkrDsYGQc+ZMQOuvz6lxfCR7f/QPsZTiIJaEFQjWvRBALDuuvDkk0VauFrSX4CVJF0EXFWksSAoJwJ9HQhRs+alXl1qtk+T9HdgHeAx24/UxXAQEF03QZB45BFYe+2UFoCkG20/ZvtK24+E3lJQT3pFi17SNsAvSfOD3yD5/RPb/+xG3ZWBVfIWg92xdRpwqO2P5tnhoPcxYADsuGNKa4iktUit+BUk7Zaz+9ECG48EvYemD/SSliEF+a/YflvSGsDvenCLlUkLU7oM9JL62D5knhwNejef/jScemoRd14OWBdYDFgv530AHFqEsSCoRtMHepKmzEW23waw/aSklwEkHQPcZnucpAuAY4DBwBHAOyRxs72BTSVtbHsbSb8kBf7ZQGke8/mkJ4Wxkr5DWm5+FLAi6Yfiuby586rAn4FXSP+438/z7oNeSEnYrK2tDf7zH5g2Ddra4BOfqJmNLEd8S9aVnw2sYPuJmhkIgm7QG/rolwde7kH5rwN72N4a+BswmvRDsY2k9Uj/aMNIS8SPyHWWBfa0/aeKez1ie1vSTIklgcOAg4Hdcp0OCVGz5mcuYbOHH4ZVVklpMewOjAWukNRX0g1FGQqCSnpDi/5lPi5GVqKa7swo4ChJC5GXi5exNjBM0sSyewNM6aBPvvRfPw1YAvgsMNX2R5I6jQi2R5N+ZGhra4v9QZuQuYTNVl8dbroppcXwPWBzYHz+/oTSXVA3ekOLfiywt6TFASStRmrlQ9KvWV6SSANeAM/b3pcUZH8KfEj74pQngJtsD8ut+u/k/I4kXSt/SJ4F1su69utUrxL0FuYSNhswALbbruaDsWWUvmOxYCqoO03forc9XdKvgRtyQH+TNJgFcDVJoXIXkggZwDGSNiL1oR9KapWfIOly23tKeiW36A1cSlK57C6/JemKv5btfThfby5oHl57Da64AvbYoyhN+rOACcBqkm7K50FQF5o+0APYHkf7lmqVVGrK/neVMluU3es4Pt6ls3fZ9WH58JiyvBEAkhayvVFu0d8OvNq190Gv4MUX4Sc/gY03LiTQ2748B/hVgWdsv1lzI0HQAb0i0DcRq0g6l7yhie1Z3akUoma9gCFD4MPiHtDyYP4uwEBgC0nY/t/CDAZBGRHoe4DtJyl7OgiCHvA3YAzwYqMdCRY8ItDXgdC6aV7mPGk99RQcfDCcfnpRM2/etj2qiBsHQVdEoA8CAAn69UtpMVyYJxU8SJ55Y/vqoowFQTkR6IMA0u5S115bpIX/Ap4n7WkMKdhHoA/qQo/n0UsaJmlU2XlPdGe6a2OFLC6GpJ9KmiTpNkmn57zbeujvKj0oX/P3EwSkjUf2tX1sfv2qqwqSTpV0a+l7X5a/bv5/uF3S+jnvAkl3S5oo6dtFvYmgdzLfC6Zs/7gWjkASFcuH2wH/kPRlYE1gy7xX67y0gIYB3Qr0WdSsZu8naH6mTZvGtGnT4L77UrfNffcVZepFST+QNETShpI27Kxwvt7f9ubAwpK+UHb518C3gD3ycYm98mLAP9fc+6BXM9+BvtS6zi2K03NL4+ict5qkmyTdIumonPfd3OqYLGn7srpnAjfm2w4DbiF9kU+2XerTvKXC9sQsdUBJ1kDS8bmlM0HSSsAI4BRJp0haRNKlksZLulxSP0kj8vEYYP2y9zNR0gmS7pX0/Zy3i6T7JI3u6qkitG56B6NHj056NyutBH/8Y0qLYSHgS8CPgB/ntDM2pn3tyDhgo7JrA22/aPslkjQHpK6gCyVdLykkkIO5qHUf/UTbB0u6m7Tj/XEkhccXc4BdEbjc9vmSlgD+QvvK1Ntt/yivfl3U9juSeipoBrApsLnt2fleF9CucPlj4Drbl0r6IUloCuAt23sCaO7BuCtIC6f+AZxHWoy1BbAkMLEzJ0LrpncwR+tm0CAoHReA7e/2sMqSwNP5eAZzS270qXJ8qO03JW0GnEL7d7sqkkYCIwFWKu7HLWgSaq11UxL6ei+nawIX5db22sAKwA75/Drg02V1S8/MGwBT8nFngmZQXdTsJOD/ch9/pXDU2sAh2f5/0a5A2dHz+sN5Q+eSTslHtv+dW1Kvd+JX0Etoa2tLejczZsDYsYVtDi5pI0l/kzRF0lRJU7uo8hZQEt4ZkM9LzK48Lq20tX0b8Kmu/LE92vZQ20MHDRrU7fcR9E5qHegrW65PAN/KsgJDgHtJ0sBfBr5KlS8ssD2pBQ2pRX1YbpkjqXKxUknUbBnav9zjbe9D0qP5Ch8XNTsp92NuBJxdYbur99NH0qKS2oBlOqgT9Eb++U8YPjylxXAGSeL6TWBPoKspPncC2+TjbYG7yq69KWnF/D2cASBpQE7XZO4fhaCJGT58OMOHDy/czrx23eyVhcO64kjgT5I+QQq4XwduACYB91D9CzmEJB6G7b9JWpu0cUMfUst7UlnZ0cD1wG3A9Jz3V7VLwH4D+H/A8ZK+BJwMnCPpANITwBH0jJOy/QcJnZvWYr31kt5NMYJmADPzpjmy/ZikjTsrbPt+Se9LupX0hPuCpCOzVtPRwGWk7/CBucolkpYiNU5+WNSbCGrL2LFj62JHeZyzaZC0h+0rGu1HNZREzWZJWgEYbbtbP8VtbW3ef//9C/YumBe6qUE036uoJJ0PHEAat/oMsHSZgF5DGTp0qCdPnjxX3okP1L5n8ueDP/4QXISdjmx1l1ILe8yYMbVyp0NKY4LzEYe79d1sugVTzRrkM7vnQdz+wEHdrRSiZr2AF16Ak06Cww+v+cyb3PV4hu33JB1GGoeK7QSblHq1sutJb9h4pGmwfZntLfMg1h2N9ieoIe+8A5MmpbTG5OnBo/LxbNsP2I45t0HdaLoWfSsSombNy5wnrc99DqZ2NRFmvnguT++9C/gIUj98kQaDoEQE+iCoD4sCG+YXpEHT7zXOnWBBIgJ9EAA89BDssgtcd12agVNj5mHBVBDUjJr20UvaMssLTJR0s6RNJM3IcgS3SVqjovwISU/k63+VtHDOfyLf425JB1WUnSjpb2X3uFzSKpIuzueb5PtVLpYqle+2IFqwALHUUrDnniktgHlYMBUENaNmLfq8aOlYYBfbMyUtDqwGPGR7qzyP/XBg34qqJ9s+V0kLZ3vSPPvptofl2Qp3kRabzClbZvOTwJyVq/mH5JTsQwx2BV0ybdo0ANpWXBFOPLFIU2eQ9ib+I2maZShM1pmeTuXsbvn5mcpZL2rZot8JuNj2TADbb9t+oOz6g8CKndRfskrewqSFVh2xBWmxFMBywIXA3ranA0i6QklQ7abSysES+cnguCw8tp+ki3Nra4d8/Yhc925Jg8vqzCV01hEhatY7mCNq9t578NhjKS2GmXkrStl+jCRaFgR1oZaBvo0sQCbp27mr5rdl17eg+tzhn0l6nDRI9fecNyjr0TwK3FxRdqKkE/L5drSLon0ReNz202XlR9jekiSlsGcV21eSRNBOBA4lSTOUVhqenuvuBRxWVucKYDOSVk6HlGuJLLpo1V6koAkYOXJkEjZ79NE08+bRR4sy9aKkRYD7JV1FjI8FdaSWX7Y5AmS2/yzpDpLy43qSJpDkDg6Q9DNgOElVEpIswYWkoDuQJC1Q3nVzudplV+fqugE+a/s5SSuTpBD6SNrP9jmS+gInS1qPJAp1TRWfH7b9oaTHbb8KkJeRA+wjaS9St5Cr1OlIHyfoRbS1Zc28xRdP8+jXWKPzCvNIaTA2FkwFjaCWLfqxwN5K8sPQ/iPykO2tbO9q+2XbJ2dRsQtKFW1/AJwF/KT8hnmhydtU6daRtCzwSkX2ftmHLYDPkzZu2CLfu9pSYVeklJU7gKSLv19F3ebSjAhqw+KLw+abp7SGSNpO0kOSnpZ0D7B+LJgK6k3NAn3uFz8WuFbSeJIy5IU9uMU/gG2yANqg3EUzCZhle0qV8tvRvjFDyYf/kDYrOQ34AFhN0o2kbp2ecg9JwCymxS0IvPoqnHJKSmvLKGBH26sC3yd1EwZBXalpP6Ht8cD4iuzK8/LyF5QdzwZK26Wt2VnZzOskKVdsP0ea0UDugiktStmsyn02y+mwyryK6/tVqTus2nFXhNZNL2DaNDjmGNh6a1huuVreubR/AbYfKk0hDoJ60msHhGz/vetSQdBNBg+Gt98u4s7rS7ouH6v83PYuRRgM5o81N9u20S7UnF4b6IOglzCk0Q4EPWPEGZc22oWaE4G+DoSoWfMz8PXX+fGjj8Lvf1/TmTe2n6/ZzYJeQxGLs+ZnYVbIFAcBMLtPH1hmGVgo2j5B69GtQN9ADZs+kgZLGpdXqd4oaVVJx0jqVkeapCUl7dbdD0TSz5V2kAoWIN4aOBAuvxxWWaXRrgRBzeky0Jdp2HwtzzT5GvAeeX48aUXp4VWqnpyvTyZp2EBeCAVsRFpxjLJGAAAP1ElEQVRxWl52mO0vZ5slDZuFgDNJsgZbAiPoeXfTkkC3Ar2kPrZPLM2SCBYgZs+G999PaRC0GN1p0TdSw2ZjYILtV7LtV2zPWVGYnwb2zcfHSBqWnzbuzk8g3wNGAtvlJ4ZBkvaVdGt+bZjrTlFSvzxc0gWSVsv3vlTS2PySpKXzU8pYSddKGtb5Rxf0FgY8/TQssgjTb7qp68J1RNKp+bt6ekX+uvlp+nZJ63eUFwTQvUDfSA2b5Uu2e8BOwH/b3ho4HxgN/CM/SRjYJfv8VeCXuc6KwP62KxezvGZ7J+AlYH2S8uYfc16n86EVoma9ignPPst/AX8a3+Gyj7qTGyL9bW8OLCzpC2WXfw18i7RA8Ned5AVBtwL9XBo2pIVJy9CuYXMAcKKkUrAekeudTAqO75A0bKC962Y1YO0KDZthto/I55/Ni6Dm2O6AatIFZwN7SLqI9gVYJVYh6YxMAK6m/WnjCdv/rnL/h3P6Ui77WaCkI/5gJ36FqFkvY/WNNuKgyZPZ55BDGu1KORvTvvp7HKnLs8RA2y/mbsYlOskLgm71d48FrpR0he0ZVGjYlJU7Ob8oBXvbH0gqadj8vFTQtiV1R8PmTuB4SZ+y/Uq+Vr4zxAygtB3QeqQA/i/bB0hqA84jtcL75jLPAvfa3j3b6pfzO+qYrfwheTbbeZT0IxaLtlqEZfv1Y8jzzzfbYOySQEmNdQawTtm1PlWOq+VVRdJIUrcm/fr1Q6omBVVbjui6yAJr64gNB82TnSQH1jVdtugbqWGTxc5+BFws6RbgIubu278Z2LFs5SHA/vn+N5AUMl8BBkq6krQp8xhJk/LTyM/pGecCP8z6OaLzcYagF7HUv/4FX/86PPNMo10p5y2S8io5favs2uwqx9XyqlL+xLn++utjO14NeJX9Pearfld0awaLG6Rhk68/AFROpTym7HjTimsTSaJm5exQdnx+fpX7UK51MyIf/rMs7xhIs3KArW3PljQGeI6gJXhlueVg+nRYstrcgYZxJ7A/aQ+EbWmX9gZ4U9KKpIA+o5O8IGi+lbFubg2bxUhPBAsD49zNaZghahbMC7bvl/S+pFuBKcALko60fRxwNHAZ6cmytFlOtbwgaL5A38w4TTHdvNF+BAXw3HMwahQcdRSsvHKjvZmD7YMrso7L+VOpUGetlhcEEBIIQZB47z2YOrXIPWODoGFEiz4IANZeG+65p9FeBEEhRIs+CIKgxYlAHwSQum3a2lIaBC1GdN0EASSJ4pEjUxoEdWKnnXaqi50I9EEAqTV/zDGN9iJYwBgzZkxd7ETXTRAEQYsTgT4IgqDFiUAfBEHQ4kSgD4IgaHHUEwW0YN7IkszVNmdpBMuQhOMaTbP4AV378rrtHevlTL2RNB14fh6r1/PvGLY+Tre+mxHo64CkybaHNtoPaB5fmsUPaC5fehv1/OzC1rwTXTdBEAQtTgT6IAiCFicCfX0Y3WgHymgWX5rFD2guX3ob9fzswtY8En30QRAELU606IMgCFqcCPRBEAQtTgT6IAiCFicCfRAEQYsTgb6OSIrPuwmRtFijfQiCIgk9+gLJgf0A4C3gdtvPSupje3aDXaNRfkgS8HvgAmCy7Vn19qHMlz7AWcC1wI0l/xxT0bpE0nrAkoAAbE8qyM72wA+AJUgNU9veuob3Px+o+ve2/b1a2cm2Vuromu0Xammrkgj0xXIB8G/gbmC0pP1tP9OIIJuD2oHAG8Adtp9rRFCzbUlfBN5Jbul+2/+ppw9l/BZ4yvaNkpYDZtue3iw/xs2KpBuAl4BpOctAIYEeOBEYbvvlgu5/TE5/CtwF3AcMBjYtwNaxOV0BWB54CFgPeAXYrgB7c4hAXyyPA2fYfkfSR8AfJP0gB/u6Bdncir6EFOSnAP8raYTtmfWwX+ZHX9sfAVcBqwAbAu9JesT2h/X0JXMv8LCkq0miXhtL2sv209Gy7xTZ3r9Otm4HPijq5rafB5D0edsH5+wnJR1QgK3vZltXA5+3/ZGkvqT/h0KJQF8AZUFiIVJLfm/bF6V4yx8lfdP2G3V0aVvgWdv/k/1bGVgEmFnPgJaDPMDVwDdIrfqvAAdL+rXtZ+rhR9l7Xhw4Dbje9mmSvgNcImnXAluQrcAMSacCD5O7PWz/qSBbQ4A7JL2Wz217iwLs3CJpLDCV1Mou6gkFUmv+C5IeBD4PfKpAW0CsjK05kpawPaPs/HhgZWCv3G1xNHCO7Wkd3aNovySdDfzd9rWSFrX9br1s5/Plgd1snyXpGmAN4Ju2HyrSjw58OYP0Q7hO/vscB1xp+4GifemtSPqvyjzb/9cIX2qFpKWBfsBnSE9379t+qyBbKwM/Iz3VPgP81vazRdiaYzMCfW3I3SPnk1rxNwH/tH1HvnY8sBJwJ/AjYCvbr9TBp8qgtpDtWfmx9EGgL/Bl4Bdlre1a2q/8TB63fU++dhTwbeA6YDJwU5FdSV38fc4AliJ1ExwMbGv7paJ8aQUkfQJYlvbB2EIGEyWtCOwPtJXZqukgabYzvnyQV9IVtveotZ1870Ntn1J2PtJ2oXo3EehrRG7lfBE4CtiZ1Ip/wPb1+fqOwMKkYPdkwb50GGDz9e2BXwMzgINsP16QH9U+k/tt35AHh79q+5oibHfTlwdtX5evb0v6+zxl+6l6+NRbkfRzYHtgLeBp4D+2ty3I1i3Ar0jf118CO9n+aQ3vvyuwG2kw9KacvRAwsNabzUhaElgauJjUyBGpsXW27RiM7SXcTwoks4Ergc2AwZImk3aQuaOOg5/fIc32KQW1HSUtY3tsvt6H1E/4rYL7xTv6TO7N9sdB3aZ6dubLMsA99R6c7sXsYnsTSRNtD5N0eYG2Ztu+WdLRtsflH5laMp703XiBdhXJD0kzYWrNlsDXSI2MX5AC/YfAHwqwNRexgKd2PJVfW5J+pe8A+gPrkgLJInX05f6cloLaXcB6kpaVtJ7tG4EhdRj87OgzWZ80B3tRgDpNZWymv09vpzQL5l1JWwCfK9DWFEmfBMZLmgC8V+P7r0raju/IfO8fAAdRwACp7WvzzJsv5u6nA4FDbBc+6yYCfY2w/T7wV9Io+s6kp6XXgC2ASbZfraM7HQW1DYDlJC1re3rRTnTxmdxaz8+kyf4+vZ2Dcx/9ocDuwGFFGbJ9iO33bR8DfB3YpcYmfkf7j8fFpLGrm0hdnzVF0h/yuNmLkvYjPU3cIOngrurOLxHoa4jt54DzSEHkd6SpgxcWMdDZhR+dBbUJtl/rpHqtfXmOJvhMms2XXs7DwK6k2Uo/Bd4sypCk7SWNk/QQaUzpjBqb+MD27Dzr5lO2L7c9njQDp9asbXtGHkM7HBiWB4B3L8DWXMRgbEFIGgj0sV2vneSr+bACsA1poGkQ8ONGDjQ2w2fSjL70NiRdCtwCfCf31Y8rcDD2NmAr0qysrSpnx9Tg/jeQFixtDLxo+9eSFgJus71RrexkW3cCm2RbP7O9a86/1fbmtbRVSQzGFoTtwlo5PfDhJeDC/GVueFBrhs+kRDP50gsZZPsPkgqZfljBR/nl3BJWje//LdLkhTuBi3JeG3B8je1Aeoq8izQ2NQJA0mokLaxCiRZ9EAQ9QtKfSUJwhwAnkRa/7VOQreGkMYC1SatWT7X9tyJstTIR6IMg6BF5Fsy+pOD7OGml9/sF2hNpZtTr9ZLraDUi0AdB0C1UR5ldSR1q59RyZaykxW2/Xav7NSvRRx8EQXd5jqT4+UhZnkjCZrWWJRhE0qAfk19FBeNrga0BJJ1n+/sF2ZmDpB2dpLFXA34CXO6C9PxLRKAPgqC7fIm0snMd4DHgKtuTizBke2elnb9KffQLAxcV3D//2QLvXc5hpI1u/gc4hzRl9AtFGoxAHwRBt7B9L6lFj6QvAedImmx7v4LsvSPpWeBl0o/LoALMrCLpV6Qnk9Jxyf4vC7AHsHjuBvvI9p2S/l2QnTlEoA+CoFtIWgb4KmltxvvAccDYTivNu63TSYqv9wD/BzyZ82uti1QuuTyuhvftjBOAUcCoPLB9V9EGYzA2aBok7Q0MI20Iso/twnYWCnqOpA+BB4AJJL2bOcGj1q3frGtTjsnjAbVcMJVtLUFaIfuepE2BgcCNrvGuZ5JWKT9l7s+vUN2pCPRBzZH0sO1156P+YaRdn56ooVvBfCJpy46u2b6lnr7UCkknkLa07Au8SNpu8y1gA9vfqLGt0kbkKksh/XjVXGO/nOi6CZoKSYsDq5Mf1YPmobcG8y7YzPbmeX+EqaUGSpUnilpwKLCwyzYdkvQpCtwTt0SImi1gSBomaYykayQ9KGkXSatIul3SRKV9M0tlT855t0n6Yq57Ztn1h3PaV9Kfc9njcl4fSedJmpTz1+jAnymSfitpfFZEPBk4MhbGBHViVl6QBTBdib4F2ToHWKwirz9wbkH25hCBfsGkjbQ59zCSpscw0lS5YSRFRyR9Geib83YFftPJ/b4KvJHL3lSW967TRs6HV6svaVFgBeD03O/6C9IUt1GSVp2fNxjUnvy0VW+bO+Z0NUlnZf37WnNzfkGSDi5qUHZp2/8sz7D9NGlMoFCi62bBZLLtWcBbkt4BLgeOknQRSYL2N6QNOXaQNDHXWZyywaOyVhCkrpbSzIG7c7pGKc/2PaXALelHJFnWB7Pda22/mMsdVeP3GdSWui8uouA557a3qpYvadla2ai476K23y07r2zhF0IE+gWTIfnxdLH8mm37CABJ/5B0DfAocF1Z/sKkPUI/ne9R/s9W2uTkEtJ2faW8rYBLJH2RtNs9ts8Ezsz3/CF1mFoWFEK9FhcVOudc0lIkBUtI6pUrAUeStHW2r6Ut0r63YyVdSFobsAKwF2lP3EKJQL9g8gpwDe1f6p0lHZivvQw8Y/tJSZvnFr1J85l/TurTnEDSIy9xLbB7LlsK3H8FdpF0a64/soofGwB/rOUbCwqlEYuLip5zfhlwIWlry9uAycAxtms+GSDvffsgabXv+qT/tW/UQz48plcuYEgaBuxu+0eN9iXoXdRzemW95pxLutn2Nvn4QWBwK04EiBZ9EATdwvYtHS0uKsDcL+hgzjm1FVBbtewpZUng2NLwU4FPKXUnWvRBEHSLOi8uGkgHc85ruTtYKy4Cq0a06IMg6C71XFx0DvDfFXn9gbOB3WplpJWCeWdEoA+CoLt8bHERxa3FqTrnPLf0gx4SgT4Igp5wM+395eOLNNSoOeetSAT6IAi6RZ0XFzVsznkrEoOxQRB0i84WF9mu9eIiJC1NmnO+PCnYj63HnPNWJAJ9EATdQtLfaV9cNJK0uOg3RSwuCmpLBPogCLrFgrK4qBWJPvogCLrLArG4qBWJFn0QBN1iQVlc1IpEoA+CIGhxYuORIAiCFicCfRAEQYsTgT4IgqDFiUAfBEHQ4vx/H4riuFstbn4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x180 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_results(loo, summary, chance_loo, vars_=None, labels=None):\n",
    "    _loo = loo.sort_values('LOO')\n",
    "    _loo[r'pseudo-$r^2$'] = 1 - _loo.LOO / chance_loo\n",
    "    _loo['r2_err'] = _loo.LOO_se / chance_loo\n",
    "    fig, axes = plt.subplots(1, 2, figsize=fig_size, gridspec_kw=dict(wspace=0.5))\n",
    "    ax = axes[0]\n",
    "    min_r2 = np.max(_loo[r'pseudo-$r^2$'])\n",
    "    with sns.axes_style('ticks'):\n",
    "        plt.sca(ax)\n",
    "        ax.barh(range(len(_loo[r'pseudo-$r^2$'].values)), \n",
    "                 _loo[r'pseudo-$r^2$'].values, color='grey', align='center')\n",
    "        ax.set_yticks(range(len(_loo.LOO.values)))\n",
    "        ax.set_yticklabels(_loo.index.tolist())\n",
    "        ax.errorbar(y = range(len(_loo.LOO.values)), \n",
    "                    x = _loo[r'pseudo-$r^2$'].values,\n",
    "                    xerr = _loo['r2_err'].values, linestyle='None',\n",
    "                    color='k',\n",
    "                   )\n",
    "        xlb, xub = plt.gca().get_xlim()\n",
    "        ylb, yub = plt.gca().get_ylim()\n",
    "\n",
    "        ax.plot([min_r2, min_r2],[ylb, yub], 'r:')\n",
    "        ax.set_xticks([0., 0.2, 0.4, 0.6])\n",
    "        sns.despine()\n",
    "\n",
    "\n",
    "        ax.set_xlabel(r'pseudo-$r^2$')\n",
    "        plt.xticks(rotation=45)\n",
    "\n",
    "        #### second plot\n",
    "        ax = axes[1]\n",
    "        plt.sca(ax)\n",
    "\n",
    "        \n",
    "        if vars_ is None:\n",
    "            vars_  = ['mu_beta_rbf_mean', 'mu_beta_cls_mean', \n",
    "                      'mu_beta_rbf_stdv', 'mu_beta_cls_stdv',]\n",
    "        if labels is None:\n",
    "            labels = ['RBF Mean', 'Cls Mean', 'RBF Std', 'Cls Std']\n",
    "\n",
    "        y = summary.loc[vars_, 'mean'].values\n",
    "        plt.bar(range(4), y, color='skyblue')\n",
    "\n",
    "        for ii, b in enumerate(vars_):\n",
    "            ylb = summary.loc[b, 'hpd_2.5']\n",
    "            yub = summary.loc[b, 'hpd_97.5']\n",
    "            ax.plot([ii,ii], [ylb, yub], 'k', linewidth=2)\n",
    "        plt.xticks(range(len(y)))\n",
    "        ax.set_xticklabels(labels)\n",
    "\n",
    "        ax.axhline(y=0, color='k', linewidth=1)\n",
    "\n",
    "        plt.ylabel('Parameter Estimate')\n",
    "        plt.xticks(rotation=90)\n",
    "        ax.xaxis.label.set_visible(False)\n",
    "        return fig\n",
    "\n",
    "_ = plot_results(loo, summary, chance_loo)\n",
    "plt.savefig('fig_model_fits_exp_linear.pdf', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Experiment Change Point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chance: 163444.105176\n"
     ]
    },
    {
     "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>LOO</th>\n",
       "      <th>dLOO</th>\n",
       "      <th>LOO_se</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Sticky-Choice</th>\n",
       "      <td>102276.618072</td>\n",
       "      <td>23921.019071</td>\n",
       "      <td>486.410436</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lin-Reg</th>\n",
       "      <td>101425.613110</td>\n",
       "      <td>23070.014110</td>\n",
       "      <td>484.558314</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF</th>\n",
       "      <td>87189.447312</td>\n",
       "      <td>8833.848311</td>\n",
       "      <td>473.710920</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bayesian-GP</th>\n",
       "      <td>83465.549982</td>\n",
       "      <td>5109.950982</td>\n",
       "      <td>462.814536</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kalman</th>\n",
       "      <td>82540.859890</td>\n",
       "      <td>4185.260889</td>\n",
       "      <td>485.059184</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Clustering</th>\n",
       "      <td>81779.400685</td>\n",
       "      <td>3423.801684</td>\n",
       "      <td>481.802894</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Linear/Clustering</th>\n",
       "      <td>81325.362616</td>\n",
       "      <td>2969.763616</td>\n",
       "      <td>479.926131</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF/Kalman</th>\n",
       "      <td>79682.879552</td>\n",
       "      <td>1327.280551</td>\n",
       "      <td>481.747704</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF/Clustering</th>\n",
       "      <td>78355.599000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>476.086496</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                             LOO          dLOO      LOO_se\n",
       "Sticky-Choice      102276.618072  23921.019071  486.410436\n",
       "Lin-Reg            101425.613110  23070.014110  484.558314\n",
       "GP-RBF              87189.447312   8833.848311  473.710920\n",
       "Bayesian-GP         83465.549982   5109.950982  462.814536\n",
       "Kalman              82540.859890   4185.260889  485.059184\n",
       "Clustering          81779.400685   3423.801684  481.802894\n",
       "Linear/Clustering   81325.362616   2969.763616  479.926131\n",
       "GP-RBF/Kalman       79682.879552   1327.280551  481.747704\n",
       "GP-RBF/Clustering   78355.599000      0.000000  476.086496"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loo = pd.read_pickle('Data/model_fits/model_fits_exp_cp.pkl')\n",
    "summary = pd.read_pickle('Data/model_fits/model_params_exp_cp_rbf_cls.pkl')\n",
    "\n",
    "# vars_  = ['mu_beta_lin_mean', 'mu_beta_cls_mean', \n",
    "#           'mu_beta_lin_stdv', 'mu_beta_cls_stdv',]\n",
    "# labels = ['Lin Mean', 'Cls Mean', 'Lin Std', 'Cls Std']\n",
    "            \n",
    "# calculate chance\n",
    "n = len(pd.read_csv('Data/exp_changepoint/changerbfpred.csv')) # n_subj * n_trails_per_subj\n",
    "n_arms = 8\n",
    "chance_loo = -2 * n * np.log(1. / n_arms)\n",
    "\n",
    "print \"Chance: {}\".format(chance_loo)\n",
    "\n",
    "loo.sort_values('LOO', ascending=False, inplace=True)\n",
    "loo['dLOO'] = loo.LOO - np.min(loo.LOO)\n",
    "loo['LOO dLOO LOO_se'.split()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAADKCAYAAACrHYtRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXe4XFXZvu8noZcYSiiHIgKCKMWQoNQQpBOlCaICGpUmICAin3wgTRAElQ4aQfgApUiRkoAYSIBQhNA7SP8RSgBJQHry/P5Ya5LJyTknp8zsPWfOe1/XXGtm7b3X+86cOe+svda7niXbBEEQBM1Lv7IdCIIgCOpLBPogCIImJwJ9EARBkxOBPgiCoMmJQB8EQdDkRKAPgiBociLQB0EQNDkR6IMgCJqcCPRBEARNTgT6Athqq60MxKP3Ppqa+H726keniEBfAG+++WbZLgRBu8T3s/mJQB8EQdDkRKAPgiBociLQB0EQNDkR6IMgaHhGjBjBiBEjynaj1zJX2Q4EQd35z3/g5pth001hkUXK9iboBmPGjCnbhV5N9OiD5ue552DnnVMZBH2Q6NEHzc9aa8Fbb8GAAWV7EgSlEIE+aH7mmgsWXbRsL4I2OPGBruXwd/b8XwxevDvuNC0xdBM0Py+8AD/4QSp7OZJOkXS7pNNa1f9R0h2SJkhasyz/gsYkAn3Q/Hz4ITzxRCp7MZLWBha0vREwj6R1qg6faHsD4AfAUaU4GDQsMXRTAJMmTeKYY44p242G46ijCopHX/gC3H13Mbbqy3rA2Px8LLAucC+A7edz/SfAtOJdCxqZCPRB0HsYCDybn08BvtTGOScAp8+pIUl7AXsBLL/88rXyr26suuFmZbvQqyls6EbSAEmjJY2XdJekoZJ2rDp+RjvXDZd0XBdtLS7p0mxrgqStO9uOpC9L+lFX7AUNzkMPwZJLprJ38w5QSR0akF/PQNJBwOO2J8ypIdujbA+1PXTQoEG197TGjDz9EkaefknZbvRaihyj/x5wle3hwEbAR8CMQG/7JzW0dQZwVrb1NVr9Q3SE7Qdtn1dDX4KyWWIJ2H//VPZu7gI2zc83A2aMR0naAlgf6FKnKOgbFBno3wfWk7S47U+B7wCb5173IEkTACRtkLMHxknapXJxviO4VtIX87H+uf5KSUtUndcfWNr27QC2P7Z9Vz68hqTrcvsLSZpL0iWSbsvlXNU9f0k/lnR3trdq9vPa/PrsQj61JmXq1KlMmjSpGGNLLw2//GUqezG27wc+lHQ7MB14SdLh+fAZwOeAcZL+WJaPQWNSZKC/CHiJ9EUcC1wL/NP2cNuTq847EdjO9ibA33Ldwvn6w2w/DtwMbCJpADC37Teqrh8EVLc3C7a/AYwh9Yx2IN3qDgMeA75ZOS//eOwMbJB9eQb4BXBCfv2upPXasyNpL0kTJU18//335/jh9DXuu+8+Ro0aVYyx//4X7rsvlb0c2wfa3sj2/rZfs318rl/V9jr5/2nvsv0MGovCAr3tT2wfa3sN4DzgoA7OfTOX03PVzsCDth/Lr/8K7EIK1FfnnvZ4SeNJQb69QcdHc/kKaWJrJeD+XDcRWLnq3M8B99ueVuXLasCJ2c6mQEsH72HGGOgCCyzQ3ml9liFDhrDXXnsVY+zJJ2Ho0FQGQR+kyMnYz0qaO798A3gX6N/GqZa0WL6m4t/5wLKStgew/RwpyH6LNO4/OfdkhufA/KqkjXIbc0tat9J2tUvAc8CQ/HooMzMayMcGV3zI5VPAwdnOUOCabn0YAQMGDKClpd3fydryhS/AxImpDII+SJFDN18GJuTe8C9Ik0aLSrpCUvX69MOA6ySNI/XkIQXovYHdJG2c68YAH9ue0oatnwD7ZVvjSL33trga+JKk24A1gCsrB/Jw0pXAndmXzwO/Bo6QdEseflquKx9AUBILLghDhqQyCPogheXR276G2XvAW1Y93zCfdwcpe6Ca8bncqbpJ0hBOW7beBL7dxqHx+fgFVXW7tHFO5byzgdaTrju0ZTNoYF59Fc49F/bYo9dPyAZBd+iVC6Yk7UsKuNuU7UtnaGlpKW4VaDA7b7wBZ54J224bgT7ok/TKQN9OTzsI2matteD118v2IghKo1cG+t5GaN20T9zpBEH9CfXKoPl58klYd91Irwz6LBHog+ZnvvlgtdVSWQKSNpf0fUn9Ja1WihNBn6apAn1r4bL2hNLauG68pFuz2FpBq3iCwlhhBTj//FQWjKRzSNll++U1HnNUlgyCWtNUgb41XRRK2xTYABhZH2+C0vj0U3j77VQWz6q2DwEq+gttLRIMgrrS1IG+SijtAkmnZcnijmb/5gXmy9dI0jl5cdRoSYtImjeLmt2YZZBHFvA2mpbChM0eeggWW6wsmeLpWTfJeWFgKb82Qd+mqQN9K8bb3pD2c+9vBiYxM23z68BLtr8GnAnsA2wP3Gl7K+A/HRkLUbM5U5iw2Yorwt/+lsriORi4FFgduBz4WRlOBH2bvpReWRE0+wBA0s2k2+jKCtpNgS+Q5BnOJQmYfVvSlqTP6S7SatyH8/kPdmTM9ihgFEBLS4s7OrevUpiw2SKLwE47zfm8+rBo7iwA0JHiaRDUi77Uo58l2NreNIuTvVZV9yjQL2dGPAVcmM/ZEPhf4HmSJg7AmgX53bQUJmz21ltw8cWpLJ4jW73+aRlOBH2bZuzR71qlVtkdzgV+DBwInC7pllx/KvB34G+S/gG8R9qIOWh0XngBdt89KVgutlghJiV9B/gusKaka3P13B1cEgR1o6kCve3xJB351vUjq54Pb+P48KrntwCV4D5b1o6kHW1/mtPmnuupz0EBDB6cNh0pNo9+DHAnqcNwWq77GHit3SuCoE40VaAviNGSFgL+XbVFYYeEqFnJ9OsHBW/+kuWzpwAHS1qStH+CgKWZudlNEBRCXxqjrwm2t7S9ge3vl+1L0Emeew6++91UFoykY0n7GowGziJtldnTNk+RdLuk01rVHy5pUvWiwSCA6NEXQoia1Y9O3Sl98glMmpTK4tnS9lfzJjib0M4eCp1F0trAgrY3yus81rF9bz58Lmm4aNMeeRw0HdGjD5qfVVeF8eNTWTxTczmNtMn9F3vY3nrA2Px8LDAj8cD267TKLgsCiEAfBPXmPEnzk7K2bqWHPXrStpiVH48pwCLdaaR6Qd/kyZO75ciIESMYMWJEt64NiqXpAr2kjbNswXhJN0taX9IUSeOyBMIqrc4fKempfPzvkubJ9U/lNv4l6YBW546XdEMZ7y/oBg88AJ/5TCoLxvaltj+wfZ3twbZ/08Mm3wEG5OcD8uvu+DXK9lDbQwcNGtQtR8aMGcOYMWO6dW1QLE0V6CUtDhwDbJ9TJrcnrYR9xPYmpOXnh7Zx6cn5+ERgi1w3ObexLrBrq3OH2966Pu8i6AqTJk2as17O0kvD4YeXso2gpO1zx+AhSQ9LenjOV3XIXcwcg98MuLuH7QV9gKYK9CQdm4ttTwWw/a7t6m7cg8CyHVw/sI26eejGwqjQuimGUaNGzVkvZ6ml4NBDU1k8xwK72V7L9pq2e7Si2vb9wIeSbgemAy9JOhxA0o+A35EWDZ7VU8eD5qHZsm5agEcAJH0X2JdZezzDSNIGrfm5pENIomaH5bpBOVNiOeDiVufuBtxl+zDaIbRuiqFTWjnvvZeUK9daCxZaqP5OzcozwKu1bND2ga2qjs/15wHn1cLGiQ+8WfNzfzF48e66E/SQZgv0r5KCPbb/KulO4GhgDUnjSOOZ+0r6OTACuCBfdzJwIXAFsCjwOnnoRpKAyyR9tnKu7XMLej/BHOiUVs5TT8GGGyYJhCFD6u/UrFwIPCPpsfzatrct2omgb9NsgX4McIWky/PKxMr7q4zRVzg5P6hoytv+ON/u/pSkYEmut6R3aXtYJ+gNfPGL8MgjsNJKZVj/NfA94OUyjNeTVTfcrGwXgk7SVIHe9mRJxwDXSJpO2uThRNI4aWf4J3CcpHmZOXTTD3jC9kOSBtfD76DOzD8/rL56Wdafsj2hLOP1ZOTpl5TtQtBJmirQw2yiZBVav64+/4Kq59OBdfLL2VbXVJ/bFULrpmReeQXOPhv23ReWWaZo6wtKuoOkb2MA2wcU7UTQt2m6QB8Es/H220mP/tvfLiPQ/7pog0HPqCwCGz16dMme1I4I9EHzs8Ya8OKLhZqU9HnbzwDvFmo46DHNuAgsAn0BhKhZcTTQENkOwEnMvqeBgR8W707Ql4lAHzQ/jz8Ou+4Kf/lLysApANsn5acX2h5XqZe0fiEOBEEVpa2MlTRc0ot5efgdeZ/WWrZ/Rg+unV/SWVX6OGfl+qPzUvY7JMXYa29hwQVh3XVTWTy/bPX6oDKcCPo2ZffoL7J9RO7l7EPadq0m2J5tG8AucCRp5et+AJKGVR37me2xkkZLWs520+VHNx2f/Sycc06hJmPP2MajK6t9u3J+b1jx2yhaNwOAqZIOk3RrVowcLGlpSZcCSJqrslG3pD3yDju3S1pb0qL5zmCcpNPzORNy+YN8bKKkLXLdBZJOy731tgZ1N7A9Q/bA9m1tnPMwUHgKR9AxbYqcffxx2njk44+LdGUMsD9pZexP8uOHJD2mICiUsgP97pJuA84HLgdOs70xSS3yENuvkvKQFyYp9o3NCpXbknRrtiP1vtcGxufVr63vCi7LKpSbAodU1Y+3vSFt/+NpxhPpJkmPSlq2qq4/8BU62Bw8RM3KoU2Rs0ceSWmVjzxSmB+2p9h+EfhbLj8B9gN6JGoWBN2hUYZuliRtg3a9pF1JqnwVIbCrSAH9a8BxwIrAWsC4qnZuBYZJ+itwA3BR1bEtJR1ICt5LVNU/mssPACTdDPQHvl3toO0tJF3AzM/qd8B/gEttv9HeGwtRs3JoU+RspZXg2mvLkkA4ARhO6pDcBpwJbFSGI0HfpexAX+Fd0vDNvsBgYCXgT/nYlcBfgLltP5d1Z+61vROApLmB/raPzK8fZNZAfxiwMTAvcEdV/SzB1/aMfTYl3Slpd9uVdqo/p5/ZHkvQkLQpcjZwIHzjG8U7k6jcHS6Qhfb2KMuRoHM0o4ZP2YF+d0kbAvOReuvbkXo9M8bEbU+V9CF5n8ysZzM6D/lMI8kb3JqzYOZm5n6aFa7P7d1D53fjOQb4ff6n/JgkNfta995iUDpvvgl//ztsvz0sXvjE2TOS7gVOzJ0SzemCoFyaUcOntEBvezzw2VbV17dz+iekMfzKteeTxvWr2bBV+xvm8lhaiZrZHln1fHgbvr1PygJqzdHt+Bc0Mi++CHvuCYMHFx7obe8haS7bn0rqR9r1LAgKpewe/RyRNAp4I0/M9kpC1Kxk1l4bPv0U+hWXeyDptKoNQn4OnGB7uqTTgJGFORIElJ91M0ds72U7FpkE3UeC/v1TWRzV2TWbVz1vfRcbBHWn4QN9EPSYZ5+FHXdMZXG0l2kVGVhB4TT80E0zEKJmxTLbMNn06fD++6ksjlXy4j21ev75njYs6RRgKHB/9f6xklYH/pDt/Nj2wz21FTQHEeiD5ufzn4cbbyza6q5Vz69o53mXkbQ2sKDtjSSdI2kd2/fmw78CvkNah3I2KYstCJor0EsaDmyWF2F9k6Q1snPeOapyzgrAcbZ3K8XJoE9g+9Y6Nb0eM1OIxwLrApVAv2hFe0nSZzpqRNJewF4Ac889Nypg/uKwultoflutsTs3EthUgb6CpA1Iy81HVAf5oI9y//2w/vpw550pA6d3MxCoTDZMAb5UdaxfO89no3rl9tChQz1x4sRa+tjn6aqAWmfoiXhaM07Grgr8FtgZ+FTSzZJuk3Rl1qiZgaS7JJ2RpYe/LekqSQ9LWisfPy2LrN0uafmqa86U9KCkrQp/d8EcmU3YbNll4aSTUlkwkpaY81ld4h3SKnJyWb0IcHo7z4M+TjMG+i2Af9h+C/gU+LrtYcATJL2cahYjrcjdGjiVNNSzD/CjfPywLLJ2DLB31TXHAiOq6mYjRM3KYzZhsyWWgAMOSGXx/KXG7d1FEugD2Ay4u+rY25KWldRC6u0HAdCcQzfnAOvn3vYEYJSkZYAlSVIGz1SdO9n26wCSnrX9oaRJwCL5+KGSNiVJKzxRdc0b+ZqB7TkRomblMZuw2dSpcO+9sM46MGBA2xfVjwclfYMUkKcB2H67u43Zvl/Sh5JuBx4CXpJ0uO3jgaOAS0lZN/v13PWgWWjGQP8pSYHyJuBq4Gnb35V0PLPrjLid55K0GDA8ZzdszswsilnOq63rQS2YTdjsmWdgs81g4kQYMqRod4bmRwUz+51ll6hOqcwcn+sfppUUSBBAcwZ6bL8taXfgOqCfpKGkW9lnOr5yFv4DvJc3O4l85N7Ml76Ugn0JY/R5jwQkLWT7vcIdCAKaLNBnobTx+fkTwMrtnLpbPmdG76dKBO2FynHg623YqL5meI+dDurPfPPByu19FeqLpG1IczwL5En+UbZ/WIozQZ+lqQJ9oxKiZiXz8stw6qlw0EGw3HJFWz+cNJxyY1awXKFoB4KgGbNugmBWpkyB665LZfFMs/0RM+d2+nd0chDUg+jRB83P6qvD00+XZf0qSX8Dlpd0EWnHtCAolAj0BRCiZsXTKENltk+V9A/SCtYnbD9Wtk9B3yOGboLm57HHYLXVUlkwkm60/YTtK2w/Juniwp0I+jy9okefFy0dScpbf4vk909t/7sT164ArGj7lk7aOpW0Afi0bjscNBYDBsBWWxW6WErSF0i9+GUk7Zir5yY2HglKoOEDvaTFSUH+67bflbQKcEYXmliBtEBljoFeUr/YzaoJWW45OOWUoq0uCawOLASskes+Bn5WtCNB0PCBnqQpc5HtdwFsPy3pVQBJRwMTbI+VdAFp8+7BJOXQ90jiZrsBG0haz/amko4kBf7pQCWf+XzSncIYSd8jaYgcASxL+qF4IW/yvBLwV+A10j/wj3LefdBgVETNWlpa4KOPYNIkaGmBeectxH6WKb41bxIyHVjG9lOFGA+CVvSGMfqlga5sDP5N4Fu2vwbcQNKbuSgH+TVI/3DDSVogFSnpJYBdbP+5VVuP2d6MlDExEDgEOBDYMV/TLiFqVi6zCJs9+iisuGIqi2cnYAxwuaT+kq4vw4mgb9MbevSvAi3tHGtLd+Y44AhJc5E1QKpYDRguaXxV2wAPtTMmX4kMk4DPAJ8DHrY9TVKHUSNEzcplFmGzz38ebroplcXzQ2Aj4Jb8vVmgDCeCvk1vCPRjgCskXZrH6Fcm9fIh6dcsrbQ9TmUDhhfzMMv6wMHAxcxcpPIUcJPtnwBImhtYhva1u1v/kDwPrCFpIrNu+BA0GLMImw0YAJtvXpYrle9WLJgKSqPhA73tyZJ+BVyfA/rbpEktgKtICpXbkkTIAI6WtC5pDP1npF75CZIus72LpNdyj97AJSSVy87yW5K++BvZ3ic9enNBMbzxBlx+OXzrW2Vo0p8FjANWlnRTfh0EhdLwgR7A9lhm7pPZmta6s//TxjnDqto6ntmHdHarOj48Pz26qm4kgKS5bK+bd6q6A3h9zt4HpfPyy/DTn8J66xUe6G1flgP8SsBzPdGiD4Lu0isCfQOxoqRzgQWB82x/2pmLQtSsZIYMgU/KufnKk/jbAosCwyRh+/elOBP0WSLQdwHbT1N1dxAEneAGYDTwctmOBH2XCPQFEFo35TDjLuqZZ+DAA+G008rIvHnX9nFFGw2CanpDHn0Q9AwJ5p47lcVzoaRfSfqmpB2r5BC6jKSFJV0n6Y68sK/18WskvSNps565HDQb0aMPmp+VV4ZrrinL+veBF0l7GUPK9rqqm23tScoUuwwYl1OOP646vg+wd3cdDZqXLvfoJQ2XdFzV667oznTWxjJZXAxJB0u6TdIESaflugld9HfFLpxf8/cT9Gmm2d7D9jH5cWwP2loPGJsX9z0ErFp90HZXVpAHfYgeD91UFh/VAkkVfzYH/ilpa9KXeeO8V2t3ekLDgU4F+ixqVrP3E5TH1KlTZ+jdcN99adjmvvvKcOVlSftIGiJpbUlr96CtgcDU/HwKsEh3G6qW6Jg8eXIPXAp6Az0O9JXetaQLJJ2We95H5bqVJd0k6VZJR+S6H0gan79kW1RdeyZwY252OHAr8C3gZNuGGUJR1bbHZ6kDKrIGkn6dxzDHSVoeGAn8TtLvJM0v6RJJt0i6TNLckkbm56OBNavez3hJJ0i6V9KPct22ku6TNGpOdxWhdVMu991330ytm+WXhz/+MZXFMxfwVWB/4Ce57BBJS+XvX/XjUuAdoKK1PCC/7ha2R9keanvooEGDuttM0Euo9Rj9eNsHSvoXcAxpYdKPbL+cA+yywGW2z5f0GeBvzFyZeoft/fPq1wVsvyepq4JmABsAG9mentu6gJkKlz8BrrV9iaQfkwSnAN6xvQuAZp2wu5y0cOqfwHmkxVjDSD2r8R05EVo35TJkyJCZejeDBkG19k2B2P5BN655jdTZmQVJBwObSroc+DJJ0iMI5kits24qQl8f5HJV4KLc216NpCuzZX59LbBc1bWV++q1SOOP0LGgGbQtanYS8H95jL+1gNRqwEHZ/veZqUDZ3j39o3lj54peyTTb/7X9CvBmB34FJTNgwICZejdTpsCYMaVsDi5pXUk3SHpI0sOSHu5Bc+cCuwK3A3+2/ZGkrSSNyLZOB74HnCSpnF+2oCGpdY++dc/1KeAg269m2QADdwMbA/OSZAQqVILpFqQeNKQe9SGS9rFtScNs31Z1TUXU7ANgqVx3i+3Rkv4X+DpJj6Za1Oxm21fCDFGzXemcqBlAPyX1wYHA4u1+CkFj8e9/w4gRMHFiWiVbLKeTJDb+COwLfLe7DdmeSvpOV9fdWPX8AOCA7rYfNC/dDfS7KgmHzYnDgT9LmpcUcL8JXA/cBtxD22OMQ0jiYdi+QdJqpA0c+pF63tWBfhRwHTABqMwo/V0zpWB3Bv4f8GtJXwVOBv4kaV/SHcBhdI2Tsv0HCZ2b3sMaayS9m+IFzQCm5s1yZPsJSeuV4UTQt+lyoLc9nqTL3rp+ZNXz4bn8N7B1q1OPzY82rwWutD296tjvgd+3On/DXI4hyRhXH9uylb3XSHcQFXZrdfxf7bQ9vPX7AcbYvlbSMuTx96AXMM88sOyyZVl/WdL8wP2SriTWrgQl0HBfOtuXl+1DB+yUJ3EXpAu3yCFqVjIvvQQnnQSHHlpo5k1OBjjd9geSDiHNP8UEalA4IYHQBWxfanvjnJZ2Z9n+BJ3kvffgtttSWSA5Lfi4/Hy67QdsR65tUDgN16NvRkLUrFyOOuooeLgnyS494oWc1ns3MA3A9v1lORP0TSLQB0F9WQBYOz8gZXL9sDx3gr5IBPqg+XnkEdh2W7j22pSBUyDdWTAVBLWmpmP0kjbO8gLjJd0saX1JU7IcwQRJq7Q6f6Skp/Lxv0uaJ9c/ldv4l6QDWp07XtINVW1cJmlFSRfn1+vn9lovlqqc32lBtKBJWGQR2GWXVBZMjRdMBUG3qFmgl7Q4SfZg+5yOuD1phewjtjchbdR9aBuXnpyPTyQtlgKYnNtYl7Sgqfrc4ba3zjbnA2asXM0/JL8DvhWTXgFkcbN+/eDEE8tKsTwdOJC0qf0uQGl6yUHfpZY9+m2Ai/PqPWy/a/uBquMPAh39pw1so24e0kKr9hhGWiwFsCRwIbCb7ckAki5XElS7SdKA6gvzncHxWXhsT0kX517Xlvn4Yfnaf0kaXHXNLEJn7RGiZo3Bfffdx5/POgueeAI++GDOF9SeqXkLStl+giQ1HASFUstA30IWIJP03TxU89uq48NoO4f455KeJE1W/SPXDcp6NI8DN7c6d7ykE/LrzZkpivYV4Enbz1adP9L2xiQphV3asH0FSQTtRNIdx9bAfvnYafnaXYFDqq65HNiQpJXTLtXqgAss0OYoUlAAQ4YMYe+NNoIvfhEef7wMF2LBVFA6tfzSzRAgs/1XSXeSlB/XkDSOJHewr6SfAyNIqpKQZAkuJAXdRUnSApNtD88LTi6T9NnKubbPrbL5OdsvSFqBJIXQT9Ketv+UtXVOlrQGSdL16jZ8ftT2J5KetP06gKTKQO7ukio6OG7jmvb0cYIGYsCAAQzaYIOUR7/KKnO+oMZUJmNjwVRQJrXs0Y8BdlOSH4aZPyKP2N7E9g62X7VdGWe/oHJh3g7tLOCn1Q3mBSfv0sawjqQlSPIG1eyZfRhGknFd0Paw3HZbG4a6VUnVefuSpGL3bHVtSA73NhZeGDbaKJUFIWlzSY9IelbSPcCasWAqaM0FB3yHCw74Tt3t1CzQ53HxY4BrJN0CnE3qqXeWf5K0tuclD91Iug341PZDbZy/OTC2lQ8fkTYrORX4GFhZ0o2kYZ2ucg9JwCzS43o7r78Ov/tdKovjOGAr2ysBPyINDwbBLDw1YSxPTRg75xN7iPLmTb2OPGl6V2Xyt5EZOnSoJ06cWLYbfZcHHoBhw9LwzeDB3WmhrbvBji+QbrH9tfZeNxLx/aw9Jz7Que0qDls77e51wv1z3s7xF4PbVEbv1Hez104M2f7HnM8KAlJwf/fdoq2uKena/FzVr21vW7QzQd+m1wb6IGhwCt/hJAjaIwJ9AYSoWbks+uab/OTxx+GccwrLvLH9YiGGgqAThExx0PRM79cPFl8c5urd/RpJC0u6TtIdkr7XxrGbJd0m6XpJxaUYBQ1PpwJ9iRo2/SQNljQ2r1K9UdJKko6WtFknfR8oacfOfiCSfqG0g1TQJLyz6KJw2WWw4oplu9JT9gQuIS0+3KPyf5X5hLQqfBhJZmFk8e4FjcocA33JGjZzAWeSvsAbk768Xe2WDQQ6Fegl9bN9ou1XumgjaGSmT4cPP0xl72Y9YKztacBDwKqVA7Y/tP1qfvkpWfs+CKBzQXM2DRvggbRoFaivhs16wDjbr2XbrwGvVWxLGgnMZftcSUcD40n586cA/wUuBlYGNs+SCjsD2zFTvuBA2/dLegh4BHhU0hdIOdAbAlsClQVgI0grd68g/dB9ApyS99ANGpilX3sN5p+fyTfcwKCttirbnZ4wEKikE08BZpPjlLQQsBez79Xc+ry98nksX+D2ikE5dGbopkxLAJTqAAAOyklEQVQNm6UrtrvANsD/5Jzl80mbeP8z30kY2Db7vB1wZL5mWWBv260XtbxhexvgFWBNYA/gj7luHjogRM0ah3cGDuTq7bbjz7fcUrYrnULSUvn/ofpxKUlGpCLONyC/rr5OwJ+Bw22/QwdUazENGjSoHm8jaCA606MvU8Pms0BHO0W0JV1wNnBEVpc8A3ij6pwVSXoj41q185Tt/7bR/qO5fIXUm/ocSVMH0p1M+47Zo0g/MrS0tPTOVWlNwgcLLMBXzz6br5btSCfJd67DW9dLOpi0evxyksRH6w7WscAdtnvHL1oT087iptk4rIvnd5fO9OjL1LC5CxguaanKMUmrVp0+hdTrh5k/CP+xvS/wP6S5hU+A/vnY88C92c/hpDsHyHr2bdD6h+T5KjtrtnNN0GDM98EHtNx9Ny3zz1+2Kz3lXNLc1u3An21/JGkrSSMktZC+8zvkO4Afl+pp0FDMsUdve7KkiobNdNJEz4mk3kNn+CdwXLWGDekH5gnbDylrvVcxQ8PG9seS9gculjQ38CFQ/QW+mTTss072C2DvnGWzEPAb0o/GopKuII1Jjs4aOtOAW4BfdfJ9QPpHu1LSD0g/Dh3NMwQNwiL/+Q9885swcSIM6b3rmPI82ddb1d1Y9bLD4cSg79JwWjeNrGEjqR+A7emSRgN7dSZDp6WlxXvvvXfd/QvaRtOmceQBB8DAgd3Npe+y1k1vIrRuyqOSWNKDONw7tW4aXMNmIdIdwTykNLdOpWG2tLRw1FFH1dezIAiCdoiVsV3A9lTbG9n+qu3Dy/Yn6CQvvAB77JHKIOiDRKAPmp8PPoCHHy5rz9ggKJ2GG7oJgpqz2mpwzz1lexEEpRE9+iAIgiYnAn3Q/Dz8MLS0pDII+iAxdBM0P4svDnvtlcogaCC22WabQuxEoA+an5YWOProsr0IgtkYPXp0IXZi6CYIgqDJiUAfBEHQ5ESgD4IgaHIi0AdBEDQ5DSdq1oxIepe2N2cpm8WBN8t2ohWN6NN8tlcv24l6IWky8GI3Ly/y7xW2ZudN23PcNi2yborhKdtDy3aiNZImNppfjepT2T7UE9vd3mKqyL9X2Oo+MXQTBEHQ5ESgD4IgaHIi0BfDqLIdaIdG9Ct86l0U+dmErW4Sk7FBEARNTvTogyAImpwI9EEQBE1OBPogCIImJwJ9EARBkxOBvkAkxefdC5G0UNk+BEFPiJWxdSQH9n2Bd4A7bD8vqZ/t6SW7NoNG8UeSgHOAC4CJtj8t16MZf7+zgGuAG3OdHKlqSFoDGAgIwPZtdbKzBbAP8BlSx9S2v1bD9s8H2vx72v5hrexkW8u3d8z2S7W01ZoI9PXlAuC/wL+AUZL2tv1cmcE1B6/9gLeAO22/0AjBy7YlfQV4L7mp+21/VKZPwG+BZ2zfKGlJYLrtyY3y41gWkq4HXgEm5SoDdQn0wInACNuv1qn9o3N5MHA3cB8wGNigDraOyeUywNLAI8AawGvA5nWwN4MI9PXlSeB02+9Jmgb8QdI+OdgXHlxzr/kvpCD/EPB7SSNtTy3Sjzb86m97GnAlsCKwNvCBpMdsf1Kia/cCj0q6iiT6tZ6kXW0/2wg/jiUi23sXZOsO4ON6NW77RQBJX7Z9YK5+WtK+dbD1g2zrKuDLtqdJ6k/63teVCPR1oCoIzEXqye9m+6IUZ/mjpG/bfqsE1zYDnrf9v9nPFYD5gallBq4c5AGuAnYm9eq/Dhwo6Ve2nyvSn6rPYmHgVOA626dK+h7wF0k71LGH2RuYIukU4FHysIftP9fJ1hDgTklv5Ne2PawOdm6VNAZ4mNTLrtcdCqTe/DqSHgS+DCxVR1tArIytOZI+Y3tK1etfAysAu+bhiaOAP9me1F4bRfkn6WzgH7avkbSA7ffL8iW/XhrY0fZZkq4GVgG+bfuREn06nfQD+aX89zseuML2A0X51GhI+n7rOtv/V4YvtULSYsDcwGdJd28f2n6nTrZWAH5Ount9Dvit7efrYWuGzQj0tSEPi5xP6sXfBPzb9p352K+B5YG7gP2BTWy/VqBvrYPXXLY/zbenDwL9ga2BX1b1ruvpT+vP6knb9+RjRwDfBa4FJgI3FTG0NIe/3+nAIqRhhAOBzWy/Um+fGhlJ8wJLMHMyti6TiZKWBfYGWqps1XSSNNu5pXqSV9Lltr9Vazu57Z/Z/l3V671s11XvJgJ9jci9nK8ARwDfIPXiH7B9XT6+FTAPKag9XZBP7QbUfHwL4FfAFOAA208W5Fdbn9X9tq/Pk8Xb2b66CF/m4NODtq/Nxzcj/f2esf1Mkb41GpJ+AWwBfAF4FvjI9mZ1snUrcCzpe3oksI3tg2vY/g7AjqTJ0Jty9VzAop3Z0KOLtgYCiwEXkzozInWyzrYdk7G9hPtJgWI6cAWwITA4b1qxOCnDpehJz++Rsn4qwWsrSYvbHpOP9yONF36n4HHw9j6re7M/Y6Hw1M+OfFocuKfsSesGYlvb60sab3u4pMvqaGu67ZslHWV7bP6RqSW3kP72LzFTRfITUiZMrdkY2J7UifglKdB/AvyhDrZmIRbw1I5n8mNj0q/0ncCCwOqkQDF/CT7dn8tK8LobWEPSEpLWsH0jMKToyU7a/6zWJOVmLwBQcApjI/79GpVKFsz7koYBX6yjrYckzQfcImkc8EGN21+JtB3f4bntfYADqMMEqe1rcubNV/Lw037AQbbrnnUTgb5G2P4Q+DtpFv0bpLulN4BhwG22Xy/BrfaC11rAkpKWsD25aKfm8FndXsZn1aB/v0blwDxG/zNgJ+CQehmyfZDtD20fDXwT2LbGJs5g5o/HxaQ5q5tIQ541RdIf8nzZy5L2JN1NXC/pwDld21Mi0NcQ2y8A55GCxBmkFMELi5jgbMefjoLXONtvdHB5vX17gQb6rBrVpwblUWAHUjbSwcDb9TIkaQtJYyU9QppLOr3GJj62PT1n3Sxl+zLbt5AycGrNaran5LmzQ4HheQJ4pzrYmoWYjK0TkhYF+tkuaif5jnxZBtiUNOE0CPhJI00oNtJnVaERfWoUJF0C3Ap8L4/Vj63jZOwEYBNS9tUmrbNjatD+9aQFS+sBL9v+laS5gAm2162VnWzrLmD9bOvntnfI9bfb3qiWtloTk7F1wnbdejldJacCXpi/1A0XvBrps6rQiD41EINs/0FSXdIPWzEtP5x7wqpx+98hJS3cBVyU61qAX9fYDqS7xLtJc1AjASStTNLCqivRow+CoEtI+itJ6O0g4CTSIrfd62RrBGkOYDXSqtVTbN9QD1vNTAT6IAi6RM6C2YMUfJ8krfT+sI72RMp8erMsmY7eTgT6IAg6hQqU2ZXUrnZOLVfGSlrY9ru1aq9RiTH6IAg6ywskRc/HqupEEjartSzBIJIG/ej8qFcwvgb4GoCk82z/qE52ZiBpKyfp65WBnwKXuU56/hUi0AdB0Fm+SlrZ+SXgCeBK2xPrYcj2N5R29qqM0c8DXFTn8fnP1bHtag4hbWTzv8CfSCmj69TTYAT6IAg6he17ST16JH0V+JOkibb3rJO99yQ9D7xK+nEZVAczK0o6lnRnUnlesX9kHewBLJyHwabZvkvSf+tkZwYR6IMg6BSSFge2I63J+BA4HhjT4UXdt3UaSfH1HuD/gKdzfa31j6oll8fWsN2OOAE4DjguT2zfXW+DMRkbNAySdgOGkzb82N123XYWCrqOpE+AB4BxJL2bGcGj1r3frGtTjcnzAbVcMJVtfYa0QvYDSRsAiwI3usa7m0lasfols35+ddWbikAf1BxJj9pevQfXH0La1empGroV9BBJG7d3zPatRfpSKySdQNq6sj/wMmmbzXeAtWzvXGNblY3IVVVC+vGqucZ+NTF0EzQUkhYGPk++VQ8ah94azOfAhrY3yvsgPFzpoLRxR1ELfgbM46pNhyQtRR33xK0QomZ9DEnDJY2WdLWkByVtK2lFSXdIGq+0b2bl3JNz3QRJX8nXnll1/NFc9pf013zu8bmun6TzJN2W61dpx5+HJP1W0i1ZEfFk4PBYGBMUxKd5QRbAZCX618nWn4CFWtUtCJxbJ3sziEDfN2khbcI9nKTpMZyUKjecpNiIpK2B/rluB+A3HbS3HfBWPvemqrr3nTZyPrSt6yUtACwDnJbHXX9JSnE7TtJKPXmDQe3Jd1tF29wqlytLOivr39eam/MDknRwvSZlF7P97+oK28+S5gTqSgzd9E0m2v4UeEfSe8BlwBGSLiJJ0P6GtOHGlpLG52sWpmryqKoXBGmopZI58K9crlKps31PJXBL2p8ky/pgtnuN7ZfzeUfU+H0GtaXwxUXUOefc9iZt1UtaolY2WrW7gO33q1637uHXhQj0fZMh+fZ0ofyYbvswAEn/lHQ18DhwbVX9PKQ9QpfLbVT/s1U2N/kLaTu+St0mwF8kfYW02z22zwTOzG3+mAJSy4K6UNTiorrmnEtahKRgCUm9cnngcJK2zha1tEXa93aMpAtJawOWAXYl7YlbVyLQ901eA65m5pf6G5L2y8deBZ6z/bSkjXKP3qR85l+QxjTHkfTIK1wD7JTPrQTuvwPbSro9X79XG36sBfyxlm8sqCtlLC6qd875pcCFpC0sJwATgaNt1zwZIO99+yBpte+apP+1nYuQDY/0yj6GpOHATrb3L9uXoHdRZHplUTnnkm62vWl+/iAwuBkTAaJHHwRBp7B9a3uLi+pg7pe0k3NObQXUVqq6SxkIHFOZfqrjXUrhRI8+CIJOUfDiokVpJ+e8lrt/NeMisLaIHn0QBJ2lyMVFfwL+p1XdgsDZwI61MtJMwbwjItAHQdBZZltcRP3W4rSZc557+kEXiUAfBEFXuJmZ4+W31NNQWTnnzUgE+iAIOkXBi4tKyzlvRmIyNgiCTtHR4iLbtV5chKTFSDnnS5OC/Zgics6bkQj0QRB0Ckn/YObior1Ii4t+U4/FRUFtiUAfBEGn6CuLi5qRGKMPgqCz9InFRc1I9OiDIOgUfWVxUTMSgT4IgqDJiY1HgiAImpwI9EEQBE1OBPogCIImJwJ9EARBk/P/ATCSOoln92poAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x180 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = plot_results(loo, summary, chance_loo)\n",
    "plt.savefig('fig_model_fits_exp_cp.pdf', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Experiment Scrambled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chance: 102308.523851\n"
     ]
    },
    {
     "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>LOO</th>\n",
       "      <th>dLOO</th>\n",
       "      <th>LOO_se</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Sticky-Choice</th>\n",
       "      <td>61623.540831</td>\n",
       "      <td>17825.538563</td>\n",
       "      <td>392.600549</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lin-Reg</th>\n",
       "      <td>59843.580073</td>\n",
       "      <td>16045.577805</td>\n",
       "      <td>389.632792</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF</th>\n",
       "      <td>47515.709879</td>\n",
       "      <td>3717.707611</td>\n",
       "      <td>384.791200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bayesian-GP</th>\n",
       "      <td>47438.683217</td>\n",
       "      <td>3640.680950</td>\n",
       "      <td>383.125824</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Clustering</th>\n",
       "      <td>45595.940702</td>\n",
       "      <td>1797.938434</td>\n",
       "      <td>378.727403</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Linear/Clustering</th>\n",
       "      <td>45248.992712</td>\n",
       "      <td>1450.990444</td>\n",
       "      <td>377.569572</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kalman</th>\n",
       "      <td>44355.167170</td>\n",
       "      <td>557.164902</td>\n",
       "      <td>375.373361</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF/Clustering</th>\n",
       "      <td>43858.175944</td>\n",
       "      <td>60.173677</td>\n",
       "      <td>377.763555</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF/Kalman</th>\n",
       "      <td>43798.002267</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>378.292263</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            LOO          dLOO      LOO_se\n",
       "Sticky-Choice      61623.540831  17825.538563  392.600549\n",
       "Lin-Reg            59843.580073  16045.577805  389.632792\n",
       "GP-RBF             47515.709879   3717.707611  384.791200\n",
       "Bayesian-GP        47438.683217   3640.680950  383.125824\n",
       "Clustering         45595.940702   1797.938434  378.727403\n",
       "Linear/Clustering  45248.992712   1450.990444  377.569572\n",
       "Kalman             44355.167170    557.164902  375.373361\n",
       "GP-RBF/Clustering  43858.175944     60.173677  377.763555\n",
       "GP-RBF/Kalman      43798.002267      0.000000  378.292263"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loo = pd.read_pickle('Data/model_fits/model_fits_exp_scram.pkl')\n",
    "summary = pd.read_pickle('Data/model_fits/model_params_exp_scram_rbf_cls.pkl')\n",
    "loo.rename(index={'Kalman w/SC':'Simple Kalman'}, inplace=True)\n",
    "\n",
    "#note the \"Simple Kalman\" is modeled without sticky choice and wasn't included in the paper\n",
    "loo.drop('Simple Kalman', inplace=True)\n",
    "loo.drop('Scrambled', inplace=True)\n",
    "\n",
    "# calculate chance\n",
    "n = len(pd.read_csv('Data/exp_scrambled/gprbfscrambled.csv')) # n_subj * n_trails_per_subj\n",
    "n_arms = 8\n",
    "chance_loo = -2 * n * np.log(1. / n_arms)\n",
    "\n",
    "print \"Chance: {}\".format(chance_loo)\n",
    "\n",
    "loo.sort_values('LOO', ascending=False, inplace=True)\n",
    "loo['dLOO'] = loo.LOO - np.min(loo.LOO)\n",
    "loo['LOO dLOO LOO_se'.split()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAADKCAYAAACrHYtRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXmYVMXV/z9fNhUUFUFlNGrUiBrR4BCDioCCK2rcNeKCiaJxw1eN0bgnuETjq+KWoEZf930HFdlEcAVEcUMjavyJC4gBUQwg5/dHVUPT9Mz0dN9epud8nqef2123btXpnp7TdatOfY/MDMdxHKd6aVFuAxzHcZzi4o7ecRynynFH7ziOU+W4o3ccx6ly3NE7juNUOe7oHcdxqhx39I7jOFWOO3rHcZwqxx294zhOleOOvgTsscceBvij6T6qGv9+NulHTrijLwGzZ88utwmOUyf+/ax+3NE7juNUOe7oHcdxqhx39I7jOFWOO3rHcSqe/v37079//3Kb0WRpVW4DHKci+OYbGD0a+vaFNdcstzVOBiNGjCi3CU0aH9E7DsCMGXDwweHoOFWGO3rHAdhmG/j663CsECRdLekFSddmlP9D0kRJEyRtHctqJI2R9KKkfuWx2KlU3NE7DkCrVtChQzhWAJK2BdqZ2U5AG0m/TDt9uZntCBwDXBjLzgbOA3aLR8dZijt6xwH4+GM45phwrAy2B0bF56OAHqkTZvZRfLoI+DE+3xp4yczmA99KWq2+xiUNkjRJ0qRZs2Yla7lTcbijdxyAH36Ad98Nx8pgDWBefD4XyLZCfBkwND5vaWbWQP2lmNkwM+tuZt07deqUhL1OBVMZ96lVzsyZM7n44ovLbUbFc+GFFzZcqVhsvjm8/HL5+l+R/wDt4/P28fVSJJ0GvGNmE2LRj2mnV6jvNG98RO84lclLQN/4vB+w9FdI0m7ADsCQtPpvStpeUjugvZnNw3EiJRvRS2oP3Au0A1YCTgE2MLNH4vnrzOyULNf1AfqZWc4LTJI6AtcD6xLe4yXAglzakfQLoNbMbs21P6cKeOMN2G03GDmyIiJvzGyKpB8kvQC8Afxb0rlmdglwHWFaZ6yk6WZ2PHAFcAewCssWaCuey19vnKBarvXP7tYxH3OqllJO3RwFPGJmt0pqBWwBHAA8ApDNyRfAdcANZvaCpDZALeHHpUHMbCowNUFbnKbA2mvDySeHY4VgZoMzii6J5V2y1P1/wC6lsMtpepRy6uZ7YHtJHc1sMfAbYFdJ4yR1kjQBQNKOMUZ4rKRDUxdLai/pCUlbxnMtY/nDktZOq9cS6GxmLwCY2UIzeyme7irpydj+qpJaSbpX0vh4bCWpj6Qhsa3fS3o59tcl2vlEfH1jST61ZsK8efOYOXNm+Qzo3BnOPz8cHafKKKWjvxP4N+F2cxTwBPCcmfUxs/T4rsuBX5vZzsCDsWy1eP05ZvYOMBrYOU4HtTazr9Ku7wTUGS9mZvsAIwjzn/sTFrR6AW8DB6bqxR+Pg4Edoy0fEGKVL4uvv5W0fV39pIevff/99w1+OM2dyZMnM2zYsPIZ8N13MHlyODpOlVEyR29mi8zsz2bWFbgVOK2eurPjcUksOhiYamZvx9f3AIcSHPWjcaQ9TtI4gpOvK17srXj8jBC+tgkwJZZNAjZNq/tTYIqZ/ZhmyxbA5bGfvkBNPe9hafha27Zt66rmRGpraxk0aFD5DHjvPejePRwdp8oomaOXtKGk1vHlV8C3QMssVU3SWvGalH23AetL2g/AzGYQnOwhhHn/WfHOoE90zJ9L2im20VpSarNJeuotATMI8/cA3YEP087PALqlbIjH6cDpsZ/uwON5fRjOCrRv356amjp/N4vP5pvDpEnh6DhVRimnbn4BTIij4bMJoWEdJD0kqUNavXOAJyWNJYzkITjo44EjJPWOZSOAhWY2N0tfpwAnxb7GEkbv2XgU+Lmk8UBX4OHUiTid9DDwYrTlZ8ClwHlRU2QU8JPGfABOBdOuHdTWhqNTcXTp2Y8uPV3CJ1+0bDNd00LSicAsM3uwwcplpqamxo4//vhym1HxlHXD1Oefwy23wLHHZluQVTlMKhXdu3e3SZMmlaXvxoZX5kozCq/M6bvZJHfGRie/P7BXuW3JhZqamvI6MadhvvoKrr8e9t3XI2+cqqNJOnozuxHw8EYnObbZBr78stxWOE5RaJKOvqnhWjeNx++AHCc5XOvGcSCEVfboUVHhlfUkHjlX0szUxr5YdrukV2KY8eGlt9apZNzROw7AyivDFluEYwFI2lXS0ZJaStqigHbqSzxyCzAgy2UDYujvPfn261QnVeXo0+UL4uvrcrxunKTnJb0kqYy7dpyysdFGcNtt4Zgnkm4CdgdOivs5hjZwSX3Ul3jkS5bfE0J8fUeU+NiwgH6dKqSqHH0mjRRK6wvsCAwsjjVORbN4McyZE47508XMzgRSOgrZNgTmSi6JR9I5w8x2AP4KXNVQ455hqnlR1Y4+TSjtdknXxmTK9a3yrQSsHK+RpJvi5qjhktaUtFIUNXtG0n2SBpbgbTQ75s2bx+TJk0srcvbGG7DWWuGYP0uiRpLFTYCF/GrUm3gkEzObE48TCPLc9eIZppoXVe3oMxhnZj2pO/Z+NDCTZWGbewP/NrNdCNr2JwD7AS+a2R7AN/V15qJm+TN58mS6d+9eWpGzjTeGBx8Mx/w5HbgP2Ap4ADijgLbqTDySjSjwh6QueHYpJ4PmFF6ZEjRbACBpNOHW+rBY3hfYnCDPcAtBwOwwSbsTPqeXCPOgb8b69WrWm9kwYBiEnbGJvYtmQG1tLVdddRWdS7lxac014aCDCm2lQxwYAFCfumlD1Jd4RNLvgBMJEiJrmtlJwN2S1iR8R39f4Ptwqozm5OiXc7ZmlhotISlV9pakFjFaYjpwh5ldFeu0JiRK6UrQ2dkaeK00pjcv2rdvT21tbcMVk+Trr+Hpp2HPPcMUTn5cAIxLe/0/hAFCXtSTeORWggJset198u3HqX6q0dEPSFOrzIdbCCOiwcBQSWNi+TXAY8CDkp4F5gOLCrLUqRw+/hiOPDIoWDbS0Uv6DXA4sLWkJ2Jx63oucZySUlWO3szGEXTkM8sHpj3vk+V8n7TnY4CUc8+Ww/YAM1scQ+lmFGqzUyF06xaSjuQXRz8CeJEwOEhtbloIfJGQdY5TEFXl6EvEcEmrAv9KS1FYLy5q1gRo0QLyTBATpbLnAqdLWoeQK0FAZ5YltnGcsuGOvpGY2e7ltsEpAjNmwHnnwZAheUfeSPozIUH3xsAnhOQ6uyVnpOPkhzv6EuCiZpVDnXdWixbBzJnhmD+7m9mvYsKbnQkpLx2n7LijdxyALl1g3LhCW0ntZP2RkNB+y0IbdJwkaE4bphyn2NwqaRVChNbz+IjeqRCqztFL6h1lC8ZJGi1pB0lzJY2NEgibZdQfKGl6PP+YpDaxfHps4xVJp2bUHSfp6XK8P6dIvP46rL56OOaJmd1nZgvM7Ekz62Zmf03QQsfJm6py9JI6AhcD+8WQyf0IO2GnmdnOhC3pZ2W59Mp4fhLLFs9mxTZ6sLwk7JVRCnbP4rwLp1jUq6HTuTOce25BaQQl7RcHAW9IelPSmw1f5TjFp6ocPUHH5i4zmwdgZt+aWfoQbSqwfj3Xr5GlrA15bIxyrZvKo14NnXXXhbPOCsf8+TNwhJltY2Zbm9nWhTTmOElRbY6+BvgcQNLhcarmb2nnexGkDTL5g6T3gG2BZ2NZpxg98Q5B8Cy97jhJl9VnSLo6YNs847OdZKmtrWXSpEkMGpQl5cD8+TBxYjjmzwfE719SNDLL1FbxOz9Rkv/IOEupNkf/OcHZE7PsHAF0BLpKGksQgrpcUspZD4zXXUnQrpkPdIhlqambTYEt0pI5pKZuzinFG3KSI6WhU1NTs+LJ6dOhZ89wzJ87gA9i8o8n0+QQ8iKPLFN/AX4DHBKfF5X+/fvTv3//YnfjJEC1hVeOAB6S9EDcrZh6f6k5+hRXxgcpZ29mCyXdQBCiOjtV0cxM0rdkn9ZxqoUtt4Rp02CTTQpp5VLgKODTZIzKmmXqNQhZprRiqsIOZvYpgKTVE7KhTkaMGFHsLpyEqKoRvZnNIizGPh7FyG4kjLJy5Tmgr6SViFM3ksYDi82soIwUToWzyiqw1VbhmD/TzWyCmX2SehRoVWOzTLWo4/kKeIap5kW1jegzRclSZL5Or3972vMlQOr2uEt9dRuDa900AT77DG68EU48EdZbL99W2kmaSNC3MQAzO7UAqxqVZQpYUsfzFUjPl9C9e3fPl1DlVJ2jd5y8mDMH7roLDjusEEd/aZImEbTsjydkq+oH3N5A/TmS1ic4+bkJ2+I0YdzROw5A167wSX4zLZJ+ZmYfEETMEiOPLFMXElIZCjipkL4vf3124nXP7tYxX3OcAnFHXwJc1KxyKNIU2v7AFayYv8CA3xbScCOzTL0J9CykP6c6cUfvOADvvAMDBsDdd4cInEZgZlfEp3eY2dhUuaQdkjTRcfKlbFE3kvpI+iRGtkzMEipWaPvXFXDtKpJuSNPHuSGWXxS3t0+UlPR8rFNO2rWDHj3CMX/Oz3h9WiGNOU5SlHtEf6eZnRdHPicQUrElgpmtkAawEVwAvBTnPZHUK+3cGWY2StJwST9JxS07TZwNN4SbbsrrUs8Z61Q65Xb0KdoD8ySdA+wBrExw/F8AV5vZYZJaASPNbBdJxwJHx2sHAx8DjxDmRKeZ2amSJphZT0nHxLqrAn8ys5GSbidEJdQCz5lZ5gT6juk7X81sfBab3wTWI7nNMU4JSAmarbA7duFCmD0bOnaENm0a22yzzBnbpWe/cpvg5Ei5N0wdGTck3UYIIbvWzHoTtnafaWafE2KTVwP6AqOiQuW+BN2aXxNG39sC4+Lu18y7gvujlEFf4My08nFm1pMghJaJlj6RRkp6K4atpcpaAttRT3JwFzWrTIYNG5Zd1GzatBBWOW1ao9s0s7lxc9SD8biIEPVS1XozA4fey8Ch95bbDCcHyj2iT03drEPQ7nhK0gBCHHBqE8cjBIe+CzCEkI9zG2BsWjvPA70k3QM8DdyZdm53SYMJznvttPK34nEBgKTRQEvgsHQDzWy3eAeQ+qyuAr4B7jOzr+p6Y+kbUmpqanxDSoWQVdAMgvTBE08UKoFwGdCHMPgYD1wP7FRIg46TBOV29Cm+JUzfnAh0AzYBbo7nHgbuBlqb2YyoO/OamR0EIKk10NLMLoivp7K8oz8H6A2sBExMK1/O+ZpZ39RzSS9KOtLMUu2kf05nmNkonCZJVkEzgDXWgH32KbT51J1gWzO7J04xOk7ZKbejP1JST8Kc/BDCyH18fABgZvMk/UAUdzKzWXEhdDwhN+cY4PkYBdOaZSJQKZ6K7b1Kw1vIU1wM/G/8R11IkJ+t6vnWZs/s2fDYY7DffmGePj8+kPQaQSG1NWlTgI5TTsrm6M1sHLBhRvFTdVRfRJjDT117G2FeP53lNorE+XfM7M+EhBDp5wamPe+TxbbvCYvBmVxUh31OU+eTT+C446Bbt7wdvZkdK6mVmS2W1IKQ4cxxyk65R/QNImkY8FVcmG2SuKhZE2DbbWHxYmjR+PgESdem7WD9A3CZmS2JyUIGJmil4+RFuaNuGsTMBpmZbzxxiosELVuGY+NJj67ZNe155h1rI02qM7vUCpmkJN2ukMh+nKTDC+nXqT4q3tE7Tkn48EM44IBwbDx1RVXlHW3VQHapujJJDYjZz+7Jt1+nOqn4qZtqwEXNKpelU2pLlsD334dj49lM0lDC4mv6858VYFqd2aXInknKgDskfQ2c3FDSE0mDgEEAG2ywQQFmOk0Bd/SOA/Czn8Ezz+R7dXru1ofqeN5Y1gBStxdzgZ+nncuWSeoMM5sTo9iuAg6qr3FPPNK8aBKOXlJfwiYUAV8T7P4fM/tXDtduBGwcM0/l0tc1hH+aH/M22GlWmNnzRWi2vuxSK2SSMrM58ThB0uVFsKfZkEp4Pnz48DJbkhwVP0cfJQ8uAPY2s17AH4HGiJFsRNhVm0tfLczsNHfyzZApU2DllcOxMniJINsBIbvUy2nn5khaX1INMZOUpPbx2IXc94s4WRgxYkTVJT6veEcP9CdIJXwLYGbvA5/DUtngfvH57ZI2krS/pFcljZG0F2Ee8sgocYCkC2JkwphYf6MoR/wQMDCeaxXbvkXSKEm3xGs3iZENj0saHe8WnCbKvHnzloqcsf76cMUV4ZgnktZuuFZumNkUIJVdagkxu1Q8ncok9WB8DnC3pAkEKZGzk7LDqQ6awtRNZ6AxSlMHAoeY2ceSBHwPzIiaOl2B9cysT9S/P4egT7I20M/MfpR0VFpbb8dNMCMlrUEQRRtMWBSbWp8R6Ytdq6++en1VnTIxefJkhg0bxkUXXQRrrw2nFpLHGwhSHbs2WCtH6skutUImKTMrWL/BqV6awoj+c6AOgZLlwtdSAdBDgPOiENmmGfW3APpIGgfcxLI50DfqmK5JCZ/NBFYHfgq8Geu+laX+MsPMhplZdzPr3rZt2/qqOmWitrZ2mcjZvHkwenQ45s9USftI6iSpg6QOSdjpOIXSFEb0I4CHJN1nZt9K2pQwyocwP9k5jtxTUQmfxFH4DsDpwF0EVUqA6QRN+1NgqSDaeiy/uJVO5g/JR0BXSZNYPgrCaYK0b99+mcjZBx9Av34waRLU1ubbZPf4SGHkuD7kOMWk4h19FDH7C0HCWMAcgtAYBAnjRwn69N/Esosk9SAkGjmDMPK+TNL9ZnaopC/iiN6Ae4GRjTDnb4Tb869if4sKenNO5fDznwdnX8AcfcyHgKRVzWx+YrY5ToFUvKMHiLLAdUkDZw6//pilztJUgGZ2CXGuM40j0s73iU8vSisbCBAFq3rExCMTgS8btt5pEqy8MmyaOdPXOOLi/8VAW0nbAMPM7LdJmOcUzuWvzy5K/bO75a12WrJQzibh6CuIjWMETjvgVjNbnMtFLmrWBPj0U7jmGjjtNPjJT/Jt5VzCIukzUcFyo6TMc6qTUoVxNoXF2IrBzN43s15mVmtmN5bbHidB5s6FJ58Mx/z50cz+y7K1nZb1VXacUuEjescB2GoreP/9Qlt5RNKDwAaS7iRkR3OcsuOOvgS4qFnlkuSUmpldI+lZQkTWu2b2dmKNO04B+NSN4wC8/TZssUU45omkZ8zsXTN7yMzelnRXghY6Tt402tFL6iNpSNrr65I1CSStF8XFkHS6pPEx0cK1sWxCI+3duBH1E38/ThOgfXvYY49wbCSSNpd0ILCepAPi41AKTDwS225M8pEVyhwHEpi6SW0+SoIoKraEsI38OUl7Al2A3mZmknrn0WwfYAIwI8f+E3s/ThPiJz+Bq6/O9+p1gK0Ieze6xrKFhH0ceZOefETSTZJ+aWYpTfpU8pElwI3Ar+socxpJl579ym1C4hTs6CVNMLOeUXJgLiGu/TkzuzjuYr0RWCmWDZF0DHA04Z/iT2Y2Ml47H9gM2I3gnE8GrgMuMTODFeVg48anfjGUbVzUsLkU6E34RzuakLNzf0mjgPOAfxL+MWcR4ucHAHtGe86VdH18P+MICoL9gL+b2a2S9iWISE0GtkwlIHeaLilRs5q11oKZM6GmBlZaqVFtxO/l85KuJjjZ9cxsegLmNTb5SLYyp5EMHHpvuU1InKQXY8eZ2WBJrxA2jlwC/M7MPpV0r6T1gfvN7Lb4RXyQZTtTJ5rZyXH3a1szmy+pM1GpshHsCOwUkzMLuB2YYGajJJ0CPGFm90r6PcuSM/zHzA4F0PI5Qx8gbJx6DriVsBmrFyEpxLj6jHBRs6bBsGHDALhon32ge/dCJRAOIgwsVo+j8cfNbO8CzGts8pFsZVlJ/362bt0683tfFM4peg9Nt698P/84Bm6QpB19SuhrQTx2Ae6Mb2INgq7MLyUNJmjHpMu6To7HbYA34vOUoNkHdfSXTdTsCuD/FFKqnZtRfwugVtLxwMoECYS5aX2v8H7MbJGklBbOj2b2HfCdpHq3zaVn8KmpqfEMPhXKUlGzVVeFkSNDpqn8+S2wEzAmKqEWqmbXqOQjdZRlJTPD1KRJkwqz1FmOXHfVnrNtJwAumzKrwbqF7MBN2tFnOrTpwGlm9nmUDTBCAoXehOmciWl1U1/M3QgjaAgj6jMlnRDn6HuZ2fi0a1KiZguAdWPZGDMbLulPwN4EPZp0UbPRZvYwLBU1G0BuomYALeI/7xpA/p+6UzEsFTUD2LVgheHU9yipDVMvAccT/g/6Ee5OU8yJd8hLiMlH6ihznLwd/YAoHNYQ5wL/lLQSweEeCDwFjAdeJXsmnFqCeBhm9nTUjX9eUgvCyDvd0Q8DniQstqZ+Eh9LG0kdDPw/4FJJvwKuBG6WdCLhDqCxd11XxP6n4jo31cVXX8EDD8AhhwRt+vy4ARgLbCppZHydN2Y2RVIq+cgbxOQjUa8plXxEwEnxkmxljoNyneMpFZIOMbMHym1HNqKo2WJJ6xEEq/rncl1NTY0df/zxRbbOyYelG6YmT4YePeDll7PN0ec8gSppTWATQrKbOUnZWUx86iZ5Sjh1k9N3s+J2xlaqk48cFBdx2wE5pyNyUbMmQG0tLCpMdVohC9m+QAeglyTM7H+TMM9xCqHiHH0lY2b3EW6NHScbTwPDgU/LbYjjpOOOvgS41k3lc+Hhh8PgwXDttYVE3nxrZkMaruY4pcUdveMASNC6dTjmzx0xG9pUYuSNmT2ShHmOUwju6B0HQnapxx8vtJWjgU+AVEIaI6S7dJyyUlWOXlIfgiTCeVFk6nDg4Kifk6qzETDEzI7I2ojj5M+PZnZsuY1wnEyqUqZY0o6EOOIj0p2842Rj3rx5vHvXXWHaZnJdm6Rz4lNJJ0iqlbRtlEFwnLJTjY6+C2HD1cHAYkmjo8zxw3F37lIkvSTpOklvSDpM0iOS3lRI7IykayU9H2ViN0i75npJUyXtUZcRkgZJmiRp0vfff1/Et+sUyuTJk+l15JE8uffesMEGhTTVCvgVQZDvlHh0nLJTVVM3kd2Aq83s6yhqtreZLVDQ0N+F5XVz1gKGELaqTwE2ArYFfkeIkz/HzL6X1I+wFf3ceM2fgdbA9cAz2YxwrZumQ21tLVdddRWdO3eGTp3ybsfMjknQLMdJjGp09DcBO8TR9gRgWNzJug7Byac7+llm9iWApA/N7AdJM4E14/mzJPUlOPV30675Kl6zRvHfjlNs2rdvT+2mm8LEidCuHeSpNhplQS4kCPEJwMw8AYhTdqpx6mYxcBhhpD4YeN/MehMSNWfGzlkdzyVpLaCPme0EnJ92bTbFTKep869/Qf/+4Zg/QwnfuTnAoUDeYTySVpP0ZMwWdVSW8wMkvSjpKUntY9l0SePiY8t8+3aqj2oc0WNmcyQdSRA8ayGpO0HNry6542x8A8yXNAZ4swhmOpVE167w6aeFCJoBzDOz9yXJzN6VtH0BbR1HkNG+Hxgr6T4zWwhLVVdPIORGOJAwrXgl4W6zTyFvwKlOqsrRm9k4YkIQM3sX2LSOqkfEOkszRKWem9nHqfMEmePMPtKv6VOw0U5l0KYNrL9+oa18KmkVYIqkhyns/2t74KSoa/8GIchgWjy3GTAtCuyNIq4FAR0kjSdMMw42sx8K6N+pIqrK0VcqLmrWBPj3v+GKK+Css/KKvIkL/0Pjwv+ZhAQ6haQTXAOYF5/PZdm6UX3nesa72T8RskcNrcfepRmmNigs0shpAlTjHL3jNJ7582H8+HDMg5jXeEh8vsTMXjezBuNqJa2bNq+eetxH/dmlsp5Lk0V+lJCsvD57h5lZdzPr3qmASCOnaeAj+hLgomZNgwvfLHgp5uOYl/hl4EcIyUPqu8DMvgD6ZJZLOh3oK+kB4Bcsf3fwPrBV3BfSD3hZUhtCfon/EvImf5jZptN8cUfvOMnRlrAPI7Uj1gh5ZPPhFuAewsarYWb23xgy3DKmyrwZeIEQNHA4YfrmaUnzY5lLfDhLySnDlKTehPjgFoSRyvkE7e0phBjz35rZ+2n1BxLS9M0kzCEeYmYLJU0nJPxeBbjbzIam1f0cWGBme8Y27gd+Q5jrvDL2s4AgbXAkMMHMRuVg+xrALrmqCEo6G7jTzD7LpX4ueIapymftL7/k908/DU88ESJwlqeqw2g9w1T5UFRLLSDTX07fzQbn6CV1BC4G9otRJvsRHO40M9sZOAM4K8ulV8bzkwi7VWFZ+FcPQlLu9Lp90pz8ysB/CXcc1xM0a3oDA2n8XcgawAG5VJTUwswuT9LJO02DBSuvDIceCmuu2XDlOpDUQ9LTUVLjTUkelutUBLk4zb2Au8xsHoCZfQu8rmW63VOB+uLSsu0ebUNIFl4XvQi7WrcHxsZ5zNR85hepvuPdQCszu0XSRYTQyoXA1cB3wF2EEMtdJY0j6N/8miAnCyEEbUoMX5sGvCVpc8KiWk9gdyC1TbI/IUXcQ4QfukUEqYVx9bwPp4nwmcTMU0+lpqamkGaGEqZM/gGcSJhScZyyk0vUTQ1hWgVJh0uaIOlvaed7kT2M7A+S3iPMVz4byzpFh/sOMDqj7jhJl8XXuwIjgc6pvhvBXsAfzWwX4DZCjPFz8U7CCDk9exEc/gXxmvWB483s8oy2vjKzvYDPgK2BY4F/xLI29RnhomZNi2mvvsojl1wCCxYU0sy8OIWpuI+jkA1TjpMYuTj6zwnOHjO7hzBi6Qh0lTSWMHK5XFLKWQ+M111JcI7zCSNhWDZ1symwhaQNU3Xj1M058fVP48alpX3XQTY5ghuBQyTdCfwyo/7GhDn/sYSEEKm7jelm9l2W9t+Kx89i3Z+ybJfs1HrsWi58rW3btvVVdSqAPTfckJNvvBHeeaeQZpLcMOU4iZGLox8BHCEpNYWR+vJOM7OdzWx/M/vczFLO+vbUhXHL9g3A/6Q3GGOOvyXLtI6ktYEv4suXgD6S1k2dk9QlrfpcwqgfILWC9o2ZnQj8kbC2sIigTgnwEfBatLMP4c4BoC7N+swfko/S+nGxqirihw02CHH0m22WdxtmdoyZLQDOJEz/7ZWUfY5TCA06ejObRXCYj0fdlxuBOxrRx3NQ1p19AAARgElEQVSEeOCViFM3cZv2YjN7I0v9XYFRse+FBE3vuyQ9D9zJ8nP7o4E9JD2RVnZ8bP8p4HbCj0YHSQ8RIoaGK+jTjwXObsT7gBDy9ntJzxAcf33rDE4TYuFKK8FOO8FqqzX6Wkm7Spom6UNJrwJb57phynFKQU7hlaVE0u7AS6nF30pCUgsIOx8lDQcG5RKh4+GVlU+7+fM5c9114YgjYJ11Mk/XG8Im6RXgADP7TFJXwlRknUlpKg0PrywfpQqvrLg5RDN7tuFaZWNVwh1BG2BUrmGYrnXTBHj9dejVC3bZJZujb4jvUt8FM5sWvx+OUzFUnKOvZOJdxk7ltsMpAt26wbff5nv11mnTh0p/bWb7JmGe4xSCO3rHKZzachvgOPXhjr4EuKhZ5dNh9mxOeecduOmmRkfemNknSdsjaTWC1k0Hwt6NOzLOPw70Bg5KSYFIGkCQCJkDHF6J61xOeXCZYscBlrRoAR07QquKGfukMkz1Ao7NMu9/AnBN6kVG1qk7CVmnHAdwR+84APynQwe4/37YeONym5Jie8KC/49AKsPUUswsc8f40qxThPDkHiWx0mkSJOroJfWWNCbGyo+WtIOkuZLGRumEzTLqD4wJjcdKeiw1aklLcvyKpFMz6o6T9HRaG/dL2ljSXfH1DrG9rNtRJU1I8j07VcKSJfDDD+FYGdSXYarg+ukSHbNmzSrIUKfySczRl1nlckl8vRlwFUEW2TerODnT/sMPYZVVmDVyZEn7zTPDVDYaVd8zTDUvkhzRr6ByaWavp50vpsolwDqEHbtHxN28SHpA0vOSRkpqn35h/Ie6JI5qjpN0V5SX3T2ePyde+4qkbmnXXCbpNUm/q8cuFzVrYoz96COOBv45ZkxJ+zWzL1KSHGmPwwjyH30VskhlZpjKxgpZp4prudOUSNLRl1PlEmA74D0zS0+hNjDq2D8AHJql74cIadcuJ9xx7EmIWgC4Nl47gKBdkuIBgoTx0dSDi5o1LX7WowenTprEkaedVm5TUtxC+O69APwzlWFKUn8ASUOBo4ArJA0ys0VAKuvU0QSpZMcBkg2vXE7lUtKLwEUsU7n8D3CipD8QtN1vj9ddSRiJP0QIJfuSOHWjsD/4/gyVy1vS+vypmX0saSPgSaCFpOPM7OY4srkybklvT0iYnMlbZrZI0ntm9iWApNTc5pExXG0Jy4ubpa6pmMlcp3DWbt2a2k8+qZjF2HhnvHdG2TNpz08FTs04fych4sZxliPJEX05VS5THBdt6EW43W1nZr1i29k0ISzjSFq9EwlJm4/LuLayxIGcRFjzm2/gwANhxoxym+I4iZOYoy+nymWaDf8FDiHEFy8ENo1Kk9s1+g3Bq8B44Jg8rnWaGF+ssw7MmgXbbFNuUxwncSpOvTJXKlnlMhNXB2zyeHJwpyiUSr2yyW6YMrNnm4KTd5oIH38Mxx4bjo5TZTRZR+84ibJgAbz5ZqE5Yx2nIqkYYQ/HKStbbAGvvlpuKxynKPiI3nEcp8pxR+84EKZtamrC0XGqDJ+6cRwIEsWDBoWj45SIvfbaqyT9uKN3HAij+YsuKrcVTjNj+PDhJenHp24cpwKRtJqkJyVNlHRUlvOPS/qPpH5pZSkZ73GStiytxU4l4yN6x6lMUhmm7gfGSrovSoWkOIEVs0il5L0dZzl8RO84lUljM0wBdJA0XtI/Yq6GOvHEI80Ld/SOU5k0NsMUQM8o4vcJMKi+ip54pHnhUzclYPLkyfMlNZQ4otR0BGaX24gMKtWmZ8xsj2I0Lmld4L6M4i9YljHqB3LLMIWZzYlPHyVDCbY+Jk+ePFvSJ7nWz6CUfzPva0Vy+m66oy8N082se7mNSEfSJLepYaJNRXHyEDJMEeSwM/s9naDm+gA5ZJiK+ZYVFVx3BD6sr36GDXkP6Uv5N/O+8scdveNUJrcA9wCnAMNSGaaAlmY2PGaY2hvYV9LfgceBpyXNB74BjiiX4U7l4Y7ecSqQfDJMEdJxOs4K+GJsaRhWbgOy4DblRiXaVEmU8vPxvvKkySYecRzHcXLDR/SO4zhVjjt6x3GcKscdveM4TpXjjt5xHKfKcUdfQiT5593EkLRquW1wnELxOPoiEh37iYTt6xPN7CNJLcxsSZlNW0ol2CNJwE3A7cAkM1tcTntg6d/uBsJGpGdimczD1ACQ1JWgxyMAMxtfpH52Iyh1rk4YmJqZ7ZJg+7cBWf+mZvbbpPqJfW1Q1zkz+3eSfWXijr643A58B7wCDJN0vJnNKKdzjQ7sJOBr4EUz+7jcDszMTNJ2wPxgoqbErfzl5G/AB2b2jKR1gCVmNqsSfhjLjaSngM+AmbHIgKI4euByoH8dap1JcFE8ng68DEwGuhFkJJLm4nhcD+gMTAO6ErSNdi1Cf0txR19c3gOGmtl8ST8Cf5d0QnT2JXeuceR8N8HJvwH8r6SBcRdmWZDUMkrxPgxsTNjduUDS22a2qFx2Aa8Bb0l6hKAGub2kAWb2Ybl/GCsAmVmmFn6xmAgsbLBWnpjZJwCSfmFmg2Px+5JOLEJfx8S+HgF+YWY/SmpJ+O4XFXf0RSDNEbQijOSPMLM7g5/lH5IOM7Ovy2BaP+AjM/tTtHMjYBVgXrmcV3TyAI8ABxNG9XsDgyX9xcxmlNKetM9hNeAa4EkzuyZmebpb0v5FHF02FeZKuhp4izjtYWb/LFJftcCLkr6Kry1KMSfN85JGAG8SRtnFukOBMJr/paSpBMG6dYvYF+A7YxNH0upmNjft9aXARsCAOEVxIXCzmc2sq41S2SfpRuBZM3tcUlsz+74cdsTXnYEDzOwGSY8CmwGHmdm0Mto0lPDj+PP4t7sEeMjMXi+VTZWIpKMzy8zs/8phS1JIWgtoDWxIuIP7wcwalIbOs6+NgD8Q7mBnAH8zs4+K0dfSPt3RJ0OcFrmNMIofCfzLzF6M5y4FNgBeAk4Gdo7ytKWyLdOBtTKzxfH2dCrQEtgTOD9thF0sWzI/p/fM7NV47jzgcOAJYBIwshTTSg387YYSkn5MBAYD/czss2LbVOlIWglYm2WLsUVZTJS0PiFlYk1aX4kuksZ+xqQv8kp6wMwOSbqf2PYZZnZV2utBZlZUvRt39AkRRznbAecB+xBG8a+b2ZPx/B5AG4Jje79ENtXpVOP53YC/EDIYnWpm75XApmyf0xQzeyouFP/azB4tth052DTVzJ6I5/sR/nYfmNkHpbStEpF0NrAbsDlB9/6/Ztav/qvy7ut54M+E7+kFwF5mdnqC7e8PHEBYDB0Zi1sBHZLOQyBpDWAt4C7CgEaEQdaNZuaLsU2EKQRnsQR4COgJdJM0iZBB5sUyLHoeRYj6STmwPSR1NLMR8XwLwnzhb0o4F17X5/RatGUUlDzssz6bOgKvlnPBugLZ18x2kDTOzPpIur+IfS0xs9GSLjSzUfFHJknGEP7+/2aZiuQiQiRM0vQG9iMMJM4nOPpFwN+L0Ndy+Aae5PggPnoTfqVfBNoBWxGcxSplsGlKPKYc2MtAV0lrS+oa9c1rS7zgWdfntDUhLrstQIlDGCvxb1fJpKJgvpfUC9iyiH29oZDofIykscCChNvfBJhtZufGtk8g6PwnvkBqZo/HyJvt4vTTScBpZlb0qBt39AlhZj8AjxFW0fch3C19BfQCxpvZl2Uwqy4Htg2wjqS1zWxWKQ1q4HN6oRyfU4X+7SqZwXGO/gzgIODMYnVkZqeZ2Q9mdhFwILBvwl1cx7Ifj7sIa1YjCVOeiSLp73G97FNJxxHuJp6SNLihawvFHX2CmNnHwK0ER3EdIUzwjmIvcNZjT30ObKyZfVXP5cW062Mq6HOqVJsqmLeA/QkRSacDc+qvnj+SdpM0StI0wlrS0IS7WGhmS2LUzbpmdr+ZjSFE4CTNFmY2N66dnQX0iQvABxWhr+XwxdgiIakD0MLMSpVJvj5b1gP6EhacOgGnVMqiYiV9Tikq0aZKQtK9wPPAUXGuflQRF2MnADsTIrB2zoyOSaD9pwgblrYHPjWzv0hqBUwwsx5J9RP7egnYIfb1BzPbP5a/YGY7JdlXJr4YWyTMrGijnMYSwwHviF/qinJglfQ5pahEmyqMTmb2d0lFCT/M4Mf4sDgSVsLt/4YQtPAScGcsqwEuTbgfCHeKLxPWoQYCSNqUoIVVVHxE7zhOo5B0D0Hs7TTgCsJGtyOL1Fd/whrAFoRdq1eb2dPF6KuacUfvOE6jiFEwxxKc73uEnd4/FLE/EaKfZpdDpqMacEfvOE5OqIQyu5Lq1M5JcmespNXM7Nuk2qtUfI7ecZxc+Zig6vl2WpkIwmZJyxJ0ImjQD4+PYjnjx4FdACTdama/K1I/S5G0hwX5602B/wHutyLp+adwR+84Tq78irCz8+fAu8DDZjapGB2Z2T4K2b1Sc/RtgDuLPD//0yK2nc6ZhGQ2fwJuJoSM/rKYHbqjdxwnJ8zsNcKIHkm/Am6WNMnMjitSf/MlfQR8Tvhx6VSEbjaW9GfCnUnqear/C4rQH8BqcRrsRzN7SdJ3RepnKe7oHcfJCUkdgV8T9mT8AFwCjKj3ovz7upag+Poq8H/A+7E8aQ2kdMnlUQm2Wx+XAUOAIXFh++Vid+iLsU7FIOkIoA8h6ceRZla0zEJO45G0CHgdGEvQu1nqPJIe/UZdm3SMuB6Q5Iap2NfqhB2yCyTtCHQAnrGEM5xJ2jj9Jct/fkXVm3JH7ySOpLfMbKsCrj+TkNlpeoJmOQUiqXdd58zs+VLakhSSLiOkr2wJfEpIs/kfYBszOzjhvlKJyJV2hPDjlbjGfjo+deNUFJJWA35GvFV3Koem6swboKeZ7RRzIbyZGqBkuaNIgjOANpaWdEjSuhQxJ24KFzVrZkjqI2m4pEclTZW0r6SNJU2UNE4hb2aq7pWxbIKk7eK116edfyseW0q6J9a9JJa1kHSrpPGxfLM67HlD0t8kjYmKiFcC5/rGGKdELI4bsgBmKdCySH3dDKyaUdYOuKVI/S3FHX3zpIaQiLsPQdOjDyFUrg9BtRFJewItY9n+wF/rae/XwNex7si0su8tJHI+K9v1ktoC6wHXxnnX8wkhbkMkbVLIG3SSJ95tlbrPPeJxU0k3RP37pBkdHxCkg4u1KLuWmf0rvcDMPiSsCRQVn7ppnkwys8XAfyTNB+4HzpN0J0GC9q+EpBu7SxoXr1mNtMWjtFEQhKmWVOTAK/G4WarMzF5NOW5JJxNkWafGfh83s09jvfMSfp9OspR8cxFFjjk3s52zlUtaO6k+Mtpta2bfp73OHOEXBXf0zZPaeHu6anwsMbNzACQ9J+lR4B3gibTyNoQcoT+JbaT/s6WSm9xNSMmXKtsZuFvSdoRs95jZ9cD1sc3fU4LQMqcolGpzUVFjziWtSVCwhKBeuQFwLkFbZ7ck+yLkvR0h6Q7C3oD1gAGEnLhFxR198+QL4FGWfan3kXRSPPc5MMPM3pe0UxzRGyGe+WzCnOZYgh55iseBg2LdlON+DNhX0gvx+kFZ7NgG+EeSb8wpKuXYXFTsmPP7gDsIaSwnAJOAi8ws8WCAmPt2KmG379aE/7WDSyEb7uGVzQxJfYCDzOzkctviNC1KGV5ZqphzSaPNrG98PhXoVo2BAD6idxwnJ8zs+bo2FxWhu/OpI+acZAXUNkm7S1kDuDi1/FTEu5SS4yN6x3FyosSbizpQR8x5khnAqnETWDZ8RO84Tq6UcnPRzcAfM8raATcCByTVSTU58/pwR+84Tq6ssLmI4u3FyRpzHkf6TiNxR+84TmMYzbL58jHF7KhcMefViDt6x3FyosSbi8oWc16N+GKs4zg5Ud/mIjNLenMRktYixJx3Jjj7EaWIOa9G3NE7jpMTkp5l2eaiQYTNRX8txuYiJ1nc0TuOkxPNZXNRNeJz9I7j5Eqz2FxUjfiI3nGcnGgum4uqEXf0juM4VY4nHnEcx6ly3NE7juNUOe7oHcdxqhx39I7jOFXO/wefmuK03F08aQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x180 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = plot_results(loo, summary, chance_loo)\n",
    "plt.savefig('fig_model_fits_exp_scrambled.pdf', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Experiment Shifted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chance: 163444.105176\n"
     ]
    },
    {
     "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>LOO</th>\n",
       "      <th>dLOO</th>\n",
       "      <th>LOO_se</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Sticky-Choice</th>\n",
       "      <td>120859.072811</td>\n",
       "      <td>27308.598199</td>\n",
       "      <td>430.277644</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lin-Reg</th>\n",
       "      <td>115895.964440</td>\n",
       "      <td>22345.489829</td>\n",
       "      <td>444.243539</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bayesian-GP</th>\n",
       "      <td>109823.165807</td>\n",
       "      <td>16272.691195</td>\n",
       "      <td>457.807441</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kalman</th>\n",
       "      <td>109613.762119</td>\n",
       "      <td>16063.287507</td>\n",
       "      <td>459.350548</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Clustering</th>\n",
       "      <td>103176.607294</td>\n",
       "      <td>9626.132682</td>\n",
       "      <td>467.884898</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF</th>\n",
       "      <td>101631.696967</td>\n",
       "      <td>8081.222356</td>\n",
       "      <td>471.337316</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Linear/Clustering</th>\n",
       "      <td>100949.175983</td>\n",
       "      <td>7398.701372</td>\n",
       "      <td>468.220908</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF/Kalman</th>\n",
       "      <td>98830.526410</td>\n",
       "      <td>5280.051799</td>\n",
       "      <td>472.622342</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF/Clustering</th>\n",
       "      <td>93550.474611</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>477.185323</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                             LOO          dLOO      LOO_se\n",
       "Sticky-Choice      120859.072811  27308.598199  430.277644\n",
       "Lin-Reg            115895.964440  22345.489829  444.243539\n",
       "Bayesian-GP        109823.165807  16272.691195  457.807441\n",
       "Kalman             109613.762119  16063.287507  459.350548\n",
       "Clustering         103176.607294   9626.132682  467.884898\n",
       "GP-RBF             101631.696967   8081.222356  471.337316\n",
       "Linear/Clustering  100949.175983   7398.701372  468.220908\n",
       "GP-RBF/Kalman       98830.526410   5280.051799  472.622342\n",
       "GP-RBF/Clustering   93550.474611      0.000000  477.185323"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loo = pd.read_pickle('Data/model_fits/model_fits_exp_shifted.pkl')\n",
    "summary = pd.read_pickle('Data/model_fits/model_params_exp_shifted_rbf_cls.pkl')\n",
    "loo.rename(index={'Kalman w/SC':'Simple Kalman'}, inplace=True)\n",
    "\n",
    "#note the \"Simple Kalman\" is modeled without sticky choice and wasn't included in the paper\n",
    "loo.drop('Scrambled', inplace=True)\n",
    "\n",
    "# calculate chance\n",
    "n = len(pd.read_csv('Data/exp_shifted/gprbfshifted.csv')) # n_subj * n_trails_per_subj\n",
    "n_arms = 8\n",
    "chance_loo = -2 * n * np.log(1. / n_arms)\n",
    "\n",
    "print \"Chance: {}\".format(chance_loo)\n",
    "\n",
    "loo.sort_values('LOO', ascending=False, inplace=True)\n",
    "loo['dLOO'] = loo.LOO - np.min(loo.LOO)\n",
    "loo['LOO dLOO LOO_se'.split()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAADKCAYAAACrHYtRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXmYnEXVvu8nC7KGAAlIQERWURYxQVlD2JEoyCYqoPDJoqiAqAiCGBABQWWRRSN88AEu7IIQMUISMCxigCzsyCY/whJAE3YIeX5/VHWm0+mZ6cn0Pue+rr7e7nqr3jrd03O63qpTz5FtgiAIgvalX6MNCIIgCGpLOPogCII2Jxx9EARBmxOOPgiCoM0JRx8EQdDmhKMPgiBoc8LRB0EQtDnh6IMgCNqccPRBEARtTjj6OrDzzjsbiEfrPtqa+H629KMiwtHXgZdffrnRJgRBp8T3s/0JRx8EQdDmhKMPgiBoc8LRB0EQtDnh6IMgCIoYPXo0o0ePbrQZVWVAow0Igi75z3/g1lthu+1gueUabU3QBxg3blyjTag6MaIPmpsnn4S9907HIAgWiRjRB83NRhvBK6/AoEGNtiQIWpZw9EFzM2AALL98o60IgpYmpm6C5ubpp+HAA9OxDyLpTEl/l3R2SflxkmZKOrmobH1JkyXdIWnD+lsbNCvh6IPm5u234eGH07GPIemTwFK2twIWk7RJ0ekLgX1LmvwE+BLwhfw8CICYuqkLM2fO5MQTT2y0Gb3ixz/+cWM6/uhH4e67G9N349kMuCU/vwXYFPgngO0XJa1XUn95288CSFq2blYGTU+M6IOgeRkMzMnPZwPdxZf26+T5Qkg6RNIUSVNmzZrVCxODVqBujl7SIEk3SZok6S5JIyTtUXT+V520G1U8D1lhX0Mk/TH3NVnSZyq9jqRPSPpaT/oLasi0abDSSunY9/gvUAg3GpRfd8W8Tp4vhO2xtkfYHjF06NBemBi0AvUc0X8FuNb2KGAr4B1gvqO3/e0q9vUr4Lzc17Z0/w8yH9tTbV9URVuC3rDiivCtb6Vj3+MuYLv8fHuguzmsVyWtKmkY6Q4gCID6Ovo3gc0kDbE9l7RotEMedQ+VNBlA0hY5amCipH0KjfMdwQ2SPpbP9c/l10hasahef2Bl238HsP2u7bvy6Q0k/Tlff2lJAyT9QdLt+TigeOQv6RuS7s79rZvtvCG/Pr8un1oTMGfOHGbOnNmYzldeGX70o3TsY9i+D3hb0t9JI/R/SzoOIN91/gLYV9J5ucmPgT8CV+XnQQDUdzH2MmBVYKKkF4HjgdVs7wcgqVDvNGA32y9L6geMBJbJ7Y+1/ZCkW4FtJN0DDLT9UlE/Q4FOJx1tfy7/s2wHLAY8ZPtLko4H9gRezPasCOwNbGH7/WzLGcCptu+S9DNJmxX9iCyApEOAQwCWXba118Xuvfdexo4dy5gxY+rf+RtvwCOPpEXZpZaqf/8NxvYRJUU/zeUXAReV1J0ObFkn01qO0+7vme5+pfWP2XjIophTV+o2orf9nu2TbG9A+oIe2UXdl/OxMM+4NzDV9oP59e+BfYDdgevySHuSpEkkJ9/ZpOMD+fgcaaFrTeC+XDYFWKuo7keA+2y/X2TLesBpuZ/tgGFdvIf5c6BLLrlkZ9VaguHDh3PIIYc0pvNHHoERI9IxCIJFop6LsR+WNDC/fAl4DehfpqolrZDbFOy7GFhV0ucBbD9JcrJfIM37z7I9Kj/eB56XtFW+xkBJmxauXWwS8CQwPL8eATxRdP5JYOOCDfn4KHBU7mcEcP0ifRgtxqBBgxg2rNPftNry0Y/ClCnpGPRZ2lFRsp7Uc47+E8DkPBo+BjgZWF7S1ZKK97gfC/xZ0kTSSB6Sgz4U2E/S1rlsHPCu7XKLTt8Gvpn7mkgavZfjOuDjkm4HNgCuKZywPSu/vjPbsjZwCnC8pAmSbgE+1JMPIFgElloKhg/vk9M2QQfjxo1rS1XJelG3OXrb17PwCHinoudb5np3AJuX1JuUj3sVX5I0hVOur5eBL5Y5NSmfv6SobJ8ydQr1zgdKF113L9dnUCOefx4uvBAOOqhPLsgGQTVoyZ2xkg4jOdxdGm1LJQwbNqxxO0tbnZdegnPPhV13DUcfBItISzr6TkbaQTuy0Ubw4ouNtiIIWpqWdPStRjto3ZQj7lKCoDUIrZuguXnkEdh00wivDIJeEI4+aG4WXxzWWy8dWxxJO0j6qqT+ZZQng6BmtNXUjaRRwPa2j8+vf1WJhk4OwxRpp+zFtsfW0s6gB6y+Olx8caOt6DWSLgDeAEba/j9J5wA7NNishtOXd6vWk7Zy9KX0UChtO5KeyGQgHH2zMHcuzJmTcsYOaOmv67q2t817MqD8ZsGgCVh3y+0bbULVaen/nO6QNNn2lpIuIan5DQf+ZruzldEPAIvntiJF9qwLvAXsRxJmu4o08v8vcHNJTH6foiB0VtNds9OmJQmEKVPSxqnWZV7WT3LeIDi30QYF5TngnD802oSq05fm6CfZ3pLOY+9vBWbSEbb5WeDftrcFzgW+DnweuNP2zsB/uuqsOLHDm2++WZU30GyMHTuWsWNrfPOzxhpw1VXp2NocRVKWXB+4EvhuY80J+hJtPaIvoSBo9hZAVsDsT8cO2u2Aj5LkGS4kCZh9UdJOpM/pLtJu3Om5/tSuOsvz/GMBhg0b5q7qtip1ETpbbjnYa6/u6zU/y+dBAwCSNqukkaQzSTpM9xUrWUpaH/g1aW3pG7anS9ob+D7pe3pK3o0eBH1qRL+As7W9XRYne6Go7AGgX46IeBS4NNfZEvgh8BRJEwdgwzrZ3bQMGzas9mJnr7wCl1+ejq3NCSWvv9NdA3WdHLxcIvDvAKPy46he2hu0Ee04ot+3SK1yUbgQ+AZwBHCOpAm5/CzgT8BVkv4KvA681ytLg+55+mnYf/80R7/CCo22psdI+hLwZWBDSTfk4oFdNCmm0+TglE8E/ihQUH8r5JrtzK75+RJWW221Cs1pHO24QFpP2srR255E0pEvLT+g6PmoMudHFT2fABSc+0JRO5L2sD03h8s92Vubg27YeOOUfKR14+jHAXeSBg5n57J3gRc6bdHBYDqks2cDHy86Vy4R+LWk/Ar9gAO7unDx1OKIESOafmqxHRdI60lbOfo6cZOkpYF/dZZdqpQQNesF/fpBCyduyTLas4GjJK1EyqMgYGU6kt50RlfJwcslAv8JHT8G44Dxi2550E70pTn6qmB7J9tb2P5qo23pEzz5JHz5y+nYwkg6iZTf4CbgPFLKzO7oKjl4uUTg75BCgN8ghQDXlEgG0jrEiL4OtKuoWVdU7Q7mvfdg5sx0bG12sv3pvAt7GzrJpVCM7fskFZKDTyMnB7f9UzoSgQv4Zm5yAXBHfl7zTX+RCKR1CEcfNDfrrguTJjXaimpQWBx9n5Ts/mOVNOoiOfhCicDz5r1LemNk0J7E1E0Q1IeLJC1Bit66jQpG9EFQLRrm6CWNkvSMpEmS7qi2mp+kX/Wi7RKSzpM0UdJkSefl8jGSpmV7T6metUGn3H8/LLtsOrYwtv9o+y3bf7a9se2fNdqmoO/Q6BH9ZTm08fskiYGq0UNBs1JOAO6yvU3eLHVF0bnv2t4C2EhSJAfvhJkzZ87XwukVK68Mxx3X8mkEJX0+D2qmSZouaXr3rYKgOjTLHP0gYI6kY4GdScJiXyfFGp9p+4uSBgDjswLgQUAh6uUI4GlSDLGBGbYPLxI0OzDXXRr4oe3xFYicbWH72MIL27eXsXk6sArwbLk3VLwhZdllly1Xpa0paOCMGTOmdxf64Afh6KN7b1DjOQnYxfb/a7QhldAT+eCQDm5+Gu3o95c0Elgb2BF4wvapktYCTrS9r6SlJC0DbA7cImkIsCswElgO+F+S6Ngk22Oy6mQxV9i+OO8evIqO2OJJto+Q9A+g1NHPv4ak8aTY552LyvoDnwLO7OyN9QWtm66omg7O668nBcuNNoKll67ONRvD48DzjTYi6Js02tFfZvv4vJHkQuBGSfuSNoAUnOO1wG7AtsDJwBrARsDEouvcBoyU9HvgL8BlRed2knQEyXmvWFTencgZALZ3zHcAhc/qFyTlyj/afmkR33fbUzUNnEcfhS23bAeZ4kuBxyU9mF/b9q6NNCjoOzTa0Rd4jTR9cxiwMbAm8Nt87hrgd8BA209Keg34p+29ACQNBPrbPiG/nsqCjv5YYGuS1vwdReULiZwVnku6U9L+tgvXKf6cvmv7FoL68LGPwYwZsOaajbakt5wCfIVOpvqCoJY02tHvL2lL0pz8yaSR++35AYDtOZLeJos72Z4l6SZJt5NikicAt+UomIF0iEAVuDFf7x4W3ELeFScCv8xrAe+Sbrsr0SYJqs0SS8D66zfaimrwqO3JjTYi6Js0zNFnAbIPlxTf2En190jJGgptLwZKE4mWbh7ZMh9PIi2EFZ87oOj5qDK2vUn5KKAxndjXJaF10wueew7OPx8OOwxWWaXR1vSGpSTdQdK3MYDtwxtrUtBXaPSIvlskjQVesh0LWX2RV19NevRf/GKrO/q223cR0sGtQ9M7ett1SGMUNC0bbADPPNNoKxYZSWvbfpy0DtVWhHRw69D0jr4d6IuiZlBFYbPWZnfgdBbObWDgf+pvTtAXCUcfNDcPPQT77gu/+12KwGkxbJ+en15qe35IsKTNG2RS0AdptARCVcn6OSfn53tKukZSv5I6q0u6vDEWBj1mqaVg003TsbX5UcnrIytpJOlMSX+XdHZJ+fpZh+kOSRvmsuUlXSlpgqTjqmR30Aa05Yhe0hYkje7Rtud1Vz9oYj78YbjggkZbscj0JmdscXJwSRdI2sR2IWdsITn4POB8Umjyj4ETbD9S1TcRtDzt6OjXJWXl+SwwN+94HQjMAr5QXFHSXcAUkpzCqfn8WsD+tqflUdQnSHc++9r+d25zLymc8xjbN9fnbbUWc+bMYebMmb3fIfvuu/DyyzBkCCxW86RJtaA3OWN7mhx8feCHWWzvh12lumy15OBB72irqZvMjsBfbb8CzAU+a3sk8DBJRqGYFUgbtT5D0gn/Mil+/mv5/LG2tyZtoDq0qM1JwOiisoWQdIikKZKmvPnmm1V5Y63EvffeO1/YrFfMmJHCKmfM6P21GoDt2bafAa7Kx/dId5sbVtB8MB0JS2aTtJ0KlEsOvjlpwPJF4Ixu7Bpre4TtEUOHDq3AlKCVaUdHfwGwuaSdgaVICR9uA/YiiZMVM8v2i7ZnkgTV3gZm0vEPdXRO43ZyUdtZtl+y/RzpH7Esxf9IS7ZwcutFZfjw4dURNltzTbjhhnaQQDg1H08AHiIJ8XVHT5ODP2b7YdsvlpwP+jjt6OjnkkY0J5Nulx/Lo/JrKFKlzLiT55K0AjDK9lakhTSVq1dNw9uJQYMGVUfYbPBg+Nzn0rG1KXxXlrT9e9LIvjt6mhz8MUkrS1qK9pyWDRaRdnT02H4V2B84EPiKpBuB1Xt4mf8Ar0uaQJqmCRrByy/DhRemY2vzuKR/AtdnIb5uBwm27wMKycHnkZOD59OF5OBX5eeFsj+Q9J9OrrL9QQvTVr/6WT9nUn7+MGlhtRz75Trz9XGKtHGeLpwnLeiW9lHcZlSvjQ665pln4OCDYeON04Jsi2L7IEkDbM/NIb+fr7BdT5KDPwSMqoK5QZvRVo6+WQlRs17wyU/C3LnQrzVvPiWdXeSsvw+cantejug6oHGWBX2J1vzvCfoOEvTvn46tSXF0zQ5Fz0uVW4OgZoSjD5qbJ56APfZIx9akszSSfS69ZNA4YuqmDvRVUbOuqHgqa948ePPNdGxN1pF0Dmnxtfj52o01K+hLhKMPmpu114abW3rz8b5Fz6/u5HkQ1JSWcPSStiNtNBHwCsnu79j+VwVtVwfWsD2hwr7OIuWFfX+RDQ6CjO3bGm1DEDT9HL2kISQnX5Ay+AHQE9GT1VlY+qCzvvrZPjKcfBNx332w+OLpGATBItH0jp60Weky268B2H4MeB5A0hhJ2+fnl2QJ4t0l3ZOlWnchCTftn8XNkHSCpEn5/Or5MVHS1cAB+dyAfO0LJd0i6cLcdk1J/5B0vaRb891C0EMKgmcVseqqcPrp6djCSFqx0TYEfZdWcPQrkx17hewJfMH2tsBfgLGkH4rtJG0ArJI3On0TODa3WRHYx/b/llzrQdvbA6tJGgx8jySrsEdu0yl9XdSsK3okeLbiinD44enY2vyu0QYEfZdWcPTPs7AYWYFyujMnA8dLuoSFd8auB4ySNIkkflYQjJrWyXTNA/k4E1gW+AgwPdd9oEz9DsP6uKhZV/RI8GzOHLj11nRsbaZK+pykoTlByPKNNijoO7TCYuw44GpJf7T9mqS1SKN8SGJOK0sS8PFc9kzebr45cBRwOdA/n3sUGG/72wBZc2QVOlf6K/0heQrYQNKUov6CHtIjwbPHH4ftt4cpU2D48NoaVltG5EcBU+HaURD0lqZ39LZnSfoJcGN26K+SEjcAXAtcB+xKEiEDGCNpU2Bp4Lukkfepkq6wvY+kF/KI3iQBqPE9MOfnpFvwl3J/lSgQBr3h4x9Pzr7F5+htbwMgaWnbrzfanqBv0fSOHsD2LXRk2imldJj3gzJ1RhZd66dkYagi9is6Pyo/HVNUdgBAFqXaVFJ/4A7gxe6tD3rF4ovDWp1p07UOOTDgRGBJSRsBY23/T4PNCvoILeHom4g1cgTOUsBFtudW0ihEzXrBs8/CWWfBkUfChz7UaGt6w3Ektcmbs4Ll6pU0knQmacrnvmIlS0nrA78mTSl+I6tZImkJ0hTjfnmAFAQtsRjbNNh+zPZI28Ntn99oe/oEs2fDn/+cjq3N+7bfoWPdp39XlWHB5ODAYpI2KTpdSA7+hfy8wCF0EygQ9D1iRB80N+uvD4891mgrqsG1kq4ihepeRsp41h09Sg4uaTHg08Dk7i4cycH7FuHo60CImnVOX5nSsn2WpL+SorUetv1gBc0GAwXZztksGOlVLjn4gcBlJGffnT1jSXtMGDFiRChptjkxdRM0Nw8+COutl44tjKSbc+Luq20/KOnyCppVnBxc0gBgJ9t/qZLJQRvRdo5e0tZZ3mBSlinYXNLsLHMwWdI6JfUPkPRoPv+nfPtLLpuUJQ8OL6k7SVL8Q9WDQYNg553TsQWR9FFJewKrSNojP/ahssQjPUkOvhLwIUk3k6LITpW0XPXeSdDKtNXUTRZAOxHY1fYcScuQdsfOsL2NpE8DRwMHlTQ9w/aFko4HdgRuBGbZHpVj9+8GzimuW5c3FKRImzPPbLQVvWElYH3Svo4Nctm7pD0eXWL7PkmF5ODTyMnBc4hwITm4gG/afg7YBJIGFDDZ9n86uXTQx2grRw/sAlxuew5AFkK7Xx1p6KYCXe28GVymbDFiY1RNKIibdblL9p13YOZMGDYMPvCB+hlXJbJM8W05THIeSWvp0R60rzg5eFGbMYtmbdCutNvUzTA6lC2/nKdqfl50fiRJBqGU70t6BPgk8NdcNjTvoH0IuLWk7iRJp3ZlSIiadU9F4mYPPABrrJGOrc1eJDmPKyX1l3Rjow0K+g7t5ujnC6DZ/j1prnIISZ9mInAYcJqkgrM+ILc7g5TE+XWgIDY1K++SXQtYT1JhTvUM26NsF5QvyxKiZt1TkbjZ2mvD+PHp2Nr8D7A18GoWxYsvRVA32s3RjwP2K8QV0zE1NcP2NrZ3t/287YKzvqTQ0Pa7wHnAd4ovaNvAa5Sf1gl6QUXiZoMGwQ47tOxibBGFKJmKN0wFQbVoK0dvexZpMfZ6SROA84FLe3CJvwHbSfoAeepG0u3AXNvTqm9x0C0vvQTnnpuOrc15wERSgvDx+XUQ1IV2W4wl54YtzQ/bab7YklH9PHLkArBuV3WDOvHss/Cd78Bmm7V08hHbV2QHvybwpO1XG21T0HdoO0ffjISoWS8YPhzea/2gp5yhbFfSGtBISdj+ZYPNCvoIbTV1EwRNzF+AD5HyKbySH0FQF2JEXwdC66Zzur3TefxxOOIIOPvsVo+8ec32yY02IuibhKMPmhsJBg5Mx9bm0pwpbSo58sb2tY01KegrhKMPmpu11oLrr2+0FdXgq8AzQCFZjUmpMIOg5vR4jl7SKEknF73+VXVNAkmrSDorPz9K0u15l+vZuaxbve0Se9foQf2qv58gICUeOcj2iflxUqMNCvoOvV6Mtf3tahgCIKlgzw7A3yR9hhTmuLXtLVm0EdAooCJHL6lfNd9P0D0zZ85k5syZnVe49940bXPvvfUzqjY8K+nrkoZL+mTOHhUEdaHXjr4wupZ0iaSz88j7x7lsLUnjJd2WlSGRdGDeiDRF0o5Fbc8Fbs6XHQXcRkqTdkbenVoQiCrue1LW4Sbr0iDpFEl3ZNnh1YADgF9I+oWkJST9IcsYXyFpYJYevkLSTcCGRe9nkqRTJf1T0tdy2a6S7pU0tru7itC6qYyxY8d2rXez2mrwm9+kY2szgJQQ5FvAt/OxSySdKenvhTvZovL18//ZHZI2zGW/ya8nF8qCoEC15+gn2T5C0j9IO1R/CnzN9rPZwa4KXGH74ixTcBUwPre9w/a3sizwkrZfl7QyWaSsB2wBbGV7Xr7WJSTJ1lskfRu4wfYfJH2DJDQF8F/b+wBowUW/K4ExpB2zFwE/IAmjDQYmdWVEcQafYcOGRQafTuhW62boUOiuTgtg+8Ce1FdRvlhJF0jaxHYhjWAhX+w80u7v3YDTbD8laW3gNGDPKpoftDjVdvQFicG38nFd4LLsPAcDqwCbSDqCpKNdvNWxcG++EUl7GzpEyh7vpL9iB1rw0KcD/yfpFeC4kvrrAcMlHQosDvyBlLShs3mBB2y/J6mgU/K+7TeANyS93EmboAd0q3UzezbccQdssQUsu2zXdZsYSZuSNOSHkb+rtrsaefcoX6ztp/K594D3q2p80PJUe8NU6cj1UeBLWQVyOOmLeizwGdIoZIF0aPm4I2kEDWlE/b08MkfSyJLrzwZWVko48sFcNsH2/sBLwGdJX/yCgNSjwOlZ0GxT0miouO/u3k8/SUsqZfUZ0kmboJr8618wenQ6tjbnAEeQNkztA3QXSjQYmJOfzwaKs0WVyxdb4FQ6kuR0SvHU4qxZs7qrHrQ4izqi3zePULrjOOB/s0jYe6TbyRuB24F7WDAHZoHhwM8BbP9F0nqkxA39SCPv24vqjgX+TMp6X/i2/klSQQJ2b+D/AacoZZc6A/itpMNIo6oupYbLcHrufyrwYg/bBovCBhskvZsW1rnJzLH9mCTZfljSZt3UrzhfbOGJpCOBh2x3G5UWycH7Fj129LYnAR8pU35A0fNR+fgv0ui9mJPyo2xb4JosLlY490vglyX1t8zHcSRp4uJzO5X09wJJB7zAfiXn/9HJtUeVvh9gnO0bJK1C/icJasxii8GqXSUFaxmelbQEcJ+ka+j+f+8u4FDSXe32pLWmAq/m9a55pNE+ObBhc9LdQhAsQNNtmLJ9ZaNt6IK98iLuUsDhlTYKUbNe8O9/w+mnw9FHt2zkTZ56PMf2W5K+R1qH6jKdYE/yxeYmvyJN9UyU9KjtQ2v1foLqMXr0aABuuummmvbTdI6+mbH9R9I/WFAvXn8dbr8dDjus0ZYsMrattMlwdL5bvb/CdhXni7W9kKx20PyMGzeu+0pVIBx9HQhRs8ooe9fzsY/B9On1N6b6PJ3De+8mR8XYvq+xJgV9hXD0QVAfliQlny/siDUpj2wQ1Jxw9EFzM2MG7Lor3HBDisBpUXq6YSoIqklV4+glbZ3lBSZJulXS5pJmZzmCyZLWKal/gKRH8/k/SVoslz+ar/EPSYeX1J0k6S9F17hC0hqSLs+vN8/XW5IydCddEDQZyy0H++yTji2MpE0l/UXSNEnTJbXFfFTQGlTN0edNSycCn8/hiJ8n7ZCdYXsb4LvA0WWanpHPTyFtlgKYla+xKbBvSd1Rtj+T+1wceIccS5x/SH4BfMF2CMy0EHPmzCkvbrbqqnDaae0QYtnTDVNBUDWqOaLfBbjc9hwA26/ZLo4umAp09d86uEzZYqSNVp0xkrRZCmAl4FJgP9uzACRdqSSoNl7SoOKG+c7gp3l34MGSLs+jrZ3y+WNz239I2riozQJCZ52hEDXrEffee295cbO33oKHH07H1maO7ccA2X6YJHEQBHWhmo5+GFmATNKX81TNz4vOj6R87PD3JT1CWqT6ay4bqqRG+RBwa0ndSZJOza93oEMU7VPAI7afKKp/gO2tSZtOym0kuZokgnYa6Y7jM3TEJZ+d2+4LfK+ozZWk0LavlrnefGyPtT3C9ogllyw7ixQUMXz48PICZw89lCJvHnqo/kZVl55umAqCqlHNL1tBgAzbv5d0J0n5cQNJE0lbuA+T9H1gNB07/c4gjcSvBpYnSQvMsj0qbzS5QtKHC3VtX1jU50dsPy1pdZIUQj9JB9v+raT+wBmSNiBtIb+ujM0F0bJHbL8IIKkwGby/pH1J00Iu06YzfZxgERg0aFB5gbN11klx9Ouss/C5FqKwGFvphqkgqCbVHNGPA/YrqOnR8SMyw/Y2tne3/bztwjz7JYWGtt8FzgO+U3zBrEP/GmWmdSStSJI3KObgbMNI4BMkmdeR+drlko665EhRvcNIuvgHl7QNXZB6sswysNVW6diCSNpB0gxJT0i6B9jQ9v2xhhTUk6o5+jwvfiJwvaQJJGXIS3twib8B22UBtKF5iuZ2YK7taWXq70CHjGvBhndIyUrOAt4F1pJ0M2lap6fcQxIwi7C4RvLii/CLX6Rja3IysLPtNYGvkaYJg6CuVHWe0PYEYEJJcenr4vqXFD2fB2ySXy60nbu4buZlkvATtp8mi5XlKZjCppQtS9p0Jlq2ZZnzB5dpO6rc8+4IrZteMHMmjBkD224LK63UaGsWhTdsPwdge0YhhDgI6knLLgjZ/mv3tYKWZ+ON4bXXGm1Fb9hQ0g35uYpf2961cWYFfYmWdfRB0CIMb7QBQRCOvg6EqNmis/zLL/Pthx6CCy5oycgb28/0pr2kM4ERwH3FapaS1gd+TbpL+Ibt6eXKetN30D6Eow+amnn9+sFhYAoFAAAQKElEQVSQITCg731V1fME4eXKggZw2v09SyldSf1jNl707KUVRd00UMOmn6SNJd2Sd6neLGlNSWMkbV+h7YMl7VHpByLpGKUMUkET8N/ll4crroA11mi0KY2gXILwAsvbfjYv9C7bRVkQdD+iL9Kw2dX2HEnLAGuR4+OVcrEeDRxU0vQM2xdKOp6kYXMjC26EupuOJMYLbIQq0rAZAJwL7Gn7BUkfpOdf4MHAHsC1FbzXfrYj/K2ZmDcP3n47pRTsV+1c9k3PYKCw03s28PGic+UShHeVNHwBJB0CHAIwcOBA0r9kbelpgua+1Nexnxy6SP2krUbdU8n98EIaNsD9RV+MWmrYbAZMtP1C7vsF4IVC35IOAAbkH5QxwCRS/PyZwBvA5aQfpR2ypMLepNvZgnzBETll2zRgBvCApI+SYp+3BHai44dlNGnn7tUksbb3gDNzDt2gRgx64glYYgmYMgWG97l1zZ4mCC+bNLwcpcnBp0yZ0jtLg0Wi4MsqddiLSiVDpEZq2Kxc6LsH7AL8wPa2wMWkL/Pfcty7gV2zzbsBJ+Q2qwKHlhnNv2R7F+A5YEPSXctvclmX8dAhalYdJj71FNftthusvnqjTWkEdwHb5efbk+6CC7wqaVVJw8gJwjspC4KKHP0CGjakjUlD6NCwOQw4TVLBWR+Q251Bco6vk0bC0CE/vBawXomGzSjbhbuTj+RNUPP77oRy0gXnA1+QdBkdG7AKrEHSGZlImsop3G08avuNMtd/IB+fy3U/AhQiGaZ2YVeImlWJtTfdlE+ffz6ssEKjTak7OdVgIUH4PHKC8Hy6kCD8qvy8s7IgqGjqZhxwtaQrbc+mRMOmqN4Z+VGYUsH2u5IKGjbHFCrmZMmVaNjcBZwi6YN5jn5FoDgDxWygkHZoA5ID/4/tw/Ko5iLSKLx/rvMU8E/be+W+Bubyzm5zS39Insr9PET6EYtNWzVmxYEDGXb33bDNNi2ffGRR6GGC8IXKggAqGNE3UsMmi519C7hc0m3AZSw4t38rsHPRzkOAQ/P1byQpZL4ALC/palJS5psk3Z7vRo6hZ1wIfCPr54iu1xmCKrDcf/4De+4JTz7ZaFOCoGVRrRcBeopS4o+7Cou/zYSkfpB0eSTdBBxS0DHpimHDhvnQQw+tuX3tiN5/nxMOPxwGD25kLH3tQ1IaSCzGNo4qLMZW9N1sul0oTa5hszTpjmAx4JZKnDyEqFkQBI2lzwUm9wbbc2xvZfvTto/rvkXQa55+Gg46KB2DIFgkwtEHzc1bb8H06e2QMzYIGkbTTd0EwQKstx7cc0+jrQiCliZG9EEQBG1OOPqguZk+HYYNS8cgCBaJmLoJmpshQ+CQQ9IxCNqMXXbZpS79hKMPmpthw1LO2CBoQ2666aa69BNTN0EQBG1OOPogCII2Jxx9EARBmxOOPgiCoM1pOlGzdiRLMpdLztJohgA9y2Jce5rRpsVtr99oI2qFpFnAM4vYvJ5/r+hrYV62vXN3lSLqpj48antEo40oRdKUZrOrWW1qtA21xHb3CUs7oZ5/r+hr0YmpmyAIgjYnHH0QBEGbE46+PoxttAGd0Ix2hU2tRT0/m+hrEYnF2CAIgjYnRvRBEARtTjj6IAiCNiccfRAEQZsTjj4IgqDNCUdfRyTF592CSFq60TYEQW+InbE1JDv2w4D/AnfYfkpSP9vzGmzafJrFHkkCLgAuAabYnttYi+b//c4DrgduzmVyhKohaQNgMCAA27fXqJ8dga8Dy5IGpra9bRWvfzFQ9u9p+3+q1U/ua7XOztn+dzX7KiUcfW25BHgD+AcwVtKhtp9spHPNzuubwCvAnbafbgbnZduSPgW8nszUfbbfaaRNwM+Bx23fLGklYJ7tWc3y49goJN0IPAfMzEUGauLogdOA0bafr9H1x+TjUcDdwL3AxsAWNejrxHxcBVgZmAFsALwA7FCD/uYTjr62PAKcY/t1Se8Dv5b09ezs6+5c86j5dyQnPw34paQDbM+ppx1l7Opv+33gGmAN4JPAW5IetP1eA037J/CApGtJol+bSdrX9hPN8OPYQGT70Dr1dQfwbq0ubvsZAEmfsH1ELn5M0mE16OvA3Ne1wCdsvy+pP+l7X1PC0deAIicwgDSS38/2ZcnP8htJX7T9SgNM2x54yvYPs52rA0sAcxrpuLKTB7gW2Js0qv8scISkn9h+sp72FH0WywBnAX+2fZakrwC/k7R7DUeYrcBsSWcCD5CnPWz/b436Gg7cKeml/Nq2R9agn9skjQOmk0bZtbpDgTSa30TSVOATwAdr2BcQO2OrjqRlbc8uen0KsDqwb56e+DHwW9szO7tGveyTdD7wV9vXS1rS9puNsiW/XhnYw/Z5kq4D1gG+aHtGA206h/QD+fH89/spcLXt++tlU7Mh6aulZbb/rxG2VAtJKwADgQ+T7t7etv3fGvW1OvB90t3rk8DPbT9Vi77m9xmOvjrkaZGLSaP48cC/bN+Zz50CrAbcBXwL2Mb2C3W0rdR5DbA9N9+eTgX6A58BflQ0uq6lPaWf1SO278nnjge+DNwATAHG12NqqZu/3znAcqRphCOA7W0/V2ubmhlJHwBWpGMxtiaLiZJWBQ4FhhX1VdVF0tzPhOJFXklX2v5CtfvJ1/6u7V8UvT7Edk31bsLRV4k8yvkUcDzwOdIo/n7bf87ndwYWIzm1x+pkU6cONZ/fEfgJMBs43PYjdbKr3Gd1n+0b82Lxbravq4ct3dg01fYN+fz2pL/f47Yfr6dtzYakY4AdgY8CTwDv2N6+Rn3dBpxE+p6eAOxi+6gqXn93YA/SYuj4XDwAWL6ShB497GswsAJwOWkwI9Ig63zbsRjbItxHchTzgKuBLYGNc9KKIaQIl3oven6FFPVTcF47Sxpie1w+3480X/ilOs+Dd/ZZ/TPbcwvUPfSzK5uGAPc0etG6idjV9uaSJtkeJemKGvY1z/atkn5s+5b8I1NNJpD+9v+mQ0XyPVIkTLXZGvg8aRDxI5Kjfw/4dQ36WoDYwFM9Hs+PrUm/0ncCSwHrkxzFEg2w6b58LDivu4ENJK0oaQPbNwPD673YSeef1Yak2OwlAeocwtiMf79mpRAF86akkcDHatjXNEmLAxMkTQTeqvL11ySl4zsuX/vrwOHUYIHU9vU58uZTefrpm8CRtmsedROOvkrYfhv4E2kV/XOku6WXgJHA7bZfbIBZnTmvjYCVJK1oe1a9jerms/p7Iz6rJv37NStH5Dn67wJ7Ad+rVUe2j7T9tu0xwJ7ArlXu4ld0/HhcTlqzGk+a8qwqkn6d18uelXQw6W7iRklHdNe2t4SjryK2nwYuIjmJX5FCBC+txwJnJ/Z05bwm2n6pi+a1tu1pmuizalabmpQHgN1J0UhHAa/WqiNJO0q6RdIM0lrSOVXu4l3b83LUzQdtX2F7AikCp9qsZ3t2Xjs7GhiVF4D3qkFfCxCLsTVC0vJAP9v1yiTflS2rANuRFpyGAt9upgXFZvqsCjSjTc2CpD8AtwFfyXP1t9RwMXYysA0p+mqb0uiYKlz/RtKGpc2AZ23/RNIAYLLtTavVT+7rLmDz3Nf3be+ey/9ue6tq9lVKLMbWCNs1G+X0lBwKeGn+Ujed82qmz6pAM9rURAy1/WtJNQk/LOH9/HAeCavK1/8SKWjhLuCyXDYMOKXK/UC6S7ybtAZ1AICktUhaWDUlRvRBEPQISb8nCb0dCZxO2uS2f436Gk1aA1iPtGv1TNt/qUVf7Uw4+iAIekSOgjmI5HwfIe30fruG/YkU+fRyo2Q6Wp1w9EEQVITqKLMrqVPtnGrujJW0jO3XqnW9ZiXm6IMgqJSnSYqeDxaViSRsVm1ZgqEkDfqb8qNWzvh6YFsASRfZ/lqN+pmPpJ2dpK/XAr4DXOEa6fkXCEcfBEGlfJq0s/PjwMPANban1KIj259TyuxVmKNfDLisxvPzH6nhtYv5HimRzQ+B35JCRjepZYfh6IMgqAjb/ySN6JH0aeC3kqbYPrhG/b0u6SngedKPy9AadLOGpJNIdyaF54X+T6hBfwDL5Gmw923fJemNGvUzn3D0QRBUhKQhwG6kPRlvAz8FxnXZaNH7Opuk+HoP8H/AY7m82vpHxZLLt1Txul1xKnAycHJe2L671h3GYmzQNEjaDxhFSvixv+2aZRYKeo6k94D7gYkkvZv5zqPao9+sa1OMyesB1dwwlftalrRD9i1JWwDLAze7ytnNJK1R/JIFP7+a6k2Fow+qjqQHbK/fi/bfI2V1erSKZgW9RNLWnZ2zfVs9bakWkk4lpa7sDzxLSrP5X2Aj23tXua9CInIVHSH9eFVdY7+YmLoJmgpJywBrk2/Vg+ahVZ15N2xpe6ucB2F6YYBS5o6iGnwXWMxFSYckfZAa5sQtEKJmfQxJoyTdJOk6SVMl7SppDUl3SJqklDezUPeMXDZZ0qdy23OLzj+Qj/0l/T7X/Wku6yfpIkm35/J1OrFnmqSfS5qQFRHPAI6LjTFBnZibN2QBzFKif436+i2wdEnZUsCFNepvPuHo+ybDSEm4R5E0PUaRQuVGkRQbkfQZoH8u2x34WRfX2w14JdcdX1T2plMi56PLtZe0JLAKcHaed/0RKcTtZElr9uYNBtUn323Vu8+d83EtSedl/ftqc2t+QJIOrtWi7Aq2/1VcYPsJ0ppATYmpm77JFNtzgf9Keh24Ajhe0mUkCdqfkRJu7CRpUm6zDEWLR0WjIEhTLYXIgX/k4zqFMtv3FBy3pG+RZFmn5n6vt/1srnd8ld9nUF3qvrmIGsec296mXLmkFavVR8l1l7T9ZtHr0hF+TQhH3zcZnm9Pl86PebaPBZD0N0nXAQ8BNxSVL0bKEfqhfI3if7ZCcpPfkdLxFcq2AX4n6VOkbPfYPhc4N1/zG9QhtCyoCfXaXFTTmHNJy5EULCGpV64GHEfS1tmxmn2R8t6Ok3QpaW/AKsC+pJy4NSUcfd/kBeA6Or7Un5P0zXzueeBJ249J2iqP6E2KZz6GNKc5kaRHXuB6YK9ct+C4/wTsKunvuf0hZezYCPhNNd9YUFMasbmo1jHnfwQuJaWwnAxMAcbYrnowQM59O5W023dD0v/a3vWQDY/wyj6GpFHAXra/1WhbgtainuGV9Yo5l3Sr7e3y86nAxu0YCBAj+iAIKsL2bZ1tLqpBdz+ik5hzqiugtmbRXcpg4MTC8lMN71LqTozogyCoiDpvLlqeTmLOq5n9qx03gZUjRvRBEFRKPTcX/Rb4QUnZUsD5wB7V6qSdnHlXhKMPgqBSFtpcRO324pSNOc8j/aCHhKMPgqAn3ErHfPmEWnbUqJjzdiQcfRAEFVHnzUUNizlvR2IxNgiCiuhqc5Htam8uQtIKpJjzlUnOflw9Ys7bkXD0QRBUhKS/0rG56BDS5qKf1WJzUVBdwtEHQVARfWVzUTsSc/RBEFRKn9hc1I7EiD4IgoroK5uL2pFw9EEQBG1OJB4JgiBoc8LRB0EQtDnh6IMgCNqccPRBEARtzv8HNB0Aeb2QyNEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x180 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = plot_results(loo, summary, chance_loo)\n",
    "plt.savefig('fig_model_fits_exp_shifted.pdf', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Experiment SRS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chance: 140986.136526\n"
     ]
    },
    {
     "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>LOO</th>\n",
       "      <th>dLOO</th>\n",
       "      <th>LOO_se</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Sticky-Choice</th>\n",
       "      <td>88261.808536</td>\n",
       "      <td>24670.740192</td>\n",
       "      <td>452.377869</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lin-Reg</th>\n",
       "      <td>86065.490931</td>\n",
       "      <td>22474.422586</td>\n",
       "      <td>451.208158</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bayesian-GP</th>\n",
       "      <td>69447.723773</td>\n",
       "      <td>5856.655429</td>\n",
       "      <td>450.577691</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF</th>\n",
       "      <td>68376.396843</td>\n",
       "      <td>4785.328499</td>\n",
       "      <td>449.838894</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Clustering</th>\n",
       "      <td>66262.397727</td>\n",
       "      <td>2671.329383</td>\n",
       "      <td>444.585014</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Linear/Clustering</th>\n",
       "      <td>65932.360612</td>\n",
       "      <td>2341.292268</td>\n",
       "      <td>443.335323</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Kalman</th>\n",
       "      <td>65121.517049</td>\n",
       "      <td>1530.448705</td>\n",
       "      <td>441.186982</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF/Kalman</th>\n",
       "      <td>64290.905294</td>\n",
       "      <td>699.836950</td>\n",
       "      <td>440.729319</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GP-RBF/Clustering</th>\n",
       "      <td>63591.068344</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>443.057297</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            LOO          dLOO      LOO_se\n",
       "Sticky-Choice      88261.808536  24670.740192  452.377869\n",
       "Lin-Reg            86065.490931  22474.422586  451.208158\n",
       "Bayesian-GP        69447.723773   5856.655429  450.577691\n",
       "GP-RBF             68376.396843   4785.328499  449.838894\n",
       "Clustering         66262.397727   2671.329383  444.585014\n",
       "Linear/Clustering  65932.360612   2341.292268  443.335323\n",
       "Kalman             65121.517049   1530.448705  441.186982\n",
       "GP-RBF/Kalman      64290.905294    699.836950  440.729319\n",
       "GP-RBF/Clustering  63591.068344      0.000000  443.057297"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loo = pd.read_pickle('Data/model_fits/model_fits_exp_srs.pkl')\n",
    "summary = pd.read_pickle('Data/model_fits/model_params_exp_srs_rbf_cls.pkl')\n",
    "\n",
    "#note the \"Simple Kalman\" is modeled without sticky choice and wasn't included in the paper\n",
    "loo.drop('Scrambled', inplace=True)\n",
    "\n",
    "\n",
    "# calculate chance\n",
    "n = len(pd.read_csv('Data/exp_srs/gprbfsrs.csv')) # n_subj * n_trails_per_subj\n",
    "n_arms = 8\n",
    "chance_loo = -2 * n * np.log(1. / n_arms)\n",
    "\n",
    "print \"Chance: {}\".format(chance_loo)\n",
    "\n",
    "loo.sort_values('LOO', ascending=False, inplace=True)\n",
    "loo['dLOO'] = loo.LOO - np.min(loo.LOO)\n",
    "loo['LOO dLOO LOO_se'.split()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAADKCAYAAACrHYtRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXmYXEXZvu8nIcgaQiAoAREBRZTFkCggEAIERCIgCKICip8ICgooiiBIQNkE+RAU0ADKx6KCLLJFxJCEHSEJkCCr7D/DEkAStrDl+f1R1ZlDZ2YymTm9TM97X9dcdbqqTlV1z8zbdareel7ZJgiCIGhd+jV6AEEQBEFtCUMfBEHQ4oShD4IgaHHC0AdBELQ4YeiDIAhanDD0QRAELU4Y+iAIghYnDH0QNCmSTpV0s6TTqvJ/J+lWSbdIWj/nDZU0UdJtkkY3ZsRBsxKGPgiaEEkbAkvb3hxYXNKnCsUn2t4U+AYwNucdBhwJbJvTIJhPGPo6sN122xmIn9770wg2ASbk6wnAxpUC24/ny7eBd/P1+sDttl8FXpG0bFc7ir/PXv3TJcLQ14EXXnih0UMIeh+DgDn5ejawfDt1TgBOz9f93aZn0lH9+UjaV9IUSVMeeOCBMsYbNDFh6IOgOXkZGJivB+bX85F0MHC/7Vty1ruF4gXqV2N7nO0RtkcMGTKkpCEHzUoY+iBoTm4Hts7Xo4E7KgWStgU+AxxbqD9d0iaSlgYG2p5DEGTC0AdBE2J7GjBX0s3APOApSUfk4l8DHwYmSfpdzjsJOI60nn98vccbNDeLNXoAQRC0j+2DqrKOy/lrt1P3/wFb1WNcjWDMmDEAXHvttQ0eSe8kDH3QN/nvf+GGG2DrrWH5TvctgyZg/PjxjR5CryaWboK+yWOPwW67pTQIWpyY0Qd9kw02gBdfhIEDF143CHo5YeiDvslii8HgwY0eRRDUhVi6CfomTzwB3/hGSoOgxQlDH/RN5s6FBx5IaRC0OLF0UwdmzpzJMccc0+hhNB1jx45tXOcf+xjcccfC6wU15cS7F00epKv1Dxu2YneG07LEjD4IgqDFqZuhlzRQ0rWSJku6XdIISbsUyn/dwX2jJB3bXlknfa0o6c+5r1skfa6r7Uj6pKRvLkp/QS/k3nvh/e9PaRC0OPWc0X8NuNz2KGBz4E1gvqG3/b0S+/o1cEbuaysWIvBUxPY9ts8tcSxBM7LSSvDd76Y0CFqcehr614FNJK1o+x3gK8A2edY9RNItAJI2zdFzJknavXJzfiK4StLHc1n/nH+ZpJUK9foDK9u+GcD2W7Zvz8XrSbo6t7+MpMUk/UnSTTldrDjzl/QdSXfk/tbO47wqvz6zLp9aizJnzhxmzpzZuAGsvDL89KcpDYIWp56G/gLgKZIQ0wTgKuAftkfZnlWodyKwk+0tgb/kvGXz/Yfbvh+4AdhS0kBggO3nC/cPAYrtvQfbOwDjScqAO5OkXkcC/wK+WKmXvzx2AzbNY3mEFMXnhPz6FUmbdNRPUe/79ddfX+iH09eYOnUq48aNa9wAXnsNpk5NaRC0OHUz9Lbftv0z2+sB5wIHd1L3hZzOy1m7AffY/ld+/Udgd5KhviLPtCdLmkwy8h0JbN+X0/+QAjusCUzLeVOAtQp1PwxMs/1uYSzrACfmfrYGhnbyHubrfS+11FIdVeuzDB8+nH333bdxA3jwQRgxIqVB07P2ZqNZe7MIhdtd6uZeKelDwEzbbwPPA68Ay7RT1ZJWsP2ipMoX0R+AVSV9wfZfbT8maSjwJeCrtmcDowp9PSNpc9s3SxoADK+0XRwS8FguuxYYAfy7UP4YMExSP9vz8lgeAi60PTX3E+6p3WTgwIEMHdrh92Tt+djHYMqUlAZNz96n/6nRQ+jV1HPp5pPALXk2fBgpaMJgSZdKKp5FPxy4WtIk0kwekoHeD9hT0hY5bzzwVjby1XwPOCD3NYk0e2+PK4BPSLoJWA+4rFKQl5MuA27LY/kISef7SEkT8/LTBxflAwiaiKWXhuHDUxoELU7dZqS2rwSurMr+bOF6s1zvVlL0nCKTc7prsUnSEk57fb0AfLmdosm5/LxC3u7t1KnUOxOo3nTdub0+g17GM8/AOefAPvs07YaspFNJT5rTitr0OQDJAcDvbR+Z884jLS2+AYyz3e7/RtA36ZVLD5L2Jxnc7Rs9lq4wdOjQxp4CDRbk+efhN7+BHXdsSkMvaUNgadubSzpL0qds35WLzwFuoy3UYIU9bP+bIKiiV56MtX2m7W3yen8QLDobbADPPZfS5mQTUlhAcrpxpcD2c7x3v4n8+vzsPvyh+gwx6C30yhl9byO0bhZOPPEswCDg0Xw9G/jEQuofYvslSZsBp/DeZc4FkLQvsC/Aaqut1sOhBs1Or5zRB0GPefBB2HjjZnavfBmoREUZyEJOd9t+Kae3AB9YWONF998hQzryRg5ahTD0Qd9kiSVgnXVSWiKStpH0dUn9Ja3Tg6Zup20NfjTQqdRmPjyIpLVZBMmPoG/QUoa+WrisI6G0du6bLOnGLLbWwFM8Qd1YfXX4wx9SWhKSziJ5kh2QD9qd3t22bE8D5kq6GZgHPJW9bciie6cAe0g6I99yUZYROYfkvhwE82npNfpFFErbmvQPdQvQwLP5QV145x2YMyfFjF2stH+DtW1vlc9dAPTvSWNFl8rMcTn/XNLp8mLdHXrSV9DatNSMvpqCUNp5kk7LksWd7fq9D1gi36Ps1jYxyysvL+l9WdTsuiyDvHcd3kbL0xCBs3vvhRVWKFumeF7WSHI+BPhOmY0HQXdpaUNfxWTbm9Gx7/0NwEzaDkh9HnjK9lbAb4BvA18AbrO9HfDfzjoLUbOu0xCBszXWgL/8JaXl8QPgz8C6wCXAIWU2HgTdpaWXbqqoCJq9ASDpBtKjdeUE7dbAx0jrm+eQThl+WdJnSZ/T7SRf5em5/j2ddWZ7HHkJaOjQodU+z0GBhgicLb887NqpB2J3GJwnBgB0pm4aBPWkL83o32NsbW+dJZKfLeTdB/TL3hIPAefnOpsBPwEeJ2niAKxfp3G3PA0ROHvxRbjwwpSWx1FVr79fZuNB0F1acUa/h6SNF16tQ84BvgMcBJwuaWLO/xXwV+Avkv4OvArEydzeyhNPwF57JQXLFVboUVOSvgJ8FVhf0lU5e0APRxgEpdFSht72ZJKOfHX+3oXrUe2UjypcTwQqxn0Brx1Ju9h+J7vSPdbTMQcNYtiwFHSkHD/68STtmYOA03LeW8CzHd4RBHWkpQx9nbhW0jLAvwshCjslRM2akH79oKSAMFkqezbwA0nvJwWkEbAybYFtgqBhhKFfRGx/duG1gqbnscfgyCPh2GNL87yR9DNSMPo1gCdJwXW2LaXxIOgBYejrQIiaNY4On6TefhtmzkxpeXzW9kY54M2WdBAvIQjqTRj6oG+y9toweXLZrc7J6bukgPYfL7uDIOgOfcm9MghqzbmSliR5aN1IzOiDJqFhhj4LkD2ZBcVu7aHSX3vtd0nQrIN7l5R0hqRJWTbhjJx/tKR783iPL2+0Qd25+25YbrmUloTtP9t+w/bVtofZ/kVpjQdBD2j0jP6C7Nr4I5LEQGksoqBZNUcBt9veMh+WurhQdojtTYENJEVw8CamUw2dlVeGI44oNYygpC/kicu9kqZLmr7wu4Kg9jTa0FcYCMyRdHiWC/6npGGSVpb0ZwBJi1UOL0naR9LN+WdDSYPzP9gkSafnOhVBs2/ksimSts15CxM529T2hZUXtm9qp850YJWO3lBo3TSeTjV0PvABOPTQlJbHz4A9bW9ge33bcXo6aAoabej3knQT8AeSCNRptrcA9gB+aPsZYGlJy5K0aCZIWhHYERgJ7ESafW9IEi3bknRopcjF+alha+CHhfzORM40/0K6XtJ9klYt5PUHPk0nB6aKEXyWKslfO1g0OtXQefVVuPXWlJbHI8AzZTYo6dQ8oTmtKv8ISTOr4i+smycvt0qKL5lgPo029BfYHgl8EjieNsN/DunQCcDlJIO+O0kZcA1gA2BSLhtE2vjqJ+mPwJ5VfXw2u7tdBRSXWhYQOcsz//dM8WxvC0yhzUPpFJLS5V9sP9/9tx7Umk41dB56CDbbLKXlcT7wiFKA7qsLcgjdQtKGwNK2NwcWl/SpQvE5pAlRkZ8DXwG+lK+DAGge98pXSMs3+wPDgDWBs3PZZcBFwADbj0l6BbjL9q4AkgYA/W0flV/fA1xQaPtwYAuS1vythfwFRM4q15Juk7SX7Uo7xc/pENsTevJmgybg4x+HGTNgzTXLbPV44GvA0yW1twlQ+VubAGwM3AVg+7l2HBgG234aQNJyJY2hQ8aMGQPAtddeW+uugh7SaEO/l1LU+iWAY0kz95vyDwC250iaS/6Dtz0rBwK5ieSvPBG4MXvBDKDtH6PCNbm9O+l6LM1jgP+VtA9Js+QRQrektVhySVh33bJbfSgH5y6LQcCj+Xo28ImF1O/XwfUCKIXM3BdgtdVW69bgxo8f3637gvrTMEOfBcg+VJV9TQfV3yat4Vfu/QNpXb/IZlXtb5bTn5E2yYplexeuR7Uzttdp3wvo6A7G1ymhddOE/Oc/cOaZsP/+sEqHe+qLytKSbiXp2xjA9oE9aO9l0pMuOV3YRGVeB9cLUIyXMGLEiIiX0OI0eka/UCSNA57PG7NBUA4vvZT06L/85TINfdlnK24H9iNNckYD5y2k/kvZaWAe6QkgCIBeYOht1zn0UNAnWG89ePLJUpqS9BHbj5D2mkrD9jRJcyXdDNwLPCXpCNvHSfomaU9rsKTlbR8AjCU5LAg4oMyxBL2bpjf0rUCImjWWOiyb7QycxILxCwz8T08atl3tLnxczj8XOLeq7nSqljC7y4l3v1B63cOGrdjd4QQ9JAx90De5/37YYw+46KLkgdMDbJ+UL8+3PamSL+kzPWo4CEqi0X70pSNpC0kTs0/8DZI+I2l2Qbfmo1X195b0UC7/q6TFc/5DuY1/Sjqwqu5kSX9rxPsLSmLppWHjjVNaHj+ten1wmY0HQXdpqRl9PjV7DLBjdstcFlgLmGF7S0kbAYcC+1TderLtcyQdSQoUcQ0wy/YoSQLuAE4v1q3LGwpqx4c+BGedVUpTipixQZPTUoaeJGdwoe05ALZfAe5OthqAe4BVO7gXkt9yNYsTQcB7LXPmzGHq1KmsvPLK7z0l+9Zb8MILsOKKsPjiPe2mT8aMXXuz0Y0eQtBFWm3pZihZa0TSV/NSzS8L5SOB9s68/0jSgyTNnL/nvCFZOuF+kuRBse5kSSd0NpAQNWsOpk6dyogRIxYUN5sxI7lVzpjR4z5sz7b9JEkW40nSxOAAoKX1ZvY+/U/sffqfGj2MoAu0mqF/hqyRY7uie7MisJ6kSSR3tBMlVYz13vm+k0n/lK8Cg3PerHyYai1gHUmVw10n2x5l+/DOBhKiZs3B8OHDmTJlyoLiZmuuCVddVbYEQuXL/yjSBOE3ZTYeBN2l1Qz9eGDPgs5HZWlqRtaW39n2M7Yrxvq8yo223wLOAL5fbNC2Sf7R7S3rBE3OwIEDGT58+ILiZoMGwQ47pLQ8KmuES+WJRiz5BU1BSxl627NIm7FXZu36M0mKgl3lH8DWkt5HXrrJmjrv2L63/BEHDeOFF+Ccc1JaHo9Iuov09zeAgtx1EDSSVtuMxfZEktBZkerXxfrnFa7nARUp2LU7qxv0cp58Er71LRg2LG3IloDtfSQtZvsdSf2AL5TScBD0kJYz9M1IiJo1IRtuCO+8A/16/lAr6bTCCdYfASfYnpeDhezd4w6CoIe01NJNEHQZCfr3T2nPKXrXbFO4rlZnDYKGEIY+6Js8+ijssktKe05HMr8h/xs0BbF0UwdC1Kx5mL+ENm8evP56SnvOR3NQelVdf6SMxoOgp4ShD/omH/kIXHddWa0VY7de2sF1EDSMXmHoJW1NOoQi4EXSuL9v+99duHd1YI3sjdOVvn5Figv7brcHHPQpbN/Y6DEEQWc0/Rp9Fio7Cvi87ZHAj0n6M11ldWCrLvbVz/bBYeT7ANOmwRJLpLRJkXSqpJuz904xf90s73GrpPVz3nlZaXWypK82ZsRBs9L0hh4YA1yQBcqw/TBtejZHSxqdr8+TtLqknSXdmaWKtycFQN5L0g253lH5n2Firr96lii+FNg7ly2W2z5H0gRJ5+R718z/TFdmCeTV6/5pBN1mzpw5zJw5M71YdVU46aSUloSklUpsa0NgadubA4tL+lSh+OfAV4Av5esKe+QT338saxxBa9AbDP3KZMPeRb4IfMn2VsDfSAGQL7C9taT1gFWyhs0BQEWvZiVgd9u/r2rrX7ZHA6tJGgT8kKRQuEu+p0NC1Kz5mDp1apu42UorwYEHprQ8LiqxrU2ACfl6ArBxoWyw7adt/weoyH0YOF/S1QVdpg4p/n3OmjWrxGEHzUhvMPTzhcraoei+VnGIPhY4UtJ5JEGyIusAo7Iq5VnAwJx/bwfLNffldCbpH+rDwPRc97526rcNLETNmo7hw4e3iZvNmQM33JDS8rhH0g6ShkgaLGnwwm/pkEFAZXCzgeULZf3auT7E9meAXwCnLKzx4t/nkCFDejDMoDfQGwx9RahsWQBJa5Fm+ZD+AVbOwUE+kfOetL0PaSb/A5KwVP9c9hBwfX68HQV8Led35GNX/UXyOEkJs3+hv6CXMHDgwDZxs0cegdGjU1oeI0h/c5cAl9Ezr5uXaZuIDMyvK8yrvrb9Uk5vAT7Qg36DFqTpvW5sz5L0c+CabNBfIgV1ALgcuALYEfhvzjta0sbAMsAhpJn3CZIutr27pGfzjN7An4DrF2E4vyQ9nj+f+wt1wt7KJz6RjHyJa/S2twSQtIztV3vY3O3AfqQvjdHAeYWylyStSjLys3OfA3NUtbV575dCEDS/oQewPYG29cpqhle9/nE7dUYW2joOOK6qfM9C+ah8eXQhb2+ALFi1cZ7R3wo8t/DRB03JEkvAWtUrez0jb/4fAywlaQNgnO3/6U5btqdJmivpZuBe4ClJR+S/37HAn0lPmQfkWy6StDxpAvOdnr6XoLXoFYa+iVgje+AsDZxr+52u3BSiZk3I00/Dr34FBx8MH/xgWa0eAWwGXJcVLFfvSWMFobQKx+X86bmfYt0detJX0Nr0hjX6psH2w7ZH2h5u+8xGjyfoAbNnw9VXp7Q83rX9Jm17O/07qxwEY8aMYcyYMTXvJ2b0Qd9k3XXh4YfLbvVySX8hueNeQNqQDYIOGT9+fF36CUNfB0LUrHmo5RKa7V9J+jvJI+sB2/+qWWdBsAjE0k3QN/nXv2CddVJaEpKus/2A7Utt/0vShaU1HgQ9YJENvaRRko4tvP51uUMCSatkcTEk/UDSTVnb47Scd8sijneNRahf+vsJmpCBA2G77VLaQyR9TNIXgVUk7ZJ/dicCj/RK6rVuXk96vHRj+3tlDATmi4rNI0Xp+Yekz5Fit25h25K26Eazo4BbgMe62H9p7ydoYj74QTj11LJaez+wLunsxno57y3SOY6gl1GvdfN60uOlm8rsOouKnZZn3mNz3lqSrpd0o6Qjc943snDYFEnbFu79DVARCB8F3EgSbTrZtmFBOdiKAFnlOqfHZ1W/SZJWI8XsPEXSKZKWlPSnLGh2saQBkvbO19cC6xfez2RJJ0i6S9I3c96OkqZKGrcoTxVB8zBz5swkbPbmm/D44yntIbZvtH0MsAFJfuDPtk+wfWePGw+CEih7jX6y7c2A7fPr44Bv2t4C+EQ+zXdxPpS0NUkkrMKttrfNp1+XyicLF1XQDGBTYPN8SvFp0onCQ2wfAuwDXJUFzyYDu+Z7XrY9xvY9VW1dQvJX/np+/WPS4atjSLO4DglRs+Zk3LhxSdjsvvtgjTVSWh67kiQ7LpHUX9I1ZTYeBN2lbK+byn/NGzldG7gg2W4GAasAn5J0EOlUX1E6cGpONyCdBIQ2QbOOBEnaEzU7Cfg/SS+SDrAUWQcYLmk/YAmSBMLsQt8LvB/bb0uqaIu8a/s14DVJL3RwTxqYPY6kt8PQoUMjdmiTMF/UbJll4PrrU6Sp8vgfYHNgou13JYWaXdAUlG3oqw3aQ8DBtp/JsgEG7gC2AN5HkhGoUDGm2wL/yNeXAD+U9O28Rj/S9k2FeyqiZm/QJuQ00fa1kn4CfJ4FRc1usH0ZgKQBpDBwXRE1A+iX/3kHASt2+CkETct8UTOAbbYpu/nK31EcmAqaiu4a+j2ycNjCOAL4vaT3kQzuF4FrgJuAO2lffGk4STwM23+TtA5wo6R+pJl30dCPA64mbbZWRLX/WphJ7Qb8P+B4SRsBJwNnS9qf9ARwOIvGSbn/ewidm97N88/DJZfAl75Upib9GcAkYC1J1+fXQZNw4t2dPoR3u/5hw5p/zrfIht72ZJIue3X+3oXrUTn9N/C5qqo/yz/t3gtclj1vKmX/C/xvVf3NcjqetCZaLPtsVX/Pkp4gKuxZVf7PDtoeVf1+gPG2r5K0CnlZJuilPP00fP/7sMkmpRl62xdnA78m8FhFOjgIGk3TnYy1fUmjx9AJu0r6DknU7MCu3hSiZk3I8OHwdrkq00pRyHYEBgMjJVUmKkHQUJrO0Dcztv9MkocNgvb4G3AtydsrCJqGMPR1ILRumouxY8emoCMHHQSnnVam580rto9deLWuI+lUUuSqaUXZYknrAr8l7TV9x/b09vLKHEvQewlDH/RNJBgwIKXlcb5SNLR7yJ43ti/vbmOSNgSWtr25pLMkfcr2Xbn458BXSJ4+ZwI7dZAXBGHogz7KWmvBlVeW3erXgSeBSkAak8JddpdNaIusNgHYGKgY+sG2nwaQtFwnecEisvZmoxs9hNJpKUMvaRQw2vaRWWTqq8BuRS+eHPXnWNvV3jdB0FPezYHpy2IQ8Gi+ns17A9L3a+e6vbx2kbQvsC/AgAEDULlPNu2yqL7Mje7r8A2H1K2v7n7+WR1mobSUoa8gaVNSLM0xRSMfBJD2TAZMn86Qz30OpkxJHjjl8LSkb5Nm3ZWlm2k9aO9loCKvOZD3njuZ1851e3ntUjy5PWLECE+ZMqUHwwyq6aoPfuXL5IRpsxZSs2f++q1o6Ncm6eh8HnhH0g3AANKBqi8VK0q6HZhC0q85IZevBexl+94si/xJ0uxoD9tP5XumkjRwDrN9He1QnDEtt1w8RTcT48aNY6nXXuPQ3/0OVlutzKYXAzbKP5CMfbeCg2duB/YjnRAfTdJtqvBS1o6aR5rtd5QXBC1p6LcFTrX9YhZI+7ztN7KG/la8VzdnBeBY0lH1acDqwIbAN0l+8ofbfl3SaNI/3BH5np+RvjyKipvvIbRumpf5ejdFOYQSsP2NktubJmmupJtJ+k9PSTrC9nHAWJKrr0hPr3SQFwQtaejPAj4jaTuSNMK4fJL1/SQjXzT0s2w/ByDpUdtzJc0Els/lh0rammTUHyjc83y+Z1Dt305QNkOHDk1BwcePh003hZKeuLIsyFiSEJ8AbK/fkzaLLpWZ43L+dNJTZbHuAnlBAK0ZSvAd4MukmfpBwMNZJvky2hQuK7iDa0laARhle3Pgp4V721PMDHob//43jBmT0vI4nfQ39xKwO1C6W08QdIdWnNFj+yVJe5EEz/pJGkFas+xI7rg9/gu8KmkiEAdPWo311kt6N+UJmgHMsf2wJNl+QNImZTYeBN2lpQx9FlybnK8fIG2stseeuc78x9yCmNkTtAmffb6dPor3jOrxoIPGsPjisOqqZbf6tKQlgWmSLqPF/r+C3kv8IdaBEDVrQp56Ck46CQ49tBTPm7zxf3re+P8hKYDOQz1uOAhKoBXX6INg4bz6Ktx0U0pLIMc1PjZfz7N9t+2IIRk0BTGjrwMhatZ8jB07FqaXvvXyhKTvkaKovQs9PjAVBKUQhj4IymMp0jmMDfPrnh6YCoJSCEMf9E1mzIAdd4SrrkoeOCVQ9oGpoPfSVbmCik5OrcMRlrpGL2kLSRMlTZZ0g6TPSJotaZKkWyR9tKr+3pIeyuV/lbR4zn8ot/FPSQdW1Z0s6W+FNi6WtIakC/Prz+T2lqIdJN1S5nsOeinLLw+7757SkpC0saS/SbpX0nRJ4ZYbNAWlGXpJKwLHAF/IbodfAN4AZtjeEjgEOLSdW0/O5VNI8gWQTp+OIsmy7lFVd5Ttz+U+lwDeJAs45S+SU4AvxUZY0Bkz+/Vj5oEHlu1iGQemgqakzBn99sCFtucA2H7F9t2F8nuAzv6r2pMTWBzoLLDnSJLMASSJg/OBPW3PApB0iaQbJV0vaWDxxvxkcJykKZK+JenCPBP7bC4/PN/7T0nDCvecIOkuSd/sZFxI2je3PeX11+M7p9n4/RlncPlxx8Ebb5TZ7BzbDwPK5zjiwFTQFJRp6IcCzwBI+mpeqvlloXwk7fsV/0jSg6QNrL/nvCGSJgP3AzdU1Z0s6YT8ehvg+nz9aeBB248W6u+d5Q8uIc2wqrkU2BQ4kfTE8TnaxKBOy/fuAfywcM8lJD2Rr7fT3nxsj7M9wvaIpZZqdxUpaCD7bb453z3zTLj//jKbjQNTQVNS5h/iMyRjj+0/SroNOBpYT9Ikkpb2/pJ+BIyhTXL1ZNJM/FJgMPAceekmH0K5WNKHKnVtn1Po88O2n8jBRCpyB9+yfbak/sDJktYjaXlf0c6Y77P9tqQHC+JmlUXbvSTtQVoWcjv3hM59L2bIppsmP/qPfnThlbtIZTM2DkwFzUaZM/rxwJ5qC2FW+RKZYXtL2zvbfsZ2ZZ39vMqNtt8CzgC+X2wwH0J5hXaWdSStBDxblf2tPIaRJB35pW2PzG23J0DmqpRCvf2BUblNtXNP0JtZdlnYfPOU9hBJ20iaIelRSXcC68eBqaCZKM3Q53XxY4ArsxDYmaSZelf5B7C1pPeRl24k3QS8Y/vedupvQ1s8zcoY3iQFD/kV8BawlqTrSMs6i8qdwE1AuMy1Is89B6ecktKecyywne01SbEMTuxpg5KWlXS1pFslfa2d8j0k3Sbpmsr+U8ErbbKkj/d0DEHroK7GHGw28qbp7ZXN32YmQrU1IXffDSNHpuWbYcMWVrtTOWpJE21v1dF2oIImAAARiUlEQVTr7iDpB6Qn1ouBSaRYyG/lsgHARGBL4IvAarZPlnRLUXSvq8TfZ+NQjhXbAzvcJan0XrtZZPvvC68VBB0wbBi88kpZra0v6ap8reJr2zt2s81NgANsvyvpXlKIzBm57KOkJdF3JE0gRzIDBuen4AeAg2zP7WbfQYvRaw19EDQRpUUXLzAIqDytzqYt6llnZZvlWAw/IcUrPr2jxlWIabxauXFzgyYkDH0dCFGz5mPwCy/wvfvvh7PO6rHnje0nu3uvpA+Q4rwWeZbkpTYQmJvTlwvllTKKZbZfynlXUOXY0M6Y58c0HjFiRO9cvw26TBj6oE8yr18/WHFFWKyx/wK2nyV5d72HvEa/taRLSB5kRVfNh4F1swvxaOCOLB+i7JCwKfBodZtB36VLXjcN1LDpJ2mYpAn5lOp1ktaUdLSk0V0c+yBJu3T1A5F0mFIw8aCFeXnwYLj4YlhjjUYPpSPOIR3Wuxn4ve03JW0naYztt4Gzc9nXgd+Rlm9uz2v0O5C83oIA6MKMXm0aNjvaniNpWVKIvhm2t5S0EUnDZp+qW0+2fY6kI0kaNtfw3oNQd9C2hvieg1Bq07BZDPgN8EXbz+bH3OVYNAYBuwCXd+G99rPdY9e4oBcwbx7MnZtCCvZrvvg72Zvs81V51xWuLwAuKBTPpk0eOQjeQ1eeWxfQsAHurrgFUVsNm02ASfnxtvKY+2ylb0l7A4vlL5SjSfFi3wJOBV4DLiR9KW2TJRV2A3aiTb7gINvTslfDDOA+SR8j+UVvBnyWti+WMaSTu5eSxNreBk7NcWqDXsbKzz4LSy7JrL/9jSHbbdfo4QRBTenKVKaRGjYrV/peBLYHfpz9mP9A2nD6R1bDNLBjHvNOwFH5nlWB/dqZzT9ve3vgP8D6pKeW3+W8xTsbhELUrKl5edAgrthpJ34/cWKjhxIENacrM/pGath8COgsKkR70gVnAkcqqUv+Gni+UGcNkgbJpKp2HrL9Wjvt35fT/5CeTD5M0tSB9CTT8cAKXg1Dhw4Nr4Ym442llmKjM89ko0YPJAjqQFdm9I3UsLkdGJXX5pG0kqS1C9Vnk2b90PaF8F/b+wM/Ju0tvA30z2WPA3flcY4iPTlA1rNvh+ovkscL/azfwT1BL2CJN95g6B13MHTJJRs9lCCoOQs19I3UsMlfFN8FLpR0I2nzqbi2fwOwXeFUIsB+uf1rSE8Xz5JODF5KCth8raSb8tPIYYvwPiB5Qnwn6+eIzvcZgiZm+f/+F774RXjssUYPJQhqTtNp3TSzho2kfgC250m6FtjX9n8Wdt/QoUO933771Xx8QdfRu+9y1IEHwqBBXfGl75KeSG8ltG4aR5/VumlyDZtlSE8EiwMTumLkAYYOHcrYsWNrO7IgCIIOaD4H4ibG9hzbm9veyPYRjR5P0AOeeAL22SelQdDihKEP+iZvvAHTp5cdMzYImpKmW7oJgrqwzjpw552NHkUQ1IWY0QdBELQ4YeiDvsn06TB0aEqDoMWJpZugb7LiirDvvikNggax/fbb16WfMPRB32ToUDj66EaPIujjXHvttXXpJ5ZugiAIWpww9EHQhEhaVtLVkm6V9LV2yq+U9HIxAI+kPSTdJukaSQOr7wn6LmHog6A5+RbwJ5Kk9j6VKG0Fvg38qvJC0oCcN5KkCRWaG8F8wtAHQXOyCUlm413gXqCo2ort6jgNHyUpyr5DEgXcuLPGi/ESZs2aVeKwg2YkNmPrwNSpU1+V1F5wlkazIvBCowdRRTOO6T7b69a5z0FARdhvNikmbGn1i/ESJM2S9GQ3x1nP31f0tSDX2V5oiLQw9PXhIdsjGj2IaiRNabZxNeuYatj2B4A/V2U/SwroMxCYm9OXF9JUpT5drD8f20O6Wreaev6+oq/uE4Y+CBpIjoM8qjpf0g9IcRwuAT5J++E6izwMrCupPzAauKPkoQa9mFijD4Lm5BxgD+Bm4Pe235S0naQxAJJOB74GnCRpX9tvA2fn+l8HftegcQdNSMzo68O4Rg+gA5pxXDEmkiQ28PmqvOsK1wcCB1aVX0DyuKkn9fxsoq9u0nQRpoIgCIJyiaWbIAiCFicMfRAEQYsThj4IgqDFCUMfBEHQ4oShryOS4vPuhUhaptFjCIKeEO6VNSQb9v1JpxRvtf24pH625zV4aPNplvFIEnAWcB4wJWu2NJT8+zsDuBK4LufJ4aqGpPVIsgsCsH1TjfrZliTWthxpYmrbW5XY/h+Adn+ftv+nrH5yX6t1VGb7qTL7qiYMfW05D3gN+CcwTtJ+th9rpHHNxusA4EXgNttPNIPxsm1JnwZeTcPUNNtvNnJMwC+BR2xfJ+n9wDzbs5rly7FRSLoG+A8wM2cZqImhB04ExrQj4lYWR+f0B6TTxFOBYcCmNejrmJyuAqwMzADWI0lebFOD/uYThr62PAicbvtVSe8Cv5X07Wzs625c86z5IpKRvxf4X0l758M5DUNS/6zSeBmwBrAh8Iakf+UTn43iLuA+SZcDTwKbSNrD9qPN8OXYQGS7XjLItwJv1apx208CSPqk7YNy9sOS9q9BX9/IfV0OfNL2u1my4rKy+6omDH0NKBiBxUgz+T1tX5DsLL+T9GXbLzZgaKOBx23/JI9zdWBJYE4jDVc28gCXA7uRZvWfBw6S9HPbj9VzPIXPYlmS5vvVtn+VA4BcJGnnGs4wewOzJZ0K3Ede9rD9+xr1NRy4TdLz+bVtj6xBPzdKGg9MJ82ya/WEAmk2/ylJ95B0jD5Qw76AOBlbOpKWsz278Pp4YHVgj7w8MRY42/bMjtqo1/gknQn83faVkpay/XqjxpJfrwzsYvsMSVeQNNa/bHtGA8d0OukL8hP593cccKntu+s1pmZD0ter82z/XyPGUhaSVgAGAB8iPb3Ntd1lBdBF7Gt14Eekp9fHgF/afrwWfc3vMwx9OeRlkT+QZvHXA/+2fVsuOx5YDbgd+C6wZVYtrNfYqo3XYrbfyY+n9wD9gc8BPy3Mrms5nurP6kHbd+ayI4GvAlcBU4Dr67G0tJDf3+kkffdbgYOA0bb/U+sxNTOS3gesRNtmbE02EyWtSoqWNbTQV6mbpLmficVNXkmX2P5S2f3ktg+xfUrh9b45PkDNCENfEnmW82ngSGAH0iz+bttX5/LtgMVJRu3hOo2pQ4Oay7cFfk4KVHGg7QfrNK72Pqtptq/Jm8U72b6iHmNZyJjusX1VLh9N+v09YvuReo6t2ZB0GLAt8DHgUeBN26M7v6vbfd0I/Iz0d3oUsL3tH5TY/s7ALqTN0Otz9mLA4K4E9FjEvgYBKwAXkiYzIk2yzrQdm7G9hGkkQzEPuBTYDBiWg1asSPJwqfem59dIXj8V47WdpBVtj8/l/UjrhV+p8zp4R5/VXXk8E6Durp+djWlF4M5Gb1o3ETva/oykybZHSbq4hn3Ns32DpLG2J+QvmTKZSPrdP0WbiuTbJE+YstkC+AJpEvFTkqF/G/htDfp6D3GApzweyT9bkL6lbwOWBtYlGYolGzCmaTmtGK87gPUkrSRpvSx7O7zem510/FmtT/LNXgqgzi6Mzfj7a1YqXjCvSxoJfLyGfd0raQlgoqRJwBslt78m8ILtI3Lb3ybJP5e+QWr7yux58+m8/HQAcLDtmnvdhKEvCdtzgb+SdtF3ID0tPQ+MBG6y/VwDhtWR8doAeL+klWzXPTL0Qj6rmxvxWTXp769ZOSiv0R8C7Ar8sFYd2T7Y9lzbRwNfBHYsuYtf0/blcSFpz+p60pJnqUj6bd4ve1rSt0hPE9dIOmhh9/aUMPQlYvsJ4FySkfg1yUXw/HpscHYwns6M1yTbz3dye63H9gRN9Fk165ialPuAnUneSD8AXqpVR5K2lTRB0gzSXtLpJXfxlu152evmA7Yvtj2R5IFTNuvYnp33zg4FRuUN4F1r0Nd7iM3YGiFpMNDPdr0iyXc2llWArUkbTkOA7zXThmIzfVYVmnFMzYKkPwE3Al/La/UTargZewuwJcn7astq75gS2r+GdGBpE+Bp2z+XtBhwi+2Ny+on93U78Jnc149s75zzb7a9eZl9VRObsTXCds1mOYtKdgU8P/9RN53xaqbPqkIzjqmJGGL7t5Jq4n5Yxbv5x3kmrJLb/wrJaeF22sIwDgWOL7kfSE+Jd5D2oPYGkLQWSQurpsSMPgiCRULSH0lCbwcDJ5EOue1Vo77GkPYA1iGdWj3V9t9q0VcrE4Y+CIJFInvB7EMyvg+STnrPrWF/Ink+vdAomY7eThj6IAi6hOoosyupQ+2cMk/GSlrW9itltdesxBp9EARd5QmSoue/CnkiCZuVLUswhKRBf23+qZUxvhLYCkDSuba/WaN+5iNpOyfp67WA7wMXu0Z6/hXC0AdB0FU2Ip3s/ATwAHCZ7Sm16Mj2DkqRvSpr9IsDF9R4ff7DNWy7yA9JgWx+ApxNchn9VC07DEMfBEGXsH0XaUaPpI2AsyVNsf2tGvX3qqTHgWdIXy5DatDNGpJ+RnoyqVxX+j+qBv0BLJuXwd61fbuk12rUz3zC0AdB0CUkrQjsRDqTMRc4Dhjf6U3d7+s0kuLrncD/AQ/n/LL1j4qSyxNKbLczTgCOBY7NG9t31LrD2IwNmgZJewKjSAE/9rJds8hCwaIj6W3gbmASSe9mvvEoe/abdW2KmLwfUOaBqdzXcqQTsm9I2hQYDFznkqObSVqj+JL3fn411ZsKQx+UjqT7bK/bg/t/SIrq9FCJwwp6iKQtOiqzfWM9x1IWkk4gha7sDzxNCrP5MrCB7d1K7qsSiFyFFNKXV+ka+0Vi6SZoKiQtC3yE/KgeNA+91ZgvhM1sb57jIEyvTFDaeaIog0OAxV0IOiTpA9QwJm6FEDXrY0gaJelaSVdIukfSjpLWkHSrpMlKcTMrdU/OebdI+nS+9zeF8vty2l/SH3Pd43JeP0nnSrop53+0g/HcK+mXkiZmRcSTgSPiYExQJ97JB7IAZinRv0Z9nQ0sU5W3NHBOjfqbTxj6vslQUhDuUSRNj1EkV7lRJMVGJH0O6J/zdgZ+0Ul7OwEv5rrXF/JedwrkfGh790taClgFOC2vu/6U5OJ2rKQ1e/IGg/LJT1v17nO7nK4l6Yysf182N+QfSNLBtdqUXcH2v4sZth8l7QnUlFi66ZtMsf0O8LKkV4GLgSMlXUCSoP0FKeDGZyVNzvcsS2HzqDALgrTUUvEc+GdOP1rJs31nxXBL+i5JlvWe3O+Vtp/O9Y4s+X0G5VL3w0XU2Ofc9pbt5Utaqaw+qtpdyvbrhdfVM/yaEIa+bzI8P54uk3/m2T4cQNI/JF0B3A9cVchfnBQj9IO5jeI/WyW4yUWkcHyVvC2BiyR9mhTtHtu/AX6T2/wOdXAtC2pCvQ4X1dTnXNLyJAVLSOqVqwFHkLR1ti2zL1Lc2/GSziedDVgF2IMUE7emhKHvmzwLXEHbH/UOkg7IZc8Aj9l+WNLmeUZvkj/zYaQ1zUkkPfIKVwK75roVw/1XYEdJN+f7921nHBsAvyvzjQU1pRGHi2rtc/5n4HxSCMtbgCnA0bZLdwbIsW/vIZ32XZ/0v7ZbPWTDw72yjyFpFLCr7e82eixB76Ke7pX18jmXdIPtrfP1PcCwVnQEiBl9EARdwvaNHR0uqkF3P6UDn3PKFVBbs/CUMgg4prL9VMOnlLoTM/ogCLpEnQ8XDaYDn/Myo3+14iGw9ogZfRAEXaWeh4vOBn5clbc0cCawS1mdtJIx74ww9EEQdJUFDhdRu7M47fqc55l+sIiEoQ+CYFG4gbb18om17KhRPuetSBj6IAi6RJ0PFzXM57wVic3YIAi6RGeHi2yXfbgISSuQfM5XJhn78fXwOW9FwtAHQdAlJP2dtsNF+5IOF/2iFoeLgnIJQx8EQZfoK4eLWpFYow+CoKv0icNFrUjM6IMg6BJ95XBRKxKGPgiCoMWJwCNBEAQtThj6IAiCFicMfRAEQYsThj4IgqDF+f8KgFtOiFs4ZwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x180 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = plot_results(loo, summary, chance_loo)\n",
    "plt.savefig('fig_model_fits_exp_srs.pdf', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.0820212806667227"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.log2(np.exp(1.5/-2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAAEWCAYAAAD1m1U5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xm8lWW5//HPl8FZnEADJ8SBUJFsYz9LRRRnc8jUOg7n6CvDNCvL8ujJU1ZmnozSsjqR5+gJsjQ1hxxSESTn2CqCKJLkkDigKTgGxvX7476XLLZrj26e9ey1v+/Xi9e91zNe64F9cT/TfSkiMDOzFfWpdwBmZmXk5GhmVoOTo5lZDU6OZmY1ODmamdXg5GhmVoOTo5lZDU6OZmY1ODmamdXQr94B9Ab77bdf3HzzzfUOw6y3UldWcs+xAC+99FK9QzCzTnJyNDOrwcnRzKyGUiZHSQMk3SBpmqR7JI1eCfsYK+mcFtOGSprcwfWndXdMZlYeZb0h86/A1RHxP5L6Aau3tqCkPhGxrLjQzKw3KGXPEXgT+KikgRHxTkS8JukcSXdJul3SupIulXQRcLOkIZKmSrpT0s/g3Z7hNZKuzdOPljQl90grd69GSbopb3P96gAkfVzSdEl3S9ovTxsv6V5J5xV6NMyscGVNjpOAp4Gpkm6T9FFgWETsAowDFuXl7oqIfYCXgL0jYldggKSt83xFxCHAjcBHImIc8CywY56/WkTsD/wCGF/ZuaQ+wFeBPYGxwNdyD/YzwK7Ade19gZxIZ0iasXDhwi4fCDOrj1Imx4hYGhHfjoiRwP8AlwF353kRy4cvb87tBsCV+TrgrsCQPH12bhe0+Hm9/PODuX0I2KoqhIHACOA24BZgMDAIeCoi3qnab1vfYWJEjI6I0YMGDerQ9zaz8ihlcpS0uaT++eOLwDxg56r5ldPiyrXGo4BrImIscBfLH/qsrgFR/fO7p9VV7RNV818CZgHj8jZHAQuBzSX1ZXnP08waVFlvyHwIuELSW8BS4HjgJEl3Af8ADmux/O3AryQd2sn9LJV0M7Aa8ElgbYCIWCbph8AUSQHMiYjPS7qE1IO9o6tfzKzXW7AAJk6E8eNhyJD2l68TucDWyjd69OiYMWNGvcMwK4dZs+Dgg+G662DkyCL22KXXB8vaczSzRjVyJPz1r/WOol2lvOZoZlZvTo5mVqw5c2CHHVJbYk6OZlastdaCMWNSW2K+5mhmxdpsM7joonpH0S73HM2sWEuWwN/+ltoSc3I0s2LNmgWbbpraEvNzjgUYMmRInHjiifUOw0rqm9/8Zr1DKNaiRXDXXbDLLrDOOkXs0c85mlkPsM46cMAB9Y6iXT6tNrNiLVyYXh8s+WhVTo5mVqynn4YTT0xtifm02syK1dQEPeBeR2E9x1p1YSQdVjX/J62s955aLx3Y10BJv837ulPS/h3djqQPSfpMZ/ZnZo2nyNPqSl2YscButBh6LCK+0I37+gnw07yvPYFXO7piRDwUEf/TjbGYWbW//AUOOSS1JVZkclyhLgzwL8DeuXc3SNKdAJJ2ybVipkr6VGXl3PO8TtK2eV7fPP0qSRtWLdcXGBwRfwKIiCURcU+ePVLS9Xn7a0nqJ+k3uVbMb/Lnd3uYkk7KNWOmShqe47wuf/5ZIUfNGtbixYtZsGABzc3NLFiwoN7hFCcCli4t/al1kclxhbowpDost0bE2Iiovm11HnBIROwB/C5PWzuvf2ZEzAGmAHtIGgD0j4gXq9YfRBq1u6aIOIhUU2Yc8AnSQLZjgEdIA94CkBPuEcAuOZZ5wBnA9/Ln13Jtm5qqa8i8+eab7R4c632am5uZOHEio0ePZuLEifUOpzhbbw033pjaEivshkxELAW+DXxb0r8ApwLvtLLsS7ldlisiHAH8MiIeyYtcBpwJbAz8XtIglifScaQEWUuljsyzwLqk2jAP5GkzgCbghfx5C+CBiPhnVSwjgPPy6OBrAfe38X0nAhMhPQTe2nLWezU1NTF8+HDGjx/P4MGD6x2OtVDkDZmWdWFeA/rWWDQkbZDXqcR3CbBJpQxCRMwnFdE6knQdc2HugY7Nyew5SbvlbfSXVKk/07KOzHxSQgQYzYp1ZOYDO1ZiyO1c4Ct5P6OBa7t0MMyAAQMGMGTIEJqamhhS4nIB3a65Gfr3T22JFXla/SHgzlwh8AzgHGB9SVe2qBl9JnC9pKmkHiOkpHYicIyk3fO0G4ElEbGI9/oC8Pm8r6mkXmItvwe2kzQdGAlcVZmRT/WvAu7OsWwNnAucletc3wZs2pkDYGak96p/9KPUlliPfbda0snAwoj4XbsL15nfrba29Lp3q4vXpXere+QbMjkxfgK4pt6xmFknLV4Mt96a2hLrkckxIn4WEXvnmzxm1pPMmwf77JPaEvPrg2ZWrO23h/nzS12zGpwczaxoq64KW2xR7yja1WNvyPQko0ePjhkzZtQ7DLNyeOYZ+OEP4StfKeqOde+5IWNmPdjixXDzzaW/IePTajMr1nbbwaOP1juKdrnnaGZWg685FsAPgVtXNOzD4bNnw2GHwdVXpzvXK5+vOZpZD7DOOnDQQUVVHuwyX3M0s2JtuilMmFDvKNrlnqOZFevtt9Mo4G+/Xe9I2uTkaGbFeuSRNNDtI4+0v2wdNVRybFlEq7WiXTXWmybpjlz4a/zKi9DM2HpruO02jwReT50s2jUOWAbcSR7B28xWggEDYNy4ekfRrobqObZUVbTrUkkX5jKtbT0fsSqwWl5Hkn6eB7a9QdJ6klbNBbZuzqVfjyvga1gvtHjx4sYtvPXii/DjH6e2xBo6ObYwLSJ2BQ5oZf4UYAFQqSr4ceDpiNgTuAj4HHAocHdE7Ae80tbOXGDL3o/m5ubGLbz1t7/B6aentsQa+rS6hUpxrbcAJE0h1bD5dJ4+DvggqYTDxcAI4NOS9iUdp3tI5Roezss/1NbOXGDL3o+mpiYmTJjQmIW3Pvzh0t+pht6VHFdIUBHx7kWPXOGQiJgtqU+uMjgX+FVETMjL9AcOI9WauRHYAfhzMaFbbzNgwACampraX9BWmkY8rT5a0m25AFZXXAycRKqrPTRfc7wd2J9UlmEXSX8EPgB4JHKzzpo3D/bbzyOBFykippHqTbecflzVz2NrzB9b9fPtwO3543vudks6LCLekfRzUvlWM+uMPn1gjTVSW2INlRwLcoOktYC/RMQ99Q7GrMfZcss06ETJOTl2UkTsW+8YzHq0CFi2LPUc1aUBcwpR7n6tmTWeBx6Afv1SW2Iez7EAriFjVuWll+Caa+DQQ2HgwCL22KXuqU+rzaxYAwfCCSfUO4p2+bTazIr16qtw/fWpLTEnRzMr1hNPwMEHp7bEfM2xAK4hYz1dt9azWbIkXXccOBBWWaX7tts6X3M0sx5glVVgyJB6R9Eun1abWbGeegpOOim1JebkaGbFeuMNuPfe1JaYT6vNrFjbbgsPPljvKNpVt55jrvfyVK7fclceJqw7t9+h+jGtrLu6pJ9KmppHD/9pnn62pJk53nO7L1ozK5t6n1ZPyiPifI000na36WT9mJa+AdwTEXvk0cMvr5p3WkTsAoyStOn7CtKsN5o1CzbfPLUlVu/kWDEAWCzpzFwF8D5JO0oaLOm3AJL65XEVkXSCpD/lPx+WtH7ugU6V9OO8TKV+zPF53gxJ++Rp7dWU2SUiJlc+RMT0Gss8DGzcvYfBrHy6vZ7N+uvDMcektsTqnRyPlTQduAS4ArgwInYHjga+GhHPAWtKWptUxuA2SQOBg4ExwCGkXt6HSTVi9gC+1GIfl+fe6Tjgq1XT26op8+5zUZJukTRb0iZV0/oCH6GN8RxdQ8YaRbfXs9l4Y/jud1NbYvW+ITMpIs6StBFpBO4/SDqaVCK18nT61aQkuCdwDjAMGAVMrdrOHcAYSZcBNwGTqubtK+lLpIS3YdX09mrKABAR+0i6lOXHagKpuNZvI6LV8mmuIWONotvr2bz1Vno7ZsstYfXVu2ebK0G9k2PFa6RT65OBHYEtgV/meVcBvwb6R8R8Sa8Bf46Iw+Hd2i59I+Ib+fNDrJgczwR2J5Vdvatqels1Ze6WdGxEVLZTfZxOi4iulmAw63G6vZ7NnDkwejTMmAElrpNT7+R4rKRdSbWizyH1EKfnPwBExGJJbwO35c8Lcx3p6cA/SSUN7sh3j/tXlqvyh7y9+4GOvun+LeCHkk4AlgDzgOe79hXNbAXDh8Odd6a2xHrEu9X5dPm0fA2yx/G71dbTdeu71cXr0rvV9b4h0y5JE4EXe2piNLMWnn8evv/91JZY6ZNjRIyPiFPrHYeZdZPnnkt3q58rd3+n3tcczay32XFHWLSo3lG0q0dcc+zpXEPGrK4a85qjmTWYuXNh7NjUlpiTo5kVq3//NNht//71jqRNvuZoZsUaNgwuu6zeUbTL1xwL4OccrdG8r+cely2Dt9+G1VaDPoWcvPqao5n1AA8+CGuuWfoBb50czaxYQ4fCpEmpLTFfczSzYm2wQRrPseTcczSzYr3yClx5ZWpLrOGSo6TdJd2eR/+eIuljkhZV1YPZpsXyx0mam+dfI2mVPH1u3sZ9kr7YYtlpkm6qx/cz6/Hmz4cjjkhtiTVUcsyjhH8LODSP/n0oaTDbWXmU8NOA02usen6ePwPYJ09bmLexM2lk8uplx0bE/ivnW5g1uFGj4OWXU1tijXbN8QBgckQsBoiI14AHpXfv5D8EbNLKugDr1pi2CrC0O4M069X69St9/RhosJ4jMAR4DkDSUfk0+gdV88cAtd5Z+pqkx0i1aP6Ypw2SNA2YA0xpsew0Sd/r9ujNeoBKwa0uF9168kk4/vjUllijJcfnSAmSiLgMOAYYCIyUNJVUhuE8SZUEd1xe73xgB+B1oPJfWuW0eitghKTNK8vm0+oz2wrEBbasUVUKbnW56Nbbb8Ojj6a2xBrttPpG4EpJV0TEIpZ/v8o1x4rz8x8qCTIilkj6KfBl4IzKghERuW5NrVPuVrnAljWqSsEtoGtFtz74Qbj33m6Oqvs1VHLM9WW+BVwraRnwDnAe8O0ObuJW4BxJq7L8tLoP8GhEzJS048qI26wn6faCWyXVUMkRICJuJxXdqtbyc/Xyl1b9vAzYKX98T/Wf6mXNrItmzoR99oFbbin1HetGu+ZoZmW34YZwyimpLbGG6zmaWckNHgz/+Z/1jqJd7jmaWbHeeAOam1NbYk6OZlasxx6D0aNTW2Ie7LYALrBlVuWNN1Ji/OAH07iOK1+XBrv1NUczK9aaa0IPeBTIp9VmVqznnoPvfCe1JebkaGbFevFFuOii1JaYrzkWwAW2rJG9r2JbxXCBLTOz7uLkaGbFeuwx2Hnn0j/K4+RoZsVabTUYMSK1JeZHecysWEOHwiWX1DuKdjVUz1HSWEnn5J8/KekqSX1aLDNU0uT6RGhmvPMO/P3vqS2xhkqOFZJ2AT4PHJOHITOzspg5M9Wunjmz3pG0qRGT43DgB8ARwDu5POv03IvsW72gpHsk/UTSTEmflnS1pIcljcrzL5R0h6Q/Sdqsap2LJD0kab/Cv51ZSSxevJgFCxZ0vpbMsGHwu9+ltsQaMTnuA/wxIl4mjQT+8YgYAzwK7Nli2Q2Ac4D9gQuAo4DPAZ/J88+MiN1J5V5PrFrn28CBVdPewzVkrNE1NzczceLEzteSWW89OPzw1JZYI96Q+TnwsdyruxOYKGljYCNgXv5TsTAiXgCQ9EREvC1pAVD5Wztd0jigPym5VtZ5Ma/Tal0Z15CxRtfU1MTw4cMZP35852rJvPwy3HQT7L9/Or0uqUZMju8AnwZuAX4PPB4RR0n6Lu99Uj5a+VmSNgDGRsRukvYGjq61XPeGbtZzDBgwoGv1ZJ58Eo49FmbMcHIsWkT8XdKxwPVAH0mjgUWs2GtszyvA65JuBx5eCWGa9U477piGLSv5c45+t7oAfrfaGpnfrTYz6w7z58NRR6W2xJwczaxYS5fCggWpLbGGvOZoZiU2fDhMm1bvKNrla44FcA0Zs7ryNUcz6wEefBDWWSe1JebkaGbFGjwYvv711JaYrzmaWbE+8AE4/fR6R9EuX3MsgJ9ztN6gw887vv56GpFn1ChYa62VG1Tia45m1gPMnQu77praEvNptZkVa9ttYdYs2HLLekfSJidHMyvW6qvD9tvXO4p2+bTazIr17LPpbvWzz9Y7kjb1iOQoaVwekXu6pN9Lul7SVh1cd6ikloPctrX8BS1HDDezbvT3v8PkyaktsdKfVksaCHyDNKL3a5K2AX7SiU0MJY0AfnsH9tUnIk7tUqBm1jEjR8JTT9U7inb1hJ7jgcCkiHgNICIeB54DkHS2pL3yz5fmXuInJN0v6XZJBwDjgWMlTcnLfUPStDx/aP4zVdKVwHF5Xr+87Ysl3Sbp4rzulpLuk3Rtrk0ztPCjYWaF6AnJcTA5GXbQJ4EjI2JP4CZSqYJJETFO0khg44gYS6pOeGZeZ0PgUxHxvy229UhE7AVslksifBX4EnBYXsfMSMW2mpubO1Zsa86cNODtnDnFBNdFPSE5PgcMaWVerZIF5wBnSboUaHldcgQwVtI0Uq2ZAXn6zIj4Z43tz87tAmAdYAvg4bzs7BrLLw/GBbasF2lubmb06NEdK7a15pqw886pLbHSX3MEbgSulPTbfM1xK1JvElLpg8GSBGyXpz0VESdI+hjwFWAyULnBMhe4JSK+ACCpP7Ax0Fpt65bJ96/ASEkzqvZXe0UX2LJepKmpiQkTJgC0X2xr883h5z8vIKr3p/TJMSIWSvoO8IecBP8OLMmzryYV0TqYVPMF4GxJOwNrAaeRenjfk3R5RHxK0vO55xjAb0iFuDrqB8CvgRfz/so9WqdZQTpVaGvJEnjpJRg4EFZZZeUG9j743epOkNQvIt7Jj/rcBewaEe+0t57frbbeoMPvVjc3w+jRqfpgZysXdo3frS7AMEnTgfuBX3UkMZpZC1tuCddd59cHG0l+jGhMveMw69HWXRcOOqjeUbTLPUczK9ZLL8HFF6e2xJwczaxYTz0Fn/1s6d+S8Q2ZArjAllmVCFi2DPr0AXXpXklndWknvuZoZsWSoG/5x3bxabWZFeuJJ+Cww1JbYk6OZlasZcvgzTdTW2K+5lgAPwRuvUmHHwYvjh8CNzPrLk6OZlasBx6A1VZLbYk5OZpZsTbZBL7//dSWmB/lMbNibbghfPGL9Y6iXZ3uOUoaK+mcqs+dqefS0X1sLOmC/PNXcmGtOyVdmKfd2cl4h3Vi+W7/PmZWZfFimDIltSX2vk+rKwPHdgdJlXj2Bm6VtD8wHNg9InYljd/YWWOBDiXHXGCr276PmdUwbx7stVdqS+x9J8dKLy4XuLow9/C+madtJemWXFb1rDzt+FzEaoakfarWvQi4OW92LHAHcCRwfuTnjSLijhb7niapX+Xn3J4r6a5cNGsz4DhggqQJklaX9JtcXOtySf0lHZd/vgHYoer7TJP0PUl/lvSZPO1gSc2SJnam92rWW3Solsx226XEuF2bg+nXXXffkJmWe3gH5M/fBT4TEbsD20naBLg8F7gaRypYVXFXROyTR/teIyJep/PFtQB2AXaLiD2AZ4BLgdMi4jTgBOC6XHxrGnB4XufViDgwIh5qsa0rgF2Bf8uf/500ZNm3gI3aCsI1ZKw36lAtmdVWg622Sm2JdfcNmUrRqbdyOxyYlPId65Lqtewk6UukBzOrK/g153YUMDP/XCmu1Vr/u1aBre8D/yfpZeDrLZYfATRJOhFYjVQmYVHVvt/zfSJiqaTKo/z/jIg3gDcktTnekmvIWG/UoVoyzzwDF1wAp54Km25aYHSd093JsWUSmAucGhHP5dICAdwL7A6sSio1UFFJQPsAt+afrwC+KulzERGSxkTE9Kp1KgW23gI+kKfdHhE3SPoP4OOkOi/VBbamRMRV8G6BraPpWIEtgD6S1iAl+oGtHgWzXqpDtWQWLYLrr4fjj2/I5Hh0LmLVnq8D/ytpVVKS+iTwB6BSauDVGus0kQpZERE3SRoB3JFv1jTndSsmAtcDdwIL87RrcgIDOAL4G3CupP8HnA/8UtLJpJ7mmXTO9/P+HwJe6OS6Zgaw/fbw+OP1jqJdpXu3WtKREXFFveOoparA1sbAxIg4sCPr+d1q6038bvVKUtbEmB0u6Q7gWtLNJjPrrEcegREjUltifkOmEyLit8Bv6x2HWY82YADst19qS8zJ0cyKtemm8KMf1TuKdpXummMjcg0Zsyr/+AcsWABDhsCqqxaxx8a45mhmDW72bBg2LLUl5uRoZsXaemu45ZbUlpivOZpZsQYMgL33rncU7fI1xwL4OUfrzd7z3OOLL8IVV8CRR6axHVc+X3M0sx7gmWfgy19ObYn5tNrMitXUBEuX1juKdrnnaGZWg5OjmRVr3jw44IDGHwnczKxTJOjfP7Ul1q3JUdLuuQTBNElTJH1M0qJcsuBOSdu0WP44SXPz/GskrZKnz83buE/SF1ssO03STVXbuFzSMEmT8+eP5e2tQQ0ub2BWZ1ttBddem9oS67bkKGkgqXzAobkMwqGkEcFn5ZIFpwGn11j1/Dx/BmmgW4CFeRs7kwajrV52bETsn/e5GvAP8mC1OflOAI6MCNcmMKuzSk2ZVuvJlFh39hwPACZHxGKAiHgtIh6smv8Q0FYV73VrTFuFNEhua8aQBrqFVNPlV8AxEbEQQNIVubjXLZJWGAIk90C/m+u8fFbSZEkzJe2b55+Z171P0o5V66xQdKs1riFjtrymzAr1ZJqb0yl1c2vVScqhO5PjEHIxLElH5dPoH1TNH0MqU9DS1yQ9BnwY+GOeNihXE5wDTGmx7DRJ38uf9wZuyT9/BHgsIp6oWv64XNzrCuBTNfZ9Jakg13mknu3+wOfzvAvzukezYiGwlkW3aoqIiRExOiJGr7FGzTN8s4bX1NTEjBkzGD9+/PKJm20Gv/hFakusO59zrBTDIiIuk3Q3cDYwUtJUUkmEkyV9DTiQVBUQUumCX5ES1fqk8gMLI2JsrkR4uaTNK8tGxMVV+9wiIp6UNJRULqGPpM9GxC9zzZrzJY0EBgC/rxFzpYDWYxHxAoCk9fK8YyVV6stEjXVaqztjZlnNmjKDBkF1siyp7uw53ggcI2md/LmSeGdFxB4R8YmIeC4iKtcNL62sGBFLgJ8CX67eYK5X/Ro1TrklbQg832LyZ3MMY4APAWtGxJi87Vq3xqJFS9VyJ5PqZ3+2xbp+39Ls/Vi0CG68MbUl1m3JMV/n+xZwraTbgZ+ReoQddSswLhfjGpRPn6cD70TEzBrL7w3c1iKGfwBHAhcAS4CtJN1MOuXurPtJxbSO78K6Ztaav/wFDjwwtSXWYweeyDdO7qncACozDzxhvdl7Bp5YsiQNPrHhhrDKKkWE0KUHKnvsu9UR8cf2lzKz0lllFdikrQdXysFvyJhZsZ5+Gk45JbUl5uRoZsV6/XWYPj21JdZjrzn2JC6wZVZXHuzWzKy7ODmaWbFmzYIttkhtiTk5mlmx1lsPPvWp1JaYrzkWwM85mq3oPc8+rly+5mhm5ddv6VJ49FF46616h9ImJ0czK9SghQth221hzpx6h9ImJ0czK9TLG2yQnnPcZpv2F66jHvv6oJn1TEtWXRV2263eYbSrQz3HOtaG6SNpR0m35VG5b5a0paSzJe3VwdjXlXRYRw+IpDMkbdzR5c2sc9Z8/XWYMAFeeKHeobSp3eRY59ow/YCLSKUPdgeOo/O93XWBDiVHSX0i4ryIeLaT+zCzDlr7tdfg7LOh5HVlOpJo3lMbBnhQy8sqrszaMB8FpkbE83nfzwPPV/Yt6TigX0RcLOlsYBppHMcfAW8Ak4GtgL1z2YUjgENYXuLgSxHxgKSZwCxgtqQPAueQSiHsC1QG7z2QNFL5laT/HJYCP4qIaW18DzNr4fE116R52jQGb7RRKh1QUh05ra5nbZjBlX13wgHAv0fEnsAlwETg1txjDeDgHPMhwDfyOpsAJ0bEeS229WJEHAA8C+wAnAD8Ik9rcyA6F9gyq61m0a0S6khyXKE2DHAMMJDltWFOBs6TVElwx+X1zicllNdJPS5Yflq9FTCiRW2YsRFxZv68RUQ8Wb3vVtQqb/Az4EhJk4CdWiw/DBgFTAWuZnmvdm5EvFFj+7Nz+2xedgvg4TztoTbicoEts1bsO2wYi3baiZPGjat3KG3qSHKsZ22Ye4Cxkj5QmSdpeNXii0i9S4CRuX0lIk4G/p10rXQp0DfP+yvw5xznWFIPFXLd6xpaJt+/Vu1nh1bWMbM2rLH22gzYYgs22rjc9z3bTY71rA2Tk+spwGRJdwCTWPFa5RRgP0nXVU07MW//D6QKh88D60u6EvgncIOk6bnXe0YnvgfAxcBJuS6NaPu6qZnV8Or668Pll8OwYfUOpU2le7e6zLVhJPUBiIhlkm4AxnfkzrbfrTarsmwZ3zzzzFQuoU8h76E0xrvVEfHHMibGbC3gDkn3AQ/5kR+zzhv8/POw+urw4IP1DqVNfkOmE3LSLv+j/WYl9uq668Kll8LQofUOpU1OjmZWqLfWWAP+7d/aX7DOSnfNsRG5hoxZlVdegalTYY89ihrwtjGuOZpZg5s/Hz75ydSWmE+rzaxYo0bBwoWwbq03i8vDydHMitWvHwwcWO8o2uXTajMr1pNPwgknpLbEnBzNrFhvvQUPP1z6GjI+rTazYo0YAfffX+8o2uWeo5lZDU6OZlashx+GIUNSW2JOjmZWrIEDYfz40t+x9jVHMyvWkCGphkzJuedoZlaDk6OZWQ1OjmZmNTg5mpnV4ORoZlaDk6OZWQ1OjmZmNTg5mpnV4DIJBZD0GjC33nH0EAOBl+odRA/g49Rxq0XE9p1dyW/IFGNuRIyudxA9gaQZPlbt83HqOEldKuDk02ozsxqcHM3ManByLMbEegfQg/hYdYyPU8d16Vj5hoyZWQ3uOZqZ1eDkaGZWg5OjmVkNTo5mZjU4ORZIko93B/lYWXeQtFZX1/U/wJVIUh9Jp0g6RtIWEbHMv/S15WP1BUlHSRqaj5XqHVfZKPlvSTtL8hturcj/nn4O7Fo1rVPET4Z6AAAFCElEQVT/nvyLunJdCmxHek1zoqRhTpDvlf/R/hoYDqwJ/FDSgPBzZu+Rj8lHgMOBnSStWueQyuoHwLyIuFnSRpIGRUR05nfP//OsXI8BP46I1yX9E/hvSZ+LiPmS5F/+d+0F/DUi/gNA0lBgdWCxj9NykvpGxD+Bq4BhwIeBtyQ9EhFL6xtd6fwZmC3pauAp4KOSjo6IJzr6b8o9mJWgqvte6TH2iYhJwCTgF5I28C/8chFxK/BfVZPWA3bOP69efETllBMjQOUX/nXg4+SzkroFViJVv3trAxcA0yPiy8DPgF9LGtzR3z0nx24maZ3KwY+IbwNPApPz/1aTgDsBnwqRjlXl54hYVHUNbTawUNJuwFmS+tYlwJKoPk7Zq8DLEfF/QBPpNHvNwgMrmRa/exOBR4Dx+XfvV8AU4AMd3p47MN0j/491Cam3eAvwl4i4O887F9gMuAc4BdgjIp6vV6z1VuNYPRYR91fN3wf4DrAI+GJEPFaXQOusreMk6SzgKOA6YAZwS0Qsrles9dTO796PSWcidwFfAvaKiGc7sl33HLvPvwJvAF/In/eSdBBAvpY2GXgGOKQ3J8as5bHaT9IBVfP7AIOBz/XWxJjVOk4fzz+fC3w9Is6IiCt7a2LMav3uHQwQEV8E/g94Gji4o4kRnBy70wO5XQZcCdwLbCtpsKSRwN0RcV1EPF63CMuj1rEaKWlDSSMj4magKSLm1y3Ccqh1nLaTtBGwA3Ab+JlQah+nEVW/e/dHxI0RMa8zG+3tB7U7zct/dgf6AneTrgNtTxrS3jcWlmvtWI0CNpK0YUQsrGN8ZdHacdoBWBdYAyAiltUrwJJYKb97To7dJCLeBq4BPgQcRLr+8SIwhnTH7IU6hlcq7RyrqRHxYh3DK412jtOf/G8qWVm/e74h080kbQyMA/YGBgFf6Gx3vrfwseoYH6eO6e7j5OS4kkhaH+gTEa4Q1w4fq47xceqY7jpOTo5mZjX4mqOZWQ1OjmZmNTg5mpnV4ORoZlaDk6OZWQ1OjmZmNXiwW7OCSDoGGEsaa/DYiFhS34isLe45mnWSpNldWS8iJkfECaRRqrfo3qisuzk5mhVI0trA1oBHZyo5n1ZbQ5E0FvgasITUO/sGaWTxScBS4M2IOCAvez6wE+n34CukUW4Oj4hT8vzZEbF9Hol8EjCENGhqZZiwX5IS3TJgfK3h6CTNBG4l1XvZHzifNA6jX00rOfccrRENAY4gXd87N7dXRcRYUs0VJO0P9M3TPsGKNWxaOoRUlmAsaaTpyrQ3I2IMcHqt9SWtAWwMXBgRewL/SUrY50ja8v18QVv53HO0RjQjIt4BXpX0OnA5qRbNJFIv8r9IY/3tK2laXmdt4N3eXIsax1uTBlAFuC+321SmRcT9lWQn6RRS2dSH8n6vjYhn8nJndfP3tJXIPUdrRE2S+ubCVGsByyLizIg4ljSE/jbAHOC6iBibe4QfBV4BNs3b2Klqe/NIRayoaueRKyRK+ggwHyAiLsrbPJU0vuC9WI/knqM1oueB35OKmn0dOEjS5/O854D5EfG4pN1yzzGA+4EzgHckTQXuqNretcDhedlKsrsGOFjSn/L642vEMQr4RXd+MSuOhyyzhpJvyLx7U8Wsq3xabWZWg3uOZmY1uOdoZlaDk6OZWQ1OjmZmNTg5mpnV4ORoZlaDk6OZWQ1OjmZmNTg5mpnV4ORoZlbD/wdhdjAyXmAfKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "loo = pd.read_pickle('Data/model_fits/model_fits_exp_cp.pkl') \\\n",
    " + pd.read_pickle('Data/model_fits/model_fits_exp_srs.pkl') \\\n",
    " + pd.read_pickle('Data/model_fits/model_fits_exp_lin.pkl') \\\n",
    " + pd.read_pickle('Data/model_fits/model_fits_exp_scram.pkl') \\\n",
    " + pd.read_pickle('Data/model_fits/model_fits_exp_shifted.pkl')\n",
    "loo.sort_values('LOO', ascending=False, inplace=True)\n",
    "\n",
    "#note the \"Simple Kalman\" is modeled without sticky choice and wasn't included in the paper\n",
    "loo.drop('Kalman w/SC', inplace=True)\n",
    "\n",
    "\n",
    "\n",
    "n = len(pd.read_csv('Data/exp_srs/gprbfsrs.csv')) \\\n",
    "    + len(pd.read_csv('Data/exp_shifted/gprbfshifted.csv')) \\\n",
    "    + len(pd.read_csv('Data/exp_scrambled/gprbfscrambled.csv')) \\\n",
    "    + len(pd.read_csv('Data/exp_changepoint/changerbfpred.csv')) \\\n",
    "    + len(pd.read_csv('Data/exp_linear/rbfpred.csv'))\n",
    "# n_subj * n_trails_per_subj\n",
    "n_arms = 8\n",
    "chance_loo = -2 * n * np.log(1. / n_arms)\n",
    "\n",
    "_loo = loo.sort_values('LOO')\n",
    "_loo[r'pseudo-$r^2$'] = 1 - _loo.LOO / chance_loo\n",
    "_loo['r2_err'] = _loo.LOO_se / chance_loo\n",
    "fig, ax =  plt.subplots(1, 1, figsize=(4, 4), gridspec_kw=dict(wspace=0.5))\n",
    "min_r2 = np.max(_loo[r'pseudo-$r^2$'])\n",
    "with sns.axes_style('ticks'):\n",
    "    plt.sca(ax)\n",
    "    ax.barh(range(len(_loo[r'pseudo-$r^2$'].values)), \n",
    "             _loo[r'pseudo-$r^2$'].values, color='grey', align='center')\n",
    "    ax.set_yticks(range(len(_loo.LOO.values)))\n",
    "    ax.set_yticklabels(_loo.index.tolist())\n",
    "    ax.errorbar(y = range(len(_loo.LOO.values)), \n",
    "                x = _loo[r'pseudo-$r^2$'].values,\n",
    "                xerr = _loo['r2_err'].values, linestyle='None',\n",
    "                color='k',\n",
    "               )\n",
    "    xlb, xub = plt.gca().get_xlim()\n",
    "    ylb, yub = plt.gca().get_ylim()\n",
    "\n",
    "    ax.plot([min_r2, min_r2],[ylb, yub], 'r:')\n",
    "    ax.set_xticks([0., 0.2, 0.4, 0.6])\n",
    "    sns.despine()\n",
    "\n",
    "\n",
    "    ax.set_xlabel(r'pseudo-$r^2$')\n",
    "    plt.xticks(rotation=45)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}