HiggsChallenge-NN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Higgs Challenge Example using Neural Networks\n",
    "In this part we will look at the **[Higgs Boson ML Challenge](https://www.kaggle.com/c/Higgs-boson)** on Kaggle and attempt a solution using neural networks (NN). The data is available from **[CERN Open Data](http://opendata.cern.ch/record/328)**. More information about the data is available from the links, and in particular at **[Documentation](http://opendata.cern.ch/record/329/files/atlas-higgs-challenge-2014.pdf)**. The general idea is that we want to extract $H\\to\\tau\\tau$ signal from background. In particular, the selection requires one of the taus to decay into an electron or muon and two neutrinos, and the other into hadrons and a neutrino. The challenge is based on Monte Carlo events processed through the **[ATLAS detector](http://atlas.cern/)** simulation and reconstruction.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background on Neural Networks\n",
    "(based on lectures from **[ML Course on Coursera](https://www.coursera.org/learn/machine-learning)**)\n",
    "\n",
    "As we saw from the logistic regression yesterday, linear classifiers are often not the best at solving complicated problems. Neural networks introduce nonlinearity. They were originally designed to mimic the brain, and were popular in the 80s and early 90s. Recently they have become popular again, especially as deep neural networks DNNs, including convolutional NNs (CNN), recurrent NNs (RNN), etc. Those are beyond the scope of this class, but we will introduce the basics of NNs.\n",
    "\n",
    "Below is a diagram of a simple NN:\n",
    "![NNFig](https://upload.wikimedia.org/wikipedia/commons/4/46/Colored_neural_network.svg)\n",
    "It is made up of \"neurons\" that get a number of inputs, processes them, and sends the output to other neurons. Mathematically, one can represent the a neuron's \"activation\" $a = g\\left(\\theta^Tx\\right)$, where $x$ are the inputs (a vector), and $\\theta$ are the parameters (weights) of the model (also a vector), and $g$ is the activation fuction. For example, if we use a logistic function as the activation function, we can have $g\\left(\\theta^Tx\\right) = \\frac{1}{1+\\mathrm{exp}\\left(-\\theta^Tx\\right)}$, or if a Rectified Linear Unit (ReLU), $g\\left(\\theta^Tx\\right) = \\mathrm{max}\\left(0, \\theta^Tx\\right)$. The NN above has an input layer (layer 1), a hidden layer (layer 2), and an output layer (layer 3). One can have more hidden layers. Let's label the activations of layer 2 as $a_i^{(2)} = g\\left(\\theta_i^{(1)T}x\\right)$, where $i$ is the index of the individual neurons. Note that the superscript of the $\\theta$ is (1). That is because these are the weights going from layer 1 to 2. Putting together all the individual weight vectors together forms a matrix $\\Theta^{(1)}$.\n",
    "\n",
    "Using matrix notation, we can define $z^{(j)} = \\Theta^{(j-1)}a^{(j-1)}$ and then $a^{(j)} = g(z^{(j)})$. Thus evaluating the NN is a series of matrix multiplications followed by activation functions.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cost function of a NN is similar to what we have for logistic regression, modified to take into account possible multiple outputs, and with more complicated regularization. In order to train the NN, we have to determine the weight matrix $\\Theta$ that minimizes the cost function. Backpropagation is the method used to do that. It calculates the partial derivatives \"errors\" for each $z_i^{(j)}$ by propagating the errors backwards. Usually something like (stochastic) gradient descent is used to solve the problem. For more details on backprorpagation, look, for example, at the **[ML course](https://www.coursera.org/learn/machine-learning)** mentioned above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's start trying to apply a NN to the Higgs Challenge data. We will start using Scikit Learn, and then try **[Keras](https://keras.io/)**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the usual setup: \n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load training data\n",
    "df = pd.read_csv('data/atlas-higgs-challenge-2014-v2.csv.gz')"
   ]
  },
  {
   "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>EventId</th>\n",
       "      <th>DER_mass_MMC</th>\n",
       "      <th>DER_mass_transverse_met_lep</th>\n",
       "      <th>DER_mass_vis</th>\n",
       "      <th>DER_pt_h</th>\n",
       "      <th>DER_deltaeta_jet_jet</th>\n",
       "      <th>DER_mass_jet_jet</th>\n",
       "      <th>DER_prodeta_jet_jet</th>\n",
       "      <th>DER_deltar_tau_lep</th>\n",
       "      <th>DER_pt_tot</th>\n",
       "      <th>...</th>\n",
       "      <th>PRI_jet_leading_eta</th>\n",
       "      <th>PRI_jet_leading_phi</th>\n",
       "      <th>PRI_jet_subleading_pt</th>\n",
       "      <th>PRI_jet_subleading_eta</th>\n",
       "      <th>PRI_jet_subleading_phi</th>\n",
       "      <th>PRI_jet_all_pt</th>\n",
       "      <th>Weight</th>\n",
       "      <th>Label</th>\n",
       "      <th>KaggleSet</th>\n",
       "      <th>KaggleWeight</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>100000</td>\n",
       "      <td>138.470</td>\n",
       "      <td>51.655</td>\n",
       "      <td>97.827</td>\n",
       "      <td>27.980</td>\n",
       "      <td>0.91</td>\n",
       "      <td>124.711</td>\n",
       "      <td>2.666</td>\n",
       "      <td>3.064</td>\n",
       "      <td>41.928</td>\n",
       "      <td>...</td>\n",
       "      <td>2.150</td>\n",
       "      <td>0.444</td>\n",
       "      <td>46.062</td>\n",
       "      <td>1.24</td>\n",
       "      <td>-2.475</td>\n",
       "      <td>113.497</td>\n",
       "      <td>0.000814</td>\n",
       "      <td>s</td>\n",
       "      <td>t</td>\n",
       "      <td>0.002653</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>100001</td>\n",
       "      <td>160.937</td>\n",
       "      <td>68.768</td>\n",
       "      <td>103.235</td>\n",
       "      <td>48.146</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>3.473</td>\n",
       "      <td>2.078</td>\n",
       "      <td>...</td>\n",
       "      <td>0.725</td>\n",
       "      <td>1.158</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>46.226</td>\n",
       "      <td>0.681042</td>\n",
       "      <td>b</td>\n",
       "      <td>t</td>\n",
       "      <td>2.233584</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>100002</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>162.172</td>\n",
       "      <td>125.953</td>\n",
       "      <td>35.635</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>3.148</td>\n",
       "      <td>9.336</td>\n",
       "      <td>...</td>\n",
       "      <td>2.053</td>\n",
       "      <td>-2.028</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>44.251</td>\n",
       "      <td>0.715742</td>\n",
       "      <td>b</td>\n",
       "      <td>t</td>\n",
       "      <td>2.347389</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>100003</td>\n",
       "      <td>143.905</td>\n",
       "      <td>81.417</td>\n",
       "      <td>80.943</td>\n",
       "      <td>0.414</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>3.310</td>\n",
       "      <td>0.414</td>\n",
       "      <td>...</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-0.000</td>\n",
       "      <td>1.660654</td>\n",
       "      <td>b</td>\n",
       "      <td>t</td>\n",
       "      <td>5.446378</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>100004</td>\n",
       "      <td>175.864</td>\n",
       "      <td>16.915</td>\n",
       "      <td>134.805</td>\n",
       "      <td>16.405</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>3.891</td>\n",
       "      <td>16.405</td>\n",
       "      <td>...</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>1.904263</td>\n",
       "      <td>b</td>\n",
       "      <td>t</td>\n",
       "      <td>6.245333</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 35 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   EventId  DER_mass_MMC  DER_mass_transverse_met_lep  DER_mass_vis  DER_pt_h  \\\n",
       "0   100000       138.470                       51.655        97.827    27.980   \n",
       "1   100001       160.937                       68.768       103.235    48.146   \n",
       "2   100002      -999.000                      162.172       125.953    35.635   \n",
       "3   100003       143.905                       81.417        80.943     0.414   \n",
       "4   100004       175.864                       16.915       134.805    16.405   \n",
       "\n",
       "   DER_deltaeta_jet_jet  DER_mass_jet_jet  DER_prodeta_jet_jet  \\\n",
       "0                  0.91           124.711                2.666   \n",
       "1               -999.00          -999.000             -999.000   \n",
       "2               -999.00          -999.000             -999.000   \n",
       "3               -999.00          -999.000             -999.000   \n",
       "4               -999.00          -999.000             -999.000   \n",
       "\n",
       "   DER_deltar_tau_lep  DER_pt_tot      ...       PRI_jet_leading_eta  \\\n",
       "0               3.064      41.928      ...                     2.150   \n",
       "1               3.473       2.078      ...                     0.725   \n",
       "2               3.148       9.336      ...                     2.053   \n",
       "3               3.310       0.414      ...                  -999.000   \n",
       "4               3.891      16.405      ...                  -999.000   \n",
       "\n",
       "   PRI_jet_leading_phi  PRI_jet_subleading_pt  PRI_jet_subleading_eta  \\\n",
       "0                0.444                 46.062                    1.24   \n",
       "1                1.158               -999.000                 -999.00   \n",
       "2               -2.028               -999.000                 -999.00   \n",
       "3             -999.000               -999.000                 -999.00   \n",
       "4             -999.000               -999.000                 -999.00   \n",
       "\n",
       "   PRI_jet_subleading_phi  PRI_jet_all_pt    Weight  Label  KaggleSet  \\\n",
       "0                  -2.475         113.497  0.000814      s          t   \n",
       "1                -999.000          46.226  0.681042      b          t   \n",
       "2                -999.000          44.251  0.715742      b          t   \n",
       "3                -999.000          -0.000  1.660654      b          t   \n",
       "4                -999.000           0.000  1.904263      b          t   \n",
       "\n",
       "   KaggleWeight  \n",
       "0      0.002653  \n",
       "1      2.233584  \n",
       "2      2.347389  \n",
       "3      5.446378  \n",
       "4      6.245333  \n",
       "\n",
       "[5 rows x 35 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.iloc[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0xa1e8209e8>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE8BJREFUeJzt3X+MHOV5wPHvU1MI5SiQkp6ojWojI1QLS01YQUir6q4N\nYFIMbYQquxaFFLCSCqk/kBqjVKL9o0rShiiJQkvdQGkq4gulNLEdR7RKcyWVEDVuU8yPOHGIE86k\nEIrq6igVcfP0j51zF2vXt3u7vt3Z9/uRTt55Z3bmfTx3z8098+47kZlIksrxQ8PugCRpeZn4Jakw\nJn5JKoyJX5IKY+KXpMKY+CWpMCZ+SSqMiV+SCmPil6TCnDLsDgCce+65uXr16rbrXn31Vc4444zl\n7dAyMK76GMeYwLjqpFNM+/btezkz39Lr/oaa+CNiI7Bx7dq1PPHEE223mZ2dZWpqaln7tRyMqz7G\nMSYwrjrpFFNEfHsp+xtqqSczd2Xm1rPOOmuY3ZCkoljjl6TCmPglqTAmfkkqjIlfkgpj4pekwgw1\n8UfExojYfuTIkWF2Q5KK4nBOSSrMSHxytx+rt32hbfuhD/3iMvdEkurBGr8kFcbEL0mFMfFLUmFM\n/JJUGBO/JBVm4Ik/IqYi4isRcU9ETA16/5Kk/nSV+CPivoh4KSKeOq59Q0QciIiDEbGtak5gHngT\nMDfY7kqS+tXtOP77gU8Cn15oiIgVwN3AFTQT/N6I2Al8JTP/MSImgY8CWwba4y45vl+S2ovM7G7D\niNXA7sy8uFq+HPj9zLyqWr4DIDM/WC2fCnwmM6/vsL+twFaAycnJS2ZmZtoed35+nomJiY792n+4\nt+ke1q8cjU8JLxZXXY1jXOMYExhXnXSKaXp6el9mNnrdXz+f3F0JPN+yPAdcFhHvBq4Czqb5V0Jb\nmbkd2A7QaDSy06PSFnuM2k0druw7ObSl876W0zg+Hg7GM65xjAmMq04GHVM/iT/atGVmPgw83NUO\nWp65K0laHv2M6pkDzm9ZXgW80MsOnKRNkpZfP4l/L3BhRKyp6vmbgJ297MBpmSVp+XU7nHMH8Bhw\nUUTMRcTNmXkUuA14BHgWeDAzn+7l4F7xS9Ly66rGn5mbO7TvAfYs9eDDqPF3GuYJDvWUVAYfxCJJ\nhXGuHkkqjM/claTCWOqRpMJY6pGkwljqkaTCWOqRpML0M1fP2HEqZ0klsNQjSYWx1CNJhXFUjyQV\nxsQvSYUx8UtSYby5K0mFGepwzszcBexqNBq3DrMfi3GYp6RxYqlHkgpj4pekwpj4JakwJn5JKoyJ\nX5IK43BOSSqMwzn74DBPSXVkqUeSCmPil6TCmPglqTAmfkkqjIlfkgpzUkb1RMQZwKPAnZm5+2Qc\nY5Q52kfSKOvqij8i7ouIlyLiqePaN0TEgYg4GBHbWla9H3hwkB2VJA1Gt6We+4ENrQ0RsQK4G7ga\nWAdsjoh1EfFO4BngxQH2U5I0IF2VejLz0YhYfVzzpcDBzHwOICJmgOuACeAMmr8MXouIPZn5g4H1\nWJLUl8jM7jZsJv7dmXlxtXw9sCEzb6mWbwAuy8zbquWbgJc71fgjYiuwFWBycvKSmZmZtsedn59n\nYmKiY7/2H67PdA/rV5517PVicdXVOMY1jjGBcdVJp5imp6f3ZWaj1/31c3M32rQd+y2Smfef6M2Z\nuR3YDtBoNHJqaqrtdrOzs3RaB3BThxupo+jQlqljrxeLq67GMa5xjAmMq04GHVM/iX8OOL9leRXw\nQi87iIiNwMa1a9f20Y36aB3tc/v6o8d+aTnaR9Jy6mcc/17gwohYExGnApuAnb3sIDN3ZebWs846\na/GNJUkD0e1wzh3AY8BFETEXETdn5lHgNuAR4Fngwcx8upeDOy2zJC2/bkf1bO7QvgfYs9SD131a\nZkmqIx/EIkmFGWrit8YvScvPSdokqTBDffRiacM5O3FSN0nLyVKPJBXGm7uSVBiv+CWpMN7claTC\nDPXmrk7Mm76STgZr/JJUGGv8klQYa/ySVBgTvyQVxsQvSYVxyoYacrSPpH54c1eSCmOpR5IKY+KX\npMKY+CWpMCZ+SSqMo3rGiKN9JHXDUT2SVBhLPZJUGBO/JBXG+fgLYO1fUiuv+CWpMCZ+SSqMiV+S\nCjPwxB8RPxUR90TEQxHxvkHvX5LUn64Sf0TcFxEvRcRTx7VviIgDEXEwIrYBZOazmfle4FeAxuC7\nLEnqR7ejeu4HPgl8eqEhIlYAdwNXAHPA3ojYmZnPRMS1wLbqPRpRjvaRytTVFX9mPgq8clzzpcDB\nzHwuM18HZoDrqu13ZuY7gC2D7KwkqX+Rmd1tGLEa2J2ZF1fL1wMbMvOWavkG4DLgIeDdwGnAk5l5\nd4f9bQW2AkxOTl4yMzPT9rjz8/NMTEx07Nf+w0e66v+omTwdXnxt2L1ob/3KpU+hsdj5qqNxjAmM\nq046xTQ9Pb0vM3suqffzAa5o05aZOQvMLvbmzNwObAdoNBo5NTXVdrvZ2Vk6rQO4qUO5YtTdvv4o\nd+0fzc/PHdoyteT3Lna+6mgcYwLjqpNBx9TPqJ454PyW5VXAC73sICI2RsT2I0fqedUuSXXUT+Lf\nC1wYEWsi4lRgE7BzMN2SJJ0sXdUaImIHMAWcGxFzwJ2ZeW9E3AY8AqwA7svMp3s5eGbuAnY1Go1b\ne+u2TqZOo33AET/SOOgq8Wfm5g7te4A9A+2RJOmkGuqUDdb4JWn5DXVYiaWe+vFDX1L9ecUvSYXx\nmbuSVBinZZakwgy1xh8RG4GNa9euHWY3NAALtf/b1x99w6eprf1Lo8dSjyQVxlKPJBXGxC9JhXE4\npyQVxhq/JBVmNCeE19jwk77S6LHGL0mFMfFLUmG8uStJhXF2Tg2FtX9peCz1SFJhTPySVBgTvyQV\nxnH8GinW/qWTz1E9klQYp2yQpMJY45ekwpj4JakwJn5JKoyJX5IK43BO1YLDPKXB8YpfkgpzUhJ/\nRPxSRPx5RHw+Iq48GceQJC1N16WeiLgPuAZ4KTMvbmnfAHwcWAF8KjM/lJmfAz4XEecAHwH+brDd\nlposAUm96+WK/35gQ2tDRKwA7gauBtYBmyNiXcsmv1etlySNiK4Tf2Y+CrxyXPOlwMHMfC4zXwdm\ngOui6cPAFzPzXwbXXUlSvyIzu984YjWwe6HUExHXAxsy85Zq+QbgMuDrwI3AXuCrmXlPm31tBbYC\nTE5OXjIzM9P2mPPz80xMTHTs0/7D9ZznZ/J0ePG1Yfdi8EYlrvUrBzcNyGLfg3VlXPXRKabp6el9\nmdnodX/9DueMNm2ZmZ8APnGiN2bm9oj4LrDxzDPPvGRqaqrtdrOzs3RaB3BThxrvqLt9/VHu2j9+\no2lHJa5DW6YGtq/FvgfryrjqY9Ax9fsTOgec37K8Cnih2zf76EWdLN70lTrrN/HvBS6MiDXAYWAT\n8KvdvjkiNgIb165d22c3pO74C0Hq4eZuROwAHgMuioi5iLg5M48CtwGPAM8CD2bm093u02mZJWn5\ndX3Fn5mbO7TvAfYs5eBe8UvS8vNBLJJUGOfqkaTC+MxdSSqMpR5JKoylHkkqzFA/YumoHo0Kx/er\nJJZ6JKkww59URRphq7d9gdvXH207J5R/DaiurPFLUmEczilJhbHGL0mFsdQjSYXx5q40YA4N1aiz\nxi9JhbHGL0mFscYvSYUx8UtSYUz8klQYE78kFcbEL0mFcTinJBXG4ZySVBhLPZJUGKdskJao09QM\n0qjzil+SCmPil6TCmPglqTADT/wRcUFE3BsRDw1635Kk/nV1czci7gOuAV7KzItb2jcAHwdWAJ/K\nzA9l5nPAzSZ+6Y2cp1+jotsr/vuBDa0NEbECuBu4GlgHbI6IdQPtnSRp4Lq64s/MRyNi9XHNlwIH\nqyt8ImIGuA54ZpAdlMadfwlouUVmdrdhM/HvXij1RMT1wIbMvKVavgG4DLgT+EPgCprlnw922N9W\nYCvA5OTkJTMzM22POz8/z8TERMd+7T9cz+keJk+HF18bdi8GbxzjGlZM61ee3E+0L/azVVfjGFen\nmKanp/dlZqPX/fXzAa5o05aZ+R/Aexd7c2ZuB7YDNBqNnJqaarvd7OwsndYB3FTTD9Hcvv4od+0f\nv8/PjWNcw4rp0Japk7r/xX626moc4xp0TP2M6pkDzm9ZXgW80MsOnKRNkpZfP4l/L3BhRKyJiFOB\nTcDOXnbgJG2StPy6Hc65A5gCzo2IOeDOzLw3Im4DHqE5nPO+zHy6l4NHxEZg49q1a3vrtVQwbwar\nX92O6tncoX0PsGepB8/MXcCuRqNx61L3IUnqjQ9ikaTC+CAWSSqMk7RJUmEs9UhSYSz1SFJhLPVI\nUmEs9UhSYSz1SFJhLPVIUmGGOo2iUzZInXWamkHql6UeSSqMpR5JKoyJX5IKY+KXpMJ4c1cq1P7D\nR9o+utR5/cefN3clqTCWeiSpMCZ+SSqMiV+SCmPil6TCmPglqTAO55TGRKe5fXodnjnIOYIGdWyH\nmA6WwzklqTCWeiSpMCZ+SSqMiV+SCmPil6TCmPglqTADH84ZEWcAfwK8Dsxm5gODPoYkaem6uuKP\niPsi4qWIeOq49g0RcSAiDkbEtqr53cBDmXkrcO2A+ytJ6lO3pZ77gQ2tDRGxArgbuBpYB2yOiHXA\nKuD5arP/HUw3JUmD0lXiz8xHgVeOa74UOJiZz2Xm68AMcB0wRzP5d71/SdLyiczsbsOI1cDuzLy4\nWr4e2JCZt1TLNwCXAe8HPgn8D/BPnWr8EbEV2AowOTl5yczMTNvjzs/PMzEx0bFf+w8f6ar/o2by\ndHjxtWH3YvDGMa5xjAmWJ671K9t/Kr/Xn9tO+2nnRDmj03F72f+JLCUfdXPsTjFNT0/vy8xGr8fs\n5+ZutGnLzHwVeM9ib87M7cB2gEajkVNTU223m52dpdM6oO2j4+rg9vVHuWv/UKdKOinGMa5xjAmW\nJ65DW6batvf6c9tpP+2cKGd0Om4v+z+RpeSjbo69WB7sVT+lmDng/JblVcALvewgIjZGxPYjR+p5\n1S5JddRP4t8LXBgRayLiVGATsLOXHThJmyQtv26Hc+4AHgMuioi5iLg5M48CtwGPAM8CD2bm070c\n3Ct+SVp+XRX4MnNzh/Y9wJ6lHjwzdwG7Go3GrUvdhySpN0MdbukVvyQtPx/EIkmF8YpfkgrjFb8k\nFabrT+6e1E5EfA/4dofV5wIvL2N3lotx1cc4xgTGVSedYvrJzHxLrzsbicR/IhHxxFI+kjzqjKs+\nxjEmMK46GXRMTqImSYUx8UtSYeqQ+LcPuwMniXHVxzjGBMZVJwONaeRr/JKkwarDFb8kaYBGOvF3\neKbvyIuI8yPiyxHxbEQ8HRG/WbW/OSL+PiK+Uf17TtUeEfGJKs4nI+Jtw43gxCJiRUT8a0TsrpbX\nRMTjVVyfrWZrJSJOq5YPVutXD7PfJxIRZ0fEQxHxteq8XV738xURv119/z0VETsi4k11PFftnvm9\nlHMTETdW238jIm4cRiytOsT1x9X34JMR8bcRcXbLujuquA5ExFUt7b3nycwcyS9gBfBN4ALgVODf\ngHXD7leXfT8PeFv1+kzg6zSfS/xHwLaqfRvw4er1u4Av0ny4zduBx4cdwyLx/Q7wGZpPZAN4ENhU\nvb4HeF/1+jeAe6rXm4DPDrvvJ4jpL4FbqtenAmfX+XwBK4FvAae3nKOb6niugJ8D3gY81dLW07kB\n3gw8V/17TvX6nBGM60rglOr1h1viWlflwNOANVVuXLHUPDn0k3qC/5TLgUdalu8A7hh2v5YYy+eB\nK4ADwHlV23nAger1nwGbW7Y/tt2ofdF84M6XgJ8Hdlc/YC+3fLMeO280p+y+vHp9SrVdDDuGNjH9\naJUk47j22p6vKvE/XyW6U6pzdVVdzxWw+rgE2dO5ATYDf9bS/obtRiWu49b9MvBA9foN+W/hfC01\nT45yqWfhG3fBXNVWK9WfzG8FHgcmM/O7ANW/P15tVqdYPwb8LvCDavnHgP/M5vMZ4I19PxZXtf5I\ntf2ouQD4HvAXVQnrUxFxBjU+X5l5GPgI8B3guzT/7/dR/3O1oNdzM/LnrI1fp/nXCww4rlFO/G2f\n6bvsvehDREwAfwP8Vmb+14k2bdM2crFGxDXAS5m5r7W5zabZxbpRcgrNP7n/NDPfCrxKs3zQycjH\nVdW8r6NZFvgJ4Azg6jab1u1cLaZTHLWKLyI+ABwFHlhoarPZkuMa5cTf9zN9hykifphm0n8gMx+u\nml+MiPOq9ecBL1XtdYn1Z4BrI+IQMEOz3PMx4OyIWHioT2vfj8VVrT8LeGU5O9ylOWAuMx+vlh+i\n+YugzufrncC3MvN7mfl94GHgHdT/XC3o9dzU4ZwBzZvQwDXAlqzqNww4rlFO/H0/03dYIiKAe4Fn\nM/OjLat2AgujCW6kWftfaP+1akTC24EjC3/GjpLMvCMzV2Xmaprn4x8ycwvwZeD6arPj41qI9/pq\n+5G7ysrMfweej4iLqqZfAJ6h3ufrO8DbI+JHqu/HhZhqfa5a9HpuHgGujIhzqr+GrqzaRkpEbADe\nD1ybmf/dsmonsKkafbUGuBD4Z5aaJ4d9c2ORGx/vojki5pvAB4bdnx76/bM0/9x6Evhq9fUumjXT\nLwHfqP59c7V9AHdXce4HGsOOoYsYp/j/UT0XVN+EB4G/Bk6r2t9ULR+s1l8w7H6fIJ6fBp6oztnn\naI78qPX5Av4A+BrwFPBXNEeE1O5cATto3qf4Ps0r3JuXcm5o1swPVl/vGdG4DtKs2S/kjXtatv9A\nFdcB4OqW9p7zpJ/claTCjHKpR5J0Epj4JakwJn5JKoyJX5IKY+KXpMKY+CWpMCZ+SSqMiV+SCvN/\nSV7l/WdbTtwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa1e70b048>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGd9JREFUeJzt3X+MVfd55/H3J/hHkNME/CN3KaCFbqerkKASZ2QjeVXd\ntVt7wNVCpFjFi2ycsjttBNtEmu4ap384tcMKr5a48cpxO4mnhigNRU4sUIqXso6vokixDXYIGFOX\nCZ4NE1izKdjxJKq94zz7x/kOuZlz78ydmftr5n5e0tU99znfc+730Zk7z/2ec+45igjMzMzKvafV\nHTAzs/bj4mBmZjkuDmZmluPiYGZmOS4OZmaW4+JgZmY5Lg5mZpbj4mBmZjkuDmZmlnNZqzswXdde\ne20sW7YMgJ/97GdcddVVre1QC3Vy/s7duXeameb+4osv/iQirpus3awtDsuWLePIkSMAlEolisVi\nazvUQp2cv3MvtrobLeHci9NeXtL/rqWddyuZmVmOi4OZmeW4OJiZWY6Lg5mZ5bg4mJlZjouDmZnl\nuDiYmVmOi4OZmeW4OJiZWc6s/YW02Zhl2/6uYnxox+1N7onZ3OGRg5mZ5dQ8cpA0DzgC/Dgifl/S\ncmAPcDXwEnBXRLwj6UpgN/Ax4J+AP4iIobSO+4DNwLvAn0TEwRTvAb4IzAO+EhE76pSfzRHVRgd9\nK0fxANis/qbyqfo0cBJ4f3r9EPBwROyR9Jdk//QfS88XI+I3JW1I7f5A0gpgA/Bh4NeB/yXpt9K6\nHgV+DxgGDkvaHxGvzDA363De3WQ2fTXtVpK0BLgd+Ep6LeBm4MnUZBewPk2vS69J829J7dcBeyLi\n7Yh4DRgEbkiPwYg4HRHvkI1G1s00MTMzm75ajzn8BfBfgF+k19cAb0TEaHo9DCxO04uBMwBp/pup\n/aX4uGWqxc3MrEUm3a0k6feB8xHxoqTiWLhC05hkXrV4pQIVFWJI6gV6AQqFAqVSCYCRkZFL052o\nE/LPji3kFeZXn1fN//javorxlYs/MOV+tVInbPdqnHup4e9TyzGHm4B/J2kt8F6yYw5/ASyQdFka\nHSwBzqb2w8BSYFjSZcAHgAtl8THly1SL/4qI6Af6Abq7u2PshhedfOMP6Iz875nggPTO4/U5ID20\nsViX9TRLJ2z3apx7seHvM+lupYi4LyKWRMQysgPK346IjcCzwCdSs03A2Nex/ek1af63IyJSfIOk\nK9OZTl3AC8BhoEvScklXpPfYX5fszMxsWmbyleteYI+kzwPfBx5P8ceBr0oaJBsxbACIiBOS9gKv\nAKPAloh4F0DSVuAg2amsAxFxYgb9slms2hlGZtZcUyoOEVECSmn6NNmZRuPb/DNwR5XltwPbK8QP\nAAem0hezevOpr2a/5F9Im5lZjouDmZnluDiYmVmOi4OZmeW4OJiZWY4vZ2k2CZ/FZJ3IIwczM8tx\ncTAzsxwXBzMzy/ExB2sJXybDrL155GBmZjkuDmZmluPiYGZmOS4OZmaW4+JgZmY5Lg5mZpYzaXGQ\n9F5JL0j6gaQTkv48xZ+Q9Jqko+mxKsUl6RFJg5KOSbq+bF2bJJ1Kj01l8Y9JOp6WeUSSGpGsmZnV\nppbfObwN3BwRI5IuB74r6ek07z9HxJPj2q8huz90F3Aj8Bhwo6SrgfuBbiCAFyXtj4iLqU0v8BzZ\nHeF6gKcxM7OWmLQ4REQAI+nl5ekREyyyDtidlntO0gJJi4AicCgiLgBIOgT0SCoB74+I76X4bmA9\nLg7W5nxBPpvLajrmIGmepKPAebJ/8M+nWdvTrqOHJV2ZYouBM2WLD6fYRPHhCnEzM2uRmi6fERHv\nAqskLQCekvQR4D7g/wBXAP3AvcADQKXjBTGNeI6kXrLdTxQKBUqlEgAjIyOXpjvRbMy/b+VoXdZT\nmF+/ddVLs7bFbNzu9eLcSw1/nyldWyki3ki7gXoi4r+n8NuS/hr40/R6GFhattgS4GyKF8fFSym+\npEL7Su/fT1aI6O7ujmIxW12pVGJsuhPNxvzvqdO1lfpWjrLzeHtdImxoY7Ep7zMbt3u9OPdiw9+n\nlrOVrksjBiTNB34X+Id0HIF0ZtF64OW0yH7g7nTW0mrgzYg4BxwEbpW0UNJC4FbgYJr3lqTVaV13\nA/vqm6aZmU1FLV+5FgG7JM0jKyZ7I+Jbkr4t6Tqy3UJHgT9O7Q8Aa4FB4OfAJwEi4oKkB4HDqd0D\nYwengU8BTwDzyQ5E+2C0mVkL1XK20jHgoxXiN1dpH8CWKvMGgIEK8SPARybri5mZNYd/IW1mZjku\nDmZmluPiYGZmOS4OZmaW4+JgZmY5Lg5mZpbj4mBmZjkuDmZmluPiYGZmOS4OZmaW016XszSbA3wT\nIJsLPHIwM7Mcjxysoap9izaz9uaRg5mZ5bg4mJlZjouDmZnl1HKb0PdKekHSDySdkPTnKb5c0vOS\nTkn6W0lXpPiV6fVgmr+sbF33pfirkm4ri/ek2KCkbfVP08zMpqKWkcPbwM0R8dvAKqAn3Rv6IeDh\niOgCLgKbU/vNwMWI+E3g4dQOSSuADcCHgR7gS5LmpduPPgqsAVYAd6a2ZmbWIpMWh8iMpJeXp0cA\nNwNPpvguYH2aXpdek+bfIkkpvici3o6I18juMX1DegxGxOmIeAfYk9qamVmL1HTMIX3DPwqcBw4B\nPwTeiIjR1GQYWJymFwNnANL8N4FryuPjlqkWNzOzFqnpdw4R8S6wStIC4CngQ5WapWdVmVctXqlA\nRYUYknqBXoBCoUCpVAJgZGTk0nQnauf8+1aOTt5oBgrzG/8e9VLvbdTO273RnHup4e8zpR/BRcQb\nkkrAamCBpMvS6GAJcDY1GwaWAsOSLgM+AFwoi48pX6ZafPz79wP9AN3d3VEsFoHsQzc23YnaOf97\nGvwjuL6Vo+w8Pjt+yzm0sVjX9bXzdm80515s+PvUcrbSdWnEgKT5wO8CJ4FngU+kZpuAfWl6f3pN\nmv/tiIgU35DOZloOdAEvAIeBrnT20xVkB6331yM5MzObnlq+ci0CdqWzit4D7I2Ib0l6Bdgj6fPA\n94HHU/vHga9KGiQbMWwAiIgTkvYCrwCjwJa0uwpJW4GDwDxgICJO1C1DMzObskmLQ0QcAz5aIX6a\n7Eyj8fF/Bu6osq7twPYK8QPAgRr6a2ZmTeBfSJuZWY6Lg5mZ5bg4mJlZjouDmZnluDiYmVmOi4OZ\nmeW4OJiZWY6Lg5mZ5cyOi9KYzQHLqlxnamjH7U3uidnkPHIwM7McFwczM8txcTAzsxwXBzMzy3Fx\nMDOzHBcHMzPLcXEwM7OcSX/nIGkpsBv4F8AvgP6I+KKkzwH/Efi/qeln0017kHQfsBl4F/iTiDiY\n4j3AF8nu+PaViNiR4suBPcDVwEvAXRHxTr2StMardg6/mc1OtYwcRoG+iPgQsBrYImlFmvdwRKxK\nj7HCsILs1qAfBnqAL0mal24z+iiwBlgB3Fm2nofSurqAi2SFxczMWmTS4hAR5yLipTT9FnASWDzB\nIuuAPRHxdkS8BgyS3U70BmAwIk6nUcEeYJ0kATcDT6bldwHrp5uQmZnN3JSOOUhaRnY/6edTaKuk\nY5IGJC1MscXAmbLFhlOsWvwa4I2IGB0XNzOzFqn52kqS3gd8A/hMRPxU0mPAg0Ck553AHwKqsHhQ\nuRDFBO0r9aEX6AUoFAqUSiUARkZGLk13onbIv2/l6OSNGqAwv3XvXS/T3XbtsN1bxbmXGv4+NRUH\nSZeTFYavRcQ3ASLi9bL5Xwa+lV4OA0vLFl8CnE3TleI/ARZIuiyNHsrb/4qI6Af6Abq7u6NYLALZ\nh2tsuhO1Q/73tOiAdN/KUXYen93XjxzaWJzWcu2w3VvFuRcb/j6T7lZKxwQeB05GxBfK4ovKmn0c\neDlN7wc2SLoynYXUBbwAHAa6JC2XdAXZQev9ERHAs8An0vKbgH0zS8vMzGailq9cNwF3AcclHU2x\nz5KdbbSKbBfQEPBHABFxQtJe4BWyM522RMS7AJK2AgfJTmUdiIgTaX33AnskfR74PlkxMjOzFpm0\nOETEd6l8XODABMtsB7ZXiB+otFxEnCY7m8nMzNqAfyFtZmY5Lg5mZpbj4mBmZjkuDmZmluPiYGZm\nOS4OZmaW4+JgZmY5Lg5mZpYzuy9KYzYHVLtR0tCO25vcE7Nf8sjBzMxyXBzMzCzHxcHMzHJcHMzM\nLMfFwczMclwczMwsx8XBzMxyXBzMzCynlntIL5X0rKSTkk5I+nSKXy3pkKRT6XlhikvSI5IGJR2T\ndH3Zujal9qckbSqLf0zS8bTMI+m+1WZm1iK1jBxGgb6I+BCwGtgiaQWwDXgmIrqAZ9JrgDVAV3r0\nAo9BVkyA+4EbyW4Jev9YQUltesuW65l5amZmNl2TFoeIOBcRL6Xpt4CTwGJgHbArNdsFrE/T64Dd\nkXkOWCBpEXAbcCgiLkTEReAQ0JPmvT8ivhcRAewuW5eZmbXAlK6tJGkZ8FHgeaAQEecgKyCSPpia\nLQbOlC02nGITxYcrxCu9fy/ZCINCoUCpVAJgZGTk0nQnaof8+1aOtuR9C/Nb996NNtk2bYft3irO\nvdTw96m5OEh6H/AN4DMR8dMJDgtUmhHTiOeDEf1AP0B3d3cUi0Ug+xCNTXeidsj/nioXj2u0vpWj\n7Dw+N68fObSxOOH8dtjureLciw1/n5rOVpJ0OVlh+FpEfDOFX0+7hEjP51N8GFhatvgS4Owk8SUV\n4mZm1iK1nK0k4HHgZER8oWzWfmDsjKNNwL6y+N3prKXVwJtp99NB4FZJC9OB6FuBg2neW5JWp/e6\nu2xdZmbWArWMx28C7gKOSzqaYp8FdgB7JW0GfgTckeYdANYCg8DPgU8CRMQFSQ8Ch1O7ByLiQpr+\nFPAEMB94Oj2sDVW794CZzS2TFoeI+C6VjwsA3FKhfQBbqqxrABioED8CfGSyvpiZWXP4F9JmZpbj\n4mBmZjkuDmZmluPiYGZmOXPz10Nmc0C1M8OGdtze5J5YJ/LIwczMclwczMwsx8XBzMxyXBzMzCzH\nxcHMzHJcHMzMLMfFwczMclwczMwsx8XBzMxyXBzMzCzHxcHMzHJquU3ogKTzkl4ui31O0o8lHU2P\ntWXz7pM0KOlVSbeVxXtSbFDStrL4cknPSzol6W8lXVHPBM3MbOpqGTk8AfRUiD8cEavS4wCApBXA\nBuDDaZkvSZonaR7wKLAGWAHcmdoCPJTW1QVcBDbPJCEzM5u5SYtDRHwHuDBZu2QdsCci3o6I18ju\nI31DegxGxOmIeAfYA6yTJOBm4Mm0/C5g/RRzMDOzOpvJJbu3SrobOAL0RcRFYDHwXFmb4RQDODMu\nfiNwDfBGRIxWaJ8jqRfoBSgUCpRKJQBGRkYuTXeiZubft3J08kZNVJjffn1qNP/dO/dm5D7d4vAY\n8CAQ6Xkn8IeAKrQNKo9QYoL2FUVEP9AP0N3dHcViEcg+LGPTnaiZ+d9T5R4DrdK3cpSdxzvrtiRD\nG4tAZ//dO/diw99nWp+qiHh9bFrSl4FvpZfDwNKypkuAs2m6UvwnwAJJl6XRQ3l7MzNrkWkVB0mL\nIuJcevlxYOxMpv3A30j6AvDrQBfwAtkIoUvScuDHZAet/31EhKRngU+QHYfYBOybbjJmnWDsDnF9\nK0dzIznfJc7qZdLiIOnrQBG4VtIwcD9QlLSKbBfQEPBHABFxQtJe4BVgFNgSEe+m9WwFDgLzgIGI\nOJHe4l5gj6TPA98HHq9bdmZmNi2TFoeIuLNCuOo/8IjYDmyvED8AHKgQP012NpOZmbUJ/0LazMxy\nOus0D6vZsjY7K8nMmssjBzMzy3FxMDOzHBcHMzPLcXEwM7McFwczM8txcTAzsxwXBzMzy3FxMDOz\nHBcHMzPLcXEwM7McXz7DbA6pdtkTX8rbpsojBzMzy3FxMDOzHBcHMzPLmbQ4SBqQdF7Sy2WxqyUd\nknQqPS9McUl6RNKgpGOSri9bZlNqf0rSprL4xyQdT8s8Ikn1TtLMzKamlpHDE0DPuNg24JmI6AKe\nSa8B1pDdN7oL6AUeg6yYkN1e9Eayu77dP1ZQUpvesuXGv5eZmTXZpMUhIr4DXBgXXgfsStO7gPVl\n8d2ReQ5YIGkRcBtwKCIuRMRF4BDQk+a9PyK+FxEB7C5bl5mZtch0jzkUIuIcQHr+YIovBs6UtRtO\nsYniwxXiZmbWQvX+nUOl4wUxjXjllUu9ZLugKBQKlEolAEZGRi5Nd6JG5N+3crSu62uUwvzZ09d6\nm0ruc+3z0cmf+WblPt3i8LqkRRFxLu0aOp/iw8DSsnZLgLMpXhwXL6X4kgrtK4qIfqAfoLu7O4rF\nbJWlUomx6U7UiPzvmSX3kO5bOcrO4535W86p5D60sdjYzjRZJ3/mm5X7dHcr7QfGzjjaBOwri9+d\nzlpaDbyZdjsdBG6VtDAdiL4VOJjmvSVpdTpL6e6ydZmZWYtM+rVD0tfJvvVfK2mY7KyjHcBeSZuB\nHwF3pOYHgLXAIPBz4JMAEXFB0oPA4dTugYgYO8j9KbIzouYDT6eHmdWRL6thUzVpcYiIO6vMuqVC\n2wC2VFnPADBQIX4E+Mhk/bDGqPZPw8w6m38hbWZmOS4OZmaW4+JgZmY5Lg5mZpbj4mBmZjkuDmZm\nluPiYGZmOS4OZmaW05kXpTEzwL+ctuo8cjAzsxwXBzMzy3FxMDOzHBcHMzPL8QHpDuGrr5rZVLg4\nmFmOz2Iy71YyM7OcGRUHSUOSjks6KulIil0t6ZCkU+l5YYpL0iOSBiUdk3R92Xo2pfanJG2q9n5m\nZtYc9Rg5/NuIWBUR3en1NuCZiOgCnkmvAdYAXenRCzwGWTEhu/XojcANwP1jBcXMzFqjEccc1pHd\ncxpgF1AC7k3x3elWos9JWiBpUWp7aOye0pIOAT3A1xvQtznPB57NrB5mWhwC+HtJAfxVRPQDhYg4\nBxAR5yR9MLVdDJwpW3Y4xarFzazN+EB155hpcbgpIs6mAnBI0j9M0FYVYjFBPL8CqZdslxSFQoFS\nqQTAyMjIpelOVJ5/38rR1namyQrzOy/nMe2Ue7M/f538mW9W7jMqDhFxNj2fl/QU2TGD1yUtSqOG\nRcD51HwYWFq2+BLgbIoXx8VLVd6vH+gH6O7ujmIxW6xUKjE23YnK87+nw3Yr9a0cZefxzjwju51y\nH9pYbOr7dfJnvlm5T/svS9JVwHsi4q00fSvwALAf2ATsSM/70iL7ga2S9pAdfH4zFZCDwH8tOwh9\nK3DfdPtlZs3n3U1zz0y+dhSApySNredvIuJ/SjoM7JW0GfgRcEdqfwBYCwwCPwc+CRARFyQ9CBxO\n7R4YOzhtZmatMe3iEBGngd+uEP8n4JYK8QC2VFnXADAw3b50ovJvan0rRztud5LNDh5RzF7+hbSZ\nmeW4OJiZWU57nOpgVflHbTYXeXdT+3NxMLO24aLRPrxbyczMcjxyMLO2N35EMXaGnkcUjePi0AZ8\nXMHM2o2Lg5nNWtP5YuXRRm1cHJrIIwQzmy1cHMyso0z1S1qnjjRcHBrAIwSzuaNTT691cTAzm4a5\nPgJxcZgBjxDMrFazbQTi4mBm1kJT/ZL5RM9VDerJr3JxqIFHCGbWaXz5DDMzy2mb4iCpR9KrkgYl\nbWt1f8zMOllb7FaSNA94FPg9YBg4LGl/RLzSzH5495GZWaZdRg43AIMRcToi3gH2AOta3Cczs47V\nLsVhMXCm7PVwipmZWQsoIlrdByTdAdwWEf8hvb4LuCEi/tO4dr1Ab3r5r4FX0/S1wE+a1N121Mn5\nO/fO5Nyn719GxHWTNWqLYw5kI4WlZa+XAGfHN4qIfqB/fFzSkYjoblz32lsn5+/cnXunaVbu7bJb\n6TDQJWm5pCuADcD+FvfJzKxjtcXIISJGJW0FDgLzgIGIONHibpmZday2KA4AEXEAODDNxXO7mjpM\nJ+fv3DuTc2+wtjggbWZm7aVdjjmYmVkbmfXFodMuuyFpSNJxSUclHUmxqyUdknQqPS9sdT/rQdKA\npPOSXi6LVcxVmUfS38ExSde3ruf1USX/z0n6cdr+RyWtLZt3X8r/VUm3tabXMydpqaRnJZ2UdELS\np1O8I7b9BPk3d9tHxKx9kB28/iHwG8AVwA+AFa3uV4NzHgKuHRf7b8C2NL0NeKjV/axTrr8DXA+8\nPFmuwFrgaUDAauD5Vve/Qfl/DvjTCm1XpL//K4Hl6XMxr9U5TDPvRcD1afrXgH9M+XXEtp8g/6Zu\n+9k+cvBlNzLrgF1pehewvoV9qZuI+A5wYVy4Wq7rgN2ReQ5YIGlRc3raGFXyr2YdsCci3o6I14BB\nss/HrBMR5yLipTT9FnCS7IoJHbHtJ8i/moZs+9leHDrxshsB/L2kF9MvxgEKEXEOsj8s4IMt613j\nVcu1k/4WtqbdJwNluxDnZP6SlgEfBZ6nA7f9uPyhidt+thcHVYjN9dOvboqI64E1wBZJv9PqDrWJ\nTvlbeAz4V8Aq4BywM8XnXP6S3gd8A/hMRPx0oqYVYrM6d6iYf1O3/WwvDjVddmMuiYiz6fk88BTZ\n8PH1sWF0ej7fuh42XLVcO+JvISJej4h3I+IXwJf55e6DOZW/pMvJ/jF+LSK+mcIds+0r5d/sbT/b\ni0NHXXZD0lWSfm1sGrgVeJks502p2SZgX2t62BTVct0P3J3OXFkNvDm2C2IuGbcv/eNk2x+y/DdI\nulLScqALeKHZ/asHSQIeB05GxBfKZnXEtq+Wf9O3fauPzNfhyP5asqP5PwT+rNX9aXCuv0F2VsIP\ngBNj+QLXAM8Ap9Lz1a3ua53y/TrZ8Pn/kX072lwtV7Kh9aPp7+A40N3q/jco/6+m/I6lfwqLytr/\nWcr/VWBNq/s/g7z/DdlukWPA0fRY2ynbfoL8m7rt/QtpMzPLme27lczMrAFcHMzMLMfFwczMclwc\nzMwsx8XBzMxyXBzMzCzHxcHMzHJcHMzMLOf/A3I4B+Wir7oqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa1fbfbe80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.PRI_jet_leading_pt[df.PRI_jet_leading_pt>0].hist(bins=50)\n",
    "plt.yscale('log')\n",
    "\n",
    "f=plt.figure()\n",
    "df.DER_mass_MMC[(df.DER_mass_MMC>0)&(df.DER_mass_MMC<250)].hist(bins=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is more information about the variables in the documentation. The variables that start with **DER** are derived quantities, determined by the physicists performing the analysis as variables that discriminate signal from backround. On the other hand, those that start with **PRI** are considered to be primary variables, from which the derived variables are calculated. They themselves generally do not provide much discrimination, but one if the ideas suggested by deep networks is that they can determine the necessary features from the primary variables, potentially even finding variables that the physicists did not consider. *EventId* identifies the event but is not a \"feature.\" The *Weight* is the event weight so that the sum of weights of all signal events should produce the signal yield expected to be observed in 2012, and the sum of weights of all background events should produce the backgroudn yield. Note that the weight varies event to event, because different background and signal processes contribute to the background and signal sets. *Label* indicates if it is a signal or background event. Ignore the *Kaggle* variables--they are only used if you want to reproduce exactly what was used in the Challenge. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# map y values to integers\n",
    "df['Label'] = df['Label'].map({'b':0, 's':1})\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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>EventId</th>\n",
       "      <th>DER_mass_MMC</th>\n",
       "      <th>DER_mass_transverse_met_lep</th>\n",
       "      <th>DER_mass_vis</th>\n",
       "      <th>DER_pt_h</th>\n",
       "      <th>DER_deltaeta_jet_jet</th>\n",
       "      <th>DER_mass_jet_jet</th>\n",
       "      <th>DER_prodeta_jet_jet</th>\n",
       "      <th>DER_deltar_tau_lep</th>\n",
       "      <th>DER_pt_tot</th>\n",
       "      <th>...</th>\n",
       "      <th>PRI_jet_leading_eta</th>\n",
       "      <th>PRI_jet_leading_phi</th>\n",
       "      <th>PRI_jet_subleading_pt</th>\n",
       "      <th>PRI_jet_subleading_eta</th>\n",
       "      <th>PRI_jet_subleading_phi</th>\n",
       "      <th>PRI_jet_all_pt</th>\n",
       "      <th>Weight</th>\n",
       "      <th>Label</th>\n",
       "      <th>KaggleSet</th>\n",
       "      <th>KaggleWeight</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>100000</td>\n",
       "      <td>138.470</td>\n",
       "      <td>51.655</td>\n",
       "      <td>97.827</td>\n",
       "      <td>27.980</td>\n",
       "      <td>0.91</td>\n",
       "      <td>124.711</td>\n",
       "      <td>2.666</td>\n",
       "      <td>3.064</td>\n",
       "      <td>41.928</td>\n",
       "      <td>...</td>\n",
       "      <td>2.150</td>\n",
       "      <td>0.444</td>\n",
       "      <td>46.062</td>\n",
       "      <td>1.24</td>\n",
       "      <td>-2.475</td>\n",
       "      <td>113.497</td>\n",
       "      <td>0.000814</td>\n",
       "      <td>1</td>\n",
       "      <td>t</td>\n",
       "      <td>0.002653</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>100001</td>\n",
       "      <td>160.937</td>\n",
       "      <td>68.768</td>\n",
       "      <td>103.235</td>\n",
       "      <td>48.146</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>3.473</td>\n",
       "      <td>2.078</td>\n",
       "      <td>...</td>\n",
       "      <td>0.725</td>\n",
       "      <td>1.158</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>46.226</td>\n",
       "      <td>0.681042</td>\n",
       "      <td>0</td>\n",
       "      <td>t</td>\n",
       "      <td>2.233584</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>100002</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>162.172</td>\n",
       "      <td>125.953</td>\n",
       "      <td>35.635</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>3.148</td>\n",
       "      <td>9.336</td>\n",
       "      <td>...</td>\n",
       "      <td>2.053</td>\n",
       "      <td>-2.028</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>44.251</td>\n",
       "      <td>0.715742</td>\n",
       "      <td>0</td>\n",
       "      <td>t</td>\n",
       "      <td>2.347389</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>100003</td>\n",
       "      <td>143.905</td>\n",
       "      <td>81.417</td>\n",
       "      <td>80.943</td>\n",
       "      <td>0.414</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>3.310</td>\n",
       "      <td>0.414</td>\n",
       "      <td>...</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-0.000</td>\n",
       "      <td>1.660654</td>\n",
       "      <td>0</td>\n",
       "      <td>t</td>\n",
       "      <td>5.446378</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>100004</td>\n",
       "      <td>175.864</td>\n",
       "      <td>16.915</td>\n",
       "      <td>134.805</td>\n",
       "      <td>16.405</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>3.891</td>\n",
       "      <td>16.405</td>\n",
       "      <td>...</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>-999.00</td>\n",
       "      <td>-999.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>1.904263</td>\n",
       "      <td>0</td>\n",
       "      <td>t</td>\n",
       "      <td>6.245333</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 35 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   EventId  DER_mass_MMC  DER_mass_transverse_met_lep  DER_mass_vis  DER_pt_h  \\\n",
       "0   100000       138.470                       51.655        97.827    27.980   \n",
       "1   100001       160.937                       68.768       103.235    48.146   \n",
       "2   100002      -999.000                      162.172       125.953    35.635   \n",
       "3   100003       143.905                       81.417        80.943     0.414   \n",
       "4   100004       175.864                       16.915       134.805    16.405   \n",
       "\n",
       "   DER_deltaeta_jet_jet  DER_mass_jet_jet  DER_prodeta_jet_jet  \\\n",
       "0                  0.91           124.711                2.666   \n",
       "1               -999.00          -999.000             -999.000   \n",
       "2               -999.00          -999.000             -999.000   \n",
       "3               -999.00          -999.000             -999.000   \n",
       "4               -999.00          -999.000             -999.000   \n",
       "\n",
       "   DER_deltar_tau_lep  DER_pt_tot      ...       PRI_jet_leading_eta  \\\n",
       "0               3.064      41.928      ...                     2.150   \n",
       "1               3.473       2.078      ...                     0.725   \n",
       "2               3.148       9.336      ...                     2.053   \n",
       "3               3.310       0.414      ...                  -999.000   \n",
       "4               3.891      16.405      ...                  -999.000   \n",
       "\n",
       "   PRI_jet_leading_phi  PRI_jet_subleading_pt  PRI_jet_subleading_eta  \\\n",
       "0                0.444                 46.062                    1.24   \n",
       "1                1.158               -999.000                 -999.00   \n",
       "2               -2.028               -999.000                 -999.00   \n",
       "3             -999.000               -999.000                 -999.00   \n",
       "4             -999.000               -999.000                 -999.00   \n",
       "\n",
       "   PRI_jet_subleading_phi  PRI_jet_all_pt    Weight  Label  KaggleSet  \\\n",
       "0                  -2.475         113.497  0.000814      1          t   \n",
       "1                -999.000          46.226  0.681042      0          t   \n",
       "2                -999.000          44.251  0.715742      0          t   \n",
       "3                -999.000          -0.000  1.660654      0          t   \n",
       "4                -999.000           0.000  1.904263      0          t   \n",
       "\n",
       "   KaggleWeight  \n",
       "0      0.002653  \n",
       "1      2.233584  \n",
       "2      2.347389  \n",
       "3      5.446378  \n",
       "4      6.245333  \n",
       "\n",
       "[5 rows x 35 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.iloc[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's create separate arrays\n",
    "eventID = df['EventId']\n",
    "X = df.loc[:,'DER_mass_MMC':'PRI_jet_all_pt']\n",
    "y = df['Label']\n",
    "weight = df['Weight']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#now split into testing and training samples\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_train, X_test, y_train, y_test, eventID_train, event_ID_test, weight_train, weight_test = train_test_split(\n",
    "    X, y, eventID, weight, test_size=0.33, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now let's first look at a NN in sklearn\n",
    "from sklearn.neural_network import MLPClassifier\n",
    "mlp = MLPClassifier(verbose=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 1, loss = 1.18776735\n",
      "Iteration 2, loss = 0.85972503\n",
      "Iteration 3, loss = 0.80940834\n",
      "Iteration 4, loss = 0.76393371\n",
      "Iteration 5, loss = 0.72086151\n",
      "Iteration 6, loss = 0.66503707\n",
      "Iteration 7, loss = 0.65368721\n",
      "Iteration 8, loss = 0.62204937\n",
      "Iteration 9, loss = 0.58840040\n",
      "Iteration 10, loss = 0.56955309\n",
      "Iteration 11, loss = 0.55154781\n",
      "Iteration 12, loss = 0.53934801\n",
      "Iteration 13, loss = 0.52807877\n",
      "Iteration 14, loss = 0.52571351\n",
      "Iteration 15, loss = 0.50227023\n",
      "Iteration 16, loss = 0.49297239\n",
      "Iteration 17, loss = 0.48187416\n",
      "Iteration 18, loss = 0.47043795\n",
      "Iteration 19, loss = 0.47198996\n",
      "Iteration 20, loss = 0.45996074\n",
      "Iteration 21, loss = 0.44496803\n",
      "Iteration 22, loss = 0.44629226\n",
      "Iteration 23, loss = 0.44094936\n",
      "Iteration 24, loss = 0.43497973\n",
      "Iteration 25, loss = 0.42827281\n",
      "Iteration 26, loss = 0.42259840\n",
      "Iteration 27, loss = 0.42137782\n",
      "Iteration 28, loss = 0.41831155\n",
      "Iteration 29, loss = 0.41348078\n",
      "Iteration 30, loss = 0.41127662\n",
      "Iteration 31, loss = 0.40912854\n",
      "Iteration 32, loss = 0.40754783\n",
      "Iteration 33, loss = 0.40573724\n",
      "Iteration 34, loss = 0.40377966\n",
      "Iteration 35, loss = 0.40157661\n",
      "Iteration 36, loss = 0.39962126\n",
      "Iteration 37, loss = 0.39998256\n",
      "Iteration 38, loss = 0.39801556\n",
      "Iteration 39, loss = 0.39762419\n",
      "Iteration 40, loss = 0.39676232\n",
      "Iteration 41, loss = 0.39584833\n",
      "Iteration 42, loss = 0.39580536\n",
      "Iteration 43, loss = 0.39460955\n",
      "Iteration 44, loss = 0.39429449\n",
      "Iteration 45, loss = 0.39432597\n",
      "Iteration 46, loss = 0.39410224\n",
      "Iteration 47, loss = 0.39287974\n",
      "Iteration 48, loss = 0.39348694\n",
      "Iteration 49, loss = 0.39295418\n",
      "Iteration 50, loss = 0.39296233\n",
      "Training loss did not improve more than tol=0.000100 for two consecutive epochs. Stopping.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "MLPClassifier(activation='relu', alpha=0.0001, batch_size='auto', beta_1=0.9,\n",
       "       beta_2=0.999, early_stopping=False, epsilon=1e-08,\n",
       "       hidden_layer_sizes=(100,), learning_rate='constant',\n",
       "       learning_rate_init=0.001, max_iter=200, momentum=0.9,\n",
       "       nesterovs_momentum=True, power_t=0.5, random_state=None,\n",
       "       shuffle=True, solver='adam', tol=0.0001, validation_fraction=0.1,\n",
       "       verbose=True, warm_start=False)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# and train\n",
    "mlp.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8210607401701362"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mlp.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Kaggle competition used the approximate median segnificance (AMS), as defined below, to determine how good a solution was. The number 10, added to the background yield, is a regularization term to decrease the variance of the AMS.\n",
    "\n",
    "Note that if you do not use the full data set (i.e. you split into training and testing) you have to reweigh the inputs so that the subsample yield matches to the toal yield, which we will do below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute approximate median segnificance (AMS)\n",
    "\n",
    "def ams(s,b):\n",
    "    from math import sqrt,log\n",
    "    if b==0:\n",
    "        return 0\n",
    "\n",
    "    return sqrt(2*((s+b+10)*log(1+float(s)/(b+10))-s))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8429549197920467"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's try a different probability cut, not the one given by default to predict().\n",
    "# We choose the top 15%, but can optimize\n",
    "y_train_prob = mlp.predict_proba(X_train)[:, 1]\n",
    "y_test_prob = mlp.predict_proba(X_test)[:, 1]\n",
    "pcut = np.percentile(y_train_prob,85)\n",
    "pcut"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add the probability to the original data frame\n",
    "df['Prob']=mlp.predict_proba(X)[:, 1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f5e8bd3a6d8>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGHJJREFUeJzt3X901PW95/HnOyEkqQm0EM0VUYPHHxSJgowWa1cnxeXQSsGe4lZt5dLTs/FHa13bPa7F41l3LWe3Pfbuskd7Vqgtukcb3du6VQ/FIxdGL7boTSy3Wn5Yq1aolAJqkvEaIPLePxJyQZKZbybznZnPzOtxDsdMvp/5zvudSV5+5jPf73fM3RERkXBUFbsAEREZHQW3iEhgFNwiIoFRcIuIBEbBLSISGAW3iEhgFNwiIoFRcIuIBEbBLSISmHFx7LSpqclbWlpyuu/777/PCSeckN+CSpx6Ln+V1i+o59Hq6ura5+4nRhkbS3C3tLTQ2dmZ031TqRTJZDK/BZU49Vz+Kq1fUM+jZWZ/ijpWSyUiIoFRcIuIBEbBLSISmFjWuEWk9B06dIhdu3bR19cXy/4nTpzItm3bYtl3qYrSc11dHVOnTqWmpibnx1Fwi1SoXbt20djYSEtLC2aW9/339vbS2NiY9/2Wsmw9uzv79+9n165dTJs2LefH0VKJSIXq6+tj8uTJsYS2DM/MmDx58phf5Si4RSqYQrvw8vEzV3CLiARGa9wiAsD6rXvyur9PnfqxrGOqq6tpbW3F3amurubee+/l05/+9Kgfa9myZSxcuJAlS5bkUmqsGhoaSKfTed1nyQV3b1//iL9Al89oLnA1IhKn+vp6tmzZAsDTTz/Nd7/7XZ599tmC1tDf38+4cSUXhRlpqURESkJPTw+f+MQnAEin08ybN48LLriA1tZWfvnLXw6Ne+ihhzjvvPM4//zzue66647bz5133smyZcs4fPgwa9euZfr06XzmM5/hW9/6FgsXLgTgrrvuor29nfnz57N06VL6+vr42te+RmtrK7Nnz2bjxo0ArFmzhm9+85tD+164cCGpVAoYmEnfcccdnH/++cydO5c9ewYmnG+88QYXX3wxF154IXfeeWcsP6uw/jcjImXlgw8+YNasWfT19bF79242bNgADBzr/PjjjzNhwgT27dvH3LlzWbRoEVu3bmXFihU8//zzNDU18c477xyzv9tuu43u7m5++tOfcuDAAa6//nqee+45pk2bxjXXXHPM2K6uLjZt2kR9fT0//OEPAXj55ZfZvn078+fP59VXX81Y+/vvv8/cuXNZsWIFt912G6tXr+aWW27hlltu4cYbb2Tp0qXcd999efxp/SvNuEWkaI4slWzfvp1169axdOlS3B13Z/ny5Zx33nlcfvnl/PnPf2bPnj1s2LCBJUuW0NTUBMCkSZOG9nX33Xfz3nvvcf/992NmbN++nTPOOGPoeOmPBveiRYuor68HYNOmTUOz9+nTp3P66adnDe7x48cPzeDnzJnDm2++CcDzzz8/9FjDvSLIB824RaQkXHzxxezbt4+9e/eydu1a9u7dS1dXFzU1NbS0tNDX14e7j3g43YUXXkhXVxfvvPMOkyZNwt0zPt7Rl18daey4ceM4fPjw0O2jj7+uqakZqqW6upr+/v6hbXEfZqkZt4iUhO3bt/Phhx8yefJkuru7Oemkk6ipqWHjxo386U8DVzydN28ejz32GPv37wc4ZqlkwYIF3H777VxxxRX09vYyffp0Xn/99aGZ8KOPPjriY1966aU8/PDDALz66qu89dZbnHPOObS0tLBlyxYOHz7Mzp07efHFF7P2cckll9DR0QEwtM9804xbRID8H7XV29ubdcyRNW4YmPU++OCDVFdX85WvfIUvfOELJBIJZs2axfTp0wE499xzueOOO7jsssuorq5m9uzZrFmzZmh/V111Fb29vSxatIi1a9fyox/9iAULFtDU1MRFF100Yh033XQTN9xwA62trYwbN441a9ZQW1vLJZdcwrRp02htbWXmzJlccMEFWXtauXIl1157LStXruRLX/pS1vG5sGwvJ3KRSCQ81w9SeHLdeupPax12W7keDqgLzpe/Uux327ZtfPKTn4xt/6VwrZJ0Ok1DQwPuzje+8Q3OOussbr311tgeL2rPw/3szazL3RNRHkdLJSJStlavXs2sWbM499xz6e7u5vrrry92SXmhpRIRKVu33nprrDPsYtGMW0QkMFmD28zOMbMtR/3rMbP/UIjiRETkeFmXStx9BzALwMyqgT8Dj8dcl4iIjGC0SyXzgD+6e+SPkRcRkfwa7ZuTVwM/i6MQESmyHb/K7/6mfCbrkBUrVvDII49QXV1NVVUV999/P6tXr+bb3/42M2bMyGs5cVxetVgiH8dtZuOBt4Fz3f24666aWTvQDtDc3DznyJlDo9Xd00vV+PphtzXWledBMEeONa0kldZzKfY7ceJEzjzzzKHb1X98Jq/7P9jyWaqrq0fc/sILL7B8+XLWrl1LbW0t+/fv5+DBg5x88sl5reOIk08+md27d8ey7yM+/PDDjD0f8dprr9Hd3X3M99ra2iIfxz2aJPwc8NJwoQ3g7quAVTBwAk6uJxtkOgEnqRNwykal9VyK/W7btu3Yk0Xqh58w5erD6uqMJ6P09PTQ3Nw8dMGoI2OTyST33HMPiUSCBx54gO9///tMmTKFs846i9raWu69916WLVvGhAkT6Ozs5C9/+Qs/+MEPWLJkCel0msWLF/Puu+9y6NAhvve977F48eKhx4z7hKCoJ+DU1dUxe/bsnB9nNGvc16BlEhHJk/nz57Nz507OPvtsbrrppuM+QOHtt9/m7rvvZvPmzTzzzDNs3779mO27d+9m06ZNPPXUU9x+++3Av14O9qWXXmLjxo185zvfyXqxqRBFCm4z+xjwb4FfxFuOiFSKhoYGurq6WLVqFSeeeCJf/vKXj7nuyIsvvshll13GpEmTqKmp4aqrrjrm/ldeeSVVVVXMmDFj6EMMRrocbLmJtFTi7v8CTI65FhGpMNXV1SSTSZLJJK2trTz44IND27LNlGtra48b+/DDDw97OdhyozMnRaQoduzYwR/+8Ieh21u2bOH0008fun3RRRfx7LPP8u6779Lf38/Pf/7zrPsc6XKw5aY8D9MQkdE753P53V+Wy7qm02luvvlm3nvvPcaNG8eZZ57JqlWrhj6p/ZRTTmH58uV86lOfYsqUKcyYMYOJEydm3OdIl4MtNwpuESmKOXPm8Otf//q47x/5MF6Aa6+9lvb2dvr7+/niF7/I/PnzAY5ZCweGjs9uamriN7/5zbCPVy7HcIOWSkSkhN11113MmjWLmTNnMm3aNK688spil1QSNOMWkZJ1zz33FLuEkqQZt0gFK8djnEtdPn7mCm6RClVXV8f+/fsV3gXk7uzfv5+6urox7UdLJSIVaurUqezatYu9e/fGsv++vr4xB1RoovRcV1fH1KlTx/Q4Cm6RClVTU8O0adNi238qlRrT9ThCVKietVQiIhIYBbeISGAU3CIigVFwi4gERsEtIhIYBbeISGAU3CIigVFwi4gERsEtIhKYqJ85+XEz+3sz225m28zs4rgLExGR4UU95X0lsM7dl5jZeOBjMdYkIiIZZA1uM5sAXAosA3D3g8DBeMsSEZGRRFkqOQPYC/zUzH5rZj82sxNirktEREZg2a7Fa2YJYDNwibu/YGYrgR53v/Mj49qBdoDm5uY5HR0dORXU3dNL1fj6Ybc11pXnxQzT6TQNDQ3FLqOgKq3nSusX1PNotbW1dbl7IsrYKMH9N8Bmd28ZvP1vgNvd/YqR7pNIJLyzszN6xUd5ct166k9rHXbb5TOac9pnqUulUiSTyWKXUVCV1nOl9QvqebTMLHJwZ10qcfe/ADvN7JzBb80DtuZUmYiIjFnUtYebgYcHjyh5HfhafCWJiEgmkYLb3bcAkabwIiISL505KSISGAW3iEhgFNwiIoFRcIuIBEbBLSISGAW3iEhgFNwiIoFRcIuIBEbBLSISGAW3iEhgFNwiIoFRcIuIBEbBLSISGAW3iEhgFNwiIoFRcIuIBEbBLSISGAW3iEhgFNwiIoGJ9JmTZvYm0At8CPRH/Qh5ERHJv6if8g7Q5u77YqtEREQi0VKJiEhgzN2zDzJ7A3gXcOB+d181zJh2oB2gubl5TkdHR04Fdff0UjW+fthtjXWjeYEQjnQ6TUNDQ7HLKKhK67nS+gX1PFptbW1dUZehowb3FHd/28xOAp4Bbnb350Yan0gkvLOzM3LBR3ty3XrqT2sddtvlM5pz2mepS6VSJJPJYpdRUJXWc6X1C+p5tMwscnBHWipx97cH//tX4HHgopwqExGRMcsa3GZ2gpk1HvkamA+8EndhIiIyvCiLxs3A42Z2ZPwj7r4u1qpERGREWYPb3V8Hzi9ALSIiEoEOBxQRCYyCW0QkMApuEZHAKLhFRAKj4BYRCYyCW0QkMApuEZHABHXVpvVb92TcXq7XMhEROZpm3CIigVFwi4gERsEtIhIYBbeISGAU3CIigVFwi4gERsEtIhIYBbeISGAU3CIigVFwi4gEJnJwm1m1mf3WzJ6KsyAREclsNDPuW4BtcRUiIiLRRApuM5sKXAH8ON5yREQkm6gz7v8J3AYcjrEWERGJwNw98wCzhcDn3f0mM0sC/9HdFw4zrh1oB2hubp7T0dGRU0HdPb1Uja/P6b6NdUFdpXZIOp2moaGh2GUUVKX1XGn9gnoerba2ti53T0QZGyW4/xtwHdAP1AETgF+4+1dHuk8ikfDOzs7oFR/lyXXrqT+tNaf7hno97lQqRTKZLHYZBVVpPVdav6CeR8vMIgd31qUSd/+uu0919xbgamBDptAWEZF46ThuEZHAjGpR2N1TQCqWSkREJBLNuEVEAqPgFhEJjIJbRCQwYR74LCJSLDt+lWFjbuegjJZm3CIigVFwi4gERsEtIhIYBbeISGAU3CIigVFwi4gERsEtIhIYBbeISGAU3CIigVFwi4gEpqxOeV+/dU/G7aF+Qo6IyNE04xYRCYyCW0QkMApuEZHAKLhFRAKTNbjNrM7MXjSzfzaz35vZfylEYSIiMrwoR5UcAD7r7mkzqwE2mdmv3H1zzLWJiMgwsga3uzuQHrxZM/jP4yxKRERGZgO5nGWQWTXQBZwJ3Ofu/2mYMe1AO0Bzc/Ocjo6OnArq7umlanw8H//TWFeah62n02kaGhqKXUZBVVrPldYvlHHPB3pG3JQ+VJVzz21tbV3unogyNlJwDw02+zjwOHCzu78y0rhEIuGdnZ2R93u0J9etp/601pzum02pnoCTSqVIJpPFLqOgKq3nSusXyrjnDJ85mdpdn3PPZhY5uEd1VIm7vwekgAU51CUiInkQ5aiSEwdn2phZPXA5sD3uwkREZHhRFn1PBh4cXOeuAh5z96fiLUtEREYS5aiS3wGzC1CLiIhEoDMnRUQCo+AWEQmMgltEJDAKbhGRwJTmqYQiImOR4SQZzvlc4eqIiWbcIiKBUXCLiARGwS0iEhitcYtIZcm0/g1BrIFXVHCv37pnxG2leuVAEZGP0lKJiEhgKmrGLSJlIttyR5nTjFtEJDAKbhGRwCi4RUQCo+AWEQmMgltEJDAKbhGRwGQ9HNDMTgUeAv4GOAyscveVcRcmIhWuwg/5yyTKcdz9wHfc/SUzawS6zOwZd98ac20iIjKMrEsl7r7b3V8a/LoX2AacEndhIiIyvFGtcZtZCwOf+P5CHMWIiEh25u7RBpo1AM8CK9z9F8NsbwfaAZqbm+d0dHTkVFB3Ty9V4+tzum+cGuviuzpAOp2moaEhtv2XokrrudL6hTz0fKAnf8WMRu2EzNsz1JU+VJVzz21tbV3unogyNlJwm1kN8BTwtLv/XbbxiUTCOzs7ozz+cZ5ct57601pzum+c4rx6YCqVIplMxrb/UlRpPVdav5CHnov15mS2y7pmqCu1uz7nns0scnBnXSoxMwMeALZFCW0REYlXlNf/lwDXAS+b2ZbB7y1397XxlSUiZU+H++Usa3C7+ybAClCLiIhEoOtxi4gcLYBXAgruiPSxZyJSKnStEhGRwGjGLSLxONATxLJDiDTjFhEJjIJbRCQwCm4RkcBojVtEcqc17KLQjFtEJDAKbhGRwCi4RUQCo+AWEQmM3pzMg0ynw4NOiZeA6c3HkqQZt4hIYBTcIiKB0VJJAWRbStGTIEWjpZAgacYtIhIYBbeISGD0Kr0E9Pb1j7icoiNSZMy0HFJ2onzK+0/M7K9m9kohChIRkcyiLJWsARbEXIeIiEQU5VPenzOzlvhLkeHo5B7JSkshFcfcPfuggeB+yt1nZhjTDrQDNDc3z+no6MipoO6eXqrG1+d031AdPvhBzj031oX5NkU6naahoaHYZRTMmPs90JO/YgokfcBpqLVil1FQ6UNVOT/PbW1tXe6eiDI2b3/17r4KWAWQSCQ8mUzmtJ8n162n/rTWfJUVhA/eejnnnvszbCvl2XgqlSLX35EQjbnfAGfVqTcPkWypKXYZBZXaXV+Q32sdDigiEpgwX2dLJFofD0iAM2opnqzBbWY/A5JAk5ntAv6zuz8Qd2EiZeVAj8JZ8ibKUSXXFKIQKTzNyPNMwSwFoqUSGVGmYFeoixSPgltyUpGzdc2opUQouKWyKHylDCi4JRa6BrlIfPT3I0UxpisiatYsFU7BLSUn6/p5dYEKESlRCm6JRdPbGzJu332oNuuYkWzJsG3WqR/PaZ8iIVFwS05yDd24bdn53pjur+CXECi4K1iphm8xZQp+hbqUCgV3RJlCbt+UzxawkmMpfAtnbLP5E/JWh4iCuwSMO9SjAC5zHxzsZ8vO90fcrtm8jIaCOw/GGrq7qc1TJRKqsczmFfqVp6KCW7NaKUcK/cpTUcEtIsfKFvoK9tJUVsGtGbVIfo1lNv/BwVqt68ckqOBWMIuUD832c1dywa0jLEQEtHafSckFt4jIWJX7iVSRgtvMFgArgWrgx+7+32OtSkQkJuWwRFOVbYCZVQP3AZ8DZgDXmNmMuAsTEZHhRZlxXwS85u6vA5hZB7AY2BpnYSIixTCmSxuMq89fIRlknXEDpwA7j7q9a/B7IiJSBFFm3DbM9/y4QWbtQPvgzbSZ7cixpiZgX473DZV6Ln+V1i+o59E6PerAKMG9Czj1qNtTgbc/OsjdVwGroj7wSMys090TY91PSNRz+au0fkE9xynKUsk/AWeZ2TQzGw9cDTwRb1kiIjKSrDNud+83s28CTzNwOOBP3P33sVcmIiLDinQct7uvBdbGXMsRY15uCZB6Ln+V1i+o59iY+3HvM4qISAmLssYtIiIlpGjBbWYLzGyHmb1mZrcPs73WzB4d3P6CmbUUvsr8idDvt81sq5n9zsz+wcwiHxpUqrL1fNS4JWbmZhb8EQhRejazfzf4XP/ezB4pdI35FuF3+zQz22hmvx38/f58MerMFzP7iZn91cxeGWG7mdn/Gvx5/M7MLsh7Ee5e8H8MvMn5R+AMYDzwz8CMj4y5Cfjfg19fDTxajFoL2G8b8LHBr28Mud+oPQ+OawSeAzYDiWLXXYDn+Szgt8AnBm+fVOy6C9DzKuDGwa9nAG8Wu+4x9nwpcAHwygjbPw/8ioFzYOYCL+S7hmLNuIdOo3f3g8CR0+iPthh4cPDrvwfmmdlwJwOFIGu/7r7R3f9l8OZmBo6XD1mU5xjgbuAHQF8hi4tJlJ7/PXCfu78L4O5/LXCN+RalZwcmDH49kWHOAwmJuz8HvJNhyGLgIR+wGfi4mZ2czxqKFdxRTqMfGuPu/UA3MLkg1eXfaC8b8HUG/o8dsqw9m9ls4FR3f6qQhcUoyvN8NnC2mT1vZpsHr7wZsig93wV81cx2MXB02s2FKa1oYr9MSLGuxx3lNPpIp9oHInIvZvZVIAFcFmtF8cvYs5lVAf8DWFaoggogyvM8joHlkiQDr6r+0cxmuvsYrmxUVFF6vgZY4+4/NLOLgf8z2PPh+Msritizq1gz7iin0Q+NMbNxDLzEyvTypJRFumyAmV0O3AEscvcDBaotLtl6bgRmAikze5OBtcAnAn+DMurv9S/d/ZC7vwHsYCDIQxWl568DjwG4+2+AOgau6VGuIv29j0WxgjvKafRPAH87+PUSYIMPrvwHKGu/g8sG9zMQ2qGve0KWnt29292b3L3F3VsYWNdf5O6dxSk3L6L8Xv8/Bt6IxsyaGFg6eb2gVeZXlJ7fAuYBmNknGQjuvQWtsrCeAJYOHl0yF+h29915fYQivjP7eeBVBt6RvmPwe/+VgT9eGHhy/y/wGvAicEax302Oud/1wB5gy+C/J4pdc9w9f2RsisCPKon4PBvwdwxcz/5l4Opi11yAnmcAzzNwxMkWYH6xax5jvz8DdgOHGJhdfx24AbjhqOf4vsGfx8tx/F7rzEkRkcDozEkRkcAouEVEAqPgFhEJjIJbRCQwCm4RkcAouEVEAqPgFhEJjIJbRCQw/x9jw/xs+hcnigAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5e7b9c8be0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "kwargs = dict(histtype='stepfilled', alpha=0.3, normed=True, bins=40)\n",
    "\n",
    "df[df.Label==0].Prob.hist(label='Background',**kwargs)\n",
    "df[df.Label==1].Prob.hist(label='Signal',**kwargs)\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's calculate the total weights (yields)\n",
    "sigall = weight.dot(y)\n",
    "backall = weight.dot(y == 0)\n",
    "\n",
    "# The training weights\n",
    "sigtrain = weight_train.dot(y_train)\n",
    "backtrain = weight_train.dot(y_train == 0)\n",
    "\n",
    "# The training weights\n",
    "sigtest = weight_test.dot(y_test)\n",
    "backtest = weight_test.dot(y_test == 0)\n",
    "\n",
    "# aside:  these can also be done by looping instead of using a dot product\n",
    "#  (Usually vectorized operations are faster for interpreted code)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's look at event yields that pass our selection\n",
    "sigtrain_sel = weight_train.dot(np.multiply(y_train, y_train_prob > pcut))\n",
    "backtrain_sel = weight_train.dot(np.multiply(y_train == 0, y_train_prob > pcut))\n",
    "\n",
    "sigtest_sel = weight_test.dot(np.multiply(y_test, y_test_prob > pcut))\n",
    "backtest_sel = weight_test.dot(np.multiply(y_test == 0, y_test_prob > pcut))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Corrected selected yields in training sample, signal = 211.87694541459166 , background = 5151.928791943724\n",
      "Corrected selected yields in test sample, signal = 210.78233635205316 , background = 5296.97764089601\n"
     ]
    }
   ],
   "source": [
    "# Now we need to correct the selected yields to be is if we used the full sample\n",
    "sigtrain_sel_corr = sigtrain_sel*sigall/sigtrain\n",
    "backtrain_sel_corr = backtrain_sel*backall/backtrain\n",
    "\n",
    "sigtest_sel_corr = sigtest_sel*sigall/sigtest\n",
    "backtest_sel_corr = backtest_sel*backall/backtest\n",
    "\n",
    "print(\"Corrected selected yields in training sample, signal =\", sigtrain_sel_corr, \", background =\",backtrain_sel_corr)\n",
    "print(\"Corrected selected yields in test sample, signal =\", sigtest_sel_corr, \", background =\",backtest_sel_corr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AMS of training sample 2.9291829786137367\n",
      "AMS of test sample 2.8745708039703985\n"
     ]
    }
   ],
   "source": [
    "print(\"AMS of training sample\", ams(sigtrain_sel_corr,backtrain_sel_corr))\n",
    "print(\"AMS of test sample\", ams(sigtest_sel_corr,backtest_sel_corr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How did we do? Worse than the BDT from yesterday.\n",
    "![Comparison with submissions](data/tr150908_davidRousseau_TMVAFuture_HiggsML.001.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Neural networks are quite sensitive to feature scaling, so let's try to scale the features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import RobustScaler\n",
    "scaler = RobustScaler()\n",
    "scaler.fit(X_train)\n",
    "X_train_scaled = scaler.transform(X_train)\n",
    "X_test_scaled = scaler.transform(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 1, loss = 0.40396788\n",
      "Iteration 2, loss = 0.37804418\n",
      "Iteration 3, loss = 0.37313823\n",
      "Iteration 4, loss = 0.37037245\n",
      "Iteration 5, loss = 0.36844364\n",
      "Iteration 6, loss = 0.36680546\n",
      "Iteration 7, loss = 0.36583018\n",
      "Iteration 8, loss = 0.36499575\n",
      "Iteration 9, loss = 0.36428769\n",
      "Iteration 10, loss = 0.36379259\n",
      "Iteration 11, loss = 0.36317937\n",
      "Iteration 12, loss = 0.36269424\n",
      "Iteration 13, loss = 0.36247411\n",
      "Iteration 14, loss = 0.36214801\n",
      "Iteration 15, loss = 0.36198055\n",
      "Iteration 16, loss = 0.36165097\n",
      "Iteration 17, loss = 0.36133456\n",
      "Iteration 18, loss = 0.36119266\n",
      "Iteration 19, loss = 0.36102897\n",
      "Iteration 20, loss = 0.36084392\n",
      "Iteration 21, loss = 0.36058338\n",
      "Iteration 22, loss = 0.36060154\n",
      "Iteration 23, loss = 0.36033204\n",
      "Iteration 24, loss = 0.36015293\n",
      "Iteration 25, loss = 0.36000207\n",
      "Iteration 26, loss = 0.35972131\n",
      "Iteration 27, loss = 0.35971010\n",
      "Iteration 28, loss = 0.35964028\n",
      "Iteration 29, loss = 0.35953327\n",
      "Iteration 30, loss = 0.35931813\n",
      "Iteration 31, loss = 0.35916346\n",
      "Iteration 32, loss = 0.35931728\n",
      "Iteration 33, loss = 0.35915913\n",
      "Iteration 34, loss = 0.35891894\n",
      "Iteration 35, loss = 0.35885837\n",
      "Iteration 36, loss = 0.35884762\n",
      "Iteration 37, loss = 0.35867448\n",
      "Iteration 38, loss = 0.35869264\n",
      "Iteration 39, loss = 0.35860610\n",
      "Iteration 40, loss = 0.35843331\n",
      "Iteration 41, loss = 0.35850879\n",
      "Iteration 42, loss = 0.35846496\n",
      "Iteration 43, loss = 0.35822736\n",
      "Iteration 44, loss = 0.35830133\n",
      "Iteration 45, loss = 0.35814101\n",
      "Iteration 46, loss = 0.35814574\n",
      "Training loss did not improve more than tol=0.000100 for two consecutive epochs. Stopping.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "MLPClassifier(activation='relu', alpha=0.0001, batch_size='auto', beta_1=0.9,\n",
       "       beta_2=0.999, early_stopping=False, epsilon=1e-08,\n",
       "       hidden_layer_sizes=(100,), learning_rate='constant',\n",
       "       learning_rate_init=0.001, max_iter=200, momentum=0.9,\n",
       "       nesterovs_momentum=True, power_t=0.5, random_state=None,\n",
       "       shuffle=True, solver='adam', tol=0.0001, validation_fraction=0.1,\n",
       "       verbose=True, warm_start=False)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# and train a new network\n",
    "mlp_scaled = MLPClassifier(verbose=True)\n",
    "mlp_scaled.fit(X_train_scaled, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8371633107299856"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mlp_scaled.score(X_test_scaled, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's try a different probability cut, not the one given by default to predict().\n",
    "# We choose the top 15%, but can optimize\n",
    "y_train_prob_scaled = mlp_scaled.predict_proba(X_train_scaled)[:, 1]\n",
    "y_test_prob_scaled = mlp_scaled.predict_proba(X_test_scaled)[:, 1]\n",
    "pcut_scaled = np.percentile(y_train_prob_scaled,85)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's look at event yields that pass our selection\n",
    "sigtrain_sel_scaled = weight_train.dot(np.multiply(y_train, y_train_prob_scaled > pcut_scaled))\n",
    "backtrain_sel_scaled = weight_train.dot(np.multiply(y_train == 0, y_train_prob_scaled > pcut_scaled))\n",
    "\n",
    "sigtest_sel_scaled = weight_test.dot(np.multiply(y_test, y_test_prob_scaled > pcut_scaled))\n",
    "backtest_sel_scaled = weight_test.dot(np.multiply(y_test == 0, y_test_prob_scaled > pcut_scaled))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Corrected selected yields in training sample, signal = 220.35738592020525 , background = 3836.136115385027\n",
      "Corrected selected yields in test sample, signal = 219.43945088666317 , background = 4125.75556096465\n"
     ]
    }
   ],
   "source": [
    "# Now we need to correct the selected yields to be is if we used the full sample\n",
    "sigtrain_sel_scaled_corr = sigtrain_sel_scaled*sigall/sigtrain\n",
    "backtrain_sel_scaled_corr = backtrain_sel_scaled*backall/backtrain\n",
    "\n",
    "sigtest_sel_scaled_corr = sigtest_sel_scaled*sigall/sigtest\n",
    "backtest_sel_scaled_corr = backtest_sel_scaled*backall/backtest\n",
    "\n",
    "print(\"Corrected selected yields in training sample, signal =\", sigtrain_sel_scaled_corr, \", background =\",backtrain_sel_scaled_corr)\n",
    "print(\"Corrected selected yields in test sample, signal =\", sigtest_sel_scaled_corr, \", background =\",backtest_sel_scaled_corr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AMS of training sample 3.5200207678782083\n",
      "AMS of test sample 3.382695240698033\n"
     ]
    }
   ],
   "source": [
    "print(\"AMS of training sample\", ams(sigtrain_sel_scaled_corr,backtrain_sel_scaled_corr))\n",
    "print(\"AMS of test sample\", ams(sigtest_sel_scaled_corr,backtest_sel_scaled_corr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We improved somewhat.\n",
    "\n",
    "SciKit Learn has simple NNs, but if you want to do deep NNs, or train on GPUs, you probalby want to use something like Keras instead. Let's try to create a simple NN using Keras."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/software/opt/xenial/x86_64/anaconda3/5.1/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n",
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(1337)  # for reproducibility\n",
    "\n",
    "from keras.models import Sequential\n",
    "from keras.layers import Dense\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the model\n",
    "from keras import regularizers\n",
    "\n",
    "model = Sequential()\n",
    "model.add(Dense(100, activation='relu', input_shape=(30,), kernel_regularizer=regularizers.l2(0.0001)))\n",
    "model.add(Dense(1, activation='sigmoid'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compile model\n",
    "model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: 1.990918653699893, 1: 1182.3150708556338}"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "class_weight = {0: y_train.shape[0]/backtrain, 1:y_train.shape[0]/sigtrain}\n",
    "class_weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1/20\n",
      " 10624/548219 [..............................] - ETA: 8s - loss: nan - acc: 0.0000e+00"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/software/opt/xenial/x86_64/anaconda3/5.1/lib/python3.6/site-packages/pandas/core/series.py:696: FutureWarning: \n",
      "Passing list-likes to .loc or [] with any missing label will raise\n",
      "KeyError in the future, you can use .reindex() as an alternative.\n",
      "\n",
      "See the documentation here:\n",
      "http://pandas.pydata.org/pandas-docs/stable/indexing.html#deprecate-loc-reindex-listlike\n",
      "  return self.loc[key]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "548219/548219 [==============================] - 8s 14us/step - loss: nan - acc: 0.0000e+00\n",
      "Epoch 2/20\n",
      "548219/548219 [==============================] - 8s 15us/step - loss: nan - acc: 0.0000e+00\n",
      "Epoch 3/20\n",
      "313216/548219 [================>.............] - ETA: 3s - loss: nan - acc: 0.0000e+00"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-32-d56e02bbe37d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_train_scaled\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my_train\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepochs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m20\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbatch_size\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m128\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msample_weight\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mweight_train\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/software/opt/xenial/x86_64/anaconda3/5.1/lib/python3.6/site-packages/keras/models.py\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, x, y, batch_size, epochs, verbose, callbacks, validation_split, validation_data, shuffle, class_weight, sample_weight, initial_epoch, steps_per_epoch, validation_steps, **kwargs)\u001b[0m\n\u001b[1;32m    963\u001b[0m                               \u001b[0minitial_epoch\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minitial_epoch\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    964\u001b[0m                               \u001b[0msteps_per_epoch\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msteps_per_epoch\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 965\u001b[0;31m                               validation_steps=validation_steps)\n\u001b[0m\u001b[1;32m    966\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    967\u001b[0m     def evaluate(self, x=None, y=None,\n",
      "\u001b[0;32m/software/opt/xenial/x86_64/anaconda3/5.1/lib/python3.6/site-packages/keras/engine/training.py\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, x, y, batch_size, epochs, verbose, callbacks, validation_split, validation_data, shuffle, class_weight, sample_weight, initial_epoch, steps_per_epoch, validation_steps, **kwargs)\u001b[0m\n\u001b[1;32m   1667\u001b[0m                               \u001b[0minitial_epoch\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minitial_epoch\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1668\u001b[0m                               \u001b[0msteps_per_epoch\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msteps_per_epoch\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1669\u001b[0;31m                               validation_steps=validation_steps)\n\u001b[0m\u001b[1;32m   1670\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1671\u001b[0m     def evaluate(self, x=None, y=None,\n",
      "\u001b[0;32m/software/opt/xenial/x86_64/anaconda3/5.1/lib/python3.6/site-packages/keras/engine/training.py\u001b[0m in \u001b[0;36m_fit_loop\u001b[0;34m(self, f, ins, out_labels, batch_size, epochs, verbose, callbacks, val_f, val_ins, shuffle, callback_metrics, initial_epoch, steps_per_epoch, validation_steps)\u001b[0m\n\u001b[1;32m   1204\u001b[0m                         \u001b[0mins_batch\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mins_batch\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtoarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1205\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1206\u001b[0;31m                     \u001b[0mouts\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mins_batch\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1207\u001b[0m                     \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mouts\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1208\u001b[0m                         \u001b[0mouts\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mouts\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/software/opt/xenial/x86_64/anaconda3/5.1/lib/python3.6/site-packages/keras/backend/tensorflow_backend.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, inputs)\u001b[0m\n\u001b[1;32m   2473\u001b[0m         \u001b[0msession\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_session\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2474\u001b[0m         updated = session.run(fetches=fetches, feed_dict=feed_dict,\n\u001b[0;32m-> 2475\u001b[0;31m                               **self.session_kwargs)\n\u001b[0m\u001b[1;32m   2476\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mupdated\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moutputs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2477\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/software/opt/xenial/x86_64/anaconda3/5.1/lib/python3.6/site-packages/tensorflow/python/client/session.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, fetches, feed_dict, options, run_metadata)\u001b[0m\n\u001b[1;32m    887\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    888\u001b[0m       result = self._run(None, fetches, feed_dict, options_ptr,\n\u001b[0;32m--> 889\u001b[0;31m                          run_metadata_ptr)\n\u001b[0m\u001b[1;32m    890\u001b[0m       \u001b[0;32mif\u001b[0m \u001b[0mrun_metadata\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    891\u001b[0m         \u001b[0mproto_data\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtf_session\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mTF_GetBuffer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrun_metadata_ptr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/software/opt/xenial/x86_64/anaconda3/5.1/lib/python3.6/site-packages/tensorflow/python/client/session.py\u001b[0m in \u001b[0;36m_run\u001b[0;34m(self, handle, fetches, feed_dict, options, run_metadata)\u001b[0m\n\u001b[1;32m   1118\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mfinal_fetches\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mfinal_targets\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mhandle\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mfeed_dict_tensor\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1119\u001b[0m       results = self._do_run(handle, final_targets, final_fetches,\n\u001b[0;32m-> 1120\u001b[0;31m                              feed_dict_tensor, options, run_metadata)\n\u001b[0m\u001b[1;32m   1121\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1122\u001b[0m       \u001b[0mresults\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/software/opt/xenial/x86_64/anaconda3/5.1/lib/python3.6/site-packages/tensorflow/python/client/session.py\u001b[0m in \u001b[0;36m_do_run\u001b[0;34m(self, handle, target_list, fetch_list, feed_dict, options, run_metadata)\u001b[0m\n\u001b[1;32m   1315\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mhandle\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1316\u001b[0m       return self._do_call(_run_fn, self._session, feeds, fetches, targets,\n\u001b[0;32m-> 1317\u001b[0;31m                            options, run_metadata)\n\u001b[0m\u001b[1;32m   1318\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1319\u001b[0m       \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_do_call\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_prun_fn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_session\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhandle\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfeeds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfetches\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/software/opt/xenial/x86_64/anaconda3/5.1/lib/python3.6/site-packages/tensorflow/python/client/session.py\u001b[0m in \u001b[0;36m_do_call\u001b[0;34m(self, fn, *args)\u001b[0m\n\u001b[1;32m   1321\u001b[0m   \u001b[0;32mdef\u001b[0m \u001b[0m_do_call\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1322\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1323\u001b[0;31m       \u001b[0;32mreturn\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1324\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0merrors\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mOpError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1325\u001b[0m       \u001b[0mmessage\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcompat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_text\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0me\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmessage\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/software/opt/xenial/x86_64/anaconda3/5.1/lib/python3.6/site-packages/tensorflow/python/client/session.py\u001b[0m in \u001b[0;36m_run_fn\u001b[0;34m(session, feed_dict, fetch_list, target_list, options, run_metadata)\u001b[0m\n\u001b[1;32m   1300\u001b[0m           return tf_session.TF_Run(session, options,\n\u001b[1;32m   1301\u001b[0m                                    \u001b[0mfeed_dict\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfetch_list\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtarget_list\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1302\u001b[0;31m                                    status, run_metadata)\n\u001b[0m\u001b[1;32m   1303\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1304\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_prun_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msession\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhandle\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfeed_dict\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfetch_list\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "model.fit(X_train_scaled, y_train, epochs=20, batch_size=128, sample_weight=weight_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's try a different probability cut, not the one given by default to predict().\n",
    "# We choose the top 15%, but can optimize\n",
    "y_train_prob_keras = model.predict(X_train_scaled)[:, 0]\n",
    "y_test_prob_keras = model.predict(X_test_scaled)[:, 0]\n",
    "pcut_keras = np.percentile(y_train_prob_keras,85)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's look at event yields that pass our selection\n",
    "sigtrain_sel_keras = weight_train.dot(np.multiply(y_train, y_train_prob_keras > pcut_keras))\n",
    "backtrain_sel_keras = weight_train.dot(np.multiply(y_train == 0, y_train_prob_keras > pcut_keras))\n",
    "\n",
    "sigtest_sel_keras = weight_test.dot(np.multiply(y_test, y_test_prob_keras > pcut_keras))\n",
    "backtest_sel_keras = weight_test.dot(np.multiply(y_test == 0, y_test_prob_keras > pcut_keras))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Corrected selected yields in training sample, signal = 228.45506778072468 , background = 4170.447242035395\n",
      "Corrected selected yields in test sample, signal = 227.24568728934773 , background = 4359.426704211413\n"
     ]
    }
   ],
   "source": [
    "# Now we need to correct the selected yields to be is if we used the full sample\n",
    "sigtrain_sel_keras_corr = sigtrain_sel_keras*sigall/sigtrain\n",
    "backtrain_sel_keras_corr = backtrain_sel_keras*backall/backtrain\n",
    "\n",
    "sigtest_sel_keras_corr = sigtest_sel_keras*sigall/sigtest\n",
    "backtest_sel_keras_corr = backtest_sel_keras*backall/backtest\n",
    "\n",
    "print(\"Corrected selected yields in training sample, signal =\", sigtrain_sel_keras_corr, \", background =\",backtrain_sel_keras_corr)\n",
    "print(\"Corrected selected yields in test sample, signal =\", sigtest_sel_keras_corr, \", background =\",backtest_sel_keras_corr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AMS of training sample 3.501901630959638\n",
      "AMS of test sample 3.4086509153663713\n"
     ]
    }
   ],
   "source": [
    "print(\"AMS of training sample\", ams(sigtrain_sel_keras_corr,backtrain_sel_keras_corr))\n",
    "print(\"AMS of test sample\", ams(sigtest_sel_keras_corr,backtest_sel_keras_corr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We only made a single layer NN in Keras. However, you can easily change the structure of the network. As an assignment, try adding an extra hidden layer and changing the number of neurons."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a few things we can easily vary: number of hidden layers, the activation function, the regularization ($\\alpha$). Let's go back to MLPClassifer (scaled) and play with some of them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 1, loss = 0.38988645\n",
      "Iteration 2, loss = 0.37209813\n",
      "Iteration 3, loss = 0.36866211\n",
      "Iteration 4, loss = 0.36662743\n",
      "Iteration 5, loss = 0.36519167\n",
      "Iteration 6, loss = 0.36420247\n",
      "Iteration 7, loss = 0.36365862\n",
      "Iteration 8, loss = 0.36301438\n",
      "Iteration 9, loss = 0.36250376\n",
      "Iteration 10, loss = 0.36234049\n",
      "Iteration 11, loss = 0.36197566\n",
      "Iteration 12, loss = 0.36163091\n",
      "Iteration 13, loss = 0.36146661\n",
      "Iteration 14, loss = 0.36108755\n",
      "Iteration 15, loss = 0.36105851\n",
      "Iteration 16, loss = 0.36086552\n",
      "Iteration 17, loss = 0.36060469\n",
      "Iteration 18, loss = 0.36051176\n",
      "Iteration 19, loss = 0.36040121\n",
      "Iteration 20, loss = 0.36036926\n",
      "Iteration 21, loss = 0.36006069\n",
      "Iteration 22, loss = 0.36024082\n",
      "Iteration 23, loss = 0.36014752\n",
      "Iteration 24, loss = 0.35992278\n",
      "Iteration 25, loss = 0.35967233\n",
      "Iteration 26, loss = 0.35987269\n",
      "Iteration 27, loss = 0.35960379\n",
      "Iteration 28, loss = 0.35965437\n",
      "Training loss did not improve more than tol=0.000100 for two consecutive epochs. Stopping.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "MLPClassifier(activation='relu', alpha=0.01, batch_size='auto', beta_1=0.9,\n",
       "       beta_2=0.999, early_stopping=False, epsilon=1e-08,\n",
       "       hidden_layer_sizes=(100, 100), learning_rate='constant',\n",
       "       learning_rate_init=0.001, max_iter=200, momentum=0.9,\n",
       "       nesterovs_momentum=True, power_t=0.5, random_state=None,\n",
       "       shuffle=True, solver='adam', tol=0.0001, validation_fraction=0.1,\n",
       "       verbose=True, warm_start=False)"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mlp_play = MLPClassifier(activation='relu', hidden_layer_sizes=(100,100), alpha=0.01, verbose=True)\n",
    "mlp_play.fit(X_train_scaled, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.83794844066528651"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mlp_play.score(X_test_scaled, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's try a different probability cut, not the one given by default to predict().\n",
    "# We choose the top 15%, but can optimize\n",
    "y_train_prob_play = mlp_play.predict_proba(X_train_scaled)[:, 1]\n",
    "y_test_prob_play = mlp_play.predict_proba(X_test_scaled)[:, 1]\n",
    "pcut_play = np.percentile(y_train_prob_scaled,85)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's look at event yields that pass our selection\n",
    "sigtrain_sel_play = weight_train.dot(np.multiply(y_train, y_train_prob_play > pcut_play))\n",
    "backtrain_sel_play = weight_train.dot(np.multiply(y_train == 0, y_train_prob_play > pcut_play))\n",
    "\n",
    "sigtest_sel_play = weight_test.dot(np.multiply(y_test, y_test_prob_play > pcut_play))\n",
    "backtest_sel_play = weight_test.dot(np.multiply(y_test == 0, y_test_prob_play > pcut_play))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Corrected selected yields in training sample, signal = 285.721029386 , background = 6262.75114743\n",
      "Corrected selected yields in test sample, signal = 283.997485301 , background = 6650.1640173\n"
     ]
    }
   ],
   "source": [
    "# Now we need to correct the selected yields to be is if we used the full sample\n",
    "sigtrain_sel_play_corr = sigtrain_sel_play*sigall/sigtrain\n",
    "backtrain_sel_play_corr = backtrain_sel_play*backall/backtrain\n",
    "\n",
    "sigtest_sel_play_corr = sigtest_sel_play*sigall/sigtest\n",
    "backtest_sel_play_corr = backtest_sel_play*backall/backtest\n",
    "\n",
    "print(\"Corrected selected yields in training sample, signal =\", sigtrain_sel_play_corr, \", background =\",backtrain_sel_play_corr)\n",
    "print(\"Corrected selected yields in test sample, signal =\", sigtest_sel_play_corr, \", background =\",backtest_sel_play_corr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AMS of training sample 3.580676739279546\n",
      "AMS of test sample 3.45564003442452\n"
     ]
    }
   ],
   "source": [
    "print(\"AMS of training sample\", ams(sigtrain_sel_play_corr,backtrain_sel_play_corr))\n",
    "print(\"AMS of test sample\", ams(sigtest_sel_play_corr,backtest_sel_play_corr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Problems (can do with either MLPClassifier or Keras):\n",
    "1. Vary the structure of the tree (number of hidden layers, number of neurons)\n",
    "1. Vary the activation. (In Keras can do it per layer, in MLPClassifier only for all)\n",
    "1. Vary the regularization. May have to do this as the structure changes.\n",
    "1. Try using derivied variables only or primary variables only.\n",
    "1. Missing data is represented by -999 before scaling. Is there a better value to use in the training?\n",
    "1. Try using the event weights to better match the background and signal shapes in the training. Note, though, that you should still treat background and signal separately; don't scale the signal down by the weight."
   ]
  },
  {
   "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.4"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}