simulated experiment- average or average, that is a questions.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn import linear_model\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "distribution_m = (4,1.2)\n",
    "distribution_noise = (3,1.1)\n",
    "distribution_experiment_noise = (0,1)\n",
    "n_subjects = 20 # subjects\n",
    "n_trials = int(100) # number of trials per subject\n",
    "responses = (0.5,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = {}\n",
    "labels = {}\n",
    "for n in range(n_subjects):\n",
    "    # for each trial, there is internel noise contributes to the response\n",
    "    internal_noise = np.random.normal(loc=distribution_noise[0],scale=distribution_noise[1],size=n_trials)\n",
    "    # for each trial, there is a \"m\" drew from a N(4,1.2) distribution\n",
    "    internal_m = np.random.normal(loc=distribution_m[0],scale=distribution_m[1],size=n_trials)\n",
    "    # subject responses is a binomial distribution\n",
    "    subject_responses = np.random.choice([-1,1],size=n_trials,p=[responses[0],responses[0]])\n",
    "    # there exists some fixed experimental noise for all the subjects drew from a unit normal distribution\n",
    "    noise = np.random.normal(loc=distribution_experiment_noise[0],scale=distribution_experiment_noise[1],size=n_trials)\n",
    "    # since the formula: mx + b = y, where x is the measure, b is the internal noise, and y is the response,\n",
    "    # the inverted measure simulation is: (y - b) / m + experimental noise\n",
    "    measures = (subject_responses - internal_noise)/internal_m + noise\n",
    "    results[n] = measures\n",
    "    labels[n] = subject_responses"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# fit a model for each subject and then plot the scores of them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import StratifiedKFold,cross_val_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_scores = []\n",
    "for n in range(n_subjects):\n",
    "    model = linear_model.LogisticRegression(class_weight='balanced',random_state=12345,max_iter=int(1e5))\n",
    "    cv = StratifiedKFold(n_splits=5,shuffle=True,random_state=12345)\n",
    "    scores = cross_val_score(model,results[n].reshape(-1,1),labels[n],cv=cv,scoring='roc_auc',)\n",
    "    result_scores.append(scores)\n",
    "result_scores = np.array(result_scores)\n",
    "result_scores = pd.DataFrame(np.array([result_scores.flatten(),np.repeat(np.arange(n_subjects),5)]).T,\n",
    "                            columns=['scores','subjects'])\n",
    "result_scores['subjects'] = result_scores['subjects'].astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "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 tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"8\" halign=\"left\">scores</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>subjects</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.650000</td>\n",
       "      <td>0.036056</td>\n",
       "      <td>0.600000</td>\n",
       "      <td>0.640000</td>\n",
       "      <td>0.650000</td>\n",
       "      <td>0.660000</td>\n",
       "      <td>0.700000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.563596</td>\n",
       "      <td>0.162300</td>\n",
       "      <td>0.409091</td>\n",
       "      <td>0.460000</td>\n",
       "      <td>0.488889</td>\n",
       "      <td>0.660000</td>\n",
       "      <td>0.800000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.711111</td>\n",
       "      <td>0.085292</td>\n",
       "      <td>0.595960</td>\n",
       "      <td>0.646465</td>\n",
       "      <td>0.747475</td>\n",
       "      <td>0.777778</td>\n",
       "      <td>0.787879</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.657172</td>\n",
       "      <td>0.143031</td>\n",
       "      <td>0.511111</td>\n",
       "      <td>0.515152</td>\n",
       "      <td>0.677778</td>\n",
       "      <td>0.745455</td>\n",
       "      <td>0.836364</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.446141</td>\n",
       "      <td>0.151034</td>\n",
       "      <td>0.281818</td>\n",
       "      <td>0.350000</td>\n",
       "      <td>0.388889</td>\n",
       "      <td>0.580000</td>\n",
       "      <td>0.630000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.628872</td>\n",
       "      <td>0.217203</td>\n",
       "      <td>0.458333</td>\n",
       "      <td>0.462963</td>\n",
       "      <td>0.511364</td>\n",
       "      <td>0.768519</td>\n",
       "      <td>0.943182</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.678000</td>\n",
       "      <td>0.117558</td>\n",
       "      <td>0.520000</td>\n",
       "      <td>0.640000</td>\n",
       "      <td>0.640000</td>\n",
       "      <td>0.790000</td>\n",
       "      <td>0.800000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.628283</td>\n",
       "      <td>0.049690</td>\n",
       "      <td>0.565657</td>\n",
       "      <td>0.595960</td>\n",
       "      <td>0.626263</td>\n",
       "      <td>0.666667</td>\n",
       "      <td>0.686869</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.546020</td>\n",
       "      <td>0.094183</td>\n",
       "      <td>0.418182</td>\n",
       "      <td>0.520000</td>\n",
       "      <td>0.536364</td>\n",
       "      <td>0.577778</td>\n",
       "      <td>0.677778</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.444000</td>\n",
       "      <td>0.105972</td>\n",
       "      <td>0.340000</td>\n",
       "      <td>0.350000</td>\n",
       "      <td>0.460000</td>\n",
       "      <td>0.470000</td>\n",
       "      <td>0.600000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.718768</td>\n",
       "      <td>0.091461</td>\n",
       "      <td>0.572727</td>\n",
       "      <td>0.711111</td>\n",
       "      <td>0.730000</td>\n",
       "      <td>0.760000</td>\n",
       "      <td>0.820000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.674000</td>\n",
       "      <td>0.059414</td>\n",
       "      <td>0.610000</td>\n",
       "      <td>0.630000</td>\n",
       "      <td>0.660000</td>\n",
       "      <td>0.720000</td>\n",
       "      <td>0.750000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.675152</td>\n",
       "      <td>0.221614</td>\n",
       "      <td>0.300000</td>\n",
       "      <td>0.676768</td>\n",
       "      <td>0.722222</td>\n",
       "      <td>0.828283</td>\n",
       "      <td>0.848485</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.536768</td>\n",
       "      <td>0.059410</td>\n",
       "      <td>0.466667</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.522222</td>\n",
       "      <td>0.585859</td>\n",
       "      <td>0.609091</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.714000</td>\n",
       "      <td>0.105972</td>\n",
       "      <td>0.650000</td>\n",
       "      <td>0.650000</td>\n",
       "      <td>0.670000</td>\n",
       "      <td>0.700000</td>\n",
       "      <td>0.900000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.692727</td>\n",
       "      <td>0.108928</td>\n",
       "      <td>0.515152</td>\n",
       "      <td>0.676768</td>\n",
       "      <td>0.718182</td>\n",
       "      <td>0.755556</td>\n",
       "      <td>0.797980</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.720626</td>\n",
       "      <td>0.166088</td>\n",
       "      <td>0.509091</td>\n",
       "      <td>0.590000</td>\n",
       "      <td>0.755556</td>\n",
       "      <td>0.866667</td>\n",
       "      <td>0.881818</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.630101</td>\n",
       "      <td>0.132697</td>\n",
       "      <td>0.494949</td>\n",
       "      <td>0.505051</td>\n",
       "      <td>0.626263</td>\n",
       "      <td>0.733333</td>\n",
       "      <td>0.790909</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.569293</td>\n",
       "      <td>0.149970</td>\n",
       "      <td>0.427273</td>\n",
       "      <td>0.433333</td>\n",
       "      <td>0.533333</td>\n",
       "      <td>0.707071</td>\n",
       "      <td>0.745455</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.568889</td>\n",
       "      <td>0.154502</td>\n",
       "      <td>0.409091</td>\n",
       "      <td>0.490909</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.644444</td>\n",
       "      <td>0.800000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         scores                                                              \\\n",
       "          count      mean       std       min       25%       50%       75%   \n",
       "subjects                                                                      \n",
       "0           5.0  0.650000  0.036056  0.600000  0.640000  0.650000  0.660000   \n",
       "1           5.0  0.563596  0.162300  0.409091  0.460000  0.488889  0.660000   \n",
       "2           5.0  0.711111  0.085292  0.595960  0.646465  0.747475  0.777778   \n",
       "3           5.0  0.657172  0.143031  0.511111  0.515152  0.677778  0.745455   \n",
       "4           5.0  0.446141  0.151034  0.281818  0.350000  0.388889  0.580000   \n",
       "5           5.0  0.628872  0.217203  0.458333  0.462963  0.511364  0.768519   \n",
       "6           5.0  0.678000  0.117558  0.520000  0.640000  0.640000  0.790000   \n",
       "7           5.0  0.628283  0.049690  0.565657  0.595960  0.626263  0.666667   \n",
       "8           5.0  0.546020  0.094183  0.418182  0.520000  0.536364  0.577778   \n",
       "9           5.0  0.444000  0.105972  0.340000  0.350000  0.460000  0.470000   \n",
       "10          5.0  0.718768  0.091461  0.572727  0.711111  0.730000  0.760000   \n",
       "11          5.0  0.674000  0.059414  0.610000  0.630000  0.660000  0.720000   \n",
       "12          5.0  0.675152  0.221614  0.300000  0.676768  0.722222  0.828283   \n",
       "13          5.0  0.536768  0.059410  0.466667  0.500000  0.522222  0.585859   \n",
       "14          5.0  0.714000  0.105972  0.650000  0.650000  0.670000  0.700000   \n",
       "15          5.0  0.692727  0.108928  0.515152  0.676768  0.718182  0.755556   \n",
       "16          5.0  0.720626  0.166088  0.509091  0.590000  0.755556  0.866667   \n",
       "17          5.0  0.630101  0.132697  0.494949  0.505051  0.626263  0.733333   \n",
       "18          5.0  0.569293  0.149970  0.427273  0.433333  0.533333  0.707071   \n",
       "19          5.0  0.568889  0.154502  0.409091  0.490909  0.500000  0.644444   \n",
       "\n",
       "                    \n",
       "               max  \n",
       "subjects            \n",
       "0         0.700000  \n",
       "1         0.800000  \n",
       "2         0.787879  \n",
       "3         0.836364  \n",
       "4         0.630000  \n",
       "5         0.943182  \n",
       "6         0.800000  \n",
       "7         0.686869  \n",
       "8         0.677778  \n",
       "9         0.600000  \n",
       "10        0.820000  \n",
       "11        0.750000  \n",
       "12        0.848485  \n",
       "13        0.609091  \n",
       "14        0.900000  \n",
       "15        0.797980  \n",
       "16        0.881818  \n",
       "17        0.790909  \n",
       "18        0.745455  \n",
       "19        0.800000  "
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_scores.groupby(['subjects']).describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'grand average AUC = 0.62')"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3X98XHWd7/HXO01LSymtTQN2WsovqyvLuqhd8LdswSyKgr8X0r3q1Ssrd0HXXVwFFVlWRKTr6uNeRMHfu42Irj8QuBAvlYX1KktRKS0FLJGWdgq0KU1/Udo0n/vHORmmk0kyM81JJpP38/HIIzPnnO/3fGYymc/5fs85368iAjMzM4CmsQ7AzMzqh5OCmZkVOCmYmVmBk4KZmRU4KZiZWYGTgpmZFTgpWN2Q9Jik08c6DrOJzEnBbIxJOlZSn6SvlCw/RlJIai5Z/m1Jny16PlfSNyRtkrRD0kOS/lHS9BGO87S07t2SfiHp6GG2/4ikP0jaJWmNpBemy8+U9J+Stkl6QtL1kmaMZKxWOycFq0npF1WjkDRpDHb7HuBp4BxJh1RTUNJs4FfANOCVETEDeAMwCzh+pAKUNAf4EfBpYDawAvj+ENv/D+ADwJnAYcCbgS3p6pnAZ4Ec8GJgPnD1SMVqB8dJwQokvUzSb9OjzR9I+n7/EamkUyVtkPRxSU8A35L0PEk3S9os6en08fyi+u6U9E+SfpnW2Zl+ufSv/2+S1knqlvTJYWI7M41tu6THJV1WtO42SReUbH+/pLenj/9I0s8lbZX0sKR3F233bUnXSrpV0i7gz4faV1rmPUVxf7q420tSk6RPSHo0XX9j+sU9lPcAnwL2AW8ZZttSfwfsAP4qIh4DiIjHI+IjEbGyyrqG8nZgdUT8ICL2AJcBfyrpj0o3lNQEfAb4aEQ8GIlHI2JrGl9HRNwWEbsj4mngeuDVIxirHQQnBQNA0hTgx8C3SY4Evwe8rWSz56frjgbOI/n8fCt9vgB4BvjfJWXagf8OHAFMAS5K93cCcC3w30iOGFtIjhgHs4vky3MWydHn+ZLemq7rAM4tei0npDHdknah/Dzd5oh0u69I+uOSGK8AZgD/OdS+0rq/AiwB5pIc9c4rquvDwFuB16ev62ngmsFelKTXpq/7BuDGdL/VOB34UUT0VVog7bYZ7OcTgxT7Y+D+/icRsQt4NF1ean76c2KaVP+QdmcN9n3zOmB1pfFbxiLCP/6B5B9zI6CiZf8JfDZ9fCqwF5g6RB0nAU8XPb8T+FTR8/8J3JY+vhS4oWjd9LT+0yuM90vAv6SPZ5B8kR+dPr8C+Gb6+C+Bu0vKfg34TPr428B3q9jXpcD3itYdWhw3sAY4rWj9XJIWQPMgdX8d+En6+JXptkekz48BorRsGnP/3+X3wIdG4fPxDeDzJct+CbyvzLavSuO+hSSxHgM8AnywzLZvIEmcLxzr/wH/JD9uKVi/HLAx0v/U1OMl22yOpOsAAEmHSvpa2pWyHbgLmFXSL/9E0ePdJP3L/fsr1B/JkWf3YMFJOiU9ublZUg/wIWBOWnYHyRfQOenm5wDL0sdHA6cUHw2THOU/f7DXOdS+ysS9uyTuo4EfF+1rDbAfOLLMa5oGvKs/1oj4FbCepOUC0Jv+nlxSdDJJ8iDd99zSujOwEzi8ZNnhJF1XpZ5Jf38hIrZF0q31NeBNxRtJegVJC+6dEfHIyIZrtXJSsH6bgHmSVLTsqJJtSofU/XvgRcApEXE4SWsDQAxvU3H9kg4l6UIaTAdwE3BURMwEvlqyn+8B50p6JclJ11+kyx8H/iMiZhX9HBYR5w/xuoba1yaKurnSL/biuB8H3liyv6kRsbHMa3obyRfrV9KrcJ4g6Yrq70LaRPLlf0xJuWOBdenj/wu8bYiumQEk7Rzi55JBiq0G/rSojukkJ7LLdfs8TNJ6GnQIZkkvJXmP3x8Rd1Qau2XPScH6/YrkiPYCSc2SzgZOHqbMDJKjwm3pydTPVLG/HwJvlvSa9HzG5Qz9eZwBbI2IPZJO5rmj6X63khylXw58P57rY78ZeGF6Unty+vNnkl5c475+CLxF0qvSuP+RA5PTV4ErlF6uKak1fS/LeS/wTeBPSLreTiI54XqSpD+JiP3Av6f1taSxnwucAPyftI4vkiSW7xTtc56kL0p6SbmdpklxsJ/PDRLrj0nOEbxD0lSSbrSVEfFQmfp3k1yZ9A+SZii5+OCDJH8LJJ0I3AZcGBE/G2R/NkacFAyAiNhLcoXJB4BtwF+R/BM/O0SxL5EclW8Bfk3yj17p/lYDf0NyVL6JpF95wxBF/idwuaQdJF9IN5bU9yzJJZOnp3X2L98BtJF0KeVJurOuAoa69HPQfaVxX0hyYngTSffJUzz3Pn2Z5Ai4My3/a+CU0h1ImgecBnwpIp4o+rmP5H18b1EsW4GV6X4uAM6MiCfTeLaS9OHvA+5J93kH0AOsHeI1ViUiNgPvIDlf83T6mvq765D0VUlfLSpyAUmXU57kgKODJAFC0sJsBb5R1ELxieY6oQO7kM2eI+ke4KsR8a2xjqVeSTqMJIkujIg/jHU8ZgfLLQUrkPR6Sc9Pu4/eC7yEKo7+JwpJb0lPsk8HlgIPAI+NbVRmI8NJwYq9iORa9B6SJv47I2LT2IZUl84m6RbJAwuBc8JNbmsQ7j4yM7MCtxTMzKxg3A1qNmfOnDjmmGPGOgwzs3Hlvvvu2xIRrcNtN+6SwjHHHMOKFSvGOgwzs3FF0rrht3L3kZmZFXFSMDOzAicFMzMrcFIwM7MCJwUzMytwUjAzswInBTMzKxh39ylYtq699lq6uroA2LgxmRdm3rxkCuLjjjuO888/f9CyZjb+OSnYoPbs2TP8RmbWUJwU7ADFLYGPfexjAFx99dVjFY6ZjTKfUzAzswInBTMzK3BSMDOzAicFMzMrcFIwM7MCJwUzMytwUjAzswInBTMbdd3d3Vx00UVs3bp1rEOxEpkmBUlnSHpY0lpJnyiz/mhJd0haKelOSfOzjMfM6kNHRwerVq1i2bJlYx2KlcgsKUiaBFwDvBE4AThX0gklmy0FvhsRLwEuB67MKh4zqw/d3d10dnYSEXR2drq1UGeybCmcDKyNiK6I2AvcAJxdss0JwB3p41+UWW9mDaajo4O+vj4A+vr63FqoM1kmhXnA40XPN6TLit0PvCN9/DZghqSW0ooknSdphaQVmzdvziRYMxsdy5cvp7e3F4De3l6WL18+xhFZsSyTgsosi5LnFwGvl/Rb4PXARqB3QKGI6yJiUUQsam1tHflIzWzULF68mObmZCzO5uZmFi9ePMYRWbEsk8IG4Kii5/OBfPEGEZGPiLdHxEuBT6bLejKMyczGWHt7O01NyVdPU1MTS5YsGeOIrFiWSeFeYKGkYyVNAc4BbireQNIcSf0xXAx8M8N4zKwOtLS00NbWhiTa2tqYPXv2WIdkRTJLChHRC1wA3A6sAW6MiNWSLpd0VrrZqcDDkh4BjgSuyCoeM6sf7e3tnHjiiW4l1KFMJ9mJiFuBW0uWXVr0+IfAD7OMwczqT0tLC0uXLh31/Q413Sx4ylnwzGtmNkF5utnynBTMbMLwdLPD89hHZmZW4JaCmTW04vMIxR599FHguRZDqYl6fsFJwcwaWldXF4+sWcvcWQsOWD6pbwoAOzbtHVBm07b1oxJbPXJSMLOGN3fWAv761E9WvP3X7py4V8f7nIKZmRU4KZiZWYGTgpmZFficglmD8128Vg0nBbMJxHfx2nCcFMwanO/itWo4KdSJ0htsSpv5buKb2WhwUqhTbuab2VhwUqgTpa0AN/PNbCw4KZg1II/3Y7VyUjBrQF1dXax86GGaWp5/wPK+SG5NWrV54FTofd1PjEpsVt+cFMwyMtYXDzS1PJ+pb3lPxdvv+dl3M4vFxg8nBbNR4osHbDxwUjDLiC8esPHIYx+ZmVlBpklB0hmSHpa0VtInyqxfIOkXkn4raaWkN2UZj5mZDS2zpCBpEnAN8EbgBOBcSSeUbPYp4MaIeClwDvCVrOIxM7PhZXlO4WRgbUR0AUi6ATgbeLBomwAOTx/PBPIZxjMojyI5soZ6P/1emh28LK9syzIpzAMeL3q+ATilZJvLgE5JFwLTgdPLVSTpPOA8gAULFpTbZMT4CpGR5ffTRspYX+Jbz0by/yzLpKAyy6Lk+bnAtyPinyW9EvhXSSdGRN8BhSKuA64DWLRoUWkdB82jSI4sv582GibyAUeWV7ZlmRQ2AEcVPZ/PwO6hDwBnAETEryRNBeYAT2UYl42wwYZUAA+rYCOnXi7xbfTu0SyTwr3AQknHAhtJTiS3l2yzHjgN+LakFwNTgc0ZxpSZRv+gDKWrq4vVD61k+uyB6/al7brHnlo5YN2urRkHZpaxRmytZJYUIqJX0gXA7cAk4JsRsVrS5cCKiLgJ+HvgekkfJelael9EjHj30GhrxA/KcKbPhj8+s7oyq2/JJhazLDV692imdzRHxK3ArSXLLi16/CDw6ixjGC2N/kExs4nBw1zUyP3oZtaInBRq1NXVxe8fXMm8mQMvsmren/SA7d74wIB1G3vGfe+YmTUwJ4WDMG+muPDVh1RV5n/98tmMoqmNWzxmVsxJYYLr6upizZqVzHzewHX707tF8k8MvHKo5+mMAzOzMeGkYMx8Hryurboyd3VmE4uZjS0nBbMR5O44G++cFMxGUDI38hrUMvBOvv5bcB7Y/OTAdd2+k6+cWpNscYLN5/Ps7NnF1+68ouL9btq2jh0xvYaIxz8nBbMRppbZNL/5L6oq03vz7RlFM751dXXx0ENraZ199MCVMQWA7qf2HbB489Z1oxFaw3JSMKtjHhkUWmcfzbve9KmKt//BrZ894Hkul2OH9vLXp36y4jq+ducVzJg7peLtG4mTQo2SJmlUfYnphp7gMCXjAo5E09gmlok4hEqjqfdE76Qwhrq6unh4zUqOnDlwXdP+5Pe2/IGXgz7ZMwqBWd2ol5FBLTv1luidFGqUy+XYHd013bx2aC5XeH7kTFhyauV/hmV39la1PzOrL/We6DObo9nMzMYftxTMzIYw0e49cVIwG0H5fJ7Y3lP1JabRvZX8vv0ZRWUHIxn88vcsOGzg/PBT9iVXKD27fuAFJ+t3rs88tiw4KZiZDWPBYQu4+OUfr6rMlfddlVE02XJSMBtBuVyO7smTarp5Ldd65IjFkc/n6du+gz0/+27FZfq6nyC/b9eIxWDjk5OCWZ3x/SvPyefzbO/ZPeCGtKFs7l7Hs72HZhhVY3NSMKszyfhJD6GWOQPW9c9g/sDmLQcu7z7weS6XY+vkHqa+5T0V73fPz75LrrXMTTM2oTgpmNUhtcxhypvfXvH2e2/+UYbRjJ1cLschzfuqHuai5YjJGUbV2BoqKXR3d3PllVdyySWXMHv2wFEqzUZDdG+l9+bbiZ4d0Ltv8A2bJ6OZMwplGMFzCma1yjQpSDoD+DIwCfh6RHy+ZP2/AH+ePj0UOCIiZtW6v46ODlatWsWyZcu48MILa63GrGbHHXdc4XF+336eeeaZQbedNm3acyeXW488oKw1lvF0niizpCBpEnAN8AZgA3CvpJsi4sH+bSLio0XbXwi8tNb9dXd309nZSUTQ2dnJkiVL3FqwUddoJ3ptZHR1dbH2wUdYcPjcAeum9E4CYO+GHQcsX79906jEVirLlsLJwNqI6AKQdANwNvDgINufC3ym1p11dHTQ15dMKtzX1+fWgtkYKz46zufzw7eaisYEa8QrqRYcPpdLTvlgxdt/7p7rM4xmcFmOfTQPeLzo+YZ02QCSjgaOBZYPsv48SSskrdi8eXPZnS1fvpze3mSwuN7eXpYvL1uVmY2Srq4uHnxoLRs372Xn7v08u69v0J+du/ezcfNeNm7ey4MPrR20q8Wyl2VLQWWWxSDbngP8MCLK3ucfEdcB1wEsWrSobB2LFy/mtttuo7e3l+bmZhYvXlxLzGY2gma2LODVZ19SVZlf/vRzGUVjlciypbABOKro+XwgP8i25wDfO5idtbe309SUvJympiaWLFlyMNWZmU1IWbYU7gUWSjoW2Ejyxd9eupGkFwHPA351MDtraWmhra2NW265hba2Np9kHsdGemYqX6psVrnMWgoR0QtcANwOrAFujIjVki6XdFbRpucCN0TEYF1LFWtvb+fEE090K6HB7Nmz56Bmpyq+VNnMhpbpfQoRcStwa8myS0ueXzZS+2tpaWHp0qUjVZ2NkZGcmcqXKptVp6HuaB5v8vk8O7ZVN8Xmk9tg96CnZqzURL5Uua/7Cfb87Lv09WyF3r2Db9g8haaZs+nrfgLqcOyjzVvX8YNbP8u27U+wb9/AeQv6TZ58CLMOfz6bt66j5YgXjGKEjcVJwRpauUuVJ0JSOPDO6l0880zfoNtOm3ZIMhBe68y6u6u6OJ5neyfBM+UuakxMmzaJliMm03LEC+rudYwnTgpjKJfLsY0tLDm18j/Dsjt7mVV0k48NbaJeqlwPN37l83l6tu+q+hLTnu51aN90oD5ex0TjpGB1IauxYdrb2+ns7AR8qbJZJZwUrC4kcwishDnlLohLuj5Wbll14OItg3eJ9POlymMnl8sRk/fWdPNarnVKRlHZcMZlUqh1TJX+o8paj0qL67AMzGmi6a2Vz5jV95PdFW3X3t7OunXrxk0rIZ/PE9u3VzVHQnRvIb9viJPJZhUal0khGXFwDQtmzmb/rt30DTFm/f79wd6NT7K+Z2tJ+dUsmHnYgO2n7E/q2rtx3YB163t2jkD0Ntp8qbJZ5SpKCpLeBdwWETskfQp4GfDZiPhNptENYcHM2XzqtW0Vb//ZuztLyh/Gxa99SVX7vPLulVVtb1aLXC5H9+QpVc+8lmsdOH2nWbUqbSl8OiJ+IOk1wF8AS4FrgVMyi2wI+XyeXT09A77oh7KuZyvTVXa8PTPLSE/3en7508+xq+dJensHvyu9uXkq02ceWSgzr9X3GYyVSpNC/7fpmcC1EfFTSZdlE9L4sbEn+F+/HHgzzeZdyYgdrdMHXlO9sSdYWHYAcbPGcuC9EpN45pnBR9WZNm1S4eTyvFbfZzCWKk0KGyV9DTgduErSIWQ7wuqQcrkce2NS1d1HU3IjNwfuUB/a3vRk9aHzjh+wbuG8ocuaNQpfkDGyRmtKz0qTwruBM4ClEbFN0lyg/OU5o2R9z1Y+e3cnT+7awZ4hTjRPbZ7MkdNnsL5nKy+YN3JJYag3+WDG6hmP8vk827vhv747cF1f2sZsmjRw3f5eyPd6yA6zSiQXyDzMgpmtA9ZN2Z/0SuzduPWA5et7yk9KNpSKkkJE7Jb0FPAa4PdAb/p7TBQfaU/K76dpiEtSJ02bxpTckbxgnidGz8rhhx8+6GXB/cunTpk2cOWUpKyZVWbBzFY++aq/rHj7K/7f96veR6VXH30GWAS8CPgWMBn4N+DVVe9xBLhZWl+uvfbaQddNtFaT2XhXaffR24CXAr8BiIi8pBmZRWU1Ke5zLJ2YBnzj3XgS3VvK3rwWPT0AaObMAdvjS1JtBFSaFPZGREgKAEnTM4zJRsDBTEozFvL5PGzvq/guZQC29JHf23jnJIbq5nx0e5IUji9NAK1z3D2akXw+z64du7jyvquqKrd+x3qm58ffV2WlSeHG9OqjWZI+CLwfuD67sKwWxa0Ad9uMX76IwcZSpSeal0p6A7Cd5LzCpRHx80wjswkll8uxZcrWqsc+ys3xMOKWrVwux7O9z3Lxyz9eVbkr77uKQ3KHZBRVdoZNCpImAbdHxOmAE4GZWQMb9ga0iNgP7JZUf/P0mZnZiKr0nMIe4AFJPwd29S+MiA9nElWNuru7ufLKK7nkkks8br41hNK7WEvvXB2tK8qGurLNV7U1lkqHqrgF+DRwF3Bf0c+QJJ0h6WFJayV9YpBt3i3pQUmrJXVUGng5HR0drFq1imXLlh1MNWZ1a+rUqUydOnVMY9izZ8+4u7rNKlfpiebvSJoCvDBd9HBEDD62BIVzEdcAbwA2APdKuikiHizaZiFwMfDqiHha0hG1vAhIWgmdnZ1EBJ2dnSxZssStBRv36uUI3Fe2TRyV3tF8KvAd4DFAwFGS3hsRdw1R7GRgbUR0pXXcAJwNPFi0zQeBayLiaYCIeKraF9Cvo6ODvr5kesa+vj6WLVvGhRdeWGt1Zu4ysQmp0nMK/wy0RcTDAJJeCHwPePkQZeYBjxc938DA+RdemNb3S2AScFlE3FZakaTzgPMAFixYUHZny5cvp7e3F4De3l6WL18+qkmh+AtkrPp9x70tg9y81pPOxTyzacD2jNJNvO4usYmi0qQwuT8hAETEI5ImD1Nm4GQCEGX2vxA4FZgP3C3pxIjYdkChiOuA6wAWLVpUWgcAixcv5rbbbqO3t5fm5mYWL148THjZGes+3/FoyLt4e5Ike/yckqHI52Q7DLm7TGwiqjQprJD0DeBf0+dLGP5E8wbgqKLn84HSMQk2AL9Oz0/8QdLDJEni3grjKmhvb6ezM5mJrampadQnaXdL4OD4Ll7L0qZt6/nanVccsKx755MAtBw2cEj9TdvWM2PuxJz9rdKkcD7wN8CHSVoAdwFfGabMvcBCSccCG4FzgPaSbX4CnAt8W9Icku6k8rNIDKOlpYW2tjZuueUW2trafJLZ7CDUOqEL1F936WCtyace3QvAjLlTBqybMXfizv5WaVJoBr4cEV+EwpVFQ96/HRG9ki4Abic5X/DNiFgt6XJgRUTclK5rk/QgyZSfH4uI7hpfC+3t7axbt27UWwlmlRhP99F0dXWx8qHfM6ll/oB1fZH0HK/ePHAOjf3dGzKPrVqDJSi3QMurNCncQTIV5870+TSgE3jVUIUi4lbg1pJllxY9DuDv0p+D1tLSwtKlS0eiqnFjsCO68XY0NxKK34t8Pj/oxD8A06ZNI5dLxk0arfei+D6a8XBl3KSW+cw462+rKrPjpi9lFI2NlkqTwtSI6E8IRMROSZWPXGaZ6erq4qE1K2mZdeDySC/Y2bxp5YAy3dsGLGoIydHtamg5DHY/A737B912V+xjy+Z90L1z0G1Gku+jsfGi0qSwS9LLIuI3AJIWAYMfhtmoapkFbz6t3MVe5d18R9kLuBpDy2FMOvslFW++/6cDk2YWfB+NjReVDnPxEeAHku6WdBdwA3BBdmGZNZZy99GY1aNKWwrHkkzHuYBkas5XMPCeA7MRUS+DwI2kerqPxmwolbYUPh0R24FZJGMZXQcMPlu72Qiqh0HgDlZ7eztNTcm/21jcR2NWqUpbCv1n7M4EvhoRP5V0WTYh2UQ33loBlfB9NDZeVJoUNqZzNJ8OXCXpECpvZZgZvo/GxodKv9jfTXKj2RnpuESzgfIXwJtZWf330biVYPWs0vkUdgM/Knq+CdiUVVBmZjY23AVkZmYFTgpmZlbgpGBmZgWVXn1kVvfy+Txs31nd0BXdO8nvK53mw2ziclIwM8tYPp9n1/adfO6e6ysus277Jqbnd2QYVXlOCtYwcrkcWybvq3pAvFxrLsOozMYXJwUzG5fG06RFuVyOvX07uOSUD1Zc5nP3XM+U3IwMoyrPScHMBsjn8+zfvqvqSXP2d28gv296RlEdaLxNWjReOCnYiCse5XTURzjtTk809ww9yQ7Nk2DmtGSSndbswrFseNKi7DgpjLEne2DZnb0Dlj+dTgj2vMMGbj9rHHWBj+bopsUTref3VTAdZ2sOWgef2H0iy+VyPD35mZqm48y1Tssoqud40qLsOCmMoQO+xErmFH5mb/K475nkH6x/TuFZuZH9Esvn82zrgbs6qyu37Wmgr/ylnGM1ymkjjq5q5ZWbtMhJYWQ4KYyh4i+x0ollNm7cCMC8efOA8TmxjFlWap20aKiuTfD/GWScFCSdAXwZmAR8PSI+X7L+fcDVwMZ00f+OiK9nGVO9GqsPYi6Xg6YtvK6tunJ3dULu+eOoH8saSnt7O52dSfO21kmLqunaXL9zPVfed9WA5U/tfgqAIw49omyZhSysOq6xlllSkDQJuIZkprYNwL2SboqIB0s2/X5EeL5nM6tYrZMW1XLwNVR37d5H9wJwyIJDBqxbyMJxeb4qy5bCycDaiOgCkHQDcDZQmhTMrA7t795Q9pLUvp7NADTNHHjZ1v7uDdA6OkfHozVp0VCJpL/r6eqrr840Bkjviu7ZzhX/7/sVl1nX8xTTtaeq/WSZFOYBjxc93wCcUma7d0h6HfAI8NGIeLx0A0nnAecBLFiwIINQzazYUEe4j27fB8Dx5a4yah29o+P+SYtG25hecj0KskwKKrMsSp7/DPheRDwr6UPAd4ABZ4wi4jrgOoBFixaV1mFmI6xejo7r3Whecp3L5Vi7rfxYSE/u2gbAkdNnHbBcKDlvWIUsk8IG4Kii5/OBA65hjIjuoqfXAwPP5GQgaYbt5Mq7qxhNE1jXs5Pp8oiaZhPZWLUEhj638TQAU+YdeG7lBfNmV91yyzIp3AsslHQsydVF5wDtxRtImptO7QlwFrAmw3jMzMat0Wq9ZZYUIqJX0gXA7SSXpH4zIlZLuhxYERE3AR+WdBbQC2wF3pdVPMVyuRx7Yx8Xv7by0TQBrrx7JVOqbIqZmY0nmd6nEBG3AreWLLu06PHFwMVZxmBmZpXzdJxmZlbgpGBmZgUe+8gsVTr+VLFy4+QUa4Tr083AScGsoKuri5UPPQgtMweujGRuhpWbNw5c192TcWRjr9Fv2LLnOCkYPU+XHzp7Z3qfzGFlZgTseRpyz882rjHRMpPms15bVZHem+7OKJj6NJo3bNnom7BJYf0gN689uSuZx+DI6QNv4V/fs5MXzMs8tFFVemNL8bwOe59Nfj/TdOCcDpAkhPE42NdQ8vk8bO+p/ku+u4f8vnI38DcOtwQmjgmZFIa+MzBpGk+Zd/SAdS+Y13hfhKX/7MXdBJ7TwWzimZBJweO6DG4if+nncjm2TI6auo9yrb6p0RqDL0k1M7MCJwUzMytwUjAzswInBTMzK3BSMDOzggl59ZGZ2Whbv30Tn7vn+gHLn9yVzDV25PSWAdu/gDJ3jmbMScEaXnd3N1deeSWXXHIJs2fPHr6A2Qgrvr+wmUJ5AAAO0klEQVSp+AZRgGf29d8smgyl0n+T6AuYMSb3RTkpWMPr6Ohg1apVLFu2jAsvvHCsw7EJqPj+n9KBF+vtJlEnBWto3d3ddHZ2EhF0dnayZMkStxZsTNX7DaI+0WwNraOjg76+PgD6+vpYtmzZGEdkVt+cFKyhLV++nN7eXgB6e3tZvnz5GEdkVt+cFKyhLV68mObmpJe0ubmZxYsXj3FEZvUt03MKks4AvgxMAr4eEZ8fZLt3Aj8A/iwiVmQZUzlDTSACY3/ix2rX3t5OZ2cyWURTUxNLliwZ44jMDl7pyeqRnPgos5aCpEnANcAbgROAcyWdUGa7GcCHgXuyiqUaU6dO9SQiDaSlpYW2tjYk0dbW5pPM1pBG8nsry5bCycDaiOgCkHQDcDbwYMl2/wR8Abgow1iG5FZAY2tvb2fdunVuJVjDyPI7K8tzCvOAx4ueb0iXFUh6KXBURNycYRw2wbW0tLB06VK3EswqkGVSKDc/YRRWSk3AvwB/P2xF0nmSVkhasXnz5hEM0czMimWZFDYARxU9nw/ki57PAE4E7pT0GPAK4CZJi0oriojrImJRRCxqbW3NMGQzs4kty3MK9wILJR0LbATOAdr7V0ZEDzCn/7mkO4GLxuLqo/Esn8+zvQduviOG3zjVvQ32RX74Dc1swsmspRARvcAFwO3AGuDGiFgt6XJJZ2W1XzMzq12m9ylExK3ArSXLLh1k21OzjKVR5XI5JmsLbz6t3Cmc8m6+I2id64nmzWwg39FsZmYFTgpmZlbgpGBmZgVOCmZmVuCkYGZmBU4KZmZW4KRgZmYFTgpmZlbgpGBmZgWZ3tFso6N728Cxj3p2Jr9nHlZ++9a5oxDYeNTdQ+9Ndw9c3rMr+T1zetkytM4buNxsHHJSGOeOO+64ssu3p9Pztc49fsC61rmDl5vIhnpPHt2evJ/Hl/vyb53n99MahiIqH12zHixatChWrPBAqsPpn6v16quvHuNIGoPfTxvvJN0XEQOmJijlloLZIIonRx/JidHN6pmTglkFRmpSdLN656RgNgi3BGwi8iWpZmZW4KRgZmYFTgpmZlbgpGBmZgVOCmZmVuCkYDaBdHd3c9FFF7F169axDsXqVKZJQdIZkh6WtFbSJ8qs/5CkByT9TtJ/Sjohy3jMJrqOjg5WrVrFsmXLxjoUq1OZJQVJk4BrgDcCJwDnlvnS74iIP4mIk4AvAF/MKh6zia67u5vOzk4igs7OTrcWrKwsWwonA2sjoisi9gI3AGcXbxAR24ueTgfG10BMZuNIR0cHfX19APT19bm1YGVlmRTmAY8XPd+QLjuApL+R9ChJS+HD5SqSdJ6kFZJWbN68OZNgzRrd8uXL6e3tBaC3t5fly5ePcURWj7JMCiqzbEBLICKuiYjjgY8DnypXUURcFxGLImJRa2vrCIdpNjEsXryY5uZkZJvm5mYWL148xhFZPcoyKWwAjip6Ph/ID7H9DcBbM4zHbEJrb2+nqSn5l29qamLJkiVjHJHVoyyTwr3AQknHSpoCnAPcVLyBpIVFT88Efp9hPGYTWktLC21tbUiira2N2bNnj3VIVocyGyU1InolXQDcDkwCvhkRqyVdDqyIiJuACySdDuwDngbem1U8Zpa0FtatW+dWgg3KM681kHKTwhx//HPTcXpiGLOJyzOvTXCeFMbMauGk0EDcCjCzg+Wxj8zMrMBJwczMCpwUzMyswEnBzMwKnBTMzKzAScHMzAqcFMzMrMBJwczMCsbdMBeSNgPrhthkDrDlIHfTKHXUQwz1Ukc9xFAvddRDDPVSRz3EMFp1HB0Rw889EBEN9UMy2J7rqJMY6qWOeoihXuqohxjqpY56iKGe6ogIdx+ZmdlznBTMzKygEZPCda6jrmKolzrqIYZ6qaMeYqiXOuohhnqqY/ydaDYzs+w0YkvBzMxq5KRgZmYFDZUUJJ0h6WFJayV9ooby35T0lKRVNe7/KEm/kLRG0mpJH6mhjqmS/kvS/Wkd/1hLLGldkyT9VtLNNZZ/TNIDkn4nqeo5UCXNkvRDSQ+l78krqyz/onTf/T/bJf1tDXF8NH0vV0n6nqSqp6WT9JG0/OpKYyj3eZI0W9LPJf0+/f28Ksu/K42hT9KwUysOUsfV6d9kpaQfS5pVQx3/lJb/naROSblq6yhad5GkkDSnyhguk7Sx6PPxplpikHRh+r2xWtIXqq1D0veLYnhM0u9qqOMkSb/u/1+TdHKV5f9U0q/S/9efSTp8qBiGNBLXtdbDDzAJeBQ4DpgC3A+cUGUdrwNeBqyqMYa5wMvSxzOAR2qIQcBh6ePJwD3AK2qM5++ADuDmGss/Bsw5iL/Jd4D/kT6eAsw6yL/vEyQ34FRTbh7wB2Ba+vxG4H1V1nEisAo4lGS2wv8LLKzl8wR8AfhE+vgTwFVVln8x8CLgTmBRjTG0Ac3p46uGimGIOg4vevxh4KvV1pEuPwq4neSG1EE/a4PEcBlwURV/x3J1/Hn69zwkfX5ELa+jaP0/A5fWEEcn8Mb08ZuAO6ssfy/w+vTx+4F/quYzXvzTSC2Fk4G1EdEVEXuBG4Czq6kgIu4CttYaQERsiojfpI93AGtIvpSqqSMiYmf6dHL6U/XVAJLmA2cCX6+27EhIj1ReB3wDICL2RsS2g6jyNODRiBjqbvbBNAPTJDWTfLHnqyz/YuDXEbE7InqB/wDeNlyhQT5PZ5MkS9Lfb62mfESsiYiHKw18kDo609cB8Gtgfg11bC96Op1hPqND/G/9C/APB1G+YoPUcT7w+Yh4Nt3mqVrjkCTg3cD3aqgjgP6j+5kM8RkdpPyLgLvSxz8H3jFUDENppKQwD3i86PkGqvxCHkmSjgFeSnKkX23ZSWkT9Cng5xFRdR3Al0j+2fpqKNsvgE5J90k6r8qyxwGbgW+lXVhflzT9IGI5h2H+2cqJiI3AUmA9sAnoiYjOKqtZBbxOUoukQ0mO5I6qNpbUkRGxKY1tE3BEjfWMlPcD/6eWgpKukPQ4sAS4tIbyZwEbI+L+WvafuiDtxvrmUF1xQ3gh8FpJ90j6D0l/dhCxvBZ4MiJ+X0PZvwWuTt/PpcDFVZZfBZyVPn4XtX8+GyopqMyyMbneVtJhwL8Df1tyRFWRiNgfESeRHMGdLOnEKvf/ZuCpiLiv2n2XeHVEvAx4I/A3kl5XRdlmkibutRHxUmAXSXdJ1SRNIfnA/6CGss8jOTo/FsgB0yX9VTV1RMQakm6WnwO3kXRN9g5ZaByQ9EmS17GslvIR8cmIOCotf0GV+z4U+CQ1JJMi1wLHAyeRJPx/rqGOZuB5wCuAjwE3pkf8tTiXGg5cUucDH03fz4+StrCr8H6S/9H7SLqu99YYR0MlhQ0cmB3nU303wUGTNJkkISyLiB8dTF1pd8udwBlVFn01cJakx0i60RZL+rca9p9Pfz8F/Jiki65SG4ANRa2cH5IkiVq8EfhNRDxZQ9nTgT9ExOaI2Af8CHhVtZVExDci4mUR8TqSpnstR4MAT0qaC5D+HrK7IiuS3gu8GVgSaUf0Qeig+u6K40kS9f3p53Q+8BtJz6+0goh4Mj2A6gOup7rPZ78NwI/Sbtv/ImlZD3rCezBp1+Tbge/XEAPAe0k+m5Ac/FT1WiLioYhoi4iXkySmR2uMo6GSwr3AQknHpkeW5wA3jWYA6RHGN4A1EfHFGuto7b8aRNI0ki+1h6qpIyIujoj5EXEMyfuwPCKqOjqWNF3SjP7HJCcnK74qKyKeAB6X9KJ00WnAg9XEUORgjsDWA6+QdGj69zmN5FxPVSQdkf5eQPLPX2s8N5F8AZD+/mmN9dRM0hnAx4GzImJ3jXUsLHp6FtV/Rh+IiCMi4pj0c7qB5CKNJ6qIYW7R07dRxeezyE+AxWl9LyS5IKKW0UpPBx6KiA01lIXkAPb16ePFVHnQUfT5bAI+BXy1xjga5+qjeO6s/SMkWfKTNZT/HkkzdB/Jh/QDVZZ/DUmX1Urgd+nPm6qs4yXAb9M6VjHMlQwV1HcqNVx9RHJO4P70Z3WN7+dJwIr0tfwEeF4NdRwKdAMzD+I9+EeSL61VwL+SXmlSZR13kyS1+4HTav08AS3AHST/9HcAs6ss/7b08bPAk8DtNcSwluT8W/9ndLgrh8rV8e/p+7kS+Bkwr9o6StY/xtBXH5WL4V+BB9IYbgLm1vA6pgD/lr6W3wCLa3kdwLeBDx3E5+I1wH3p5+se4OVVlv8IyXffI8DnSUerqOXHw1yYmVlBI3UfmZnZQXJSMDOzAicFMzMrcFIwM7MCJwUzMytwUjCrQDoi50Vlluck/bDGOt833OiiZqPNScHsIEREPiLeWWPx95EMvWFWN5wUbMJK79q+RcncFask/WU6Hv6cdP0iSXcWFflTScuVzIXwwXSbY/rHtU8HMrxa0r3pIG1/XbSvf0jHur9f0uclvRNYBCxLx9Cfli5/MC27dPTeCbPnNI91AGZj6AwgHxFnAkiaSTLw3WBeQjJw2nTgt5JuKVn/AZJRWP9M0iHALyV1An9EMkT2KRGxW9LsiNgq6QKS+QBWSJpNcrfyH0VEaJiJb8yy4paCTWQPAKdLukrSayOiZ5jtfxoRz0TEFuAXDBy0rA14Tzrs+T0kQ1osJBkX51uRjjMUEeXG498O7AG+LuntQE1jEpkdLCcFm7Ai4hHg5STJ4UpJl5IMJd3/f1E6bWfpmDClzwVcGBEnpT/HRjJ3g8psWxpLL0mS+XeSVsVt1b4es5HgpGATVnrlz+6I+DeSiU1eRjIw28vTTUqHgz5byRzaLSQDDd5bsv524Px0+HQkvTAdYbYTeH86hwBpVxHADpKx7/vn4JgZEbeSTLhy0ki9TrNq+JyCTWR/QjLbVR/JiJPnA9OAb0i6hIGz5v0XcAuwgGQO3Hw6w15/K+DrwDEk8wKIZOa5t0bEbZJOAlZI2gvcClxCMrLmVyU9QzJnxE8lTSVpWXw0k1dsNgyPkmp2ECS9HPhiRLx+2I3NxgF3H5nVSNIikrHtvzzWsZiNFLcUzMyswC0FMzMrcFIwM7MCJwUzMytwUjAzswInBTMzK/j//bX4uRxD0b8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2e763726da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.boxplot(x='subjects',y='scores',data=result_scores)\n",
    "plt.title('grand average AUC = %.2f'%(result_scores['scores'].mean()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we see from above, most of the subjects have a average AUC score above chance. Each subject has 5 model scores, which represent our estimate of the model mx + b = y. Thus, the grand average AUC (0.67) is how much the between-subject variability we estimate from the current simulated experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAD0NJREFUeJzt3X2MHPV9x/HPp74AhdAC9ZGYh+ZsiaJA/iE5RUloqghDSmgapw2tjJTWbSpZaUVLqlatEVUSNf80fYjaP9ogl9CgFkESBxJEnQZKIFGrxu35gQfnoLYJTQyOfSQqeWgFRfn2j50z473Z29mZ2b3zl/dLWt/uzG9mvvebuY/n5nbn54gQAODk9yMrXQAAoBsEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBJTk9zY2rVrY2ZmZpKbBICT3u7du5+NiOlh7SYa6DMzM5qbm5vkJgHgpGf7v+q045ILACRBoANAEgQ6ACRBoANAEgQ6ACQxNNBt32r7mO3HStPOsX2/7QPF17PHWyYAYJg6Z+iflHR137Rtkh6IiIskPVC8BgCsoKGBHhFfkfSdvsmbJN1WPL9N0rs7rgsAMKKm19BfFRFHJKn4em53JQEAmhj7H0Vtb7U9Z3tuYWFh3JsDME720scobZdrj9aaBvpR2+skqfh6bFDDiNgeEbMRMTs9PfRWBACAhpoG+j2SthTPt0j6fDflAACaqvO2xTsk/Zuki20ftv0bkv5E0lW2D0i6qngNAFhBQ++2GBHXDZi1seNaAAAt8ElRAEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEhi6P3QAaDTsUAHrSuiu228THGGDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkESrQLf9u7b3237M9h22T+uqMADAaBoHuu3zJf2OpNmIeJ2kNZI2d1UYAGA0bS+5TEn6UdtTkk6X9Ez7kgAATTQO9Ih4WtKfS/qGpCOSnouI+7oqDAAwmjaXXM6WtEnSeknnSTrD9nsr2m21PWd7bmFhoXmlAIBltbnkcqWkr0fEQkT8n6S7JL2lv1FEbI+I2YiYnZ6ebrE5AMBy2gT6NyS9yfbpti1po6T5bsoCAIyqzTX0XZJ2SNoj6dFiXds7qgsAMKKpNgtHxIckfaijWgAALfBJUQBIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCRaDXABALJXugIUOEMHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIolWg2z7L9g7bj9uet/3mrgoDAIym7YhFfyXpnyLiWtunSDq9g5oAAA00DnTbPybpZyT9miRFxAuSXuimLADAqNpcctkgaUHS39nea/sW22d0VBcAYERtAn1K0uslfTwiLpP0A0nb+hvZ3mp7zvbcwsJCi80BmBj7xAdOCm0C/bCkwxGxq3i9Q72AP0FEbI+I2YiYnZ6ebrE5AMByGgd6RHxL0jdtX1xM2ijpa51UBQAYWdt3ufy2pNuLd7g8KenX25cEAGiiVaBHxD5Jsx3VAgBogU+KAkASBDoAJEGgA0ASBDoAJEGgA0ASBDoAJEGgA0ASBDoAJEGgA0ASBDoAJEGgA0ASBDoAJEGgA0ASBDoAJEGgA0ASbQe4AJBBedzQiGarUG+5kJdM7582tIaq9TYr62WFM3QASIJAB4AkCHQASIJAB4AkCHQASIJAB4AkCHQASIJAB4AkCHQASIJAB4AkCHQASIJAB4AkCHQASIJAB4AkCHQASIJAB4AkWge67TW299q+t4uCAADNdHGGfoOk+Q7WAwBooVWg275A0s9JuqWbcgAATbU9Q/9LSX8g6Ycd1AIAaKFxoNt+p6RjEbF7SLuttudszy0sLDTdHIA27OrHoLZ1V6s4Poiz9NIA0YvTyvP6X5eX7Z9eFkVL2b3Slqu77veYVJsz9Mslvcv2U5LulHSF7X/obxQR2yNiNiJmp6enW2wOALCcxoEeETdGxAURMSNps6QvRcR7O6sMADAS3ocOAElMdbGSiHhI0kNdrAsA0Axn6ACQBIEOAEkQ6ACQBIEOAEkQ6ACQBIEOAEkQ6ACQBIEOAEkQ6ACQBIEOAEkQ6ACQBIEOAEkQ6ACQBIEOAEkQ6ACQBIEOAEl0MsAFgFVmcXDkiOUHSh4yiLIVxwd+rjJo3uL0/kGjh22rDisUfmm95RqqBqyWTvw2o95mTkqcoQNAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACTRONBtX2j7QdvztvfbvqHLwgAAo2kzYtGLkn4vIvbYPlPSbtv3R8TXOqoNADCCxmfoEXEkIvYUz78naV7S+V0VBgAYTSfX0G3PSLpM0q6KeVttz9meW1hY6GJzQF720HE+K5uXlxuw/LAxO/vn12lfftRRHudzccny9PK0Qcv3zy9vv3990tLuWOyq8qOqXdUyq13rQLf9SkmflfSBiPhu//yI2B4RsxExOz093XZzAIABWgW67VeoF+a3R8Rd3ZQEAGiizbtcLOkTkuYj4mPdlQQAaKLNGfrlkn5F0hW29xWPazqqCwAwosZvW4yIf5GW+esFAGCi+KQoACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACTR+H7oEzdohNaoGJy2qm1VO+S1eAyMab/XWv0oNQwboXhxPUvaRXXbvtdWnDhwct+gyv3Pq14PEvLQQaLrrKtq8OfFaVX1Dlq2f5mqOvv7omqg6Krvq7976+7aScUPZ+gAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJtAp021fbfsL2QdvbuioKADC6xoFue42kv5b0DkmXSLrO9iVdFQYAGE2bM/Q3SjoYEU9GxAuS7pS0qZuyAACjahPo50v6Zun14WIaAGAFtBkkumq01iVDodreKmlr8fL7tp9osc2KKuoNZLtMu7WSnu2omq5RWzMv1Vb3+Gio1upPbFSv36pWXLkxl/5ddoVLhmCuWG7EfeoB/3alXLPXSnp26fqX2+7g77hqzihrPt6mN3tov3VwGL6mTqM2gX5Y0oWl1xdIeqa/UURsl7S9xXbGyvZcRMyudB1VqK0ZamuG2ppZTbW1ueTyH5Iusr3e9imSNku6p5uyAACjanyGHhEv2r5e0hclrZF0a0Ts76wyAMBI2lxyUUTslLSzo1pWyqq9HCRqa4ramqG2ZlZNbY5Y8ndMAMBJiI/+A0AS6QPd9i/Z3m/7h7ZnS9NnbP+v7X3F4+YBy59j+37bB4qvZ0+gtqts77b9aPH1igHLf9j206Xv4ZqualuuvmLejcUtH56w/bMDll9ve1fRd58q/njeuWLdi33wlO19A9o9VfTpPttz46ilYpu19tFK3EbD9p/Zftz2I7bvtn3WgHYT67dh/WD71GJ/HyyOrZlx1lPa7oW2H7Q9X/xM3FDR5m22nyvt6w9OorYTRETqh6TXSrpY0kOSZkvTZyQ9VmP5P5W0rXi+TdJHJ1DbZZLOK56/TtLTA5b/sKTfX4G+u0TSw5JOlbRe0iFJayqW/7SkzcXzmyX95gT2919I+uCAeU9JWjvh42/oPlLvTQWHJG2QdErRt5dMoLa3S5oqnn900LE9qX6r0w+SfkvSzcXzzZI+NaH9uE7S64vnZ0r6z4ra3ibp3kkeX/2P9GfoETEfEW0+zLRJ0m3F89skvbt9VT2DaouIvRGx+J7+/ZJOs31qV9uta5m+2yTpzoh4PiK+LumgereCOM62JV0haUcxqdO+q1Js85cl3THO7YzBitxGIyLui4gXi5dfVe+zJCupTj+Ufx53SNpY7PexiogjEbGneP49SfNahZ+MTx/oQ6y3vdf2l22/dUCbV0XEEam3UyWdO7nyJEnvkbQ3Ip4fMP/64lfmW7u8HDREnds+/ISk/y4FxiRuDfFWSUcj4sCA+SHpvuIy1tYBbcZh2D5aDbfReJ+kLwyYN6l+q9MPx9sUx9Zz6h1rE1Nc5rlM0q6K2W+2/bDtL9i+dJJ1SS3ftrha2P5nSa+umHVTRHx+wGJHJP1kRHzb9hskfc72pRHx3VVQ2+Kyl6r3q/DbBzT5uKSPqPcD9xH1Lje8bwL11bntQ61bQ9RVs87rtPzZ+eUR8YztcyXdb/vxiPhK05rq1KZ6+6jTvqpb22K/2b5J0ouSbh+wmrH0W1W5FdPGelyNyvYrJX1W0gcqsmKPpNdExPeLv5V8TtJFk6pNShLoEXFlg2Wel/R88Xy37UOSfkpS/x99jtpeFxFHbK+TdGzctUmS7Qsk3S3pVyPi0IB1Hy21/1tJ9466nYb11bntw7OSzrI9VZxJVd4aoq5hddqekvSLkt6wzDqeKb4es323er/itw6mun24zD6qdRuNJmr02xZJ75S0MYoLwRXrGEu/VajTD4ttDhf7/MclfWcMtSxh+xXqhfntEXFX//xywEfETtt/Y3ttREzsnkcv20sutqfdu6e7bG9Q73/SJyua3iNpS/F8i6Rlz6o7qu0sSf8o6caI+Ndl2q0rvfwFSY+Nu7bCPZI2F+84WK9e3/17uUERDg9KuraYNO6+u1LS4xFxuGqm7TNsn7n4XL3fesbeXzX30YrcRsP21ZL+UNK7IuJ/BrSZZL/V6Yfyz+O1kr406D+iLhXX6T8haT4iPjagzasXr+fbfqN6+frtcdd2gpX8i+wkHur9EB1W72z8qKQvFtPfo94fHB9W71elny8tc4uKd3Wod33uAUkHiq/nTKC2P5L0A0n7So9zK2r7e0mPSnpEvQN93ST6rph3k3rvSHhC0jtK03fqpXfobFAv6A9K+oykU8e4nz8p6f19086TtLNUy8PFY796lxwmcfxV7qNybcXra9R758ShCdZ2UL3r0YvH2M39tU2636r6QdIfq/efjiSdVhxLB4tja8OE+uqn1bu080ipv66R9P7F407S9aVM+aqkt0yitvKDT4oCQBIv20suAJANgQ4ASRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJAEgQ4ASfw/N5WUcHunyfgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2e763e19b00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for (_,k),(_,l) in zip(results.items(),labels.items()): \n",
    "    idx_positive = l == 1\n",
    "    idx_negative = l == -1\n",
    "    plt.hist(k[idx_positive],color='blue',bins=50)\n",
    "    plt.hist(k[idx_negative],color='red',bins=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown above, it the two distributions of measures for the positive responses (blue) and the negative responses (red). They overlap a lot, which suggests low signal-to-noise ratio for the data as a whole."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# average measures within each subject and fit a model for the averaged data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = []\n",
    "for (_,k),(_,l) in zip(results.items(),labels.items()):\n",
    "    idx_positive = l == 1\n",
    "    idx_negative = l == -1\n",
    "    data.append([k[idx_positive].mean(),1])\n",
    "    data.append([k[idx_negative].mean(),-1])\n",
    "data = np.array(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.     1.     1.     0.9375 1.    ]\n"
     ]
    }
   ],
   "source": [
    "model = linear_model.LogisticRegression(class_weight='balanced',random_state=12345,max_iter=int(1e5))\n",
    "cv = StratifiedKFold(n_splits=5,shuffle=True,random_state=12345)\n",
    "scores = cross_val_score(model,data[:,0].reshape(-1,1),data[:,1],cv=cv,scoring='roc_auc',)\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have a perfect score when we use the averaged measures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x2e761be3e80>"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAD8CAYAAACB3pQWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAE3NJREFUeJzt3X9w1PWdx/HXuyEaqlAsxl4xxqQ9lB8hhpJQ7nKciBzl0KmnQ8QW72A8G5meHTtTcbRVpzfXmXash3jVUhl1UExJEFA7tr3qRXuOFmqTElESsaIpjWGAQqWmGA/kfX/smuFHYL+72f1uPsnzMbOT3eS738/7vZu89ruf73e/MXcXACBMH8t3AQCAzBHiABAwQhwAAkaIA0DACHEACBghDgABI8QBIGCEOAAEjBAHgICNyMVKzz77bC8rK8vFqgFgSGptbf2juxene7+chHhZWZlaWlpysWoAGJLM7PeZ3I/pFAAIGCEOAAEjxAEgYDmZE+/PoUOH1NXVpd7e3riGHBKKiopUUlKiwsLCfJcCYBCKLcS7uro0atQolZWVycziGjZo7q59+/apq6tL5eXl+S4HwCAUKcTNrFPSe5I+lHTY3avTHai3t5cAT5OZaezYsdq7d2++SwEwSKWzJX6Ju/9xIIMR4OnjMQNwKuzYBICARQ1xl/SMmbWaWX02BjbL7iUOP/rRj/Too49KklavXq3u7u6+n11//fVqb2+PpxAASIo6nVLr7t1mdo6kZ83sdXd/4egFkuFeL0mlpaVZLnNwWLp0ad/11atXq6KiQuPGjZMkPfjgg/kqCxj0srmhxf92P1akLXF3705+3SPpCUnT+1lmlbtXu3t1cXHaH//Puc7OTk2YMEGLFy9WZWWlFixYoIMHD6q5uVlTp07VlClTdN111+mDDz6QJN16662aNGmSKisrdfPNN0uSvv3tb+vuu+/W+vXr1dLSokWLFqmqqkrvv/++Zs2apZaWFq1cuVK33HJL37irV6/W1772NUnSY489punTp6uqqko33HCDPvzww/gfCABDSsoQN7MzzGzUR9clzZX0Wq4Ly4Xt27ervr5eW7du1ejRo7V8+XItWbJETU1NevXVV3X48GGtXLlS+/fv1xNPPKFt27Zp69atuv32249Zz4IFC1RdXa2Ghga1tbVp5MiRx/xs48aNfbebmpq0cOFCdXR0qKmpSS+99JLa2tpUUFCghoaG2HoHMDRF2RL/lKQXzewVSS9L+qm7/3duy8qN8847T7W1tZKka6+9Vs3NzSovL9cFF1wgSVq8eLFeeOEFjR49WkVFRbr++uu1ceNGffzjH488RnFxsT7zmc9o8+bN2rdvn7Zv367a2lo1NzertbVVNTU1qqqqUnNzs956662c9Alg+Eg5J+7ub0m6KIZaci7q4XojRozQyy+/rObmZjU2Nuq+++7Tc889F3mchQsXat26dZowYYKuvPJKmZncXYsXL9Z3v/vdTMsHgBMMq0MMd+7cqU2bNkmS1q5dqzlz5qizs1NvvvmmJGnNmjW6+OKL1dPTowMHDmj+/PlasWKF2traTljXqFGj9N577/U7zlVXXaUnn3xSa9eu1cKFCyVJl156qdavX689e/ZIkvbv36/f/z6jM08CQJ/YPnZ/vHzsYZ44caIeeeQR3XDDDRo/frzuvfdezZgxQ3V1dTp8+LBqamq0dOlS7d+/X1dccYV6e3vl7rrnnntOWNeSJUu0dOlSjRw5su+F4SNnnXWWJk2apPb2dk2fntgHPGnSJH3nO9/R3LlzdeTIERUWFur+++/X+eefH0vvAIYm8xykaXV1tR//TyE6Ojo0ceLErI8VVWdnpy6//HK99lp4+2Tz/dgBA8UhhqmZWWsmpzQZVtMpADDUDJsQLysrC3IrHABOZdiEOAAMRYQ4AASMEAeAgBHiABCw/IV4iOeiPYV3331XP/zhD/tud3d3a8GCBXmsCMBwwJZ4lhwf4uPGjdP69evzWBGA4WDYhHhnZ6cmTpyor3zlK5o8ebLmzp2r999/Xzt27NC8efM0bdo0zZw5U6+//rokaceOHZoxY4Zqamp055136swzz5Qk9fT06NJLL9XnPvc5TZkyRU899ZSkxKlrd+zYoaqqKi1btkydnZ2qqKiQJH3+85/Xtm3b+mqZNWuWWltb9Ze//EXXXXedampqNHXq1L51AUBk7p71y7Rp0/x47e3tx34j8cGr7F1SePvtt72goMC3bNni7u51dXW+Zs0anz17tr/xxhvu7r5582a/5JJL3N39sssu8x//+Mfu7r5y5Uo/44wz3N390KFDfuDAAXd337t3r3/2s5/1I0eO+Ntvv+2TJ08+ZryPbi9fvtzvvPNOd3fv7u728ePHu7v7bbfd5mvWrHF39z/96U8+fvx47+npSf3YAYGJ8U89WJJaPIO8HTZb4pJUXl6uqqoqSdK0adPU2dmpX/3qV6qrq+v7Rw27du2SJG3atEl1dXWSpC9/+ct963B3ffOb31RlZaXmzJmjd955R7t37z7luFdffbUef/xxSdK6dev61vvMM8/oe9/7nqqqqjRr1iz19vZq586dWe8bwNCVtxNg5cPpp5/ed72goEC7d+/WmDFj+j1L4ck0NDRo7969am1tVWFhocrKytTb23vK+5x77rkaO3astm7dqqamJj3wwAOSEi8IGzZs0IUXXphZQwCGvWG1JX680aNHq7y8vG8r2d31yiuvSJJmzJihDRs2SJIaGxv77nPgwAGdc845Kiws1PPPP993OtlTnZpWkq655hrdddddOnDggKZMmSJJ+sIXvqAf/OAH8uQZfbZs2ZL9JgEMafkL8WzPimeooaFBDz30kC666CJNnjy5b+fiihUrtHz5ck2fPl27du3SJz7xCUnSokWL1NLS0vfv2SZMmCBJGjt2rGpra1VRUaFly5adMM6CBQvU2Nioq6++uu97d9xxhw4dOqTKykpVVFTojjvuyLgPAMPTsDkVbboOHjyokSNHyszU2NiotWvX5u3okdAeO+B4nIo2tUxPRTus5sTT0draqhtvvFHurjFjxujhhx/Od0kAcAJC/CRmzpzZNz8OAINVrHPiuZi6Gep4zACcSmwhXlRUpH379hFKaXB37du3T0VFRfkuBcAgFdt0SklJibq6urR37964hhwSioqKVFJSku8yAAxSsYV4YWGhysvL4xoOAIaFYf1hHwAIHSEOAAEjxAEgYIQ4AASMEAeAgBHiABAwQhwAAkaIA0DAIoe4mRWY2RYzezqXBQEAoktnS/wmSR25KgQAkL5IIW5mJZIuk/RgbssBAKQj6pb4Ckm3SDqSw1oAAGlKGeJmdrmkPe7emmK5ejNrMbOWAZ2p0CyzCwAMQ1G2xGslfdHMOiU1SpptZo8dv5C7r3L3anevLi4uznKZAID+pAxxd7/N3UvcvUzSNZKec/drc14ZACAljhMHgICl9U8h3P2Xkn6Zk0oAAGljSxwAAkaIA0DACHEACBghDgABI8QBIGCEOAAEjBAHgIAR4gAQMEIcAAJGiANAwAhxAAgYIQ4AASPEASBghDgABIwQB4CAEeIAEDBCHAACRogDQMAIcQAIGCEOAAEjxAEgYIQ4AASMEAeAgBHiABAwQhwAAkaIA0DACHEACBghDgABI8QBIGCEOAAEjBAHgIAR4gAQMEIcAAKWMsTNrMjMXjazV8xsm5n9exyFAQBSGxFhmQ8kzXb3HjMrlPSimf3c3TfnuDYAQAopQ9zdXVJP8mZh8uK5LAoAEE2kOXEzKzCzNkl7JD3r7r/ObVkAgCiiTKfI3T+UVGVmYyQ9YWYV7v7a0cuYWb2kekkqLS3NeqEY5MziHc/z82Yw221ms43BXBtyJ62jU9z9XUm/lDSvn5+tcvdqd68uLi7OUnkAgFOJcnRKcXILXGY2UtIcSa/nujAAQGpRplM+LekRMytQIvTXufvTuS0LABBFlKNTtkqaGkMtAIA08YlNAAgYIQ4AASPEASBghDgABIwQB4CAEeIAEDBCHAACRogDQMAIcQAIGCEOAAEjxAEgYIQ4AASMEAeAgBHiABAwQhwAAkaIA0DACHEACBghDgABI8QBIGCEOAAEjBAHgIAR4gAQMEIcAAJGiANAwAhxAAgYIQ4AASPEASBghDgABIwQB4CAEeIAEDBCHAACRogDQMBShriZnWdmz5tZh5ltM7Ob4igMAJDaiAjLHJb0DXf/rZmNktRqZs+6e3uOawMApJByS9zdd7n7b5PX35PUIencXBcGAEgtrTlxMyuTNFXSr3NRDAAgPVGmUyRJZnampA2Svu7uf+7n5/WS6iWptLQ0awXmlFlm93OPdzwMyGB+2KkNAxVpS9zMCpUI8AZ339jfMu6+yt2r3b26uLg4mzUCAE4iytEpJukhSR3uvjz3JQEAooqyJV4r6Z8lzTaztuRlfo7rAgBEkHJO3N1flMTsGAAMQnxiEwACRogDQMAIcQAIGCEOAAEjxAEgYIQ4AASMEAeAgBHiABAwQhwAAkaIA0DACHEACBghDgABI8QBIGCEOAAEjBAHgIAR4gAQMEIcAAJGiANAwAhxAAgYIQ4AASPEASBghDgABIwQB4CAEeIAEDBCHAACRogDQMAIcQAIGCEOAAEjxAEgYIQ4AASMEAeAgBHiABAwQhwAApYyxM3sYTPbY2avxVEQACC6KFviqyXNy3EdAIAMpAxxd39B0v4YagEApGlEtlZkZvWS6iWptLQ0W6sFgGOY5buCk3OPf8ys7dh091XuXu3u1cXFxdlaLQDgFDg6BQACRogDQMCiHGK4VtImSReaWZeZ/WvuywIARJFyx6a7fymOQgAA6WM6BQACRogDQMAIcQAIGCEOAAEjxAEgYIQ4AASMEAeAgBHiABAwQhwAAkaIA0DACHEACBghDgABI8QBIGCEOAAEjBAHgIAR4gAQMEIcAAJGiANAwAhxAAgYIQ4AASPEASBghDgABIwQB4CAEeIAEDBCHAACRogDQMAIcQAIGCEOAAEjxAEgYIQ4AASMEAeAgBHiABCwSCFuZvPMbLuZvWlmt+a6KABANClD3MwKJN0v6R8lTZL0JTOblOvCAACpRdkSny7pTXd/y93/T1KjpCtyWxYAIIooIX6upD8cdbsr+T0AQJ6NiLCM9fM9P2Ehs3pJ9cmbPWa2fSCFRXC2pD8eVUCOhztKdsc6to+wxddL7p/vofK8DJU+pAB6SePXsr9ezs9kzCgh3iXpvKNul0jqPn4hd18laVUmRWTCzFrcvTqu8XJlqPQh0ctgNFT6kOjlZKJMp/xG0ngzKzez0yRdI+kn2RgcADAwKbfE3f2wmd0o6ReSCiQ97O7bcl4ZACClKNMpcvefSfpZjmtJV2xTNzk2VPqQ6GUwGip9SPTSL3M/YR8lACAQfOweAAIWTIibWZ2ZbTOzI2Z2yr26ZlZgZlvM7Om46osqSh9mdp6ZPW9mHcllb4q7ziiiPichnLbBzD5pZs+a2e+SX886yXJ3JXvuMLP/Movz2NbU0uij1MyeSfbRbmZl8VaaWtReksuONrN3zOy+OGuMKkovZlZlZpuSv19bzWxhlHUHE+KSXpN0laQXIix7k6SO3JaTsSh9HJb0DXefKGmGpH8bpKc6SNlLQKdtuFVSs7uPl9ScvH0MM/tbSbWSKiVVSKqRdHGcRUaQso+kRyV9P/k7Nl3SnpjqS0fUXiTpPyT9byxVZSZKLwcl/Yu7T5Y0T9IKMxuTasXBhLi7d7h7yg8QmVmJpMskPZj7qtIXpQ933+Xuv01ef0+JF6RB9ynZiM9JKKdtuELSI8nrj0j6p36WcUlFkk6TdLqkQkm7Y6kuupR9JF9ER7j7s5Lk7j3ufjC+EiOL8pzIzKZJ+pSkZ2KqKxMpe3H3N9z9d8nr3Uq8sBanWnEwIZ6GFZJukXQk34VkQ/Jt7lRJv85vJRkL5bQNn3L3XVLiRVTSOccv4O6bJD0vaVfy8gt3H2zv+FL2IekCSe+a2cbktOP3k++YBpuUvZjZxyT9p6RlMdeWrijPSx8zm67ExsKOVCuOdIhhXMzsfyT9VT8/+pa7PxXh/pdL2uPurWY2K9v1RTXQPo5az5mSNkj6urv/OVv1pSMLvUQ6bUMcTtVLxPv/taSJSnxqWZKeNbO/d/coU3xZM9A+lPi7n6nExsFOSU2Slkh6KBv1pSMLvXxV0s/c/Q/53j2RhV4+Ws+nJa2RtNjdU26MDqoQd/c5A1xFraQvmtl8Jd72jjazx9z92oFXF10W+pCZFSoR4A3uvnHgVWUmC71EOm1DHE7Vi5ntNrNPu/uu5B9Rf3PEV0ra7O49yfv8XIl9FrGGeBb66JK0xd3fSt7nSSX6iD3Es9DL30iaaWZflXSmpNPMrMfdY9+BnoVeZGajJf1U0u3uvjnKuENqOsXdb3P3EncvU+L0AM/FHeDZkDzi4SFJHe6+PN/1DFAop234iaTFyeuLJfX3LmOnpIvNbETyRfZiDb4d6FH6+I2ks8zso/nW2ZLaY6gtXSl7cfdF7l6a/Ju/WdKj+QjwCFL2kvz7eEKJHh6PvGZ3D+KixFZQl6QPlNiZ9Ivk98cp8Xbq+OVnSXo633Vn0oekv1NiymGrpLbkZX6+a8/0OZE0X9IbSszvfSvfdZ+kl7FKHDXwu+TXTya/Xy3pweT1AkkPKBHc7ZKW57vuTPpI3v6H5O/Xq5JWSzot37Vn2stRyy+RdF++6x7A79e1kg4d9TffJqkq1br5xCYABGxITacAwHBDiANAwAhxAAgYIQ4AASPEASBghDgABIwQB4CAEeIAELD/BydpKStGsTqIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2e764033710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(data[data[:,1]==1,0],color='blue',label='positive')\n",
    "plt.hist(data[data[:,1]==-1,0],color='red',label='negative')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown above, after we average the measure within each subject, we do have better representation (signal) of the data distribution. There is a clear boundary between positive and negative responses. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# concatenate all the subjects but don't average over trials. Carefully cross validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = linear_model.LogisticRegression(class_weight='balanced',random_state=12345,max_iter=int(1e5))\n",
    "cv = StratifiedKFold(n_splits=5,shuffle=True,random_state=12345)\n",
    "data = []\n",
    "for n,measure_ in results.items():\n",
    "    temp = pd.Series(measure_)\n",
    "    temp = temp.to_frame()\n",
    "    temp['subject'] = n\n",
    "    temp['label'] = labels[n]\n",
    "    temp.columns = ['measures','subject','label']\n",
    "    data.append(temp)\n",
    "data = pd.concat(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.62526745 0.66935363 0.59995999 0.62929085 0.63665377] average = 0.63\n"
     ]
    }
   ],
   "source": [
    "scores = cross_val_score(model,data['measures'].values.reshape(-1,1),data['label'],\n",
    "                         groups=data['subject'],# this parameter controls the split by subject not the labels\n",
    "                         scoring='roc_auc',cv=cv,\n",
    "                        )\n",
    "print(scores,'average = %.2f'%(np.mean(scores)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown above, the average score is similar to method 1"
   ]
  },
  {
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}