F0D-15-Lasso.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FLIP (00): Data Science \n",
    "**(Module 03: Linear Algebra)**\n",
    "\n",
    "---\n",
    "- Materials in this module include resources collected from various open-source online repositories.\n",
    "- You are free to use,but NOT allowed to change and distribute this package.\n",
    "\n",
    "Prepared by and for \n",
    "**Student Members** |\n",
    "2006-2023 [TULIP Lab](http://www.tulip.org.au), Australia\n",
    "\n",
    "---\n",
    "\n",
    "## Session 22 Lasso"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## General Description\n",
    "This notebook will be divided into 3 parts, each part will be corresponded to one of the three question given in the this assigment, and each part will also be broke into several subsection in order to clarify my code and thoughts. \n",
    "\n",
    "Now, we will walk through the whole assignment in the following order:\n",
    "* Implementation of LassoRegressor\n",
    "* Tune & Train the LassoRegressor\n",
    "* Feature Selection and XRegressor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part I. Implementation of LassoRegressor\n",
    "\n",
    "The LassoRegressor is a type of regressor that deal with the regression problem by solving the cost function:\n",
    "$$\\arg\\min_{\\boldsymbol{w}\\in\\mathbb{R}^{D},b}\\frac{1}{2N}\\Vert\\boldsymbol{y}-(\\boldsymbol{X}\\boldsymbol{w}-b\\boldsymbol{1})\\Vert^{2}+\\alpha\\Vert\\boldsymbol{w}\\Vert_{1},$$\n",
    " \n",
    "Where $\\alpha$ is a constant that will determine how much the **1-norm** term will influence the cost function, and the term $N$ is the number of samples that used to train the regressor.\n",
    "\n",
    "To lower the result of the cost function above, one frequently used method is **gradient descent**:\n",
    "$$x^{(t+1)}=x^{(t)}-\\eta g$$\n",
    "Where t is the t-th iteration, $\\eta$ is learning rate, or stepsize, and $g$ is the gradient. By solving this equation iteratively, we can reach a **local optimum point** with a proper **stepsize**.\n",
    "\n",
    "\n",
    "Since we need to compute gradient, we have to apply **partial differentiation** to the cost function. However, one may found that is that the **1-norm** term might be **non-differentiable** at some points. This phenomenon will cause a big problem when we are training the regressor because that we want to use gradient descent method to train, but **we cannot obtain a gradient from a non-differentiable point.**\n",
    "\n",
    "To cope with that annoying non-differentialble problem, I choose **subgradient method** as a method to tackle the problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subgradient Method\n",
    "The idea of subgradient method is simple :\n",
    "* When the point is differentialble, we just simply compute the gradient\n",
    "* If it's not differentiable, we randomly assign a value $[-1, 1]$ to it, as a **subgradient** and use it to do the descent task.\n",
    "\n",
    "However, a **negative subgradient may to point to a descent direction.** Some constraint have to be imposed on it to make sure that our algorithm will descent.\n",
    "\n",
    "Thus, I impose the following constraint to the subgradient:\n",
    "$$\\eta=\\frac{f(x_k) - f^{lev}_{k}}{\\Vert s_k\\Vert ^2}$$\n",
    "\n",
    "Where $\\eta$ is a **modified** learning rate, $s_k$ is a subgradient, $f(x_k)$ is the cost function at iteration $k$, and $f^{lev}_{k}=\\min_{0\\leq i\\leq k}f(x_{i})$.\n",
    "\n",
    "This constraint can be simply interpreted as the description below:\n",
    "* If the regressor is converged to some point, the  $\\eta$ will be closed or equal to $0$.\n",
    "  - that is, when the regressor is closed to some local optimum points, it shouldn't take a big step size. This can prevent it from leaving the local optimum point since it has limited effect on the descent direction .\n",
    "* If the regressor is far from some points with low cost, the $\\eta$ will increase.\n",
    "  - that is, since you are not closed to any local optimum, you can take a big step.\n",
    "  \n",
    "By applying this constraint, we can make sure that the non-differential point won't drag the descent direction away from local optimum point, and thus the whole process will descent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation\n",
    "\n",
    "Now let's take a look at my implementation.\n",
    "The comments in the code below will explain the details of my implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pylab import *\n",
    "from copy import deepcopy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LassoRegressor(object):\n",
    "    \n",
    "    def __init__(self, eta=0.001, alpha=0.001, n_iter=40, random_state=1, zero=1e-10):\n",
    "        \"\"\" Initialize LassoRegressor.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        eta : float\n",
    "              Learing rate.\n",
    "        alpha : float\n",
    "              The coefficient of the 1-norm term in the cost function\n",
    "        n_iter : int\n",
    "              Number of iterations.\n",
    "        random_state : int\n",
    "              Random seed.\n",
    "        zero : float\n",
    "              Value that smaller than zero will be considered as zero.\n",
    "        \"\"\"\n",
    "        self.eta = eta\n",
    "        self.alpha = alpha\n",
    "        self.n_iter = n_iter\n",
    "        self.zero = zero\n",
    "        self.random_state = random_state\n",
    "        \n",
    "    def fit(self, X, y, stop_crit=1e+15):\n",
    "        \"\"\" Fit the model with given traning set.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        X : ndarray\n",
    "            Training data.\n",
    "        y : ndarray\n",
    "            Label of training data.\n",
    "        stop_crit : float\n",
    "            If the result of cost function exceeds this value,\n",
    "            the training process will be stop since it cannot\n",
    "            converge.\n",
    "        \n",
    "        Returns\n",
    "        -------\n",
    "        self :\n",
    "               Trained model.\n",
    "        \"\"\"\n",
    "        \n",
    "        # First, we initialize the weights with random value.\n",
    "        rgen = np.random.RandomState(self.random_state)\n",
    "        self.w_ = rgen.normal(loc=0.0, scale=0.01, size=1+X.shape[1])\n",
    "        self.cost_ = []\n",
    "        N = X.shape[0]\n",
    "        \n",
    "\n",
    "        # After initializing the weights, we start iterating to\n",
    "        # solve the cost function. The self.w_[0] indicates the\n",
    "        # intercept/bias of the function, so we don't need to \n",
    "        # do gradient descent on it.\n",
    "        for i in range(self.n_iter):\n",
    "            output = self.predict(X)\n",
    "            # this term is equal to :\n",
    "            # y - (Xw - b1)\n",
    "            errors = (y - output)\n",
    "            L1_norm = np.fabs(self.w_[1:]).sum()\n",
    "            \n",
    "            # calculate the \"SSE\" term.\n",
    "            # This is equal to the gradient of\n",
    "            # (1/2N)||y - (Wx - b1)||^2\n",
    "            SSE = self.eta * X.T.dot(errors) / N\n",
    "\n",
    "            # subgradient_descent_ will help us get\n",
    "            # the gradient of a||w||_1.\n",
    "            # You can see the detailed explanation in\n",
    "            # its function defintion.\n",
    "            L1_sub = self.subgradient_descent_()\n",
    "            \n",
    "            # update weights.\n",
    "            # the reason why 'SSE' we add SSE to weights instead\n",
    "            # of subtracting is that the partial derivative term\n",
    "            # will have a negative sign. So the updating rule:\n",
    "            # x(t+1) = x(t) - g(x) = x(t) - (-eta*(y-(Wx-b)))\n",
    "            #        = x(t) + (eta*(y-(Wx-b)))\n",
    "            # and the term 'L1_sub' do not have a negative sign,\n",
    "            # so it's still a subtrahend.\n",
    "            self.w_[1:] += SSE - L1_sub\n",
    "            \n",
    "            # update bias by simply adding the errors to i.\n",
    "            self.w_[0] += self.eta * errors.sum()\n",
    "            \n",
    "            # update the cost \n",
    "            cost = (errors**2).sum() / (2.0 * N) + self.alpha * L1_norm\n",
    "            self.cost_.append(cost)\n",
    "        \n",
    "            # stop if the errors exceed stop criteria.\n",
    "            if errors.sum() > stop_crit:\n",
    "                print(\"Error! Sum of errors exceeds stop criteria\")\n",
    "                break\n",
    "        return self\n",
    "    \n",
    "    def subgradient_descent_(self):\n",
    "        \"\"\" Computing the gradient for the differential point in 1-norm\n",
    "            as well as the subgradient for the non-differetiable point.\n",
    "        \"\"\"\n",
    "        # the derivative of 1-norm will be 1 or -1 according\n",
    "        # to its sign. \n",
    "        # (ex. -5 will have a derivative 1, while 0.5 will have 1.)\n",
    "        subgrad = np.sign(self.w_[1:])\n",
    "        n_eta = deepcopy(self.eta)\n",
    "        \n",
    "        # The 1-norm is non-differetiable when the value is equal\n",
    "        # to 0. Find this point by using np. operation\n",
    "        non_diff = (np.fabs(self.w_[1:]) <= self.zero)\n",
    "        can_diff = np.logical_not(non_diff)\n",
    "\n",
    "        if len(non_diff[non_diff == True]) > 0:\n",
    "            \n",
    "            subgrad[non_diff] = np.random.uniform(-1, 1)\n",
    "            \n",
    "            # apply the constraint to it.\n",
    "            # but we have to make sure we already have costs first.\n",
    "            if len(self.cost_) > 0:\n",
    "                \n",
    "                # get f(x_k), and find the minimum cost in the previous\n",
    "                # cost results.\n",
    "                cost = self.cost_[-1]\n",
    "                min_cost = min(self.cost_)\n",
    "                \n",
    "                # adjust the learning rate.\n",
    "                # if the current cost is equal to the minimum one,\n",
    "                # then the subgradient won't influence the descent direction;\n",
    "                # otherwise, it will change the direction and stepsize based\n",
    "                # on the difference between cost and min_cost.\n",
    "                n_eta = (cost - min_cost) / ((np.fabs(subgrad).sum())**2)\n",
    "                \n",
    "        # apply different learning rate to gradient and subgradient\n",
    "        subgrad[non_diff] *= n_eta\n",
    "        subgrad[can_diff] *= self.eta\n",
    "        \n",
    "        # multiply alpha\n",
    "        subgrad *= self.alpha\n",
    "\n",
    "        return subgrad\n",
    "    \n",
    "    def predict(self, X):\n",
    "        \"\"\" Predict the regression result of X by solving\n",
    "            y = w0 + w0x0 + w1x1 ... + wnxn\n",
    "        \"\"\"\n",
    "        return np.dot(X, self.w_[1:]) + self.w_[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we need to make sure that our LassoRegressor can have the **tendency of convergence**.\n",
    "\n",
    "We test it by using small $\\eta$ and $\\alpha$, and a very small iterations number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.linear_model import Lasso\n",
    "\n",
    "# No data file\n",
    "train_set = np.loadtxt('F:\\\\data\\\\x_train.csv', delimiter=',')\n",
    "y_train = np.loadtxt('F:\\\\data\\\\y_train.csv', delimiter=',')\n",
    "\n",
    "\n",
    "sc_x = StandardScaler()\n",
    "X_std = sc_x.fit_transform(train_set)\n",
    "#y_train_std = sc_y.fit_transform(y_train)\n",
    "\n",
    "llr = LassoRegressor(eta=0.001,alpha=0.001, n_iter=20)\n",
    "llr = llr.fit(X_std, y_train)\n",
    "\n",
    "\n",
    "plt.plot(range(1, llr.n_iter+1), llr.cost_)\n",
    "plt.ylabel('cost')\n",
    "plt.xlabel('Epoch')\n",
    "plt.tight_layout()\n",
    "plt.savefig('F:\\\\data\\\\fig-lasso-cost.png', dpi=300)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The figure above shows that our regressor does regress in only just 5 iterations!\n",
    "So why do I choose small $\\eta$ and $\\alpha$ ?\n",
    "\n",
    "Since we want to verify our implementation, we need to make sure it can work properly, that is, it can converge in finite iterations. However, some other factors may also let regressor be unable to converge:\n",
    "* Overshooting, which happens when the stepsize($\\eta$) is too big\n",
    "* The regressor cannot find a set of $w$ to solve the problem\n",
    "\n",
    "When we choose a $\\alpha$ too big, the coefficients of each features will be pressed to zero, and thus resulting large errors. So, to verify our implementation, we have to rule out these two situations by selecting small $\\eta$ and $\\alpha$, and thus we can identify the problem of our LassoRegressor more easily.\n",
    "\n",
    "\n",
    "Now, let's try it with more iterations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import r2_score\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "llr = LassoRegressor(eta=0.001,alpha=0.001, n_iter=500)\n",
    "llr = llr.fit(X_std, y_train)\n",
    "\n",
    "y_train_pred = llr.predict(X_std)\n",
    "\n",
    "print('MSE train: %.2f' % (mean_squared_error(y_train, y_train_pred)))\n",
    "print('R^2 train: %.2f' % (r2_score(y_train, y_train_pred)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "By examining the $R^2$ and $MSE$, we can see the regressor works fine.\n",
    "\n",
    "Knowing that our LassoRegressor is usable, we can begin our task 2: **Tune & Train the LassoRegressor**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part II. Tune & Train the LassoRegressor\n",
    "\n",
    "The purpose of this part is to tune and train our homemade LassoRegressor. Through the process of tuning and training, we can know how each hyper-parameter can affect the behavior and performance of our regressor.\n",
    "\n",
    "But there are a few things we have to do before we directly jump into the tuning process.\n",
    "\n",
    "### Inspecting data\n",
    "Inspecting data can help us:\n",
    "* Know if our data is learnable\n",
    "* Know the distribution of each data points/features\n",
    "\n",
    "Before starting the tuning/training tasks, we should first inspect, we should first ensure that our data is learnable. If the data is just a composition of random noise without any pattern, than it is unlikely to learn anything from it, and the regressor will never converge.\n",
    "\n",
    "Besides, even if we know that our data can converge(we know it, since we have made it converge in part I), we still need to inspect our data. By inspecting data, we can learn the distribution of data points, and can have a rough idea and direction about what we can try to tune our regressor.\n",
    "\n",
    "So, inpsecting data is an essential task for all who are trying to solve a problem through machine learning techniques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let us plot the data points of our training set, and check how each feature is related to the label."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.concatenate([X_std, y_train.reshape(-1, 1)], axis=1)\n",
    "df = pd.DataFrame(data)\n",
    "sns.set(style='whitegrid', context='notebook')\n",
    "\n",
    "# plot it separately for the ease of reading\n",
    "# notice that '13' the the label.\n",
    "sns.pairplot(df, x_vars=df.columns[0:int(X_std.shape[1]/2)], y_vars=df.columns[-1], size=2.5)\n",
    "sns.pairplot(df, x_vars=df.columns[int(X_std.shape[1]/2):X_std.shape[1]], y_vars=df.columns[-1], size=2.5)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('F:\\\\data\\\\fig-pairwise-dist.png', dpi=300)\n",
    "plt.show()\n",
    "sns.reset_orig()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By observing the figure above, one will notice that some features might be useless, such as feature 3 and 8, since they are not seen to be related to the labels, while features like 5 and 12 have some obvious patterns.\n",
    "\n",
    "And since most of the features have some patterns that related to our labels, I beleive that it is possible for us to learn things from this data set. Now, let us do some testing to support my assumption."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mse_lr(model, X, y):\n",
    "    return ((model.predict(X) - y)**2).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.learning_curve import learning_curve\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "test_set = np.loadtxt('F:\\\\data\\\\X_test.csv', delimiter=',')\n",
    "y_test = np.loadtxt('F:\\\\data\\\\y_test.csv', delimiter=',')\n",
    "\n",
    "num_data = X_std.shape[0]\n",
    "lr = LinearRegression()\n",
    "\n",
    "train_sizes, train_scores, test_scores = learning_curve(estimator=lr, X=X_std, y=y_train, scoring=mse_lr)\n",
    "\n",
    "train_mean = np.mean(train_scores, axis=1)\n",
    "test_mean = np.mean(test_scores, axis=1)\n",
    "\n",
    "plt.plot(train_sizes, train_mean, color='green', marker='o', markersize=5, label='training error')\n",
    "plt.plot(train_sizes, test_mean, color='blue', marker='s', markersize=5, label='testing error')\n",
    "\n",
    "plt.hlines(y=1, xmin=0, xmax=num_data, color='red', linewidth=2, linestyle='--')\n",
    "plt.title('Learning curve of our data')\n",
    "plt.grid()\n",
    "plt.xlabel('#traning samples')\n",
    "plt.ylabel('MSE')\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig('F:\\\\data\\\\fig-learning-curve.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the figure above shows, we can surely learn something with this data set. But, one thing we need to know is that if we just use this data to train our linear regressor without processing our data, we might be unable to reach a very low $MSE$, as you can see that the training and testing curve is converged at around 22-24 $MSE$. We will leave this problem to part III.\n",
    "\n",
    "As for now, we can start tuning, since we've confirmed that the data is learnable.\n",
    "\n",
    "### Tuning & Training\n",
    "Our tasks in this section is to :\n",
    "* tune $\\eta$ and $\\alpha$\n",
    "* tune $\\eta$ and $\\alpha$ to **make $w$ sparse**\n",
    "\n",
    "First, I will explain how I find $\\eta$ and $\\alpha$ to get a fine prediction result, and then I will show how I manage to tune an $\\alpha$ such that the resulting weights $w$ will be sparse.\n",
    "\n",
    "Now, let's do the first task. We have two hyper-parameters to tune: $\\eta$ and $\\alpha$. Before tuning alpha, we should first fine a proper $\\eta$, and then find $\\alpha$ by using fixed $\\eta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "etas = [0.00001, 0.0001, 0.001, 0.01, 0.1]\n",
    "max_iter = 30\n",
    "\n",
    "for eta in etas:\n",
    "    print('---- testing eta %.5f ----' % (eta))\n",
    "    llr = LassoRegressor(eta=eta, n_iter=max_iter)\n",
    "    llr = llr.fit(X_std, y_train)\n",
    "    \n",
    "    print(llr.cost_[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, as the printed results show, 0.001 might be the best choice for us, since values below it might descent too slow, and values greater than it will cause overshooting, so the sum of errors will eventually exceed stop criteria.\n",
    "\n",
    "Choosing $\\eta=0.001$, we can start seaching for a proper $\\alpha$ now. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cross_validation import train_test_split\n",
    "\n",
    "\"\"\"    \n",
    "    We have learned that eta = 0.001 is our best choice,\n",
    "    now let us find an appropriat alpha value\n",
    "\"\"\"\n",
    "\n",
    "a = 0.001\n",
    "eta = 0.001\n",
    "alphas = [0.001, 0.01, 0.1, 1]\n",
    "max_iter = 3000\n",
    "\n",
    "candidates = []\n",
    "r2_lists = []\n",
    "mse_lists = []\n",
    "a_lists = []\n",
    "\n",
    "for a in alphas:\n",
    "    print(\"running with init alpha : %.3f\" % (a))\n",
    "\n",
    "    a_list = []\n",
    "    r2_list = []\n",
    "    mse_list = []\n",
    "    inc = a\n",
    "    for i in range(0, 10):\n",
    "        # Examine the result of choosen alpha through training-validating process.\n",
    "        X_tune_train, X_val, y_tune_train, y_val = train_test_split(\n",
    "        X_std, y_train, test_size=0.4, random_state=1)\n",
    "\n",
    "        llr = LassoRegressor(eta=eta, alpha=a, n_iter=max_iter)\n",
    "        llr = llr.fit(X_tune_train, y_tune_train)\n",
    "        \n",
    "        y_val_pred = llr.predict(X_val)\n",
    "        r2  = r2_score(y_val, y_val_pred)        \n",
    "        mse = mean_squared_error(y_val, y_val_pred)\n",
    "        \n",
    "        a_list.append(a)\n",
    "        r2_list.append(r2)\n",
    "        mse_list.append(mse)\n",
    "    \n",
    "        a += inc\n",
    "    # collect testing result\n",
    "    a_lists.append(a_list)\n",
    "    r2_lists.append(r2_list)\n",
    "    mse_lists.append(mse_list)\n",
    "    \n",
    "print('Finished searching alpha')\n",
    "\n",
    "\n",
    "\n",
    "a = 0.001\n",
    "# now, we plot the resulting R^2 value of a given set of alpha values.\n",
    "# each line in the figure indicates the resulting value of R^2 from 'a' to 'a * 9'\n",
    "# (ex. if a = 0.001, then the line of 0.001 will be the R^2 value from 0.001 to 0.009)\n",
    "for l in r2_lists:\n",
    "    a_alpha = \"%.3f\" % a\n",
    "    plt.plot(range(0,10), l, label=a_alpha)\n",
    "    a *= 10\n",
    "\n",
    "plt.title('resulting R^2 of a given set of alpha')\n",
    "plt.ylabel('R2')\n",
    "plt.xlabel('alpha(x)')\n",
    "plt.legend(loc='best')\n",
    "plt.savefig('F:\\\\data\\\\fig-alpha-plot-r2', dpi=300)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# same as above, but this time, we plot the MSE value.\n",
    "a = 0.001\n",
    "for l in mse_lists:\n",
    "    a_alpha = \"%.3f\" % a\n",
    "    plt.plot(range(0,10), l, label=a_alpha)\n",
    "    a *= 10\n",
    "    \n",
    "plt.title('resulting MSE of a given set of alpha')\n",
    "plt.ylabel('MSE')\n",
    "plt.xlabel('alpha(x)')\n",
    "plt.legend(loc='best')\n",
    "plt.savefig('F:\\\\data\\\\fig-alpha-plot-mse', dpi=300)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the figure shows above, if we want to have a good $R^2$ or $MSE$ result, we might need to pick a small $\\alpha$. We will explain this later. Now let's compare our LassoRegressor with sklearn LinearRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "X_test = np.loadtxt('F:\\\\data\\\\x_test.csv', delimiter=',')\n",
    "y_test = np.loadtxt('F:\\\\data\\\\y_test.csv', delimiter=',')\n",
    "\n",
    "X_test = sc_x.transform(X_test)\n",
    "\n",
    "olr = LinearRegression()\n",
    "olr = olr.fit(X_std, y_train)\n",
    "\n",
    "llr = LassoRegressor(alpha=0.001, n_iter=max_iter)\n",
    "llr = llr.fit(X_std, y_train)\n",
    "\n",
    "y_train_pred_o = olr.predict(X_std)\n",
    "y_test_pred_o = olr.predict(X_test)\n",
    "\n",
    "y_train_pred = llr.predict(X_std)\n",
    "y_test_pred = llr.predict(X_test)\n",
    "\n",
    "\n",
    "print(\"Ordinary least square regressor\")\n",
    "print(\" R2 train: %.2f, test: %.2f\" %(\n",
    "        r2_score(y_train, y_train_pred_o), \n",
    "        r2_score(y_test, y_test_pred_o))\n",
    "     )\n",
    "print(\"MSE train: %.2f, test: %.2f\" %(\n",
    "        mean_squared_error(y_train, y_train_pred_o), \n",
    "        mean_squared_error(y_test, y_test_pred_o))   \n",
    "    )\n",
    "\n",
    "print(\"\\nMy LassoRegressor\")\n",
    "print(\" R2 train: %.2f, test: %.2f\" %(\n",
    "        r2_score(y_train, y_train_pred), \n",
    "        r2_score(y_test, y_test_pred))\n",
    "     )\n",
    "print(\"MSE train: %.2f, test: %.2f\" %(\n",
    "        mean_squared_error(y_train, y_train_pred), \n",
    "        mean_squared_error(y_test, y_test_pred))    \n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's good to know that our LassoRegressor works very well, and resulting in a little bit higher $R^2$ and $MSE$ result in test set. But the werid thing is: our regressor does not do any better than sklearn in the training phase. \n",
    "\n",
    "So why did it can surpass the scikit-Linear regressor(ordinary least square regressor)? One possible reason is that the testing data is too small, and the testing data may just happen to match the pattern our regressor has learned. Once the testing data become bigger and bigger, the scikit one may eventually surpass our regressor.\n",
    "\n",
    "\n",
    "### Tuning $\\alpha$ to make $w$ sparse\n",
    "Now, we have to tune the $\\alpha$ to make $w$ sparse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def plot_compress(coef, max_alpha, title='LassoRegressor'):\n",
    "    \"\"\" A helper function that help us plot our coefficients(weights)\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    coef : ndarray\n",
    "           List of coefficients(weights) with respects to some alpha value.\n",
    "    max_alpha : int\n",
    "           Maximum alpha value.\n",
    "    title : str\n",
    "           The title of resulting figure.\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    plt.hlines(y=0, xmin=0, xmax=max_alpha, color='red', linewidth = 2, linestyle = '--')\n",
    "\n",
    "    for i in range(coef.shape[1]):\n",
    "        plt.plot(range(max_alpha),coef[:,i])\n",
    "    \n",
    "    plt.title(title)\n",
    "    plt.ylabel('Coefficients')\n",
    "    plt.xlabel('Alpha')\n",
    "    plt.tight_layout()\n",
    "    plt.savefig('F:\\\\data\\\\fig-%s-decay.png' % title, dpi=300)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_alpha = 10\n",
    "\n",
    "coef_lr  = np.zeros((max_alpha, X_std.shape[1]))\n",
    "coef_llr = np.zeros((max_alpha, X_std.shape[1]))\n",
    "\n",
    "for a in range(max_alpha):\n",
    "    lr = LinearRegression()\n",
    "    lr.fit(X_std, y_train)\n",
    "    \n",
    "    llr = LassoRegressor(alpha=a, n_iter=max_iter)\n",
    "    llr.fit(X_std, y_train)\n",
    "    \n",
    "    coef_lr[a,:] = lr.coef_.reshape(1, -1)\n",
    "    coef_llr[a,:] = llr.w_[1:].reshape(1, -1)\n",
    "\n",
    "# plt.savefig('F:\\\\data\\\\fig-Scikit-LinearRegressor-decay.png' % title, dpi=300)\n",
    "plot_compress(coef_lr, max_alpha, title='Scikit-LinearRegressor')\n",
    "plot_compress(coef_llr, max_alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ordinary least square regressor will just simply compute:\n",
    "$$min_{w}\\Vert Xw - y\\Vert_{2}^{2}$$\n",
    "\n",
    "So it's not surprise to have a flat $w$.\n",
    "\n",
    "However, the $\\alpha$ term will work as a compressor in **LassoRegressor** : The larger the $\\alpha$, the smaller the coefficients, and **the less useful a feature is, the faster it will be compressed to zero**. \n",
    "\n",
    "And in our LassoRegressor, the $w$ will be become sparse will $\\alpha$ is closed to 1, and will have only 3 featues left when $\\alpha=2$.\n",
    "\n",
    "This figure can also explain when $R^2$ and $MSE$ escalated quick when we picking a $\\alpha$ larger than 1. When $\\alpha$ reaches 1, many of the weights have been pressed to zero, that is, some possibly useful features might be eliminated from the regressor. and when $\\alpha$ gets to 6, all the weights have been pressed to zero, which means that the regressor make prediction based on **only on bias, since other features have been totally eliminated**.\n",
    "\n",
    "So if you check the figure of $MSE$ and $R^2$, you will find that $R^2$ reaches 0.0 and the curve of $MSE$ is flatten after setting $\\alpha=6$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusion\n",
    "\n",
    "In this section, we have found that:\n",
    "* $\\eta=0.001$ with $\\alpha=0.001$ can make our regressor work well\n",
    "* when $\\alpha$ reaches 1, it will start becoming sparse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Part III. Feature Selecetion and XRegressor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quick summary \n",
    "* Best $\\overline{R}^2$ : $0.79$ (using 9 features), with $R^2$ : $0.86$ and $MSE$ : $7.96$\n",
    "\n",
    "In this section, we have to apply what we have learned during the past few weeks to make $\\overline{R}^2$ score as higher as possible. To complete this task, we have to develop our **XRegressor** to fulfill this task. I will discuss this in the following order:\n",
    "* Design of XRegressor\n",
    "* Feature selection\n",
    "* Removing outlier\n",
    "* Generating features\n",
    "* Final result\n",
    "\n",
    "Now, let us begin!\n",
    "### Design of XRegressor\n",
    "My XRegressor is still a lasso regressor. Or, more precisely, it's a lasso based regressor, with **Coordinate Descent** method and a pipepine designed for **data processing**. The core idea of  **Coordinate Descent** is simple:\n",
    "We perform line search procedure along **one coordinate direction at a time**, instead of computing the gradient of all weights at once. In each iteration, we will minimize the cost along one direction at a time. One benefit of this method is that since we only care about **one weight** at a time, so no partial differentiation needs to be calculated, and thus we don't need to worry about the subgradient problems that we encountered in part I. The update rules of my XRegressor is as follow:\n",
    "$$w_{j}=\\frac{S(\\sum^{N}_{i=1}x_{ij}(y_{i}-\\tilde{y_{i}}), \\alpha)}{\\sum_{i=1}^{N}x_{ij}^2}, \\text{where}$$\n",
    "$$S(w,\\alpha)=\\begin{cases}\n",
    "sign(w)(|w| - \\alpha) & \\text{if w - $\\alpha$} > 0\\\\ \n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}$$\n",
    "\n",
    "$w_j$ is the $j\\text{th}$ weights, $N$ is the number of training samples, $x_{ij}$ is the $j\\text{th}$ feature of $i\\text{th}$ sample, and $S$ is called  **soft thresholding**.\n",
    "\n",
    "More details will be given in the comments of **XRegressor** below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import RANSACRegressor\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "class XRegressor(object):\n",
    "    \"\"\"\n",
    "        A lasso regressor that uses coordinate descent as a descent method\n",
    "    \"\"\"\n",
    "    def __init__(self, alpha=0.001, n_iter=40, random_state=1, zero=1e-10):\n",
    "        \"\"\" Initilialize XRegressor.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        alpha : float\n",
    "                The coefficient of the 1-norm term in the cost function\n",
    "        n_iter : int\n",
    "                Number of iterations.\n",
    "        random_state : int\n",
    "                Randome seed.\n",
    "        zero : float\n",
    "                Value that smaller than zero will be considered as zero.              \n",
    "        \"\"\"\n",
    "        self.eta = eta\n",
    "        self.alpha = alpha\n",
    "        self.n_iter = n_iter\n",
    "        self.zero = zero\n",
    "        self.random_state = random_state\n",
    "        self.poly = None\n",
    "        self.e_f_ = False\n",
    "        self.r_f_ = False\n",
    "        self.p_f_ = False\n",
    "        \n",
    "    def fit(self, X, y, n_e_iter=3000, n_alpha=2.0, n_feat=5, rf_thresh=90.0, degree=2,\n",
    "                        extract_feats=False, remove_outlier=False, augment=False, stop_crit=1e+15):\n",
    "        \"\"\" Fitting model, with lots of option that can be used to process data.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        X : ndarray\n",
    "            Training data.\n",
    "        y : ndarray\n",
    "            Label of data.\n",
    "        n_e_iter : int\n",
    "            Number of iterations for Lasso feature selector.\n",
    "        n_alpha : float\n",
    "            Alpha value for Lasso feature selector.\n",
    "        n_feat : int\n",
    "            Number of selected features that the Lasso feature selector should return.\n",
    "        rf_thresh : float\n",
    "            Residual threshold for RANSAC outlier remover.\n",
    "        degree : int\n",
    "            Degrees of Polyfeatures.\n",
    "        extract_feats : boolean\n",
    "            If XRegressor should extract feature or not. (by using Lasso feature selector)\n",
    "        remove_outlier : boolean\n",
    "            If XRegressor should remove outlier or not. (by using RANSAC)\n",
    "        augment : boolean\n",
    "            If XRegressor should augment data or not. (by using PolyFeatures)\n",
    "        stop_crit : float\n",
    "            If the result of cost function exceeds this value,\n",
    "            the training process will be stop since it cannot\n",
    "            converge.\n",
    "        \n",
    "        Returns\n",
    "        -------\n",
    "        self : XRegressor\n",
    "               Trained model.\n",
    "        \"\"\"\n",
    "        \n",
    "        # selecting features using LassoRegressor in partI\n",
    "        # the pre-defined alpha value is 2, which is based on \n",
    "        # our experiments of alpha tuning in part 2.\n",
    "        if extract_feats:\n",
    "            self.e_f_ = True\n",
    "            llr = LassoRegressor(alpha=n_alpha, n_iter=n_e_iter)\n",
    "            llr = llr.fit(X, y)\n",
    "            # argsort will return an ascending list,\n",
    "            # use [::-1] to reverse it to get the first\n",
    "            # n_feat features.\n",
    "            feat = np.argsort(np.fabs(llr.w_[1:]))[::-1]\n",
    "            self.selected = feat[0:n_feat]\n",
    "        \n",
    "        if remove_outlier:\n",
    "            self.r_f_= True\n",
    "            \n",
    "        if augment:\n",
    "            self.p_f_ = True\n",
    "            self.p_degree = degree\n",
    "        \n",
    "        # process data.\n",
    "        X, y = self.process_data_(X, y, rf_thresh, training=True)\n",
    "        \n",
    "        rgen = np.random.RandomState(self.random_state)\n",
    "        self.w_ = rgen.normal(loc=0.0, scale=0.01, size=1+X.shape[1])\n",
    "        self.cost_ = []\n",
    "        N = X.shape[0]\n",
    "        # compute the term sum(X_ij)**2\n",
    "        norm_X = (X**2).sum(axis=0)\n",
    "        \n",
    "        # same as LassoRegressor, self.w_[0] is still the bias/intercept.\n",
    "        # in coordinate descent method, we need to perform the descent\n",
    "        # \"one weight at a time\", that is why we have 2 for-loop here.\n",
    "        # (the whole process is similar to the one used in sklearn:enet_coordinate_descent)\n",
    "        for i in range(self.n_iter):\n",
    "            for j in range(1, len(self.w_)):\n",
    "                # predict the result by fixing all value,\n",
    "                # but remove the term xj_wj from prediction.\n",
    "                r_resp = self.predict(X, ignore=j)\n",
    "                \n",
    "                # compute partial residual, that is, the cost \n",
    "                # along the coordinate.\n",
    "                pr = y - r_resp\n",
    "                    \n",
    "                # Trying to minimize the cost along coordinate j.\n",
    "                # the j-1 term is the weight of featue 0 is at index 1 of w_[],\n",
    "                # so subtract 1 to get correct index.\n",
    "                r_pred = X[:, j-1].dot(pr).sum()\n",
    "                \n",
    "                # perform thresholding to remove weights that is too low\n",
    "                self.w_[j] = self.soft_threshold_(r_pred, N) / norm_X[j-1]\n",
    "            output = self.predict(X)\n",
    "            errors = y - output\n",
    "\n",
    "            # compute the cost, same as LassoRegressor.\n",
    "            L1_norm = np.abs(self.w_[1:]).sum()\n",
    "            self.w_[0] += self.eta * errors.sum()\n",
    "            cost = (errors**2).sum() / (2.0 * N) + self.alpha * L1_norm\n",
    "            self.cost_.append(cost)\n",
    "            \n",
    "            \n",
    "            if errors.sum() >= stop_crit:\n",
    "                print(\"Error! Sum of errors exceeds stop criteria\")\n",
    "                break\n",
    "        \n",
    "        return self\n",
    "    \n",
    "    def soft_threshold_(self, r_pred, N):\n",
    "        \"\"\" Sort threshold to supress coefficient that is to small to zero.\n",
    "        \n",
    "        Parameters:\n",
    "        -----------\n",
    "        r_pred : ndarray\n",
    "                 A \"predicted\" value of new coefficient.\n",
    "        N : int\n",
    "                 Number of samples\n",
    "                     \n",
    "        Returns\n",
    "        -------\n",
    "            updated coefficient.\n",
    "        \"\"\"\n",
    "        if np.fabs(r_pred) < (self.alpha * N):\n",
    "            return 0\n",
    "        return np.sign(r_pred) * (np.fabs(r_pred) - (self.alpha * N))\n",
    "\n",
    "    def process_data_(self, X, y=None, thresh=90.0, training=False):\n",
    "        \"\"\" Pre-process the data to make it meet the format of regressor\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        X : ndarray\n",
    "            Data to be processed.\n",
    "        y : ndarray\n",
    "            Labels of data.\n",
    "        thresh : float\n",
    "            Residual threshold for RANSAC outlier remover.\n",
    "        training : boolean\n",
    "            In training phase of not.\n",
    "        \n",
    "        Returns\n",
    "        -------\n",
    "        X : ndarray\n",
    "            Processed training data.\n",
    "        y : ndarray\n",
    "            Processed training labels.\n",
    "        \"\"\"\n",
    "        \n",
    "        # process data through a pipeline according optional flags.\n",
    "        if self.e_f_:\n",
    "            X = self.extract_feats_(X)\n",
    "        \n",
    "        # we can only remove outlier in training set.\n",
    "        # It makes zero sense to remove data from testing set.\n",
    "        if self.r_f_ and training:\n",
    "            X, y = self.ransac_feats_(X, y, thresh)\n",
    "        \n",
    "        # augment data using PolynomailFeatures.\n",
    "        # we only need to init one instance for augmenting data\n",
    "        if self.p_f_:\n",
    "            if training:\n",
    "                self.poly = PolynomialFeatures(degree=self.p_degree)\n",
    "                self.poly.fit(X)\n",
    "            # we only take column 1:\n",
    "            # since that the column 0 is always equal to 1,\n",
    "            # which will be useless for predicting.\n",
    "            X = self.data_aug_(X)[:, 1:]\n",
    "        return X, y\n",
    "    \n",
    "    def extract_feats_(self, X):\n",
    "        \"\"\" Extraced feature with pre-computed result.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        X : ndarray\n",
    "            Data to be extraced.\n",
    "        \n",
    "        Returns:\n",
    "        X : ndarray\n",
    "            Selected features.\n",
    "        \"\"\"\n",
    "        return X[:, self.selected]\n",
    "    \n",
    "    def ransac_feats_(self, X, y, thresh=90.0):\n",
    "        \"\"\" Remove outlier with a given residual threshold.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        X : ndarray\n",
    "            Data to be processed.\n",
    "        y : ndarray\n",
    "            Labels to be processed.\n",
    "        thresh : float\n",
    "            Residual threshold.\n",
    "        \n",
    "        Returns:\n",
    "        X_filt : ndarray\n",
    "                 Inlier.\n",
    "        y_filt : ndarray\n",
    "                 Labels of inlier.\n",
    "        \"\"\"\n",
    "        # use RANSAC to remove outlier.\n",
    "        ransac = RANSACRegressor(LinearRegression(), \n",
    "                max_trials=100, \n",
    "                min_samples=50, \n",
    "                residual_threshold=thresh, \n",
    "                random_state=0)\n",
    "        ransac = ransac.fit(X, y)\n",
    "                   \n",
    "        inlier = ransac.inlier_mask_\n",
    "        X_filt = X[inlier, :]\n",
    "        y_filt = y[inlier]\n",
    "        return X_filt, y_filt\n",
    "    \n",
    "    def data_aug_(self, X):\n",
    "        \"\"\" Using pre-fitted PolyFeatures to do data augmentation.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        X : ndarray\n",
    "            Data to be augmented.\n",
    "        \n",
    "        Returns:\n",
    "            Augmented data.\n",
    "        \"\"\"\n",
    "        return self.poly.transform(X)\n",
    "    \n",
    "    def predict(self, X, ignore=None):\n",
    "        # for x^0 to x^n \n",
    "        \"\"\" Predict the result of given input.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        X : ndarray\n",
    "            Data to be predicted.\n",
    "        ignore : int\n",
    "            The feature that will be ignored during the predicting process.\n",
    "        \n",
    "        Returns\n",
    "        -------\n",
    "            Result of prediction.\n",
    "        \n",
    "        \"\"\"\n",
    "        if ignore == None:\n",
    "            return np.dot(X, self.w_[1:]) + self.w_[0]\n",
    "        \n",
    "        ign_ = deepcopy(self.w_[1:])\n",
    "        ign_[ignore-1] = 0\n",
    "        \n",
    "        return np.dot(X, ign_) + self.w_[0]\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we need to verify our design.\n",
    "\n",
    "Let's see if it can converge:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_set = np.loadtxt('F:\\\\data\\\\x_train.csv', delimiter=',')\n",
    "y_train = np.loadtxt('F:\\\\data\\\\y_train.csv', delimiter=',')\n",
    "\n",
    "sc_x = StandardScaler()\n",
    "X_std = sc_x.fit_transform(train_set)\n",
    "\n",
    "xr = XRegressor()\n",
    "xr = xr.fit(X_std, y_train, extract_feats=False)\n",
    "\n",
    "plt.plot(range(1, xr.n_iter+1), xr.cost_)\n",
    "plt.ylabel('cost')\n",
    "plt.xlabel('Epoch')\n",
    "plt.tight_layout()\n",
    "plt.savefig('F:\\\\data\\\\fig-lasso-cd-cost.png', dpi=300)\n",
    "plt.show()\n",
    "\n",
    "y_train_pred = xr.predict(X_std)\n",
    "print('MSE train: %.2f' % (mean_squared_error(y_train, y_train_pred)))\n",
    "print('R^2 train: %.2f' % (r2_score(y_train, y_train_pred)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result shows **XRegressor** does converge! And it also converges faster than the **subgradient descent** method! It can reach a very high $R^2$ in just a few iterations, but why? One possible reason is that we will always use \"fresh information\" to update the weight of one coordinate. For example, after completing updating the weight of $w_j$, $w_{j+1}$ will utilize this new $w_{j}$ to do the descent, and thus leads to a faster convergence rates.\n",
    "\n",
    "(referecnce: On the Finite Time Convergence of Cyclic Coordinate Descent Methods)\n",
    "\n",
    "Now, we also need to check if the coordinate descent Lasso can perform the feature selection task:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_alpha = 10\n",
    "coef = np.zeros((max_alpha, X_std.shape[1]))\n",
    "\n",
    "for a in range(max_alpha):\n",
    "    xr = XRegressor(alpha=a, n_iter=max_iter)\n",
    "    xr.fit(X_std, y_train)\n",
    "    coef[a,:] = xr.w_[1:].reshape(1, -1)\n",
    "    \n",
    "plot_compress(coef, max_alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **XRegressor** can successfully select features, according to the resulting figure.\n",
    "\n",
    "Now let's see how much $\\overline{R}^2$ we can get from testing set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def r2_bar_score(y, y_pred, N, P):\n",
    "    \"\"\" Compute adjusted R^2\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : ndarray\n",
    "        Labels.\n",
    "    y_pred : ndarray\n",
    "        Predictions.\n",
    "    N : int\n",
    "        Number of predictions.\n",
    "    P : int\n",
    "        Number of selected features.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "        adjusted R^2.\n",
    "    \"\"\"\n",
    "    return 1 - ((((y - y_pred)**2).sum()) / (N - P - 1)) / ((((y - np.mean(y))**2).sum()) / (N - 1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xr = XRegressor(alpha=0.001, n_iter=max_iter)\n",
    "xr = xr.fit(X_std, y_train)\n",
    "y_test_pred = xr.predict(X_test)\n",
    "\n",
    "print(\"R2_bar_score: %.2f\" %(r2_bar_score(y_test, y_test_pred, X_test.shape[0], xr.w_.shape[0]-1)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is extremely bad! So we have to take some action to improve our performance. \n",
    "\n",
    "In the following section, you will see how I deal with all those data, and improve the $\\overline{R}^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature selection\n",
    "\n",
    "Our **LassoRegressor** will act as the feature selector of our **XRegressor**. We choose the same $\\alpha$ value, and plot the selected features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "llr = LassoRegressor(alpha=2, n_iter=3000)\n",
    "llr = llr.fit(X_std, y_train)\n",
    "feat = np.argsort(np.fabs(llr.w_[1:]))[::-1]\n",
    "\n",
    "for i in range(len(llr.w_[1:])):\n",
    "    print(\"Feature #%d : %.5f\" %(feat[i], llr.w_[feat[i]+1]))\n",
    "\n",
    "n_sel = 6\n",
    "selected = feat[0:n_sel]\n",
    "selected_data = np.concatenate([X_std, y_train.reshape(-1, 1)], axis=1)\n",
    "df = pd.DataFrame(selected_data)\n",
    "sns.set(style='whitegrid', context='notebook')\n",
    "sns.pairplot(df, x_vars=df.columns[selected], y_vars=df.columns[-1], size=2.5)\n",
    "plt.tight_layout()\n",
    "plt.savefig('F:\\\\data\\\\fig-feat-best.png', dpi=300)\n",
    "plt.show()\n",
    "sns.reset_orig()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reason why we plot out selected features is that:\n",
    "* we want to know what features is selected\n",
    "* we can examine if the feature is truly meaningful\n",
    "\n",
    "As we have mentioned in previous section, it's important for us to know what kind of data we are dealing with, and plotting out these figures do help us to understand our data. Although we can know if a feature is important or not through the process of tuning and traning, sometimes we have to rely on our own judgement, and decide if we truly need these selected features or not.\n",
    "\n",
    "And the result do meet what we have expected in part II : feature 5 and 12 are both great for prediction, and feature 3 and 8 are useless as expected, since those features are not seen to related to the labels.\n",
    "\n",
    "You can see there is actually a large gap between $3\\text{th}$ and $4\\text{th}$ features, so here we can make an assumptions that **we do not really need more than 3 features**, so let's do some simple experiment to verify this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_train_val(y_train_pred, y_train, y_valid_pred, y_valid, train_size, valid_size, P):\n",
    "    \"\"\" Print out training/validating result\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y_train_pred : ndarray\n",
    "                   Prediction of training data.\n",
    "    y_train : ndarray\n",
    "                   Label of training data.\n",
    "    y_valid_pred : ndarray\n",
    "                   Prediction of validation data.\n",
    "    y_valid : ndarray\n",
    "                   Label of validation data.\n",
    "    train_size : int\n",
    "                   Number of samples of training data.\n",
    "    valid_size : int\n",
    "                   Number of samples of validation data.\n",
    "    P : int\n",
    "                   Number of selected features.\n",
    "    \"\"\"\n",
    "    print(\"      Training R2_bar: %.2f, R2: %.2f, MSE: %.2f\" % (\n",
    "        r2_bar_score(y_train, y_train_pred, train_size, P),\n",
    "        r2_score(y_train, y_train_pred),\n",
    "        mean_squared_error(y_train, y_train_pred)\n",
    "    ))\n",
    "    \n",
    "    print(\"      Testing(Validating) R2_bar: %.2f, R2: %.2f, MSE: %.2f\" % (\n",
    "        r2_bar_score(y_valid, y_valid_pred, valid_size, P),\n",
    "        r2_score(y_valid, y_valid_pred),\n",
    "        mean_squared_error(y_valid, y_valid_pred)\n",
    "    ))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_set = np.loadtxt('F:\\\\data\\\\x_train.csv', delimiter=',')\n",
    "y_train = np.loadtxt('F:\\\\data\\\\y_train.csv', delimiter=',')\n",
    "X_std = sc_x.fit_transform(train_set)\n",
    "\n",
    "max_feats = 6\n",
    "mse_best = 100\n",
    "r2b_best = 0\n",
    "\n",
    "selected_list = []\n",
    "\n",
    "for i in range(1, max_feats+1):\n",
    "    sc_x = StandardScaler()\n",
    "\n",
    "    X_tune_train, X_valid, y_tune_train, y_valid = train_test_split(\n",
    "        X_std, y_train, test_size=0.5, random_state=0)\n",
    "    \n",
    "    xr = XRegressor(alpha=0.001, n_iter=max_iter)\n",
    "    xr = xr.fit(X_tune_train, y_tune_train, n_feat=i, extract_feats=True)\n",
    "\n",
    "    y_train_pred = xr.predict(X_tune_train[:, xr.selected])\n",
    "    y_valid_pred = xr.predict(X_valid[:, xr.selected])\n",
    "    print(\"=== %d features ===\" %(i))\n",
    "    print(\"Selected : %s\" %(xr.selected))\n",
    "    print_train_val(y_train_pred, y_tune_train, y_valid_pred, y_valid,\n",
    "                    X_tune_train.shape[0], X_valid.shape[0], xr.w_.shape[0]-1)\n",
    "    \n",
    "    selected_list.append(xr.selected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results indicates that using lesser than 2 features might not be a smart choice, while using 3 to 6 features may result in similar performance, which match our assumption.\n",
    "\n",
    "Since there exists only a very small difference between using 3 features and 6 features, it's wise to choose 3 features, since that:\n",
    "* the model will be simpler\n",
    "* those 3 are the most significant features\n",
    "\n",
    "Those 3 features together increase the peformance significantly, while other features contribute only a little, so we can regard these features as major factors of the regression prediction.\n",
    "\n",
    "Now, let's further inspect these features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "selected = selected_list[2]\n",
    "\n",
    "data = np.concatenate([X_std[:, selected], y_train.reshape(-1, 1)], axis=1)\n",
    "df = pd.DataFrame(data)\n",
    "sns.set(style='whitegrid', context='notebook')\n",
    "sns.pairplot(df, x_vars=df.columns[0:X_std[:,selected].shape[0]], y_vars=df.columns[-1], size=4)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('F:\\\\data\\\\fig-selected-3feats', dpi=300)\n",
    "plt.show()\n",
    "sns.reset_orig()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Closely examine these features, one may find that they do have some pattern. However, one may also notice that there exists many **outliers** that may influence the prediction result, these outlier may possibly drag our regressor away from a good regression line. Thus, removing outlier is a must if we want to improve the accuracy.\n",
    "\n",
    "### Removing outlier\n",
    "To remove outlier, we apply **RANSACRegressor** to selected features, mask out those outliers, and use only inliers to train our regressor.\n",
    "\n",
    "The idea of ransac is simple: we just randomly select some data points, and see if we can use some linear lines to fit those points. For all those points that are within the **residual threshold**, that is, is close enough to the line, will be marked as an inlier, and others will be treated as outliers. This whole process will iterate $T$ times(in my code, I set $T=100$), and the iteration that has the most inliers will be regarded as the **correct answer**, and use the result of that iteration to mask out outliers.\n",
    "\n",
    "Let's see what we will get after removing outliers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import RANSACRegressor\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "ransac = RANSACRegressor(LinearRegression(), \n",
    "            max_trials=100, \n",
    "            min_samples=50, \n",
    "            residual_threshold=90.0, \n",
    "            random_state=0)\n",
    "\n",
    "ransac.fit(X_std[:, selected], y_train)\n",
    "\n",
    "inlier_mask = ransac.inlier_mask_\n",
    "outlier_mask = np.logical_not(inlier_mask)\n",
    "for i in selected:\n",
    "    X_sel = X_std[:, [i]]\n",
    "    line_X = np.arange(-4, 5, 1)\n",
    "    #line_y_ransac = ransac.predict(line_X[:, np.newaxis])\n",
    "    plt.scatter(X_sel[inlier_mask], y_train[inlier_mask], \n",
    "                c='blue', marker='o', label='Inliers')\n",
    "    plt.scatter(X_sel[outlier_mask], y_train[outlier_mask],\n",
    "                c='lightgreen', marker='s', label='Outliers')\n",
    "    #plt.plot(line_X, line_y_ransac, color='red')\n",
    "    plt.xlabel('feature %d' %(i))\n",
    "    plt.ylabel('Label %d' %(i))\n",
    "    plt.legend(loc='upper left')\n",
    "    plt.tight_layout()\n",
    "    plt.show()  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the whole process is proceed in **3-dimension**, it's hardly for us to know which point will be consider outliers. But as you can see in the figures above, most of the points that seen to related to labels have been selected, so we can expect that the result may be improved, since some annoying outliers are gone.\n",
    "\n",
    "One thing that worth mentioning is the value of **residual threshold**. This value will largely affect the result of removing outliers. If you set it too high, then nothing will be filtered out; and if it's too low, you will only have few data points that met the RANSAC line, and since number of points is too small, it's unlikely to make our regressor learn a correct pattern.\n",
    "\n",
    "I choose the value 4.0 as my threshold based on some experiments. The 4.0 works good for RANSAC to keep all the meaningful data, and remove most of the outliers. So let's see if removing outliers really help to improve $R^2$, and see if we can get a nice $\\overline{R}^2$ scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xr = XRegressor(alpha=0.001, n_iter=max_iter)\n",
    "xr_nofilt = XRegressor(alpha=0.001, n_iter=max_iter)\n",
    "\n",
    "xr = xr.fit(X_std, y_train, n_feat=3, rf_thresh=90.0, extract_feats=True, remove_outlier=True)\n",
    "xr_nofilt = xr_nofilt.fit(X_std, y_train, n_feat=3, extract_feats=True)\n",
    "\n",
    "\n",
    "\n",
    "y_test_pred = xr.predict(X_test[:, xr.selected])\n",
    "\n",
    "print(\"XRegressor(removing outliers) Testing R2_bar: %.2f, R2: %.2f, MSE: %.2f\" % (\n",
    "    r2_bar_score(y_test, y_test_pred, X_test.shape[0], xr.w_.shape[0]-1),\n",
    "    r2_score(y_test, y_test_pred),\n",
    "    mean_squared_error(y_test, y_test_pred)\n",
    "    ))\n",
    "\n",
    "y_test_pred = xr_nofilt.predict(X_test[:, xr_nofilt.selected])\n",
    "\n",
    "print(\"XRegressor(without removing outliers) Testing R2_bar: %.2f, R2: %.2f, MSE: %.2f\" % (\n",
    "    r2_bar_score(y_test, y_test_pred, X_test.shape[0], xr.w_.shape[0]-1),\n",
    "    r2_score(y_test, y_test_pred),\n",
    "    mean_squared_error(y_test, y_test_pred)\n",
    "    ))\n",
    "\n",
    "selected = xr.selected\n",
    "X_filt, y_filt = xr.ransac_feats_(X_std[:, selected], y_train)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results indicates that by removing outliers, we can achieve a great performance boost! And not only do we have a better $\\overline{R}^2$, we also get a good $R^2$. \n",
    "\n",
    "After selecting features and removing outliers, now we have a good set of data points for us to traing our regressor. Knowing this data are good, we should futher \"exploit\" these data points, to see if we can push our performance to the limit!\n",
    "\n",
    "### Generating features\n",
    "Our selected features are good, and since it's good, it's reasonable for us to think that their polynomial features are useful. So, we use **PolynomialFeatures**(which we learned from Lab) to generate more new features. For example, if we have features $a$ and $b$, their degree-2 polynomial features will be $[1, a, b, a^2, ab, b^2]$. **The column 0 should be removed before we start training, since that it will always equal to 1**, which is meaningless in during the training phase, and will influence our prediction result.\n",
    "\n",
    "\n",
    "Since we want to have a high $\\overline{R}^2$ score, it's important to control the total number of features, since the more features we use, the more penalty we will get, and that is why I just use max degree = 3 in the code below. So let's see if generating new features can help us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "ds = [1, 2, 3]\n",
    "for d in ds:\n",
    "    poly = PolynomialFeatures(degree=d)\n",
    "    X = poly.fit_transform(X_filt)[:, 1:] # remove column 0, since col 0 will always equal to 1.\n",
    "    X_tune_train, X_valid, y_tune_train, y_valid = train_test_split(\n",
    "        X, y_filt, test_size=0.3, random_state=0)\n",
    "    \n",
    "    xr = XRegressor(alpha=0.001, n_iter=max_iter)\n",
    "    xr = xr.fit(X_tune_train, y_tune_train)\n",
    "    \n",
    "    y_train_pred = xr.predict(X_tune_train)\n",
    "    y_valid_pred = xr.predict(X_valid)\n",
    "\n",
    "    print(\" Train-Val with degree %d (%d features)\" %(d, X_tune_train.shape[1]))\n",
    "    print_train_val(y_train_pred, y_tune_train, y_valid_pred, y_valid,\n",
    "                    X_tune_train.shape[0], X_valid.shape[0], xr.w_.shape[0]-1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, we do get a super huge performance boost, and we even get $R^2$ score that close to $0.9$! And the $MSE$ values are also low, and break the barrier we have mentioned in partII(we can only have 22-24 $MSE$ in training phase). \n",
    "\n",
    "But now, we encounter a big issue: which degree we should choose? Degree 2 and 3 and result in a very closed performance, and degree 3 has even higher $MSE$ and $R^2$ performance. However, I finally decide to pick degree 2 based on the following reasons:\n",
    "* adding 10 features only result in a 0.19 $MSE$ increase, indicating that most poly-features are useless in degree 3\n",
    "* it doing good in both trainig and validating\n",
    "* 19 features is too much(we have 13 features originally)\n",
    "* the model complexity will be much smaller if we pick degree 2\n",
    "\n",
    "By using 9 features, we have reached a good $\\overline{R}^2$, with a much simpler model in contrast to using 19 features. Now we can finally see how much we can get from testing set!\n",
    "\n",
    "### Final result\n",
    "* iteraions : 3000\n",
    "* $\\alpha$ : 0.001\n",
    "* number of features : 9 (3 orignal, 6 polynomial)\n",
    "* residual threshold : 4.0 (for removing outlier)\n",
    "\n",
    "result : **$\\overline{R}^2 : 0.79$,  $R^2 : 0.86$**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xr = XRegressor(alpha=0.001, n_iter=max_iter)\n",
    "xr.fit(X_std, y_train, n_feat=3, degree=2, rf_thresh=90.0,\n",
    "       extract_feats=True, remove_outlier=True, augment=True)\n",
    "\n",
    "X_poly_test, _ = xr.process_data_(X_test)\n",
    "y_test_pred = xr.predict(X_poly_test)\n",
    "\n",
    "\n",
    "print(\"XRegressor (Final result) Testing R2_bar: %.2f, R2: %.2f, MSE: %.2f\" % (\n",
    "    r2_bar_score(y_test, y_test_pred, X_poly_test.shape[0], xr.w_.shape[0]-1),\n",
    "    r2_score(y_test, y_test_pred),\n",
    "    mean_squared_error(y_test, y_test_pred)\n",
    "    ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Conclusion\n",
    "In this homework, we need to use **everything** that we have learned from previous class and lab. Lots of inspection and observation also take a big part: we need to interpred data, and making our own judgement, since there does not exist a single answer for this assigmemt, **blindly tune the model will not really increase the result**. Both **coding** and **interpreting** data are important when we want to solve machine learning problem!\n",
    "\n",
    "And it's really exciting to see that we can improve the $\\overline{R}^2$ from 0.44 on testing set, to the final result 0.79."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}