Cochrane COVID-19 classifier.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# To recreate the model, calibration and results reported in the paper: Machine learning reduced workload with minimal risk of missing studies: development and evaluation of a COVID-19 research classifier for the Cochrane COVID-19 Study Register\n",
    "N.B. very slight variation is to be expected due to running in a different environment and the fact that RIS files are used here rather than reading text from the database.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import ris files\n",
    "\n",
    "from RISparser import readris\n",
    "\n",
    "TrainingIncludes = readris(open('1_Training_Included_20878.ris.txt', 'r', encoding=\"utf-8\"))\n",
    "TrainingExcludes = readris(open('2_Training_Excluded_38635.ris.txt', 'r', encoding=\"utf-8\"))\n",
    "CalibrationIncludes = readris(open('3_Calibration_Included_6005.ris.txt', 'r', encoding=\"utf-8\"))\n",
    "CalibrationExcludes = readris(open('4_Calibration_Excluded_10118.ris.txt', 'r', encoding=\"utf-8\"))\n",
    "EvaluationIncludes = readris(open('5_Evaluation_Included_2310.ris.txt', 'r', encoding=\"utf-8\"))\n",
    "EvaluationExcludes = readris(open('6_Evaluation_Excluded_2412.ris.txt', 'r', encoding=\"utf-8\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read data into dataframes (this can take a few minutes)\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "def parseData(incl, excl):\n",
    "    data = pd.DataFrame(columns=['title', 'abstract', 'label'])\n",
    "    for l in incl:\n",
    "        title = ''\n",
    "        abstract = ''\n",
    "        if 'abstract' in l:\n",
    "            abstract = l['abstract']\n",
    "        if 'primary_title' in l:\n",
    "            title = l['primary_title']\n",
    "        data = data.append({'title': title, 'abstract': abstract, 'label': 1}, ignore_index=True)\n",
    "    for l in excl:\n",
    "        title = ''\n",
    "        abstract = ''\n",
    "        if 'abstract' in l:\n",
    "            abstract = l['abstract']\n",
    "        if 'primary_title' in l:\n",
    "            title = l['primary_title']\n",
    "        data = data.append({'title': title, 'abstract': abstract, 'label': 0}, ignore_index=True)\n",
    "    data.abstract = data.abstract.fillna('none')\n",
    "    data.title = data.title.fillna('none')\n",
    "    data.reset_index(inplace=True)\n",
    "    return data\n",
    "\n",
    "TrainingDf = parseData(TrainingIncludes, TrainingExcludes)\n",
    "CalibrationDf = parseData(CalibrationIncludes, CalibrationExcludes)\n",
    "EvaluationDf = parseData(EvaluationIncludes, EvaluationExcludes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GridSearchCV(cv=StratifiedKFold(n_splits=2, random_state=None, shuffle=False),\n",
       "             error_score='raise-deprecating',\n",
       "             estimator=SGDClassifier(alpha=0.0001, average=False,\n",
       "                                     class_weight='balanced',\n",
       "                                     early_stopping=False, epsilon=0.1,\n",
       "                                     eta0=0.0, fit_intercept=True,\n",
       "                                     l1_ratio=0.15, learning_rate='optimal',\n",
       "                                     loss='log', max_iter=1000,\n",
       "                                     n_iter_no_change=5, n_jobs=None,\n",
       "                                     penalty='l2', power_t=0.5,\n",
       "                                     random_state=None, shuffle=True, tol=0.001,\n",
       "                                     validation_fraction=0.1, verbose=0,\n",
       "                                     warm_start=False),\n",
       "             iid='warn', n_jobs=None,\n",
       "             param_grid={'alpha': array([1.e-01, 1.e-02, 1.e-03, 1.e-04, 1.e-05, 1.e-06])},\n",
       "             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,\n",
       "             scoring='roc_auc', verbose=0)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# build model\n",
    "\n",
    "import numpy as np\n",
    "from sklearn.feature_extraction.text import TfidfVectorizer\n",
    "from sklearn.linear_model import SGDClassifier\n",
    "from sklearn.model_selection import GridSearchCV, GroupKFold, StratifiedKFold\n",
    "\n",
    "def getFeatures(data):\n",
    "    all_features = []\n",
    "    for index in range(len(data)):\n",
    "        title_tokens = data.title[index]\n",
    "        abstract_tokens = data.abstract[index]        \n",
    "        features = [\"t_%s\" % t for t in title_tokens.split(\" \") if not t in stopwords]      \n",
    "        features2 = [\"a_%s\" % t for t in abstract_tokens.split(\" \") if not t in stopwords]\n",
    "        features.extend(features2)\n",
    "        all_features.append(\" \".join(features))\n",
    "    return all_features\n",
    "\n",
    "stopwords = pd.read_csv('pubmed.stoplist.csv')\n",
    "ids = TrainingDf.index.tolist()\n",
    "labels = TrainingDf.label.astype(int)\n",
    "training_features = getFeatures(TrainingDf)\n",
    "\n",
    "vec = TfidfVectorizer(ngram_range=(1,3), max_features=75000, min_df=3, strip_accents='unicode')\n",
    "X = vec.fit_transform(training_features)\n",
    "\n",
    "sgd = SGDClassifier(class_weight=\"balanced\", loss=\"log\")\n",
    "parameters = {'alpha':10.0**-np.arange(1,7)}\n",
    "clf = GridSearchCV(sgd, parameters, scoring=\"roc_auc\", cv=StratifiedKFold(n_splits=2))\n",
    "clf.fit(X, labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# score calibration and test data\n",
    "\n",
    "#vec2 = TfidfVectorizer(vocabulary=vec.get_feature_names(), ngram_range=(1,3), max_features=75000, min_df=3)\n",
    "\n",
    "calibration_features = getFeatures(CalibrationDf)\n",
    "evaluation_features = getFeatures(EvaluationDf)\n",
    "\n",
    "X_calibration = vec.transform(calibration_features)\n",
    "X_evaluation = vec.transform(evaluation_features)\n",
    "\n",
    "calibration_preds = clf.predict_proba(X_calibration)\n",
    "calibration_probabilities = calibration_preds[:,1]\n",
    "\n",
    "evaluation_preds = clf.predict_proba(X_evaluation)\n",
    "evaluation_probabilities = evaluation_preds[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAG5CAYAAAADAnkFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xm8HUWd9/HPlyQkLAkBBBJBNAgiIyAoCigCKkFEx5HhmXFlUYRxQURZHHCQzLiBCiKCG/gQFjceF0ZQWSWMCuiAomxuGJQthDUJS8JWzx/dJzk5uTe599w1nc/79erXvae6uru66yy/U9VVJ6UUJEmStPJbbaQLIEmSpMFhYCdJktQQBnaSJEkNYWAnSZLUEAZ2kiRJDWFgJ0mS1BAGdpIkSQ1hYCdJktQQBnaSJEkNYWCnQZHkwCSll+XzI12+JkpyZ5IL+pBvj7oedunn/g9Nsn/3JVy1JBmf5OtJ7knyTJLrhum4y9Rvkk8meaojX5+eL4NUprWTzEiyaw/r3lOXd5PhKMtIaauXNw/iPlvXbrs+5D0vyV860u5Mcmbb483r/b2zLW2Xuu4m9bDPXyS5fKDnoaE1dqQLoMZ5F/CHjrS7R6IgWuzXwM7Azf3c7lDgTuCcQS9RMx0KHAy8H7gBWDCCZfkqcOEIHn9t4HjgKeB/Otb9N3ATMHe4C7WKOR6YuII8d1C9N7QHgLvU254JzO/Ifwjg75COcgZ2Gmw3lVL63FKRJMD4UsrCISzTKq2UMh+4dqTL0V9J1iylPDbS5eiHrYEFpZSvjHRBSil3UgXlgyLJ6sAzpZSnVph5BUop9wH3DbxUw2swr8FwKKXc1oc8i+jHe0Mp5ZYBFUrDwq5YDZskY+tm/1OSvD/JH4AngHfW68cn+XiSPyZZlGRukm8keVbHflZP8vkkc5I8luTnSXbooZthme6oOr3HrqAkb0tybb3PBUkuTvLijjznJXk4yQvq9Y8m+XuSz9Zv/O15J9RdGn+oz+eBJD9LslO9/qokN/VQvtWSzE7y3328rnsn+W2Sx5PcmuSAjvU9ddVtnuT8uttwUX0tL0+yTb3+TmBL4LVtXep/adv+uUm+meS+evtbkny4DtTbj71pkh/U1/OhJOcm2amH7p/WdX1xksuSLAAurde9LsmP6vp9PMmfk3w5yfodx/pkvd8XJflekvn1Nf98kjFJtkpyaZJHktye5Mg+Xt81kpxYb/NEXY4vJVmnXj82SQEOBCa2Xa93rmC/eye5Msm8+jl3S5Kj29a/PMl36+M+Xv/9ZpJN+1DmHp/79br/k+TGut5uS/KBjvWt58vbk3whyd3AQuB5STZK8pX6efZIknuTXJHklW3bbw7cUz/8RNv1OLNe39vr7+Akv0+yMMmD9fNmy448fX799XLudya5YCDXoF6/Tf2cfLgu72+T7NfLYddI9Z53b12PV2bZ95X+1vV6Sc6ur9OjSf47yfN6uFZ/6XnzxXmW6opN8kngM/XqO9rqbpd6/TJdsen7+/Yeqd7zHqzP8W+pXqcTlldG9Z8tdhpsY5Is9bzq4Rvu/6HqhpkB3Avcm2QMcBGwE3Ai1bfI5wH/BbwsycvbWvX+L/A24HPAFcC2wAXAMveE9FWSj9flORP4BDAeOBr4RZIdSil/bMs+nqo76Yy6DLsD/wE8DHy63t844BKqbo4vAFcCq9fn95z6/L4IfD/J7qWUWW37f2N97v/Wh6K/BPgscALVNf03YGaSP5dSrl7Odj8FngGOAv4OPAt4JbBuvf4fgR/W+zysTltYn9tGwDVUXww/BvwN+CfgZGBaK3+StevzXofqWt4GvAH4di9lmkB1Xb9C9eEypk7fHPgl1fWeVx/jCOB/kry4h+fX96m6j78KvA44EhhX/39afb32Az6X5E+llB/1dpGSrAb8CNiNqm5/CWxH9VzZKckrSylPJNm5TnslML3evNcP1ST/VpfvZ8B7qa7zlsBWbdmmAbdSXa8HgalU3b2/TrJVKeWh3va/HDu0lf9equtwWpJxpZRTOvJ+tj7fQ+rHD9RleJqqq+5eqq6+fanqYvdSys+puvfeAPwY+Bows96+167XJMdRvda/Cfw71fNxBnBt/fprb31a4etvKK9Bkn8ArgbmUNXHQ8D+wDlJNiilnNyxjxOB64B3U72+/hO4qn7u/q3O09+6nglcDLwd2BT4FDCr3ue8PlyD3ny1LuP7gTexpGW1x9s4+vq+neT5db4rqb4AzQc2Bvaiem3aYzOYSikuLgNeqF6spZdlbJ1nbP34QWCdju3fWa97U0f6jnX6wfXjrevHn+3It3+dfmZb2ieBp3oo63vqvJvUj59HdS/QyR35JlK98X+zLe28ett9OvJeQtUN3Xr8rjrfgcu5ZmOA24HvdaRfCvwRyAqu+Z3Ao8DGbWlrUH3AndaWtkddll3qxxvVjz+wgv3/Abi8h/TPUQWFL+1I/3qd/vz68WH1cfboyHdmnf7OHq7rfisoU+rn0WZ1/r076rsAh3Vsc2Od/o9taatTBSrfXcHx3lBv++GO9LfX6e/qOIeH+/BamUT1wXbliuq4Y7uxVPeuPQa8v7f67e25Xz9fngG27kj/GVVwMqFjf1f0oUxj6nLNAs5vS59S7+M/+vD6W5/qg/2/O/I9D1gEnN3f198KXjMDugbA/wMep+1111aGBcDEjn38qr2e6+fuk8BXuqjr1rU7vyP/rnX6Rzuu1V96OP/298jNWfa1+O/t9dOx/S9oe0+g7+/bb6kfv6ivz3eX7he7YjXY9gde1r6UZVtULi/Lfqt8I9UH7U9SdW2NrVv+rgfup/pWDvDq+u83O7b/DtUbdjf2ovqAOqfj2I8DP287dsvTVK0R7X4PPLft8eupgq6zeztoKeVp4HTgn5JsDJDkBVQfCKeX+h1xBX5TSrmrbZ+PA3/uKEun+6gCyn9PcniS7eqWqb56DXBjKeX6jvSZVIFXq452Ax4qpXSOouutxQ7gB50Jqbr/vp6qe/gpqg/FVgvOVp35qVoG2t1KVWcXtxJKKU/U+1jedYLqXGFJq1PLd6iCkdeuYPue7EL1peHLy6vjJBOTfK7uKmyd9wKq4L2n8+6L35VSOrv/vwVMpmrFavf9Xsr1gbrrcSFL6mO3AZTpFVStcDPbE0sptwNXsew17svrb3kGeg1eA1za/rqrnU0VjO3Yue/2ei6l/JUq2Gu9Trqp66Xe/0op/wPc1b7PYdLX9+3fUp3TmUn2TzJtmMu5SjGw02C7tZRyXfvSQ557ekjbiOqb+5M9LM+qF+o8UHWDLFZ/UHfTNdU6Nix582lf9m07dssj9fHaLaJ6E27ZALirD8HZmVT3Gba6XQ+l+pY+s49lf6CHtM6yLKWU8gzVB8BlwDFU5z23vg9o7T4cc316rsO729a3/t7bQ76e0gDml1IebU+ou3oup+oWOoHqQ/7lVMER9HyeD3Y8foKqzp7sIX1F9/esDywqHV1h9TW8lyXn2h8b1H9XNLjhu8D7qLozX0d13i+jOr9e63cF5iwnrfNclqnjVPcAnkbVFbkvVRfcy6ieS92WqXXc3p5TneXqy+tveQZ0Dai6Kvvy/F/R8drz9beu+7LP4dCn9+1Syp+oblF4gOpWi78m+UuSQ4e5vKsE77HTSOgp2Lmf6oPyjb1s0xp23wpkptAWIKS6cXrdjm0WAqslGdvRatgZqN1f/30z1bfevpR3Re4DdkiS5QV3pZSHkpwLHJLkFOAA4NxSjWQdMnVryLsBUt2g/haq+6bGUgWXy9O616rTs+u/97fl27aHfFN6K1YPaS+m6n5/ZyllcStFkheuoIyD5QFgfJJ124O7uoVzI6oW3f5q3bfU6zxuqQaGvJ6qK/OzbelrsOzzvD96uvattM4vCT3VxzupWtw7Bxt0fX9r23F7e07d30P6QAz0GjxE357/KzreA9B1Xfe2z2UGYw2xvr5vU0q5iurewrFU9zl+CPhSkjmllO8NeUlXIbbYabS4iPrer84Wv3r5U53vyvrvOzq2fyvLPp9vp+oa3KYj/R87Hl9M1b3z/F6O3dnl2Bc/peqW6W2kXLtTqc79e1T3X53WxfG6Vkr5Yynlv4BbqAZjtPTWCnIFsE2SzqBtf6ru8FYdXQWsm2R6R7639qd4bWVp15eBJYPhivpv5wjXf6Vq7buC/vsFVTfbe5eTp3VbQed5H0L1nO7Wi5O8qCPt7VT3Zd7Qh+1LZ5mSbE/VutSulacvrWi/rPMvdY3rEaG70d01Xp6BXoMrgD3qQUTt9gceoZo3snPfiyXZjKq7dlad1E1dL/X+l2oi6I3b9jkQ/am7vr5vL1ZKeaqUci1LvkC+pDOPBsYWO40W36R6A7ykbrm6jur+nY2p7mn5XinlwlLKTUm+AxyR5Gmqm563AT7CshPCXkT1Zn1WkuOp3kDfTce37VLKbUn+Ezgh1VQNl9TbTaHqEplXBz79cR7VgJIzkmxFFeSMoeq6+n0p5f+1Hf+WJJdRdVXMKqX0dyLhfknyEqoRrN+juh/vSar7+l5ENSK45UZg3yT/CswGHq/vTTqJ6kP4p/Vo4juoguV/A06t7yGCavTyh4BvJ/kP4K/A3iy5Z6ov90TeTBWgf7b+pv8wVbdsN/e2deOnVF3Bn08ymWo08IupRjZeR3VvVr+UUuYnOQr4al3vZ1KNGN2C6qb+w+qW3KuBjyZ5kGrk8qtZMqKwW3cDP65fD3Opvni8Gjii9G0uyYuAY+p6/wXwQuDjVHXUfo4PJbkL2CfJLKpWrvvKklGg7XkfTPIp4L+SnEXVLbkB1ajVx1j6OTkYBnoNZlC1sM1K8gmq5+R+VF2oR5RSOt+Hpib5AfANqvv4/ovqvE6Exdeqv3W9U5KvU90D+FyqwTJ3UI1qHagb67+HJzmP6v3hD6WUR3rI26f37VTTybyK6vX0d6qg8aB6H/6SxWAb6dEbLs1YWDIqdofl5GmNij2ll/XjqKbf+B3VwIUFVDe+f4V6pGWdbzxVYDK3zvdLqgBsqRFfdd6dqD6MH6F64zuOJbOnb9KRdx+qb7zzqLpxZwPnA69uy9PjyEd6HoU4gepD6c9U93PdT/Um9vIetj+oLtO+/bjmdwIX9JDeOXKtc1TsFKp7+P5QX5f5VC0VhwGrtW03jWqE7oJ6+7+0rXseVVBzf31uf6AKrtNRlk2ppk15pL6u51N123SOaO11RClVwHlZXc4HqQYuPJeOUZcsGRU7uWP73ursF8ANfbjOa1JNe/G3+lzvompV7RzZ3adRsW3530j1qwyPUA20uZkqMGitfw7VB/eD9bn/hOpG+s6Rjf0ZFXsBVWvjzVQtM39l2VHErf29uYcyj6cK7O+ieu1dRxXU9zQCc3r9vFpI24h1OkbFtuU/hGoQxCKqQPCHwAv7WJc9joDv7TUzkGtQr9+W6pc9Wu8Vv6VjRHfbPt4GfInq/Woh1XvM9h15+1rXrWv32ta1qJ87FwKb9XCtuhkVG6qg826qnoz2946l3lvqtBW+b1MNkPkh1WtoIdX7xs9oew9wGbwl9UWXVnr1qMmLSynvGemy9FeqyYi3p3pzXilmtu9W3dpzPNUHe083oUtDon6PuK6UMmi/3yqNNnbFSiMkyXiq+0t2pmr1OKxpQV2Sw6i+9f+Rau641wIfpJqbzKBOkgaZgZ00cp5DNW3EfODL9dI0j1PdZ/dcqm68v1H9OsCnRrJQktRUdsVKkiQ1hNOdSJIkNcSIB3ZJNk5yXpIHkjyW5IYkL21bnyQzktyd5PEkszrnIEqybpJzk8yrl3PrqQna82yT5Kp6H3cl+XiSgcwHJUmSNKqMaGCXZF2qqSqepJoX6B+AI6iGcLccTTWNwqFUk2DOAS5LMrEtz7eofuNvr3rZDji37TiTqKZLuLvexweBI+v9SpIkNcKI3mOX5ATglaWUV/WyPlTB2CmllBPrtPFUP2Hy0VLK1+rJX28Bdiql/KrO05q77IWllD8meR/wGWCjUsqiOs+/UwV4m5QVXIS6HM9m2QlwJUmShspE4O4VxSntRnpU7JuoZqz+f1Q/HXMX8OVSyhn1+mlUk6le2tqglLIoyVVUEx5+jWqqiHmtoK7Oc22SeXWeP9Z5rmoFdbVLqIK951FNRLtYHTyOb0uaSjUBqyRJ0nDahJ5/x7xHIx3YbQa8j+pXBD5N9esBpyZZVEo5hyU/dHxvx3b3Uk2fQJ1nbg/7ntu2/RQ6fvKmbZ9T6AjsgGOoJlBdyh133MGkSQP5rWtJkqQVmz9/Ps95znOgn72FIx3YrUY1C/ix9ePf1gMj3gec05avswkyHWk9NVGuKE96SYeqJe/ktscTgTsnTZpkYCdJkkatkR4Vew/V/XHtbqX6fUmoBkrAkpa3lg1Z0uI2B9ioh31v0JGnp33Asq2BlFIWlVLmtxa8t06SJK0ERjqw+yWwZUfaC6hmp4eqi3QO1Y9JA5Bkdar78a6uk64B1kny8rY8OwLrdOTZtd62ZU+qgRm3D8aJSJIkjbSRDuy+AOyU5Ngkmyd5O3AIcDpAPQrkFODYJPsk2RqYCTxGNcUJpZRbgYuBM5LsVI+IPQO4qJTyx/o43wIWATOTbJ1kH+BY4OT+jDSRJEkazUb0HrtSyv/WQdZngI9TtdAdXkr5Zlu2zwJrUP2O5rrAr4A9Synt3aPvAE5lyejZH1HNe9c6zrwk06kCxuuAh6juoWu/j06SpEZ6+umnefLJJ0e6GGozbtw4xowZM+j79bdi+6Ce4HjevHnzHDwhSVpplFKYM2cODz/88Ioza9hNnjyZKVOm0NMPYc2fP5911lkHYJ36fv8+GelRsZIkaYi0groNN9yQNddcs8cAQsOvlMJjjz3G3LnVbG1Tp04dtH0b2EmS1EBPP/304qBu/fXXH+niqMMaa6wBwNy5c9lwww0HrVt2pAdPSJKkIdC6p27NNdcc4ZKoN626Gcz7Hw3sJElqMLtfR6+hqBsDO0mSpIYwsJMkSY01c+ZMJk+ePOD97L777hx++OGDUKKh5eAJSZJWMTNmzRje4+3ev+MdeOCBPPzww1xwwQVDU6AGs8VOkiSpIQzsJEnSqLX77rtz2GGHcfTRR7PeeusxZcoUZsyYsVSehx9+mEMOOYSNNtqICRMmsPXWW3PRRRf1uL8DDzyQN7/5zUulHX744ey+++6LHz/66KPsv//+rL322kydOpWTTjppmf088cQTHH300Wy88castdZa7LjjjsyaNWugpztgdsVKkqRR7eyzz+YjH/kIv/rVr7jmmms48MADeeUrX8n06dN55plneP3rX8+CBQs477zzeP7zn88tt9wyoHnhjjrqKK688kp++MMfMmXKFI499liuv/56tttuu8V53vWud3H77bfzne98h2c/+9n88Ic/ZK+99uLGG29kiy22GIzT7oqB3WjS8Q1kpdSEc5AkjSrbbrstxx9/PABbbLEFp512GldccQXTp0/n8ssv59e//jW33norL3jBCwDYbLPNuj7WI488wje+8Q3OOeccpk+fDlSB5SabbLI4z2233ca3v/1t7rzzTp797GcDcOSRR3LxxRdz1lln8elPf7rr4w+UgZ0kSRrVtt1226UeT506dfHPcd1www1ssskmi4O6gbrtttt44okn2HnnnRenrbfeemy55ZaLH//mN7+hlLLMMRctWjTiv/JhYCdJkka1cePGLfU4Cc888wyw5Ke5+mq11VajlLJUWvsvP3Su68kzzzzDmDFjuP7665fp8l177bX7VZ7B5uAJSZK00tp222258847+dOf/tSn/BtssAH33HPPUmk33HDD4v8333xzxo0bx7XXXrs47aGHHlpq/9tvvz1PP/00c+fOZfPNN19qmTJlygDPaGAM7CRJ0kprt912Y9ddd2XfffflsssuY/bs2fz0pz/l4osv7jH/a17zGq677jrOOecc/vznP3P88cdz0003LV6/9tprc9BBB3HUUUdxxRVXcNNNN3HggQey2mpLQqYXvOAFvOMd72D//ffnBz/4AbNnz+Z///d/OfHEE/nJT34y5Oe8PHbFSpK0iunvhMGj3fe//32OPPJI3va2t/Hoo4+y+eabc8IJJ/SY93Wvex3HHXccRx99NAsXLuTd7343+++/PzfeeOPiPJ/73Od45JFHeNOb3sTEiRM54ogjmDdv3lL7Oeuss/jkJz/JEUccwV133cX666/PzjvvzN577z2k57oi6Utf8qouySRg3rx585g0adLQHagJI0qbcA6S1AALFy5k9uzZTJs2jQkTJox0cdSD5dXR/PnzWWeddQDWKaXM7+s+7YqVJElqCAM7SZKkhjCwkyRJaggDO0mSpIYwsJMkSWoIAztJkqSGMLCTJElqCAM7SZKkhjCwkyRJaggDO0mStMq4/fbbScINN9wwoP3MmDGD7bbbbpBKNXj8rVhJklY1w/3zj/083oEHHsjZZ5+9TPrrXvc6Lr744kEqVDMZ2EmSpFFnr7324qyzzloqbfz48SNUmpWHXbGSJGnUGT9+PFOmTFlqWXfddZk1axarr746P//5zxfnPemkk3jWs57FPffcA8AzzzzDiSeeyOabb8748ePZdNNN+dSnPtXjcWbOnMnkyZOXSrvgggtIslTaCSecwEYbbcTEiRM56KCDWLhw4TL7Ouuss9hqq62YMGECL3zhC/nyl7880MvQb7bYSZKklcbuu+/O4Ycfzn777cfvfvc7br/9dj72sY/x7W9/m6lTpwJwzDHHcMYZZ/CFL3yBXXbZhXvuuYc//OEPXR/z/PPP5/jjj+f000/nVa96Feeeey6nnnoqm2222eI8Z5xxBscffzynnXYa22+/Pb/97W85+OCDWWuttTjggAMGfN59ZWAnSZJGnYsuuoi11157qbSPfvSjHHfccXzyk5/k8ssv55BDDuHmm29mv/32Y5999gFgwYIFfPGLX+S0005bHFA9//nPZ5dddum6LKeccgrvfve7ec973gOw+PjtrXaf+MQnOOmkk/jnf/5nAKZNm8Ytt9zC1772NQM7SZK0anv1q1/NV77ylaXS1ltvPQBWX311zjvvPLbddlue+9zncsoppyzOc+utt7Jo0SJe+9rXDlpZbr31Vt773vculbbzzjtz5ZVXAnDfffdxxx13cNBBB3HwwQcvzvPUU0+xzjrrDFo5+sLATpIkjTprrbUWm2++ea/rr776agAefPBBHnzwQdZaay0A1lhjjX4dZ7XVVqOUslTak08+2a99PPPMM0DVHbvjjjsutW7MmDH92tdAOXhCkiStVG677TY+/OEPc8YZZ7DTTjux//77Lw6utthiC9ZYYw2uuOKKPu1rgw02YMGCBTz66KOL0zrnuNtqq6249tprl0prf7zRRhux8cYb89e//pXNN998qWXatGndnmZXbLGTJEmjzqJFi5gzZ85SaWPHjmXddddlv/32Y8899+Rd73oXr3/969lmm2046aSTOOqoo5gwYQIf/ehHOfroo1l99dV55StfyX333cfNN9/MQQcdtMxxdtxxR9Zcc02OPfZYPvjBD/LrX/+amTNnLpXnQx/6EAcccAA77LADu+yyC9/85je5+eablxo8MWPGDA477DAmTZrE61//ehYtWsR1113HQw89xEc+8pEhuUY9MbCTJGlVM9wTFHfh4osvXjzKtWXLLbfk7W9/O7fffjsXXnghAFOmTOHMM8/kX//1X5k+fTrbbbcdxx13HGPHjuXjH/84d999N1OnTl3mHrmW9dZbj/POO4+jjjqKr3/96+yxxx7MmDGDQw45ZHGet7zlLdx222189KMfZeHChey77768733v45JLLlmc5z3veQ9rrrkmn/vc5zj66KNZa6212GabbTj88MOH4Or0Lp39ylpWkknAvHnz5jFp0qShO9BK8EJboSacgyQ1wMKFC5k9ezbTpk1jwoQJI10c9WB5dTR//vzWwIt1Sinz+7pP77GTJElqCAM7SZKkhjCwkyRJaggDO0mSpIYwsJMkqcEcJDl6DUXdGNhJktRA48aNA+Cxxx4b4ZKoN626adXVYHAeO0mSGmjMmDFMnjyZuXPnArDmmmuSZIRLJaha6h577DHmzp3L5MmTB/VnxwzsJElqqClTpgAsDu40ukyePHlxHQ0WAztJkhoqCVOnTmXDDTfs9w/ba2iNGzduUFvqWgzsJElquDFjxgxJEKHRx8ETkiRJDWFgJ0mS1BAGdpIkSQ1hYCdJktQQIxrYJZmRpHQsc9rWp85zd5LHk8xK8qKOfayb5Nwk8+rl3CSTO/Jsk+Sqeh93Jfl4nMxHkiQ1zGhosbsZmNq2bNO27mjgI8ChwMuAOcBlSSa25fkWsB2wV71sB5zbWplkEnAZcHe9jw8CR9b7lSRJaozRMN3JU6WUOZ2JdYva4cCnSik/qNMOAO4F3g58LclWVMHcTqWUX9V5DgauSbJlKeWPwDuACcCBpZRFwE1JXgB8JMnJxR/RkyRJDTEaWuy2qLtaZyf5TpLN6vRpwBTg0lbGOjC7CnhFnbQzMK8V1NV5rgXmdeS5qt625RLg2cDzeipQkvFJJrUWYGJP+SRJkkaTkQ7sfgXsD7wOOJgqkLs6yfr1/1C10LW7t23dFKCn30mZ25Gnp33QlqfTMVTBYWu5c0UnIkmSNNJGtCu2lPLTtoc3JrkGuA04ALi2la1js3Sk9dSVuqI86SW95TPAyW2PJ2JwJ0mSRrmRbrFbSinlUeBGYAuqgRKwbKvahixpcZsDbNTDrjboyNPTPmDZlrxWORaVUua3FmBBn09CkiRphIyqwC7JeGAr4B5gNlVQNr1t/erAbsDVddI1wDpJXt6WZ0dgnY48u9bbtuxJNUr29iE5EUmSpBEw0vPYfT7Jbkmm1QHZ94BJwNn1aNVTgGOT7JNka2Am8BjVFCeUUm4FLgbOSLJTkp2AM4CL6hGx1HkXATOTbJ1kH+BYwBGxkiSpUUZ6upNNgG8DzwLuo7qvbqdSyt/q9Z8F1gC+DKxLNdhiz1JKe9foO4BTWTJ69kdU894BUEqZl2Q6cDpwHfAQ1f1z7ffQSZIkrfRGevDEW1ewvgAz6qW3PA8C71zBfm4Edu1/CSVJklYeo+oXmjEBAAAcrUlEQVQeO0mSJHXPwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIUZNYJfkmCQlySltaeOTfCnJ/UkeTfKjJJt0bLdpkgvr9fcnOTXJ6h15dktyfZKFSf6a5L3DdV6SJEnDZVQEdkleBhwC/L5j1SnAPsBbgV2AtYGLkoyptxsD/BhYq17/VmBf4KS2fU8DfgL8HNge+DRwapJ9h/CUJEmSht3YkS5AkrWBbwIHA//Rlr4OcBCwXynl8jrtncAdwB7AJcCewD8Azyml3F3nOQKYmeRjpZT5wHuBv5dSDq93fWuSHYAjge8PwylKkiQNi9HQYnc68ONW8NbmpcA44NJWQh283QS8ok7aGbipFdTVLgHG19u38lzK0i4BdkgyrqcC1V3Ak1oLMLH/pyVJkjS8RjSwS/JWqgDsmB5WTwGeKKU81JF+b72ulefe9pV1/ieWl6d+PBZ4Vi9FOwaY17bcuaJzkSRJGmkjFtgleQ7wReAdpZSF/dkUKG2PSxd5spxtAT4DrNO2bNJLPkmSpFFjJFvsXgpsCFyf5KkkTwG7AYfV/98LrJ5k3Y7tNmRJC9wclrTMAVDnH7e8PPU+ngIe6KlgpZRFpZT5rQVY0M0JSpIkDaeRDOyuALYBtmtbrqMaSNH6/0lgemuDJFOBrYGr66RrgK3r9JY9gUXA9W15prO0PYHrSilPDuL5SJIkjagRGxVbSllANRBisSSPAg+UUm6qH38DOCnJA8CDwOeBG4HWQItLgVuAc5McBaxX5zmjbmkD+CpwaJKTgTOoBlMcBLxtCE9PkiRp2I34dCcr8GGqLtPzgTWoWvkOLKU8DVBKeTrJG4AvA78EHge+RTWVCXWe2Un2Br4AfAC4GzislOJUJ5IkqVFGVWBXStm94/FC4IP10ts2fwfeuIL9XgW8ZBCKKEmSNGqNhnnsJEmSNAgM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIaoqvALsk7k0wY7MJIkiSpe9222J0MzEnytSQvH8wCSZIkqTvdBnbPBt4NTAV+keTmJEck2WDwiiZJkqT+6CqwK6U8VUr5QSnlTcCmwNlUgd6dSX6Q5A1JMpgFlSRJ0vINePBEKWUOcAUwCyjADsC3gD8nedVA9y9JkqS+6TqwS/KsJIcn+R3wS2BD4M3Ac4GNgYuAcwallJIkSVqhsd1slOSHwN7AbOBM4OxSyn1tWR5J8lngsIEXUZIkSX3RVWAHzAf2KKX8fDl57gG26HL/kiRJ6qeuArtSygF9yFOA27rZvyRJkvqv2wmKv5Dk0B7SP5DkpIEXS5IkSf3V7eCJfwGu7SH9GuAt3RdHkiRJ3eo2sHsW8FAP6fPrdZIkSRpm3QZ2twGv6yH9dVQjZSVJkjTMuh0VewpwSpL1gZ/Vaa8FjgaOHIyCSZIkqX+6HRV7RpIJwLHAf9bJdwKHlVL+72AVTpIkSX3XbYsdpZQvAV9KMhV4vJTy8OAVS5IkSf3VdWDXUkq5ZzAKIkmSpIHpdh67DZKcleTvSRYmeaJ9GexCSpIkacW6bbGbCTwf+BzVT4eVwSqQJEmSutNtYLcrsGsp5beDWRhJkiR1r9t57O7EVjpJkqRRpdvA7sPAZ5JsMpiFkSRJUve67Yo9F5gI/C3JfODJ9pWllA0HWjBJkiT1T7eB3b8PaikkSZI0YN3+8sQ3BrsgkiRJGphu77EjyfOSzEhybpIN67Q9k2w1eMWTJElSX3XVYpfkVcDFwK+BVwDHA3OBlwAHA/8yWAWUJEnNMmPWjJEuwoDN2H3GSBehR9222J0IzCilvBpo/6WJnwE7DbhUkiRJ6rduA7ttge/1kD4X2KD74kiSJKlb3QZ2DwNTekjfDrirrztJ8r4kv08yv16uSfL6tvXjk3wpyf1JHk3yo86585JsmuTCev39SU5NsnpHnt2SXF//ru1fk7y3n+crSZI06nUb2H0HOCHJBtS/QJFkR+DzwHn92M+dVFOn7FAvPwP+O8mL6vWnAPsAbwV2AdYGLkoypj7mGODHwFr1+rcC+wIntQ6QZBrwE+DnwPbAp4FTk+zb77OWJEkaxbqdx+5YqkmK76EKDm8BxgHnA5/o605KKRd2JH0syfuAnZLcCRwE7FdKuRwgyTuBO4A9gEuAPYF/AJ5TSrm7znMEMDPJx0op84H3An8vpRxeH+PWJDsARwLf7/eZS5IkjVLdzmP3BPCWJC+gGgm7GvCbUsofui1I3fr2L1Stb9cAL6UKFi9tO+7dSW6iGol7CbAzcFMrqKtdAoyvt7+yznMpS7sEOCjJuFLKkx3rSDK+3kfLxG7PS5Ikabh022IHQCnlT8CfBrKPJNtQBXITgEeAfUoptyTZDniilPJQxyb3suT+vin14/YyPZTkieXlqR+PBZ5F1erY6RiqKVwkSZJWGt3OY/f15a0vpRzSj939kWrQxWSq++POTrLb8g5PfV9f63Bd5MlytgX4DHBy2+OJVPcDSpIkjVrdtthN7Xg8DngRVQD0P/3ZUd2t+5f64XVJXgZ8CPgusHqSdTta7TYErq7/nwPs2L6/JOvW5bm3LU/nCN4NgaeAB3op0yJgUds++3NKkiRJI6Lbe+z+sTMtyVjgK8CtAyxTqO5vux54EphONSiDJFOBrYGj67zXUA24mFpKaXWp7kkVlF3flqezvHsC1/V0f50kSdLKakD32LUrpTyV5PPALJbuxuxVkk8DP6Ua6TqRarqS3YG9SinzknwDOCnJA8CDVNOp3AhcXu/iUqoRuecmOQpYr85zRj0iFuCrwKFJTgbOoBpMcRDwtgGdsCRJ0igzaIFdbTOqbtC+2ohq2pSpwDzg91RB3WX1+g9TdZmeD6wBXAEcWEp5GqCU8nSSNwBfBn4JPA58i2oqE+o8s5PsDXwB+ABwN3BYKcWpTiRJUqN0O3jis51JVMHZm4Bv9nU/pZSDVrB+IfDBeuktz9+BN65gP1dRTcsiSZLUWN222O3c8fgZ4D6qX5E4Y0AlkiRJUle6HTzxqsEuiCRJkgam29+KlSRJ0ijT7T12/0vvk/supZTy8m6OIUmSpP7p9h67K4F/o/o5sWvqtJ2ALYGv0Ta5ryRJkoZHt4HdZOD0Usqx7YlJPgVsVEp5z4BLJkmSpH7p9h67fwXO6iF9JvAvXZdGkiRJXes2sFsEvKKH9FdgN6wkSdKI6LYr9lTgq0m2B66t03YCDgY+PRgFkyRJUv90O4/dp5LMBj4EvLtOvhU4uJTyrcEqnCRJkvqu69+KrQM4gzhJkqRRousJipNMSnJgkv9Ksm6d9uIkUweveJIkSeqrbico3hq4HHgMeA7VaNiHqEbLbgIcMEjlkyRJUh9122L3Bapu2OcDC9vSfwzsOtBCSZIkqf+6DexeBny5lNL5s2J3AXbFSpIkjYBuA7sngLV7SN8CuL/74kiSJKlb3QZ2PwKOS9K6R68k2Rg4AfjBoJRMkiRJ/dJtYHcE8GxgDrAG8DPgr1T32x27nO0kSZI0RLqdoHheklcA04GXUAWIvwEu6eG+O0mSJA2Dfgd2ScYBPwHeX0q5FLh00EslSZKkfut3V2wp5Ulge8CWOUmSpFGk23vszgPeNZgFkSRJ0sB0/VuxwKFJ9gCuAx5tX1FKOXpApZIkSVK/dRvYvRT4ff3/th3r7KKVJEkaAf0K7JJsBswupbxqiMojSZKkLvX3Hrs/Axu0HiT5bpKNBrdIkiRJ6kZ/A7t0PN4bWGuQyiJJkqQB6HZUrCRJkkaZ/gZ2hWUHRzhYQpIkaRTo76jYADOTLKofTwC+mqRzupN/HozCSZIkqe/6G9id3fH4vMEqiCRJkgamX4FdKcVfm5AkSRqlHDwhSZLUEAZ2kiRJDWFgJ0mS1BAGdpIkSQ1hYCdJktQQBnaSJEkNYWAnSZLUEAZ2kiRJDWFgJ0mS1BAGdpIkSQ1hYCdJktQQBnaSJEkNYWAnSZLUEAZ2kiRJDWFgJ0mS1BAGdpIkSQ1hYCdJktQQBnaSJEkNYWAnSZLUEAZ2kiRJDTGigV2SY5L8b5IFSeYmuSDJlh15xif5UpL7kzya5EdJNunIs2mSC+v19yc5NcnqHXl2S3J9koVJ/prkvcNxjpIkScNlpFvsdgNOB3YCpgNjgUuTrNWW5xRgH+CtwC7A2sBFScYA1H9/DKxVr38rsC9wUmsHSaYBPwF+DmwPfBo4Ncm+Q3lykiRJw2nsSB68lLJX++Mk7wLmAi8F/ifJOsBBwH6llMvrPO8E7gD2AC4B9gT+AXhOKeXuOs8RwMwkHyulzAfeC/y9lHJ4fahbk+wAHAl8f4hPU5IkaViMdItdp3Xqvw/Wf18KjAMubWWog7ebgFfUSTsDN7WCutolwPh6+1aeS1naJcAOScZ1FqLu/p3UWoCJ3Z+SJEnS8Bg1gV2SACcDvyil3FQnTwGeKKU81JH93npdK8+97Svr/E8sL0/9eCzwrB6Kcwwwr225s7/nI0mSNNxGTWAHnAZsC7ytD3kDlLbHpYs8Wc62n6FqPWwtm/SQR5IkaVQZFYFdki8BbwJeXUppbx2bA6yeZN2OTTZkSQvcHJa0zLX2ty5VF26veep9PAU80FmeUsqiUsr81gIs6P9ZSZIkDa+Rnu4kSU4D/hl4TSlldkeW64EnqUbMtraZCmwNXF0nXQNsXae37Aksqrdv5ZnO0vYEriulPDkY5yJJkjTSRnRULNVUJ28H/glYkKTVqjavlPJ4KWVekm8AJyV5gGpQxeeBG4HL67yXArcA5yY5ClivznNG3doG8FXg0CQnA2dQDaY4iL51+0qSJK0URror9n1U97DNAu5pW97SlufDwAXA+cAvgceAfyylPA1Q/30DsLBef36d/8jWDuqWwL2B3YEbgOOAw0opTnUiSZIaY6TnsUsf8iwEPlgvveX5O/DGFeznKuAl/S2jJEnSymKkW+wkSZI0SAzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7CRJkhrCwE6SJKkhxo50AdQwM2aMdAkGrgnnIElaJdliJ0mS1BAGdpIkSQ1hYCdJktQQBnaSJEkNYWAnSZLUEAZ2kiRJDWFgJ0mS1BAGdpIkSQ1hYCdJktQQBnaSJEkNYWAnSZLUEAZ2kiRJDWFgJ0mS1BAGdpIkSQ1hYCdJktQQBnaSJEkNYWAnSZLUEAZ2kiRJDWFgJ0mS1BAGdpIkSQ1hYCdJktQQBnaSJEkNYWAnSZLUEAZ2kiRJDWFgJ0mS1BAGdpIkSQ1hYCdJktQQBnaSJEkNYWAnSZLUEAZ2kiRJDWFgJ0mS1BAGdpIkSQ1hYCdJktQQBnaSJEkNMaKBXZJdk1yY5O4kJcmbO9YnyYx6/eNJZiV5UUeedZOcm2RevZybZHJHnm2SXFXv464kH0+S4ThHSZKk4TLSLXZrAb8DDu1l/dHAR+r1LwPmAJclmdiW51vAdsBe9bIdcG5rZZJJwGXA3fU+PggcWe9XkiSpMcaO5MFLKT8FfgrQ2YBWt6gdDnyqlPKDOu0A4F7g7cDXkmxFFcztVEr5VZ3nYOCaJFuWUv4IvAOYABxYSlkE3JTkBcBHkpxcSinDcKqSJElDbkQDuxWYBkwBLm0llFIWJbkKeAXwNWBnYF4rqKvzXJtkXp3nj3Weq+qgruUS4DPA84DZnQdOMh4Y35Y0sTOPJEkjYcasGSNdBI1iI90VuzxT6r/3dqTf27ZuCjC3h23nduTpaR/tx+h0DDCvbbmzb0WWJEkaOaM5sGvp7CpNR1pPXakrypNe0ls+A6zTtmzSp5JKkiSNoNHcFTun/jsFuKctfUOWtLjNATbqYdsNOvJ0tsxtWP/tbMkDqi5fYHHXrQNoJUnSymA0t9jNpgrKprcSkqwO7AZcXSddA6yT5OVteXakamVrz7NrvW3LnlSjZG8fqsJLkiQNt5Gex27tJNsl2a5OmlY/3rQerXoKcGySfZJsDcwEHqOa4oRSyq3AxcAZSXZKshNwBnBRPSKWOu8iYGaSrZPsAxwLOCJWkiQ1ykh3xe4AXNn2+OT679nAgcBngTWALwPrAr8C9iylLGjb5h3AqSwZPfsj2ubFK6XMSzIdOB24DnioPs7JSJIkNchIz2M3iyUDGXpaX4AZ9dJbngeBd67gODcCu3ZTRkmSpJXFaL7HTpIkSf1gYCdJktQQI32PnTT6zJgx0iUYuCacgySp32yxkyRJaggDO0mSpIYwsJMkSWoIAztJkqSGMLCTJElqCAM7SZKkhjCwkyRJaggDO0mSpIYwsJMkSWoIAztJkqSGMLCTJElqCAM7SZKkhjCwkyRJaggDO0mSpIYwsJMkSWoIAztJkqSGMLCTJElqCAM7SZKkhhg70gWQJGk4zZg1Y6SLIA0ZW+wkSZIawhY7qYlmzBjpEgzMyl5+SRohtthJkiQ1hIGdJElSQxjYSZIkNYSBnSRJUkMY2EmSJDWEgZ0kSVJDGNhJkiQ1hPPYSRp9mjCPXRPOQdJKx8BOktRn/hyXNLrZFStJktQQBnaSJEkNYWAnSZLUEAZ2kiRJDWFgJ0mS1BCOipWkodCE6U6acA7SKsYWO0mSpIawxU6Shsms22eNdBH6ZZZz1kkrHVvsJEmSGsIWO0lSj3afOWukizBgsw7cfaSLIA0rAztJK4WVrRtTo0MTgtMmMMAePnbFSpIkNYQtdtIK2FIkSVpZ2GInSZLUEAZ2kiRJDWFXrIaU3ZiSpCYMYllZBoDYYidJktQQttiNYrZ2SZKk/lilWuySvD/J7CQLk1yf5FUjXSZJkqTBssoEdkneApwCfArYHvg58NMkm45owSRJkgbJKhPYAR8BvlFKObOUcmsp5XDgDuB9I1wuSZKkQbFK3GOXZHXgpcAJHasuBV7RQ/7xwPi2pIkA8+fPH6oiVhYtWurho088NbTHkyRJfbLo0aU/o4c6Juh2/ymlDHJRRp8kzwbuAl5ZSrm6Lf1Y4IBSypYd+WcAxw9rISVJkpa1SSnlrr5mXiVa7Np0RrHpIQ3gM8DJHWnrAQ8ORaFqE4E7gU2ABUN4HK2YdTE6WA+jh3UxelgXo8Nw1cNE4O7+bLCqBHb3A08DUzrSNwTu7cxcSlkELOpIHtI21yStfxeUUoa4z1fLY12MDtbD6GFdjB7WxegwjPXQ732vEoMnSilPANcD0ztWTQeuXnYLSZKklc+q0mIHVdfquUmuA64BDgE2Bb46oqWSJEkaJKtMYFdK+W6S9YGPA1OBm4C9Syl/G9mSLbYI+E+W7QLW8LMuRgfrYfSwLkYP62J0GLX1sEqMipUkSVoVrBL32EmSJK0KDOwkSZIawsBOkiSpIQzsJEmSGsLATpIkqSEM7IZRkvcnmZ1kYZLrk7xqBfn3TXJLkkX1332Gq6xN15+6SHJwkp8neaheLk/y8uEsb1P19zXRtt1bk5QkFwx1GVcVXbw/TU5yepJ76m1uTbL3cJW3ybqoi8OT/DHJ40nuSPKFJBOGq7xNlGTXJBcmubt+r3lzH7bZra6vhUn+muS9w1HWTgZ2wyTJW4BTgE8B2wM/B36aZNNe8u8MfBc4F3hx/ff8JDsOT4mbq791AewOfBt4NbAz8Hfg0iQbD31pm6uLemht91zg83V+DYIu3p9WBy4Dngf8H2BL4GCgzz9Urp51URfvAE6gmlNtK+Ag4C1Uv3mu7q0F/A44tC+Zk0wDfkJVX9sDnwZOTbLvkJWwt7I4j93wSPIr4DellPe1pd0KXFBKOaaH/N8FJpVSXt+WdjHwUCnlbcNR5qbqb130sP0Y4CHg0FLKOUNX0mbrph7qa38VcBbwKmByKWWF36S1fF28P70XOAp4YSnlyeErafN1URenAVuVUl7blnYS8PJSSp9awLV8SQqwTyml1x6CJCcCbyqlbNWW9lXgxaWUnYehmIvZYjcM6m+3LwUu7Vh1KfCKXjbbuYf8lywnv/qgy7rotCYwDnhwEIu2ShlAPXwcuK+U8o2hKtuqpsu6eBPVTzOenuTeJDclObYOvNWlLuviF8BLW7eHJNkM2Bv48VCVUz3q7TN7hyTjhrMgq8xPio2wZwFjgHs70u8FpvSyzZR+5lffdFMXnU6g6nK6fBDLtarpdz0keSVVN9N2Q1u0VU43r4nNgNcA36QKIrYATqf6TPmvoSnmKqHfdVFK+U6SDYBfJAlVHXyllHLCkJZUnXr7zB5LVa/3DFdBDOyGV2e/d3pIG0h+9V1X1zbJ0cDbgN1LKQuHomCrmD7VQ5KJwHnAwaWU+4ejYKug/rwmVgPmAoeUUp4Grk/ybKruWQO7getzXSTZHfgY8H7gV8DmwBeT3FNK+cRQFlLL6KneekofUgZ2w+N+4GmW/ca1IctG+C1z+plffdNNXQCQ5EjgWGCPUsrvh6Z4q4z+1sPzqW7Uv7BqlADqW0mSPAVsWUq5bUhK2nzdvCbuAZ6sg7qWW4EpSVYvpTwx+MVcJXRTF58Azi2lnFk/vjHJWsDXk3yqlPLM0BRVHXr7zH4KeGA4C+I9dsOgfpO7HpjesWo6cHUvm13TQ/49l5NffdBlXZDkKOA4YK9SynVDV8JVQxf18AdgG6pu2NbyI+DK+v87hqywDdfla+KXwOZJ2j9DXgDcY1DXvS7rYk2gM3h7mqq1KMtm1xDp7TP7umEfYFRKcRmGhWr4+RPAu6mGpH8BeAR4br3+HOAzbflfQRXpfxR4Yf33SWDHkT6XlX3poi6OBhYB+1J9I2sta4/0uazMS3/roYftZ1KNFBzxc1nZly5eE88BFgBfogro3kDVovSxkT6XlX3poi5mAPOBtwLTqIKLvwDfHelzWZkXYG2WfIkswIfr/zet138GOKct/zTgUeDkut7eXdfjvsNddrtih0kp5btJ1qca1TcVuAnYu5TytzrLprR96yqlXJ3krcAnqZrabwPeUkr51fCWvHn6WxdU966sDnyvY1f/SfWmqi50UQ8aIl28P92RZE+qoOP3VIOJvgicOKwFb6AuXhefpAo8PglsDNwHXEh13526twNVj0DLyfXfs4EDqepm8dyCpZTZ9QTdXwA+ANwNHFZK+f6wlLaN89hJkiQ1hPfYSZIkNYSBnSRJUkMY2EmSJDWEgZ0kSVJDGNhJkiQ1hIGdJElSQxjYSZIkNYSB3f9vtw5kAAAAAAb5W9/jK4oAACbEDgBgQuwAACbEDgBgInCNXU2YSWNWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot calibration distribution\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams.update({'figure.figsize':(7,5), 'figure.dpi':100})\n",
    "\n",
    "CalibrationDf['probabilities'] = calibration_probabilities\n",
    "\n",
    "x1 = CalibrationDf.loc[CalibrationDf.label==1, 'probabilities']\n",
    "x2 = CalibrationDf.loc[CalibrationDf.label==0, 'probabilities']\n",
    "\n",
    "kwargs = dict(alpha=0.5, bins=10, stacked=True)\n",
    "\n",
    "plt.hist(x1, **kwargs, color='g', label='Include')\n",
    "plt.hist(x2, **kwargs, color='r', label='Exclude')\n",
    "plt.gca().set(title='Frequency histogram of calibration probabilities', ylabel='Frequency')\n",
    "\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEICAYAAABI7RO5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xu8VXWd//HXW0BBRUnABjnooQkRL2mIZBeNES11SrTRxLFfoA78SrNxrEbTJo+lPWrGMnvYZHj54S1AbUx0nMxEtMkreEkNKUqSE15OggiJyuXz+2N9D26O+5yz1jln730u7+fjcR5nXb5rrc937ctnf7/ftddWRGBmZpbXNrUOwMzMehYnDjMzK8SJw8zMCnHiMDOzQpw4zMysECcOMzMrxInDOkVSSHpvK+tOlvSLasfUnUgaJOl2SWsk3VzF4zZIuqFC+z5E0tJK7LsrSZokqbGD29an53b/VtafJ+mqcmUl/Y+kaW3s+wpJ/9aRuLqLsifF2iZpOfBuYFPJ4j0jYmVtIuqeIuJG4Mb2ykmaDTRGxNcqHlT1HU/2XBkaERtrHUxHSApgTEQsA4iIXwFjaxtVbUXEt9pYd1TztKTpwD9FxEdK1n+ustFVnlscHffJiNix5O8dSaO1TytWPd3gMdgD+F1PTRrdWTd4bPuuiPBfwT9gOXB4meX1QACnAc8D96flBwMPAK8CTwKTSrYZDdwHrAXuBi4HbkjrJpF9Ei97bLLEfy7wB+AV4CZglxaxTEux/AU4v2Q//YDz0rZrgcXAKOCHwHdbHPN24KxWzkUAnwN+D6xO2yutmw78b5oWcCnwMrAG+A2wLzAT2AC8BawDbk/lxwEL0zl7Bjim5JhDU0yvAY8CFzUfpySmM1JMz6VllwEr0jaLgUNKyjcANwM3pHPxFLAn8NUU7wrgY208H8rGClyY6rUh1e20Mtu29Rj+HPhCi/JPAp/KWae8z6OJwIMp/hfInoPbpnX3p/P511SHE1vur53HanZ6Tvx3OrcPA3/bynmsT8eaCaxMsXypRZ1uSY/Ta8A/AdsB30/lV6bp7UrrTfY8/0uq88kl+/t74PG0rxVAQ8FYbmhRtn+aX5hiGwe8QdYzsQ54teScXFSyr08AT6Tz9wDwvpJ15wB/TuduKTC51u9/EeHE0aGT1n7iuA7YARgEjCR7Qzia7E3iiDQ/PG3zIPC99AI4ND1B8r7gzwIeAurS9j8G5rSI5coUx/7Am8C4tP4rZG+QY8ne1Pcne0OemF4o26Ryw4DXgXe3ci4CuAMYAuwONAFHpnXTeTtxfJzszW1IOt44YERa1/KFNABYRvaC3xY4LJ2XsWn93PS3PbA32Yu+ZeK4G9gFGJSWfSbVrz/wJeBFYGBa10D2Av94Wn8d8BxwfoplBikBlal/e7E2ND+erWzf1mP4WeDXJWX3Jntz2S5nnfI+jw4k+3DTn+x5s4SSDwrpfL63ZH7L/nLUfzawiux51Z+s63JuK+eiPh1rDtnrZz+y59PhJXXaABxL9loaBHwjnb9dgeFkb7zfLIlzI2+/vj5KlgDHlqzfL+3rfcBLwLEFYmkzcbR8DZTUczbp+Q6MJ/tw8gGyD3PT0mOzHdlrcwWwW8lxyibdqr8H1jqAnviXHth1ZC/iV4GftXgCvaek7DnA9S22vys9QXZPT+wdStb9hPwv+CWUfAIBRqQXVv+SWOpK1j8CTE3TS4EprdRvCXBEmv4CcGcb5yKAj5TM3wScm6a3vGjI3lB+R/YGtU2LfWx5IaX5Q8jeBLcpWTYnvVj7pTqOLVlXrsVxWDuP4Wpg/zTdANxdsu6T6fHtl+YHp30OKbOfVmMt2XdbiaOtx3Aw2RvdHmndxcA1BeqU63lUZj9nAbe2OJ+tJY726j8buKpk3dHAs60ctz4da6+SZf8OXF1Sp/tbbPMH4OiS+Y8Dy0vibPn6ugn4t1aO/33g0gKxdEXi+BEp0ZWsX0qW5N5LllQOBwa09Xyu9p/HODru2IgYkv6ObbFuRcn0HsAJkl5t/gM+QvYGsRuwOiL+WlL+TwVi2AO4tWS/S8iaxe8uKfNiyfTrwI5pehTZi66ca8k+zZL+X99OHK0dY4uIWEDWBfJD4CVJsyTt1Mr+dgNWRMTmkmV/Imu9DSd7Uy09x6XTZZdJ+pKkJenqpleBnclaU81eKpleD/wlIjaVzFOuXu3Emkerj2FErCXr4pmayk6l5GKDHHXKRdKeku6Q9KKk14BvFdhPnvq3+/xoofSx+1M6Rrl1zccvfc20LF/u9bUbgKQPSLpXUpOkNWRdri3r3VYsXWEP4Est3h9GkbUylpEl8QbgZUlzJXX18TvEiaMyomR6BVmLY0jJ3w4R8W2yftN3SdqhpPzuJdN/JeuOAUBSP7I3ztJ9H9Vi3wMj4s85YlwB/G0r624Apkjan6xL6Wc59teuiPhBRBwI7EM2hvCV5lUtiq4ERkkqfX7uTtbX20T2KbKuZN2ocodrnpB0CFnL79PAuyJiCNk4izpem1yx5tHeYzgHOEnSB8m6Zu6FwnVq73n0I+BZsiundiLrdsp7bjpb/3JKH8/d0zGalXuu7NFG+XKvr+b1PwHmA6MiYmfgCt5Z77ZiyaNlvC2tAC5u8fhvHxFzACLiJ5FdkbVH2td3Ch6/Ipw4Ku8G4JOSPi6pn6SB6fryuoj4E7AIuFDStpI+QtZN0ux3wEBJfy9pAPA1sr7PZlcAF0vaA0DScElTcsZ1FfBNSWOUeZ+koQAR0Ug26Hw98NOIWN/WjvKQdFD6hDeA7I2sedAQsk/77ykp/nAq86+SBkiaRHZe5qZWwH8BDZK2l7QX2VhAWwaTJZsmoL+krwOttXaKajXWnNu39xjeSfam8Q1gXskn+yJ1au95NJhsgHhdOp+fb7F9y8enVGfrX86/pcd2H+AUYF4bZecAX0vnbRjwdbLXXKnm19chZAPRzd+nGQysiog3JE0E/rGTsZTzElAnadtW1l8JfC69NiRph/Q4DZY0VtJhkrYje72sZ+uvANSME0eFRcQKYArZp7gmsk8YX+Htc/+PZANjq4ALyAZmm7ddA5xO9ib/Z7IXaOkXmi4j+8T0C0lryQYJP5AztO+R9ff+guxN42qyT7TNriUbEGyvmyqvncheJKvJmvyvAJekdVcDe6em+s8i4i3gGOAosqth/hP4bEQ8m8p/gaxb5sUU3xyygf/W3AX8D9kb6J/IXoTlurcKyxFre9p8DCPiTbJEeTjZJ+RmueuU43n0ZbLn4Vqyx6jlm2MDcG16fD7dYt+drX8595ENuN8DXBIRbX2J9CKyD1+/IbvY47G0rNmLZM+5lWTdfJ8rie104BvpvH+d7PXQmVjKWUB2pdmLkv7ScmVELCK7+OLyFOcysnERyJL7t8nO64tkFwCcV/D4FdF82aR1E5IayAYiP9Ne2QrHcSjZJ7f6Fv3X3Y6k7wB/ExHTah2LdZykerKr2QaEv/fSrbnFYe+QujP+mexqmG6XNCTtlbrWlLoYTgNurXVcZn2FE4dtRdI4skuMR5BdntgdDSbrvvkrWffCd4HbahqRWR/iriozMyvELQ4zMyukV94kbNiwYVFfX1/rMMzMepTFixf/JSKGt1euVyaO+vp6Fi1aVOswzMx6FEm57lzhriozMyvEicPMzApx4jAzs0J65RiHGcCGDRtobGzkjTfeqHUo3d7AgQOpq6tjwIABtQ7FegAnDuu1GhsbGTx4MPX19UhdcSPc3ikieOWVV2hsbGT06NG1Dsd6AHdVWa/1xhtvMHToUCeNdkhi6NChbplZbk4c1qs5aeTj82RFOHGYmVkhHuOwPqNhYUPX7m9S+/vbcccdWbduXeF9L1y4kEsuuYQ77rgjfzwNDey44458+ctfLnw8syKcOMppaOhbxzUzK8BdVWZVsHDhQiZNmsTxxx/PXnvtxcknn0zznakfffRRPvShD7H//vszceJE1q5du9W2DQ0NXHLJJVvm9913X5YvXw7AxRdfzNixYzn88MNZunTpljJ/+MMfOPLIIznwwAM55JBDePbZzvwgn9nW3OIwq5LHH3+cZ555ht12240Pf/jD/PrXv2bixImceOKJzJs3j4MOOojXXnuNQYMGtb8zYPHixcydO5fHH3+cjRs3Mn78eA488EAAZs6cyRVXXMGYMWN4+OGHOf3001mwYEElq2d9iBOHWZVMnDiRuro6AA444ACWL1/OzjvvzIgRIzjooIMA2GmnnXLv71e/+hXHHXcc22+/PQDHHHMMAOvWreOBBx7ghBNO2FL2zTfb+kl2s2KcOMyqZLvtttsy3a9fPzZu3EhEtHspbP/+/dm8+e1f8C39vkW5bTdv3syQIUN44oknuiBqs3fyGIdZDe21116sXLmSRx99FIC1a9eycePGrcrU19fz2GOPAfDYY4/x3HPPAXDooYdy6623sn79etauXcvtt98OZK2W0aNHc/PNNwPZN8OffPLJalXJ+gC3OKzPyHP5bLVtu+22zJs3jzPPPJP169czaNAgfvnLX25V5h/+4R+47rrrOOCAAzjooIPYc889ARg/fjwnnngiBxxwAHvssQeHHHLIlm1uvPFGPv/5z3PRRRexYcMGpk6dyv7771/Vulnv1St/c3zChAnRqR9y8uW4vcKSJUsYN25crcPoMXy+TNLiiJjQXjl3VZmZWSFOHGZmVogTh5mZFVKxxCHpGkkvS3q6ZNl/SHpW0m8k3SppSMm6r0paJmmppI+XLD8yLVsm6dxKxWtmZvlUssUxGziyxbK7gX0j4n3A74CvAkjaG5gK7JO2+U9J/ST1A34IHAXsDZyUypqZWY1ULHFExP3AqhbLfhERzRepPwTUpekpwNyIeDMingOWARPT37KI+GNEvAXMTWXNzKxGavk9jlOBeWl6JFkiadaYlgGsaLH8A+V2JmkmMBNg991379JArZfo6sudc+yvX79+7Lffflvmp06dyrnnFutxnT17NosWLeLyyy/Pvc306dP5xCc+wfHHH1/oWGZ51CRxSDof2Ajc2LyoTLGgfIuo7BdPImIWMAuy73F0QZhmnTZo0CDf+sN6naonDknTgE8Ak+Ptbx82AqNKitUBK9N0a8vNeqQ1a9YwceJE5s+fz9ixYznppJM47LDDmDFjBj//+c8577zz2LRpE8OGDeOee+7ZatuWLYnmH4qKCM4880wWLFjA6NGjKf1i7+LFizn77LNZt24dw4YNY/bs2YwYMaKqde6LuvqHw3Iftwp3SKjq5biSjgTOAY6JiNdLVs0HpkraTtJoYAzwCPAoMEbSaEnbkg2gz69mzGadsX79eg444IAtf/PmzWPnnXfm8ssvZ/r06cydO5fVq1czY8YMmpqamDFjBj/96U958sknt9xrKo9bb72VpUuX8tRTT3HllVfywAMPALBhwwbOPPNMbrnlFhYvXsypp57K+eefX6nqWh9RsRaHpDnAJGCYpEbgArKrqLYD7k539XwoIj4XEc9Iugn4LVkX1hkRsSnt5wvAXUA/4JqIeKZSMZt1tda6qo444ghuvvlmzjjjjC03IHzooYc49NBDGT16NAC77LJL7uPcf//9nHTSSfTr14/ddtuNww47DIClS5fy9NNPc8QRRwCwadMmtzas0yqWOCLipDKLr26j/MXAxWWW3wnc2YWhmdXc5s2bWbJkCYMGDWLVqlXU1dUVvsV6RPDWW29tWVdu24hgn3324cEHH+zaClif5m+Om9XApZdeyrhx45gzZw6nnnoqGzZs4IMf/CD33Xffltumr1q16h3b1dfXs3jxYgBuu+02NmzYAGS3WJ87dy6bNm3ihRde4N577wVg7NixNDU1bUkcGzZs4Jln3Gi3zvFt1a3vqMHdh5vHOJodeeSRnHrqqVx11VU88sgjDB48mEMPPZSLLrqICy+8kFmzZvGpT32KzZs3s+uuu3L33Xdvtb8ZM2YwZcoUJk6cyOTJk9lhhx0AOO6441iwYAH77bcfe+65Jx/96EeB7Lbtt9xyC1/84hdZs2YNGzdu5KyzzmKfffap3kmwXse3VS/Ht1XvFXyb8GJ8vrpWT7yqyrdVNzOzinDiMDOzQpw4rFfrjV2xleDzZEU4cVivNXDgQF555RW/KbYjInjllVcYOHBgrUOxHsJXVVmvVVdXR2NjI01NTbUOpdsbOHAgdXV17Rc0w4nDerEBAwZs+Ra2mXUdd1WZmVkhThxmZlaIE4eZmRXixGFmZoU4cZiZWSFOHGZmVogTh5mZFeLEYWZmhThxmJlZIU4cZmZWiBOHmZkV4sRhZmaFOHGYmVkhThxmZlZIxRKHpGskvSzp6ZJlu0i6W9Lv0/93peWS9ANJyyT9RtL4km2mpfK/lzStUvGamVk+lWxxzAaObLHsXOCeiBgD3JPmAY4CxqS/mcCPIEs0wAXAB4CJwAXNycbMzGqjYokjIu4HVrVYPAW4Nk1fCxxbsvy6yDwEDJE0Avg4cHdErIqI1cDdvDMZmZlZFVV7jOPdEfECQPq/a1o+ElhRUq4xLWtt+TtImilpkaRF/qlQM7PK6S6D4yqzLNpY/s6FEbMiYkJETBg+fHiXBmdmZm+rduJ4KXVBkf6/nJY3AqNKytUBK9tYbmZmNVLtxDEfaL4yahpwW8nyz6arqw4G1qSurLuAj0l6VxoU/1haZmZmNdK/UjuWNAeYBAyT1Eh2ddS3gZsknQY8D5yQit8JHA0sA14HTgGIiFWSvgk8msp9IyJaDribmVkVVSxxRMRJrayaXKZsAGe0sp9rgGu6MDQzM+uE7jI4bmZmPYQTh5mZFeLEYWZmhThxmJlZIU4cZmZWiBOHmZkV4sRhZmaFOHGYmVkhThxmZlaIE4eZmRXixGFmZoU4cZiZWSFOHGZmVogTh5mZFeLEYWZmhThxmJlZIU4cZmZWiBOHmZkV4sRhZmaFOHGYmVkhThxmZlaIE4eZmRWSK3FI2rfSgZiZWc+Qt8VxhaRHJJ0uaUhnDyrpXyQ9I+lpSXMkDZQ0WtLDkn4vaZ6kbVPZ7dL8srS+vrPHNzOzjsuVOCLiI8DJwChgkaSfSDqiIweUNBL4IjAhIvYF+gFTge8Al0bEGGA1cFra5DRgdUS8F7g0lTMzsxrJPcYREb8HvgacA3wU+IGkZyV9qgPH7Q8MktQf2B54ATgMuCWtvxY4Nk1PSfOk9ZMlqQPHNDOzLpB3jON9ki4FlpC9wX8yIsal6UuLHDAi/gxcAjxPljDWAIuBVyNiYyrWCIxM0yOBFWnbjan80DIxzpS0SNKipqamIiGZmVkBeVsclwOPAftHxBkR8RhARKwka4XkJuldZK2I0cBuwA7AUWWKRvMmbax7e0HErIiYEBEThg8fXiQkMzMroH/OckcD6yNiE4CkbYCBEfF6RFxf8JiHA89FRFPa138BHwKGSOqfWhV1wMpUvpFsbKUxdW3tDKwqeEwzM+sieVscvwQGlcxvn5Z1xPPAwZK2T2MVk4HfAvcCx6cy04Db0vT8NE9avyAi3tHiMDOz6sibOAZGxLrmmTS9fUcOGBEPkw1yPwY8lWKYRTbofrakZWRjGFenTa4GhqblZwPnduS4ZmbWNfJ2Vf1V0vjmsQ1JBwLrO3rQiLgAuKDF4j8CE8uUfQM4oaPHMjOzrpU3cZwF3CypedxhBHBiZUIyM7PuLFfiiIhHJe0FjCW7yunZiNhQ0cjMzKxbytviADgIqE/bvF8SEXFdRaIyM7NuK1fikHQ98LfAE8CmtDgAJw4zsz4mb4tjArC3L4M1M7O8l+M+DfxNJQMxM7OeIW+LYxjwW0mPAG82L4yIYyoSlZmZdVt5E0dDJYMwM7OeI+/luPdJ2gMYExG/lLQ92e9omJlZH5P3tuozyG4T8uO0aCTws0oFZWZm3VfewfEzgA8Dr8GWH3XatVJBmZlZ95U3cbwZEW81z6Tbm/vSXDOzPihv4rhP0nlkP/d6BHAzcHvlwjIzs+4qb+I4F2giuw36/wXupOAv/5mZWe+Q96qqzcCV6c/MzPqwvPeqeo7yv/P9ni6PyMzMurUi96pqNpDsh5V26fpwzMysu8s1xhERr5T8/Tkivg8cVuHYzMysG8rbVTW+ZHYbshbI4IpEZGZm3VrerqrvlkxvBJYDn+7yaMzMrNvLe1XV31U6EDMz6xnydlWd3db6iPhe14RjZmbdXZGrqg4C5qf5TwL3AysqEZSZmXVfRX7IaXxErAWQ1ADcHBH/1JGDShoCXAXsS/b9kFOBpcA8oJ40hhIRqyUJuAw4GngdmB4Rj3XkuGZm1nl5bzmyO/BWyfxbZG/wHXUZ8POI2AvYH1hCdluTeyJiDHBPmgc4ChiT/mYCP+rEcc3MrJPytjiuBx6RdCtZC+E44LqOHFDSTsChwHSAdNfdtyRNASalYtcCC4FzgCnAdRERwEOShkgaEREvdOT4ZmbWOXm/AHgxcAqwGngVOCUivtXBY76H7IaJ/0/S45KukrQD8O7mZJD+N//ex0i2HktpTMu2ImmmpEWSFjU1NXUwNDMza0/eriqA7YHXIuIyoFHS6A4esz8wHvhRRLwf+Ctvd0uVozLLyt03a1ZETIiICcOHD+9gaGZm1p68Px17AVm30VfTogHADR08ZiPQGBEPp/lbyBLJS5JGpOONAF4uKT+qZPs6YGUHj21mZp2Ut8VxHHAMWeuAiFhJB285EhEvAiskjU2LJgO/JbvUd1paNg24LU3PBz6rzMHAGo9vmJnVTt7B8bciIiQFQBqT6IwzgRslbQv8kWz8ZBvgJkmnAc+T3YEXsh+NOhpYRnY57imdPLaZmXVC3sRxk6QfA0MkzSD73kWHf9QpIp5g61u1N5tcpmwAZ3T0WGZm1rXy3qvqkvRb468BY4GvR8TdFY3MzMy6pXYTh6R+wF0RcTjgZGFm1se1OzgeEZuA1yXtXIV4zMysm8s7xvEG8JSku0lXVgFExBcrEpWZmXVbeRPHf6c/MzPr49pMHJJ2j4jnI+LaagVkZmbdW3tjHD9rnpD00wrHYmZmPUB7iaP0PlHvqWQgZmbWM7SXOKKVaTMz66PaGxzfX9JrZC2PQWmaNB8RsVNFozMzs26nzcQREf2qFYiZmfUMRX6Pw8zMzInDzMyKceIwM7NCnDjMzKwQJw4zMyvEicPMzApx4jAzs0KcOMzMrBAnDjMzK8SJw8zMCnHiMDOzQpw4zMyskJolDkn9JD0u6Y40P1rSw5J+L2mepG3T8u3S/LK0vr5WMZuZWW1bHP8MLCmZ/w5waUSMAVYDp6XlpwGrI+K9wKWpnJmZ1UhNEoekOuDvgavSvIDDgFtSkWuBY9P0lDRPWj85lTczsxqoVYvj+8C/ApvT/FDg1YjYmOYbgZFpeiSwAiCtX5PKb0XSTEmLJC1qamqqZOxmZn1ae78A2OUkfQJ4OSIWS5rUvLhM0cix7u0FEbOAWQATJkzwz9yaGQANCxtqHUKvU/XEAXwYOEbS0cBAYCeyFsgQSf1Tq6IOWJnKNwKjgEZJ/YGdgVXVD9vMzKAGiSMivgp8FSC1OL4cESdLuhk4HpgLTANuS5vMT/MPpvULIqJ3tigaGvrWcc2sR+pO3+M4Bzhb0jKyMYyr0/KrgaFp+dnAuTWKz8zMqE1X1RYRsRBYmKb/CEwsU+YN4ISqBmZmZq3qTi0OMzPrAZw4zMysECcOMzMrxInDzMwKceIwM7NCnDjMzKwQJw4zMyvEicPMzApx4jAzs0KcOMzMrBAnDjMzK8SJw8zMCnHiMDOzQpw4zMysECcOMzMrxInDzMwKceIwM7NCnDjMzKwQJw4zMyvEicPMzApx4jAzs0KcOMzMrBAnDjMzK6R/tQ8oaRRwHfA3wGZgVkRcJmkXYB5QDywHPh0RqyUJuAw4GngdmB4Rj1U7bjPrnIaFDbUOwbpILVocG4EvRcQ44GDgDEl7A+cC90TEGOCeNA9wFDAm/c0EflT9kM3MrFnVWxwR8QLwQppeK2kJMBKYAkxKxa4FFgLnpOXXRUQAD0kaImlE2o91hYaGvnVcM+uUmo5xSKoH3g88DLy7ORmk/7umYiOBFSWbNaZlZmZWAzVLHJJ2BH4KnBURr7VVtMyyKLO/mZIWSVrU1NTUVWGamVkLNUkckgaQJY0bI+K/0uKXJI1I60cAL6fljcCoks3rgJUt9xkRsyJiQkRMGD58eOWCNzPr46qeONJVUlcDSyLieyWr5gPT0vQ04LaS5Z9V5mBgjcc3zMxqp+qD48CHgf8DPCXpibTsPODbwE2STgOeB05I6+4kuxR3GdnluKdUN1wzMytVi6uq/pfy4xYAk8uUD+CMigZlZma5+ZvjZmZWiBOHmZkV4sRhZmaF1GJw3Czjb6yb9UhucZiZWSFOHGZmVoi7qsz6GN/e3DrLLQ4zMyvEicPMzApx4jAzs0I8xmFWIx5rsJ7KLQ4zMyvEicPMzApxV5X1Pf7GulmnuMVhZmaFOHGYmVkh7qqybmPh8oW1DqGypk/aanZS2ULWWQtbnGfrek4c9g69/g3czDrFXVVmZlaIWxzdmD/5mxU3afbCWodQVbXomnOLw8zMCnHiMDOzQpw4zMysEI9x5OTxBjOzTI9pcUg6UtJSScsknVvreMzM+qoekTgk9QN+CBwF7A2cJGnv2kZlZtY39YjEAUwElkXEHyPiLWAuMKXGMZmZ9Uk9ZYxjJLCiZL4R+EBpAUkzgZlpdp2kpZ043jDgL53YvidynfuGvlbn3l/fa+/bavZCLuxMnffIU6inJA6VWRZbzUTMAmZ1ycGkRRExoSv21VO4zn1DX6tzX6svVKfOPaWrqhEYVTJfB6ysUSxmZn1aT0kcjwJjJI2WtC0wFZhf45jMzPqkHtFVFREbJX0BuAvoB1wTEc9U8JBd0uXVw7jOfUNfq3Nfqy9Uoc6KiPZLmZmZJT2lq8rMzLoJJw4zMyukzyaO9m5hImk7SfPS+ocl1Vc/yq6Vo85nS/qtpN9IukdSrmu6u7O8t6qRdLykkNTjL93MU2dJn06P9TOSflLtGLtajuf27pLulfR4en4fXYs4u5KkayS9LOnpVtZL0g/SOfmNpPFddvCI6HN/ZAPsfwDeA2wLPAns3aLM6cAVaXq69iuIAAACyklEQVQqMK/WcVehzn8HbJ+mP98X6pzKDQbuBx4CJtQ67io8zmOAx4F3pfldax13Feo8C/h8mt4bWF7ruLug3ocC44GnW1l/NPA/ZN+DOxh4uKuO3VdbHHluYTIFuDZN3wJMllTui4g9Rbt1joh7I+L1NPsQ2fdlerK8t6r5JvDvwBvVDK5C8tR5BvDDiFgNEBEvVznGrpanzgHslKZ3phd8Dywi7gdWtVFkCnBdZB4Chkga0RXH7quJo9wtTEa2ViYiNgJrgKFVia4y8tS51Glkn1Z6snbrLOn9wKiIuKOagVVQnsd5T2BPSb+W9JCkI6sWXWXkqXMD8BlJjcCdwJnVCa2mir7mc+sR3+OogHZvYZKzTE+Suz6SPgNMAD5a0Ygqr806S9oGuBSYXq2AqiDP49yfrLtqElmr8leS9o2IVyscW6XkqfNJwOyI+K6kDwLXpzpvrnx4NVOx97C+2uLIcwuTLWUk9Sdr3rbVLOzuct22RdLhwPnAMRHxZpViq5T26jwY2BdYKGk5WT/w/B4+QJ73uX1bRGyIiOeApWSJpKfKU+fTgJsAIuJBYCDZDRB7s4rdqqmvJo48tzCZD0xL08cDCyKNOPVQ7dY5ddv8mCxp9PR+b2inzhGxJiKGRUR9RNSTjescExGLahNul8jz3P4Z2YUQSBpG1nX1x6pG2bXy1Pl5YDKApHFkiaOpqlFW33zgs+nqqoOBNRHxQlfsuE92VUUrtzCR9A1gUUTMB64ma84uI2tpTK1dxJ2Xs87/AewI3JyuA3g+Io6pWdCdlLPOvUrOOt8FfEzSb4FNwFci4pXaRd05Oev8JeBKSf9C1l0zvYd/EETSHLLuxmFp7OYCYABARFxBNpZzNLAMeB04pcuO3cPPnZmZVVlf7aoyM7MOcuIwM7NCnDjMzKwQJw4zMyvEicPMzApx4jAzs0KcOMzMrJD/Dwy0eGxIaQfLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot evaluation distribution\n",
    "\n",
    "EvaluationDf['probabilities'] = evaluation_probabilities\n",
    "\n",
    "x1 = EvaluationDf.loc[EvaluationDf.label==1, 'probabilities']\n",
    "x2 = EvaluationDf.loc[EvaluationDf.label==0, 'probabilities']\n",
    "\n",
    "kwargs = dict(alpha=0.5, bins=10, stacked=True)\n",
    "\n",
    "plt.hist(x1, **kwargs, color='g', label='Include')\n",
    "plt.hist(x2, **kwargs, color='r', label='Exclude')\n",
    "plt.gca().set(title='Frequency histogram of evaluation probabilities', ylabel='Frequency')\n",
    "\n",
    "plt.legend();"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}