{ "cells": [ { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "UQ9Cw8pMSu8u" }, "source": [ "# Multivariate Linear Regression in Python with Numpy\n", "\n", "Welcome to one more tutorial! In the last post (see <a href=\"/20180809linearRegressionScratch/#more\">here</a>) we saw how to do a linear regression on Python using barely no library but native functions (except for visualization).\n", "\n", "In this exercise, we will see how to implement a linear regression with multiple inputs using Numpy. We will also use the Gradient Descent algorithm to train our model.\n", "\n", "The first step is to import all the necessary libraries. The ones we will use are:\n", "- Numpy - for numerical calculations;\n", "- Pandas - to read csv and data processing;\n", "- Matplotlib - for visualization" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": {}, "colab_type": "code", "id": "GpQjC6oYzgGo" }, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "tHhIs2BNUbq_" }, "source": [ "Now import the data that we want to work on. This data (hypothetical) consists in the following information from real state properties:\n", "- Size (the size of the housing in square meters);\n", "- Nr Bedrooms -> the number of bedrooms;\n", "- Nr Bathrooms -> the number of bathrooms\n", "- Price -> The price of the house, in terms of thousands of dollars (or any other currency since the data is hypothetical)\n", "\n", "**Hypothesis**\n", "The price is linearly correlated with the size, nr of bedrooms and nr of bathrooms of a housing.\n", "\n", "We will check validity of the above hypothesis through linear regression.\n", "\n", " Pandas function `read_csv()` is used to read the csv file 'housingprices.csv' and place it as a dataframe." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 195 }, "colab_type": "code", "id": "pibqK2ndzqUv", "outputId": "df5e4e94-dea9-44ac-9282-05e83397c299" }, "outputs": [ { "data": { "text/html": [ "<div>\n", "<style scoped>\n", " .dataframe tbody tr th:only-of-type {\n", " vertical-align: middle;\n", " }\n", "\n", " .dataframe tbody tr th {\n", " vertical-align: top;\n", " }\n", "\n", " .dataframe thead th {\n", " text-align: right;\n", " }\n", "</style>\n", "<table border=\"1\" class=\"dataframe\">\n", " <thead>\n", " <tr style=\"text-align: right;\">\n", " <th></th>\n", " <th>Size</th>\n", " <th>Nr Bedrooms</th>\n", " <th>Nr Bathrooms</th>\n", " <th>Price</th>\n", " </tr>\n", " </thead>\n", " <tbody>\n", " <tr>\n", " <th>0</th>\n", " <td>84</td>\n", " <td>1</td>\n", " <td>1</td>\n", " <td>43.747</td>\n", " </tr>\n", " <tr>\n", " <th>1</th>\n", " <td>45</td>\n", " <td>4</td>\n", " <td>1</td>\n", " <td>30.100</td>\n", " </tr>\n", " <tr>\n", " <th>2</th>\n", " <td>73</td>\n", " <td>1</td>\n", " <td>3</td>\n", " <td>39.481</td>\n", " </tr>\n", " <tr>\n", " <th>3</th>\n", " <td>34</td>\n", " <td>2</td>\n", " <td>3</td>\n", " <td>23.908</td>\n", " </tr>\n", " <tr>\n", " <th>4</th>\n", " <td>31</td>\n", " <td>4</td>\n", " <td>3</td>\n", " <td>24.144</td>\n", " </tr>\n", " </tbody>\n", "</table>\n", "</div>" ], "text/plain": [ " Size Nr Bedrooms Nr Bathrooms Price\n", "0 84 1 1 43.747\n", "1 45 4 1 30.100\n", "2 73 1 3 39.481\n", "3 34 2 3 23.908\n", "4 31 4 3 24.144" ] }, "execution_count": 3, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "df= pd.read_csv('housingprices.csv')\n", "df.head()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 34 }, "colab_type": "code", "id": "16jQsby5z2mS", "outputId": "d97157d1-ee47-422e-8a96-49fd12d39ae5" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Index(['Size', 'Nr Bedrooms', 'Nr Bathrooms', 'Price'], dtype='object')\n" ] } ], "source": [ "print(df.columns)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "MOQ2MXGZVO8N" }, "source": [ "A good practice, before performing any computation, is to check wheter the data contains invalued values (such as NaNs - not a number). This can be done using `pd.isna()` function, which returns a dataframe of True or False values. Since we want to summarize the results for each column initially and know wheter there is AT LEAST one invalid value, we can use the `any()` function, which returns True if there is any invalid number, otherwise False." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 101 }, "colab_type": "code", "id": "-EhYjCnMz72D", "outputId": "68ef50fa-6650-44c4-f46c-435562d9785c" }, "outputs": [ { "data": { "text/plain": [ "Size False\n", "Nr Bedrooms False\n", "Nr Bathrooms False\n", "Price False\n", "dtype: bool" ] }, "execution_count": 5, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "pd.isna(df).any()" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "4P-mDCSsVx6e" }, "source": [ "No invalid value was found in the dataframe.\n", "\n", "Split the dataset into inputs (x) and output(y). Use the method `values` to transform from a DataFrame object to an array object, which can efficiently managed by Numpy library." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 34 }, "colab_type": "code", "id": "ludmulkp0AYb", "outputId": "e2a06ad0-563f-4b91-8502-f7887f57016b" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "99\n" ] } ], "source": [ "x = df[['Size','Nr Bedrooms','Nr Bathrooms']].values\n", "y = df['Price'].values.reshape(-1,1)\n", "m = len(y)\n", "print(m)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "qkkpM-cQWIrv" }, "source": [ "Let's generate a simple visualization the price in relation to each input variable." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 286 }, "colab_type": "code", "id": "fdvXkeW50L8j", "outputId": "2cdb5b40-3e09-486a-cace-46606011e6ef" }, "outputs": [ { "data": { "text/plain": [ "[<matplotlib.lines.Line2D at 0x7f428efc30b8>]" ] }, "execution_count": 7, "metadata": { "tags": [] }, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD8CAYAAAB9y7/cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnX+QHPV14D9vRyMYYR8rjM4Hi0Cy\nTaQgC1izAXxK5QzYCIMRW4Ax2CQ4xYVyEjvYuDaWLhSSCAny7Tk2qfLlrNjEJHaEQOjWEpDoCIhy\nhViClVdCLKBYyIAYwMgRK//QAqvdd390z2pmtnumZ7Z7unv6faqmdvpNz3zffuc7r7/9vu/7nqgq\nhmEYRvvREbcChmEYRjSYgTcMw2hTzMAbhmG0KWbgDcMw2hQz8IZhGG2KGXjDMIw2xQy8YRhGmxLI\nwItIp4hsEJHnReQ5EfmwiJwgIo+IyE/cv7OjVtYwwsbGttHOBJ3B3wX8s6ouBM4CngOWA4+q6unA\no+6xYaQNG9tG2yL1drKKyPHATuB9WnayiOwBPqKqr4nIScDjqrqg1medeOKJOm/evOlrbRge7Nix\n4+eqOifo+WGNbRvXRtQ0OrZLzAhwznzgAPB3InIWsAO4GXivqr7mnvM68N56HzRv3jwGBwcb1dEw\nAiEiLzX4llDGto1rI2qaGNtAMAM/A/gQ8AVV3S4id1F1y6qqKiKetwIichNwE8Cpp57ajI5GBhgY\nKtK/ZQ+vjoxycmeBvqUL6O3uirrZpse2jWsjLKIc+0F88K8Ar6jqdvd4A86P4mfu7Svu3ze83qyq\na1W1R1V75sxp+A7DyAADQ0VWbNxNcWQUBYojo6zYuJuBoWLUTTc9tm1cG2EQ9diva+BV9XVgv4iU\nfJAXAc8Cm4AbXNkNwA9C0cjIHP1b9jA6Nl4hGx0bp3/LnkjbtbFtxE3UYz+IiwbgC8D3RWQmsA/4\nfZyLw30iciPwEnBNKBoZmePVkdGG5CFjY9uIjajHfiADr6o7gR6Ply4KRQsj05zcWaDoMaA7RBgY\nKkbqi7exbcSJ39g/ubMQyufbTlYjdvqWLqCQz02Rj6u2yhdvGLFwwULv9Rs/eaMEddEYxrSoFSlQ\n+vvl+3YxXrUvo9wfGUOUjZEiYorEmhZbnz/QkLxRzMAbkTIwVGT15mHePDw2KSuOjNJ3/y6ACiP/\npfU7PT+jFFlQWowqHZe/38g2pWiUtI0RL/dMLXmjmIvGiIzSj67cuJcYm1BWbRqukPn5HXMisUTZ\nGOkhrkispGMG3ogMrx9dOSOjlYbfyxdfyOemuG1KtCjKpikGhoosWfMY85c/xJI1j9k6QsTEHImV\nWMzAG5HRyI+r5D8dHRsnJwJAV2eBO69czOxZec/3dPrI4ybGjVvTJq0XJr+7v7CiUdKKGXgjMur9\nuEqGe2CoSN+GXZN+x3FV8jmZXCTzy4dXJ09ebKTVXZDmC1Pf0gXkO6RClu9wxlCSqVK5rrzhzw/n\nYwxjKn7hjwD5nLDy8kUArN48zNh4pbUeG1dWb3Z89NWunBJ+8riJeuEsKtJ6YZqk2iiGZCSj5NPn\neecx8pM3ihl4IzJ6u7u488rFdLkz+XLXS//VZ01GN3gtwpbLS++rxk8eN2nTt0Sa/dj9W/Z4ThKS\nfnG6o3cx159/6uTYyIlw/fmnckfv4lA+38IkjUjp7e6adpia3yKrnzxu0qZviah3VUZJmi9Od/Qu\nDs2gV2MzeCN2Ogs+i6iuvMvHwPjJ4yZt+pbwi2JKuh8bbJHVD5vBGzUpRbcUR0bJiTCuSlfZLsEw\ndg+uWraIvvt3MTZxdIab7xBWLXN89H1LF1RsYoFkG5606Vui9L2lbTcopLfPAT7ztz/iiRcOTh4v\nef8JfP8PPhzKZ5uBN3yp3h1YcjGUoisGXzrIAzuK0949WM+wpM3wpE3fcsJwqcVBWvu82rgDPPHC\nQT7ztz8KxcjXrckaJj09PWqlzdLBwFDRMzdMOaUZfTVdnQWeWH5hlOp5IiI7VNUrM2Sk2Lg2mmXe\n8od8X3txzWWTz5sd2+aDN6ZQmrnXWxRM4w5Tw8gS5qLJOF4+9HopBkqIeG82Kl/Yqv78CxbOYevz\nB6bcRqcxE2Atbh3Yzbrt+xlXJSfCdefNjSxSIkza7XvIOmbgM4xXBr7qxc5aeBn38oUtr8//3raX\nJ88N25efFG4d2F3xf46rTh4n2cinNSNjCbs4TcVcNBnGa6Ye1Lh7MXtWnjuvXDw5I//yfbvq3gmM\njo2zbvv+dO+grGLd9v0NyZNCmneylqe7KKVZ6NuwKxVpFqLEDHyGCXvr/KyZMyaNexAffol28+Wn\ndaNTWlMsQP10F1nFXDQZxi8KplmKI6M1owIa1SOtm1T8/p+kpypIq95QP91FVrEZfEYZGCo2bdzD\n/L0X8jmuO29uandQevG+ObMakieFtN55GP6Ygc8gJRdKsxx/bN43S2QtSomUujoLCEfzvd/Ru3gy\nKVm5PK0LZPsOHG5InhTSmiQN6qe7SCrHzfT+HfnJG8VcNBkkaBikH4dGx/j6p86ejFgIMr8r5HM1\njXZad1B6kdaZcFr1hvrpLpJKPtcBTP0tOvLpYwY+g9RbvKznmz+5s8DgSwd5/dBbdY27uOdnKWQt\nrb7stOoN6U1VcMinpoGfvFHMwGcQv7SwpRQD82sslOZzwhu/qIxn9yPMvNZp4vz3zZ6SX6QkTzJp\nnsFDOu8Co07RbD74DNK3dAH5XFV5s9zR8mZ+g6tDAIWxidqfH3bRgrQx/OovG5InhTT74NNK1Cma\nbQafVaomZeNuzPCX1u/k+EKefE4q4ooL+RzHzOioWyavPEFSVklbicESaZ/Bp5GoXUtm4DNI/5Y9\nU3asTnA0ZnhkdIx8hzB7Vp6Rw2OTg+5L63fW/Fyb6aUbv9xC9rVGS5SuJTPwGSTIDtGxCWXWzBkM\n3XbxpKxU+MOP686bG4p+aWf2rLznBpvZs5Idsuc3UbcJfHoxH3wGCbqAU30huGDhHM/zhOwuqHpx\n2ZknNSQ3jKgwA59BvBZ2vKi+EGx9/oDveWbcj+LXT37ypOB3h5H0Ow/DH3PRJISwU53W+rzqhZ3j\nC3l+/c6RKYuq1Sv5aa5c30rS2k8rL19E34ZdFeMgnxNWXp7szUIl0pouOEq9zcAngLDzcHvmed+w\ni1Wbhjk0enTRtLysXpBBFnXMbruQ1n5K62YhSG8u+6j1NgOfAGrl4W7mS/bM8z6uk2F6XoMoyEp+\nmivXt5ILFs7x3Ajmt4aRJNK4WQjC/w21iqj1NgOfAMK8pR8YKgbK393MIErzDK+VpNUHD+l1c6TV\nLRa13mbgE0Ajt/S1foCNZolsZhCldYbXStJqbNLq5oD0usUSk6pARHIiMiQiD7rH80Vku4jsFZH1\nIjIzFI0ySNDtyqUfYHVZsrNX/z/mL38oUIm8cpI++FtBFOPar1+T3t9pLtnn5/5Kulss6lQFjYRJ\n3gw8V3b8VeDrqvoB4E3gxlA0yiC93V2++dAHhoosWfOYrwEv+daV2lvK8x2V2xHNdz5J6OM66h9t\nVKS5ZF9a3WK1fvthEMhFIyKnAJcBfwHcIiICXAh82j3lHmAV8DehaJVimvVherk+qm+Zm80JkhPh\nU+fOZevzB1LnW42SqMZ1b3cXgy8dZN32/YyrkhPhqnOS79pKc7rgNF+coiSoD/4bwJ8C73aP3wOM\nqOoR9/gVINmjtwV4hifev4vVm4cnc7pcsHBOIEM7MFTky/ftCiXR07gqD+woprpKUkREMq4Hhoo8\nsONoScRS//ecdkKi+z/NycbSenEaGCpWFCop2QwIZ92jrotGRD4BvKGqO5ppQERuEpFBERk8cCDZ\nt0vTxTM8cUJ58/DYpM/8e9tervChr9i4m4GhYsV7SheKMH9YafGltooox3VafdlpThec1ovTqk3D\nUxL/jU0oqzYNh/L5QXzwS4BlIvIicC/OLexdQKeIlO4ATgGKXm9W1bWq2qOqPXPmJHvBY7o0EyXh\n9cOfTkm9WikIsn67WkVk4zqt7oK0Gsk0E3Vq6boGXlVXqOopqjoPuBZ4TFU/A2wFrnZPuwH4QSga\npYzyRdCOJmc61ReGRi8UOZGKBZo0z8RaRZTjOq393+UT5eMnN5LPdJKNfQVnYWovju/yO+GolB6q\nwxabnel0ViVzaiScrpDP8bVrzuKnay7jieUX0tvdZTOx6THtcZ3W/q9X6SvJ+F06k31JjZ6GNjqp\n6uPA4+7zfcC54auUHvxcKTkRJlQ9k3h5MXJ4jPnLH5pcdPVKCeDF7Fl5Vl6+aMpiTFeNmqvGVMIe\n16nu/+qhmuxr0iR+aiZdfbcKpqc8DCxdcIOUu2T8fKoTqvx0zWXsXHkx/VefVfcz1X2U7xwsj431\nu7WfNXOG50p7WuOw24W09r9Xpa+xCU384nCaifrCZKkKfPCKZwcCzazLXSy93V11KyGVU1p0Lblb\nAOYvf8jzXD9fveWMiZe09n9aUyykmairf5mB98AvJ8cxMzrqGnevmVpQl0uJV0dGKy4wHT4xvrV8\n9ZYzxmiU4wt5z+iN4wtW8CMq3vKxCX7yRjED74FfHHMtAy3gO1MrHQfduNQ5K193B2sabvmzSlqT\ndvkF+SQ8+AdI77rH6NhEQ/JGMR+8B43eknZ1FiqiWLzo7e5iIoBxL+RzqOK7eBtFvgojXNK60WnE\nw1VQS54k0rruETU2g/fAL4Xn7Fl53hqbqPjx5juEw+8cqYiC8TO8fp9biropvf9L63d6vr+0eGsk\nm7T6stOachfSu+4RtQ/eZvAe+BWlVoWrzumajG7pLORBqEhFUJ16oDzq5tdvH5kSZ+wVx57WdLOG\nQ1q/v7Sm3E0zKy9f5Ln3IKw6uDaD96B01V+9ebji6joyOsb6J/fzrmOdbvvlW0em+MdHx8b58n27\nJo/LEwmNjI7RgXN1LiUf85plWGm8dJPWkn1pTbkLzkTqlvU7KXmuiyOj3OLeCSd5Fh/1nYcZeB9K\n4Y3Vt0+l5GFQe8eis6imU+KKJ3DuBGq5WtJ6u2k4PLjrNV/5Hb2LW6xNcNLqWgJYsfFpqpclJ1x5\n0n83UUa8mYGvwXQGdq2ImyCJhCzMMb1EnUAqKjp9/MHVqTSSSNTRKGnFfPA1SLrP1DDC5JBPtIyf\n3Eg+mZ/B16rA5OdLnS5hrZAbySTqyIio8JvrZnsOHD3NVoELQqZn8F5FrMujYMJYXIpyhdxIJlFH\nRhhTOW6mdx0EP3lSqGeDpkumDXy9DSnTLdDQ1Vmg/+qzKgrq9l99lvnW25ze7i773lvM4Xe817z8\n5Ekh6k1xmXXR3Dqw29eAl3LBTIdSWKMtlmYT+95bS1o3aUUduZTJGfytA7tr+tZFmCx82wydhbyl\nEjCMFpLWTVp+EUphRS5lcga/bvv+mq9PKIHyxoCTquBdx86ouXHJMIxoSesmLT8zE1bxr0wa+OmU\nTjODbhjJI62btA757I3wkzdKagx8mKFEOZ/86kHo/6QtlhlG0kirDz5qvVPhgw8zlGhgqMgxM5pP\ncG3G3TCSR1rTBUe9dpAKAx9WKFHpQnHYY/tyrqO+0U968QDDmA5pjSUHZ+JVXsc4LTUTol47SIWB\n9/OjFRsMZ/S6UIAzGL72ybMmDbiXqU/DbMBIDuVpopeseSy0jStRks95mwM/uTF9Mh8mOTBUpKNG\nzbBGXDW1OrO3u2vyNq/aOz97loU9GsGJendiVES94Bclae3zqGsHJNrAl760WguiJVdNkBlTvc70\nm+HPmjnDjLsRmLSW7EtroRJIb59n2gfvZ3CrKV2t612963VmWkOtjGTht0N6uqkvoiatm4Ugvb/d\nTPvgg345OZFAV+96nRn1rjIjG/it1wdYx4+VtG4WgvTefWTaBx/kyynkc74unOpOqteZUe8qM7LB\nhM948ZMnhbTOgsEJk/TK4Jn0wIjM+OC9fOh+t4bHzcxVhEL5hS9Wd1K9GXqaF5kMY7qkdRY8SfUF\nNOEXVIj+wpQIA++1Av7F9Tv5x+3eCcHeqopjD7LJYWCoyIhPZZrSDD31A9xIBH6emIR7aOhbumDK\nfpBcR/JnweCs11XXPx6b0MQvsgKRXpgSYeD9FlP9bmnHVSsWUwGuOqeLnBtOmRPhqnOOpmstXUD8\n+q00Q0/rbjgjWfiNs6RPKAdfOsh41Y9ufEIZfOlgTBoFJ60L21FfmBJh4KfzJYyOjbN68zAP7ChO\n+uLHVXlgR3EyiqZeNE5php7W3XBGssj57NvwkyeF7/vcMfvJk0RaF7ajXvdIRLKx6ST/AjzrX5ai\naHq7u2p2VvUM3Qo1GNPFbyxPZ4y3gjQHGaR1YTsTycaiGvglw+7XWTkRm6EboeO36G+5jIxqonYL\nJ8LAT2fgF/I5Ogve0TElw+7XiV+7xlL/GuGT1rWcWXlvc+AnTxJ+NsBPnhSidgvH7qIZGCpy+J0j\nnq/NnpXnsjNPYuvzBybzwF+wcE7FcelHs2Lj7go/e/kPqtRZYeWTN4xapHW8/eWVZ/LF9Ts95Uln\n1bJF9N2/q2LBMt8hrFq2KEatghGlWzhWA1+KbqleAO0s5Fm1bFHD/3StH5T51o1Wktbxls8JY+Na\ncZwG0npRjZq6Bl5E5gJ/D7wXJ9JrrareJSInAOuBecCLwDWq+mYjjftFtxx3TOPJvdL6gzLiI8qx\nnUb6t+ypMO4AY+M6GayQdMwGTCWIc+0I8GVVPQM4H/hjETkDWA48qqqnA4+6xw2R5q3RRlsQ2dhO\nI/Z7bD/qGnhVfU1Vf+w+/yXwHNAFXAHc4552D9DbaOO2c9SIkyjHdhqx32P70dDyuIjMA7qB7cB7\nVfU196XXcW5zvd5zk4gMisjggQOVWen6li4gX7UTIZ+SrdFGe9Ho2K41riGdFZ3SGv1j+BPYwIvI\nu4AHgC+q6i/KX1NVxWcntqquVdUeVe2ZM8cjeVjVGs7YhLJq03AqfhBGe9DM2K41rgeGivRt2FWR\nW6lvw67Ej2nbyd1+BIqiEZE8zg/g+6q60RX/TEROUtXXROQk4I1GG/da1AEYGR2bzDFjg8uIkijG\n9urNw56Llas3Dyd+PNtCZXtRdwYvIgJ8B3hOVf+q7KVNwA3u8xuAHzTaeK3FmzSU2zLSTVRj2yt1\nRi25YURFEBfNEuB3gQtFZKf7uBRYA3xMRH4CfNQ9boh6ize2em9ETGRj2zCSQF0Xjar+K/6prC+a\nTuN9Sxd4bnQqYav3RpRENbY7C3lGPIrEJH3bvNF+xJpkorSoM9uj0pKt3htpZdWyRZ7RYWnYNm+0\nF7Hnoikt6gwMFW2bsdEW2LZ5IymItjDZs4gcAF5qUXMnAj9vUVvNYPpNDy/9TlNV70K+EVJnXCe9\nH/1Iq96QXt1r6d3U2G6pgW8lIjKoqj1x6+GH6Tc9kq5fibToWU1a9Yb06h6F3slP9GwYhmE0hRl4\nwzCMNqWdDfzauBWog+k3PZKuX4m06FlNWvWG9Ooeut5t64M3DMPIOu08gzcMw8g0ZuANwzDalLYx\n8CKSE5EhEXnQPZ4vIttFZK+IrBeRmTHq1ikiG0TkeRF5TkQ+LCIniMgjIvIT9+/sGPX7kogMi8gz\nIrJORI6Nu/9E5G4ReUNEnimTefaZOPy1q+vTIvKhVurqhZf+aUBE5orIVhF51h0TN8etUxDcMfuk\niOxy9V4dt06NUG2/wqJtDDxwM05FnhJfBb6uqh8A3gRujEUrh7uAf1bVhcBZOHomoiyciHQBfwL0\nqOoHgRxwLfH333eBS6pkfn32ceB093ET8Dct0rEW32Wq/mnAr4xh0nkbuFBVzwLOBi4RkfNj1qkR\nqu1XKLSFgReRU4DLgG+7xwJcCGxwT4mt7JqIHA/8Dk5aWlT1HVUdIVll4WYABRGZAcwCXiPm/lPV\nHwIHq8R+fXYF8PfqsA3odPO4x4aP/omnRhnDRON+979yD/PuIxURJNX2K0zawsAD3wD+FJhwj98D\njKjqEff4FeIbpPOBA8Dfubdg3xaR4whY8jBqVLUI/C/gZRzDfgjYQXL6rxy/PusC9pedlxR9U01V\nGcPE47o5duIUaHlEVVOhN1PtV2ik3sCLyCeAN1R1R9y6+DAD+BDwN6raDfyaKndMrZKHUeP6sa/A\nuRCdDBxHClwLcfZZFqhVxjCpqOq4qp4NnAKcKyIfjFunekRtv1Jv4HGKNiwTkReBe3FcC3fh3KaX\nsmWeAsRVEPMV4JWy2cQGHIP/s5IbodmShyHxUeCnqnpAVceAjTh9mpT+K8evz4rA3LLzkqJvKvEp\nY5gaXBfoVlIwUcHDfonI98L68NQbeFVdoaqnqOo8nMXBx1T1Mzhf8NXuaU2VFAxJv9eB/SJSSm5/\nEfAsIZQ8DImXgfNFZJa7dlHSLxH9V4Vfn20Cfs+NpjkfOFTmyjEaoEYZw0QjInNEpNN9XgA+Bjwf\nr1b18bFf14fZQNs8gI8AD7rP3wc8CewF7geOiVGvs4FB4GlgAJiNs07wKPAT4F+AE2LUbzXOj+EZ\n4B+AY+LuP2AdzprAGM5d0I1+fYZTlembwAvAbpyIoLjH4hT949YpoN6/jeP6ehrY6T4ujVuvAHqf\nCQy5ej8D3Ba3Tk38D5P2K6yHpSow2g4RuRso+Tan+GHdWepdwKXAYeCz6kaOiMgNwK3uqXeo6j3V\n7zeMtJB6F41hePBdavtfPePmReQEYCVwHnAusDLODWiGMV0CGfik78Q0jHK0fgy6X9z8UpzwuoOq\n+ibwCOlYqDMMT4LO4BO7E9MwmsAvbt7i6Y22om7R7bKdmJ8FZycm8I6IXIGzKADOrsLHga/U+qwT\nTzxR582b17SyhlGLHTt2/FxbVJNVRG7Cce9w3HHHnbNw4cJWNGtklGbHdl0DT+VOzLNwdjneTMCd\nmOU/hFNPPZXBwcFGdTQMBoaK9G/Zw6sjo5zcWaBv6QJ6uysn1yIStKC7X9x8kaOTlpL8ca8PUNW1\nuAUaenp61Ma1ESUNjO0KgrhoprUTU1XXqmqPqvbMmdPygvdGGzAwVGTFxt0UR0ZRoDgyyoqNuxkY\nanovk1/c/BbgYhGZ7a4pXezKDCOVBDHwSd+JabQ5/Vv2MDo2XiEbHRunf8sez/NFZB3wI2CBiLwi\nIjeKyOdE5HPuKQ8D+3Bi/P8W+CMAVT0I/DnwlPu43ZUZRiqp66JR1ddFZL+ILFDVPRzd6fgszm7C\nNSRnp6PRhrw6MtqQXFWvq/V57h3nH/u8djdwd2MaGkYyCeKDB/gC8H236MM+4PdxZv/3iciNwEvA\nNdGoaGSdkzsLFD2M+cmdhRi0MYz0EMjAq+pOoMfjpYvCVccwptK3dAF9G3YxNn50mSefE/qWLqjx\nLsMwbCerkQ6ql/Atw4Zh1CWoi8ZoIUFCArNE/5Y9jE1UWvSxCaV/y55M94th1MMMfMIohQSWokZK\nIYFAZo2Zl/+9ltwwDAdz0SSMRkMCDcMw/DADnzAaDQk0DMPwI3YXjfmbK7GQQMMwwiLWGfzAUJG+\n+3dVbEHvu3/XdLagp56+pQso5HMVskI+l+mQQL9BarefhlGbWH8jqzYNe0ZHrNo0HJNG8dPb3cWd\nVy6mq7OAAF2dBe68cnGm72omGpQbhuEQq4tmZHSsIXlW6O3uyrRBNwwjHOwu10g8s2flG5IbhuEQ\nq4G3H64RhJWXLyKfkwpZPiesvHxRTBoZRjqI1UWz8vJF3HLfTsrd8B1C5n+4tw7sZt32/YyrkhPh\nuvPmckfv4rjVio2Su8qirQyjMWIPk8x1CBNlSaRyHVLj7Pbn1oHdfG/by5PH46qTx1k38mbQDaMx\nYnXR9G/ZU5EhEGBsXDO9a3Pd9v0NyQ3DMPyI1cDbrs2pjKt3mkQ/ueGNiFwiIntEZK+ILPd4/esi\nstN9/LuIjJS9Nl722qbWam4Y4RGri2bWzBy/fmfcU24YzSIiOeCbwMdwSk4+JSKbVPXZ0jmq+qWy\n878AdJd9xKiqnt0qfQ0jKmI18F7GvZbcyC4NprQ4F9irqvsARORe4AqcMpNeXAesDF1pw4gZi4M3\nEk8phXJ5SosVG3fXSmnRBZQvWrziyqYgIqcB84HHysTHisigiGwTkd4Q/gXDiAUz8EbiiTiF8rXA\nBlUtb+A0Ve0BPg18Q0TeX/0mEbnJvQgMHjhwIAw9DCN0zMAbiaeJxfgiMLfs+BRX5sW1wLpygaoW\n3b/7gMep9M+Xzlmrqj2q2jNnzpwa2htGfJiBNxLP8QXvnc1+cuAp4HQRmS8iM3GM+JRoGBFZCMwG\nflQmmy0ix7jPTwSW4O+7N4xEE+siayHfwejY1JyAhbxdd4yjiM/eNz+5qh4Rkc8DW4AccLeqDovI\n7cCgqpaM/bXAvaoVMai/CXxLRCZwJkBryqNvDCNNxGrg3/Iw7rXkRjYZOeyTddRHDqCqDwMPV8lu\nqzpe5fG+fwOyu2XYaCtinSo3cettZBC/alZW5cowahOrgW/01tvIJn1LF3hmk8xylSvDCEKsBr6Z\nW+92J+dzdfOTZ4bqTA2WucEw6hKrge/0yfvuJ88ClotmKv1b9niWdsxyUjrDCEKsBt7PZmXYltHp\ns/7gJ88ClpTOMJojVgN/yKf2qp88C9i6xFRskdUwmsNcNAnD1iWm0rd0AYV8ZYbRQj5ni6yGUQdz\n0SQMm61Opbe7izuvXExXZwEBujoL3HnlYqvwZBh1iHWjk7lopnLBwjkVJfvK5VnGSvYZRuPEOoO3\n2epUHnr6tYbkhmEYfsRq4P1mpVmerb7p42v3kxuGYfgRq4F/cJf3rNRPbhiGYQQnsIEXkZyIDInI\ng+7xfBHZ7hY1Xu+mZW2IER9fu588C1gcvGEYYdHIDP5m4Lmy468CX1fVDwBvAjeGqVhWWXTyuxuS\nZ4WBoSJL1jzG/OUPsWTNY7XK9RmG4RLIwIvIKcBlwLfdYwEuBDa4p9wDNFy7crZPvLufPAts2/dm\nQ/IsMDBUpG/DroqarH0bdpmRN4w6BJ3BfwP4U6CUqP09wIiqHnGPfYsa12Ll5YvIdVRu0cx1CCsv\nX9ToR7UNlotmKqs3DzM2XpWLZlxZvXk4Jo0MIx3UNfAi8gngDVXd0UwD9YoTa5Xhqj7OGn4ZCTKc\nqaCpyCIRuURE9rhrRMs9Xv8otH8CAAAN80lEQVSsiBwQkZ3u47+XvXaDiPzEfdwQxv9gGHEQZAa/\nBFgmIi8C9+K4Zu4COkWktFHKt6hxreLEqzcPU5UkkAkl0zMzv8tbti97jSEiOeCbwMeBM4DrROQM\nj1PXq+rZ7qPkfjwBWAmcB5wLrBSR2S1S3TBCpa6BV9UVqnqKqs7DqWH5mKp+BtgKXO2edgPwg0Yb\nt5hvIyLOBfaq6j5VfQdnYnJFwPcuBR5R1YOq+ibwCHBJRHoaRqRMJw7+K8AtIrIXxyf/nXBUyjYd\nPr4YP7nhSRewv+zYb43oKhF5WkQ2iMjcRt5bz/VoGEmgIQOvqo+r6ifc5/tU9VxV/YCqflJV3260\n8bxP637yLJDzMeR+cqNpNgPzVPVMnFn6PY28uZbr0TCSQqymdGyiMXkWsD6ZShPhtEVgbtnxlDUi\nVf2PsknJt4Fzgr7XMNJChufKRlpYefkiz6LbNcJpnwJOd3dbz8RZO9pUfoKInFR2uIyjm/i2ABeL\nyGx3cfViV2YYqSPWdMHGVGblOzjsMV2flWG/VW93F4MvHWTd9v2Mq5IT4VO/Ndc3fbCqHhGRz+MY\n5hxwt6oOi8jtwKCqbgL+RESWAUeAg8Bn3fceFJE/x7lIANyuqgej/Q8NIxpiNfCCd/hflt3NM2fk\nPA38zBk5j7OzwcBQkQd2FCc3e42r8sCOIj2nnVDLyD8MPFwlu63s+Qpghc977wbuDkl9w4iNeCs6\nNSjPAlYEZSr9W/YwOjZeIRsdG6d/y56YNDKMdBCrgc/5VJL2k2cBK4IylVdHRhuSG4bhEKuBt7wr\nU+lbusBzQTHLBabtomcYzRGrge/y+YH6ybPCeFVirerjrGGVvwyjOaxkX8JYtWmY6iXWCVeeVbY+\n771T1E9uGIZDrAbeCkxPxapcTcV88IbRHLEaeEs2ZgSh02fHqp/cMAyH7O6eMVKD35p7htfiDSMQ\nsRp4K25hBMHcVobRHLbRyUg8tl/CMJrDNjoZicf2SxhGc9hGJyPx2H4Jw2gOm8Ebicf2SxhGc9gM\n3kg8D+7y3hfhJzcMw8Fm8AnD3BFTsSgaw2gOm8EnjL6lC8hXVdjOd2Q72VgziMglIrJHRPaKyHKP\n128RkWfdotuPishpZa+Ni8hO97Gp+r2GkRYs2VgC8cpFk2UarckqIjngm8DHgTOA60TkjKrThoAe\nt+j2BuB/lr02qqpnu49l01TfMGIjVgPft3QBhXxlpaJCPpfp2erqzcOMT1Rlk5xQVm/ObrKxJmqy\nngvsVdV9qvoOcC9wRfkJqrpVVQ+7h9twimsbRlsRq4Hv7e7iqnO6Jn3uORGuOqfLtwxbFrD8PFPp\n7e7iU781t2Kc1KrJCnQB+8uOX3FlftwI/FPZ8bEiMigi20SkdxqqG0asxGrgB4aKrH9qf0WtzfVP\n7WdgqBinWkbCGBgqsv7JqnHyZDjjRESuB3qA/jLxaaraA3wa+IaIvN/jfTe5F4HBAwcsbbGRTGI1\n8Ks3DzNWVcxibDzb7ojOgk/mRB95Fli1aZixKrfV2ITWypFfBOaWHZ/iyioQkY8CfwYsU9W3S3JV\nLbp/9wGPA93V71XVtarao6o9c+ZYPL6RTCxdcMJYtWyRZxTNqmW+/ua2p4kwyaeA00VkvojMBK4F\nKqJhRKQb+BaOcX+jTD5bRI5xn58ILAGene7/YBhxMCNuBYxKSn7l/i17eHVklJM7C/QtXZDpdYlG\nUdUjIvJ5YAuQA+5W1WERuR0YVNVNOC6ZdwH3i+Pbf9mNmPlN4FsiMoEzAVqjqmbgjVQSq4EXvDNH\nZnebkxEWqvow8HCV7Lay5x/1ed+/AYuj1c4wWkOsBt7SBU9lYKjIio27GR0bB6A4MsqKjbsBMjuL\n7yzkPd0xWV6XMIwgWMGPhNG/Zc+kcS8xOjZO/5Y9MWkUP4tOfndDcsMwHKzgR8Io+hSS9pNngW37\n3mxIbhiGg9VkTRgdPrcvfvIsYDmLDKM5zMAnjAkfm+UnzwKWddQwmsN88Ebiue68uQ3JDcNwiNXA\nf+b8UxuSZwHbyTqVO3oXc/35p1bkorn+/FO5o9eiGQ2jFrGGSd7Ru5ifHvgVT7xwcFK25P0nZPqH\nu2rZIvru31WxNT/rO1nBGStZHheG0Qx1Z/AiMldEtrrFEYZF5GZXfoKIPCIiP3H/zm608YGhIj9+\n+VCF7McvH8p0srHe7i76P3kWXZ0FBCc3fv8nz8psDLxhGM0TZAZ/BPiyqv5YRN4N7BCRR4DPAo+q\n6hq3Ys5y4CuNNF4r5jvLBq23O9spkw3DCIe6M3hVfU1Vf+w+/yXwHE5u7SuAe9zT7gEazpv9qk9s\nt5/cMAzDCE5Di6wiMg8ndep24L2qWipr/zrw3kYbP9mnNJ+f3DAMwwhOYAMvIu8CHgC+qKq/KH9N\nVRWfDai1CiNYgWlvBoaKLFnzGPOXP8SSNY9lek3CMIzmCWTgRSSPY9y/r6obXfHPROQk9/WTgDe8\n3lu3MEJ10HvGg+BLycaKI6MoR5ONmZE3DKNRgkTRCPAd4DlV/auylzYBN7jPbwB+0Gjj/Vv2eFZ0\nynJiLUs2ZhhGWASJolkC/C6wW0R2urL/AawB7hORG4GXgGsabdwWWadifWIYRlgEiaL5V1UVVT1T\nVc92Hw+r6n+o6kWqerqqflRVD9b7rGpskXUq1ifhICKXiMgeEdnrhvFWv36MiKx3X9/uBhCUXlvh\nyveIyNJW6m0YYRJrqoK+pQso5HMVskI+l+lFVuuT6SMiOeCbwMeBM4DrROSMqtNuBN5U1Q8AXwe+\n6r73DJwarouAS4D/7X6eYaSOWA18b3cXd165uGLX5p1XLs70Jh/rk1A4F9irqvtU9R3gXpx9G+WU\n7+PYAFzkrjddAdyrqm+r6k+Bve7nGUbqiL3otu3anIr1ybTpAvaXHb8CnOd3jluk+xDwHle+req9\n9mUYqaSlBn7Hjh0/F5GXfF4+Efh5K/XxISl6gOniRS09TmuVEiJyE3CTe/i2iDzTqrariOt7iXM8\nZPF/bspH21IDr6oegfAOIjKoqj2t1CfJeoDpMg09ikB5svhTXJnXOa+IyAzgeOA/Ar4XVV0LrG1A\np0iIq237n1vfdjPvs4pORjvyFHC6iMwXkZk4i6abqs4p38dxNfCYuyN7E3CtG2UzHzgdeLJFehtG\nqMTugzeMsHF96p8HtgA54G5VHRaR24FBVd2Es3nvH0RkL3AQ5yKAe959wLM4mVT/WFXHPRsyjIST\nJAO/Nm4FXJKiB5guXgTSQ1UfBh6ukt1W9vwt4JM+7/0L4C/C1iki4mrb/ucUtC1qlekNwzDaEvPB\nG4ZhtCktNfAicreIvOEXUiYOf+1uE39aRD4Ukx4fEZFDIrLTfdzmdV5IuniWRKw6J/J+CahHS/pF\nRI4VkSdFZJery2qPc3xTDUTFdNIftKDtW9zv7mkReVREQgkZrddu2XlXiYiKSChRJkHaFZFrysbr\nP4bRbpC2ReRU97cy5Pb3pSG1G759VNWWPYDfAT4EPOPz+qXAP+EkDT4f2B6THh8BHmxRn5wEfMh9\n/m7g34EzWt0vAfVoSb+4/+e73Od5nAIz51ed80fA/3GfXwusj1inHPAC8D5gJrDLo38i0Slg2xcA\ns9znfxhG20HaLRsvP8TZINbTov/3dGAImO0e/+cW9vVa4A/d52cAL4bUduj2saUzeFX9IU7Egh9X\nAH+vDtuATnFzzrdYj5ah/iURy4m8XwLq0RLc//NX7mHefVQvFvmlGoiK6aQ/iLxtVd2qqofdw204\n8fuRt+vy5zi5fN4Koc2g7f4B8E1VfRNAVT3rUUTUtgL/yX1+PPBqGA1HYR+T5oP32mIe1zbxD7su\ngn8SkUWtaFAqSyKW09J+qaEHtKhfRCQnTnrqN4BHVNW3T1T1CFBKNRAVQb6DqHRq9Pu/EWemF3m7\nrptgrqo+FEJ7gdsFfgP4DRF5QkS2icglLWx7FXC9iLyCE6n1hZDarkfDdiBJYZJJ4sfAaar6K9e/\nNoBzSxgZUqMkYiupo0fL+kWd2POzRaQT+L8i8kFVjSsdQGoQkeuBHuC/taCtDuCvgM9G3ZYHM3DG\n3kdw7lZ+KCKLVXWkBW1fB3xXVb8mIh/G2U/xQVWdaEHbDZG0GXygbeJRo6q/KLkI1ImnzovIiVG1\nJ94lEctpSb/U06PV/eK2MwJsxUndW85kn0hlqoGoaCT9Qdg6Bfr+ReSjwJ8By1T17Ra0+27gg8Dj\nIvIijl94UwgLrUH+31eATao6pk7Wz38nnMlGkLZvBO4DUNUfAcfi5KmJmobtQNIM/Cbg99zV4vOB\nQ6r6WquVEJH/UvKdisi5OP0UifFw2/EqiVhO5P0SRI9W9YuIzHFn7ohIAfgY8HzVaX6pBqJiOukP\nIm9bRLqBb+EY97D80TXbVdVDqnqiqs5T1Xk4vv9lqtpU3pSg7boM4MzecScZvwHsm2a7Qdt+GbjI\nbfs3cQz8gRDarkfjdiCM1d8GVonXAa8BYzhX4BuBzwGfc18XnEINLwC7CWFFvkk9Pg8M46ygbwP+\na4R98ts4izZPAzvdx6Wt7peAerSkX4AzcSIkngaeAW5z5bfjGBBwflT34+RrfxJ4XwvG76U4M8UX\ngD9rpU4B2v4X4Gdl392mVrRbde7jYY3NAP+v4LiHnnV/E9e2sK/PAJ5wfwc7gYtDajd0+2g7WQ3D\nMNqUpLloDMMwjJAwA28YhtGmmIE3DMNoU8zAG4ZhtClm4A3DMNoUM/CGYRhtihl4wzCMNsUMvGEY\nRpvy/wEQI3K7R7SsrwAAAABJRU5ErkJggg==\n", "text/plain": [ "<Figure size 432x288 with 4 Axes>" ] }, "metadata": { "tags": [] }, "output_type": "display_data" } ], "source": [ "fig,axs = plt.subplots(2,2)\n", "axs[0, 0].plot(x[:,0],y,'o')\n", "axs[0, 1].plot(x[:,1],y,'o')\n", "axs[1, 0].plot(x[:,2],y,'o')" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "G9gEqC2nWugI" }, "source": [ "Linear correlation can be evaluated through Pearson's coefficient, which returns a value between 0 and 1. 0 means there is no correlation, while 1 means perfect correlation. Everything in betweeen indicates that the data is somehow correlated, though usually a correlation of more than 0.8 is expected for a variable to be considered a predictor, i.e an input to a Machine Learning model." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 67 }, "colab_type": "code", "id": "33fF8V9f0T9w", "outputId": "1af07804-d0a3-4e4e-a8ad-f469bec8491a" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Correlation between x1 and y = 0.97\n", "Correlation between x2 and y = 0.24\n", "Correlation between x3 and y = 0.11\n" ] } ], "source": [ "# Check correlation of each input with the input\n", "from scipy.stats import pearsonr\n", "print(f'Correlation between x1 and y = {pearsonr(x[:,0],y[:,0])[0]:.2f}')\n", "print(f'Correlation between x2 and y = {pearsonr(x[:,1],y[:,0])[0]:.2f}')\n", "print(f'Correlation between x3 and y = {pearsonr(x[:,2],y[:,0])[0]:.2f}')" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "wWSIcpExXLlY" }, "source": [ "Results above shows that only the size shows high correlation with the price. Even though, we will keep the other variables as predictor, for the sake of this exercise of a multivariate linear regression.\n", "\n", "Add a bias column to the input vector. This is a column of ones so when we calibrate the parameters it will also multiply such bias." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 34 }, "colab_type": "code", "id": "8DAvADy6DqF_", "outputId": "89d1e214-8769-43b9-9466-70efc48ea146" }, "outputs": [ { "data": { "text/plain": [ "(99, 4)" ] }, "execution_count": 9, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "# Add a bias to the input vector\n", "X = np.concatenate((np.ones((len(x),1)),x),axis=1)\n", "X.shape" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "mtlwbE6kYgDO" }, "source": [ "Another important pre-processing is data normalization. In multivariate regression, the difference in the scale of each variable may cause difficulties for the optimization algorithm to converge, i.e to find the best optimum according the model structure. This procedure is also known as **Feature Scaling**." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 185 }, "colab_type": "code", "id": "Zn3E8AsONg-j", "outputId": "a2e58710-a10e-4623-afc6-879378c80fb3" }, "outputs": [ { "data": { "text/plain": [ "array([[1. , 0.84693878, 0. , 0. ],\n", " [1. , 0.44897959, 0.03061224, 0. ],\n", " [1. , 0.73469388, 0. , 0.02040816],\n", " [1. , 0.33673469, 0.01020408, 0.02040816],\n", " [1. , 0.30612245, 0.03061224, 0.02040816],\n", " [1. , 0.30612245, 0.02040816, 0.01020408],\n", " [1. , 0.91836735, 0.03061224, 0.02040816],\n", " [1. , 0.33673469, 0.02040816, 0.01020408],\n", " [1. , 0.33673469, 0.03061224, 0. ],\n", " [1. , 0.79591837, 0.02040816, 0.02040816]])" ] }, "execution_count": 10, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "Xnorm = X.copy()\n", "minx = np.min(X[:,1:])\n", "maxx = np.max(X[:,1:])\n", "Xnorm[:,1:] = (X[:,1:]-minx)/(maxx-minx)\n", "Xnorm[:10,:]" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 50 }, "colab_type": "code", "id": "spvruJQ8N7uD", "outputId": "79ebe17d-a1ec-49d3-ef01-f82a8af0ba32" }, "outputs": [ { "data": { "text/plain": [ "array([0.62943377, 0.28902469, 0.5230232 , 0.13457221, 0.14045897,\n", " 0.11666251, 0.92968321, 0.15153405, 0.19299077, 0.75762035])" ] }, "execution_count": 11, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "ynorm = y.copy()\n", "maxy = np.max(y)\n", "miny = np.min(y)\n", "ynorm = (y-miny)/(maxy - miny) \n", "ynorm[:10,0]" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "1cCUls3IY6QF" }, "source": [ "Now we are ready to start working on the model. It is reasonable to give an initial guess on the parameters and check this initial guess makes sense or is it complete nonsense. In this exercise, the choosen vector of parameters $\\theta$ is \n", "\n", "$$\n", "\\theta = \\begin{bmatrix}0.4 \\\\ 0.4 \\\\ 0.4 \\\\0.4 \\end{bmatrix}\n", "$$\n", "\n", "Notice that there are 4 parameters, which corresponds to the bias + the input variables from the data. The linear regression model works according the following formula.\n", "\n", "$$\n", "Y = X\\cdot \\theta\n", "$$\n", "\n", "Thus, $X$ is the input matrix with dimension (99,4), while the vector $theta$ is a vector of $(4,1)$, thus the resultant matrix has dimension $(99,1)$, which indicates that our calculation process is correct.\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 286 }, "colab_type": "code", "id": "bwx_TwkqDLu8", "outputId": "6e7bbba7-d8e4-4afc-bcc4-f4b1b94556c0" }, "outputs": [ { "data": { "text/plain": [ "[<matplotlib.lines.Line2D at 0x7f428ef8aef0>]" ] }, "execution_count": 12, "metadata": { "tags": [] }, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmYlNWVx/Hv6aaABpRmi7IKIqII\nCqRVDOMWNYg6QIwRNcYlKiaRic4YMzqZJEZjRBnXSFRi3IiRcWEQlYgLxAWX0IgBQVHEjXZDpFGh\nobc7f9wqurqo5a3uarrrrd/nefrprrfeqrpF6enb5z33XHPOISIi4VLU2gMQEZHcU3AXEQkhBXcR\nkRBScBcRCSEFdxGREFJwFxEJIQV3EZEQUnAXEQkhBXcRkRBq11ov3LNnTzdw4MDWenkRkby0dOnS\nz51zvTKd12rBfeDAgZSXl7fWy4uI5CUzez/IeUrLiIiEkIK7iEgIKbiLiISQgruISAgpuIuIhFDG\n4G5md5rZZ2b2eor7zcxuNrM1ZrbczEbnfpgiIvlv7rIKxk5byKBLH2fstIXMXVbRYq8VZOZ+N3Bs\nmvvHA0OiX1OAW5s/LBGRcJm7rILL5qygorIKB1RUVnHZnBUtFuAzBnfn3HPAF2lOmQjc67yXgVIz\n652rAYqIhMH0BaupqqlrdKyqpo7pC1a3yOvlIufeF/gw7va66LEdmNkUMys3s/L169fn4KVFRPLD\nR5VVWR1vrp16QdU5N9M5V+acK+vVK+PqWRGR0OhTWpLV8ebKRXCvAPrH3e4XPSYiIlGXjBtKSaS4\n0bGSSDGXjBvaIq+Xi+A+DzgjWjUzBtjknPs4B88rIhIak0b15eoTR9C3tAQD+paWcPWJI5g0KmkW\nu9kyNg4zs/uBI4CeZrYO+A0QAXDO3QbMB44D1gBbgLNbZKQiInlu0qi+LRbME2UM7s65UzPc74AL\ncjYiERFpNq1QFREJoVbr5y4iUgjmLqtg+oLVfFRZRZ/SEi4ZN3SnpGYU3EVEWkhsVWps8VJsVSrQ\n4gFeaRkRkRays1elxlNwFxFpITt7VWo8pWVERJopPq/etSSCGVRuqaHIjDrndji/pValxlNwFxFp\nglhAr6iswoBYCK+sqtl+TrLA3pKrUuMpuIuIBJQqoO8YwhsrNqPeOVXLiIi0NYmVL5kCerx653h3\n2vEtM7AUdEFVRCSAZJUvQe2MHHsiBXcRkQCaWuGys3LsiRTcRUQCyDT7tuj30pII3TpFdkrnx3SU\ncxcRSSPVRVRg++2+O/FCaVAK7iIiKSS7iNqWA3o8BXcRkQTxs/VEscC++NJv7/yBZUHBXUSE9OmX\nRDujfUBzKbiLSMHLtoa9NUobs6XgLiIFK136JZXWKm3MloK7iBSkxNl6EG39Imo8BXcRKUjZrDgt\niRS3Wr16Uym4i0hBynRRNF9KHlNRcBeRgtSntCRlrj1fA3o8tR8QkYJ0ybihlESKGx0riRRz4+SR\nLL7023kd2EEzdxEpULHgHdtBaWf2Wt8ZFNxFpGBNGtU3NME8kdIyIiIhpOAuIhJCSsuISKjFVqHG\n8upH7tOLRW+uD2WePZ6Cu4iETqomYBWVVfzl5Q+2n1dRWcVlc1YAhC7AKy0jIqESaysQq2HP1ASs\nqqaO6QtWt/zAdrJAwd3MjjWz1Wa2xswuTXL/ADNbZGbLzGy5mR2X+6GKiGTWlI2sc97Ct3YbrJwL\nG94Bl+nXS8vImJYxs2JgBnAMsA5YYmbznHOr4k77b+AB59ytZjYMmA8MbIHxioik1ZRAnbGF77vP\nwbavYJ/jgz3hwt/Bizf7n3ftB8O/C9/5Xdbjao4gM/eDgDXOubXOuWpgNjAx4RwH7Br9uSvwUe6G\nKCKS2dxlFYydtjBjGiZRoBa+NVUw+zRY/1bmJ9zwDrx0C4w4GY6/HvqVQW21v885eOK/shxh0wS5\noNoX+DDu9jrg4IRzLgeeNLN/AzoDR+dkdCIiAWRq3xvfBCxltUzFUnj4PGjXEQYdBv0Pgo//CUf9\nBvqMgo6l8OjP4Kz5UJRmXtx9T/j+3TDocCgphQPPabiv+mvovX8u33pKuaqWORW42zl3nZkdAswy\ns+HOufr4k8xsCjAFYMCAATl6aREpVEE228jYBGzLF/57t0FQ2t/PrpfeBa/cCh12hVGnQ88hMO73\n8MhPYemdcOC5yZ9r29fQoQsMS0xuRHXYBQ44JYt32HRBgnsF0D/udr/osXjnAMcCOOdeMrOOQE/g\ns/iTnHMzgZkAZWVlrXOVQUTyTnyteteSCGawcUtNxr1ODdJvZP35Gnj4HPjRE9CpO5zxiD9esxU+\netUH/F17+2MjT4MVD8BTl8Pe46GkG7z+EKyPVtq4enjtrzDpj8Fz8y0oSHBfAgwxs0H4oH4KcFrC\nOR8ARwF3m9m+QEdgfS4HKiLhFySIf1m1jd5sYCO9mr/X6Wt/gU9WwNZNEIk7N9IR9vhW43PN4IQb\n4cGzoGojdOoBj/8crAiKot0lu/aDPqODv+EWlDG4O+dqzWwqsAAoBu50zq00syuAcufcPOBi4E9m\n9u/4z+As51qp/kdE8lJi3ryyqmb7fbFgUspX3BSZwd5F6zhm27V8TaeUz5fxQmldLbx2Pwz5Duyy\ne7BBdh8EU/7uAz3Av5VD1/4Nt9uQQDl359x8fHlj/LFfx/28Chib26GJSFilmqGn0o5axhS9wbTI\nn+hFJZfXnkk9RYyyt1nmhuxwfqDNNtY8DV9/4nPq2YgP5KVt99qh2g+ISItIFsArt9TQtSTC5upa\naur8fDx+hp5MZ6p4ucNUdrEqKlwPvl/9G5a7wdwYuYXDipZzyLZb2EZ7IMu9TpfNgs7fgL3HNfu9\ntkUK7iKSM4k9XYbYh1zR7imurPoh1USAzME80WZKuL32BN52fXmhfgSb8bnx/607kknFLzKh+EUe\nrDsi+63xBh0GAw+F4khW48kX1lqp8bKyMldeXt4qry0izbDkDvjyI1/HPWAMtOsAJK81vz1yPeOK\ny7mi5ofcWTc+q5f5afFc9iz6hF/UTKGeou0XVUu3/xVQzdMll9Fj106UXvTyjnnv1U9An5HB8+l5\nwsyWOufKMp2nmbuIZOejZf5C5PPXwTeGwZmPQueeO/R06WfrObpoKdtcO37c7lFm1R1DzQ4hx+EL\nFhtrTw0/iixgFYNxNUWpZ+Xl6+Gxi+CDl2GPQxqOv/k4PHAGjPg+fPc2+OoT6LKb/wWw7StY9Qjs\n911o3zl3/y5tjIK7SCHZ8gW8eo+vFIkZMAYGHep/rqtJnaZYeg/sPhwmzoBxV8NbC2DeVJg1Cc58\ndIeeLt8qep06ijir5j+pcD23B/Z+9hlTih/nmtpT+FW7WSx3g/lr3VHbH1cSKWbWgR/Q89VNHHb6\nL3l3rzQL3vc/GZ6+3F8cjQX3NU/7csXeB8D4a2Hrl/DnY3zNek0VrFsCOCiKwAGTs/v3yyMK7iKF\n5O/T4B+3Nz526M99cF+3FB7+EZz4J7/0Pt7mDTD/Ehj1A+j7Tea+8RXTF/Rgry0XMv2TmUy99iEc\nfRo95IG6I1lYN5rP6br9WBe2cGvkRgbYeh6OnECv2k38rt2dFLcv4S9V39reDqCs/FroMQT2TLMA\nCfzM+4JXGlIvD54Nq+dDr6Fw+sPQcVe/4vTgH8Oz1/jnPPQ/YM8jG36hhZSCu0hYPH+9X5Bz/HV+\ntWUyuw+HQy+GI+KaV8Vy1e18xQl3HQfHXu2X2Mfue/VuqNvGD5aPZPELj2/Pf1dwAIduvSFareIo\n5Wsq2YX21FBNhA3RwF7asR3TuZ6D3HK62hZePngGj4w/HWpOgfsnc+W7f+TKo9vD2J/5xlsVS2H8\n9PQ9XGJigd05iHTydesn3OBXkMbe3yEX+K8CoguqImHw5cdw0/5QV+0X1Zx8L/RtwkrJqo0w53x4\newHsP5lH+1/C9KffY/bWKayt783pNb9M+dDzih/jvHbzuaD6Z/w2cg+vRkbT5YSrGvLkT18OL9wA\n//LvcPTlDQ+s3gyPXuSX9vcZBZPv839dHHaJ78UijQS9oKrgLhIGr93vOxaeOBOe/JWvYLngHw3L\n4qs2wutzYOQP/NL6dOrr4fnrcIuu4ur6M/iwtpRb29/EudUX83T9N1M+bKh9wG2RG9jDPqPIHEy4\nBUb/sOGE2m2w5hk/sy5OkjT45HU/zpCnS5pL1TIihWTkqbDX0dClFww8DDZ/5gN7bTW4OnhlJvz9\n99D/YJ+aSaeoCA6/hPMXd+HpL/sxuXgRK+oHsrB+VNqHrXYDmFB9FddGbmd4u3X0H3FS4xPadYB9\n0mzSlmlckhXN3EXy3ZYvUufY//af8MFLsPF93wjr1PszPl3yNrrJSxaTKYkUc/V3hzNpdL9A50t2\ngs7ctUG2SD7bugluHgmLb05+/55HwBfvwdZKXxWTQeLm0g12DOyxI6UlEbp1imD4ni5XnzhCgb0N\nUFpGJJ/9408+wO95ePL7h46HHz8Hn66Efqnz5TGZNpeO39Eoq6X+stMpuIvsDB8v95ssH3JB9u1h\nv/oEyu/yKyq/sU/D8dpt8MptsNcxfsFOKt0G+q8A0m0urYCeXxTcRXaGJy6F9xf7jR0O+Wmwx2ze\nAItv8LPz2q3+8Wc91nD/yrmweX3w54sT37GxT9y+oqmuwPUtLUm/o5G0Ocq5i6Sy8T24fhi890Lz\nn+v0OdBzb3jqV74PSoxz/mJnfUIqpL7OL5l/8RYYNskvPOo+yFe/xCy9yz/nnkdmNZT4vLoDKiqr\n+MvLH6TchzTjphfSJmnmLpLKS3+ELytysyFDpCOc8xTMPML3PTn/edj2JTx0Nnz8T+hY6uu7e+7t\nV48Wt4Px1/gFSfGpmHiT74NNH2Sd5smUV4+nVEz+UnAXSWbb1/DP+2H4SdkF97efAgyGxDW7enUW\nfLEWjvo1TJ4F906ET1dA/zHQriMc/VvY8DasfQ7eeBQGHAJDjvFfidaVQ4ddodfe0LmH/wooeYlj\nahk3l5Y2TcFdJJnls/3M+uDzff/ybV/5ZfPpbP0S7jsJxvy0cXBf9heo2wZH/wZ2HwEXrWhoNXvO\nkw3nOedfp+OuyZ+/ejPMOhF67+/TNuOuCtxiIFmv9Uwybi4tbZpy7hJedTXw8Hnwvz/0S+qDcs6v\n6OwzCvodCB/+A567zgfvmK1f7vicS+7w3/c/Oe68Tb7F7OC4GXCqHuJmqQN77HEHT4H3nvcLkzp2\nTX1ugmxSMaA8exgouEs41dXCw+f6ZlRvzPM9zINyzje2OurXPuAedD5Uf+XTNACfvQHXDYU55/pz\nAaq3wEszoi0AdvP5eud8+aOrg8FHpXq17Bz8E4h09v1ZegzOePrcZRWMnbYwbSqmb2kJp48ZQN/S\nksYLkZRnz2tKy0jbVrMVyu+Eb54F7TsFf1zl+/Dus3DMlfD2k74j4YjvQ4cumR9bVNS4B0q/b0Lf\nMvjHTDjwPOi1D3TqCa8/7Gf2Y34Cr94LWz73q0DffhIWXOa3eHtnIbTv4s/Lhc494Nynocs3Mp4a\nJBWjEsfwUnCXtm3Jn+DJ//a7Ax10XvDH9RgMU8uhc0/Y9wTYsjFzYK+vh9cfgvVvwpgLGl+sPPh8\nmHMevPecX9J/0XKY/QM/tj6jfapkj7F+N6DeB8BTv/ELjLrsDvsc39ArPRd2G5b27qAXTpV6CTcF\nd2m7nPN7XQKseDBYcP9irX/Mgef5wA7QfU+I9dXaumnHXLVzfmu2hVf6ssRe+/ide+INm+RTLLF6\ndDOY9MeG0sapS6C+xt/XvhOMPsOnaS5aDl13Xp+VoBdOVeIYfsq5S9tlBmf/zVeffPiKX+wT7/O3\nfdom3gs3wqKrofrrHZ9v8U1w82hYv7rhmHN+9n3fSVBVCZNug5+8uGPao117mHgL7BWXOy8p9Zti\njPud/6sgtvMP+F2MXB0s+XPT3nsTBblwGkvFKLCHm4K7tE11NT5wF0d8SmTfCb6XSsynq2DGwTD7\n1IbZ9KYKeO2vMOr0hq3X4u1zgl/+f88Ev5Ub+Jn86vk+GE8t933RYxtcBNF7fxj+vR2Pd9vD5+nX\nPBX8uZohyIVTUCqmkCi4S8va8kXjEsJ4a5+FhVc1VJzEW/4A3HQAVH7gm15NnuUX7sQsvtEH6ncW\nwrPX+mMv/gFcPYy9MPnr9RgMZzzit6K7d6J/7pJSOPcZv19nLvPi4PvAnPV4bp8zidRtehtTFUxh\nUXCXYOJ7mgSxdZNPj9y4P7wY7TWeGMTrquG5a/0in3j1dfD8dX5Xoa79G45veMfvFQp+A+Sz58MB\np/pd7Vc8BEvvhgNO8bPmVHYbBmfM9Z0W753oSyY7dQ+2EXO2IiVZ1aI3VaZUTEmkmBsnj1QqpsAo\nuEuDrV/CpnWNj9XV+j05r+63Y847mZf+CHcdB9cOhmenweAjfdri87fhznG+GVddjQ/0g4/yFSZP\n/tIH25gXrocv3vFlhbG+KVu+gFsO9Bsn11b7BT39D4Ljr4dv/ZuflQ84OPMqUvDVLOc9A6N+mPnc\nPJCpTa9m64VJ1TLil7W/cru/4Hjyvb66o77e120/9CNf5jfi5IYZcV0NWDF8Eu1R/vlb/mIj+Auf\n1V/72u/h3/O13gBfvAufvQkPnAE9h/qUyndvg3+9CW4dC3/7hX/tl2+Fhb/zNen7TmgYY6fufkOK\nl2b4apgzH/XjbN8JvnOlP+eMR4K/594HpO+Bnkf6lJYkTcmohr2wKbgXutV/g3k/8xsqDxnnywDB\nz6aX/cWnTibd5i80gs9xPzIVarb4nerBdzLc9hV02AVOujP5Bcnug3wwn32qLzc84jI/K+85BA7/\nhS9DfPNxv6nFvhP8ayamSkZ8379+h11glz4t92+SZy4ZN3SH8kddOJVAwd3MjgVuAoqBO5xz05Kc\nczJwOX4Xrn86507L4TgFfCoj2118Yurr/WPjH//ZG35m3n2wv2A5YEzDfbv09jPbcb/3FSEx7bv4\nuvHSATDocN+mdte4QJuu0mSf4/zzff4WHPaLhuNjL/SVMAMOgb3HQ32tb3mbaN9/9bn1I/+rZXLk\neSZ+w42uJRE6Roqo3FJDH9WwC2AuWaVC/AlmxcBbwDHAOmAJcKpzblXcOUOAB4BvO+c2mtk3nHOf\npXvesrIyV15e3tzxF453n4fHLoIpz6Zfafn52/4iZKRjw7GqjXD3Cf4C30l3QWn0IuXKuX7ziB89\nCbv2btnxS07Erz6N7WcaUxIpVn69AJjZUudcWabzgkx/DgLWOOfWOueqgdnAxIRzzgNmOOc2AmQK\n7NIEr90HG9ZAuw6pz3EOHjgT/ny0z3HHdOgKpXv47oS3HwbvLPLH95vka7sV2PNCYslj4rSsqqaO\n6QtW7/hAKUhBgntf4MO42+uix+LtDextZovN7OVoGmcHZjbFzMrNrHz9+vVNG3Ehqq+HNc/AfidC\nUTt4brrfqm1TBbx2v29rW1XpUy6HXezrt2ceDq/P8dUvRUVw6l9h6lLfsXDWJHjuf/xzp/tlIW1K\nkNWn6SpnpLDk6oJqO2AIcATQD3jOzEY45yrjT3LOzQRmgk/L5Oi1W1/1Ft9edtCh8M2zG6dEcuHj\nZf6C597H+guX/5zta8hd9H/0Tj38Rco9D49WqIyGB37ot3Drvif89BW/QKfnXr4E8IlL/S+A5uTw\nZafJZgclbbAhMUGCewUQt5KEftFj8dYBrzjnaoB3zewtfLBfkpNRtnXtO8F+3/X9vV/8g6/+GHl6\n8ouCADVVvgqluL3f1T6+J0l9nV8A1Kl7w7G3FgDme4V33BXOmOfLFksH+ID+jf0aX2DsPsjv1/ns\ntf4iafzKy/adYcIfcvr2peVks4OSKmQkXpC0zBJgiJkNMrP2wCnAvIRz5uJn7ZhZT3yaZm0Ox9n2\n7f99X2e9S2949EIf6BN3tI9Z9QhMGwBX7Q6Lft/4vofPhZtHwddxaauSbjDipIYWtF37wnHXwrem\n+m3bklWOREr8tm57j8vN+5NWkSkVE/u7S4uVJFHG4O6cqwWmAguAN4AHnHMrzewKM4utMlkAbDCz\nVcAi4BLn3IaWGnSr2bTO57vrahuOfbEW/jwumhY5wm+kcMwVsPL/YO2i5M+z4iHoOsA3uHr1Xvjq\nU3/8jcdg5RzYWgkvz2g4f8xP4Ht3tNS7kjYo6A5KN0weyXvTjldrAdlBoJy7c24+MD/h2K/jfnbA\nf0S/8sez1/rWrt88K9j5rh5evMXvwlN2tj+24mH48GWf9wafwx57IQw6zO/BmWjz534hztif+eXv\nr/0VXroFDr0YHr8YdhsBR/3KPx78DL6kW+oUj4SOdlCSXCjclSCfroRFV/mLk0F17e9XcP79atj2\ntb8gueJB3x8lcUOGWGB//0V4/vqG4yv/z18IHX6S74ey34l+G7n1q30Qn3CzT6VESvxfCI9eCH86\notlvV/JHkEZgyq1LJoUR3Nct9VUm8Z6/zq+2POWvwZ5j0e/9L4LvXAlff+pn25++Dp+vTt7PO2bF\ng/DMb6MXRfH7bvbaF3bbz98+9GLfi+XT1/0mEX1H++OfroQ/jPL7cfY/OLv3K3kpaCpGuXUJIvx/\n61dvgbuPg15D4dyFPr2x4R0/gz5kqq9KefNx2LLBb42WTOWH/pdB2Tm+x8qwibD4Zt/hsKid34It\nlXG/94uH5kyB85+Df70ZNq9vKEHcbRj8+AXYbXjjssRuA333w/oa3/NFQiexfcDm6lpq6lJXCCsV\nI9kI/8y9fScYf42/4LnoKn/shRt8GeIhU31q5dVZMP8XjVd1xlt8E2A+Tw5w1G98IK7e7J8jfiPl\nRJES3+3QOd8RsXQADBzb+JzdR+xYb96+M4yfBr1H+vp5CZX41aYOqKyqSRvYlYqRbGXsLdNSdkpv\nmVXz/IXJklLf+fDVe+DU2X4WX/GqLy8Ev9JzxsFQt80H/aJiuPQDf9/aZ+HeCf7iZ6ytLSTfaDmd\nNx+H2afBpFthpHqqFapsFiTFaDNriRe0t0x40zIfvQYPnuk3Vx53FYy/Fj5+Df7vfDj/+YbADr5u\n/AcP+AAMjWfRT1/uvyduApHtDjv7HA+T7/MLjKQgZbMgKUapGGmqcAb36s0wbyp07gWH/dwfi3T0\n6ZGX/rjjzvYAe3zLfyWa8Ae/orTH4OaPa98Tmv8ckreC9IaJp1SMNEf4gvvna3xflfVvwsmzGi/t\n7zbQr+zMxu7Dczo8KQyJF0vNYOOWmrSPiRQZXTq2U092yYn8DO4v3uLTIgecuuPinvkX+/04T38Y\nBuvPWcm9ZIG7cktNoyAe32u9sip9UAfl1SX38i+4OwdvPeH39Vx8o+/CuP4Nf8FzwBiYOMOvJC0d\n0NojlRBKzJvHB+74n4OWKWiDDWkp+RfczfzmyKv/5jdSfvKX0LEUBh7qg3viSlGRHGhKlUsmmq1L\nS8q/4A4+wO9znO9vXvm+n6Wn27tTpAnSbWnXXKqCkZaWn8E9pqhIpYWSM+lWjOYysKsKRnaG/A7u\nIjmSLpfeVLHZfmncRVdVwcjOouAuQvY16DGlSaplFMSlLVBwFyH7jaVV5SJtnYK7FJxkdepBcuqx\nNIuqXCQfKLhLQUhV+ZIut64Vo5LPFNwltFIF9CCzdM3OJd8puEuoNCegxxioBl3ynoK7hEZiOWNT\na9P7lJbkblAirST8OzFJwWhqOWM8LTCSsNDMXUIjm3JGLTCSsFNwl7wXy7NnSsOolFEKiYK75KWg\nTb0U0KVQKbhL3si2EkYBXQqZgrvkhWwrYVTOKIVOwV3atKZukqFyRil0Cu7SZiXO1oNSOaOIgru0\nYdnUrevCqUhjgYK7mR0L3AQUA3c456alOO97wEPAgc658pyNUgpSprp1BXSR1DIGdzMrBmYAxwDr\ngCVmNs85tyrhvF2AC4FXWmKgUnj6lJakzLUroIukF6T9wEHAGufcWudcNTAbmJjkvCuBa4CtORyf\nFLBLxg2lJNJ44/OSSDE3Th7J4ku/rcAukkaQ4N4X+DDu9rrose3MbDTQ3zn3eA7HJgVu0qi+XH3i\nCPqWlmD42bp2PxIJptkXVM2sCLgeOCvAuVOAKQADBgxo7ktLCMXvkhTr9aJ6dZHsBQnuFUD/uNv9\nosdidgGGA383M4DdgXlmNiHxoqpzbiYwE6CsrKypHVklZFKtPK2orOKyOSsANFsXyVKQ4L4EGGJm\ng/BB/RTgtNidzrlNQM/YbTP7O/BzVctIOkFbCVTV1DF9wWoFd5EsZQzuzrlaM5sKLMCXQt7pnFtp\nZlcA5c65eS09SAmXbFsJZNPKV0S8QDl359x8YH7CsV+nOPeI5g9L8ll83rxrXL/02M8bt6TelDoZ\ntRIQyZ5WqEpOJc7KK6saAnn8z0GplYBI0yi4S07lYqs7rTwVaT4Fd8mJpnZvjFFAF8ktBXfJSrJ8\n+sYtNWl3Q8pEAV0k9xTcJa3EYL65upaaOh/G43PoTQnsJZFirTgVaSEK7pJSuoujQZUmqZap3FKz\nffWpArtIy1BwFyB1uqU5+paWqHWASCtRcJeczNATqYRRpHUpuBeYlpihx6jiRaTtUHAvAKn6uDRl\nhh4pMrp0bKccukgbp+Aectn2cUlHM3KR/KHgHnK5WDGqkkWR/KPgHnJN6ahYqnSLSN5TcA+5dJtM\nJ9IMXSQ8FNxDJGhrgNhtzdBFwkvBPc8FqYRxqExRpNAouOehoFvUxYsFdq0YFSkMCu55pjmljdqu\nTqRwFLX2ACQ7zSlt1HZ1IoVDM/c8EH+htKmLkNTrRaSwKLi3cYlpmExUCSMioODe5gVJw6gSRkQS\nKbi3cekughpoVi4iSSm4t1GxPHuqHLvKGkUkHQX3NihTnl0XR0UkEwX3Nihdnl15dREJQsG9DYlf\neZqMgVIxIhKIgnsbEaTkUYuQRCQorVBtIzKVPCrPLiLZ0My9jUhX8qg8u4hkS8G9FSTru66SRxHJ\npUBpGTM71sxWm9kaM7s0yf3/YWarzGy5mT1jZnvkfqjhEMutV0T7xFRW1bBxS03Sc5WKEZGmyhjc\nzawYmAGMB4YBp5rZsITTlgFlzrn9gYeAa3M90LAI2tWxb2mJtrwTkSYLkpY5CFjjnFsLYGazgYnA\nqtgJzrlFcee/DJyey0Hmq1Tryw7TAAAH10lEQVTb3mWikkcRaa4gwb0v8GHc7XXAwWnOPwf4W7I7\nzGwKMAVgwIABAYeYX4Jse5eJSh5FpLlyekHVzE4HyoDDk93vnJsJzAQoKytramvyNqs5uyTFKM8u\nIrkQJLhXAP3jbveLHmvEzI4Gfgkc7pzblpvh5YdMK0vTUd91EWkJQYL7EmCImQ3CB/VTgNPiTzCz\nUcDtwLHOuc9yPso2LNvNNOKpzFFEWkrGahnnXC0wFVgAvAE84JxbaWZXmNmE6GnTgS7Ag2b2mpnN\na7ERtzFN3dNU6RcRaUmBcu7OufnA/IRjv477+egcjytvpFtZCtr2TkRah1aoNlOf0pKUuXa1DRCR\n1qLg3gSJ9euRYqOmrqE2piRSrAVIItKqFNwDSle/HikyunWKKOUiIm2GgnsaqQJ6Yv16Tb2jU/t2\nLPv1d3byCEVEklNwJ3WbgHQBPVGmC6siIjtTwQb3IG0CsllhqpYBItKWhD6452JWnolq1kWkrQl1\ncE9cPdrUWXkysV8OKncUkbYolMG9Ob1e0lFAF5F8Ebrg3pxeL8kooItIPgpdcG9qrxe1CRCRMAld\ncM+mJFGzchEJq9AF93S9XjQrF5FCEbrgfsm4oTvk3NXrRUQKTeiCeyyAx2rbNUMXkUIUuuAOPsAr\nmItIIcu4E5OIiOQfBXcRkRBScBcRCaHQ5NzjG4TpIqqIFLpQBPfElgMVlVVcNmcFgAK8iBSkUKRl\nkrUcqKqpY/qC1a00IhGR1hWK4J6q5YB2RxKRQhWK4J5qFyTtjiQihSoUwf2ScUMpiRQ3OqbdkUSk\nkIXigqpaDoiINBaK4A5qOSAiEi8UaRkREWlMwV1EJIQU3EVEQihQcDezY81stZmtMbNLk9zfwcz+\nN3r/K2Y2MNcDTTR3WQVjpy1k0KWPM3baQuYuq2jplxQRyRsZL6iaWTEwAzgGWAcsMbN5zrlVcaed\nA2x0zu1lZqcA1wCTcz3YWP+Yisqq7fufgtoNiIgkCjJzPwhY45xb65yrBmYDExPOmQjcE/35IeAo\nM7PcDbOhf0xsf1SXcL/aDYiINAgS3PsCH8bdXhc9lvQc51wtsAnokYsBxiTrH5NI7QZERLydekHV\nzKaYWbmZla9fvz6rxwYJ3Go3ICLiBQnuFUD/uNv9oseSnmNm7YCuwIbEJ3LOzXTOlTnnynr16pXV\nQDMFbrUbEBFpECS4LwGGmNkgM2sPnALMSzhnHnBm9OeTgIXOucS0eLMk6x8TS+r3LS3h6hNH6GKq\niEhUxmoZ51ytmU0FFgDFwJ3OuZVmdgVQ7pybB/wZmGVma4Av8L8Ackr9Y0REgrMcT7ADKysrc+Xl\n5a3y2iIi+crMljrnyjKdpxWqIiIhpOAuIhJCCu4iIiGk4C4iEkIK7iIiIdRq1TJmth54v4kP7wl8\nnsPh5ItCfd9QuO9d77uwBHnfezjnMq4CbbXg3hxmVh6kFChsCvV9Q+G+d73vwpLL9620jIhICCm4\ni4iEUL4G95mtPYBWUqjvGwr3vet9F5acve+8zLmLiEh6+TpzFxGRNPIuuGfarDsszKy/mS0ys1Vm\nttLMLowe725mT5nZ29Hv3Vp7rC3BzIrNbJmZPRa9PSi6+fqa6Gbs7Vt7jLlmZqVm9pCZvWlmb5jZ\nIYXweZvZv0f/G3/dzO43s45h/bzN7E4z+8zMXo87lvQzNu/m6L/BcjMbnc1r5VVwj9usezwwDDjV\nzIa17qhaTC1wsXNuGDAGuCD6Xi8FnnHODQGeid4OowuBN+JuXwPc4JzbC9iI35Q9bG4CnnDO7QMc\ngH//of68zawv8DOgzDk3HN9W/BTC+3nfDRybcCzVZzweGBL9mgLcms0L5VVwJ9hm3aHgnPvYOfdq\n9Oev8P+j96XxZuT3AJNaZ4Qtx8z6AccDd0RvG/Bt/ObrEML3bWZdgcPweyPgnKt2zlVSAJ83fl+J\nkugubp2Ajwnp5+2cew6/50W8VJ/xROBe570MlJpZ76CvlW/BPchm3aFjZgOBUcArwG7OuY+jd30C\n7NZKw2pJNwK/AOqjt3sAldHN1yGcn/sgYD1wVzQddYeZdSbkn7dzrgL4H+ADfFDfBCwl/J93vFSf\ncbPiXb4F94JjZl2Ah4GLnHNfxt8X3cowVOVOZnYC8Jlzbmlrj2UnaweMBm51zo0CNpOQggnp590N\nP0MdBPQBOrNj2qJg5PIzzrfgHmSz7tAwswg+sN/nnJsTPfxp7E+z6PfPWmt8LWQsMMHM3sOn3b6N\nz0WXRv9sh3B+7uuAdc65V6K3H8IH+7B/3kcD7zrn1jvnaoA5+P8Gwv55x0v1GTcr3uVbcA+yWXco\nRPPMfwbecM5dH3dX/GbkZwKP7OyxtSTn3GXOuX7OuYH4z3ehc+4HwCL85usQzvf9CfChmQ2NHjoK\nWEXIP298OmaMmXWK/jcfe9+h/rwTpPqM5wFnRKtmxgCb4tI3mTnn8uoLOA54C3gH+GVrj6cF3+e/\n4P88Ww68Fv06Dp9/fgZ4G3ga6N7aY23Bf4MjgMeiP+8J/ANYAzwIdGjt8bXA+x0JlEc/87lAt0L4\nvIHfAm8CrwOzgA5h/byB+/HXFmrwf62dk+ozBgxfHfgOsAJfURT4tbRCVUQkhPItLSMiIgEouIuI\nhJCCu4hICCm4i4iEkIK7iEgIKbiLiISQgruISAgpuIuIhND/A6dgFwIMbh0lAAAAAElFTkSuQmCC\n", "text/plain": [ "<Figure size 432x288 with 1 Axes>" ] }, "metadata": { "tags": [] }, "output_type": "display_data" } ], "source": [ "# Initial estimate of parameters\n", "theta0 = np.zeros((X.shape[1],1))+0.4\n", "#theta0 = np.array([[0],[0.5],[2],[0.5]])\n", "\n", "ypred = Xnorm.dot(theta0)\n", "\n", "sortidx = np.argsort(ynorm[:,0]) # sort the values for better visualization\n", "plt.plot(ynorm[sortidx,0],'o')\n", "plt.plot(ypred[sortidx,0],'--')" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "CT8J4nDxhC6t" }, "source": [ "Create a function `grad()` to compute the necessary gradients of the cost function. This is explained with higher details in the other post I mentioned in the beginning of this one. But only to remember, the gradient can be calculated as:\n", "\n", "$$\n", "\\frac{\\partial J}{\\partial \\theta_j} = \\frac{1}{m}\\sum(f(\\theta) - y)x_j\n", "$$\n", "\n", "Where $j=0,1,..,3$ since there are four predictors. This will produce a derivative for each input $j$, thus the complete gradient consists on the vector:\n", "\n", "$$\n", "\\nabla J = [\\frac{\\partial J}{\\partial \\theta_0},...,\\frac{\\partial J}{\\partial \\theta_3}]\n", "$$" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": {}, "colab_type": "code", "id": "e6yprDM1E1NT" }, "outputs": [], "source": [ "# calculate gradient\n", "def grad(theta):\n", " dJ = 1/m*np.sum((Xnorm.dot(theta)-ynorm)*Xnorm,axis=0).reshape(-1,1)\n", " return dJ" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 84 }, "colab_type": "code", "id": "SLJGzF32Gj1D", "outputId": "2f29ace3-b653-46b6-8772-06f06f926ff9" }, "outputs": [ { "data": { "text/plain": [ "array([[0.16830492],\n", " [0.0739841 ],\n", " [0.00197455],\n", " [0.00134667]])" ] }, "execution_count": 14, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "grad(theta0)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "n6sdJkuLiTE7" }, "source": [ "Similarly, calculate the cost function, also known as objective function, which can be expressed as the sum of the squared errors, as follows.\n", "\n", "$$\n", "J = \\sum(f(\\theta)-y)^2\n", "$$" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": {}, "colab_type": "code", "id": "Lf2S5qAwIpV1" }, "outputs": [], "source": [ "def cost(theta):\n", " J = np.sum((Xnorm.dot(theta)-ynorm)**2,axis=0)[0]\n", " return J" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 34 }, "colab_type": "code", "id": "ubJCF8eqIvti", "outputId": "8b16be10-ccf2-44b8-925a-7f522b0dcac0" }, "outputs": [ { "data": { "text/plain": [ "5.923489427685076" ] }, "execution_count": 16, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "cost(theta0)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "-M8pp59Wiqr8" }, "source": [ "We are ready to implement the Gradient Descent algorithm! The steps of this algorithm consists of:\n", "- Obtain the gradients of the cost function according the actual value of the parameters;\n", "- Calculate the cost to keep track of it;\n", "- Update the parameters according the following schedule:\n", "\n", "$$\n", "\\theta^{i+1} = \\theta^i - \\alpha\\times\\nabla J(\\theta^i)\n", "$$\n", "\n", "Where the superscript $i$ refers to the current iteration. Then the iterative steps are repeated until the algorithm converges." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": {}, "colab_type": "code", "id": "Oi7jm_89GlAf" }, "outputs": [], "source": [ "def GD(theta0,learning_rate = 0.5,epochs=1000,TOL=1e-7):\n", " \n", " theta_history = [theta0]\n", " J_history = [cost(theta0)]\n", " \n", " thetanew = theta0*10000\n", " print(f'epoch \\t Cost(J) \\t')\n", " for epoch in range(epochs):\n", " if epoch%100 == 0:\n", " print(f'{epoch:5d}\\t{J_history[-1]:7.4f}\\t')\n", " dJ = grad(theta0)\n", " J = cost(theta0)\n", " \n", " thetanew = theta0 - learning_rate*dJ\n", " theta_history.append(thetanew)\n", " J_history.append(J)\n", " \n", " if np.sum((thetanew - theta0)**2) < TOL:\n", " print('Convergence achieved.')\n", " break\n", " theta0 = thetanew\n", "\n", " return thetanew,theta_history,J_history\n", " \n", " \n", " " ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "H5xd0xF61E67" }, "source": [ "In the above code, a maximum number of iterations is fixed. This avoid the algorithm to loop infinitely. Additionally, the iteration process is stopped if the following criteria is met at any point.\n", "\n", "$$\n", "\\sum(\\theta^{i+1}-\\theta^{i})^2 < \\text{TOL}\n", "$$\n", "\n", "Where TOL is a tolerance, i.e a maximum difference between the values of the parameters between iterations so it can be stated that the values converged.\n", "\n", "Next, evaluate the Gradient Descent to determine the optimum set of parameters for the linear regression." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 118 }, "colab_type": "code", "id": "BzIc5-TVKZs6", "outputId": "7557efe1-48ec-4628-ecb3-3b32643fbaaa" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "epoch \t Cost(J) \t\n", " 0\t 5.9235\t\n", " 100\t 0.5097\t\n", " 200\t 0.3353\t\n", " 300\t 0.3226\t\n", "Convergence achieved.\n" ] } ], "source": [ "theta,theta_history,J_history = GD(theta0)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 286 }, "colab_type": "code", "id": "VNMFGZcKKf4P", "outputId": "00dcbeac-3c24-4a1f-f413-ebe5b7ac6026" }, "outputs": [ { "data": { "text/plain": [ "[<matplotlib.lines.Line2D at 0x7f4283fff2e8>]" ] }, "execution_count": 19, "metadata": { "tags": [] }, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGetJREFUeJzt3Xl0XOWd5vHvr6pUZW22LKu8L7Lx\nAjYYL8LYhLjJ5rB0FhJ6oBOyTHfGk5DJSc8kZ5ok032Scybp6XR30qEnp4lDQ9IJBBICE0IHAm0M\nJjQYZDDed4w32ZaN90VL6Z0/6sqWRVWpbFS675Wezzl16tatW6Xn3CM/fvXWrXvNOYeIiERHLOwA\nIiJyYVTcIiIRo+IWEYkYFbeISMSouEVEIkbFLSISMSpuEZGIUXGLiESMiltEJGISpXjTuro6V19f\nX4q3FhHpl1auXHnQOZcuZtuSFHd9fT2NjY2leGsRkX7JzN4sdltNlYiIRIyKW0QkYooqbjOrMbOH\nzWyjmW0wswWlDiYiIrkVO8f9A+BJ59wtZpYEKkqYSURECuixuM1sCLAQ+CyAc64VaC1tLBERyaeY\nqZKJQDNwn5m9Zmb3mFlliXOJiEgexRR3ApgD/LNzbjZwEriz+0ZmttjMGs2ssbm5uZdjiohIp2KK\nezew2zm3Inj8MNkiP49zbolzrsE515BOF3UM+dvctXQLdy3dwqETLRf1ehGRgaDHOW7n3D4z22Vm\n05xzm4D3AetLEeZ7T28GoLYyye3zJ5TiR4iIRF6xx3F/CbjfzFYDs4DvlCLMa3/1AQDaMh2leHsR\nkX6hqMMBnXOrgIYSZyEWMwA6dOF5EZG8vPrmZNDbdKi5RUTy8qq440FzZ5yKW0QkH6+KO2ZBcWvE\nLSKSl1fF3Tni1lSJiEh+fhW3aapERKQnXhV3TCNuEZEeeVXckJ0u0YhbRCQ//4rbDH3/RkQkP++K\nOxaDDo24RUTy8q6442aa4xYRKcC74o6Z5rhFRArxr7hjGnGLiBTiXXHrqBIRkcK8K+6YjioRESnI\nu+KOx/QFHBGRQvwrbn04KSJSkHfFrQ8nRUQK86649eGkiEhh/hW3mS5dJiJSgHfFrakSEZHCvCvu\n7EmmVNwiIvl4V9xmupCCiEgh3hV3XFMlIiIFeVncGnGLiOTnXXHHNMctIlKQd8Udj5kupCAiUkCi\nmI3MbAdwHMgA7c65hlIF0lElIiKFFVXcgfc45w6WLEkge+myUv8UEZHo8nOqRM0tIpJXscXtgKfM\nbKWZLS5pIJ0dUESkoGKnSq51zu0xs+HA02a20Tm3vOsGQaEvBhg/fvxFB9KIW0SksKJG3M65PcH9\nAeBRYF6ObZY45xqccw3pdPriA2nELSJSUI/FbWaVZlbduQwsAtaWLJAuXSYiUlAxUyUjgEfNrHP7\nB5xzT5YqkC5dJiJSWI/F7ZzbDlzZB1kAfeVdRKQn3h0OGDN9OCkiUoh3xa2vvIuIFOZfceuoEhGR\ngrwr7uyly8JOISLiL++KWyeZEhEpzLvijumoEhGRgvwrbtNx3CIihXhX3DqOW0SkMO+KW5cuExEp\nzLvi1tkBRUQK87O41dsiInl5V9w6rauISGHeFbfODigiUph/xa0Rt4hIQd4VdyxmOAdO5S0ikpN3\nxR3PXrBBhwSKiOThXXHHYkFxa8QtIpKTf8UdjLh1hkARkdy8K+54kEgjbhGR3Lwr7rMjbhW3iEhO\n3hV3PNY5VaLiFhHJxdvi1lElIiK5eVfcnVMlmuMWEcnNu+I+N1USchAREU/5V9wacYuIFORdccf0\n4aSISEFFF7eZxc3sNTN7vKSBsr2tDydFRPK4kBH3l4ENpQrSKa6vvIuIFFRUcZvZWOAm4J7Sxun6\nlXcVt4hILsWOuP8R+J9AyY/1OHtUiXpbRCSnHovbzP4YOOCcW9nDdovNrNHMGpubmy8+kE7rKiJS\nUDEj7ncBHzazHcCDwHvN7OfdN3LOLXHONTjnGtLp9EUHOjfiVnGLiOTSY3E7577mnBvrnKsHbgOe\ncc7dXqpAZ88OqBG3iEhO/h3HrS/giIgUlLiQjZ1zzwLPliRJQGcHFBEpzLsRt645KSJSmHfFbZoq\nEREpyLvi1tkBRUQK87C4s/c6HFBEJDfvijsZjwPQ0q4ht4hILt4Vd21VEoBDJ1pCTiIi4ifvirsu\nKO6DKm4RkZy8K+5UIs6Q8jKaj6u4RURy8a64AdLVKZo14hYRycnP4q5KacQtIpKHn8VdreIWEclH\nxS0iEjFeFnddVYqTrRlOtrSHHUVExDteFne6OgXA/mNnQk4iIuIfL4v70pHVAKzZczTkJCIi/vGy\nuC8bNZiqVIJXdrwVdhQREe94WdzxmDFnwlBeeeNw2FFERLzjZXEDzKsfyqb9x3XOEhGRbrwt7oVT\ns1eKX76lOeQkIiJ+8ba4Lx89hLqqJMs2qrhFRLrytrhjMWPh1DTLtzTr+pMiIl14W9wA75k2nCOn\n2li160jYUUREvOF1cS+ckiZm8OymA2FHERHxhtfFPaSijLkThrJMxS0icpbXxQ1w3bThrN1zjAPH\n9fV3ERGIRHFnDwt8bpOOLhERgQgU9/RRgxkxOMXSDZouERGBIorbzAaZ2ctm9rqZrTOzb/VFsC4/\nn/dfNoLnNjdzpi3Tlz9aRMRLxYy4W4D3OueuBGYB15vZ/NLGOt8HZ4zkdFuG57cc7MsfKyLipR6L\n22WdCB6WBbc+/UbM/EnDqB6U4Pfr9vXljxUR8VJRc9xmFjezVcAB4Gnn3Ioc2yw2s0Yza2xu7t0P\nEpOJGO+7dDhLN+ynPdPRq+8tIhI1RRW3cy7jnJsFjAXmmdnlObZZ4pxrcM41pNPp3s7JohkjOXyq\njZd1jm4RGeAu6KgS59wRYBlwfWni5PdHU9MkEzGeWre/r3+0iIhXijmqJG1mNcFyOfABYGOpg3VX\nmUqwcEodT6/fj3M66ZSIDFzFjLhHAcvMbDXwCtk57sdLGyu3RTNGsufIaV2LUkQGtERPGzjnVgOz\n+yBLjxZNH8HXY8bjq5uYObYm7DgiIqHw/puTXdVUJFk4Nc3jr++lQ+foFpEBKlLFDfDhK0ez9+gZ\nVu7UhYRFZGCKXHG/f/oIUokYv319b9hRRERCEbnirkoleP9lI/jdmiZ9GUdEBqTIFTfAh64cxcET\nrby4/VDYUURE+lwki/u6acOpSiU0XSIiA1Iki3tQWZxFM0bwxJp9OtWriAw4kSxugFvmjOV4S7vO\nGCgiA05ki3v+pGGMHVrOwyt3hx1FRKRPRba4YzHj43PG8oetB9lz5HTYcURE+kxkixvglrljcQ4e\nfVWjbhEZOCJd3ONqK5g/qZaHV+7WGQNFZMCIdHED/Mnccew4dIpXdugr8CIyMES+uG+4YiRVqQQP\nvrwz7CgiIn0i8sVdkUxw8+wxPL6micMnW8OOIyJScpEvboDb50+gtb2DX63cFXYUEZGS6xfFPW1k\nNfPqa7l/xU6dp1tE+r1+UdwAn5w/njcPneL5rQfDjiIiUlL9privv3wkdVVJfv7Sm2FHEREpqX5T\n3KlEnP/UMI6lG/az+/CpsOOIiJRMvylugE8tmEDMjPte2BF2FBGRkulXxT1qSDl/PHMUD72yi2Nn\n2sKOIyJSEv2quAE+9+5JnGhp1xdyRKTf6nfFffmYISyYNIz7XthBm65JKSL9UL8rboD/snAiTUfP\n8G+rm8KOIiLS6/plcV83dTiTh1dx93Pb9IUcEel3eixuMxtnZsvMbL2ZrTOzL/dFsHciFjPuuO4S\nNu47zr9v2B92HBGRXlXMiLsd+IpzbjowH/iimU0vbax37sNXjmbCsAruemaLztUtIv1Kj8XtnGty\nzr0aLB8HNgBjSh3snUrEY3zxPZNZu+cYz25qDjuOiEivuaA5bjOrB2YDK0oRprfdPHsMY4eW84Ol\nGnWLSP9RdHGbWRXwa+AvnHPHcjy/2MwazayxudmPEW5ZMOpetesIz232I5OIyDtVVHGbWRnZ0r7f\nOfdIrm2cc0uccw3OuYZ0Ot2bGd+Rj88Zy7jacr775CYdYSIi/UIxR5UY8C/ABufc90ofqXclEzG+\numga65uO8dvVe8OOIyLyjhUz4n4X8CngvWa2KrjdWOJcvepDM0czfdRg/v6pTbS269uUIhJtxRxV\n8gfnnDnnZjrnZgW33/VFuN4Sixl/ecOl7HrrNA+s0Pm6RSTa+uU3J3NZOKWOBZOGcdczWzl6WmcO\nFJHoGjDFbWZ846bLOHyqlbuWbgk7jojIRRswxQ3ZMwfedtV4fvofO9iy/3jYcURELsqAKm6Ary6a\nSkUyzrd+u15fyhGRSBpwxT2sKsX/+MBU/rD1IL9fpxNQiUj0DLjiBrh9/gQuHVnNNx9bx3Fd4kxE\nImZAFnciHuNvPnYF+4+f4e9+vynsOCIiF2RAFjfA7PFD+ew19fzspTdp3PFW2HFERIo2YIsb4KuL\npjF6SDl3PrKGlvZM2HFERIoyoIu7MpXg2zdfztYDJ/iHpzaHHUdEpCgDurgBrps2nE9cPZ4fP7+d\nF7cdCjuOiEiPBnxxA/yvmy6jflglX/nlKn0dXkS8p+IGKpIJvn/rLPYfb+Gv/t9afTFHRLym4g7M\nGlfDf3//FB57fS8PvLwz7DgiInmpuLu447rJLJya5luPrWftnqNhxxERyUnF3UUsZvzjrbMYVpXk\nC/ev5OgpzXeLiH9U3N3UVib54Sfn0HTkDF968DXaM7pijoj4RcWdw5zxQ/nfH72c5Zub+fbvNoQd\nR0TkPImwA/jqtnnj2bz/BPe+8AaTh1fxyasnhB1JRARQcRf09RsvZVvzCf76N+sYPaSc91w6POxI\nIiKaKikkEY/xfz8xm8tGVfOF+1ey8s3DYUcSEVFx96R6UBn3fXYeIwYP4s9+8gqbdckzEQmZirsI\n6eoUP/uzq0kmYnzynhVsPaDyFpHwqLiLNH5YBQ987mqcg9uWrNDFhkUkNCruCzBlRDUPLp6PGfzp\nj1/StImIhELFfYEmD6/iwcXziZnxp0teYvXuI2FHEpEBpsfiNrN7zeyAma3ti0BRcEm6iof+6wLK\nk3Fu/dFLLNt4IOxIIjKAFDPi/glwfYlzRM7EukoeueMaLhleyef+tZFf6IyCItJHeixu59xyQFfT\nzWF49SAeWryAayfX8bVH1vA3T2wg06FzeYtIaWmO+x2qTCW45zMNfOLq8fzoue185t6Xeetka9ix\nRKQf67XiNrPFZtZoZo3Nzc299baRUBaP8Z2br+BvP34FL+94iw/90x94fZc+tBSR0ui14nbOLXHO\nNTjnGtLpdG+9baTcetV4Hv78AgD+5O4X+fHy7XRo6kREepmmSnrZzLE1PP6la7luWppv/24Dn7jn\nJfYcOR12LBHpR4o5HPAXwIvANDPbbWZ/XvpY0Ta0MsmPPjWX794ykzW7j3L995fzy8ZdugixiPQK\nK0WZNDQ0uMbGxl5/3yjaeegUX/nVKl7ZcZh5E2v5zs2XM3l4ddixRMQzZrbSOddQzLaaKimx8cMq\neGjxAv7Px65g077j3PCD5/m732/kZEt72NFEJKJU3H0gFjNumzeepV/5Iz40czQ/XLaN6/7+WR5Y\nsVPXtBSRC6bi7kN1VSm+d+ssHrnjGibUVvD1R9dw/Q+e58m1TTr6RESKpuIOwZzxQ/nV5xdw9+1z\n6ehwfP7nr3LjXc/z+Oq9+ualiPRIH06GrD3TweOrm/inZ7awrfkkl6Qr+c/vmsjH5oyhIqlLgooM\nFBfy4aSK2xOZDscTa5u4+7ltrN1zjMGDEtx61Tg+vaCecbUVYccTkRJTcUeYc45Xdx7mvhd28MTa\nfXQ4x7unpPn4nDF8cMZIBpXFw44oIiVwIcWtv8U9Y2bMnVDL3Am1NB09zS9W7OTXr+7hyw+uojqV\n4KaZo/jo7DFcVV9LPGZhxxWREGjEHQEdHY6X3jjEr1fu4Ym1TZxqzTCsMsmiGSNYNGMk11wyjFRC\nI3GRKNNUST92sqWdZzc18+S6fTyzYT8nWzNUpRJcc8kw3j01zcIpdUwYVhl2TBG5QJoq6ccqg+mS\nm2aO4kxbhv/YdpCn1x9g+eZmnlq/H4DxtRVcO6WOqyfWMnfCUMbUlGOmaRWR/kIj7n7COceOQ6d4\nfkszz285yIvbDnEi+Fr9yMGDmFs/lLnjh3LluBouG1WtQw1FPKMR9wBkZkysq2RiXSWfXlBPe6aD\njfuOs/LNw2dv/7a6Kdg2e83MGaOHMH3UYGaMHsy0kdUMr05pZC4SARpxDyBNR0+zZvdR1jcdY93e\nY6zfe+y8c4VXpRJMrKtkUrqSSXVV2ft0JfXDKqlM6f94kVLSiFtyGjWknFFDylk0Y+TZdUdOtbJ+\n7zG2HDjB9uYTbD94ksYdh/nNqr3nvbamoowxNeXZ29Ds/dih5YwcUk66OkVdVVJHtoj0ERX3AFdT\nkeSayXVcM7nuvPWnWzO8cfAk2w+eYOdbp9hz+DR7jpxmx6GTvLD1ICdbM297r8GDEqSrU8FtEOmq\nFHXVSYZWJKkpL2NIeRlDKsqoqUgypLyMymRcUzMiF0HFLTmVJ+NMHz2Y6aMHv+055xxHT7ex+/Bp\nDhw/Q/PxlnO3E9n7NbuPcPBE69kPSHNJxIyaijIGl5cxeFAZVakElak4lckElanglowHy8F98Nyg\nshipRJxUIsagsux9KlinLyZJf6filgtmZtRUJKmpSAJDCm57qrWdI6faOHq6LbhvPbt8JLg/djr7\n/MnWdg4cP8PJlgwnW9s51ZKh9SLOV14WN1KJ+HnlngrKvSxuJGIxEnGjLB4jEQvug/Vlceu2HKMs\nlr1PxI2yWOe2hpkRjxlxM8zILneuNyMeo8vyuW1ilr1ll7Pna48H62Kx87eJWfY9jOyHykb2feh8\nHDwXCzIYQJftjG6vtzzru7xvLN82+uvIGypuKamKZIKKZILRNeUX9frW9g5OtbZzoqWdU62Z7H1L\nhjNtGVraO2hpz96ffdzWwZn2DC1t2efOtJ2/TXvG0d7Rwem27H17xtGW6aC9w5233JbpOLttW0an\n2u2qs9TPFbwF/1lw7j+V7H8hXR53vtbOe0y+53t4nXV7g7dvX1wOum9fxOsKZa+tSPLLzy+g1FTc\n4rVkIkYy0Tm6D4dzjkyH61boLrveOTpc9rQEmQ5Hh8veMh0E953rsmeA7HyvjHO4YF12OfuaTI73\ndWS3dQ5ckMcBOM49R/bndS4TbJN93bnljuAoMtfttV0fQ+fPfft2570v2YydywSvdV32W+fPOn99\n5+Pzn6f763rYvvvzdH++yNe5bsHflr/L9vme61yoHtQ3lariFumBWTB9EkdnZxQv6Ao4IiIRo+IW\nEYkYFbeISMSouEVEIkbFLSISMSpuEZGIUXGLiESMiltEJGJKcj5uM2sG3rzIl9cBB3sxTilFKSso\nbylFKStEK2+UssLF553gnEsXs2FJivudMLPGYk8mHrYoZQXlLaUoZYVo5Y1SVuibvJoqERGJGBW3\niEjE+FjcS8IOcAGilBWUt5SilBWilTdKWaEP8no3xy0iIoX5OOIWEZECvCluM7vezDaZ2VYzuzPs\nPLmY2Q4zW2Nmq8ysMVhXa2ZPm9mW4H5oiPnuNbMDZra2y7qc+SzrrmB/rzazOR5k/aaZ7Qn27yoz\nu7HLc18Lsm4ysw/2Zdbg548zs2Vmtt7M1pnZl4P13u3fAlm93L9mNsjMXjaz14O83wrWTzSzFUGu\nh8wsGaxPBY+3Bs/Xe5D1J2b2Rpd9OytYX5rfAxdcfSPMGxAHtgGTgCTwOjA97Fw5cu4A6rqt+y5w\nZ7B8J/C3IeZbCMwB1vaUD7gReILslZfmAys8yPpN4Ks5tp0e/E6kgInB70q8j/OOAuYEy9XA5iCX\nd/u3QFYv92+wj6qC5TJgRbDPfgncFqy/G/hCsHwHcHewfBvwkAdZfwLckmP7kvwe+DLingdsdc5t\nd861Ag8CHwk5U7E+Avw0WP4p8NGwgjjnlgNvdVudL99HgH91WS8BNWY2qm+S5s2az0eAB51zLc65\nN4CtZH9n+oxzrsk592qwfBzYAIzBw/1bIGs+oe7fYB+dCB6WBTcHvBd4OFjffd927vOHgfdZH13J\nuEDWfErye+BLcY8BdnV5vJvCv2hhccBTZrbSzBYH60Y455qC5X3AiHCi5ZUvn6/7/L8Ff1Le22Xa\nyauswZ/ms8mOtrzev92ygqf718ziZrYKOAA8TXbUf8Q5154j09m8wfNHgWFhZXXOde7bbwf79vtm\nluqeNdAr+9aX4o6Ka51zc4AbgC+a2cKuT7rs30beHqbjez7gn4FLgFlAE/AP4cZ5OzOrAn4N/IVz\n7ljX53zbvzmyert/nXMZ59wsYCzZ0f6lIUfKq3tWM7sc+BrZzFcBtcBfljKDL8W9BxjX5fHYYJ1X\nnHN7gvsDwKNkf8H2d/7pE9wfCC9hTvnyebfPnXP7g38UHcCPOffnuhdZzayMbBHe75x7JFjt5f7N\nldX3/QvgnDsCLAMWkJ1W6LygeddMZ/MGzw8BDvVx1K5Zrw+mp5xzrgW4jxLvW1+K+xVgSvApcpLs\nBw6PhZzpPGZWaWbVncvAImAt2ZyfCTb7DPCbcBLmlS/fY8Cng0+95wNHu/zJH4puc383k92/kM16\nW3A0wURgCvByH2cz4F+ADc6573V5yrv9my+rr/vXzNJmVhMslwMfIDsvvwy4Jdis+77t3Oe3AM8E\nf+2ElXVjl/+8jexcfNd92/u/B6X8BPZCbmQ/fd1Mdm7rG2HnyZFvEtlP3l8H1nVmJDu3thTYAvw7\nUBtixl+Q/RO4jexc2p/ny0f2U+4fBvt7DdDgQdafBVlWB7/wo7ps/40g6ybghhD27bVkp0FWA6uC\n240+7t8CWb3cv8BM4LUg11rgr4P1k8j+B7IV+BWQCtYPCh5vDZ6f5EHWZ4J9uxb4OeeOPCnJ74G+\nOSkiEjG+TJWIiEiRVNwiIhGj4hYRiRgVt4hIxKi4RUQiRsUtIhIxKm4RkYhRcYuIRMz/B4ok9iPV\nm1dkAAAAAElFTkSuQmCC\n", "text/plain": [ "<Figure size 432x288 with 1 Axes>" ] }, "metadata": { "tags": [] }, "output_type": "display_data" } ], "source": [ "plt.plot(J_history)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "H5dy43rk1_hn" }, "source": [ "Observe in the plot above how the cost function J drastically reduces at the initial iterations, converging to a value much smaller than the initial one. By using an appropriate Tolerance TOL, the iteration process was halted at less than 350 iterations, though the maximum number was initially fixed to 1000.\n", "\n", "We can perform predictions on the training set using the following code." ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 286 }, "colab_type": "code", "id": "4Ax6CZ0pK3sq", "outputId": "39965f10-6cb2-49a7-d647-a8422eea92dd" }, "outputs": [ { "data": { "text/plain": [ "[<matplotlib.lines.Line2D at 0x7f4283fb4908>]" ] }, "execution_count": 20, "metadata": { "tags": [] }, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd8VFXawPHfM5NJoySUUBJKkCqC\ngsSC2MCCWBF7W13dxXV1dV2XFXbdrmvBsuvua8G+LmtZxaCCYgEFUZASpIdekgAJkEpCMpk57x9n\nJjOTOiEJpDzfzyd7Z+7ce+eMYZ85Ofc5zxFjDEoppVo+x7FugFJKqcahAV0ppVoJDehKKdVKaEBX\nSqlWQgO6Ukq1EhrQlVKqldCArpRSrYQGdKWUaiU0oCulVCsREc5BIhIPvAwMAwxwO5AOvAMkAzuA\na40xubVdp2vXriY5OfnIW6uUUm3QihUr9htjEuo6TsKZ+i8ibwCLjDEvi0gkEAv8FjhojHlMRKYC\nnYwxD9Z2nZSUFLN8+fLwPoFSSikARGSFMSalruPqHHIRkTjgbOAVAGNMmTEmD7gCeMN32BvAxCNv\nrlJKqYYKZwy9H5ADvCYiaSLysoi0A7obY/b4jtkLdK/uZBGZLCLLRWR5Tk5O47RaKaVUFeEE9Ajg\nZOB5Y8xI4BAwNfgAY8dtqh27McbMMMakGGNSEhLqHAJSSil1hMIJ6BlAhjFmqe/5e9gAv09EegL4\nttlN00SllFLhqDOgG2P2ArtFZLBv13nAeuBD4FbfvluB2U3SQqWUUmEJK20R+AUw05fhsg34MfbL\n4F0RuQPYCVzbNE1USqmWKzUtk+nz0snKKyExPoYp4wczcWRSk7xXWAHdGLMKqC5l5rzGbY5SSrUe\nqWmZTJu1hhK3B4DMvBKmzVoD0CRBXWeKKqVUE5k+L70imPuVuD1Mn5feJO+nAV0ppZpIVl5JvfY3\nlAZ0pZRqIonxMfXa31Aa0JVSqolMGT+YGJczZF+My8mU8YNrOKNhws1yUUopVU/+G5/NKstFKaXU\nkZk4MqnJAnhlOuSilFKthAZ0pZRqJXTIRSmlGtHRnBlamQZ0pZRqJEd7ZmhlOuSilFKN5GjPDK1M\nA7pSSjWSoz0ztDIdclFKqSMQPFYeF+NCpIZVfmi6maGVaUBXSqkw+YN4Zl4JQiCA55W4Ge1YxwOR\n/+P2sl9TQPuKc5pyZmhlGtCVUqoWNQXxyr3xm52fk+LYxL0RH/Bw+S0AJDXHLBcR2QEUAh6g3BiT\nIiJ/An6KXUAa4LfGmLlN0UillDqqjAERZq/YybTUDRU3OmsaUgHoSDEAq7wDABBg8dRxTdzQUPXp\noY81xuyvtO8ZY8yTjdkgpZQ6aooPQuYK6D8OHEFFtL6fAZ/8hrG0o8T9UliXGuLYzbvl5/CxdzRw\n9MbNg2mWi1Kq7Vr1X5h5NezfFLp/9TsAdOQQHSmq8zKCl3+UTyLVO4ZI3NwdOYdHUw41RYtrFW5A\nN8BnIrJCRCYH7b9HRFaLyKsi0qkJ2qeUUk3DGFjxOvQ6FRDYs9ruP7DV9tr7ngnAAMmq8RLi28bF\nRDEn6mK+8w6jV1wUd8d+ydlbnwRPeZN+hMrCHXI50xiTKSLdgM9FZCPwPPBXbLD/K/AUcHvlE31f\nAJMB+vTp0yiNVkqpBtv1HRzYDJf/E/4zCboOhB/NhrXvYxDu3nMRz/ENAx2ZrPQMqjjNf2M05Ibn\nga1gvNBlAIjA6iKY9ROYeRVc9Sq063JUPlK4i0Rn+rbZIvIBcKoxZqH/dRF5Cfi4hnNnADMAUlJS\narunoJRSR8+K1yGqIwy7Cg7lwJd/gewNFC77LxvM8XxakExJVCQDJLP6IB5s8d9h4xyYstU+P/Ea\nKD8Mcx6AGefAdf+BxBFN/pHqDOgi0g5wGGMKfY8vBP4iIj2NMXt8h10JrG3CdiqlVOPxlEPGcjjx\nWohsB6N+jGfB43zy4m9ZfvgsMkwCXhw8X345G03vimBeY9bKvvXQbajtnfudfAt0PwHeux3czWem\naHfgA7ENjQD+a4z5VETeFJER2C+uHcCdTdZKpZRqTM4IuGcZHy/fzKOPzSczr4S/RYxhknMRv/f8\nk1w6AvCsZ1LFKTVO3/d6IWcjjLip6mtJJ8M9y8DpaopPUUWdAd0Ysw04qZr9tzRJi5RSqikZA95y\nUldnM+2jHRU55q95LuJC53IGOzJY4h0K2OyVXpLDXtOFbvEdqr9e/i4oK4LuQ6t//SgFc9C0RaVU\nG/PFokXs/+tA3nj3fyGVETebXjzqvpH9pmPFvgsdy1kUdT8nuTJqnr6/b73ddjuhKZsdFp36r5Rq\nM1LTMvnk88Wc78xFqpn3+b737JDnW4y9+TntFAejapq+3/s0uG5mzT30o0h76EqplmH397D5iwZd\nYvq8dDp7bbWSLFN7KmGMy8m9V18IDhejYvfVfGC7LnD8pfbm6jGmAV0p1TIsehq++GODLpGVV0JP\nOUC5cZBN1bmQ/hyVpPgYHp00nCtGJUOX/pBTywIVq96CvWsa1K7GokMuSqnmb+F0yM+A/N0Nukxi\nfAyJhw6STTzeSv3ZGnPMEwbD3hqyssvL4MN74Ix7ocfwBrWtMWhAV0o1bwe2wvyHwRUL7mI4XADR\nHes+rxpTxg9m4wfJ5HkCwyMxLiePThpec4nbU34KJbk1tG0zeMttDnozoEMuSqnmbdkr4IiAcQ/Z\n5wWZR3ypiSOTGHLlVF5rPxkhMLRSa73yfmfB0Mur7ncfhp3f2sfN4IYoaA9dKdWclR2CtP/A0Csg\nKcXuy8+Abscf2fWMYeJJPeu34ITHDRnLoEMPECdEREOH7rDjG5j7a3BGQZeBR9aeRqYBXSnVfK1+\nF0rz4dTJ0GMY3LnQFsA6UiW58OQguOQpGHVreOd4y+H1S2zN9MwVNk3xxneg2xC47Fl70zQi8sjb\n1Ig0oCulmi8RGHC+DaIi0LPKpPX6yc8Arxti4sM/xxUDnZJhyxd2rHz83+z+uF7hfykcJRrQlVLN\n16jb7I/f2lkgDjhhYq2n+dcBzcorIS7GhQjkFbu5psNangDoWM81PlNuh4Pb4IK/QlT7uo8/RjSg\nK6Wan+Wv2sAdHMwBvn/J9tSrCeg1LeacV+KuOCayeA+44NNdTi7qVY/2nPGL+n6CY0KzXJRSzcuu\nJTB3CmycaysZBovrVW0uempaJtNmrSHTVxGxpoUXEuUAbuPkka8rL4/cOmhAV0o1H16PrR8e3xcm\nzQBHpRAV1wsKsuxxQabPSw8ptFWTVd7+vOy5mIz8ssZsdbOhQy5KqeajaJ/NM7/k6epvXMb1slkn\nRfugY2LF7hprlVfymfcUPvOeQlJ8TGO1uFnRHrpSqnHl7oT0T47s3PwMu42vYf3huN52m2eHXVLT\nMhnz2Pwah1gq60I+sS6puRRuCxdWD11EdgCFgAcoN8akiEhn4B0gGbti0bXGmBrmxyql2oyXz4dD\n2fD7/fVf3KFdVzjrAUgYUv3r/c6CKdsgtnPFuHlNQy3+G6PxviyX/OJSvov+BTsG3MqgkRfXr10t\nRH2GXMYaY4LvJEwFvjTGPCYiU33PH2zU1imlWp6yIrvN22Un3dRH5+PgvD/U/LorhtS1B5k+b0HF\nDdDqVFtoqygbnixn0MAavixagYaMoV8BnOt7/AbwFRrQlVI3z4LXLoKD2+sf0Av32qn1MfHV5pLn\nFrv5mfMjTjGdyOTMai8hUP1izv4aMEFj761NuAHdAJ+JiAFeNMbMALobY/b4Xt+LXUxaKdXa5Gyy\niyp3Pi684/3HHdwKnF/n4cGB+9XoZ0hiHxcefqzGXPIrnIvJMAmkeqsP6Ik13fDM9wX0uHpOKmpB\nwg3oZxpjMkWkG/C5iGwMftEYY3zBvgoRmQxMBujTp4YbHUqp5itjGcz/K0REwXFjYcx90Llf9cce\nOgDv3AwTpkPKj+u8dOVx8ARvDrt9KwnVdKMz03QhSarPI491CVPP6139iRU99PrMKGpZwgroxphM\n3zZbRD4ATgX2iUhPY8weEekJZNdw7gxgBkBKSkq4N6OVUs3F4Am2DvnW+fDDW3Zs/JZZ1R+btwMy\nvocz7w+5IRrcC0+Mj2HskAQWbMypMg6eKPtJ89ZefCvLdCXFsanK/qT4GF7p+xlD5twAw3dXrZne\n6xQY+zt747WVqjOgi0g7wGGMKfQ9vhD4C/AhcCvwmG87uykbqpQ6Bg4dsBkro26DU38Kn0yFFa+B\nu8QWrark+7RVnAo8/uZs4l3v84LjOnKL3SHDJ5l5Jfxnya4q58ZwmM5SRJapPeBmmS7EyyHaUcIh\nYkIXqPjz1fag9Llw0vWhJyadbH9asXB66N2BD0TEf/x/jTGfisgy4F0RuQPYCVzbdM1USh0T6XPg\nw1/AfT/YioMnTLSlYoMCenANlTudyzjVBT3kADea+TxRfAXgDCtPPFEOAHZIpTZZpiulxkWC5BEf\n1zmQzVJaaA8QJ/Q5veqJOekQ26Vt99CNMduAKjUrjTEHgPOaolFKqWYid4cNkP5x5z6nhwTLymPg\nibKffBPLWtMPl3hIkv3sMjZfYoBkIMBmk0RgOeaAg6YDv3PfzkoTWCyici55XrGbVR3O5ZPxd/HV\nyZXGwnd8A8YDt35kv3wq+8/V0He0LSnQSunUf6Xakqw0WPu+XZez7JC90XnJ0+CKrv743B0Q39tm\nufiVl8KeH6D3qVVqqBwwcSzyDmeHtwcAybK3IqDfE5HKGMc6Lit9mL1U7YXn0pE5kRMQASl2k1jT\nos01GXAB3P6ZrZm+cY5dg7T/WPua1wOFWfUvm9vCaEBXqi0oL4P3fgwbP7ZLpsV0sj9Xv1pzMAcb\n0Cv3dpe+AJ//AX61sUoNlWc9k8ADCdhJ48myl4W+P/BPls0kSD6vRz7BRWWPh5wX43Ly7IUduGBQ\np/DW5/zwXrvoxcibAvucEdDnNDAGvvgTtOsWCOiHcmwNmFacsghay0WptiEi0vZYz/0tTNkCv06H\nu5fUHTxzd0Anm6Lor5ty8cc2e+UPzzxb49h4DvEUmWi6Sj4ACeTRx5FDjoljiGM3d42KJSk+JmSh\n5gsOzISZV4f3ebbOh+0LA8/zM+DTaXYykwgMuwp2LraVGQH2b7bbVt5D14CuVEt26IAdPqmN8YXd\nq16Ccx8MTefLWgWvXxooilXpvCXD/8Ida4aSPHUO97+zisy8EjaY3uSYOEaVrwo5PJ5Cvo26h0sd\n3wHCWHmV1yJvQoDzO+wEIOHSPwLw4IAsFk8dx/bHLmHx1HF2WCV/t62mGI64XnBgc+CzbfkSljwH\n5Yft82FXAQbWfWCff/MMRMdDoma5KKWaq+nHQXQcTK2aBljh1YvsAsuXPFX1tagOsGORDXy+VXlC\nV/7phKETEEg7NDhY6B3OuY4fELwYX7+wl+SQKAe5ZcxA/nXJJaHv89kSWOKCETfCV4/ZHnbwcAnY\nL5XEkeF97kEXwRd/hJVv2JTKbQugQ89AUa+uA6HHcFj5Joy+Gy7/J3jKoEPrntCuPXSlWip/7/Rw\nPpQWVX+M1wN7Vtn6KNXp0h96joA17wGhK//0ln2c5fiBKKouBrHIcyJdpJBhsqNiXy/f7M3TTh5h\nd2z6DN6+CTzlMOaX8KPZNtWx/zgbgINXI/J67dT8cHvoZ/zCXmfub2DfOtj2lZ3FKkHZMydebydE\neXxj5zXNbm1FNKAr1VIVHwg8Prit+mMObrfDEN1qGSsffrUN+ge2hmStXORYxr8jHycKd5VT5ntH\ncFXpH1lv+lbsGxqbZx/4a5kX7bU3YQsyoF0XSB5j95/+M5j0EiGT+4v3g6c0UO+8Lg4nTHoZzp5i\ns3VKcm2AD3baz+xNX4czvGu2AhrQlWqpjBdO+Snc8Tn0PLH6Y7LX2W1tNz9PmATAS89PD5mK30ey\nyTPtKKBdlVMKaM8KMxgPNljGuJxc0tsNUR3tWDVAZ1+lxU3zYNHTUJRjnyeOhAHnhQbayPZwwzsw\n8IK6P7dfuy5wzhS7elG7BDju3NDXnRHQKyW0197K6Ri6Ui1V+25wyZP2sX/4pXLwyt4ACHSteYWe\n1G1Q4r2AtcWdQvb3kWx2mW4h+4In+vQkhxvcHzA79mpumXAW/T1Z0CU60AZ/1cXvZ8CBLXDidYEL\nZaXZG7L+Al6RsTD4ovA/e7DjL4Mhl7apwF0TDehKtVSH8yEixk7yefcWuP6/VWuVdDseTrvTBswa\nTJ+XTmZZ1cqIvSWb9aZvRRCvsmhE3i74+8/50Zmnw8gk4PbQC3ToYVMlD2yxM02Dc8DXpcJ3/7LD\nPVEdYO9am2I44PyqC0OHQ4M5oEMuSrVcn/8Rnhlqe8JF++zQRmVDr4AJj1fdH8Q/Oag9xQwTOxbv\nwEsvySE3MolnrhvBjuD0Qr/4PtDnDFj9jv0LwVMeemERSBplH/c+JfS1gRfYiT6r3rLPV/0X/neb\nBuYG0h66Ui1V3k6bFdKuC/Q6FTZ9AmOnBV73lMPhvIpiVNWtAJRX7MYhgscYnnK9wHDHNsaU2glD\nk6Om8/pd50GnWibjnHgNfHw/bP8a3pxkUyOD66BPmgFPH29ndQbrO8ZmpXz5Fzj+0kAOugb0BtEe\nulItVe7OQEbJoPF26MU/MxLsDdHp/WHDRyHpiAa7AlBusRsDeHzj7x97TidRDnKaYyPRLhcTJ1xU\nfZGrYEMngsMFC5+0hbFiK9Vo2b8JnJH2CyeYCFz6NHjdMHeKzUFv5dPyjwYN6Eq1RF6v7dXG+9IG\nB/luKG7+LHDMvvV223VQlSJa1ZlvUigy0dwUs4TnxwoTvV/YQly1ie0MJ15rh3wg8AXjd9y5MHU3\nJI6oem7n42Dc7+04v/+vDdUgGtCVaomK9tqZj518Ab3b8XDGvdB9WOCY7PV4HC7OfnlXlZWBqlNs\nImk/4koui1jGue6FMOdXtnRuXSY+Z2drQtWADrb4V0254GfcA2f92ubUh5uDrmqkY+hKtUQR0XD+\nn+1YNNghjAv/GnLIvi0rOehJZFdx1Zme1UmMj7G97R/egrT/2ADrDDNE5PlKD8R0qv246jhd8LPF\nEBNf/3NViLB76CLiFJE0EfnY9/x1EdkuIqt8P9X8TaWUahKxneHMX9qeebD9W+DrJ8AYJHsDG7zh\n9XpjXE6mjB8M/c6Bny+1wyF1jZ8Hi2xn65EfyU1Nh9PWmtEhlwarTw/9PmADELzy6hRjzHuN2ySl\nVJ3ydgFiF58ItuULWPAI93/jJMJ9FRkmodrTg1cAqrKQRLchtmzu0CvCb895fziij6EaV1gBXUR6\nAZcAjwC/atIWKaXqtuBRW+DqgY0huz+MGM9I8xQ/KX2DSz2PVFRCDJYUH8PiqeOq7K9wOB9KDtrV\njFSLEu6Qy9+B3wDeSvsfEZHVIvKMiFT72xeRySKyXESW5+TkNKStSim/vF2BDBcCi0/c+956pruv\n4QTHTqZGvAWVlqCoGFqpTXQcXP8WnDut9uNUs1NnQBeRS4FsY8yKSi9NA4YApwCdgQerO98YM8MY\nk2KMSUlIqP7PP6VUPeXtrMhwCc4xB/jIOxqAOyPmhJziXxkorDU6h1ysNylboHCGXMYAl4vIxUA0\n0FFE/mOMudn3eqmIvAb8uqkaqZQK4nFDQWZFD71yjrnBwVmlz9CTg9hyWmEMs6hWoc4eujFmmjGm\nlzEmGbgemG+MuVlEegKIiAATgbVN2lKllJW/G4yXlQUdGfPY/GpzzHeb7nxvbAZMWMMsqlVoSB76\nTBFJwHYBVgE/a5wmKaVqFduFZaOmM/X7aDLdtU8YqlIhUbVq9QroxpivgK98j/XvN6WOheg4frlu\nQK3BPMblDH+8XLUaOvVfqebqcEFg4Ypge1bTI/+HGk+r181P1apoQFeqOcpaBY/1hvWzq762+B88\nG/18taf5b35qMG+bNKAr1Rx1GWC3Gcuqvpa7g8guycS4Qgte6c1PpQFdqeYoqr1d+Li0IGR3alom\nBzM3M39vDFERDjrFuhB0mEVZWm1Rqaby9XTbw77p3fqd9+k0W6gqrjfk7a5YaSgzr4RYDrM+Op+d\npht5JW5iXE6euW6EBnIFaA9dqar8pWAbasHDsHkelIdXvhaAg9tg6Yt25aH43hTu2x4yCzRJ9gOw\n23QDoMTtYfq89MZpr2rxNKArFWz1u/D34bDz28a7ZllR+McuehocEXDGL2DUj5ledlXILNDdJoGr\nS//AYm9gIYusMBavUG2DDrkoFWznYrvduxb6ntGwa8X1huSzbO3ycORn2MUlRv2Y1C0eps8zZBaO\nCjnkMFEsN0NC9iXGxzSsnarV0B66UsH8S67VZ3GHmpxyBwyeAEXZ4R2/8k3wevgs/lqmzVrDgbw8\nTpZNdCJwY/RMxxrGO76veK6ZLSqYBnSlgsX3gT6jYdCFDb/WmffbHvebk8I7vt/ZMO53/HlRESVu\nD8myj1lRf2K0Y33FIbc5P+WXEbMAzWxRVemQi1LBzvyl/Wmo8jK78HHXQbBpHrhLwFX70EhqbjLT\nvy2tuAGaaboCgRuhAMdF7Ce6+yB2/OyShrdRtTraQ1eqsg9+BjOvadg1stfB00OgeD8YD+xdE3ht\n3zr4ZKoN8j4rZj/Hi7PmhVROLCSWAhNbEdCT4qI5zrmfxH6V1hFVykcDulLBXhpnh0lyGpgKeOiA\n3Q64wG6z0gKvzX8Ylj4Pn/zGPi/K4cS03zPJfF7lMpmmK0mynxiXk4fO7QLlJSErFSkVTAO6Un7u\nEsj0LcxVuAe8lVdcrIdiX0DvPgzad4fMlb7r7rVDMIkjWdDxCsY8Np9HHv0jLsp5x3NulctkmK4k\nRxzk0UnDmdDLl8/eGDdsVaukY+hK+RXusdvEkbZHXbwf2nc7smsV+8a923WBC/4KHRMBWDf3BU4w\nHsZuv5kd20oxGO6Kms0K70C2mF5VLvNOzHW8fPNIBvZJAm8P+OUaiO1yZG1SrV7YPXQRcYpImoh8\n7HveT0SWisgWEXlHRCKbrplKHQUFvoCelOJ7nnnk1zq036ZARsfDSddBv7NITcvko3X7mes5le2m\nJwb4S8TrdJYi3vaMrXKJGJeTSydcBn1OszscTpuFE9nuyNulWrX69NDvAzYAHX3PHweeMca8LSIv\nAHcA1df0VKol8PfQjzsHDudDRJgTdr78K3QfCsOuCuwbeKEdahHhw+Xb+WzeR6wsjCOLCcCEisNe\n81xECZF85BldsU+wk4WmjB/MxMExsPZ9O0Fpy5d2DD3l9oZ/VtUqhRXQRaQXcAnwCPAr3zqi44Ab\nfYe8AfwJDeiqJYuOs7ng/c6B4y8L7xxjYNGT9nFwQO87GvqOtoW1Zn/PYufveTfiHB4s/ykm6A/j\n7aYnj5bfVPG8ymLOGSvgvdvhhrch7U0wXg3oqkbhDrn8HfgN4L9L1AXIM8aU+55nADq7QbVsAy+A\nWz+CaN8foR533eeUl1a/P2cTFO61VRLd7SkwMVwb8TVTImquvFjtrM/43nabtwtyd+oNUVWrOnvo\nInIpkG2MWSEi59b3DURkMjAZoE+fPvVuoFLHxHOjocdwmDSj9uNc0XDKT2DNexVlbrPySpgf/Ws2\n05fMw78AwOnrC83xnBZyugCGWhZzbpcAzig4sMWO6WtAV7UIZ8hlDHC5iFwMRGPH0P8BxItIhK+X\n3guo9g6SMWYGMAMgJSWlmgUSlWpk8x+Bwiy44v/qd96bk6BdVxvEo+NsCVuf4GAdF+NCBPKK3SRG\nu7lKyjjP05kH31lGKS4A4kwB+zyBm5d3un/F2Y7VrDP9KvbVGMSDidja6Du/A4wGdFWrOgO6MWYa\nMA3A10P/tTHmJhH5H3A18DZwK1DN4odKHQOHsmHDx3D5v2xADNeBLYHKiB0TK3LHU9MymTZrTUUZ\n27ySwFDMdeUfcJfzI44ve41y3/+dHHiJ5xAHK/IH4BvvcL7xDgfs0Eq9arDE94ZtX4E4dFKRqlVD\n8tAfBN4WkYeBNOCVxmmSUg2w/kN7I/Fwnh137hRmADTGTvrp0NM+75iIZ/1HnP3ol2TmH67xtL6S\nzV7TuSKYA3SiEIcYDpgOVY4Pq1de2YTpEBFl2yY6F1DVrF7/OowxXxljLvU93maMOdUYM8AYc40x\npoa7Q0odRUtfhH2+uil7V4d/XvFB8JTy9++L6Dd1Do9/W4TTW0Zxvi19K3i5wzmH7hwMOa2v7COf\ndsyOfIgrHN8A0EkKAcitFND9GSz1ro6YMMh+MUVEglPnAqqa6de9aj1KC2H3Ejj1TtuT3VN3QE9N\ny2TMY/OZ8PA7AKQXd8AAi0uP47nyyzHYIZuRsoXfu2bysOu1kPP7yD7WepM5QXYwwGHH3HNMPPeX\n3cVK78CK4xpUtzxvFzyeDJ88eGTnqzZDA7pqfooPwjfP1L+Wyo5vwFsOQy6xeeR1TJH3j41n5pVw\nmEje95zFFmN7z6tNf54ov548bC/7XOcqADpLYLGJDhTTWYrYZnqyl84kyn4EyKc9C6LGURybiNAI\ndcvzdkFJLix7+cjOV22G/v2mmp+5v7azI5NSoN9Z4Z+35UtwxUKf0+1szzpMn5decaNzu+nJA+67\ngl41dKAEwUsB7TlZNgOw3DuIeF+WS3kxPMmt/BB5Ans8K+nnyuOZSSOYmFxuZ50mpTTOEEmcLxfd\nof93VbXTfyGq+UkaZQN6TZN2anIo206Rj4iyz42xMysdzmoPD15c2UU5bpzgG2JxYEiLmswLnst4\nsvw6fuSeRhKHeOCqs1gV0tP21U1/fxlkLOPkkUnw9XRY8DA8lE2j/F+sYxLEdoUL/9rwa6lWTQO6\nan4Gjod5vw1ULAzXtf8Gj2/yctYqeONyuOY1GHBetXnkwZMi/hzxGmc51nJW2T8A8OIgm070jchD\nyqFnfDseGH8yE0ckBlYfyt1hvzA6Hwe9TwOv7e1TvB+iOga+WBrKGQG/2do411KtmgZ01fyU5Npt\nuIsrB/MPccT3gdJ8nnv7A54oOlwxIxNC88j9ekgu+bQLmbkZFduHy+IMl3VdYMvonngWPHGcraVy\n3u9h4ZO2tvmUzXDqT+0P2Fro/nx2pY4iDeiq+XnrOjvdfdBF4Z/z/k8hIpLUPr+19VPySvgmqiu9\nSzcDE6hrinIPyeVgRALPXDmuz4M3AAAgAElEQVQicPPyf/+2Oe0ZK+Ck68HpshOO9tgbpOTusL3z\nyooPaM1ydUxolotqXkqLbEA8d6rNvw5D6soM8tfM5Z1lGdz/zqqKdTnXeZMZKjvDukZ3OciZIytl\nonRMgvxd4D4Eg8bbfYkj7OIXxsDBbdDZN5X/wFZ4chCsn21rocd2DfsjK9VYNKCr5iXPF4DdxbD7\n+zoPT03L5PEPviOOIjaZpJCe+FpvMv1kL7HUPNMTIBI3XaQwMEvUb8ilNrPE4bI3W8GuZlR8wFcs\nKws6+QJ6bBco2gd5u+HiJ+Gc34T5gZVqPDrkopqXvF12u+J1SP8U7vqm1sOnz0unW3kWOGGX6R7y\n2iLviUR7yojETTHRNV6jgwvSB93J4OQzQ1/oc7otjNV1EETG2n09R9rt+lTABHro0XEQ2QHyM+CM\ne8L8sEo1Lg3oqnnJ9fXQe50SWLC5Fll5JYxy2JunOysF9FVmAKvKB1Q899/wjA+ulhgfw5TxIxg8\ncmLVi5cfhr5n2kUv/LqfAGc9YHvs17xh2wmBqogHt8LaWdArxd6YVeoo0oCumpcB58Flz9qhl03z\n7GxRR9WRQX8aogFyac8XnpHsMlUXdG7HYbpIPp645JqLYmVvtLNTYzqFVmd0xcDESiV4XdFw3h+q\nb3tcL8hYDps/s5UeT76lHh9cqYbTgK6al64D7c+SF8B4oOSgrVGODeKzP5nL1wU9MDgqxssXeU9k\nkffEiksEpx7O6vicLaj18++qvlfxQZj/MCx/1Qbju78PDK3UpuwQLHvF9s77BtYCZfAE+85bvtAs\nF3VMaEBXzcu2r+xQRfsEAL5cvpY/fOclM6+EwbKbeVEP8mzERJ4uv7biFCcePNjZoFXK0369AhY8\nYoN3cG74rqXw1vW2zO5pd8K508IL5gCr34XPfw9RcTBtV2D/KXfYtm/5ouJLSKmjSQO6aj6MgXdu\nsTnfZ/2axaNf5jdf5HLAHQlAL7Fj5ac5NoactijqPj7xnMbD5beELrAM0PcMwMDupb4etM+S52xF\nxjsXQY9h9Wtnou/GaGl+1df8qxxpD10dAxrQVfNRkgulBTy7soxnFi7HIe3wmEAiYi+xpQDuLruv\nYl8UZSTKQfJMOxLjY6peM2kUOCNh5+LQgN7teEgYUv9gDtBtqN2eMCl0/64l8NG99rEGdHUMhLNI\ndDSwEIjyHf+eMeaPIvI6cA7g76bcZoxZ1VQNVa3fgqUrGAusK45H8HKhLGMHPdhg7KpDfSSbYhPF\n/qCl3fr4eu17nD2rrzfuirZBfWelMfRzpx55QyMi4VcbIKbS9P4OPex25M02jVGpoyycHnopMM4Y\nUyQiLuAbEfnE99oUY8x7Tdc81WLlZ0JUB4juGLq/JM9WUezQvcopX3z7PWOBDNMNL8I/XP/iFc/F\nbCi3Af1dzzkUEc2zrn/xqPtG9tKFvrIPgMvOOYNzaqo3ft4fQgtlFWRBuwQ7lf9IdUysuq9DIiC2\n3G191jJVqpHUOVPUWEW+py7fT12lMVRb9+p4+N+tVff/36nwVPVT+tsVZwCw2yQAwgHi6EpgnDrd\n9OErzwgud37H2R0yeea6Ebx8mR3aOOf002puS98zbC/d738/hn9Xk3feUBGRgLEZMEodA2FN/RcR\np4isArKBz40xS30vPSIiq0XkGRFppFqhqlVIHAn71ofuM8ZOj/c/rmRJ+/O4pWwqBbQDYL/pSFfJ\nxymCYLi5w0runnAyIDw+RmwmS/dhcPrP665uuHEObP4cDudDxjLoU8sXQEMdOoIqkUo1grBuihpj\nPMAIEYkHPhCRYcA0YC8QCcwAHgT+UvlcEZkMTAbo00dnzrV6pYU2RbDP6bDhQyjYAx19NVIObLHb\ny56tdkji9otGM21W+4q64vtNHN0cBTx17UlMHOiCJ2+CqO7QpT/s9S0Efdw5Ya1OxNdP2CGg0+60\n+e39z2uMT1vVA+k2e0apY6Be//KMMXnAAuAiY8we33BMKfAacGoN58wwxqQYY1ISEhIa3mLVvC19\nEf45ylYqBJsu6Od/3CvFBv5KJjoW8fxYm0suQHFkF46LKbY98dwd9qBOybZXvm+tfZ6fEVjUojZ9\nx9ieefqnENkeelf7z7XhOvSwtdOVOgbqDOgikuDrmSMiMcAFwEYR6enbJ8BEYG1TNlS1AIcL4Lt/\n2en7gy+GiOjQionGCz2Gw6sXwVePVexOTctkzKNfUjzrXvYsnsmU8YPZ/tglXHLPM8Tc8bE96OB2\nu+2UbMfDozralYP+McIu91aXvqNtbZZV/7G1WRpyQ1SpZiqcIZeewBsi4sR+AbxrjPlYROaLiL17\nBauAnzVhO1VL8P0Mm0t+zoP2BuGZ99t8b7+Tf2R/XjybfVvSmJQ2n8y8EgToTD6x0aVsPNyJd2fZ\n4ZSJI5MD5/p76PF9YMy99id3J3jdgRK2tenjm6Lf/zzbLqVaoToDujFmNTCymv3jqjlctVUFWfDt\nP+16oEkn233Bud5eL6mrspj+2SbuPxTH2Y7VZJbahSgM0FtyAJvhUuL2MH1eOhP7eezNzGFX2YDe\nIdHmlfvl+nrtncMI6O262olEMfFNN9yi1DGmd29Uzfatg1X/De/Ynd/azJULg4Y/jLH1zYuy+W7e\n24yZPYaY/M1s9Pamm+TRiYKKQ3v7Jgjt9lVMzMorsSsCfToVctJh3ENw/X8C1377Jnj7Zvs4nB46\nwG1z4SpNKVStl079b+sObIXCvZA8JrDP64GXxtpMEofLjjnH9ar1Mqnlo3mBf5L+1GbiYnYgAqY4\nl1XRk3lWbsTpLiLFWUiGSSDd9AZgiGM333lPAKCf7MVjhAxji1olxsdAO9/NxUPZ0O8siAuaOHQ4\nH8oKbfuqm+RTnXY6HV+1btpDb8vcJTDzavtzONBbJmMZ7PkBzv+zfb7g0ZqvsfNbvvv0v0ybtYaN\n+REYIK/ETW6xmzzas8WbyNDyDYxybGad6ctholjrTeYJ97UVwRvgX56JXFr2N0qIJsbltNP4/dki\n+Rl2OCcnPfC+PXzlcsf/DRzOxvnvoVQLpz30tmzhk3ZYA2DdLBh1m32c/oldS3PUrXYi0JLnYPTd\n0N1XlKqsGHZ9C1sXwA9vkVgSjdv9KNX9c1rhHcTFzqW4KGem53wAcunIc57QmZoGBxtM39Dyt14v\niNN+wWz4yPbYE3z1WvxFtcLJQVeqjdCA3lYd2AqL/w4nXg8DL4CBFwZe2zTPTpePjrPLra18E778\nC9z4tg2s790BnlI8jkhWcjy/PXwT5TX8U1phBnKdfGUfewdW7O9MAX0lmzQzgLvaL+LqPkX0v/nZ\n0NWJHA57M3P3Mvu8U3Lgte6+gL59YSDIK9XGaUBvqzofZ2dsDroodGw5dwfkbKhYPi01vYSdnquI\nWL+ftx/9nInHuUjiQuaWDWWZdwgl1F7xYYXX1m0pNDEs8wYC77R2H3El84n4XRa89ASURlS71Bw/\nnW8XlPjyz6HZLAmDbYGtmE5H/J9AqdZGA3pbJQIjbwo8X/G6XVpt8AR29pnElAXxfD97jm85NztU\nQn4Z/0wDuCHst9lqErmn7Bcs8Q5lP3HEuJw8Omk4E7058NEc2PkN7FkVGK+vLK4XHMoBV6wN4H4R\nUTBlSz0/tFKtm94UbYvevgm+/Vfovm1fw8LpfLgNLtp+Hd8XxANHVlYzPsZFp1gXAsTHRLI4+mwO\nEEdSfIwN5iOToJvNbuHrJ+z2+Muqv9jmz+0YfqdkLUmrVB20h97WHC6wk3W6DyM1LZPp89LJyivh\nwugTeJFZ7E79IyXl12EnANdfUnxM1WXgqtNtiN3uWGSDe5f+1R+3/Wu7vXnWEbVHqbZEe+htTeZy\nwLC4bADTZq0hM68EA3x22I5v3x3xIadIeq2XqElFumE4ojpAnK/65ogbaz7On4se2e6I2qRUW6I9\n9DYguCf+YPQH/BThzgVQgqfiGIODR9w38iPn56SZAWFf246xE5puGK4r/mVzzYPrvVTmL6K1bQEM\nvSL8ayvVBmlAb6X8Qdxf/Mo/Fj7Us5F06UMRsVXOeclzKS95Lq2yPzhojx2SwIKNOWTllZB4JEE8\nWDg55C5fOwv2HNl7KNWGaEBvhVLTMnl41lK6lu8DemOCxsO3mkS+9w4J+1pH1PNuTCNvtr30E687\nNu+vVAuiAb05y91hF1XueVLtGR6HDoTkkk+fl84ZnuU8G/V/vF5+Ic+VX0E2Nl/7z+XVrPNZjYr0\nwmMVyP0cztrH2JVSFfSmaHO27GW7GERt9q6F6cdBWqASYVZeCec4fwDgtojPmOj8BoBYDlNTImJw\nqmFIeqFSqsXQHnpzVpQN7RNgxzd22TZXTNVjMpfb7fcz7PAEkBQXxdmH1zDbcwZ9JJsrnYuZ4bmM\nR1yvMFgyuLgsUGyr2fTElVINFs4SdNEi8r2I/CAi60Tkz779/URkqYhsEZF3RCSy6ZvbxhTtg/xM\neONSWDWz2kNSHRfwrvMSPFmrOf9P7zDyL58RV5BOguSz0HMiH3jGcLxjF0NkFymyiUxHT+2JK9VK\nhdNDLwXGGWOKRMQFfCMinwC/Ap4xxrwtIi8AdwDPN2Fb256ibBg8wQb2xc/CybeRunpfRQpil2jI\ndwt9vOdybdQcznEv5BXPxVzvXA3AIu9wynHyB9ebzOj/Lb0zcug9/n7SRl9YxxsrpVqiOnvoxiry\nPXX5fgwwDnjPt/8N7ELRbceBrfDaxXZxiHDNvgdW/jv844v2QftuLEn8EeTt5N7f/4H731lVMRno\nac+jPOp4nq0miX+WTyTNa/PH/+25gBvLfks2nYiJ707EwPPpk/GhvWbv08N/f6VUixLWTVERcYrI\nKiAb+BzYCuQZY8p9h2QA1f7dLiKTRWS5iCzPyclpjDY3D+s+gJ2LYcHfwjveGEh7Ez78RfjvceUM\nFrS7mNu/68ombxJ3R8zGif1P3kuyOdu5hp1eO5PyqfJrWWlsZcNDxPCt15aXzcorgWvegNN/DhHR\n0GN4+O+vlGpRwgroxhiPMWYE0As4FQg7kdkYM8MYk2KMSUlISKj7hBbDly2y5QsbrOtSWlj/txh4\nPg8tdVLsNvyr/Eq6SAFxHALgWudXeI3wnicwOWew7OJe5yzucX5AO+wCzInxMRAZC4Mvhgv+ChF6\nq0Op1qpeWS7GmDwRWQCMBuJFJMLXS+8FZDZFA5utfudC4hxbBbC0wC4GURv/0Mz48Hr0c5es4esv\n51BY2A9oz4feM/iu9HgOEIeLcu6NSGWB5yT2EMg//13ETM52ruGQieJFz2WhtVX6nWV/lFKtVjhZ\nLgkiEu97HANcAGwAFgBX+w67FZjdVI1slnqfApO/gmterzuYg12c4RcrYcRNdR6ampbJ7Dkf87j7\nUY6TwBh9jm9y0HXOBQC85RmHyyEVWSvznLa3/p33BLrFd9AMFqXamHB66D2BN0TEif0CeNcY87GI\nrAfeFpGHgTTglSZsZ/OTvQFiu9jiUjmb7MrzUe1rPt7psut0Ln3RLiwR16vGQ6fPS+cM70FwQo6p\n+mXxluc80rwDyOt4PNMvGhII2qVnw4tzOX/cPZw/LIwStkqpVqXOgG6MWQ2MrGb/Nux4ets081ro\nOxrOuBdeGAOX/xNO/lHNx29fZCf/bPgQEkfUGtCz8kro5swDYD+hAT1QW6WaBSGiOsC9aUf0cZRS\nLZ/OFD0SxkDRXmjfHbqfAF0H26n3tQX0TZ/aYA6Qt6vaQ/wVEg2QIHnkm1hKCdzEDHvxCKVUm6S1\nXI5ESS54yqBDD9/anDfD7qWwf3PN5xTugfi+4IyE/N1VXk5Ny6xYcAJsQN8fNNxSr8UjlFJtkgb0\nI1G0z2479LDb4b57w+lzaz6ncB90TLI/+Rmhr618kz4fXUuJu7xi1xPl13O/++eATtFXSoVHA3pd\nvnoM/hQXmmvuT0Fs7wvoHROh21DYOr/m6xTthQ7dIb43FGRV7J69YidffvwWJ3vX0plArvpO04PV\npj8CLJ46ToO5UqpOOoZelyXP2W354UC1w27Hw5Uvhi6ddtXLNrDXpHAfDLgALnsWIm02TGpaJjNS\nv2CuczEA/SWLg6YjADc6vyTNO4CCuPAXo1BKtW0a0OvS8yRwl4SWru3QA066PvS47ifUfp0pW+y4\ne3THil3T56VzkmcHOO3z4xx7WOYZQgyH+ZvrFZ7y3kD/8bqOplIqPG1vyGXVW7D63fCPz8+EDj1t\nUPfbuwYyVlQ9dumLsPzV6q/jirbBfO8aW8+lIIusvBIGOwLj6f3FDsV0lXwAxo4apkMtSqmwtb0e\n+qqZdvjkxGvrPtYYm5FycCt83R/O/5Pdv/BJG5jvXRl6/KZP7fh4yu2h+3PS2frJP5mSeSbRhbv4\nb+S/mZyWjGEIg2U3W7092WD6cNB0AGBoh8NQBicP1eEWpVT42lYPfcdi2LEIMpaFV1DLGLjR15s/\nsDWwv2if7bVX1n8c5GysksWybOki+m97k8KCfDJMVwDiyuyN1cGyi3TTm3vc9/GC53JiXE5uP8m3\n0n37bvX+iEqptqttBfSDlYJyXRwO6D8WBo6Hg9sD+wt9GSuV9T/PbrcuCNn93aq1AGSbePaazniN\nkCT7Afi1+2e8UG5nffrTE09LcNsT21fzHkopVYO2NeSSuzPwOCc9kEdek4Pb7NBKhx52XU9/r76m\nHnq34ymJ7sa3H7/FT97tTFyMCxG4syyHUqeLfNoBQg5xJGED+gpjJwud7VjNv52vQNJHEHcDJJ8F\n7bo2wodWSrUVbauHnrfTztSEGqffh9j8Obz7IxvQ3YdsIC8tAHdxSO85NS2TMY/NJ3naXFKLTsDt\nLsMAeSVucovddJdcsk08IADsMD1wSTlDZQcTHEuJoJyI9l3hUI6dbRrVHroNAYez8f8bKKVarbbX\nQ+99Gtzwli1kVZf8DPsFcPxlIE5bLTEiBm79yE7jJzBlv8TtAWBa+U/wB26/GMrY6yt9C3Bd2e8B\n4XcR/+EW5+eM8v6bSeefA3OB/Zvgh7fte/lnoCqlVBjaVkDvfSq0SwgvmAMUZNrJQj2GhyzdlprX\nn+nvpJOVtxaHCJ6QG6zi+18vxvcH0M/c9yN4K46Ij4lEBIa4d7PD0ZtHJp7EpSOTYFEiHNhiSwhE\nx2tAV0rVS9sK6OMfsdt1qbBpHlz5fO3H52dCR1+Z26IcKD/MF6t3svCLz9hfNhJDZKVgbj0R8SLd\nJZdb3VMr9vmDe1J8DIuvj4RFT8O2dXDi9Qzx55p3HWB76EXZtoKjUkrVQzgrFvUWkQUisl5E1onI\nfb79fxKRTBFZ5fu5uOmbGyZjYPX/YNvXgX1eb+CmZu52+OG/cDi/9usUZAbqlr92Ecz7LRsW/o+n\nHf8gkvKaTyOW0x0biKKMKMp4zvV3znSsCVRMLDsEW78E44XuQwMnDr3CZsoU7dOURaVUvYVzU7Qc\neMAYMxQ4HbhbRPxR6BljzAjfTy2lBo+yvWtg1k9gedAiSlu/hL8lwp4fAr3fnE21X+eWVDj3Qfu4\nc384uJ3owzmUmEgKianxtG+9w4gSN2dHb2NQTCEXO7/n+NiiQMXEuN6Bg7sFBfRTfgKjf25LBGjK\nolKqnuoM6MaYPcaYlb7Hhdj1RJv3fPTvX7TbrFWBfbk7AtkpCb6Avj+99ut0HQCdj7OPu/SHg9vo\nG1kQkrHi5xRBsEMqkyZeA+LkpbOK+ei2gQD87vqxgWn88b6AnnI79Bkd+p7+7BvtoSul6qleY+gi\nkoxdjm4pMAa4R0R+BCzH9uJzG7uBR2TfervN2wnFByG2s30cEW0Dutdjs1dyagnouTth4xwYdhV0\n6M4PhzpzkvsQ/b1byCY+5NAYl7NqvfLVo2D7QugxzD4PznmP6mBveooDImMD+wuy4MVz4IK/wpBL\nGvgfQSnV1oSdhy4i7YH3gV8aYwqA54H+wAhgD/BUDedNFpHlIrI8JyenEZpcB6/XTr/3D2Xs8fXS\nc3faoQ4RcEZA0ig7hl2TzBUwbxq3/HMuyVPn8NQKO2be37GHbBNf0T+vcfGJU35ivwwq1073i+pY\ndbZq+x72S6dwb2h1R6WUCkNYPXQRcWGD+UxjzCwAY8y+oNdfAj6u7lxjzAxgBkBKSkoYBVQaKG+H\nHVo56Qb4/Pe2t95/nO2hd+obOO72T6s93b+u54TCL3jIBT8UtANgjTeZKe7JbPP2JId4DHWs8XnS\ndXa75HnokGj/Sgh28/v2iyWYwwHlJbDk/2whsIhIlFIqXHUGdBER4BVggzHm6aD9PY0xe3xPrwTW\nNk0T68k/3NL3DHhgU2As+oQroV1gXNofuLPySiqm6OcWuxHAAIkRBygy0RRgh0Ry6cj/POeGvFVW\nXgm1OrQfks+E0++q+lrCoNrPdbpqf10ppSoJp4c+BrgFWCMi/ruMvwVuEJER2Pi3A7izSVpYX/3H\nwU/m2wUnXNGB/Wfeb4P4Y/PJzCthpGzhOdfr/Fp+xuaSXox2rGNSxCKmlN8JCD3lAHtMF4Jvfp7u\nWM9tznk8Un4ju013EuPrGBZ578dQnAt3fRN++2+ba6tBitR9rFJKBakzoBtjvqFySofVfNIUg0XG\nkprTg+n/+ZYu+eu4K+oTnpWb2FPiJJ/2GN9HOYyLkxzbmORcxKvlE+jJAa6JWMjbnrGsMIN9AT10\nmGSG62k6SjEzPeexPyLR5pTXpt/ZMP9h+OwhuPDh8NqfPMb+KKVUPbW64lxr//cws2a9Q2ZeCbFy\nmAks5hb3e6yKvpMTJFACd5vpidcId0V8xHTXi8zznkKJiWSS0/ambyh7iF+5fx5y7UPYHr+07179\njdDK+p1rt6vearTPp5RSNWkVAd1f7XDw1FSGrH2KU8wPAKz1JgNwsXMpABkmoeKcUiLZZnqSabrw\nUPntHCKGed4ULnEuIRI3JUSzn7iQbJbcobcA8O9fTgxvabjEkTZrZdzvGu2zKqVUTVp8LZfgaofH\nSxYR4iXdayfuFBHLVm9P+jv2UGhiyKN9yLk/dk/hkInhIHbh5lTPmUx0fsu1zq/oK/v4KvYirplw\nfiB4m7FQNiX84l7OCHgojIU0lFKqEbT4gD59XnpF6drBshuAdBOYWr/G9KM/e3y989BbARmmOwaI\n92W5fFM8nJWOYfx4cDn9t83lp9dNhv5BPXGR8IO5UkodZS0+oAenDg5x7KbMONluApN40rwDmej8\nlv3G9sL9aYlJ8TFMGT+4mqGTy2098m1vBgpzKaVUC9DiA3pifAyZvqDeR/ax1SRSTgROEbzGMDvy\nUgztyfTE1BLEKzmwxW47Nu+SNUopFazFB/Qp4wdXjKH/3H0fHSmuprbKhfW76MLpdhtcZ0UppZq5\nFh/Q/UHbP+uzQ3zX8Hrhtbn5/dAFpZVSqgVo8QEdbFCf2DEdfvgQzp0KnRs4VDLg/MZpmFJKHUWt\nIqBz6ACk/hwi24eWqVVKqTakZUwsMsauwekvRVv5tQ/vgeIDcPUrWnZWKdVmtYyAfnAbLHwS3rzS\nLlgRbNnLkD7XlpvtedKxaJ1SSjULLSOgd+kP18+06YQzr4HSQlLTMjnn0Xns+Hg63zlGkhp1+bFu\npVJKHVMtZwy9/1i4+jV490eUPjWcl0oeYqe7J1fyZ6QMSj5YB+JoWHaLUkq1YC2jh+53/KUwaQbr\nypOIKi8E7MITB+lIidvD9Hl1LPqslFKtWMvpofsNv5qrZsZQ3Vp2da4gpJRSrVidPXQR6S0iC0Rk\nvYisE5H7fPs7i8jnIrLZt+3U9M21alopqM4VhJRSqhULZ8ilHHjAGDMUOB24W0SGAlOBL40xA4Ev\nfc+PiinjBxPjcobsi3E5615BSCmlWrE6A7oxZo8xZqXvcSGwAUgCrgDe8B32BjCxqRpZ2cSRSTw6\naThJ8TEItnJiWCsIKaVUK1avMXQRSQZGAkuB7saYPb6X9gLdG7VldZg4MkkDuFJKBQk7y0VE2gPv\nA780xhQEv2aMMVDtfUpEZLKILBeR5Tk5OQ1qrFJKqZqFFdBFxIUN5jONMbN8u/eJSE/f6z2B7OrO\nNcbMMMakGGNSEhISqjtEKaVUIwgny0WAV4ANxping176ELjV9/hWYHbjN08ppVS4whlDHwPcAqwR\nkVW+fb8FHgPeFZE7gJ3AtU3TxIDUtMyKuueJ4a4+pJRSbUSdAd0Y8w2VV1cOOK9xm1OVP4hn5pVU\nrAcKkJlXwrRZawA0qCulFM186n9qWibTZq2pWDO08l1Xne6vlFIBzTqgT5+XTonbU+sxOt1fKaWs\nZh3QwwnWOt1fKaWsZh3Q6wrWOt1fKaUCmnVAr65mi//urE73V0qpUM26fK4/WGuqolJK1a1ZB3TQ\nmi1KKRWuZj3kopRSKnwa0JVSqpXQgK6UUq2EBnSllGolNKArpVQrIXZtiqP0ZiI52MqMR6IrsL8R\nm9NS6Odue9rqZ9fPXbO+xpg6F5Q4qgG9IURkuTEm5Vi342jTz932tNXPrp+74XTIRSmlWgkN6Eop\n1Uq0pIA+41g34BjRz932tNXPrp+7gVrMGLpSSqnataQeulJKqVq0iIAuIheJSLqIbBGRqce6PU1F\nRHqLyAIRWS8i60TkPt/+ziLyuYhs9m07Heu2NgURcYpImoh87HveT0SW+n7v74hI5LFuY2MTkXgR\neU9ENorIBhEZ3RZ+3yJyv+/f+FoReUtEolvj71tEXhWRbBFZG7Sv2t+vWM/6Pv9qETm5vu/X7AO6\niDiB/wMmAEOBG0Rk6LFtVZMpBx4wxgwFTgfu9n3WqcCXxpiBwJe+563RfcCGoOePA88YYwYAucAd\nx6RVTesfwKfGmCHASdjP36p/3yKSBNwLpBhjhgFO4Hpa5+/7deCiSvtq+v1OAAb6fiYDz9f3zZp9\nQAdOBbYYY7YZY8qAt4ErjnGbmoQxZo8xZqXvcSH2/9xJ2M/7hu+wN4CJx6aFTUdEegGXAC/7ngsw\nDnjPd0ir+9wiEgecDZwxs6gAAAJ5SURBVLwCYIwpM8bk0QZ+39jS3TEiEgHEAntohb9vY8xC4GCl\n3TX9fq8A/m2sJUC8iPSsz/u1hICeBOwOep7h29eqiUgyMBJYCnQ3xuzxvbQX6H6MmtWU/g78BvD6\nnncB8owx5b7nrfH33g/IAV7zDTW9LCLtaOW/b2NMJvAksAsbyPOBFbT+37dfTb/fBse6lhDQ2xwR\naQ+8D/zSGFMQ/JqxaUmtKjVJRC4Fso0xK451W46yCOBk4HljzEjgEJWGV1rp77sTtjfaD0gE2lF1\nWKJNaOzfb0sI6JlA76DnvXz7WiURcWGD+UxjzCzf7n3+P7182+xj1b4mMga4XER2YIfUxmHHluN9\nf5JD6/y9ZwAZxpilvufvYQN8a/99nw9sN8bkGGPcwCzsv4HW/vv2q+n32+BY1xIC+jJgoO8OeCT2\n5smHx7hNTcI3bvwKsMEY83TQSx8Ct/oe3wrMPtpta0rGmGnGmF7GmGTs73e+MeYmYAFwte+w1vi5\n9wK7/799O0aJGIjCAPxNtWCnR7CxtdzCQrDbQ9h4DCvPYmFhI2KpFxALURFR9xDWFrGYEbZZWIsl\nOv4fDIQ0mckfHuRlUkrZaacO8KzzvNVWy7SUstGe+e91d533gmX5XuKw7XaZ4mOhNbOaYRh+/cAM\nr5jjeOz5rHGde+rr1wPu25ip/eQbvOEaW2PPdY33YB9X7Xgbt3jHOSZjz28N693FXcv8Apv/IW+c\n4AVPOMWkx7xxpn4n+FTfyI6W5Yui7uib41HdBfSj6+VP0YiITvyFlktERKwgBT0iohMp6BERnUhB\nj4joRAp6REQnUtAjIjqRgh4R0YkU9IiITnwBdQCMUhE2xnYAAAAASUVORK5CYII=\n", "text/plain": [ "<Figure size 432x288 with 1 Axes>" ] }, "metadata": { "tags": [] }, "output_type": "display_data" } ], "source": [ "yprednorm = Xnorm.dot(theta)\n", "\n", "ypred = yprednorm*(maxy-miny) + miny\n", "plt.plot(y[sortidx,0],'o')\n", "plt.plot(ypred[sortidx,0],'--')" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "H7z0pVOS2mhw" }, "source": [ "The following function is to used to get an input, normalize it, perform predictions using the values of $\\theta$ encountered, and denormalizing the output to come back to the original range." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": {}, "colab_type": "code", "id": "90SPW1f9Rmkh" }, "outputs": [], "source": [ "def predict(x,theta):\n", " xnorm = (x-minx)/(maxx-minx)\n", " yprednorm = xnorm.dot(theta)\n", " ypred = yprednorm*(maxy - miny) + miny\n", " return ypred" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "QDxDwAa521v3" }, "source": [ "Let's use our model to predict the price of a house with 73 square meters, 1 bedroom and 1 bathroom." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 34 }, "colab_type": "code", "id": "iaRkEiDBdNat", "outputId": "b43e5487-7642-4a1f-d5da-4462fb3b7a94" }, "outputs": [ { "data": { "text/plain": [ "array([54.35879773])" ] }, "execution_count": 22, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "x = np.array([1,73,1,1])\n", "\n", "predict(x,theta)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "azSumKQw3AaA" }, "source": [ "To confirm that the model is statistically significant, use the Pearson correlation and evaluate the predicted output against the observed one." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 34 }, "colab_type": "code", "id": "hOtZRXKRda8C", "outputId": "72eb9fc9-f8db-491b-884f-eeb6b3620446" }, "outputs": [ { "data": { "text/plain": [ "(0.9753195824331272, 1.930066280643269e-65)" ] }, "execution_count": 23, "metadata": { "tags": [] }, "output_type": "execute_result" } ], "source": [ "pearsonr(ypred.reshape(-1),y.reshape(-1))" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "ajpwpjJh3MLZ" }, "source": [ "Notice here that the value of pearson correlation is 0.97 which indicates high correlation.\n", "The second value (1.93e-65) is the p-value. It must be a very low value for us to reject the null hypothesis, i,e there is no initial hypothesis of linear correlation. Assuming 1% significance level, we have:\n", "0.01 >> 1.93e-65 thus we can reject the null hypothesis thus indicating linear correlation between the prediction and observations.\n", "\n", "See you next post!" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": {}, "colab_type": "code", "id": "XDsWEeE74R1v" }, "outputs": [], "source": [] } ], "metadata": { "colab": { "collapsed_sections": [], "name": "linearRegressionNumpy.ipynb", "provenance": [], "version": "0.3.2" }, "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.2" } }, "nbformat": 4, "nbformat_minor": 1 }