{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "2ef4894b-b542-48c9-82aa-2f678750fbaf"
      },
      "source": [
        "###Airbnb New User Bookings Competition\n",
        "*Author*: **Sandro Vega Pons** (sv.pons@gmail.com)\n",
        "\n",
        "The main 3 points of this notebook are:\n",
        "\n",
        " 1. Source code of the ensemble techniques I used in my solution.\n",
        " 2. Example of how to use them in a 3-layer learning architecture.\n",
        " 3. Analysis of the performance of the methods on problems with different number of classes. \n",
        "    Comparison with stack generalization based on LogisticRegression (sklearn implementation) and GradientBoosting (XGBoost implementation)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "cf8c4ab5-6499-48d8-a911-f62dd8d7f37b"
      },
      "source": [
        "# 1- Ensemble techniques based on scipy.optimize package. \n",
        "\n",
        "## Source code of the two ensemble techniques I used in my solution (EN_optA and EN_optB). \n",
        "Given a set of predictions (e.g. predictions obtained with different or the same classifier with different parameters values),\n",
        "the two ensemblers define two different linear problems and find the optimal coefficients that minimize an objective function \n",
        "(in this case the multi-class logloss).\n",
        "\n",
        "Useful ideas were taken from this [discussion](https://www.kaggle.com/c/otto-group-product-classification-challenge/forums/t/13868/ensamble-weights) \n",
        "on the *Otto Group Product Classification Challenge* forum. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "6acfa7ab-a9f3-4f30-8849-882af0145bfa"
      },
      "outputs": [],
      "source": [
        "import numpy as  np\n",
        "from sklearn.metrics import log_loss\n",
        "from sklearn.base import BaseEstimator\n",
        "from scipy.optimize import minimize"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "c2ae75ea-66a6-4a82-9b8c-589402f59d04"
      },
      "source": [
        "## First ensemble technique (EN_optA)\n",
        "Given a set of predictions $X_1, X_2, ..., X_n$,  it computes the optimal set of weights\n",
        "$w_1, w_2, ..., w_n$; such that minimizes $log\\_loss(y_T, y_E)$, \n",
        "where $y_E = X_1*w_1 + X_2*w_2 +...+ X_n*w_n$ and $y_T$ is the true solution."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "512e8a1f-0a1c-4096-8df4-b7c31cff65e2"
      },
      "outputs": [],
      "source": [
        "def objf_ens_optA(w, Xs, y, n_class=12):\n",
        "    \"\"\"\n",
        "    Function to be minimized in the EN_optA ensembler.\n",
        "    \n",
        "    Parameters:\n",
        "    ----------\n",
        "    w: array-like, shape=(n_preds)\n",
        "       Candidate solution to the optimization problem (vector of weights).\n",
        "    Xs: list of predictions to combine\n",
        "       Each prediction is the solution of an individual classifier and has a\n",
        "       shape=(n_samples, n_classes).\n",
        "    y: array-like sahpe=(n_samples,)\n",
        "       Class labels\n",
        "    n_class: int\n",
        "       Number of classes in the problem (12 in Airbnb competition)\n",
        "    \n",
        "    Return:\n",
        "    ------\n",
        "    score: Score of the candidate solution.\n",
        "    \"\"\"\n",
        "    w = np.abs(w)\n",
        "    sol = np.zeros(Xs[0].shape)\n",
        "    for i in range(len(w)):\n",
        "        sol += Xs[i] * w[i]\n",
        "    #Using log-loss as objective function (different objective functions can be used here). \n",
        "    score = log_loss(y, sol)   \n",
        "    return score\n",
        "        \n",
        "\n",
        "class EN_optA(BaseEstimator):\n",
        "    \"\"\"\n",
        "    Given a set of predictions $X_1, X_2, ..., X_n$,  it computes the optimal set of weights\n",
        "    $w_1, w_2, ..., w_n$; such that minimizes $log\\_loss(y_T, y_E)$, \n",
        "    where $y_E = X_1*w_1 + X_2*w_2 +...+ X_n*w_n$ and $y_T$ is the true solution.\n",
        "    \"\"\"\n",
        "    def __init__(self, n_class=12):\n",
        "        super(EN_optA, self).__init__()\n",
        "        self.n_class = n_class\n",
        "        \n",
        "    def fit(self, X, y):\n",
        "        \"\"\"\n",
        "        Learn the optimal weights by solving an optimization problem.\n",
        "        \n",
        "        Parameters:\n",
        "        ----------\n",
        "        Xs: list of predictions to be ensembled\n",
        "           Each prediction is the solution of an individual classifier and has \n",
        "           shape=(n_samples, n_classes).\n",
        "        y: array-like\n",
        "           Class labels\n",
        "        \"\"\"\n",
        "        Xs = np.hsplit(X, X.shape[1]/self.n_class)\n",
        "        #Initial solution has equal weight for all individual predictions.\n",
        "        x0 = np.ones(len(Xs)) / float(len(Xs)) \n",
        "        #Weights must be bounded in [0, 1]\n",
        "        bounds = [(0,1)]*len(x0)   \n",
        "        #All weights must sum to 1\n",
        "        cons = ({'type':'eq','fun':lambda w: 1-sum(w)})\n",
        "        #Calling the solver\n",
        "        res = minimize(objf_ens_optA, x0, args=(Xs, y, self.n_class), \n",
        "                       method='SLSQP', \n",
        "                       bounds=bounds,\n",
        "                       constraints=cons\n",
        "                       )\n",
        "        self.w = res.x\n",
        "        return self\n",
        "    \n",
        "    def predict_proba(self, X):\n",
        "        \"\"\"\n",
        "        Use the weights learned in training to predict class probabilities.\n",
        "        \n",
        "        Parameters:\n",
        "        ----------\n",
        "        Xs: list of predictions to be blended.\n",
        "            Each prediction is the solution of an individual classifier and has \n",
        "            shape=(n_samples, n_classes).\n",
        "            \n",
        "        Return:\n",
        "        ------\n",
        "        y_pred: array_like, shape=(n_samples, n_class)\n",
        "                The blended prediction.\n",
        "        \"\"\"\n",
        "        Xs = np.hsplit(X, X.shape[1]/self.n_class)\n",
        "        y_pred = np.zeros(Xs[0].shape)\n",
        "        for i in range(len(self.w)):\n",
        "            y_pred += Xs[i] * self.w[i] \n",
        "        return y_pred  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "38eb0cc8-f44d-4231-9e7f-11aefbe4f6e8"
      },
      "source": [
        "## Second ensemble technique (EN_optB)\n",
        "Given a set of predictions $X_1, X_2, ..., X_n$, where each $X_i$ has\n",
        "$m=12$ clases, i.e. $X_i = X_{i1}, X_{i2},...,X_{im}$. The algorithm finds the optimal \n",
        "set of weights $w_{11}, w_{12}, ..., w_{nm}$; such that minimizes \n",
        "$log\\_loss(y_T, y_E)$, where $y_E = X_{11}*w_{11} +... + X_{21}*w_{21} + ... + X_{nm}*w_{nm}$ \n",
        "and and $y_T$ is the true solution."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "a22d495e-fad9-46c4-99bc-ba8d5ac3ff31"
      },
      "outputs": [],
      "source": [
        "def objf_ens_optB(w, Xs, y, n_class=12):\n",
        "    \"\"\"\n",
        "    Function to be minimized in the EN_optB ensembler.\n",
        "    \n",
        "    Parameters:\n",
        "    ----------\n",
        "    w: array-like, shape=(n_preds)\n",
        "       Candidate solution to the optimization problem (vector of weights).\n",
        "    Xs: list of predictions to combine\n",
        "       Each prediction is the solution of an individual classifier and has a\n",
        "       shape=(n_samples, n_classes).\n",
        "    y: array-like sahpe=(n_samples,)\n",
        "       Class labels\n",
        "    n_class: int\n",
        "       Number of classes in the problem, i.e. = 12\n",
        "    \n",
        "    Return:\n",
        "    ------\n",
        "    score: Score of the candidate solution.\n",
        "    \"\"\"\n",
        "    #Constraining the weights for each class to sum up to 1.\n",
        "    #This constraint can be defined in the scipy.minimize function, but doing \n",
        "    #it here gives more flexibility to the scipy.minimize function \n",
        "    #(e.g. more solvers are allowed).\n",
        "    w_range = np.arange(len(w))%n_class \n",
        "    for i in range(n_class): \n",
        "        w[w_range==i] = w[w_range==i] / np.sum(w[w_range==i])\n",
        "        \n",
        "    sol = np.zeros(Xs[0].shape)\n",
        "    for i in range(len(w)):\n",
        "        sol[:, i % n_class] += Xs[int(i / n_class)][:, i % n_class] * w[i] \n",
        "        \n",
        "    #Using log-loss as objective function (different objective functions can be used here). \n",
        "    score = log_loss(y, sol)   \n",
        "    return score\n",
        "    \n",
        "\n",
        "class EN_optB(BaseEstimator):\n",
        "    \"\"\"\n",
        "    Given a set of predictions $X_1, X_2, ..., X_n$, where each $X_i$ has\n",
        "    $m=12$ clases, i.e. $X_i = X_{i1}, X_{i2},...,X_{im}$. The algorithm finds the optimal \n",
        "    set of weights $w_{11}, w_{12}, ..., w_{nm}$; such that minimizes \n",
        "    $log\\_loss(y_T, y_E)$, where $y_E = X_{11}*w_{11} +... + X_{21}*w_{21} + ... \n",
        "    + X_{nm}*w_{nm}$ and and $y_T$ is the true solution.\n",
        "    \"\"\"\n",
        "    def __init__(self, n_class=12):\n",
        "        super(EN_optB, self).__init__()\n",
        "        self.n_class = n_class\n",
        "        \n",
        "    def fit(self, X, y):\n",
        "        \"\"\"\n",
        "        Learn the optimal weights by solving an optimization problem.\n",
        "        \n",
        "        Parameters:\n",
        "        ----------\n",
        "        Xs: list of predictions to be ensembled\n",
        "           Each prediction is the solution of an individual classifier and has \n",
        "           shape=(n_samples, n_classes).\n",
        "        y: array-like\n",
        "           Class labels\n",
        "        \"\"\"\n",
        "        Xs = np.hsplit(X, X.shape[1]/self.n_class)\n",
        "        #Initial solution has equal weight for all individual predictions.\n",
        "        x0 = np.ones(self.n_class * len(Xs)) / float(len(Xs)) \n",
        "        #Weights must be bounded in [0, 1]\n",
        "        bounds = [(0,1)]*len(x0)   \n",
        "        #Calling the solver (constraints are directly defined in the objective\n",
        "        #function)\n",
        "        res = minimize(objf_ens_optB, x0, args=(Xs, y, self.n_class), \n",
        "                       method='L-BFGS-B', \n",
        "                       bounds=bounds, \n",
        "                       )\n",
        "        self.w = res.x\n",
        "        return self\n",
        "    \n",
        "    def predict_proba(self, X):\n",
        "        \"\"\"\n",
        "        Use the weights learned in training to predict class probabilities.\n",
        "        \n",
        "        Parameters:\n",
        "        ----------\n",
        "        Xs: list of predictions to be ensembled\n",
        "            Each prediction is the solution of an individual classifier and has \n",
        "            shape=(n_samples, n_classes).\n",
        "            \n",
        "        Return:\n",
        "        ------\n",
        "        y_pred: array_like, shape=(n_samples, n_class)\n",
        "                The ensembled prediction.\n",
        "        \"\"\"\n",
        "        Xs = np.hsplit(X, X.shape[1]/self.n_class)\n",
        "        y_pred = np.zeros(Xs[0].shape)\n",
        "        for i in range(len(self.w)):\n",
        "            y_pred[:, i % self.n_class] += \\\n",
        "                   Xs[int(i / self.n_class)][:, i % self.n_class] * self.w[i]  \n",
        "        return y_pred      "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "92a9fa9b-6bf6-4dc9-9080-448915d911aa"
      },
      "source": [
        "# 2- How to use EN_optA and EN_optB in a 3-layers classification architecture.\n",
        "\n",
        "Somehow similar 3-layers architectures have been previously used in Kaggle competitions\n",
        "(e.g. [here](https://www.kaggle.com/c/otto-group-product-classification-challenge/forums/t/14335/1st-place-winner-solution-gilberto-titericz-stanislav-semenov))\n",
        "\n",
        "### Data\n",
        "For simplicity I am using here synthetic data instead of the original data from Airbnb competition. \n",
        "All the feature engineering step is avoided and it is also easier to play with the number of classes.\n",
        "Moreover, it is easier to change the parameters of the data generation function to study the performance of \n",
        "the algorithms on different types of data.\n",
        "\n",
        "Once the data is generated it is splitted into:\n",
        "\n",
        "- training set: (X_train, y_train)\n",
        "- validation set: (X_valid, y_valid)\n",
        "- test set: (X_test, y_test)\n",
        "\n",
        "### Learning architecture\n",
        "\n",
        " * First layer: I am using 6 classifiers from scikit-learn (Support_Vector_Machines, Logistic_Regression, \n",
        "   Random_Forest, Gradient_Boosting, Extra_Trees_Classifier, K_Nearest_Neighbors). All classifiers are used with \n",
        "   (almost) default parameters. At this level, many other classifiers can be used. \n",
        "   All classifiers are applied twice:\n",
        "     1. Classifiers are trained on (X_train, y_train) and used to predict the class probabilities of (X_valid).\n",
        "     2. Classifiers are trained on (X = (X_train + X_valid), y = (y_train + y_valid)) and used to predict \n",
        "        the class probabilities of (X_test)\n",
        " * Second layer: The predictions from the previous layer on X_valid are concatenated and used to create a new \n",
        "   training set (XV, y_valid). The predictions on X_test are concatenated to create a new test set (XT, y_test). \n",
        "   The two proposed ensemble methods (EN_optA and EN_optB) and their calibrated versions are trained on \n",
        "   (XV, y_valid) and used to predict the class probabilites of (XT).\n",
        " * Third layer: The four prediction from the previous layer are linearly combined using fixed weights."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "51191efb-a210-4d14-a67c-6f6d7a2db5ed"
      },
      "outputs": [],
      "source": [
        "from sklearn.datasets import make_classification\n",
        "from sklearn.linear_model import LogisticRegression\n",
        "from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier, ExtraTreesClassifier\n",
        "from sklearn.svm import SVC\n",
        "from sklearn.neighbors import KNeighborsClassifier\n",
        "from sklearn.calibration import CalibratedClassifierCV\n",
        "from sklearn.cross_validation import train_test_split\n",
        "from sklearn.linear_model import LogisticRegressionCV\n",
        "from xgboost.sklearn import XGBClassifier\n",
        "\n",
        "#fixing random state\n",
        "random_state=1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "c23c08cc-a81c-445b-9d78-f7c8c25524cb"
      },
      "source": [
        "## Generating dataset\n",
        "    \n",
        "Parameters can be changed to explore different types of synthetic data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "f23d5ab3-dc4c-44d7-b69f-8215d4072353"
      },
      "outputs": [],
      "source": [
        "n_classes = 12  # Same number of classes as in Airbnb competition.\n",
        "data, labels = make_classification(n_samples=2000, n_features=100, \n",
        "                                   n_informative=50, n_classes=n_classes, \n",
        "                                   random_state=random_state)\n",
        "\n",
        "#Spliting data into train and test sets.\n",
        "X, X_test, y, y_test = train_test_split(data, labels, test_size=0.2, \n",
        "                                        random_state=random_state)\n",
        "    \n",
        "#Spliting train data into training and validation sets.\n",
        "X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.25, \n",
        "                                                      random_state=random_state)\n",
        "\n",
        "print('Data shape:')\n",
        "print('X_train: %s, X_valid: %s, X_test: %s \\n' %(X_train.shape, X_valid.shape, \n",
        "                                                  X_test.shape))\n",
        "    "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "1ba82de4-ad82-4b16-8ea1-bce1a8793410"
      },
      "source": [
        "## First layer (individual classifiers)\n",
        "All classifiers are applied twice:\n",
        "\n",
        " - Training on (X_train, y_train) and predicting on (X_valid)\n",
        " - Training on (X, y) and predicting on (X_test)\n",
        "    \n",
        "You can add / remove classifiers or change parameter values to see the effect on final results."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "70ef6605-40ba-4c4f-ac4d-edf116917dcc"
      },
      "outputs": [],
      "source": [
        "#Defining the classifiers\n",
        "clfs = {'LR'  : LogisticRegression(random_state=random_state), \n",
        "        'SVM' : SVC(probability=True, random_state=random_state), \n",
        "        'RF'  : RandomForestClassifier(n_estimators=100, n_jobs=-1, \n",
        "                                       random_state=random_state), \n",
        "        'GBM' : GradientBoostingClassifier(n_estimators=50, \n",
        "                                           random_state=random_state), \n",
        "        'ETC' : ExtraTreesClassifier(n_estimators=100, n_jobs=-1, \n",
        "                                     random_state=random_state),\n",
        "        'KNN' : KNeighborsClassifier(n_neighbors=30)}\n",
        "    \n",
        "#predictions on the validation and test sets\n",
        "p_valid = []\n",
        "p_test = []\n",
        "   \n",
        "print('Performance of individual classifiers (1st layer) on X_test')   \n",
        "print('------------------------------------------------------------')\n",
        "   \n",
        "for nm, clf in clfs.items():\n",
        "    #First run. Training on (X_train, y_train) and predicting on X_valid.\n",
        "    clf.fit(X_train, y_train)\n",
        "    yv = clf.predict_proba(X_valid)\n",
        "    p_valid.append(yv)\n",
        "        \n",
        "    #Second run. Training on (X, y) and predicting on X_test.\n",
        "    clf.fit(X, y)\n",
        "    yt = clf.predict_proba(X_test)\n",
        "    p_test.append(yt)\n",
        "       \n",
        "    #Printing out the performance of the classifier\n",
        "    print('{:10s} {:2s} {:1.7f}'.format('%s: ' %(nm), 'logloss  =>', log_loss(y_test, yt)))\n",
        "print('')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "c5911fab-e208-417a-87c3-8bdf361e2060"
      },
      "source": [
        "## Second layer (optimization based ensembles)\n",
        "Predictions on X_valid are used as training set (XV) and predictions on X_test are used as test set (XT). \n",
        "EN_optA, EN_optB and their calibrated versions are applied."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "0ff51335-04ac-4361-89f5-ff6a76fbae67"
      },
      "outputs": [],
      "source": [
        "print('Performance of optimization based ensemblers (2nd layer) on X_test')   \n",
        "print('------------------------------------------------------------')\n",
        "    \n",
        "#Creating the data for the 2nd layer.\n",
        "XV = np.hstack(p_valid)\n",
        "XT = np.hstack(p_test)  \n",
        "        \n",
        "#EN_optA\n",
        "enA = EN_optA(n_classes)\n",
        "enA.fit(XV, y_valid)\n",
        "w_enA = enA.w\n",
        "y_enA = enA.predict_proba(XT)\n",
        "print('{:20s} {:2s} {:1.7f}'.format('EN_optA:', 'logloss  =>', log_loss(y_test, y_enA)))\n",
        "    \n",
        "#Calibrated version of EN_optA \n",
        "cc_optA = CalibratedClassifierCV(enA, method='isotonic')\n",
        "cc_optA.fit(XV, y_valid)\n",
        "y_ccA = cc_optA.predict_proba(XT)\n",
        "print('{:20s} {:2s} {:1.7f}'.format('Calibrated_EN_optA:', 'logloss  =>', log_loss(y_test, y_ccA)))\n",
        "        \n",
        "#EN_optB\n",
        "enB = EN_optB(n_classes) \n",
        "enB.fit(XV, y_valid)\n",
        "w_enB = enB.w\n",
        "y_enB = enB.predict_proba(XT)\n",
        "print('{:20s} {:2s} {:1.7f}'.format('EN_optB:', 'logloss  =>', log_loss(y_test, y_enB)))\n",
        "\n",
        "#Calibrated version of EN_optB\n",
        "cc_optB = CalibratedClassifierCV(enB, method='isotonic')\n",
        "cc_optB.fit(XV, y_valid)\n",
        "y_ccB = cc_optB.predict_proba(XT)  \n",
        "print('{:20s} {:2s} {:1.7f}'.format('Calibrated_EN_optB:', 'logloss  =>', log_loss(y_test, y_ccB)))\n",
        "print('')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "700b7d8c-59c7-4cd4-af69-bc6793207c8c"
      },
      "source": [
        "## Third layer (weighted average)\n",
        "Simple weighted average of the previous 4 predictions."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "1a380d87-4cc9-4ca2-bb3e-3417c13043de"
      },
      "outputs": [],
      "source": [
        "y_3l = (y_enA * 4./9.) + (y_ccA * 2./9.) + (y_enB * 2./9.) + (y_ccB * 1./9.)\n",
        "print('{:20s} {:2s} {:1.7f}'.format('3rd_layer:', 'logloss  =>', log_loss(y_test, y_3l)))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "5971bb93-56e8-4c71-9db8-7644e35eb90d"
      },
      "source": [
        "### Plotting the weights of each ensemble\n",
        "In the case of EN_optA, there is a weight for each prediction and in the case of EN_optB there is \n",
        "a weight for each class for each prediction."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "1c47ff5e-f272-49e9-90d2-db6c20749fb2"
      },
      "outputs": [],
      "source": [
        "from tabulate import tabulate\n",
        "print('               Weights of EN_optA:')\n",
        "print('|---------------------------------------------|')\n",
        "wA = np.round(w_enA, decimals=2).reshape(1,-1)\n",
        "print(tabulate(wA, headers=clfs.keys(), tablefmt=\"orgtbl\"))\n",
        "print('')\n",
        "print('                                    Weights of EN_optB:')\n",
        "print('|-------------------------------------------------------------------------------------------|')\n",
        "wB = np.round(w_enB.reshape((-1,n_classes)), decimals=2)\n",
        "wB = np.hstack((np.array(list(clfs.keys()), dtype=str).reshape(-1,1), wB))\n",
        "print(tabulate(wB, headers=['y%s'%(i) for i in range(n_classes)], tablefmt=\"orgtbl\"))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "24982acf-809b-473f-a399-5b352f8d596d"
      },
      "source": [
        "### Comparing our ensemble results with sklearn LogisticRegression based stacking of classifiers.\n",
        "Both techniques *EN_optA* and *EN_optB* optimizes an objective function. In this experiment I am using the multi-class \n",
        "logloss as objective function. Therefore, the two proposed methods basically become implementations of LogisticRegression.\n",
        "The following code allows to compare the results of sklearn implementation of LogisticRegression with the proposed ensembles."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "33483233-097b-49b7-993e-76ff2c6965d8"
      },
      "outputs": [],
      "source": [
        "#By default the best C parameter is obtained with a cross-validation approach, doing grid search with\n",
        "#10 values defined in a logarithmic scale between 1e-4 and 1e4.\n",
        "#Change parameters to see how they affect the final results.\n",
        "lr = LogisticRegressionCV(Cs=10, dual=False, fit_intercept=True, \n",
        "                          intercept_scaling=1.0, max_iter=100,\n",
        "                          multi_class='ovr', n_jobs=1, penalty='l2', \n",
        "                          random_state=random_state,\n",
        "                          solver='lbfgs', tol=0.0001)\n",
        "\n",
        "lr.fit(XV, y_valid)\n",
        "y_lr = lr.predict_proba(XT)\n",
        "print('{:20s} {:2s} {:1.7f}'.format('Log_Reg:', 'logloss  =>', log_loss(y_test, y_lr)))\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "02a25c52-1d0f-4d4f-bb49-6f9c78666d8b"
      },
      "source": [
        "### Is there any parameters configuration for LogisticRegression that produces better results than the proposed ensemble techniques?\n",
        "I wasn't able to find such parameter configuration for a problem with 12 number of classes."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "0a290560-536c-4213-9886-e8707cfdc810"
      },
      "source": [
        "# 3- Comparison of the ensemble techniques on problems with different number of classes\n",
        "Let's explore how the different ensemble techniques perform according to the number of classes in the problem.\n",
        "We generate different dataset with different number of classes (e.g. from 3 to 15 classes) and compare the result of \n",
        "the different ensembling methods."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "f72442e8-7d56-4bf4-ba29-88712c915a4f"
      },
      "outputs": [],
      "source": [
        "#For each value in classes, a dataset with that number of classes will be created. \n",
        "classes = range(3, 15)\n",
        "\n",
        "ll_sc = []  #to store logloss of individual classifiers\n",
        "ll_eA = []  #to store logloss of EN_optA ensembler\n",
        "ll_eB = []  #to store logloss of EN_optB ensembler\n",
        "ll_e3 = []  #to store logloss of the third-layer ensembler (method used for submission in the competition).\n",
        "ll_lr = []  #to store logloss of LogisticRegression as 2nd layer ensembler.\n",
        "ll_gb = []  #to store logloss of GradientBoosting as 2nd layer ensembler.\n",
        "\n",
        "#Same code as above for generating the dataset, applying the 3-layer learning architecture and copmparing with  \n",
        "#LogisticRegression and GradientBoosting based ensembles. \n",
        "#The code is applied to each independent problem/dataset (each dataset with a different number of classes).\n",
        "for i in classes:\n",
        "    print('Working on dataset with n_classes: %s' %(i))\n",
        "    n_classes=i\n",
        "    \n",
        "    #Generating the data\n",
        "    data, labels = make_classification(n_samples=2000, n_features=100, \n",
        "                                       n_informative=50, n_classes=n_classes,\n",
        "                                       random_state=random_state)\n",
        "    X, X_test, y, y_test = train_test_split(data, labels, test_size=0.2, \n",
        "                                            random_state=random_state)\n",
        "    X_train, X_valid, y_train, y_valid = train_test_split(X, y, \n",
        "                                              test_size=0.25, \n",
        "                                              random_state=random_state)\n",
        "    \n",
        "    #First layer\n",
        "    clfs = [LogisticRegression(random_state=random_state), \n",
        "            SVC(probability=True, random_state=random_state), \n",
        "            RandomForestClassifier(n_estimators=100, n_jobs=-1, \n",
        "                                   random_state=random_state), \n",
        "            GradientBoostingClassifier(n_estimators=50, \n",
        "                                       random_state=random_state), \n",
        "            ExtraTreesClassifier(n_estimators=100, n_jobs=-1, \n",
        "                                 random_state=random_state), \n",
        "            KNeighborsClassifier(n_neighbors=30, n_jobs=-1)]\n",
        "    p_valid = []\n",
        "    p_test = []\n",
        "    for clf in clfs:\n",
        "        #First run\n",
        "        clf.fit(X_train, y_train)\n",
        "        yv = clf.predict_proba(X_valid)\n",
        "        p_valid.append(yv)\n",
        "        #Second run\n",
        "        clf.fit(X, y)\n",
        "        yt = clf.predict_proba(X_test)\n",
        "        p_test.append(yt)\n",
        "        #Saving the logloss score\n",
        "        ll_sc.append(log_loss(y_test, yt))\n",
        "\n",
        "    #Second layer\n",
        "    XV = np.hstack(p_valid)\n",
        "    XT = np.hstack(p_test)  \n",
        "    \n",
        "    enA = EN_optA(n_classes)   #EN_optA\n",
        "    enA.fit(XV, y_valid)\n",
        "    y_enA = enA.predict_proba(XT)    \n",
        "    ll_eA.append(log_loss(y_test, y_enA))  #Saving the logloss score\n",
        "    \n",
        "    cc_optA = CalibratedClassifierCV(enA, method='isotonic') #Calibrated version of EN_optA \n",
        "    cc_optA.fit(XV, y_valid)\n",
        "    y_ccA = cc_optA.predict_proba(XT)\n",
        "    \n",
        "    enB = EN_optB(n_classes)   #EN_optB\n",
        "    enB.fit(XV, y_valid)\n",
        "    y_enB = enB.predict_proba(XT)   #Saving the logloss score\n",
        "    ll_eB.append(log_loss(y_test, y_enB))\n",
        "    \n",
        "    cc_optB = CalibratedClassifierCV(enB, method='isotonic') #Calibrated version of EN_optB \n",
        "    cc_optB.fit(XV, y_valid)\n",
        "    y_ccB = cc_optB.predict_proba(XT) \n",
        "    \n",
        "    #Third layer\n",
        "    y_3l = (y_enA * 4./9.) + (y_ccA * 2./9.) + (y_enB * 2./9.) + (y_ccB * 1./9.)\n",
        "    ll_e3.append(log_loss(y_test, y_3l))   #Saving the logloss score\n",
        "    \n",
        "    #Logistic regresson\n",
        "    lr = LogisticRegressionCV(Cs=10, dual=False, fit_intercept=True,\n",
        "                              intercept_scaling=1.0, max_iter=100,\n",
        "                              multi_class='ovr', n_jobs=1, penalty='l2',\n",
        "                              random_state=random_state,\n",
        "                              solver='lbfgs', tol=0.0001)\n",
        "    lr.fit(XV, y_valid)\n",
        "    y_lr = lr.predict_proba(XT)\n",
        "    ll_lr.append(log_loss(y_test, y_lr))   #Saving the logloss score\n",
        "    \n",
        "    #Gradient boosting\n",
        "    xgb = XGBClassifier(max_depth=5, learning_rate=0.1, \n",
        "                        n_estimators=10000, objective='multi:softprob', \n",
        "                        seed=random_state)\n",
        "    #Computing best number of iterations on an internal validation set\n",
        "    XV_train, XV_valid, yv_train, yv_valid = train_test_split(XV, y_valid,\n",
        "                                             test_size=0.15, random_state=random_state)\n",
        "    xgb.fit(XV_train, yv_train, eval_set=[(XV_valid, yv_valid)], \n",
        "            eval_metric='mlogloss', \n",
        "            early_stopping_rounds=15, verbose=False)\n",
        "    xgb.n_estimators = xgb.best_iteration\n",
        "    \n",
        "    xgb.fit(XV, y_valid)\n",
        "    y_gb = xgb.predict_proba(XT)\n",
        "    ll_gb.append(log_loss(y_test, y_gb)) #Saving the logloss score \n",
        "\n",
        "ll_sc = np.array(ll_sc).reshape(-1, len(clfs)).T \n",
        "ll_eA = np.array(ll_eA) \n",
        "ll_eB = np.array(ll_eB) \n",
        "ll_e3 = np.array(ll_e3)\n",
        "ll_lr = np.array(ll_lr) \n",
        "ll_gb = np.array(ll_gb)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "b9156d0f-8f4d-4bfd-be3b-92d9307914ad"
      },
      "source": [
        "## Plotting the results\n",
        "Notice that sklearn LogisticRegression and XGBoost produce better results for problems with few classes, but as the number of classes increases\n",
        "the proposed ensembling methods outperform LogisticRegression and XGBoost. Again the question here is whether it is possible to fine-tune\n",
        "LogisticRegression (or XGBoost) to produce better (or comparable) results than the ones produced by EN_optA, EN_optB on problems with high number of clases.\n",
        "It can also be noticed that the *3rd-layer* ensemble always produces better results than the 2nd-layer ensemblers\n",
        "(e.g. EN_optA, EN_optB)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "7e9e0a28-db58-4904-b40d-6c800d6a11aa"
      },
      "outputs": [],
      "source": [
        "import matplotlib.pylab as plt\n",
        "plt.figure(figsize=(10,7))\n",
        "\n",
        "plt.plot(classes, ll_sc[0], color='black', label='Single_Classifiers')\n",
        "for i in range(1, 6):\n",
        "    plt.plot(classes, ll_sc[i], color='black')\n",
        "    \n",
        "plt.plot(classes, ll_lr, 'bo-', label='EN_LogisticRegression')\n",
        "plt.plot(classes, ll_gb, 'mo-', label='EN_XGBoost')\n",
        "plt.plot(classes, ll_eA, 'yo-', label='EN_optA')\n",
        "plt.plot(classes, ll_eB, 'go-', label='EN_optB')\n",
        "plt.plot(classes, ll_e3, 'ro-', label='EN_3rd_layer')\n",
        "\n",
        "plt.title('Log-loss of the different models for different number of classes.')\n",
        "plt.xlabel('Number of classes')\n",
        "plt.ylabel('Log-loss')\n",
        "plt.grid(True)\n",
        "plt.legend(loc=4)\n",
        "plt.show()"
      ]
    }
  ],
  "metadata": {
    "_change_revision": 0,
    "_is_fork": false,
    "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.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}