Distribution regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from numpy import exp, sqrt, dot\n",
    "from scipy.spatial.distance import cdist\n",
    "\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.preprocessing import scale\n",
    "from sklearn.model_selection import KFold\n",
    "from sklearn.svm import SVR\n",
    "\n",
    "from hyperopt import STATUS_OK, hp, fmin, tpe, Trials, space_eval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Import data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def load_data():\n",
    "    full_data = pd.read_csv(\"Data/X.csv\")\n",
    "    train_y = pd.read_csv(\"Data/y_train.csv\")\n",
    "    # Rename columns to something more interpretable\n",
    "    columns = ([\"reflectance_\" + str(i) for i in range(7)]\n",
    "               + [\"solar_\" + str(i) for i in range(5)] + [\"id\"])\n",
    "    full_data.columns = columns\n",
    "    \n",
    "    # Move ID column to the beginning\n",
    "    id_column = full_data[\"id\"]\n",
    "    full_data.drop(labels=[\"id\"], axis=1, inplace = True)\n",
    "    full_data.insert(0, \"id\", id_column)\n",
    "    \n",
    "    # Add the target value column to the training part\n",
    "    # in full_data\n",
    "    split = 98000\n",
    "    y_id_dict = train_y.set_index(\"id\")[\"y\"].to_dict()\n",
    "    full_data.loc[:(split-1), \"y\"] = full_data.loc[:(split-1), \"id\"].map(y_id_dict)\n",
    "    \n",
    "    # Split into training and testing data\n",
    "    train, test = full_data[:split], full_data[split:]\n",
    "    return (train, test)\n",
    "\n",
    "train, test = load_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Set parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "random_seed = 8\n",
    "# Set folds for out-of-fold prediction\n",
    "n_folds = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "cols_excl = [\"id\", \"y\"]\n",
    "cols_orig = [c for c in train.columns if c not in cols_excl]\n",
    "\n",
    "# Standardise data can make training faster and reduce\n",
    "# the chances of getting stuck in local optima\n",
    "train[cols_orig] = scale(train[cols_orig])\n",
    "test[cols_orig] = scale(test[cols_orig])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 1. Distribution regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Kernel mean embedding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Instead of fitting a model to the instances, the idea of distribution regression is to find a regression on the underlying probability distributions the instances come from. It is based on the assumption that the data is $\\{(x_i, y_i)\\}_{i=1}^{n}$ with:\n",
    "\n",
    "* $n$ the number of bags in the dataset ;\n",
    "* $x_i$ the probability distribution of bag $i$ ;\n",
    "* $y_i$ is the aerosol optical depth of bag $i$.\n",
    "\n",
    "However, $x_i$ is not observed: for each bag $i$, the $100$ instances $x_{i,l}$, $l=1,...,100$, are samples from the distribution $x_i$. Our dataset is thus  $\\{(\\{x_{i,l}\\}_{l=1}^{100}, y_i)\\}_{i=1}^{n}$ and we want to find a mapping $\\hat{f}$ that will best predict unseen bags.\n",
    "\n",
    "The mapping  $\\hat{f}$ on $\\{(\\{x_{i,l}\\}_{l=1}^{100}, y_i)\\}_{i=1}^{n}$ will try to learn the relationship between the true distributions $\\{x_i\\}_{i=1}^{n}$ and the target values $\\{y_i\\}_{i=1}^{n}$. To achieve that, the information of the 100 instances in each bag has to be summarised whilst losing the less information possible. The aggregated approach that simply computes the mean of the features for each bag is an example of information summary, yet plenty of data is lost that way.\n",
    "A better way to represent each bag is via kernel mean embedding:\n",
    "$$\\mu_{\\hat{x}_i} = \\frac{1}{100}\\sum_{l=1}^{100} k(\\cdot, x_{i,l})$$\n",
    "\n",
    "Each bag is represented as a linear combination of kernels, and with the right choice of kernel, the lost information can be very negligible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Kernel Ridge Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We now want to find $\\hat{f}$ that minimises the following regularised least square problem:\n",
    "$$ \\underset{f}{arg min} \\sum_{i=1}^{n} (f(\\mu_{\\hat{x}_i}) - y_i)^2 + \\lambda \\Vert f \\Vert^2$$\n",
    "\n",
    "with $\\lambda>0$ the L2 regularisation parameter.\n",
    "\n",
    "In kernel ridge regression, $f$ is interpreted as a linear combination of feature space mappings $\\phi$ of the data points $\\mu_{\\hat{x}_i}$:\n",
    "$$ f = \\sum_{i=1}^{n} \\alpha_i \\phi(\\mu_{\\hat{x}_i} ) $$\n",
    "\n",
    "The equation thus becomes:\n",
    "$$ \\underset{\\alpha}{arg min} (\\Vert y -K\\alpha \\Vert^2 + \\lambda \\alpha^T K \\alpha)$$\n",
    "with :\n",
    "* $K(i,j) = k'(\\mu_{\\hat{x}_i} , \\mu_{\\hat{x}_j})$ for $i,j=1..n$ ;\n",
    "* $k'$ another kernel.\n",
    "\n",
    "By differentiating with respect to $\\alpha$ and setting it to zero:\n",
    "$$ \\alpha^{*} = (K + \\lambda I_n)^{-1}y $$\n",
    "\n",
    "For the sake of simplicity and because the results proved to be reasonably good, we set $k'$ as the linear kernel and as a result:\n",
    "$$ K(i,j) = \\frac{1}{100^2} \\sum_{l,k=1}^{100} k(x_{i,l} , x_{j,k})$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "class Kernel(object):\n",
    "    \"\"\" Kernel class from Zoltan Szabo\n",
    "        giving the kernel mean embedding.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, par=None):\n",
    "        \"\"\" Initialization.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        par : dictionary, optional\n",
    "              Name of the kernel and its parameters (default is\n",
    "              {\"name\": \"RBF\", \"sigma\": 1}). The name of the kernel comes\n",
    "              from \"RBF\", \"exponential\", \"Cauchy\", \"student\", \"Matern3p2\",\n",
    "              \"Matern5p2\", \"polynomial\", \"ratquadr\" (rational quadratic),\n",
    "              \"invmquadr\" (inverse multiquadr).\n",
    "        \"\"\"\n",
    "        if par is None:\n",
    "            par = {\"name\": \"RBF\", \"sigma\": 1}\n",
    "\n",
    "        name = par[\"name\"]\n",
    "        self.name = name\n",
    "\n",
    "        # other attributes:\n",
    "        if name == \"RBF\" or name == \"exponential\" or name == \"Cauchy\":\n",
    "            self.sigma = par[\"sigma\"]\n",
    "        elif name == \"student\":\n",
    "            self.d = par[\"d\"]\n",
    "        elif name == \"Matern3p2\" or name == \"Matern5p2\":\n",
    "            self.l = par[\"l\"]\n",
    "        elif name == \"polynomial\":\n",
    "            self.c = par[\"c\"]\n",
    "            self.exponent = par[\"exponent\"]\n",
    "        elif name == \"ratquadr\" or name == \"invmquadr\":\n",
    "            self.c = par[\"c\"]\n",
    "        else:\n",
    "            raise Exception(\"kernel=?\")\n",
    "\n",
    "    def gram_matrix(self, y1, y2):\n",
    "        \"\"\"  Compute the Gram matrix = [k(y1[i,:], y2[j,:])]; i, j: running.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        y1 : (number of samples1, dimension)-ndarray\n",
    "             One row of y1 corresponds to one sample.\n",
    "        y2 : (number of samples2, dimension)-ndarray\n",
    "             One row of y2 corresponds to one sample.\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        g : ndarray.\n",
    "            Gram matrix of y1 and y2.\n",
    "        \"\"\"\n",
    "\n",
    "        if self.name == \"RBF\":\n",
    "            sigma = self.sigma\n",
    "            g = cdist(y1, y2)\n",
    "            g = exp(-g ** 2 / (2 * sigma ** 2))\n",
    "        elif self.name == \"exponential\":\n",
    "            sigma = self.sigma\n",
    "            g = cdist(y1, y2)\n",
    "            g = exp(-g / (2 * sigma ** 2))\n",
    "        elif self.name == \"Cauchy\":\n",
    "            sigma = self.sigma\n",
    "            g = cdist(y1, y2)\n",
    "            g = 1 / (1 + g ** 2 / sigma ** 2)\n",
    "        elif self.name == \"student\":\n",
    "            d = self.d\n",
    "            g = cdist(y1, y2)\n",
    "            g = 1 / (1 + g ** d)\n",
    "        elif self.name == \"Matern3p2\":\n",
    "            l = self.l\n",
    "            g = cdist(y1, y2) \n",
    "            g = (1 + sqrt(3) * g / l) * exp(-sqrt(3) * g / l)\n",
    "        elif self.name == \"Matern5p2\":\n",
    "            l = self.l\n",
    "            g = cdist(y1, y2)\n",
    "            g = (1 + sqrt(5) * g / l + 5 * g ** 2 / (3 * l ** 2)) * \\\n",
    "                exp(-sqrt(5) * g / l)\n",
    "        elif self.name == \"polynomial\":\n",
    "            c = self.c\n",
    "            exponent = self.exponent\n",
    "            g = (dot(y1, y2.T) + c) ** exponent\n",
    "        elif self.name == \"ratquadr\":\n",
    "            c = self.c\n",
    "            g = cdist(y1, y2) ** 2\n",
    "            g = 1 - g / (g + c)\n",
    "        elif self.name == \"invmquadr\":\n",
    "            c = self.c\n",
    "            g = cdist(y1, y2)\n",
    "            g = 1 / sqrt(g ** 2 + c ** 2)\n",
    "        else:\n",
    "            raise Exception(\"kernel=?\")\n",
    "\n",
    "        return g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Compute the linear kernel product of \n",
    "# the mean embedding of X1 and X2\n",
    "# denoted as K(i, j) above\n",
    "def mean_embedding(X1, X2, kernel):\n",
    "    k = Kernel(kernel)\n",
    "    gram_mat = k.gram_matrix(X1, X2)\n",
    "    # Number of instances in the bag\n",
    "    N = float(gram_mat.shape[0])\n",
    "    mu_X1_X2 = gram_mat.ravel().sum() / N**2\n",
    "    return (mu_X1_X2)\n",
    "\n",
    "# Return a symmetrised matrix\n",
    "def symmetrise(A):\n",
    "    return(A + A.T - np.diag(A.diagonal()))\n",
    "\n",
    "# Compute the Gram matrix K given the kernel and \n",
    "# the smoothing parameter theta\n",
    "def compute_gram(df, kernel, theta):\n",
    "    nb_bag = df[\"id\"].nunique()\n",
    "    K_matrix = np.zeros((nb_bag, nb_bag))\n",
    "    print(\"Computing {0} Gram matrix for theta={1}:\".format(kernel, theta))\n",
    "    for i in range(nb_bag):\n",
    "        if (i%50 == 0):\n",
    "            print(\"Bag number: {0}\". format(i))\n",
    "        \n",
    "        for j in range(i+1):\n",
    "            # Compute mean embedding\n",
    "            X1 = df.loc[train[\"id\"] == (i+1), cols_orig].values\n",
    "            X2 = df.loc[train[\"id\"] == (j+1), cols_orig].values\n",
    "\n",
    "            K_matrix[i, j] = mean_embedding(X1, X2, {'name': kernel, 'sigma': theta})\n",
    "            \n",
    "    return symmetrise(K_matrix)\n",
    "        \n",
    "#K_cauchy = compute_gram(train, \"Cauchy\", 2**4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Class for kernel ridge regression\n",
    "class RidgeRegression(object):\n",
    "    def __init__(self, l2_reg):\n",
    "        self.l2_reg = l2_reg\n",
    "\n",
    "    def fit(self, G, y):\n",
    "        # Train size\n",
    "        n_train = G.shape[0]\n",
    "        ridge_mat = G + (self.l2_reg * n_train) * np.identity(n_train)\n",
    "        self.ridge_mat = ridge_mat\n",
    "        # Shape of y_train is (1, n_train)\n",
    "        self.y_train = y\n",
    "\n",
    "    def predict(self, G_test):\n",
    "        y_test_hat = self.y_train.dot(np.linalg.solve(self.ridge_mat, G_test))\n",
    "        return y_test_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Kernel selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "A kernel is characterised by a parameter we will call $\\theta$ and the ridge regression depends on the L2 regularisation $\\lambda$. Through cross-validation, we selected the kernels giving the most stable validation loss. They are given below with their associated parameters:\n",
    "\n",
    "* __Cauchy__: \n",
    "$$k_C(a,b) = \\dfrac{1}{1 + \\dfrac{\\Vert a-b\\Vert_2^2}{\\theta^2}}, \\quad\\theta_C = 16, \\quad\\lambda_C = 2^{-23} $$\n",
    "* __Matérn 5/2__: \n",
    "$$k_M(a,b) = \\left(1 + \\dfrac{\\sqrt{5}\\Vert a-b\\Vert_2^2}{\\theta} + \\dfrac{5\\Vert a-b\\Vert_2^2}{3\\theta^2} \\right)e^{-\\dfrac{\\sqrt{5}\\Vert a-b\\Vert_2^2}{\\theta}}, \\quad\\theta_M = 64, \\quad\\lambda_M = 2^{-31} $$\n",
    "* __Rational quadratic__: \n",
    "$$k_r(a,b) = 1 - \\dfrac{\\Vert a-b\\Vert_2^2}{\\Vert a-b\\Vert_2^2 + \\theta}, \\quad\\theta_r = 512, \\quad\\lambda_r = 2^{-26}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 2. Stacking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We will then map the features in the three spaces that describes the data in different ways. Each kernel ridge regression gives a prediction of the labels and combining them might give a better result for three reasons:\n",
    "\n",
    "* __Statistical reason__: we might not have enough data and even if each model $h_i$ performs well on the training set, the true model $f$ might still be not reached ;\n",
    "* __Computational reason__: each model $h_i$ only finds a local optima ;\n",
    "* __Representational reason__: the true model is out of the representation of functions we're considering.\n",
    "\n",
    "Combining our model might take us a step closer to finding the true model $f$. The ensembling technique we used was __out-of-fold stacking__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Out-of-fold prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "In the __first stage__, out-of-fold prediction is applied to ensure that each first-layer regressor does not overfit by predicting on data already seen. For each regressor, we iteratively separate the training data in $N$ folds ($N=5$ in our model), and then use N-1 folds to train the model and then predict the target value of the remaining fold. To create the new testing set, the average of the predictions of each fold is taken."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# G_train and G_test are pandas dataframes\n",
    "# krr is a kernel ridge regression\n",
    "def oof_prediction(krr, G_train, y_train, G_test, n_folds, random_seed):\n",
    "    kf = KFold(n_splits=n_folds, shuffle=True, random_state=random_seed)\n",
    "    n_train = G_train.shape[0]\n",
    "    n_test = G_test.shape[1]\n",
    "    oof_train = np.zeros(n_train)\n",
    "    oof_test = np.zeros(n_test)\n",
    "    oof_test_folds = np.zeros((n_test, n_folds))\n",
    "\n",
    "    for i, (train_index, test_index) in enumerate(kf.split(G_train)):\n",
    "        G_tr = G_train.loc[train_index, train_index].values\n",
    "        y_tr = y_train[train_index].reshape((1, -1))\n",
    "        G_te = G_train.loc[train_index, test_index].values\n",
    "\n",
    "        krr.fit(G_tr, y_tr)\n",
    "        oof_train[test_index] = krr.predict(G_te)\n",
    "        G_test_partial = G_test.loc[train_index, :]\n",
    "        oof_test_folds[:, i] = krr.predict(G_test_partial.values)\n",
    "\n",
    "    oof_test = oof_test_folds.mean(axis=1)\n",
    "    return oof_train, oof_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "nb_bags_train = 980\n",
    "# Create a vector with the unique values of y for each ID.\n",
    "y_train = train.groupby(\"id\")[\"y\"].median().values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Load Gram matrices\n",
    "def load_gram(csv_file, nb_bags_train):\n",
    "    # Import data\n",
    "    G = pd.read_csv(csv_file, header=None)\n",
    "    idx_train = nb_bags_train - 1\n",
    "    idx_test = nb_bags_train\n",
    "    G_train = G.loc[:idx_train, :idx_train]\n",
    "    G_test = G.loc[:idx_train, idx_test:]\n",
    "    return (G_train, G_test)\n",
    "\n",
    "# Define models and import Gram matrices\n",
    "# Cauchy\n",
    "l2_reg_cauchy = 2**(-23)\n",
    "cauchy = RidgeRegression(l2_reg_cauchy)\n",
    "G_train_cauchy, G_test_cauchy = load_gram(\"kernels_me/Cauchy_16.csv\", nb_bags_train)\n",
    "\n",
    "# Matern 5/2\n",
    "l2_reg_matern_52 = 2**(-31)\n",
    "matern_52 = RidgeRegression(l2_reg_matern_52)\n",
    "G_train_matern_52, G_test_matern_52 = load_gram(\"kernels_me/Matern_52_64.csv\", nb_bags_train)\n",
    "\n",
    "# Rational quadratic\n",
    "l2_reg_rquadr = 2**(-26)\n",
    "rquadr = RidgeRegression(l2_reg_rquadr)\n",
    "G_train_rquadr, G_test_rquadr = load_gram(\"kernels_me/rquadr_512.csv\", nb_bags_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training is finished.\n"
     ]
    }
   ],
   "source": [
    "# Create OOF train and test predictions\n",
    "# Cauchy\n",
    "cauchy_oof_train, cauchy_oof_test = oof_prediction(cauchy, G_train_cauchy, \n",
    "                                                   y_train, G_test_cauchy,\n",
    "                                                   n_folds, random_seed)\n",
    "# Matern 5/2\n",
    "matern_52_oof_train, matern_52_oof_test = oof_prediction(matern_52, G_train_matern_52, \n",
    "                                                         y_train, G_test_matern_52,\n",
    "                                                         n_folds, random_seed)\n",
    "# Rational quadratic\n",
    "rquadr_oof_train, rquadr_oof_test = oof_prediction(rquadr, G_train_rquadr, \n",
    "                                                   y_train, G_test_rquadr,\n",
    "                                                   n_folds, random_seed)\n",
    "print(\"Training is finished.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Second stage prediction with SVR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "In the __second stage__, a Support Vector Regression is fed with the different predictions from the kernel ridge regression to predict the target value $y$.\n",
    "The SVR uses these predictions to compute the optimal weights assigned to each kernel regression and we might hope to find a better optimum to approximate the true regression function $f$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>cauchy</th>\n",
       "      <th>matern_52</th>\n",
       "      <th>rquadr</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-3.553306</td>\n",
       "      <td>-3.497803</td>\n",
       "      <td>-3.584285</td>\n",
       "      <td>-3.998082</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-4.238140</td>\n",
       "      <td>-4.246409</td>\n",
       "      <td>-4.246125</td>\n",
       "      <td>-4.137141</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-2.585482</td>\n",
       "      <td>-2.645614</td>\n",
       "      <td>-2.563139</td>\n",
       "      <td>-2.694732</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-3.893719</td>\n",
       "      <td>-3.965724</td>\n",
       "      <td>-3.902051</td>\n",
       "      <td>-3.296275</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-3.684191</td>\n",
       "      <td>-3.718706</td>\n",
       "      <td>-3.654680</td>\n",
       "      <td>-3.181391</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     cauchy  matern_52    rquadr         y\n",
       "0 -3.553306  -3.497803 -3.584285 -3.998082\n",
       "1 -4.238140  -4.246409 -4.246125 -4.137141\n",
       "2 -2.585482  -2.645614 -2.563139 -2.694732\n",
       "3 -3.893719  -3.965724 -3.902051 -3.296275\n",
       "4 -3.684191  -3.718706 -3.654680 -3.181391"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Building the new data frames using the\n",
    "# of out-of-fold predictions\n",
    "kernel_train = pd.DataFrame({'cauchy': cauchy_oof_train,\n",
    "                             'matern_52': matern_52_oof_train,\n",
    "                             'rquadr': rquadr_oof_train})\n",
    "\n",
    "kernel_train[\"y\"] = y_train\n",
    "\n",
    "kernel_test = pd.DataFrame({'cauchy': cauchy_oof_test,\n",
    "                            'matern_52': matern_52_oof_test,\n",
    "                            'rquadr': rquadr_oof_test})\n",
    "\n",
    "cols_excl_kernel = [\"y\"]\n",
    "cols_kernel = [c for c in kernel_train.columns if c not in cols_excl_kernel]\n",
    "\n",
    "kernel_train.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Tuning SVR's parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Root mean squared error metric\n",
    "def RMSE(y, y_hat):\n",
    "    out = np.sqrt(mean_squared_error(y.reshape((-1,)), y_hat.reshape((-1,))))\n",
    "    return (out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def scoring_function(parameters):\n",
    "    print(\"Training the model with parameters: \")\n",
    "    print(parameters)\n",
    "    \n",
    "    # Run several KFold shuffles and take the mean RMSE\n",
    "    average_RMSE = []\n",
    "    nb_run = 10\n",
    "    \n",
    "    for m in range(nb_run):\n",
    "        KFold_RMSE = 0.0\n",
    "        n_splits = 5\n",
    "\n",
    "        kf = KFold(n_splits=n_splits, shuffle=True, random_state=(random_seed+m))\n",
    "        nb_fold = 0\n",
    "        for train_index, validation_index in kf.split(kernel_train):\n",
    "            nb_fold += 1\n",
    "            train_fold, validation_fold = kernel_train.loc[train_index], kernel_train.loc[validation_index]\n",
    "\n",
    "            svr = SVR(C=parameters[\"C\"], epsilon=parameters[\"epsilon\"])\n",
    "            svr.fit(train_fold[cols_kernel], train_fold[\"y\"])\n",
    "\n",
    "            y_hat_test = svr.predict(validation_fold[cols_kernel])\n",
    "            RMSE_test = RMSE(y_hat_test, validation_fold[\"y\"].values)\n",
    "\n",
    "            KFold_RMSE += RMSE_test\n",
    "\n",
    "        KFold_RMSE /= n_splits\n",
    "        \n",
    "        average_RMSE.append(KFold_RMSE)\n",
    "    \n",
    "    average_RMSE = np.array(average_RMSE)\n",
    "\n",
    "    print(\"Cross-validation score: {0} +/- {1}\\n\".format(average_RMSE.mean(),\n",
    "                                                         2*average_RMSE.std()))\n",
    "    \n",
    "    return {\"loss\": average_RMSE.mean(), \"status\": STATUS_OK}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the model with parameters: \n",
      "{'epsilon': 0.050000000000000003, 'C': 1.5}\n",
      "Cross-validation score: 0.686196784675 +/- 0.00310294556159\n",
      "\n",
      "Training the model with parameters: \n",
      "{'epsilon': 0.10000000000000001, 'C': 0.5}\n",
      "Cross-validation score: 0.688370690304 +/- 0.00292268709048\n",
      "\n",
      "Training the model with parameters: \n",
      "{'epsilon': 0.050000000000000003, 'C': 2.0}\n",
      "Cross-validation score: 0.685749746794 +/- 0.00324723589177\n",
      "\n",
      "Training the model with parameters: \n",
      "{'epsilon': 0.20000000000000001, 'C': 1.5}\n",
      "Cross-validation score: 0.684769970815 +/- 0.00224621885902\n",
      "\n",
      "Training the model with parameters: \n",
      "{'epsilon': 0.15000000000000002, 'C': 2.5}\n",
      "Cross-validation score: 0.685005070367 +/- 0.00345812467647\n",
      "\n",
      "Training the model with parameters: \n",
      "{'epsilon': 0.10000000000000001, 'C': 1.0}\n",
      "Cross-validation score: 0.686645044299 +/- 0.00267387859763\n",
      "\n",
      "Training the model with parameters: \n",
      "{'epsilon': 0.15000000000000002, 'C': 1.0}\n",
      "Cross-validation score: 0.685922402595 +/- 0.0022856147877\n",
      "\n",
      "Training the model with parameters: \n",
      "{'epsilon': 0.15000000000000002, 'C': 0.5}\n",
      "Cross-validation score: 0.687533030786 +/- 0.00241636467058\n",
      "\n",
      "Training the model with parameters: \n",
      "{'epsilon': 0.10000000000000001, 'C': 1.5}\n",
      "Cross-validation score: 0.686081435086 +/- 0.00327535796236\n",
      "\n",
      "Training the model with parameters: \n",
      "{'epsilon': 0.050000000000000003, 'C': 2.0}\n",
      "Cross-validation score: 0.685749746794 +/- 0.00324723589177\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Grid to pick parameters from.\n",
    "parameters_grid = {\"C\":       hp.choice(\"C\", np.arange(0.5, 3, 0.5)),\n",
    "                   \"epsilon\": hp.choice(\"epsilon\", np.arange(0.05, 0.25, 0.05))\n",
    "                  }\n",
    "# Record the information about the cross-validation.\n",
    "trials = Trials()\n",
    "\n",
    "best = fmin(scoring_function, parameters_grid, algo=tpe.suggest, max_evals=10, \n",
    "            trials=trials)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.684769970815494"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(trials.losses())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>C</th>\n",
       "      <th>epsilon</th>\n",
       "      <th>score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.5</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.68477</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     C  epsilon    score\n",
       "0  1.5      0.2  0.68477"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Save the best parameters as a csv.\n",
    "best_parameters = pd.DataFrame({key: [value] for (key, value) in \n",
    "                                zip(space_eval(parameters_grid, best).keys(),\n",
    "                                    space_eval(parameters_grid, best).values())})\n",
    "# Add the corresponding score.\n",
    "best_parameters[\"score\"] = min(trials.losses())\n",
    "best_parameters.to_csv(\"Parameters/best_parameters_SVR.csv\", encoding=\"utf-8\", index=False)\n",
    "\n",
    "best_parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>C</th>\n",
       "      <th>epsilon</th>\n",
       "      <th>score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.5</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.68477</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     C  epsilon    score\n",
       "0  1.5      0.2  0.68477"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "best_parameters = pd.read_csv(\"Parameters/best_parameters_SVR.csv\", encoding=\"utf-8\")\n",
    "best_parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SVR(C=1.5, cache_size=200, coef0=0.0, degree=3, epsilon=0.20000000000000001,\n",
       "  gamma='auto', kernel='rbf', max_iter=-1, shrinking=True, tol=0.001,\n",
       "  verbose=False)"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "svr = SVR(C=best_parameters[\"C\"][0],\n",
    "          epsilon=best_parameters[\"epsilon\"][0])\n",
    "svr.fit(kernel_train[cols_kernel], y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.67393208259018023"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Training error\n",
    "RMSE(svr.predict(kernel_train[cols_kernel]), y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Prediction\n",
    "y_hat_test = svr.predict(kernel_test[cols_kernel])\n",
    "\n",
    "test_pred = test.groupby(\"id\")[[\"y\"]].mean().reset_index()\n",
    "test_pred[\"y\"] = y_hat_test\n",
    "test_pred.columns = [\"Id\", \"y\"]\n",
    "\n",
    "# Save as a .csv\n",
    "test_pred.to_csv(\"Predictions/Prediction_SVR.csv\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}