{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import sys\n",
"sys.path.append('../..')\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"import os\n",
"\n",
"from definitions import RESULTS_FOLDER, FIGURE_FOLDER\n",
"import statsmodels.formula.api as smf\n",
"\n",
"\n",
"figure_location = os.path.join(FIGURE_FOLDER,'twostep')\n",
"if not os.path.exists(figure_location):\n",
" os.makedirs(figure_location)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"data_lesion_hpc = pd.read_csv(os.path.join(RESULTS_FOLDER, 'twostep', 'results_lesion_hpc.csv'))\n",
"data_lesion_dls = pd.read_csv(os.path.join(RESULTS_FOLDER, 'twostep', 'results_lesion_dls.csv'))\n",
"data_control = pd.read_csv(os.path.join(RESULTS_FOLDER, 'twostep', 'results_control.csv'))"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Unnamed: 0</th>\n",
" <th>Action1</th>\n",
" <th>Action2</th>\n",
" <th>Agent</th>\n",
" <th>DLS reliability</th>\n",
" <th>HPC reliability</th>\n",
" <th>P(SR)</th>\n",
" <th>Reward</th>\n",
" <th>StartState</th>\n",
" <th>State2</th>\n",
" <th>Terminus</th>\n",
" <th>Trial</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.057327</td>\n",
" <td>0.087327</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>3.0</td>\n",
" <td>5.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.115198</td>\n",
" <td>0.167028</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>3.0</td>\n",
" <td>5.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.175993</td>\n",
" <td>0.239769</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>3.0</td>\n",
" <td>5.0</td>\n",
" <td>2.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>3</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.224014</td>\n",
" <td>0.306158</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>3.0</td>\n",
" <td>5.0</td>\n",
" <td>3.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>4</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.252047</td>\n",
" <td>0.366749</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>4.0</td>\n",
" <td>7.0</td>\n",
" <td>4.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Unnamed: 0 Action1 Action2 Agent DLS reliability HPC reliability \\\n",
"0 0 0.0 0.0 0.0 0.057327 0.087327 \n",
"1 1 0.0 0.0 0.0 0.115198 0.167028 \n",
"2 2 0.0 0.0 0.0 0.175993 0.239769 \n",
"3 3 0.0 0.0 0.0 0.224014 0.306158 \n",
"4 4 0.0 0.0 0.0 0.252047 0.366749 \n",
"\n",
" P(SR) Reward StartState State2 Terminus Trial \n",
"0 1.0 1.0 0.0 3.0 5.0 0.0 \n",
"1 1.0 1.0 0.0 3.0 5.0 1.0 \n",
"2 1.0 1.0 0.0 3.0 5.0 2.0 \n",
"3 1.0 0.0 0.0 3.0 5.0 3.0 \n",
"4 1.0 1.0 0.0 4.0 7.0 4.0 "
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_lesion_dls.head()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def is_common_or_rare(action, out):\n",
" left_outcomes = (5, 6)\n",
" right_outcomes = (7, 8)\n",
" if action == 0 and out in left_outcomes:\n",
" return 'common'\n",
" elif action == 0 and out in right_outcomes:\n",
" return 'rare'\n",
" elif action == 1 and out in left_outcomes:\n",
" return 'rare'\n",
" elif action == 1 and out in right_outcomes:\n",
" return 'common'\n",
" else:\n",
" raise ValueError('The combination of action and outcome does not make sense')"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def add_relevant_columns(dataframe):\n",
" dataframe['PreviousAction'] = dataframe.groupby(['Agent'])['Action1'].shift(1)\n",
" dataframe['PreviousStart'] = dataframe.groupby(['Agent'])['StartState'].shift(1)\n",
" dataframe['PreviousReward'] = dataframe.groupby(['Agent'])['Reward'].shift(1)\n",
" dataframe['Stay'] = (dataframe.PreviousAction == dataframe.Action1)\n",
" dataframe['Transition'] = np.vectorize(is_common_or_rare)(dataframe['Action1'], dataframe['Terminus'])\n",
" dataframe['PreviousTransition'] = dataframe.groupby(['Agent'])['Transition'].shift(1)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"add_relevant_columns(data_control)\n",
"add_relevant_columns(data_lesion_dls)\n",
"add_relevant_columns(data_lesion_hpc)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Unnamed: 0</th>\n",
" <th>Action1</th>\n",
" <th>Action2</th>\n",
" <th>Agent</th>\n",
" <th>DLS reliability</th>\n",
" <th>HPC reliability</th>\n",
" <th>P(SR)</th>\n",
" <th>Reward</th>\n",
" <th>StartState</th>\n",
" <th>State2</th>\n",
" <th>Terminus</th>\n",
" <th>Trial</th>\n",
" <th>PreviousAction</th>\n",
" <th>PreviousStart</th>\n",
" <th>PreviousReward</th>\n",
" <th>Stay</th>\n",
" <th>Transition</th>\n",
" <th>PreviousTransition</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.057327</td>\n",
" <td>0.087327</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>3.0</td>\n",
" <td>5.0</td>\n",
" <td>0.0</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>False</td>\n",
" <td>common</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.115198</td>\n",
" <td>0.167028</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>3.0</td>\n",
" <td>5.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>True</td>\n",
" <td>common</td>\n",
" <td>common</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.175993</td>\n",
" <td>0.239769</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>3.0</td>\n",
" <td>5.0</td>\n",
" <td>2.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>True</td>\n",
" <td>common</td>\n",
" <td>common</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>3</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.224014</td>\n",
" <td>0.306158</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>3.0</td>\n",
" <td>5.0</td>\n",
" <td>3.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>True</td>\n",
" <td>common</td>\n",
" <td>common</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>4</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.252047</td>\n",
" <td>0.366749</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>4.0</td>\n",
" <td>7.0</td>\n",
" <td>4.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>True</td>\n",
" <td>rare</td>\n",
" <td>common</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Unnamed: 0 Action1 Action2 Agent DLS reliability HPC reliability \\\n",
"0 0 0.0 0.0 0.0 0.057327 0.087327 \n",
"1 1 0.0 0.0 0.0 0.115198 0.167028 \n",
"2 2 0.0 0.0 0.0 0.175993 0.239769 \n",
"3 3 0.0 0.0 0.0 0.224014 0.306158 \n",
"4 4 0.0 0.0 0.0 0.252047 0.366749 \n",
"\n",
" P(SR) Reward StartState State2 Terminus Trial PreviousAction \\\n",
"0 1.0 1.0 0.0 3.0 5.0 0.0 NaN \n",
"1 1.0 1.0 0.0 3.0 5.0 1.0 0.0 \n",
"2 1.0 1.0 0.0 3.0 5.0 2.0 0.0 \n",
"3 1.0 0.0 0.0 3.0 5.0 3.0 0.0 \n",
"4 1.0 1.0 0.0 4.0 7.0 4.0 0.0 \n",
"\n",
" PreviousStart PreviousReward Stay Transition PreviousTransition \n",
"0 NaN NaN False common NaN \n",
"1 0.0 1.0 True common common \n",
"2 0.0 1.0 True common common \n",
"3 0.0 1.0 True common common \n",
"4 0.0 0.0 True rare common "
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_lesion_dls.head()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def compute_mean_stay_prob(data):\n",
" means = data[data['Trial']>0].groupby(['PreviousTransition', 'PreviousReward'])['Stay'].mean()\n",
" sems = data.groupby(['PreviousTransition', 'PreviousReward'])['Stay'].sem()\n",
" return means, sems"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"mean_lesion_hpc, sem_lesion_hpc = compute_mean_stay_prob(data_lesion_hpc)\n",
"mean_lesion_dls, sem_lesion_dls = compute_mean_stay_prob(data_lesion_dls)\n",
"mean_full, sem_full = compute_mean_stay_prob(data_control)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"PreviousTransition PreviousReward\n",
"common 0.0 0.558473\n",
" 1.0 0.836690\n",
"rare 0.0 0.563003\n",
" 1.0 0.823360\n",
"Name: Stay, dtype: float64"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mean_lesion_hpc"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"def plot_daw_style(ax, data, yerr=None, title=''):\n",
" lightgray = '#d1d1d1'\n",
" darkgray = '#929292'\n",
"\n",
" bar_width= 0.2\n",
"\n",
" bars1 = data[:2][::-1]\n",
" bars2 = data[2:][::-1]\n",
" if yerr is not None:\n",
" errs1 = yerr[:2][::-1]\n",
" errs2 = yerr[2:][::-1]\n",
" else:\n",
" errs1 = yerr\n",
" errs2 = yerr\n",
" \n",
" # The x position of bars\n",
" r1 = np.array([0.125, 0.625]) \n",
" r2 = [x + bar_width + .05 for x in r1]\n",
" list(sem_full),\n",
" plt.sca(ax)\n",
" \n",
" plt.bar(r1, bars1, width=bar_width, color='blue', capsize=4,yerr=errs1)\n",
" plt.bar(r2, bars2, width=bar_width, color='red', capsize=4,yerr=errs1)\n",
" plt.xticks([r+ bar_width/2 +.025 for r in r1], ['Rewarded', 'Unrewarded'], fontsize=12)\n",
" plt.yticks(fontsize=12)\n",
" plt.title(title, fontsize=16)\n",
" plt.ylim([0.4, 1])\n",
" plt.xlim([0, 1])"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAC8CAYAAACNKy/mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucVVX9//HXm4uBA4gkoiiCpXxNURHJawiYlzTF1PAGKt80U+xbVvZTy+AMWpaa9u2bl7wUXr6SaYoY3vICQkqKBn3FlEwBFS8gggyCCnx+f+w9w+FwzswenDkzw7yfj8d+zDlrr73XOvusc2adtddFEYGZmZlZS9WmqTNgZmZm9mm4MmNmZmYtmiszZmZm1qK5MmNmZmYtmiszZmZm1qK5MmNmZmYtmiszZSTpa5KelPSupJWS5kuaKOkreXGGSMpJyvzeSBolKST1qWd+uqZpDajPcdZy5JWNnYrsa5fuyxXE7VPmbJo1qLyyXGw7ZCPOV/M5SZ/nJDXreU0kTZE0ZSOOa5HfA67MlImk7wD3Av8CzgC+Clya7j44L+oQYCz1e28mA/sDb9UzW13TtFyZMdj4cmTWXA0nKdP52zNNmiNrFO2aOgOtyPnAxIg4Iy/sceDG+rTC5JPUHlgdEYuARQ2QR2vFXI5sEzQrIl5p6kxY43PLTPl0A94utiMi1kLSdEnSUgLwSXWzaLqvT/p8tKTLJS0EPgK6FmsWlHSSpMclLZJUJenvkk7P298HeC19emNeE+yodP88SeML81qquVXSLpIelrRC0gJJ/5nuP1XSS2kenpD0+XpfOSuLEuVonqTbJX1T0iuSVkl6XtLQgmPHS3pD0gGSnk3jzZP0X0XS2UfSo2mZWCHpMUn7FIk3WNJfJC1L482WdEbe/lrLeF68kHSppB+kt3ZXSJosaet0+2OaxuuSLihxTQ5KbwlXSXpP0jWSOubFG5LGG5Lhmp6S5rUqTff/JH2r1jfHGlyp2ykNeQsp7/NzqqSXlXQvmCZpZ0kVkn6blqd3JP1SUruC4/9D0r2SlqbHzlBet4S8eCel37MfSZoj6dgS+dlK0nWS3kzjviTprIZ4rU3NLTPl8wxwuqRXgfsiYm6RODcB25PchvoSsKZInB8DzwJnAW2BVSXS+xxwN/BzYC1wEHCTpI4RcT3JrYTjgHuAy4BJ6XH/rv9LA+Au4EbgSmA08DtJO5PcNrsQaA/8N3AHsO9GpmEbr23hFyVJ+cliMLA3Sdn7CLgAeFDSnhHxcl68LsCdwC+AV4CTgF9LWh4R4wEk7QFMBV4ERgFBUj6mStovIman8Y4B/gT8FfgWsBjYDeidl15dZTzfqcALJGWzB/Ar4FagM/AgcAPJLYmfS/q/iHig4PjbgT8C1wL7AGOAivQ1ZCbpS+m5fg38kOQH5S4kt3yt4RWW+4iIYt+rjekg4PMkn5vNSMren4BXWfc5OQi4mOT791oAST2B6cBy4NvAMuBcYLKkoyLiwTTeISTfq5OBHwDdSb5r2wM1n09JXUg+Tx2BHMmP2cOB6yR9JiL+p7EuQFlEhLcybEBf4B8kX95B8uU8ATisIF4u3d+uILxPGv48oIJ9o9J9fUqk3Yak4nojMLvIOc8scsw8YHyR8AByRfJ7Wl7YlsBq4D2gS174d9K4vZv6/WgtW17ZqG3LlSpHaTn4GNghL6wzsAS4LS9sfHrsSQXp/wWYX11mSSofS4GueXG6pOe7J32uNN2ZQJuMr7NoGc8rs3PzP1PAVWn4xXlh7YB3gd8XuX7XF5zzxyQ/Nvqmz4ek8YaUuP590ufnA0uaulxs6lst5X56qfcmLzxHUukpLEO52uKUyMe8tGxvkRdW/T14U0Hc54En8p5fSfI9ulNeWFuSCsrzeWF/Jflx0CYvbN80jSl5YT8h+fG7c0G6N5L8P2pX23Vp7ptvM5VJJC0xe5H8yv0pMAs4FnhY0sX1ONXESEtcbdJmzAmS3gQ+Sbczgf+od+azebD6QUS8T/JPYUZEfJAX56X0b69GyoOVdizwxYJtv4zHzoiIBdVPImI56zoL51tD8osz3x+AHYDt0ucHAX+OiKV55/uApGVwcBr0HyQtMDdFegu2mHqW8b9ExOq859Vl8eG8fKwm+aVcrHz+scjrakPSSlMfzwJbprcejpLkFpnGVVjuz6g9eqN4OiKW5T3foOzlheeXvYNIPns1fX4iaVWaAPSX1EVSW5LXdXf+ZyUi/kZSkcr3FeBvwGtKRjK2S1utHgY+C+y6sS+wOfBtpjJKC+KT6VbdjPgQMFbSNWkloC51jjSR1InkF/GHJE34/yb5dX0O8I2Ny32dCvP+cYkwgA6NlAcr7YUo6AhZ5LZTKe+UCNuuIOz9iPikxLHbAW+Q9B0rVobfJmnRg+SLlTR+URtRxkuVxWLhxcpn4TXIf12ZRcRUScOB/yIZ3YikqcD3I+If9TmXZbJBuW8CG1v2ugF/L3K+t0laL7ckuWXUntKf0XxbAzuRVPqL+WyJ8BbBlZkmFBELJd1Ecn9zZ7INGczSMW1/kl+2gyJienVgPf55QdIcuVl+gKRu9TjeNh09SoS9WRC2paT2BRWa6mOr4y4Btilyvm3SfZA0eUPtFYWGKOP10QOYU/Ac1r2u6r5r631mKPIPIiLuBu5OK2RDSPoYPSRp+9paoqzBZX7Pmkhtn5VI939IUjkp9Rmdn/f8PZIW8++WSO/lEuEtgm8zlYmkUrdWdkn/Vo90+ij927FI3Kw2T//W/FORtCVwTEG82tKaD/QrCDvqU+TJWq798suvpM4k8yQ9XRCvLXB8QdhJwALW/dOfCnw1PUf++Y5O90HSv2UecKYklchT1jLeUE4oeH4SSafj6h8g1f80Cj8zR5Y6YURURcSfgd8C29J8/om2Fhu8Z2ll+LCmyc4GppJ89vpUB6S3lU4E/h4Ry9PW/meBrytvig9J+5L0icz3EMn/mwURMbPItrxxX07jcstM+bwg6QmSpuXXSDo9HgmcDfwxr0/Ci+nfH0h6EFgTETPrmdZTwAfANZLGkoy6uJjkF+8WefHeIamtnyTpH8AK4LWIeI+kT8DvJF0N/BnYk3qO3LBNxjvAI0qG5FePZqoALimItxy4XNJWJJNDngwcAozK6+d1CUml+DFJvyD5hXkBSeVkHCS9KiWdRzLS7nFJ15PMf/MFYOuIGEv2Mt5QjpR0BfAIST+ZscCtaV84IuKt9HbRRZIWk/wCHkkyiqWGpHEkv5ifABaSjF78Dsl8KJ7jp7yeJbk9eUVaEfiIZLTbZ5o0V+tcTfKd+5e0jH9Akr++JD8mqo0lKZcTJf2WZDRTJRtOBXI1SUVoWvq9/jLJ52YXkhbOxvohUBZumSmfC0iu9ziSgncnSVP5hSTDRqv9mWRo3miSX77P1jeh9EvxWJJfyneTDL2+iWRIaH68tSQdJrcEHk3TOjrdfQvJh+Q44H6SIXxF5y6wTd5U4JfAz0jKbQfgiNhweoEPSFosTgfuA4YC342IW6ojpP1ChqRxbwFuA6qAwZEOy07j3Qccmj69maSD8FmknRqzlvEGNJLkn8i9JMNfbyT5jBbGmUEy7Ho8SYvUpQVx/kbyi/lqkj4/vyBtrWqcbFspaYfvY4DXSd6va0jek/FNl6t1ImIhyRQdc4DrSMp5N+CrEfFQXrxHgREkHd/vIRnyfx4Ft43STsgHAA+Q/D96GPgdyTV4opFfTqNT1D0wxsxaKUnzSIazjqwj3njgkIjYvhz5Khclk0j+nmQ4a1N3JDWzEtwyY2ZmZi1a2Sozkr4taWY6hfL4OuJ+T9LbSqb6/p2k5nIP08zMzJqZst1mknQcSe//w4GOETGqRLzDSaYZP5ikg9y9JBMHXViWjJqZmVmLUraWmYi4JyImkoyeqc3pwM0RMSedRO4SPIrGzMzMSmiOfWZ2A2bnPZ8N9JDkORjMzMxsA81xnplOJKuDVqt+3JmCVp106fKzACoqKvbeZZddMGsIzz333OKI6N7Y6bgMW2MpVxkGl2NrHPUpw82xMlNFMqFcterHG8xOGBE3ADcADBw4MGbOrO/ccmbFSZpfd6xPz2XYGku5yjC4HFvjqE8Zbo63meaQzDZbbU/gnXRWWjMzM7P1lHNodjtJHUhm7GwrqUOJReFuBc6QtGu61srFNJMZGc3MzKz5KWfLzMXASpLp+0emjy+WtIOkKkk7AKTTNF9OMr3y/HQbW8Z8mpmZWQtStj4zEZEDciV2dyqIexVwVSNnyczMzDYBzbHPjJmZmVlmrsyYmZlZi+bKjJmZmbVorsyYmZlZi+bKjJmZmbVorsyY2UZr27Yt/fv3p1+/fgwfPpwPP/ywQc575JFHsnTp0gY5V6FZs2bxwAMPNMq5811//fXceuutAIwfP56FCxfW7DvzzDN58cUXGz0PZq2FKzNmttE6duzIrFmzeOGFF9hss824/vrr19sfEaxdu7be533ggQfo2rVrQ2VzPbVVZlavXt1g6Zx99tmcdtppwIaVmZtuuoldd921wdIya+1cmTGzBjFo0CBeeeUV5s2bxxe+8AVGjx7NgAEDeP3113nkkUfYf//9GTBgAMOHD6eqqooHH3yQE044oeb4KVOmcPTRRwPQp08fFi9eDMBVV11Fv3796NevH7/61a8AmDdvHv369as59sorrySXywHw61//ml133ZU99tiDk046ab08fvzxx4wZM4Y777yT/v37c+edd5LL5TjrrLM47LDDOO2005g3bx6DBg1iwIABDBgwgKeeeqomf0OGDOHrX/86u+yyCyNGjCAiALjwwgtr0jz//PMByOVyXHnlldx9993MnDmTESNG0L9/f1auXMmQIUOoXr9owoQJ7L777vTr148LLrigJq+dOnXixz/+MXvuuSf77bcf77zzToO9V2abnIjYJLa99947zBoKMDNchutUUVERERGffPJJDBs2LK699tp47bXXQlI8/fTTERGxaNGiGDRoUFRVVUVExM9//vOorKyMTz75JHr16lUTfvbZZ8dtt90WERG9e/eORYsWxcyZM6Nfv35RVVUVy5cvj1133TWef/75eO2112K33XaryccVV1wRY8eOjYiIbbfdNlatWhUREe+///4Gef79738f5557bs3zsWPHxoABA+LDDz+MiIgVK1bEypUrIyJi7ty5Uf2+PPHEE9GlS5d4/fXXY82aNbHffvvFtGnT4r333ou+ffvG2rVr10tz7NixccUVV0RExODBg+PZZ5+tSbP6+Ztvvhm9evWKd999Nz755JMYOnRo3HvvvRERAcSkSZMiIuKHP/xhXHLJJfV6b5qiDEcLLcfWPNWnDLtlxsw22sqVK+nfvz8DBw5khx124IwzzgCgd+/e7LfffgDMmDGDF198kQMPPJD+/ftzyy23MH/+fNq1a8dXvvIV7r//flavXs3kyZM55phj1jv/9OnTOfbYY6moqKBTp04cd9xxTJs2rdY87bHHHowYMYLbb7+ddu2yTXI+bNgwOnbsCMAnn3zCN7/5TXbffXeGDx++Xt+WffbZh+233542bdrQv39/5s2bR5cuXejQoQNnnnkm99xzD5tvvnnm6/fss88yZMgQunfvTrt27RgxYgRPPvkkAJttthlHHXUUAHvvvTfz5s3LfF6z1qZsyxmY2aanus9MoYqKiprHEcGhhx7KhAkTNoh34okncs0119CtWze++MUv0rlz5/X2R3obp1C7du3W64uzatWqmseTJ0/mySefZNKkSVxyySXMmTOnzkpNfn6vvvpqevTowezZs1m7di0dOnSo2feZz3ym5nHbtm1ZvXo17dq145lnnuGxxx7jD3/4A7/5zW94/PHHa02vrtcH0L59eyStl5aZFeeWGTNrVPvttx9//etfeeWVVwD48MMPmTt3LgBDhgzh+eef58Ybb+TEE0/c4NiDDjqIiRMn8uGHH7JixQruvfdeBg0aRI8ePXj33Xd57733+Oijj/jzn/8MwNq1a3n99dcZOnQol19+OUuXLqWqqmq9c3bu3Jnly5eXzO+yZcvYdtttadOmDbfddhtr1qyp9fVVVVWxbNkyjjzySH71q18VrdyVSnPfffdl6tSpLF68mDVr1jBhwgQGDx5ca3pmtiG3zJhZo+revTvjx4/n5JNP5qOPPgLg0ksvpW/fvrRt25ajjjqK8ePHc8stt2xw7IABAxg1ahT77LMPkAxp3muvvQAYM2YM++67LzvuuCO77LILAGvWrGHkyJEsW7aMiOB73/veBqOihg4dys9//nP69+/PRRddtEGao0eP5vjjj+euu+5i6NCh67XaFLN8+XKOOeYYVq1aRURw9dVXbxBn1KhRnH322XTs2JGnn366JnzbbbflsssuY+jQoUQERx555Aa32sysbqqtmbNBE5K6ATcDhwGLgYsi4o4i8boC/w0ckQZdG8mK27UaOHBgVI8OMPu0JD0XEQPLmabLsDWkpijD4HJsDac+ZbicLTPXAB8DPYD+wGRJsyNiTkG8q4HNgT7A1sBjkuZHxO/LmFczMzNrIcrSZ0ZSBXA88JOIqIqI6cAk4NQi0Y8GLo+IDyNiHklrzjfKkU8zMzNreTJVZiT9XdJ5knpsZDp9gTURMTcvbDawW6kkCx73KxHPzMzMWrmsLTOXAgcBr0p6UNIpkjrWI51OwLKCsGVA5yJxHwIulNRZ0k4krTJFJ26QdJakmZJmLlq0qB7ZMWseXIZtU+BybE0tU2UmIv4UEccBvYD7gNHAW5J+J+ngDKeoAroUhHUBio2P/A6wEvhXmtYE4I0S+bohIgZGxMDu3btneSlmzYrLsG0KXI6tqdWrz0xELAFuBa4HFpD0g7lB0lxJh9Ry6FygnaSd88L2BAo7/xIRSyJiRERsExG7pXl8pj75NDMzs9Yja5+ZNpIOl3Q7sBAYAfwc2CYidgIuAm4vdXxErADuAcZJqpB0IHAMcFuRtD4v6bOS2ko6AjiL5DaXmZmZ2QaytswsBH4J/APYNSKOiIg7ImIlJLehgH/WcY7RQEfgXZJbR+dExBxJgyTlT9G5N/B/JLegLgNGFBm+bWZmZgZkr8wcFRH9IuLyiFhYLEJEDK3tBOnto69FREVE7FA9YV5ETIuITnnx/hgRPSNi84joHxEPZ341ZtYsSbVv9XXHHXcwcOBAOnXqxLbbbssRRxzB9OnTGz7jZtYiZK3MPFIsUNK7DZgXM7M6XXXVVZx33nn86Ec/4p133mHBggWMHj2a++67r6mzZmZNJGtlpn1hgKT2QNuGzY6ZWWnLli1jzJgxXHPNNRx33HFUVFTQvn17jj76aK644go++ugjzjvvPHr27EnPnj0577zzataDmjJlCttvvz2XX345W2+9Ndtuuy0TJ07kgQceoG/fvnTr1o2f/exnNWnlcjmGDx/OyJEj6dy5M7vvvjtz587lsssuY+utt6ZXr1488si633kLFy5k2LBhdOvWjZ122okbb7xxvXOdcMIJnHbaaXTu3JnddtsNT/lv1nBqrcxImibpSaCDpCfzN+Bl4Kmy5NLMDHj66adZtWoVxx57bNH9P/3pT5kxYwazZs1i9uzZPPPMM1x66brxA2+//TarVq3izTffZNy4cXzzm9/k9ttv57nnnmPatGmMGzeOV199tSb+/fffz6mnnsr777/PXnvtxeGHH87atWt58803GTNmDN/61rdq4p588slsv/32LFy4kLvvvpsf/ehHPPbYYzX7J02axEknncTSpUsZNmwY3/72txvhCpm1UhFRcgNOB0aRzPtyet52GnA40L6248u57b333mHWUICZ4TLcYKD2Lavbb789evToUXL/5z73uZg8eXLN84ceeih69+4dERFPPPFEdOjQIVavXh0RER988EEAMWPGjJr4AwYMiHvvvTciIsaOHRuHHHJIzb5JkyZFRUXFBse///77sWDBgmjTpk188MEHNfEvvPDCOP3002vO9eUvf7lm35w5c6JDhw7ZX/hGaIoyHJt4Obbyqk8ZrnWhyYi4BUDSjIh4qXGqU2Zm2Xz2s59l8eLFrF69mnbtNvz6WrhwIb1796553rt3bxYuXLje8W3bJnfHO3ZMJjHv0WPdKi0dO3akqmrd4MrCfVtttdUGx1dVVbFw4UK6detG587rJjXv3bv3ereSttlmm5rHm2++OatWrSr5OsysfkreZpKUvwjkAZK+UWwrQx7NzADYf//96dChAxMnTiy6v2fPnsyfP7/m+YIFC+jZs2ej56tnz54sWbKE5cvXTWq+YMECtttuu0ZP28yotWXmZNZNaldsdWuAAH7XoDkyMythiy22YNy4cZx77rm0a9eOww47jPbt2/Poo4/yxBNPcPLJJ3PppZfyxS9+EUmMGzeOkSNHNnq+evXqxQEHHMBFF13ElVdeydy5c7n55pu5/faSc4maWQMqWZmJiCPzHtc6h0xLl8vlqKysLLl/7Nix5HK58mXIzEr6/ve/T48ePbj00ksZMWIEnTt3Zu+99+bHP/4xAwYM4IMPPmCPPfYAYPjw4Vx88cVlydeECRM4++yz6dmzJ1tuuSWVlZUceuihZUnbrLVT0semyA4p6yKUaxs0Rxtp4MCBkX9/emMm4koMSf9OqfeRJS6ltUCSnouIgeVMs7AMm30aTVGGweXYGk59ynBtFZbVwCe1bNX7NwE5QOk2Nd2Ut+WaKmNmZq1CLpdDUsnNreNWm9oqMzsCn6tlq96/CciRdP8pteXKnyN/sM2shatrGYv8rbIyx7rv3MHptu57uLIyl/lc1vrU1mdmfql9tpHq8SnLsa4KNST9OyU/QmVlstXF977MrEXIAYXfafnfmWNxK7mVUtvQ7BvyHt8m6dZiW9aEJHWTdK+kFZLmSzqlRLzPSLpe0juSlki6X1KrG9+Ywze+zKw1ydHcWsit5ahtaPZreY9faYC0rgE+BnoA/YHJkmZHxJyCeN8F9gf2AJYBNwL/AxzXAHloMXL4o2tmZpZFbbeZLst7nOF+RmmSKoDjgX4RUQVMlzSJZP6aCwui7wg8HBHvpMf+Abjq06RvZmZmm66sq2Yj6WBJN0qanP79cj3S6QusiYi5eWGzgd2KxL0ZOFBST0mbAyOAB+uRllmr4w7jZtaaZarMSPo+8AdgCTAZeA+4Q9IPMqbTieSWUb5lQOcicecCC4A3gQ+ALwDjSuTrLEkzJc1ctGhRxqyYNR+1lWGPBLGWwt/F1tSytsz8ADg4Ii6IiGsj4kLg4DQ8iyqgS0FYF2B5kbjXAR2AzwIVwD2UaJmJiBsiYmBEDOzevXvGrJg1Hw1XhnO4y7g1FX8XW1PLfJuJDTsBv0rysy+LuUA7STvnhe0JFHb+rQ4fHxFLIuIjks6/+0jaqh55NWtlcjTrkSAN2BzUp08fOnbsSKdOndhmm20YNWrUeitdm1nrU9vQ7DbVG8k34c2SdpbUUVJf4AaSgf91iogVJC0s4yRVSDoQOIZ1C1nmexY4TdIWktoDo4GFEbG4Xq/MzDZZ999/P1VVVcyaNYu///3vXHbZZXUfVGD16tWNkDMzawpZlzP4Lckq2i+T3DL6J0nH3N/WI63RQEfgXWACcE5EzJE0SFL+z6rzgVXAv4BFwJHAsfVIx8xaiW222YbDDz+cWbNmATB58mT22msvunTpQq9evdbr+Dxv3jwkcfPNN7PDDjtw8MEHAzBjxgwOOOAAunbtyp577smUKVOa4JWY2aeRdTmDHfO2jVrOIL1t9LWIqIiIHSLijjR8WkR0yov3XkSMiIitI6JrRHwpIp6p/0szs03dG2+8wYMPPshOO+0EQEVFBbfeeitLly5l8uTJXHfddUycOHG9Y6ZOnco///lPHn74Yd58802++tWvcvHFF7NkyRKuvPJKjj/+eNyJ1cCjBFuSkqtmtzQNt2r2xqvzUjbLTFkxzWHV7KYaWdQoRaauF1OPRPv06cPixYuRRFVVFQcffDB/+tOf6Nq16wZxzzvvPCRx9dVXM2/ePHbccUf+/e9/87nPJb/DfvGLX/DCCy9w223r7ngffvjhnHLKKZx++umZ89QcNYdVs5tlGd7ITA1J/07ZqKPxd/FGaKhVswtPOkzSLyXdsjHLGZiZNZSJEyeyfPlypkyZwksvvcTixUmXur/97W8MHTqU7t27s8UWW3D99dfX7KvWq1evmsfz58/nrrvuomvXrjXb9OnTeeutt8r6eqx5yuExgi1F1nlmxpL0j2kDDCeZZ+ZwYGnjZc3MrHaDBw9m1KhRnH/++QCccsopDBs2jNdff51ly5Zx9tlnU9j6rLxf5r169eLUU09l6dKlNduKFSu48MLCicmtNcrRrMcIWp6sLTPfAA6NiO8BH6d/jwb6NFbGzMyyOO+88/jLX/7CrFmzWL58Od26daNDhw4888wz3HHHHbUeO3LkSO6//34efvhh1qxZw6pVq5gyZQpvvPFGmXJvZg0ha2Wma0S8kD7+WFL7tFPu4EbKl5lZJt27d+e0007jkksu4dprr2XMmDF07tyZcePGccIJJ9R6bK9evbjvvvv42c9+Rvfu3enVqxdXXHEFa9euLVPuzawh1LZqdr5/S9otXeH6BeAcSe8D7zde1sxsk9GAnR/nzZu3Qdh1111X8/jrX/960eP69OmzwS0ngH333ZepU6c2WP7MrPyyVmYuJlleAOAi4H9J1lsa3RiZMjMzM8sqU2UmIh7Ie/w3YKdGy5GZmZlZPWRtmSFdV+kEoCewEPhjRPyrsTJmZmZmlkXWodmnAH8H9gBWALsDz6fhZmZmZk0ma8vMpcCREfFkdYCkQSQLRdY+9tHMzMysEWUdmt0ZeLogbAZQ0bDZMTNrXbz+j9mnl7UycxXwM0kdACR1BH6ahpuZWT4p85arrKyZUXZwuq03y2xlZbZzmbViJW8zSXqd5LMEyTIU2wDfTeeX2TINewu4LEtCkroBNwOHAYuBi6pXzi6I9yAwKC9oM+DliNg9SzpmZmbWutTWZ2ZkA6d1DfAx0APoD0yWNDudiK9GRByR/1zSFODxBs6LmVmzkAMqC8Ly21nG4jWAzOpSsjITEQ02JaakCuB4oF9EVAHTJU0CTgVKrugmqQ9JK83YEd1bAAAO2ElEQVR/NlRezMyakxyurJh9WlmHZreXVCnpVUmr0r+VkjbLmE5fYE1EzM0Lmw3sVsdxpwHTIuK1jOmYmZlZK5O1A/DlwCHA2cCe6d+DgV9kPL4TsKwgbBnJKKnanAaML7VT0lmSZkqauWjRooxZMWs+XIZtU+BybE0ta2VmODAsIh6JiJcj4hHgWJIZgbOoAroUhHUBlpc6QNKXSDod310qTkTcEBEDI2Jg9+7dM2bFrPlwGbZNgcuxNbWslZlS4/6yjgecC7RLl0Soticwp0R8gNOBe9I+NmZmZmZFZa3M3AXcL+lwSV+Q9BVgIvDHLAdHxArgHmCcpApJBwLHkMwgvIF0Hpvh1HKLyczMzAyyV2b+H/AoyfDq54D/AZ4ALqhHWqOBjsC7wATgnIiYI2mQpMLWl6+R9Kl5oh7nNzMzs1aozrWZJLUlmXPmZxExZmMTioglJJWUwvBpJB2E88MmkFR4zMzMzGpVZ8tMRKwBroqIVWXIj5mZmVm9ZL3NdL+koxs1J2ZmZmYboc7bTKkOwN2Sngby12wiIk5rjIyZmZmZZZG1MvNCupmZmZk1K5kqMxFRuA6amZmZWbOQtWUGSQcDJwM9gYXAHyLiscbKmJmZmVkWWRea/D7wB2AJMBl4D7hD0g8aMW9mZmZmdcraMvMD4OCIqOk3I+k24C/ALxsjY2ZmZmZZZB2aDfBKwfNXyRvVZGZmZtYUslZmcsDNknaW1FFSX+AGYKykNtVbo+XSzMzMrISst5l+m/49maQ1pnq17BHpPqXhbRs0d2ZmZmZ1yFqZ2bFRc2FmZma2kbLOMzO/sTNiZmZmtjHK1s9FUjdJ90paIWm+pFNqiTtA0pOSqiS9I+m75cqnmZmZtSyZJ81rANcAHwM9gP7AZEmzI2JOfiRJWwEPAd8D7gY2A7YvYz7NzMysBSlLy4ykCuB44CcRURUR04FJwKlFon8feDgi/jciPoqI5RHxz3Lk08zMzOonl8shqeSWy+UaPQ9ZZwAeJunTtOL0BdZExNy8sNnAbkXi7gcskfSUpHcl3S9phxL5OkvSTEkzFy1a9CmyZ9Y0XIZtU+ByvAmSMm+5ykqCZEjz4HSLvC1XWZn9fBspa8vMJcBbkn4jad+NSKcTsKwgbBnQuUjc7YHTge8COwCvAROKnTQiboiIgRExsHv37huRLbOm5TJsmwKX49YtRzI/i4Cp6aa8LVeGPGSqzETEnsAhwErgT5JelnSxpD4Z06kCuhSEdQGWF4m7Erg3Ip6NiFVAJXCApC0ypmVmZmZlkmP9lpjCLVeGPGTuMxMRsyPih0Av4FxgOPDvdNTRiDpmAJ4LtJO0c17YnsCcInH/wfrLJFQ/3vj2JzMzM9tk1asDsKTPA2OA64AO6eMbgW+TjDwqKiJWAPcA4yRVSDoQOAa4rUj03wPHSuovqT3wE2B6RCytT17NzMysdcjUqVfSuSQjj3YC/gicGhEz8vb/CXi3jtOMBn6XxnsPOCci5kgaBDwYEZ0AIuJxST8CJgObA9OBknPSmJmZWeuWdYTSEcAvgfsi4uPCnRHxoaTjajtBRCwBvlYkfBpJB+H8sOtIWn/MzMzMapV1OYOjMsR55NNnx8zMzKx+MveZSeea+aWkWyTdWr01ZubMatMcJmoyM7Oml3XSvLHAb9P4w0n6vBwOuFNuK1K2ykMLmqjJzMyaXtaWmW8Ah0bE94CP079HA30aK2NWHvWY5JHKyhzUUn2orMyVtd6Qo+knajIzs6aXtTLTNSJeSB9/LKl9RDxD8t/MWo0czan6kKPpJ2oyM7Oml3U0078l7ZaucP0CcI6k94H3Gy9r1vzkcBXBzMyam6yVmYuBz6aPLwTuIBlOfW5jZMrMzMwsq6xDsx/Ie/wMyeR5ZmZmZk0u62imJSXC65r118zMzKxRZe0A3L4wIF03qW3DZsfMzMysfmq9zSRpGsnAkA6SnizYvT3wVGNlzMzMzCyLuvrM3EQy5vaLwM154QG8AzzeSPkyMzMzy6TWykxE3AIgaUZEvPRpEpLUjaRCdBiwGLgoIu4oEi8H/Bj4KC94j4h49dOkb2ZmZpumWvvMSNpbUr/qioyk7pL+V9JsSddL6lTb8QWuAT4GegAjgOsk7VYi7p0R0Slvc0XGzMzMiqqrA/CvgG3ynt8E9AVuAPoBl2dJRFIFcDzwk4ioiojpwCTg1Hrn2MzMrJXyArvF1VWZ+QIwDUBSV+AIYEREXAOcTLI+UxZ9gTURMTcvbDZQqmXmaElLJM2RdE7GNMzMzFqcplgjb1NbX7euykw7kltDAPsBb1dXSCLidaBrxnQ6AcsKwpYBnYvE/SNJJao78E1gjKSTi51U0lmSZkqauWjRooxZMWs+XIZtU+BybE2trsrMHGB4+vgk4NHqHZK2Y8MKSilVQJeCsC7A8sKIEfFiRCyMiDUR8RTw38DXi500Im6IiIERMbB79+4Zs2LWfLgM26bA5biccjSnBX+bi7qGZl8A3C/pemAN8KW8fScCf82YzlygnaSdI+JfadieJJWlugTJO2RmZtbK5WitFZba1Noyk3bU3QE4FPhcRLyct3sy8L0siUTECuAeYJykCkkHAscAtxXGlXSMpC2V2Af4DnBfpldjZmZmrU6dyxlExPKIeC4ilheEvxwRC+uR1migI/AuMAE4JyLmSBokqSov3knAKyS3oG4FflE9342ZmZlZoUyrZjeEiFgCfK1I+DSSDsLVz4t29jUzMzMrJutCk2ZmZmbNkiszZmZm1qK5MmNmZmYtmiszZmZm1qK5MmNmZmYtmiszZtYovCCemZVL2YZmm9kmoB6r0+VYN0/pkPTvlPwIlZXJlkVE5nTNrPVxy4yZNYocXkHGzMrDLTNm1ihyuMJiZuXhlhkzMzNr0VyZMTMzsxbNlRkzMzNr0VyZMbNWw8PFzTZNZavMSOom6V5JKyTNl3RKHfE3k/SSpDfKlUcza3mk7FtlZQ6IdBucblGzVVbmMp/LzJqPcrbMXAN8DPQARgDXSdqtlvg/BN4tR8bMrLXI4QHjZpueslRmJFUAxwM/iYiqiJgOTAJOLRF/R2AkcFk58mdmrUWO/JaYDbdcU2XMzD6FcrXM9AXWRMTcvLDZQKmWmf8BfgSsbOyMmZmZWctWrknzOgHLCsKWAZ0LI0o6FmgXEfdKGlLbSSWdBZyVPq2S9HID5HWjFbmPvhWwuPw5ydMMb+63kOvUuzzJNq8yDBtciqZ/b6AllOOmv05NVIaTpJtXOW4h3zFNrtmVYSjMVOYyrCjDmieS9gL+GhGb54X9ABgSEUfnhVUAs4AjI+JfaWXm9ojYvtEz2QgkzYyIgU2dj+bO16n58nuTja9T8+b3p24t/RqVq2VmLtBO0s4R8a80bE9gTkG8nYE+wDQltbPNgC0kvQ3sFxHzypNdMzMzaynKUpmJiBWS7gHGSToT6A8cAxxQEPUFoFfe8wOA3wADgEXlyKuZmZm1LOUcmj0a6Egy3HoCcE5EzJE0SFIVQESsjoi3qzdgCbA2fb6mjHltKDc0dQZaCF+n5svvTTa+Ts2b35+6tehrVJY+M2ZmZmaNxcsZmJmZWYvmykwzImm8pEvLfeymxNewafn6f3q+hk3L179hlPs6tprKjKR5klZKqpL0dnqxOjV1vpqCpJC0U0FYTtLtTZUnq5vL8Douwy2Ty/D6XI4bTqupzKSOjohOJKOp9gIuaopMSCrXkPiykNS2CdLcpK5hPbgMNwKX4bJyGW4krbkct7bKDADpSKmHST5MSPqMpCslLZD0jqTrJXVM902VdHz6+EtpTfrI9Pkhkmaljz8v6XFJ70laLOl/JXWtTjP9RXKBpH8AKyS1k7SXpOclLZd0J9AhP5+SjpI0S9JSSU9J2iNvX63HfhqShkh6Q9IPJL0r6S1J/5m3f7yk6yQ9IGkFMNTXsLxchmvnMtz8uQzXzeU4u1ZZmZG0PXAE8Eoa9AuS9aP6AzsB2wFj0n1TgSHp44OAV4HBec+nVp+WZGHMnsAXSObLyRUkfTLwVaArybWfCNwGdAPuIlmMszqPA4DfAd8CPgv8FpiUFtTNaju2gWwDbEFyLc4ArpG0Zd7+U4CfkixJMR1fw7JyGc7EZbgZcxnOzOU4i4hoFRswD6gClpMsj/tY+iYIWAF8Pi/u/sBr6eMvA/9IHz8EnAnMSJ9PBY4rkd7XgL8XpP+NvOcHAQtJh8enYU8Bl6aPrwMuKTjnyyQFr9ZjM1yLAHYqCMuRLB0BSWFfSbJGVvX+d0lmYQYYD9yat6/VXUOXYZfhln4NXYab/vq7HDdcOW4W97rK6GsR8aikwcAdJAtrbQZsDjyndQtcCai+9/g00FdSD5Ka7jCgUtJWwD7AkwCStgZ+DQwiqSG3Ad4vSP/1vMc9gTcjfedS8/Me9wZOl/RfeWGbpcdFHcfWZQ3QviCsPfBJ3vP3ImJ13vMPSRYMrZb/WrrT+q5hU3EZTrgMuwxvCtff5biBynGrvM0UEVNJarRXkqwSuhLYLSK6ptsWkXRQIyI+BJ4Dvgu8EBEfk9Qavw/8OyKqVxm9jORN2SMiugAjSQrReknnPX4L2E5ab4nQHfIevw78NC9PXSNi84iYkOHYuiwgWQMr347UrwDlv5bWeA2blMuwy3CGY5s1l2HA5bjBynGrrMykfgUcCuwB3AhcndZGkbSdpMPz4k4Fvs26+4lTCp5DUnutApZK2g74YR3pPw2sBr6Tdp46jqRmXO1G4GxJ+ypRIemrkjpnOLYudwIXS9peUhtJhwBHA3fX4xw1ImItre8aNgcuwy7DLsMt+/q7HDdQOW61lZmIWATcCvwEuICkE9oMSR8AjwL/kRd9Ksmb+2SJ5wCVJAtiLgMmA/fUkf7HwHHAKJLmuxPzj4mImcA3SRbafD/N36gsx2YwjqQ2Pj09/nJgRES8UI9zFGpt17DJuQy7DNd2bEvQysswuBw3WDn22kxmZmbWorXalhkzMzPbNLgyY2ZmZi2aKzNmZmbWorkyY2ZmZi2aKzNmZmbWorkyY2ZmZi2aKzNmZmbWorkyY2ZmZi2aKzNmZmbWov1/4AraKI7bcy8AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 648x180 with 3 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(1,3, figsize= (9.5,2.5), sharey=True)\n",
"\n",
"plot_daw_style(axes[0], list(mean_lesion_hpc), yerr=sem_lesion_hpc, title='Striatum')\n",
"plot_daw_style(axes[1], list(mean_lesion_dls), yerr=sem_lesion_dls, title='Hippocampus')\n",
"plot_daw_style(axes[2], list(mean_full), yerr=sem_full, title='Full model')\n",
"\n",
"\n",
"leg = axes[1].legend(['Common', 'Rare'], fontsize=12, frameon=False, handlelength=0.7, title='Previous transition')\n",
"plt.sca(axes[0])\n",
"plt.ylabel('Stay probability', fontsize=12)\n",
"\n",
"plt.savefig(os.path.join(figure_location, 'StayProbability.pdf'))\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"plt.close()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Doll et al analysis"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"figure_location = os.path.join(FIGURE_FOLDER,'twostep_deterministic')\n",
"if not os.path.exists(figure_location):\n",
" os.makedirs(figure_location)\n",
"\n",
"# load data\n",
"data_lesion_hpc = pd.read_csv(os.path.join(RESULTS_FOLDER, 'twostep_deterministic', 'results_lesion_hpc.csv'))\n",
"data_lesion_dls = pd.read_csv(os.path.join(RESULTS_FOLDER, 'twostep_deterministic', 'results_lesion_dls.csv'))\n",
"data_control = pd.read_csv(os.path.join(RESULTS_FOLDER, 'twostep_deterministic', 'results_control.csv'))\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Unnamed: 0</th>\n",
" <th>Action1</th>\n",
" <th>Action2</th>\n",
" <th>Agent</th>\n",
" <th>DLS reliability</th>\n",
" <th>HPC reliability</th>\n",
" <th>P(SR)</th>\n",
" <th>Reward</th>\n",
" <th>StartState</th>\n",
" <th>State2</th>\n",
" <th>Terminus</th>\n",
" <th>Trial</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.000000</td>\n",
" <td>0.030000</td>\n",
" <td>0.681250</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>2.0</td>\n",
" <td>4.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.000000</td>\n",
" <td>0.059100</td>\n",
" <td>0.671966</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>3.0</td>\n",
" <td>7.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.021000</td>\n",
" <td>0.087327</td>\n",
" <td>0.675625</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>3.0</td>\n",
" <td>7.0</td>\n",
" <td>2.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>3</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.029370</td>\n",
" <td>0.114707</td>\n",
" <td>0.670047</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>2.0</td>\n",
" <td>4.0</td>\n",
" <td>3.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>4</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.055789</td>\n",
" <td>0.141266</td>\n",
" <td>0.669577</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>3.0</td>\n",
" <td>7.0</td>\n",
" <td>4.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Unnamed: 0 Action1 Action2 Agent DLS reliability HPC reliability \\\n",
"0 0 0.0 0.0 0.0 0.000000 0.030000 \n",
"1 1 1.0 1.0 0.0 0.000000 0.059100 \n",
"2 2 1.0 1.0 0.0 0.021000 0.087327 \n",
"3 3 0.0 0.0 0.0 0.029370 0.114707 \n",
"4 4 1.0 1.0 0.0 0.055789 0.141266 \n",
"\n",
" P(SR) Reward StartState State2 Terminus Trial \n",
"0 0.681250 1.0 1.0 2.0 4.0 0.0 \n",
"1 0.671966 1.0 0.0 3.0 7.0 1.0 \n",
"2 0.675625 1.0 1.0 3.0 7.0 2.0 \n",
"3 0.670047 0.0 1.0 2.0 4.0 3.0 \n",
"4 0.669577 1.0 0.0 3.0 7.0 4.0 "
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_control.head()"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"def add_relevant_columns(dataframe):\n",
" dataframe['PreviousAction'] = dataframe.groupby(['Agent'])['Action1'].shift(1)\n",
" dataframe['PreviousStart'] = dataframe.groupby(['Agent'])['StartState'].shift(1)\n",
" dataframe['PreviousReward'] = dataframe.groupby(['Agent'])['Reward'].shift(1)\n",
" dataframe['Stay'] = (dataframe.PreviousAction == dataframe.Action1)\n",
" dataframe['SameStart'] = (dataframe.StartState == dataframe.PreviousStart)\n",
"\n",
"\n",
"add_relevant_columns(data_control)\n",
"add_relevant_columns(data_lesion_dls)\n",
"add_relevant_columns(data_lesion_hpc)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def compute_mean_stay_prob(data):\n",
" means = data[data['Trial']>0].groupby(['PreviousReward', 'SameStart'])['Stay'].mean()\n",
" sems = data.groupby(['PreviousReward', 'SameStart'])['Stay'].sem()\n",
" return means[::-1], sems[::-1]\n"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"mean_lesion_hpc, sem_lesion_hpc = compute_mean_stay_prob(data_lesion_hpc)\n",
"mean_lesion_dls, sem_lesion_dls = compute_mean_stay_prob(data_lesion_dls)\n",
"mean_full, sem_full = compute_mean_stay_prob(data_control)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"PreviousReward SameStart\n",
"1.0 True 0.944564\n",
" False 0.953787\n",
"0.0 True 0.484103\n",
" False 0.484945\n",
"Name: Stay, dtype: float64"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mean_lesion_dls"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"def plot_doll_style(ax, data, yerr=None, title=''):\n",
" lightgray = '#d1d1d1'\n",
" darkgray = '#929292'\n",
"\n",
" bar_width = 0.2\n",
"\n",
" bars1 = data[:2]\n",
" bars2 = data[2:]\n",
" if yerr is not None:\n",
" errs1 = yerr[:2]\n",
" errs2 = yerr[2:]\n",
" else:\n",
" errs1 = None\n",
" errs2 = None \n",
" \n",
" # The x position of bars\n",
" r1 = np.arange(len(bars1)) * .8 + 1.5 * bar_width\n",
" r2 = [x + bar_width for x in r1]\n",
"\n",
" plt.sca(ax)\n",
"\n",
" handle1 = plt.bar(r1, bars1, width=bar_width, color=lightgray, yerr=errs1, capsize=4)\n",
" handle2 = plt.bar(r2, bars2, width=bar_width, color=darkgray, yerr=errs2, capsize=4)\n",
" plt.ylabel('Stay probability', fontsize=15)\n",
" plt.xticks([r + bar_width / 2 for r in r1], ['same', 'different'], fontsize=15)\n",
" plt.yticks(fontsize=15)\n",
" plt.title(title, fontsize=18)\n",
" plt.ylim([0.45, 1.])\n",
" plt.xlim([0, 1.6])\n",
"\n",
" ax.spines['right'].set_visible(False)\n",
" ax.spines['top'].set_visible(False)\n",
" return handle1, handle2"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/matplotlib/tight_layout.py:176: UserWarning: Tight layout not applied. The left and right margins cannot be made large enough to accommodate all axes decorations. \n",
" warnings.warn('Tight layout not applied. The left and right margins '\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAC2CAYAAADEHHvuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XecVNX9//HXh7qUBekgRVBAQKQINoyNYCVGBBVjj6DGHk3QGDGW71djiZKIWAgaIbFEjdhiiYqiye+rERRQBBEFkSZtaVKX/fz+OHdwHGZ3Z3dn2Zmd9/PxuI/dOffce87Mfs7dM/eee665OyIiIiIiEtSo6gqIiIiIiGQSdZBFREREROKogywiIiIiEkcdZBERERGROOogi4iIiIjEUQdZRERERCSOOshZzMzONzM3s6Oqui5S/ZjZO2a2sKrrIZLNzOxvZlaYkPa/0bG7XVXVqzzM7N9mNr8C2+/yWYhkKnWQK5GZ7W1m481srpltMrMCM/vMzCaa2dFx+W42syFVUL+OUdl9dnfZsnuZ2VHRP+Rfl5DHzezl3VkvkUwR10aKWw6p6jqKyO5Tq6orUF2ZWX9gKrAdmATMBuoBXYGTgA3A21H2m4CJwPNlLOavwFPAtnJWs2NU9kJgRjn3IdXXsYBVdSVEdrMngVeSpJf7zKmIZB91kCvPTUB9oK+7/6DzaWaXA63Lu2Mzy3f3De6+A9hRsWqKJOfu5f3iJZLNPnL3v1V1JUSkammIReXpAqxO7BwDuHuRuy+NhjjEnvV9XvzlvFje6PVjZvbjaPzXRuClaN0uY5DNLD8a3/aBma0ys61mNt/M7jCz+nH5zuf7M9h/iSv7neL2HbftLmNTzWxhlN7bzN40s41mtsLM/mBmtcwsL/p9iZltMbN3zax7eT5Y2T2K+Tu/E/2t9zazF8xsnZmtN7PJZrZ3Qt7YJevzzewKM5sX/e3nmdkVxZR5hJm9Ee13s5l9ZGYjisnb2cz+YmaLzWybmS2N6tQvLs+xZvZ3M/sq2t9aM/uXmR1Z3PuN2uXkKG9B1P4amlkNM/utmS2I3sdHZnZYRd5zrN0kSd+5n7i0PAtDoj63MGRrrZl9YmZ3J/t8pHKY2aDob3N2knVpH2MbHffnm1mnuDa3xsweMbP6ZlbTzEbHxeV0SzIcJIrhO6K2sNXMlkex3T5J3qbR/ldHx/IpZta3hDoeFNVtdbTvz83sejOrmc7PQmR30hnkyvMlsK+ZDXX354rJsxI4hzBU4j1gfDH5+gPDgD8ThmKUpC0wEvgH8ARQCBwJXAv0BY6L8r0L3A78Nir3vSj921L2X5J2wBvA34FnCZfof0U4y70fYYjJHUBz4NfA82bW3d2LKlCmlE19M2tewX00IHy5+i9wPeHL4KXAIWbW192XJ+S/gnDF5GHC0KKfAfeZWVN3vyWWycxOAiYDy4F7orxnABPMbG93vyEub3/gLaA28AjwKdCUEOsDgOlR1vOj9EnAYr5vH2+Z2dHuHov7+Pc2hdA+fgMcCFwA5AGrgYOBsVG5vwZeMrO93H1Ded5zGY2L6jIJGAPUJHz2A8u5P0kuWRvZmuRvvDvlE9rcFOA6QhxeANQFvgMOAO6LXo8CXjazju6+EcDMahOOzYcATwN/IAz3uwQ41sz6u/vSKG+dKO8BhP83/41+fwtYC/zgeG1mPyUc7z8H7gYKgMOA24BehNgXyT7urqUSFuBQwthgB+YBjxIORt2T5HXgsWL249EyKMm686N1R8Wl1QFqJ8n7P1Heg+LSjorSzk9l33Hr3gEWJqQtjPKflpA+nXBAfQGwuPQro/zHVfXfKheWuL91acvLpfyd34ny/TEh/ZQo/aEkZW4A2iXE6H8J4/PbRWk1ga8J/4D3TMj7H8KXrC5RmhE6xFuAXknea4243xskWd8KWAW8Usx7G5WQ/lwUw9Pi2xbw0yj/xeV5z1H6QuCdEv5e58elrUmss5bd1kaeiss3KEo7O8k+/gYUppD2v9E+2qVQr39Hea9OSH8xissPgFpx6UOj/CPi0i6J0m5P2MfJUfpf4tIujdJuTMj76yh9flxafcKJnreBmgn5R0X5f1TSZ6FFS6YuGmJRSdz9/4B+hG/gjYGfAw8An5nZe4mXo0sx093fTLHcbe6+HcDC0IYm0dmQ2PYHl6Hcslri7s8kpP2b0KEZ6+4elx47c9elEusjuxoPHFPMUhZ3xL9w98mEM0jJZmN53N0Xx+XdRjgDWotwwyqEttIBeNSjM1lxee8mDAc7OUruQ7gi8Rd3n5VYmMddkXD372K/R5eYmxE62x+QvC3sIJwhjvceIYYfirWtuHRIHsOpvOeyWgfsZ2Y9y7m9pCZZG/nfKq1RuBI4LiEtFpcPunthQjr8MC5PifZxZ/wO3P0FwpfNIWYWuyF3COGL3JiE8u4nnK2OdxzhiuCjQBMzax5b+P5Gx2NLf3simUdDLCqRu39COBOLme1FuPw7EjgceMHM+nlqN0LNK0u5ZnYp8AtCJyLxS1CTsuyrjBYkSSsoZl0svVnlVUeS+KK4L1vf/38s1VrfdRgFwBzCP9oG8R3TKD3RZ9HP2BfFTtHP2UnyfpqQN/aP/+PSKmpm+xAu9R4H7JGw2nfdgmXuviUhLWkMu3tB9Jkli+FU3nNZ/ZIwHOsTM/uKcNbuJeAl1zCldCq2jVShxUn+V5Tl2Nop2se6JPueDfQk/G9YQ4jPJR4Nz4hx9y1mtoAwVC4mdh/JpBLq3qqEdSIZSx3k3cTdvwYmmVlsvPFhwEGEM6yl2ZRqOWZ2DWH85r8IY9KWEoZ6tAUeI/UbM5N1HmKKi5uSZtQobp2mEcs+xcVGcX/LZPkT85YlDmJ5S4pRzKwhYSxxA+CPwCeEoQ9FhLHTycbupiuGU3nPxeWDJG3M3V8ws47AiYQv24OAEcB7ZjYoxS/bUnHlOTZWVEXjsqztK9U2Hnt9DaF9JbOkDGWLZAx1kHczd3cz+4DQQW5bCUWcQxjXeEL8WSUzOz5ZdUrYz5roZ9Mk6zoRLsFJbmpiZq2TnEXuBqxIOHsM0CPJPmJnnr6Kfn4Z/dwvSd4eCXk/j34We1d95MfAnsAF7v6X+BVmVtmXzFN5zxDaWbI2lvQss7uvIYzj/Ft0SfwOwg24JwOJw5ukcpR0bCzv1YHK9iUw0Mwaufv6hHU9CGP/C+LyHmVmDePPIptZHmHu/Ph2/0X0c2MGnnUXqRCNQa4kZnaMme3yBcTM6vH9mKzYJdeNJD/YlscOQsd35zf9qB6/SZI3dvBLVnZsWMeg+EQz+xmh0yG57QfxZGanAPuS/GE3Z1ncI3Wju+SvJsRq7Ml9HwGLgJ+bWeu4vLX5/mafF6LkmYTLwheY2S4d6rixlLEza5aw/lgqdyw+pPaeIbSzbmbWNi5vXeCy+J1FU3n9YIhINKY/NswkXccPKd1XhL9j4rHxcMKMQ5noecIJsWvjE6OZY/YHno+7R+QFwiwtVyfs43LC1Zh4rxBmd7nezHYZvmdm9cwsv+LVF9n9dAa58owBmpnZi4RLT5uA9sCZhOl1JkVjlAHeBwaZ2XWEToK7+1PlLPdZ4PfAq2b2HNAoKjPZGd/PCJecLzWzTYSzCCvcfYq7f25mbwIXRx2OGYSbo04hPFGqdjnrJ9lvFTDUzPYkzPwQm+btW+DmJPnnAR+Y2UOEeDuTMH3a/7j7NwDuvsPCA3QmAx+a2fgo73DC1FS3u/sXUV43s58Tpp36r5nFpnnbgzD04DXCjXb/JpoyLhqasJgQw+cQ2uT+aftEyvGeI/cTprJ7M8pbJ6pf4rCqfGBZdDz5GFhBuJJzCeHM30uV+F4kjruvi4bKnW9mfyMMmetKuN/kE5JfBalqjwDnAjdEN4jH6nwpsAy4IS7vBOBC4NZoDP8HhJtoh7LrOPyNZnYuYaaXz83sL4T/D00IV5SGAj8htaGEIhlFHeTKcw3hsuePCHMY70G4C30W4U7ix+LyXkq4Q/kGwj9CCI+QLo+7CWfMRgB/InQQ/g78he/PWAPg7pvN7AzCHdp/JMyhOZUw1yaEf9RjgbOi398DjgYeJFxqk9z0HWH87hjCJX4jdEp/5e7LkuQfS/iidgVhpopFwC/d/U/xmdz9JTP7MTCacNa4DuFmtwvdfUJC3g/N7EDgRuB0wk2pqwhTqf0nyrPWzI4D7orKrkWYdvBEQvuozA5yqu/5PxYeBvJbQttdQmhf0whfAGI2EdrojwlnLhsSOjYvAr+Pn/lDdourCFc1hhA6gdOAwYSzrBnXQXb3bWZ2DN+3l1MJX6yeIkznFj9zzFYzG0SIx5OB0wid5B8T4rp1wr5fidribwj/J5pH+54f7SPZjbciGc9+OPPWbijQrDPhn98hhDtn33P3o1LYrjHhH8QQwtCQl4Er3X115dVWROJZeOpbR3fvmELeowgzLfzc3R+r1IpliFx8zyIi1VFVnEHej3AG533CGaJU/Z0wxnEk4S70Ownjqg5PdwVFREREJHdVRQf5pWhycszsWcLlmBKZ2aGEeUyPdPd3o7QlhDF+g3T3rIiIiIiky26fxaKcE9qfAHwb6xxH+/kv4YaBE9JVNxERERGRbLlJrxswN0n6nGidiOwGqdwvEJf3HXLsQTC5+J5FRKqjbJkHuQlhCrJEBVTuo5NFREREJMdkSwcZin90a9JpOMzsIjObZmbT9ttvP4/yadESW6oFxbmWFJZqQbGupZRFJK2ypYNcQJhHONEeJD+zjLuPd/f+7t6/Xr16lVo5kaqiOJdcoVgXkd0pWzrIc0k+1ri4sckiIiIiIuWSLR3kV4HWZvajWIKZ9Qf2jtaJiIiIiKTFbp/FwszqEx4UAtAWaGRmp0avX3H3TWY2H5jq7iMA3P3/zOx1YJKZ/ZrvHxTy7901B/KsWbMqtH2vXr3SVBMRERERqUwpnUGOboy41MzSMWNES+CZaDkE6BH3umWUpxZQM2G7M4CpwKPAJGA6cEoa6iNS5dLcxkQylmJdRLJBqkMsZhPO2C41s7+b2bFmVq65Pt19obtbMcvCKE9Hdz8/Ybu17v5zd9/D3Ru5+5nuvqo8dRDJQGlrYyIZTrEuIhkvpQ6yu58HtAYui36+Biwys9vMrEsl1k8kJ6iNSa5QrItINkh5DLK7f0cY3vCome0DnAecC/zGzP4TrXvK3bdUSk1Fqjm1MckVmRDruq9EREpS3lksivh+Yu4dhAd2PAAsNLNj0lExkRynNia5QrEuIhkn5TPI0ewTpwHnA4cD8wkHsYnu/q2ZNQXuBx4mTL8mImWgNpZZdIax8ijWRSTTpTqLxSPAcmAc8DVwtLt3c/e73P1bAHdfA/wJ6FhJdRWpttTGJFco1kUkG6R6Bnl/4NfAk+6+oYR8s4GjK1wrkdyTEW2somdNQWdOpVQZEesiIiVJtYN8KrDM3bcnrjCzWsCe7r7I3TcS5ioWkbJRG5NcoVgXkYyXagd5AXAo8N8k63pH6YkP9hCR1KmNSa5QrEuFbd++ncWLF7NlS5joZNu2bQurtkZSTkXAp4WFhSP79eu3oqorEy/VDnJJk7jnAVvTUBeRXKY2JrlCsS4VtnjxYvLz8+nYsSPRc2b04LAsVFRUZCtXruyxfPnyCcBPq7o+8YrtIJtZL6BPXNKJZtYtIVsecDowrxLqJlKtqY1JrlCsS7pt2bIlvnMsWapGjRreokWLdcuXL+9Z1XVJVNIZ5FOAm6LfHfhdMfkWABenWqCZ9QDGEi6xrQUmALe4+45SttsPGAP8CNgEPAOMisapiWSjSmljIhlIsS5pp85x9VCjRg2n/M/lqDQlVeh2IB9oRLgkNjB6Hb/Udfd93P3NVAozsybAm4QD5MnArcCvgFtK2a4xMAWoBwwn3AE9DPhbKuWKZKi0tzGRDKVYl2rjtddea9ixY8eMO+Mp6VXsGeToDuPYXcbp6tn/gtDJHeru64E3zKwRcLOZ3RWlJXNptN1J7r4WwMzWAC+YWX93n5am+onsNpXUxkQyjmJddoe2bdvuv3r16to1atTwevXqFR199NHrHnnkkUWNGzcuSmc5xx9//MaFCxd+ms59Vqa2bdvuP27cuIVDhgwpaVpFSVDSGOQewJfuvjX6vUTu/lkK5Z0AvJ7QEX4KuBM4EnipmO36ANNinePIvwhnogcD6iBL1qmkNiaScRTrsrs89dRTXwwZMmTDggULah977LFdr7/++jYPPPDAkvg8RUVFuDs1a2qyFCleSd/kPyVMuRP7/ZNilti6VHQD5sYnuPsiwpjixBs24uUB2xLSCgnTg3RPsWyRTFMZbUwkEynWZbfq1KnT9oEDB66bM2dOPYCDDjpo3yuuuKLtAQcc0K1+/foHzJkzp+7q1atrnn766Xu1aNGiV8uWLXtdeeWVexYWFrJ582bLz8/v8+GHH+bF9rd06dJaeXl5ByxZsqTWyy+/nN+qVaudT0T66KOP8g466KB98/Pz+3Tu3Hm/xx9/vHFs3UEHHbTvvffe2zz2+r777mvWr1+/fSF01EeMGNG+adOmvfPz8/t07dq1R3yZ8RYuXFh74MCBnRs3btynQ4cOPe+5556d+xw2bFjHK6+8cs/Y6/j6DRkypNOyZcvqnHHGGV3q16/fd/To0a0AXn/99YZ9+/btlp+f36d169a97rvvvmYAq1evrnnKKad0bNKkSe8999xz/2uvvbbNjh07dtb9gAMO6DZixIj2+fn5fdq1a7f/G2+80eC+++5r1rp1615NmzbtPXbs2GaxemzevNkuuuiidm3atNm/WbNmvc8888wOGzduzJqB4yXdpHc08Fnc7+nQhHBjXqKCaF1x5gNnmlntuMnl+xHmymyaprqJ7G6V0cZEMpFiXXar+fPn137rrbcaDx48uCCW9uyzzzZ96aWXvujdu/eWoqIiGzx48N4tW7Ys/PLLLz/dsGFDjeOPP77LmDFjto0aNWrV8ccfv3bSpEnNDjzwwCUAEydObHLggQduaNu2beHHH3+8s5ytW7fakCFDOp955pmr3n333Xn/+te/Gv7sZz/r3LNnz8969+5d4pSFkydPbvT+++83nDdv3qdNmzbdMWPGjLxmzZolnbDgtNNO23vffffd/PLLL8+cMWNG3uDBg7t27tx568knn1zisInnn39+Qdu2bRvGD7H44osv6gwdOrTLvffe+/X5559fUFBQUOOrr76qAzBy5Mj269evr/nVV199smLFilrHHXdc1zZt2my/+uqrVwHMmjWrwXnnnbfy4Ycf/uaaa67Z89xzz9170KBB6xYsWPDJq6++mn/OOefsc+655xY0bty46LLLLmu3cOHCujNmzPisTp06PmzYsL2vu+66PceNG7ekpDpnimLPILv71NgMEdHvJS5lKNOTpFkx6TF/BloAY82sdTSjxQPAjmjZdYdmF5nZNDObtnLlyjJUT2T3SEcbU5xLNlCsy+5y5plnds7Pz+9zxBFHdDvkkEM23Hbbbcti64YPH766f//+W2rXrs2KFStqvvvuu43Hjx+/qFGjRkVt27YtvPzyy7999tlnmwKcddZZqydPnrzzBNwzzzzTbPjw4WsSy3v77bcbbNq0qeZtt922PC8vz3/6059uGDhw4NqJEyc2S8ybqHbt2v7dd9/VnDlzZp67c8ABB2zZa6+9dnnC5Pz582t/9NFHDceOHbu4fv36PmDAgM1nnnnmqkmTJpVaRjKPPvpo0wEDBqy/+OKL19StW9dbt269Y8CAAZsLCwv55z//2fSuu+5a0qRJk6J9991322WXXbb8ySef3FlO27Ztt1511VWra9Wqxdlnn12wfPnyOrfffvvSevXq+dChQ9fXrl3bZ8+eXbeoqIgnn3yy+dixY79p1arVjiZNmhT99re/Xfb8889nzUnNVB8Uki4FwB5J0huT/MwyAO4+18wuIkzzdjFhaMV4Qqf622K2GR/loX///iV1vkWyluJccoViXVLxxBNPzC/uZrT27dvvHKo5f/78OoWFhdamTZvY0B/c3Vq3br0N4KSTTtpw4YUX2pQpUxq0a9du+5w5c+qdddZZBYn7/Oabb2q3bt16W/x45vbt229bunRp7dLq+tOf/nTDrFmzVlx55ZUdli5dWuf4449fO27cuG+aNm36g5sKFy1aVKdRo0aFTZo02Zm+1157bfv444/rl1ZGMt98802dTp067XJ2e9myZbW2b99uXbp02fk5derUadu333678700b958Zwe+fv36RQDt27cvjKXVrVu3aMOGDTWXLVtWa8uWLTUOOeSQHwyD3bFjR/YPsTCzlZR8VvcH3L1lCtnmkjDW2MzaAw1IGJucZP+PmtkTQBdgBeGpOasJ8yiLZJ1KamMiGUexLpkgft7kvffee3udOnV8zZo1M2rX3rUvW7NmTX7yk58U/O1vf2vaqlWr7QMHDlwX30GNad++/fbly5fX2bFjx86b/r755ps6Xbp02QpQv379HZs2bdp5tX758uU/KGz06NErRo8evWLJkiW1TjnllH1uueWW1n/605+Wxufp0KHDtvXr19cqKCioEavDokWL6rRp02Z7VEZRfBlLly4t8eRn+/btt02bNq1BYnqbNm0Ka9Wq5V988UWdfv36bQFYuHBhnVatWu1yVrs0rVu3LszLyyuaNWvW7E6dOpV5+0xQ0oc4jjIc0FL0KjDKzPLdPfYNbziwGSh1mIa7byG6gcPMziMMEXk6zXUU2V0qo42JZCLFumSUvfbaa/thhx227qKLLmp/zz33LGncuHHR3Llz6y5cuLD24MGDNwKcc845a4YPH77PHnvsseOmm25KOm72qKOO+q5evXo7brzxxtY33XTTt2+88UbDKVOm7HHrrbfOAdh///03v/jii02uuuqqVV9//XXtxx9/vHnsLOzUqVPr79ixww477LBN+fn5RXXr1i1KNrNG586dt/fp02fjVVdd1e6hhx765pNPPsl78sknm0+YMGEBQJ8+fTbdf//9rb799ttlW7dutXHjxrWK37558+bb58+fXxfYAHDBBResue+++9pMmDChyXnnnVewevXqml999VWdAQMGbD7xxBMLfvOb37R9+umnF6xcubLWuHHjWl1xxRVJr9SXpGbNmpxxxhmrLrvssvZ//vOfF7Vt27ZwwYIFtT/66KN6w4YNK25K34xS0jzIN1dCeQ8BVwLPmdmdwN7AzcC98VO/mdl8YKq7j4heNwJuAN4lzF5xNOEBIxe6+y5jgkSyQSW1MZGMo1iXTPT0008vvPLKK9t2796956ZNm2q0a9du29VXX71zzPLAgQO/q1evXtGKFStqn3rqqeuS7SMvL88nT548/5JLLtlr7NixrVu2bLn9wQcfXNC3b98tAL/97W+/PfXUU+u3bt2697777rt52LBha6ZOnZoPsHbt2pqjRo1qv3jx4rp169YtOvzww9ffdNNNy4up61cjR47cq02bNr0bNWpUeN111y095ZRT1gNccsklq6dMmdJon3326bXnnntuPeuss1Y98MADrWPbjho1avmoUaPa33LLLe2uvvrqZbfeeuu3//jHP7649tpr2/3yl7/s2LBhwx033HDDkgEDBmyeMGHCopEjR3bYe++9969bt66fffbZK6+66qpV5fl8x40bt/jaa6/d8+CDD+6+du3aWi1bttx2wQUXrASyooNs7rv3S300B+b9/PBR0zfHP2razBYC77j7+dHrBsBkoD/hgSGfAre5+/OplNm/f3+fNq1iUyXPmjWrQtv36tWr9EyyO2XNOKhUZUKcQ/WJ9WrU5hXrSVSjv29OmjNnDt27/2B46/SqqotU3MyZM5v37t27Y1XXI15JY5CfBq539y+j30vk7qenUmA0AfzAUvJ0THj9HXBsKvsXyRaV1cZEMo1iXUSyTUljkFsAscHkLdH4MZF0UxuTXKFYF5GsUtIY5KPjfj9qt9RGJIeojUmuUKyLSLYp6VHTIiIiIiI5J+UOspntb2ZPmNl8M/su+vmEmelOBZE0UBuTXKFYF5FMl9KT9MxsCGG+4S+BZwkP6mgJnAxMM7PTU51RQkR2pTYmuUKxLiLZINVHTd8JvACc7nHzwpnZ9YQD3F2ADmgi5ac2JrlCsS4iGS/VIRbtgQmeMGly9Ho80C7dFRPJMWpjkisU65LTatas2a9bt249unTpst/AgQM7r1q1atfH5+0Gn3/+eZ0uXbrsVxVlZ4NUzyBPA/YDXk+yrifwUdpqJJKb1MYkVyjWJWPMmjWrXzr316tXr1IfWFK3bt2iuXPnfgYwdOjQjnfffXeLO++8M+kT9NKpsLCQWrVS7fZJSQ8KqR/38hrgKTOrTbj0FRszdgowEjijMispUh2pjUmuUKyLJHfIIYd8N2vWrHqx1zfeeGOryZMnN922bZsNHjx47ZgxY5aOHj26VV5eno8ePXrFiBEj2s+ePbve+++/P++FF17If/TRR5u/8MILC84666wOM2fObLBly5YaJ510UsGYMWOWArRt23b/n/3sZ6vefvvtRhdffPGK7t27bx05cmTHevXqFR188MEbq+6dZ76Svkps5IeTuRvwe+D2hDSAD4AquUQgksXUxiRXKNZFEhQWFvL222/njxgxYhXAc88912j+/Pl5s2bNmuPuDBo0qPOrr77a8Oijj974hz/8oRWwYsaMGfW3bdtWY+vWrfbuu+82/NGPfrQB4N57713SqlWrHYWFhQwYMGDfDz74oN7BBx+8GSAvL69o+vTpnwN07dq1x5gxYxYNHjx448UXX6zhTCUoqYN8AZXwtCMz6wGMBQ4F1gITgFvcfUcp2/UnHEz7EQ6kHwE3uPsH6a6jyG5SKW1MJAMp1kUiW7durdGtW7ceS5YsqdOzZ89NQ4YMWQ/w2muvNXr33Xcb9ejRowfApk2basydOzfv0ksvXX3eeec1KCgoqFG3bl3v1avXxvfee6/+//3f/+WPHTt2EcDEiRObPvbYY80LCwtt5cqVtWfOnJkX6yCfe+65BQCrV6+uuWHDhpqDBw/eCHDBBResnjJlSuOq+RQyX0lP0nss3YWZWRPgTeAzwpQ++wD3EG4WHF3Cdu2j7T4Czo2SRwH/MrNe7v51uusqUtkqo40D8MKjAAAdFUlEQVSJZCLFusj3YmOQV69eXfPYY4/tfMcdd7QcPXr0Cnfnl7/85bJRo0atStymXbt2W8eNG9f8oIMO2ti7d+/Nb775Zv7XX39dt2/fvlvmzp1b5/777281ffr0OS1atNgxbNiwjlu2bNk5CUN+fn4RgLtjZom7lmLs7ifp/QKoBwx19zfc/SHgFuAaM2tUwnaDgfxou3+6+z8J49UaAidWdqVFRERE0qlZs2Y77rvvvkXjxo1rtXXrVjvhhBPW//Wvf22+bt26GgALFiyovWTJkloAAwYM2Dhu3LhWRx111IZBgwZtmDhxYosePXpsqlGjBgUFBTXr1atX1LRp0x3ffPNNrXfeeSfpWeHmzZvvaNiw4Y7XX3+9IcBjjz3WdPe92+yT8u2MZjYcuBDoCuQlrnf3lins5gTgdXdfH5f2FGFezCOBl4rZrjZQSBjHFrMxStPXIakW0tTGRDKeYl0kOOywwzZ3795984QJE5pcdtlla2bPnp134IEHdgOoX79+0eOPP76gbdu2hUceeeSG++67r/XAgQO/a9SoUVHdunX9sMMO2whw6KGHbu7Zs+emLl267NehQ4et/fr1K/bmu0ceeWRh7Ca9gQMHri8un6T+JL0zgUeBx4CB0e81gJ8SxhFPSrG8bsCU+AR3X2Rmm6J1xXWQ/wHcCtxjZrdFab8DCoBnUixbJGOlsY2JZDTFemaZNWtWhffRq1f2PiE8lWnZ0m3Tpk0fx7+eMmXK/NjvN95444obb7xxReI2J5988obCwsKdUyAuXLjw0/j1//jHPxYmK2vJkiWfxL8+/PDDN33++eefxV7fe++9S8v8BnJEqkMsRgH/A1wWvX7A3S8AOgGrgE0p7qcJ4QCYqCBal5S7LwWOBoYB30bLUOA4d1+ZYtkimSxdbUwk0ynWRSTjpdpB7gL8J5ppYgfQCMDdNxCGR1xehjKT3clsxaSHlWZtCI8gnU4YpnFC9Ps/zaxDMdtcZGbTzGzaypXqQ0vGK1cbU5xLFlKsi0jGS7WDvA6oG/2+BOget86AZinupwDYI0l6Y5KfWY4ZRRgOcqq7v+burxHOJu8Afp1sA3cf7+793b1/ixYtUqyeSJUpVxtTnEsWUqyLSMYry6OmexEeDfoi8DszKwS2EcYCpzoX8VzCWOOdoincGkTritMNmO3u22MJ7r7NzGYTpooTyXbpamMimU6xLiIZL9UO8u+BvaLffxf9/gDhaUcfAheluJ9XgVFmlh9dTgMYDmwGppaw3dfAiWZWx923AZhZXaAnxd/YJ5JN0tXGRDKdYl1EMl5KHWR3fx94P/p9LXBy1EGtmzBlW2keAq4EnjOzO4G9gZuBe+P3Y2bzganuPiJKmgCMBCab2QOEy3CXAW2A8WUoXyQjpbGNiWQ0xbqIZIMyPyjEghbAtrIezNy9APgx4UzBS4SHhIwBbkrIWivKE9tuOnA84WEhfyVMA1QfOMbdZ5b1PYhksoq0MZFsoliXXGRm/S688MJ2sde/+93vWl1zzTV7VmWdyuLll1/OP/rooztXdT0qW1keFHIi4XHQ/aLtCs1sOnBb9GS7lLj7Z4S5L0vK0zFJ2lvAW6mWI5Jt0tXGRDKdYl0yxdixY/ulc39XXHFFqfMq16lTx1955ZUmy5YtW96mTZvCipS3fft2ateuXZFdlKqwsJBatVLuLlYbKZ1BNrOLCWd8NwJXAadFPzcCL0brRaSc1MYkVyjWJdfVrFnTzz333JW33357q8R18+bNq3PooYd27dq1a49DDz206xdffFEnMc8111yz589+9rO9DjvssC5Dhw7tVFhYyMUXX9yuZ8+e3bt27drj7rvvbg5w9tlnd3j88ccbAxxzzDH7nHbaaR0BxowZ0/zKK6/cE2DQoEH77Lffft07d+683x/+8IfmsTLq16/f95e//OWevXr16vbWW281fPbZZxt16tRpv379+u377LPPJpuNrNpJdYjFb4Hx7n6suz/k7s9FP48F/gzcUHlVFMkJamOSKxTrkvNGjRq14rnnnmu6evXqmvHpv/jFLzqceeaZq+fNm/fZ8OHDV19yySXtk20/a9as+q+//vr8l156acEf//jH5o0bN97x6aefzpk5c+aciRMntpg7d26dI444YsO7776bD7B8+fI68+bNywP4z3/+0/DII4/cCPD4448vnD179pwZM2Z89vDDD7davnx5TYDNmzfX6Nmz5+ZZs2bNPfzww7+7/PLLO7744ovzP/zww89XrFhRuaesM0SqHeRmwHPFrPsH0DQ91RHJWWpjkisU65LzmjZtWnTaaaetvuOOO1rGp3/88ccNLrroojUAl1xyyZrp06c3TLb98ccfv7Zhw4YO8OabbzZ6+umnm3Xr1q1H3759uxcUFNT67LPP8o455piN77//fsPp06fnde3adXPz5s23f/3117WnT5/eYODAgRsB7rzzzlb77rtvj379+nVfvnx57dmzZ+cB1KxZk/PPP78AYMaMGXnt2rXbuv/++2+tUaMGZ5111urK/GwyRaqDSt4GjgTeSLLuSODdtNVIJDepjUmuUKyLANdff/23BxxwQI8zzjhjVVm3bdCgQVHsd3e3e+65Z9GwYcN2udF13bp1tV566aXGhx9++IY1a9bUmjRpUpMGDRoUNWnSpOjll1/Onzp1av60adPm5ufnFx100EH7bt68uQZAnTp1iuLHHZtZOd9l9ir2DLKZ9YgtwH3AOWb2oJkdZ2Z9o58PAecQZqIQkTJQG5NcoVgX2VWrVq12nHTSSQVPPPHEzrG/ffv2/W7ChAlNAB5++OGm/fv331jafo455ph1Dz74YIutW7cawKxZs+quX7++BkC/fv02Pvzwwy0HDRq08aijjto4bty41gcffPBGgLVr19Zs3Ljxjvz8/KKPP/44b+bMmQ2S7b9Pnz5bFi9eXGf27Nl1AZ566qmcuMpT0hnkTwGPe23AxdHi0euY14iblk1EUqI2JrlCsS6SxA033LB84sSJO5+d/uCDDy4677zzOv7pT39q3axZs8JJkyYtLG0fV1999aqFCxfW3X///bu7uzVt2nT7K6+88iXAj370o43vvfdeo549e27dunXrtnXr1tU84ogjNgAMGzZs3fjx41t07dq1xz777LOld+/e3yXbf/369X3s2LFf/+QnP+nctGnTwoMPPnjjnDlz6qXpI8hY5u7JV5gdWZYduXtJT8KrUv379/dp06ZVaB+zZs2q0Pa9evWq0PaSdlV+vSjdbSwT4hyqT6xXozavWE+iGv19KyRb2/ycOXPo3r17fFKp06tJ5po5c2bz3r17d6zqesQr9gxyJnd4RaoDtTHJFYp1Eck2ZZr52cwOBn5EuMt4DfBvd/+gMiomkovUxiRXKNZFJJOl1EE2swbAM4THPRcCqwlT9dQ0s9eA09x9U6XVUqSaUxuTXKFYF5FskOo8yHcBhwLDgTx3bwPkAWdE6XemWmB0J/NbZrbJzJaa2a1mVuINGWZ2s5l5Mcv1qZYtksHS1sZEMpxiXdKiuHuoJLsUFRUZUFRqxt0s1Q7yMOA6d3/G3YsA3L3I3Z8BfkN4VGipzKwJ8CbhruWTgVuBXwG3lLLpBMKBM36JHURfTfE9iGSytLQxkSygWJcKy8vLY/Xq1eokZ7mioiJbuXJlY8JMNxkl1THIjYFviln3DdAoxf38AqgHDHX39cAbZtYIuNnM7orSduHui4HF8WlmdiMw191npFi2SCZLVxsTyXSKdamwdu3asXjxYlauXAnAtm3bmpeyiWSmIuDTwsLCkVVdkUSpdpBnApeY2Wse93XNwqNVLonWp+IE4PWEjvBThLPBRwIvpbITM2sKHAP8b4rlimS6dLUxkUynWJcKq127Np06dYpP6lhFVZFqKtUO8m8JQxnmmtlk4FugJXAKIShPSHE/3YAp8QnuvsjMNkXrUuogA6cCtQmda5HqIF1tTCTTKdZFJOOl1EF29ylm1hf4HWF8WBtgGfABYbjEZymW1wRYmyS9IFqXqjOAj9x9Xhm2EclYaWxjIhlNsS4i2aDUDrKZ1SAcwBa5+xlpKDPZiHorJj1ZfdoQhmNcV0q+i4CLADp06FDGKorsPhVpY4pzySaKdRHJFqnMYlEDWEiY0L2iCoA9kqQ3JvmZ5WROJ3So/15SJncf7+793b1/ixYtSsoqUtXK3cYU55JlFOsikhVK7SC7eyHwNVA/DeXNJYw13snM2gMNonWpOIPwxKXi7oIWySppbmMiGUuxLiLZItV5kO8EbjCzin5tfxU4zszy49KGA5uBqaVtbGYdgUOAJytYD5FMk642JpLpFOsikvFSncXiWMK4sYVmNp1w13H8mGF39+Ep7Och4ErgOTO7E9gbuBm4N37qNzObD0x19xEJ259BeDTpsynWWyRbpKuNiWQ6xbqIZLxUO8jNgc8TXpeZuxeY2Y+B+wlTuq0FxhA6yYn1Svb46TOAt9x9ZXnKF8lgaWljIllAsS4iGS/Vad6OTleB0RQ+A0vJ07GY9D7pqodIJklnGxPJZIp1EckGqY5BFhERERHJCSl3kM1sfzN7wszmm9l30c8nzKxXZVZQJFeojUmuUKyLSKZLaYiFmQ0Bnga+JNwgt4LwaNCTgWlmdrq7P19ptRSp5tTGJFco1kUkG6R6k96dwAvA6e6+825jM7uecIC7C9ABTaT81MYkVyjWRSTjpTrEoj0wIf5gBmEuHmA80C7dFRPJMWpjkisU6yKS8VLtIE8D9itmXU/go/RURyRnqY1JrlCsi0jGS3WIxTXAU2ZWm3DpKzZm7BRgJHCGme18dKi7b0p3RUWqObUxyRWKdRHJeKl2kP8b/fw9cHtcukU/P0jIn+whHyJSPLUxyRWKdRHJeKl2kC/gh48CFZH0UhuTXKFYF5GMl+qT9B6r5HqI5DS1MckVinURyQa7/Ul6ZtbDzN4ys01mttTMbjWzlC6hmdlQM/vQzDab2Woze83MGlR2nUVEREQkd+zWDrKZNQHeJFxeOxm4FfgVcEsK244EngBeBU4g3MzxBakPExERERERKdXu7lz+AqgHDHX39cAbZtYIuNnM7orSdmFmzYExwBXu/ue4VZMrvcYiIiIiklN29xCLE4DXEzrCTxE6zUeWsN3p0c+JlVUxERERERFIsYOc6hjhFHQD5sYnuPsiYFO0rjgHA58DI8xssZltN7MPzGxAmuolUqXS2MZEMppiXUSyQapnkJeY2V1m1r2C5TUB1iZJL4jWFac1sC8wGrgOOAn4DnjNzFol28DMLjKzaWY2beXKlRWrtUjlK1cbU5xLFlKsi0jGS3UM8sPAOcCvzGwa8AjwVHFjhkuRbP5LKyY9pgbQEDjN3V8DMLP/B3wNXA7cuEsh7uOB8QD9+/ev8jk3x44dW+F9XHHFFWmoiWSocrWxTItzkRQo1kUk46V0Btndb3L3vYFjCEMd7gWWmdnjZjaoDOUVAHskSW9M8jPLMWuin+/E1Wk9MB3oUYbyRTJSGtuYSEZTrItINijTTXruPsXdzyUMebiCMOzhdTNbaGY3m9mepexiLgljjc2sPdCAhLHJCeYQzjBbQroBRWV4CyIZLQ1tTCQrKNZFJJOVdxaL/sARhM5uAfAeYV7i+WZ2dgnbvQocZ2b5cWnDgc3A1BK2e5nQGT46lmBmjYF+wMzyvAGRDFfeNiaSbRTrIpJxUu4gm9leZnaTmX0JvAW0AS4A9nT3c4C9CGPL7i5hNw8BW4HnzGyQmV0E3AzcGz/+zMzmm9kjsdfuPg14AXjEzM4zs8HAi8B2YFyq70Ekk6WpjYlkPMW6iGS6lG7SM7MphG/4i4HHgL+4+9fxedx9h5k9AVxV3H7cvcDMfgzcD7xEGHc8htBJTqxX4lRAZxMOlvcC9YH/AAPdvSCV9yCSydLVxkQynWK9+qnoTei6AV0yUaqzWKwCTgTecPeS7h6eAXQqaUfu/hkwsJQ8HZOkbQQuiRaR6iZtbUwkwynWRSTjpdRBdvfTS88F7r6dMPWaiJSB2pjkCsW6iGSDVM8gA2Bm7YCuQF7iOnd/JV2VEslVamOSKxTrIpLJUh2DnA88DRwbS4p+xl8e0+NDRcpJbUxyhWJdRLJBqrNY/B7oABxOOJidAhxFeALSAuCQyqicSA5RG5NcoVgXkYyXagf5ROA24IPo9VJ3f9fdLyJMvzaqMionkkPUxiRXKNZFJOOl2kFuBXzj7juA74Cmcete4ftLZSJSPmpjkisU6yKS8VLtIH8DNI9+/wL4Sdy6g4Et6ayUSA5SG5NcoVgXkYyX6iwWbwCDgMmEB3tMNLN+hKfiHQHcUznVE8kZamOSKxTrIpLxUu0gX0d4eh3u/lcz2wicCtQDLic8ElREyk9tTHKFYl1EMl6qDwrZBGyKez2Z8O1fRNJAbUxyhWJdRLJBSmOQzWyHmR1UzLp+ZrYj1QLNrIeZvWVmm8xsqZndamYlznlpZh3NzJMsT6VarkgmS2cbE8lkinURyQapDrGwEtbVBgpT2olZE+BN4DPgZGAfwnizGsDoFHbxa+A/ca9XpVKuSBZISxsTyQKKdRHJeMV2kM2sA9AxLqmvmSU+EjQPOI8wuXsqfkEYZzbU3dcDb5hZI+BmM7srSivJ5+7+fopliWS0SmpjIhlHsS4i2aakM8g/B24iPP7TgQeLybcZGJlieScAryd0hJ8C7gSOBF5KcT8i1UFltDGRTKRYF5GsUlIH+QHgWcLlsFnAWdHPeNuARe6+NcXyugFT4hPcfZGZbYrWldZB/ouZNQVWAE8CN7j75hTLljSYNSsxBMpu6tSpFd7HFVdcUeF9ZIDKaGMimUixLiJZpdgOsruvBFYCmFknYJm7b6tgeU2AtUnSC6J1xdkKjAP+BawHjiJMFbQPYSzzLszsIuAigA4dOpS7wiKVJR1tTHEu2UCxLiLZJtVp3r6O/W5m9YERhDO+y4FJ8etT2V2SNCsmPVb+MsL8mDHvmNm3wANm1sfdZyTZZjwwHqB///7F7lskE5S3jSnOJdso1kUkGxQ7zZuZ3WNm8xLS8oGPgD8Cw4HfATPNrGuK5RUAeyRJb0zyM8sleTb6eUAZtxPJCJXUxkQyjmJdRLJNSfMgHw38LSHt10BX4EJ3bw7sCSwEbkyxvLmEMwU7mVl7oEG0riw84adItqmMNiaSiRTrIpJVShpi0RGYnpA2DPjM3R+FMK7MzO4BbkmxvFeBUWaW7+4borThhDuXy3rn1qnRz8Q6imSLjqS/jVW5sWPHVmj7anIDZoU/B6g+nwXVMNb19xWp3ko6g1wL2BJ7Ec0e0Z2EWSgI3/hbp1jeQ4Qb7p4zs0HRTRc3A/fGT/1mZvPN7JG41zdHl+iGRtvdCowBnnP3ik+rIFI1KqONiWQixbqIZJWSOsjzCLNFxPwk+vl6Qr6WwJpUCnP3AuDHQE3ClG63EDq6NyVkrRXliZlLmCf5L8ArwJnA3dFPkWyV9jYmkqEU6yKSVUoaYnE/8Gczawx8C1xJeMLRvxLyHQt8mmqB7v4ZMLCUPB0TXj9FeKCISHVSKW1MJAMp1kUkq5Q0D/JjZtYGuIww88RHwGXuvj2Wx8xaEOYhzooxYyKZRG1McoViXUSyTYnzILv774Hfl7B+JRovJlJuamOSKxTrIpJNShqDLCIiIiKSc9RBFhERERGJY+7V/zkbZrYSKMvjsCtDc2BVFdchU2TCZ7HK3Y+v4jqkVYbEOWTG3zcTZMrnoFivHJny980EmfBZVLs4l6qVEx3kTGBm09y9f1XXIxPos6je9PcN9DlUb/r7fk+fhVRHGmIhIiIiIhJHHWQRERERkTjqIO8+46u6AhlEn0X1pr9voM+hetPf93v6LKTa0RhkEREREZE4OoMsIiIiIhJHHWTZycx6mpmb2VHRazezy+PW1zCzcWb2bbTu5ij9ZDObY2bbzGxhlVS+GGZ2upmdX9X1kMyiWJdcoDgXKb8SHzUtOe9QYEHc66HApcAI4DNgsZnVBCYBrwIXAt/t7kqW4nTCHJ2PVXE9JLMp1iUXKM5FUqQOshTL3d9PSOoGFLj7o7EEM2sHNAKecPd/V6Q8M6sNFLn7jorsR6SsFOuSCxTnImXg7lpKWID9gNeANYRv0nOAy6J1g4E3gBXAeuB94NiE7W8mPGHoYGAasBn4N9AJaAk8D2yM9jswSfkjgdnAVsKTo65N43u7FPgmel8vAccADhwVrXfg8uj3d6LX8cv5SdJujvLXAH4DzI/qPg84L6H8d4BngYuAL4EdQPtoXU/gn8CGaHkGaB237VGxukbrNgJfAZfG5XmsuPppyZ04V6xryZVYV5xr0ZK+pcorkOlL1Mj/CZwI/Dg6AP0mWnc5cCVwXHQgujc6IBwWt/3NwCZgJnAWMARYFB1Q3wJ+DRwLvAmsBurHbTsK2A7cFu3/N9GB6fI0vK+To4PLg1H9b48OrMUdTHsAE4C1wCHR0go4Jcr3qyitXZR/XHSAuxYYBNwZfTY/iavDO8Ay4GPg1OgzbgR0BtZFn88QYBjh8t+HfD/zSuxg+gUwOvp8Ho3SDory7ANMAT6Kq3O7qo6pTFyqa5wr1qs+tjJtqa6xrjjXoiW9S5VXIJMXwjgnB/ZPIW8NwpCV14FH49JvjvZxZFzapVHa7+LSekRpJ0SvG0UHo5sSyrkVWA7UrOB7+y/wakLan4s7mMa9l1UJ23SM8sUfJDsDRex6dmES8GHc63cIZ19aJ+T7K/A5UCcurUt0MB4cvY4dTG+Ny1MbWAncEZf2LPBOVcdSJi/VOc6jfSnWtVT7WFeca9GS3kWzWJRsDeEb+ENmNtzMWsavNLN2ZjbRzJYAhYQzA8cCXRP2sw14L+71/OjnlCRpbaOfhwINgGfMrFZsibZpBbQr75uKbsLoC7yQsOq58u4zwY8JB9PJCXV/C+gTlR8z3d2XJ2w/CJgMFMVtuwBYCPRPyPuv2C/uvp1w9qHcn02OqpZxHtVdsS7xqmWsK85F0k8d5BK4exHh4LiccKlnuZm9Z2Z9zawG8CIwAPgdcDRwIOHO37yEXW2I9hWzLfq5Nq6sWFps2+bRz9mEg3RseTtKb1+Bt9aCcGZkRUJ64uvyag7UJFxSi6/7Y1G5beLyflvM9tclbLsd2Jtd3/fahNfb2PXzlxJU4zgHxbrEqcaxrjgXSTPNYlEKd58LDIvuxj2cMO7qn4TLQX0Jl89ei+U3s3ppKnpN9PMnJD/gfF6Bfa8knB1pmZCe+Lq81kT7P4xw1iFR/EHbi9l+MmF8XKJVFa6d7KKaxjko1iVBNY11xblImqmDnKLoUs8UM7sXeILvvzFvjeUxs70IB5BZaSjy/whjufZ093+mYX87ufsOM5tBuKnjobhVQ9NUxBTC2YbG7v5GObZ/i3DH83R3T3awLQudfSiD6hTnoFiX4lWnWFeci6SfOsglMLNewB+AvxOmm2lCuEw0kzD9z2LgHjO7EcgHbgGWpKNsd18bPdXoT9FB+l3CkJiuwNHufkoFi7gdeM7MHiR8sz8SOL6C+wTA3T83s4eAp8zsLsJUSHmE6ZW6uvvIUnZxM+GGk3+a2aOEMwxtCXc1P+bu75ShOnOBk81sCOHvtdTdl5bl/VR31TzOQbEukWoe64pzkTTSGOSSLSdcCruBMA7tAcLclj91962Eb+eFhLtq/wf4PTA1XYW7+12E+SRPINx88SRhWqH3StouxX1PBq4ATiLM29mX8DSldLmM8JmcC7xCGKs2mPBPobS6zSNM37MJGE/47G8hnNmZX8KmyTxAuOnjUcKUQheVcftcUG3jPNq/Yl1iqm2sK85F0ssqfrVDRERERKT60BlkEREREZE46iCLiIiIiMRRB1lEREREJI46yCIiIiIicdRBFhERERGJow6yiIiIiEgcdZBFREREROKogywiIiIiEkcdZBERERGROP8f1oaq1ZWgZCsAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 266.4x180 with 3 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(1,3, figsize= (3.7,2.5), sharey=True)\n",
"\n",
"h1, h2 = plot_doll_style(axes[0], list(mean_lesion_hpc), title='Striatum')\n",
"h1, h2 = plot_doll_style(axes[1], list(mean_lesion_dls), title='Hippocampus')\n",
"h1, h2 = plot_doll_style(axes[2], list(mean_full), title='Full model')\n",
"\n",
"\n",
"axes[0].set_position([0.1,0.1,0.5,0.7])\n",
"axes[1].set_position([0.8,0.1,0.5,0.7])\n",
"axes[2].set_position([1.5,0.1,0.5,0.7])\n",
"\n",
"#leg = axes[2].legend(['Reward', 'No reward'], fontsize=12, frameon=False, handlelength=.7)\n",
"\n",
"#plt.subplots_adjust(left=0.07, right=.93, wspace=0.25, hspace=0.35)\n",
"leg = fig.legend([h1, h2], ['Reward', 'No reward'], bbox_to_anchor=(2.1, .5), loc = (1,.5), title=\"Previous outcome\")\n",
"leg.set_title('Previous outcome', prop={'size': 12})\n",
"\n",
"plt.tight_layout()\n",
"\n",
"plt.savefig(os.path.join(figure_location, 'StayProbability_DeterministicTask.pdf'), bbox_inches='tight')\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"plt.close()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"SameStart PreviousReward\n",
"False 0.0 0.492189\n",
" 1.0 0.493106\n",
"True 0.0 0.491967\n",
" 1.0 0.492369\n",
"Name: P(SR), dtype: float64"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_control.groupby([\"SameStart\", \"PreviousReward\"])['P(SR)'].mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Regression analysis"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Unnamed: 0</th>\n",
" <th>Action1</th>\n",
" <th>Action2</th>\n",
" <th>Agent</th>\n",
" <th>DLS reliability</th>\n",
" <th>HPC reliability</th>\n",
" <th>P(SR)</th>\n",
" <th>Reward</th>\n",
" <th>StartState</th>\n",
" <th>State2</th>\n",
" <th>Terminus</th>\n",
" <th>Trial</th>\n",
" <th>PreviousAction</th>\n",
" <th>PreviousStart</th>\n",
" <th>PreviousReward</th>\n",
" <th>Stay</th>\n",
" <th>SameStart</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.000000</td>\n",
" <td>0.030000</td>\n",
" <td>0.681250</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>2.0</td>\n",
" <td>4.0</td>\n",
" <td>0.0</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>False</td>\n",
" <td>False</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.000000</td>\n",
" <td>0.059100</td>\n",
" <td>0.671966</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>3.0</td>\n",
" <td>7.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>False</td>\n",
" <td>False</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.021000</td>\n",
" <td>0.087327</td>\n",
" <td>0.675625</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>3.0</td>\n",
" <td>7.0</td>\n",
" <td>2.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>True</td>\n",
" <td>False</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>3</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.029370</td>\n",
" <td>0.114707</td>\n",
" <td>0.670047</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>2.0</td>\n",
" <td>4.0</td>\n",
" <td>3.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>False</td>\n",
" <td>True</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>4</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.055789</td>\n",
" <td>0.141266</td>\n",
" <td>0.669577</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>3.0</td>\n",
" <td>7.0</td>\n",
" <td>4.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>False</td>\n",
" <td>False</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Unnamed: 0 Action1 Action2 Agent DLS reliability HPC reliability \\\n",
"0 0 0.0 0.0 0.0 0.000000 0.030000 \n",
"1 1 1.0 1.0 0.0 0.000000 0.059100 \n",
"2 2 1.0 1.0 0.0 0.021000 0.087327 \n",
"3 3 0.0 0.0 0.0 0.029370 0.114707 \n",
"4 4 1.0 1.0 0.0 0.055789 0.141266 \n",
"\n",
" P(SR) Reward StartState State2 Terminus Trial PreviousAction \\\n",
"0 0.681250 1.0 1.0 2.0 4.0 0.0 NaN \n",
"1 0.671966 1.0 0.0 3.0 7.0 1.0 0.0 \n",
"2 0.675625 1.0 1.0 3.0 7.0 2.0 1.0 \n",
"3 0.670047 0.0 1.0 2.0 4.0 3.0 1.0 \n",
"4 0.669577 1.0 0.0 3.0 7.0 4.0 0.0 \n",
"\n",
" PreviousStart PreviousReward Stay SameStart \n",
"0 NaN NaN False False \n",
"1 1.0 1.0 False False \n",
"2 0.0 1.0 True False \n",
"3 1.0 1.0 False True \n",
"4 1.0 0.0 False False "
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_control.head()"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"data = data_control[['Agent', 'PreviousReward', 'PreviousAction', 'SameStart', 'Action1']]"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.499005\n",
" Iterations 7\n"
]
}
],
"source": [
"mod = smf.logit(formula='Action1 ~ PreviousReward * PreviousAction * SameStart', data=data)\n",
"res = mod.fit()"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.489139\n",
" Iterations 6\n",
"Optimization terminated successfully.\n",
" Current function value: 0.538574\n",
" Iterations 8\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.412727\n",
" Iterations: 35\n",
"Optimization terminated successfully.\n",
" Current function value: 0.507848\n",
" Iterations 7\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.496648\n",
" Iterations: 35\n",
"Optimization terminated successfully.\n",
" Current function value: 0.518046\n",
" Iterations 8\n",
"Optimization terminated successfully.\n",
" Current function value: 0.498005\n",
" Iterations 7\n",
"Optimization terminated successfully.\n",
" Current function value: 0.425808\n",
" Iterations 8\n",
"Optimization terminated successfully.\n",
" Current function value: 0.549229\n",
" Iterations 6\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.444553\n",
" Iterations: 35\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.450759\n",
" Iterations: 35\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.473348\n",
" Iterations: 35\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.456640\n",
" Iterations: 35\n",
"Optimization terminated successfully.\n",
" Current function value: 0.476260\n",
" Iterations 8\n",
"Optimization terminated successfully.\n",
" Current function value: 0.442589\n",
" Iterations 8\n",
"Optimization terminated successfully.\n",
" Current function value: 0.523259\n",
" Iterations 7\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.478907\n",
" Iterations: 35\n",
"Optimization terminated successfully.\n",
" Current function value: 0.471507\n",
" Iterations 21\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.507018\n",
" Iterations: 35\n",
"Optimization terminated successfully.\n",
" Current function value: 0.506375\n",
" Iterations 7\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.453712\n",
" Iterations: 35\n",
"Optimization terminated successfully.\n",
" Current function value: 0.497200\n",
" Iterations 7\n",
"Optimization terminated successfully.\n",
" Current function value: 0.486236\n",
" Iterations 8\n",
"Optimization terminated successfully.\n",
" Current function value: 0.529446\n",
" Iterations 7\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.483552\n",
" Iterations: 35\n",
"Optimization terminated successfully.\n",
" Current function value: 0.525430\n",
" Iterations 7\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.426731\n",
" Iterations: 35\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.492921\n",
" Iterations: 35\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.514083\n",
" Iterations: 35\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.473534\n",
" Iterations: 35\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n",
"/Users/jessegeerts/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py:508: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
" \"Check mle_retvals\", ConvergenceWarning)\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.445715\n",
" Iterations: 35\n",
"Warning: Maximum number of iterations has been exceeded.\n",
" Current function value: 0.490803\n",
" Iterations: 35\n",
"Optimization terminated successfully.\n",
" Current function value: 0.428061\n",
" Iterations 24\n"
]
},
{
"ename": "LinAlgError",
"evalue": "Singular matrix",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mLinAlgError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-27-f72b44884927>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0magent\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mAgent\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munique\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mmod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msmf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlogit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mformula\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Action1 ~ PreviousReward * PreviousAction * SameStart'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mAgent\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0magent\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmod\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 5\u001b[0m \u001b[0mdf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mres\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_dict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mignore_index\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/discrete/discrete_model.py\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, start_params, method, maxiter, full_output, disp, callback, **kwargs)\u001b[0m\n\u001b[1;32m 1832\u001b[0m bnryfit = super(Logit, self).fit(start_params=start_params,\n\u001b[1;32m 1833\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxiter\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaxiter\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfull_output\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1834\u001b[0;31m disp=disp, callback=callback, **kwargs)\n\u001b[0m\u001b[1;32m 1835\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1836\u001b[0m \u001b[0mdiscretefit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mLogitResults\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbnryfit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/discrete/discrete_model.py\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, start_params, method, maxiter, full_output, disp, callback, **kwargs)\u001b[0m\n\u001b[1;32m 218\u001b[0m mlefit = super(DiscreteModel, self).fit(start_params=start_params,\n\u001b[1;32m 219\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxiter\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaxiter\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfull_output\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 220\u001b[0;31m disp=disp, callback=callback, **kwargs)\n\u001b[0m\u001b[1;32m 221\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 222\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mmlefit\u001b[0m \u001b[0;31m# up to subclasses to wrap results\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/models/lib/python3.6/site-packages/statsmodels/base/model.py\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, start_params, method, maxiter, full_output, disp, fargs, callback, retall, skip_hessian, **kwargs)\u001b[0m\n\u001b[1;32m 471\u001b[0m \u001b[0mHinv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcov_params_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mxopt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mretvals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 472\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'newton'\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 473\u001b[0;31m \u001b[0mHinv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mretvals\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Hessian'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mnobs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 474\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mskip_hessian\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 475\u001b[0m \u001b[0mH\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhessian\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxopt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/models/lib/python3.6/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36minv\u001b[0;34m(a)\u001b[0m\n\u001b[1;32m 549\u001b[0m \u001b[0msignature\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'D->D'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m'd->d'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 550\u001b[0m \u001b[0mextobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 551\u001b[0;31m \u001b[0mainv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_umath_linalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 552\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mainv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 553\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/models/lib/python3.6/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36m_raise_linalgerror_singular\u001b[0;34m(err, flag)\u001b[0m\n\u001b[1;32m 95\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 96\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 97\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Singular matrix\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 98\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 99\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_nonposdef\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mLinAlgError\u001b[0m: Singular matrix"
]
}
],
"source": [
"df = pd.DataFrame({}, columns=res.params.keys())\n",
"for agent in data.Agent.unique():\n",
" mod = smf.logit(formula='Action1 ~ PreviousReward * PreviousAction * SameStart', data=data[data.Agent==agent])\n",
" res = mod.fit()\n",
" df = df.append(res.params.to_dict(), ignore_index=True)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 -1.215337\n",
"1 -0.680725\n",
"2 -1.007641\n",
"3 -0.481838\n",
"4 -1.121497\n",
"5 -0.757686\n",
"6 -0.867501\n",
"7 -1.722767\n",
"8 -1.446227\n",
"9 -1.963610\n",
"10 -0.885519\n",
"11 -0.344840\n",
"12 -1.178655\n",
"13 -1.174598\n",
"14 -1.515912\n",
"15 -2.018533\n",
"16 -1.317301\n",
"17 -1.163151\n",
"18 -0.597837\n",
"19 -0.346625\n",
"20 -0.938596\n",
"21 -1.128959\n",
"22 -1.376725\n",
"23 -0.385504\n",
"24 -0.693147\n",
"25 -1.402043\n",
"26 -1.709521\n",
"27 -1.317707\n",
"28 -0.435318\n",
"29 -1.410987\n",
"30 -1.544083\n",
"31 -0.976010\n",
"Name: PreviousReward, dtype: float64"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df[\"PreviousReward\"]"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>Logit Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>Action1</td> <th> No. Observations: </th> <td> 271</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>Logit</td> <th> Df Residuals: </th> <td> 263</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>MLE</td> <th> Df Model: </th> <td> 7</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Wed, 29 Jul 2020</td> <th> Pseudo R-squ.: </th> <td>0.2475</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>15:23:33</td> <th> Log-Likelihood: </th> <td> -133.01</td> \n",
"</tr>\n",
"<tr>\n",
" <th>converged:</th> <td>False</td> <th> LL-Null: </th> <td> -176.75</td> \n",
"</tr>\n",
"<tr>\n",
" <th> </th> <td> </td> <th> LLR p-value: </th> <td>4.054e-16</td>\n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> -0.7732</td> <td> 0.349</td> <td> -2.215</td> <td> 0.027</td> <td> -1.457</td> <td> -0.089</td>\n",
"</tr>\n",
"<tr>\n",
" <th>SameStart[T.True]</th> <td> 0.6678</td> <td> 0.477</td> <td> 1.401</td> <td> 0.161</td> <td> -0.267</td> <td> 1.602</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousReward</th> <td> -0.9760</td> <td> 0.518</td> <td> -1.883</td> <td> 0.060</td> <td> -1.992</td> <td> 0.040</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousReward:SameStart[T.True]</th> <td> -1.9391</td> <td> 0.948</td> <td> -2.046</td> <td> 0.041</td> <td> -3.797</td> <td> -0.081</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousAction</th> <td> 0.2426</td> <td> 0.530</td> <td> 0.458</td> <td> 0.647</td> <td> -0.796</td> <td> 1.281</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousAction:SameStart[T.True]</th> <td> -0.4474</td> <td> 0.737</td> <td> -0.607</td> <td> 0.544</td> <td> -1.893</td> <td> 0.998</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousReward:PreviousAction</th> <td> 2.1426</td> <td> 0.773</td> <td> 2.772</td> <td> 0.006</td> <td> 0.628</td> <td> 3.657</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousReward:PreviousAction:SameStart[T.True]</th> <td> 22.1332</td> <td> 8544.817</td> <td> 0.003</td> <td> 0.998</td> <td>-1.67e+04</td> <td> 1.68e+04</td>\n",
"</tr>\n",
"</table>"
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" Logit Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Action1 No. Observations: 271\n",
"Model: Logit Df Residuals: 263\n",
"Method: MLE Df Model: 7\n",
"Date: Wed, 29 Jul 2020 Pseudo R-squ.: 0.2475\n",
"Time: 15:23:33 Log-Likelihood: -133.01\n",
"converged: False LL-Null: -176.75\n",
" LLR p-value: 4.054e-16\n",
"===================================================================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"-------------------------------------------------------------------------------------------------------------------\n",
"Intercept -0.7732 0.349 -2.215 0.027 -1.457 -0.089\n",
"SameStart[T.True] 0.6678 0.477 1.401 0.161 -0.267 1.602\n",
"PreviousReward -0.9760 0.518 -1.883 0.060 -1.992 0.040\n",
"PreviousReward:SameStart[T.True] -1.9391 0.948 -2.046 0.041 -3.797 -0.081\n",
"PreviousAction 0.2426 0.530 0.458 0.647 -0.796 1.281\n",
"PreviousAction:SameStart[T.True] -0.4474 0.737 -0.607 0.544 -1.893 0.998\n",
"PreviousReward:PreviousAction 2.1426 0.773 2.772 0.006 0.628 3.657\n",
"PreviousReward:PreviousAction:SameStart[T.True] 22.1332 8544.817 0.003 0.998 -1.67e+04 1.68e+04\n",
"===================================================================================================================\n",
"\"\"\""
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res.summary()"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 -1.215337\n",
"1 -0.680725\n",
"2 -1.007641\n",
"3 -0.481838\n",
"4 -1.121497\n",
"5 -0.757686\n",
"6 -0.867501\n",
"7 -1.722767\n",
"8 -1.446227\n",
"9 -1.963610\n",
"10 -0.885519\n",
"11 -0.344840\n",
"12 -1.178655\n",
"13 -1.174598\n",
"14 -1.515912\n",
"15 -2.018533\n",
"16 -1.317301\n",
"17 -1.163151\n",
"18 -0.597837\n",
"19 -0.346625\n",
"20 -0.938596\n",
"21 -1.128959\n",
"22 -1.376725\n",
"23 -0.385504\n",
"24 -0.693147\n",
"25 -1.402043\n",
"26 -1.709521\n",
"27 -1.317707\n",
"28 -0.435318\n",
"29 -1.410987\n",
"30 -1.544083\n",
"31 -0.976010\n",
"Name: PreviousReward, dtype: float64"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df['PreviousReward']"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Ttest_1sampResult(statistic=-13.497324667045547, pvalue=1.6042574528289812e-14)"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from scipy.stats import ttest_ind\n",
"from scipy.stats import ttest_1samp\n",
"\n",
"ttest_1samp(df['PreviousReward'],0)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Ttest_1sampResult(statistic=-4.5382898405157395, pvalue=8.029457107185251e-05)"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ttest_1samp(df['PreviousReward:SameStart[T.True]'], 0)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Help on function ttest_1samp in module scipy.stats.stats:\n",
"\n",
"ttest_1samp(a, popmean, axis=0, nan_policy='propagate')\n",
" Calculate the T-test for the mean of ONE group of scores.\n",
" \n",
" This is a two-sided test for the null hypothesis that the expected value\n",
" (mean) of a sample of independent observations `a` is equal to the given\n",
" population mean, `popmean`.\n",
" \n",
" Parameters\n",
" ----------\n",
" a : array_like\n",
" sample observation\n",
" popmean : float or array_like\n",
" expected value in null hypothesis. If array_like, then it must have the\n",
" same shape as `a` excluding the axis dimension\n",
" axis : int or None, optional\n",
" Axis along which to compute test. If None, compute over the whole\n",
" array `a`.\n",
" nan_policy : {'propagate', 'raise', 'omit'}, optional\n",
" Defines how to handle when input contains nan. 'propagate' returns nan,\n",
" 'raise' throws an error, 'omit' performs the calculations ignoring nan\n",
" values. Default is 'propagate'.\n",
" \n",
" Returns\n",
" -------\n",
" statistic : float or array\n",
" t-statistic\n",
" pvalue : float or array\n",
" two-tailed p-value\n",
" \n",
" Examples\n",
" --------\n",
" >>> from scipy import stats\n",
" \n",
" >>> np.random.seed(7654567) # fix seed to get the same result\n",
" >>> rvs = stats.norm.rvs(loc=5, scale=10, size=(50,2))\n",
" \n",
" Test if mean of random sample is equal to true mean, and different mean.\n",
" We reject the null hypothesis in the second case and don't reject it in\n",
" the first case.\n",
" \n",
" >>> stats.ttest_1samp(rvs,5.0)\n",
" (array([-0.68014479, -0.04323899]), array([ 0.49961383, 0.96568674]))\n",
" >>> stats.ttest_1samp(rvs,0.0)\n",
" (array([ 2.77025808, 4.11038784]), array([ 0.00789095, 0.00014999]))\n",
" \n",
" Examples using axis and non-scalar dimension for population mean.\n",
" \n",
" >>> stats.ttest_1samp(rvs,[5.0,0.0])\n",
" (array([-0.68014479, 4.11038784]), array([ 4.99613833e-01, 1.49986458e-04]))\n",
" >>> stats.ttest_1samp(rvs.T,[5.0,0.0],axis=1)\n",
" (array([-0.68014479, 4.11038784]), array([ 4.99613833e-01, 1.49986458e-04]))\n",
" >>> stats.ttest_1samp(rvs,[[5.0],[0.0]])\n",
" (array([[-0.68014479, -0.04323899],\n",
" [ 2.77025808, 4.11038784]]), array([[ 4.99613833e-01, 9.65686743e-01],\n",
" [ 7.89094663e-03, 1.49986458e-04]]))\n",
"\n"
]
}
],
"source": [
"help(ttest_1samp)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.498717\n",
" Iterations 7\n"
]
}
],
"source": [
"mod = smf.logit(formula='Action1 ~ PreviousReward * PreviousAction * SameStart * Agent', data=data)\n",
"res = mod.fit()"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>Logit Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>Action1</td> <th> No. Observations: </th> <td> 10840</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>Logit</td> <th> Df Residuals: </th> <td> 10824</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>MLE</td> <th> Df Model: </th> <td> 15</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Wed, 29 Jul 2020</td> <th> Pseudo R-squ.: </th> <td>0.2800</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>15:23:35</td> <th> Log-Likelihood: </th> <td> -5406.1</td>\n",
"</tr>\n",
"<tr>\n",
" <th>converged:</th> <td>True</td> <th> LL-Null: </th> <td> -7508.1</td>\n",
"</tr>\n",
"<tr>\n",
" <th> </th> <td> </td> <th> LLR p-value: </th> <td> 0.000</td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> -0.4366</td> <td> 0.105</td> <td> -4.157</td> <td> 0.000</td> <td> -0.642</td> <td> -0.231</td>\n",
"</tr>\n",
"<tr>\n",
" <th>SameStart[T.True]</th> <td> 0.4108</td> <td> 0.149</td> <td> 2.760</td> <td> 0.006</td> <td> 0.119</td> <td> 0.702</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousReward</th> <td> -1.0693</td> <td> 0.166</td> <td> -6.456</td> <td> 0.000</td> <td> -1.394</td> <td> -0.745</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousReward:SameStart[T.True]</th> <td> -2.3645</td> <td> 0.356</td> <td> -6.638</td> <td> 0.000</td> <td> -3.063</td> <td> -1.666</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousAction</th> <td> 0.5712</td> <td> 0.154</td> <td> 3.699</td> <td> 0.000</td> <td> 0.269</td> <td> 0.874</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousAction:SameStart[T.True]</th> <td> -0.8381</td> <td> 0.218</td> <td> -3.843</td> <td> 0.000</td> <td> -1.265</td> <td> -0.411</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousReward:PreviousAction</th> <td> 2.3226</td> <td> 0.243</td> <td> 9.567</td> <td> 0.000</td> <td> 1.847</td> <td> 2.798</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousReward:PreviousAction:SameStart[T.True]</th> <td> 4.4519</td> <td> 0.502</td> <td> 8.867</td> <td> 0.000</td> <td> 3.468</td> <td> 5.436</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Agent</th> <td> 0.0032</td> <td> 0.005</td> <td> 0.683</td> <td> 0.495</td> <td> -0.006</td> <td> 0.012</td>\n",
"</tr>\n",
"<tr>\n",
" <th>SameStart[T.True]:Agent</th> <td> -0.0003</td> <td> 0.007</td> <td> -0.043</td> <td> 0.966</td> <td> -0.013</td> <td> 0.013</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousReward:Agent</th> <td> -0.0011</td> <td> 0.007</td> <td> -0.144</td> <td> 0.885</td> <td> -0.015</td> <td> 0.013</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousReward:SameStart[T.True]:Agent</th> <td> -0.0007</td> <td> 0.016</td> <td> -0.045</td> <td> 0.964</td> <td> -0.031</td> <td> 0.030</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousAction:Agent</th> <td> 0.0009</td> <td> 0.007</td> <td> 0.128</td> <td> 0.898</td> <td> -0.013</td> <td> 0.014</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousAction:SameStart[T.True]:Agent</th> <td> 0.0057</td> <td> 0.010</td> <td> 0.590</td> <td> 0.555</td> <td> -0.013</td> <td> 0.025</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousReward:PreviousAction:Agent</th> <td> -0.0020</td> <td> 0.011</td> <td> -0.188</td> <td> 0.851</td> <td> -0.023</td> <td> 0.019</td>\n",
"</tr>\n",
"<tr>\n",
" <th>PreviousReward:PreviousAction:SameStart[T.True]:Agent</th> <td> 0.0063</td> <td> 0.022</td> <td> 0.282</td> <td> 0.778</td> <td> -0.037</td> <td> 0.050</td>\n",
"</tr>\n",
"</table>"
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" Logit Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Action1 No. Observations: 10840\n",
"Model: Logit Df Residuals: 10824\n",
"Method: MLE Df Model: 15\n",
"Date: Wed, 29 Jul 2020 Pseudo R-squ.: 0.2800\n",
"Time: 15:23:35 Log-Likelihood: -5406.1\n",
"converged: True LL-Null: -7508.1\n",
" LLR p-value: 0.000\n",
"=========================================================================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"-------------------------------------------------------------------------------------------------------------------------\n",
"Intercept -0.4366 0.105 -4.157 0.000 -0.642 -0.231\n",
"SameStart[T.True] 0.4108 0.149 2.760 0.006 0.119 0.702\n",
"PreviousReward -1.0693 0.166 -6.456 0.000 -1.394 -0.745\n",
"PreviousReward:SameStart[T.True] -2.3645 0.356 -6.638 0.000 -3.063 -1.666\n",
"PreviousAction 0.5712 0.154 3.699 0.000 0.269 0.874\n",
"PreviousAction:SameStart[T.True] -0.8381 0.218 -3.843 0.000 -1.265 -0.411\n",
"PreviousReward:PreviousAction 2.3226 0.243 9.567 0.000 1.847 2.798\n",
"PreviousReward:PreviousAction:SameStart[T.True] 4.4519 0.502 8.867 0.000 3.468 5.436\n",
"Agent 0.0032 0.005 0.683 0.495 -0.006 0.012\n",
"SameStart[T.True]:Agent -0.0003 0.007 -0.043 0.966 -0.013 0.013\n",
"PreviousReward:Agent -0.0011 0.007 -0.144 0.885 -0.015 0.013\n",
"PreviousReward:SameStart[T.True]:Agent -0.0007 0.016 -0.045 0.964 -0.031 0.030\n",
"PreviousAction:Agent 0.0009 0.007 0.128 0.898 -0.013 0.014\n",
"PreviousAction:SameStart[T.True]:Agent 0.0057 0.010 0.590 0.555 -0.013 0.025\n",
"PreviousReward:PreviousAction:Agent -0.0020 0.011 -0.188 0.851 -0.023 0.019\n",
"PreviousReward:PreviousAction:SameStart[T.True]:Agent 0.0063 0.022 0.282 0.778 -0.037 0.050\n",
"=========================================================================================================================\n",
"\"\"\""
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res.summary()"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.17537174578342e-11"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res.pvalues['PreviousReward:SameStart[T.True]']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}