{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework 6: Prediction on Housing Prices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from sklearn import linear_model as lm\n",
    "\n",
    "from IPython.display import display, Latex, Markdown"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -U okpy\n",
    "from client.api.notebook import Notebook\n",
    "ok = Notebook('hw6.ok')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kaggle \n",
    "\n",
    "This assignment is purposefully left nearly open-ended.  The Ames data in your possession comes from a larger data set.  Your goal is to provide a linear model that accurately predicts the prices of the held-out homes, measured by root mean square error.  That is, the score you will see on the Kaggle leaderboard is calculated as follows:\n",
    "\n",
    "$$score = \\sqrt{\\dfrac{\\sum_{\\text{houses in public test set}}(\\text{actual price for house} - \\text{predicted price for house})^2}{\\text{# of houses}}}$$\n",
    "\n",
    "Perfect prediction of house prices would have a score of 0, so you want your score to be as low as possible!\n",
    "\n",
    "**Kaggle Submission Site:** https://inclass.kaggle.com/c/ds100-2017-hw6  \n",
    "**Max number of submissions per day:** 2  \n",
    "**Max number of final submissions:** 1 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Ames data set consists of 2930 records taken from the Ames Assessor’s Office.  The data set has 23 nominal, 23 ordinal, 14 discrete, and 20 continuous variables (and 2 additional observation identifiers) --- 82 features in total.  An explanation of each variable can be found in the included `README.txt` file.  The information was used in computing assessed values for individual residential properties sold in Ames, Iowa from 2006 to 2010.  Since the data is publicly available, we have injected noise into all the sale prices to remove the temptation to do \"oracle learning.\"\n",
    "\n",
    "The data are split into training and test sets with 2000 and 930 observations, respectively.  The actual sale price is withheld from you in the test set.  In addition, the test data are further split into public and private test sets.  When you upload a test set prediction onto Kaggle for validation, the score you receive will be calculated using the public test set.  The private test set will be used in the final evaluation of this homework assignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "raw_data = pd.read_csv(\"ames_train.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example Data\n",
    "\n",
    "Throughout this assignment, we will use this reduced data set for examples. This is only for demonstration; in your final submission you'll want to use more features than just these."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "small_data = (\n",
    "    raw_data[[\"SalePrice\", \"Gr_Liv_Area\", \"Lot_Area\", \"Bedroom_AbvGr\"]]\n",
    "    .rename(columns = {\n",
    "        \"SalePrice\": \"price\",\n",
    "        \"Gr_Liv_Area\": \"sqft\",\n",
    "        \"Lot_Area\": \"lotsize\",\n",
    "        \"Bedroom_AbvGr\": \"bedrooms\"\n",
    "    })\n",
    ")\n",
    "\n",
    "small_data.iloc[1:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Grading\n",
    "Grading will be based on a number of set criteria, enumerated below:\n",
    "\n",
    "Task | Description\n",
    "--- | ---\n",
    "EDA | You create exploratory plots for at least 3 (basic) features to motivate your work.  The minimal 3 should cover each of the 3 variable types: categorical, discrete, continuous.\n",
    "Transformations | Your final model includes transformations of the data.\n",
    "Diagnostics | You have diagnostic checks with commentary for your model\n",
    "RMSE | Your model beats the RMSE threshold of **$30,000**.  This should be attainable with a well-thought-out model.\n",
    "Model | Your modeling pipeline is encapsulated in a pipeline object called `final_pipeline`.\n",
    "Written Questions | Your submission should include answers to the written questions at the bottom of this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Submission\n",
    "\n",
    "This assignment requires a Kaggle submission in addition to the usual okpy one.  To submit to Kaggle, you should create a `csv` file with 930 rows---one for each house in the test data---and 2 columns:\n",
    "\n",
    "* `PID` The house identification number\n",
    "* `SalePrice` Your estimate for the sale price of the house\n",
    "\n",
    "An example kaggle submission file has been included with this assignment.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Restrictions\n",
    "\n",
    "While we want you to be creative with your models, we want to make it fair to students who are seeing these techniques for the first time.  As such, you are only allowed to train linear models and their regularized forms (e.g. ridge and lasso).  This means no random forest, CART, neural nets, etc.  However, you are free to feature engineer to your heart's content.  Remember that domain knowledge is the third component of data science...\n",
    "\n",
    "That being said, you may want to explore the [sklearn API](http://scikit-learn.org/stable/modules/classes.html#module-sklearn.linear_model) for more information on Lasso, Ridge, and ElasticNet."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prizes\n",
    "The top 10 students, evaluated by their score in the private test set will: \n",
    "1. Have bragging rights \n",
    "2. Be invited to attend a lunch at the Faculty Club, hosted by Professor Yu."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EDA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make plots to explore the data. You may create as many plots as you wish and you will choose three of them to be graded.\n",
    "\n",
    "The 3 plots you submit should cover each of the three variable types (categorical, discrete, continuous).\n",
    "\n",
    "Insert this comment at the top of the three cells you want to submit for EDA:\n",
    "\n",
    "    # EDA_SUBMIT\n",
    "    \n",
    "We will use this tag to grade your 3 submitted plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# EDA_SUBMIT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# EDA_SUBMIT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# EDA_SUBMIT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transformations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have already encountered one of `sklearn`'s `Transformer` classes: the `DictVectorizer` from lab 10.  A transformer is an object that cleans, reduces, expands, or generates features.  \n",
    "\n",
    "A transformer's `fit` method, will learn parameters. For the `DictVectorizer`, the parameters are the allowed values for a categorical variable.\n",
    "\n",
    "The `transform` method takes the learned parameters and transforms any inputted new data.  For `DictVectorizer`, this means taking a vector of categorical values and transforming the data into a matrix where each row has at most one non-zero value (they may all be zero if this vector contains a category that was previously unseen).\n",
    "\n",
    "`fit_transform` simply learns from and transforms the input data all in one go\n",
    "\n",
    "Since we might want to perform different transformations on different columns, we've provided you with a `ColumnSelector` class.  You may want to use our code as a template for your own custom transformers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.base import BaseEstimator, TransformerMixin\n",
    "\n",
    "class ColumnSelector(BaseEstimator, TransformerMixin):\n",
    "    \"\"\"\n",
    "    Transformer that extracts a column of a data frame\n",
    "    \n",
    "    Example Usage\n",
    "    >> data = pd.DataFrame({'a': [1, 2, 3, 4],\n",
    "                            'b': [5, 6, 7, 8],\n",
    "                            'c': [9, 10, 11, 12]})\n",
    "    >> cs = ColumnSelector(cols=['a', 'b'])\n",
    "    >> data['a'] == cs.transform(data)\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    col : list of strings, required\n",
    "        The name(s) corresponding to the desired column of a DataFrame.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, cols):\n",
    "        self.cols = cols\n",
    "\n",
    "    def fit(self, X, y=None):\n",
    "        \"\"\"\n",
    "        Returns itself, nothing to be fit\n",
    "        \"\"\"\n",
    "        return self\n",
    "\n",
    "    def transform(self, X, y=None):\n",
    "        \"\"\"\n",
    "        Returns the desired column as a matrix\n",
    "        \"\"\"\n",
    "        return X.as_matrix(self.cols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have also seen another transformation in class: the polynomial transformation.  In practice, you would use `sklearn`'s nice `PolynomialFeatures`.  To give you experience implementing your own transformer class, write a bivariate (exactly 2 input features) `BiPolyTrans` transformer class that, given two features, $W$ and $Z$ of a matrix $X$, calculates all powers up to a given `degree`.  That is for every record (row) $x_i = \\begin{bmatrix} w_i & z_i \\end{bmatrix}$, \n",
    "$$\\phi_{degree}(x_i) = \\begin{bmatrix} 1 & w_i & z_i & w_iz_i & w_i^2z_i & w_iz_i^2 & \\dots & w_iz_i^{degree-1} & w_i^{degree} & z_i^{degree} \\end{bmatrix} $$\n",
    "\n",
    "If you are worried about efficiency, you may want to make use of Python's `itertools`.  Namely, `chain` and `combinations_with_replacement` should be helpful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "for_assignment_type": "solution"
   },
   "outputs": [],
   "source": [
    "from itertools import chain, combinations_with_replacement\n",
    "\n",
    "class BiPolyTrans(BaseEstimator, TransformerMixin):\n",
    "    \"\"\"\n",
    "    Transforms the data from a n x 2 matrix to a matrix with\n",
    "    polynomial features up to the specified degree.\n",
    "    \n",
    "    Example Usage\n",
    "    data = np.array([[1, 2], [3, 4]])\n",
    "    d3polytrans = BiPolyTrans(2)\n",
    "    d3polytrans.fit_transform(data) == np.array([[1, 1, 2, 1, 2, 4], [1, 3, 4, 9, 12,16]])\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    degree : integer, required\n",
    "        largest polynomial degree to calculate with the two features\n",
    "    \"\"\"\n",
    "    def __init__(self, degree):\n",
    "        self.degree = degree\n",
    "    \n",
    "    def fit(self, X, y=None):\n",
    "        \"\"\"\n",
    "        Calculates the number of input and output features\n",
    "        \"\"\"\n",
    "        self.n_input_features = 2\n",
    "        self.n_output_features = (self.degree + 1) * (self.degree + 2) // 2\n",
    "        return self\n",
    "    \n",
    "    def transform(self, X, y=None):\n",
    "        \"\"\"\n",
    "        Transforms the data into polynomial features\n",
    "        \n",
    "        Input\n",
    "        -----\n",
    "        X : an n x 2 matrix, required.\n",
    "            \n",
    "        Output\n",
    "        ------\n",
    "        A higher-dimensional matrix with polynomial features up to the specified degree\n",
    "        \"\"\"\n",
    "        n_records = X.shape[0]\n",
    "        output = np.empty((n_records, self.n_output_features), dtype=X.dtype)\n",
    "        \n",
    "        # Get all combinations up to specified degree\n",
    "        combs = chain.from_iterable(\n",
    "            combinations_with_replacement(range(self.n_input_features), i)\n",
    "            for i in range(self.degree + 1)\n",
    "        )\n",
    "        \n",
    "        for colnum, comb in enumerate(combs):\n",
    "            output[:, colnum] = X[:, comb].prod(axis = 1)\n",
    "        \n",
    "        return(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.grade('qtransform')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pipelines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.pipeline import Pipeline, FeatureUnion\n",
    "from sklearn.preprocessing import OneHotEncoder, PolynomialFeatures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, we will formalize our data cleaning, extraction, transformation, and training all into an abstraction called a Pipeline.  In a nutshell, a pipeline is the recipe for going from a clean but untransformed data set to a trained model.  For more information, see [sklearn's docs](http://scikit-learn.org/stable/modules/pipeline.html).  In the example below, we extract polynomial features from each home's square footage and then fit a linear model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex_pipeline1 = Pipeline([\n",
    "    ('selector', ColumnSelector(['sqft'])), \n",
    "    ('poly_feats', PolynomialFeatures(3, include_bias=False)),\n",
    "    ('lm', lm.LinearRegression(fit_intercept=False))\n",
    "])\n",
    "\n",
    "ex_pipeline1.fit(small_data, small_data[['price']])\n",
    "\n",
    "print(\"Training RMSE:\",\n",
    "      (np.mean((ex_pipeline1.predict(small_data) \n",
    "                - small_data[['price']])**2)**(.5)).item())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we've learned, training error definitely isn't everything! In addition to our training error, we want to be able to calculate validation error using cross validation. Luckily, sklearn makes this quite easy for us. Below, we calculate validation error for our initial pipeline, using 10-fold cross-validation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import make_scorer\n",
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "def score_func(y, y_pred, **kwargs): \n",
    "    return np.mean((y-y_pred)**2)**0.5\n",
    "scorer = make_scorer(score_func)\n",
    "\n",
    "cv_scores = cross_val_score(ex_pipeline1, small_data, small_data[['price']], cv=10, scoring=scorer)\n",
    "print(\"Validation RMSE:\", np.mean(cv_scores))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course we wouldn't want to just use one predictor given how rich our dataset is.  To append more features, we can use `FeatureUnion`, which combines several transformers into a MEGATRANSFORMER, which outputs a concatenation of the output of its constituents.  Note: `FeatureUnion` does NOT check if the transformations create linearly independent features.  In the example below, we combine the polynomial lift of square footage and lot size.\n",
    "\n",
    "<img src=\"pipeline.png\" style=\"height: 3in;\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex_pipeline2 = Pipeline([\n",
    "    ('union', FeatureUnion(n_jobs=1,  transformer_list=[\n",
    "        ('poly_lotsize', Pipeline([\n",
    "            ('selector', ColumnSelector(['lotsize'])), \n",
    "            ('poly_feats', PolynomialFeatures(3, include_bias=False))\n",
    "        ])),\n",
    "        \n",
    "        ('poly_sqft', Pipeline([\n",
    "            ('selector', ColumnSelector(['sqft'])), \n",
    "            ('poly_feats', PolynomialFeatures(3, include_bias=False)),\n",
    "        ]))\n",
    "    ])),\n",
    "        \n",
    "    ('lm', lm.LinearRegression(fit_intercept=False))\n",
    "])\n",
    "\n",
    "ex_pipeline2.fit(small_data, small_data[['price']])\n",
    "\n",
    "print(\"Training RMSE:\",\n",
    "      (np.mean((ex_pipeline2.predict(small_data) \n",
    "                - small_data[['price']])**2)**(.5)).item())\n",
    "\n",
    "cv_scores = cross_val_score(ex_pipeline2, small_data, small_data[['price']], cv=10, scoring=scorer)\n",
    "print(\"Validation RMSE:\", np.mean(cv_scores))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your final model should be presented as a data pipeline.  We should be able to train the pipeline on a new (clean) data set without any issues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_pipeline = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Submitting to Kaggle\n",
    "\n",
    "The following code will write your predictions on the test dataset to a CSV, which you can submit to kaggle.  You may need to modify it a little to suit your needs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "\n",
    "test_data = pd.read_csv(\"ames_test.csv\")\n",
    "\n",
    "submission_df = pd.DataFrame(\n",
    "    {\n",
    "    \"PID\": test_data[\"PID\"], \n",
    "    \"SalePrice\": final_pipeline.predict(test_data).reshape(-1,)\n",
    "    }\n",
    ")\n",
    "\n",
    "timestamp = datetime.isoformat(datetime.now()).split(\".\")[0]\n",
    "\n",
    "submission_df.to_csv(\"submission_{}\".format(timestamp), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your final kaggle submission should achieve a test-set RMSE threshold of 30,000 or lower. Write your best test-set RMSE (as shown on kaggle) here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_test_RMSE = 53240\n",
    "\n",
    "_ = ok.grade('qkaggle')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diagnostics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make some plots to investigate how well your models fit the data.  Pick an intermediate (not final) model for your diagnostics submission.  Provide commentary about patterns you notice and how you addressed them.  Include this comment on top of the cells you would like us to grade.\n",
    "\n",
    "    # DIAGNOSTIC_SUBMIT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DIAGNOSTIC_SUBMIT\n",
    "\n",
    "# Code for plot\n",
    "\n",
    "# Commentary\n",
    "diagnostic_commentary = r\"\"\"\n",
    "\n",
    "Put your commentary about diagnostics here, replacing this text.\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "display(Markdown(diagnostic_commentary))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**FYI:** Recall from lecture that stability is a measure of how robust your modeling procedure is to perturbations of the data.  While the formal definition is a little technical, the concept is intuitive: if you create pseudoreplicates of the data, the coefficients of your model shouldn't change too much since that would mean that your model is too sensitive to small changes in the training data.  Below, we use our pipeline to do a five-fold stability check.  This method is really a heuristic (as easily noted by the arbitrary choice of 5 folds).  To get a better assessment of your model, you could carry out a bootstrap analysis.  For this particular model, it would seem that the coefficients are not changing too crazily relative to the magnitude of their impact on home prices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import KFold\n",
    "\n",
    "fivefold = KFold(n_splits=5, shuffle=True)\n",
    "\n",
    "def calc_coefs(X, y, modeler):\n",
    "    model = modeler\n",
    "    model.fit(X, y)\n",
    "\n",
    "    return(model.steps[1][1].coef_[0])\n",
    "    \n",
    "np.vstack(calc_coefs(small_data.iloc[fold,:], small_data.iloc[fold, :][['price']], ex_pipeline2) \n",
    "          for (fold, _) in fivefold.split(small_data))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Properties of Least Squares\n",
    "Here we ask you to prove some basic properties about least squares.  While the focus of the class isn't the mathematical machinery behind data science, we want to at least motivate how theory can inform application."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 1\n",
    "Let $y$ be an $n \\times 1$ response vector and $X$ be an $n \\times p$ full rank design matrix **with a column of 1s**.  We use the least squares procedure to fit $y$ on $X$:\n",
    "$$\\hat y = X\\hat\\theta$$ where $\\hat\\theta = (X^TX)^{-1}X^Ty$.  The residuals are given by $e = y - \\hat y$.\n",
    "\n",
    "##### Part a\n",
    "Show that $\\sum_{i=1}^n e_i = 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Part b\n",
    "Show that $e$ is in the null space of $X^T$.  In other words, prove that $X^Te=0$.  Note that this is property is where the name \"Normal Equations\" come from: $e$ must be normal (orthogonal) to the space spanned by the columns of $X$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Part c\n",
    "Your friend fits a linear model of sale price on home square footage with an intercept as shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "model = lm.LinearRegression(fit_intercept = True)\n",
    "model.fit(small_data[[\"sqft\"]], small_data[[\"price\"]])\n",
    "print(\"Intercept:\", model.intercept_[0])\n",
    "print(\"Slope:\", model.coef_[0,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "She wants to know if the stochastic model $Y = X\\theta + \\epsilon$, where $\\epsilon$ is a mean 0 vector independent of the columns of the design matrix $X$ is plausible.  One assumption is that `sqft` must be independent of the noise term $\\epsilon$.  To test for this, your friend writes the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_independent(variable, error):\n",
    "    # Inputs\n",
    "    #   variable: n x 1 numpy array with variable of interest\n",
    "    #   error: n x 1 numpy array estimates of the error term epsilon given by y - y_fitted\n",
    "    # Outputs\n",
    "    #   boolean, True if the variable passes test for independence\n",
    "    return sum(variable * error)[0]\n",
    "\n",
    "n = small_data.shape[0]\n",
    "sqft = small_data[[\"sqft\"]].values.reshape(n, 1)\n",
    "fitted = (small_data[[\"price\"]] - model.predict(sqft)).values.reshape(n, 1)\n",
    "test_independent(sqft, fitted)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "She concludes that since this value is very small, `sqft` and the noise are most likely independent of each other.  Is this a reasonable conclusion? Why or why not?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**SOLUTION:** Her test will almost always computationally yield small values since by part (b), the inner product should theoretically yield 0.  In fact, it is only due to numerical reasons that her code does not produce 0. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 2\n",
    "\n",
    "Centering takes every data point and subtracts the overall mean from it.  We can write the transformation function $\\phi$ as:\n",
    "\n",
    "$$\\begin{align}\\phi(X) &= \\left[\\begin{array}{c|c|c|c} X_1 - \\bar{X}_1 & X_2 - \\bar{X}_2 & \\dots & X_d - \\bar{X}_d \\end{array}\\right] \\\\\n",
    "\\phi(y) &= y - \\bar{y} \\end{align}$$\n",
    "\n",
    "where $\\bar{X}_j$ is the arithmetic mean of the $j^{th}$ column of $X$ and $\\bar{y}$ is the average of the responses.  Show that if a bias/intercept term is included in a regression after centering, then it will always be 0.  This, of course, means that adding a column of 1s to your design matrix after centering your data might be a little silly.\n",
    "\n",
    "**Hint:** You will want to use what we've proved in Question 1a."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "for_assignment_type": "solution"
   },
   "source": [
    "**SOLUTION:**\n",
    "Let $X$ be an $n \\times d$ design matrix with an intercept column.\n",
    "\n",
    "**Claim:** \n",
    "Centering a variable makes it mean 0.  \n",
    "**Proof:** \n",
    "$n^{-1}\\sum_i\\left(x_i - \\bar{x}\\right) \n",
    "= n^{-1}\\sum_i \\left(x_i\\right) - n^{-1}n\\bar{x} \n",
    "= \\bar{x} - \\bar{x} = 0$\n",
    "\n",
    "**Proof of main problem:**\n",
    "Let $y$ and $x_1, \\dots, x_j$ be centered variables.  From question 1, we have that the sum of the residuals sums to 0.\n",
    "$$\\sum_{i=1}^n e_i = \\sum_{i=1}^n \\left( y_i - \\theta_0 - \\sum_{j=1}^d \\theta_j x_{ij} \\right) = 0$$\n",
    "Where $\\theta_0$ is the intercept and the remaining $\\theta_j$ are the other regression coefficients.  Dividing by $n$ on both sides,\n",
    "$$n^{-1}\\sum_{i=1}^n \\left( y_i - \\theta_0 - \\sum_{j=1}^d \\theta_j x_{ij} \\right) = 0$$\n",
    "We then use commutivity of addition to get\n",
    "$$\\begin{align} \n",
    "0 = \n",
    "n^{-1}\\sum_{i=1}^n y_i - n^{-1}\\sum_{i=1}^n \\theta_0 - n^{-1}\\sum_{i=1}^n\\sum_{j=1}^d \\theta_j x_{ij} \n",
    "&= \\bar{y} - \\theta_0 - \\sum_{j=1}^d \\theta_j n^{-1}\\sum_{i=1}^n  x_{ij} \\\\\n",
    "&= 0 - \\theta_0 - \\sum_{j=1}^d \\theta_j \\bar{x}_j \\\\\n",
    "&= \\theta_0 - \\sum_{j=1}^d \\theta_j \\cdot 0 \\\\\n",
    "&= \\theta_0 \\end{align}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Submitting your assignment\n",
    "Congratulations, you're done with this homework!\n",
    "Run the next cell to submit the assignment to OkPy so that the staff will know to grade it. You can submit as many times as you want, and you can choose which submission you want us to grade by going to https://okpy.org/cal/data100/sp17/.  After you've done that, make sure you've pushed your changes to Github as well!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = ok.submit()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
