548 lines
65 KiB
Plaintext
548 lines
65 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### CT4101 Machine learning, Semester 1 2024-2025\n",
|
|
"# Linear regression using numpy and scikit-learn\n",
|
|
"\n",
|
|
"#### 18 September 2024"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Machine learning (ML) algorithms allow a computer program to learn to improve its performance at a specified task with experience. ML algorithms may be applied to tasks such as classification, clustering, regression, transcription, translation, anomaly detection and sequential decision making (e.g. in Markov decision processes or Markov games). The benefits of ML research are beginning to be felt in society; outputs of ML research include useful everyday services that can recognise handwriting, protect users from email spam and financial fraud, and recommend suitable products or services according to a customer's preferences."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Linear regression in one or more variables is one of the most common tasks which must be solved by machine learning practitioners. Given a set of measurements for $j$ independent scalar variables $\\{x_{1},x_{2},...,x_{j}\\}$ and a corresponding set of observations for a dependent scalar variable $y$, the goal of linear regression is to fit a model which may be used to predict $y$ for any given values of $\\{x_{1},x_{2},...,x_{j}\\}$. Formally (Eqn. 1): \n",
|
|
"\n",
|
|
"\\begin{equation}\n",
|
|
" y = \\theta_{0} + \\theta_{1} x_1 + \\theta_{2} x_2 + ... + \\theta_{j} x_j\n",
|
|
"\\end{equation}\n",
|
|
"\n",
|
|
"where $\\{\\theta_{1},\\theta_{2},...,\\theta_{j}\\}$ are the weights for each independent variable and $\\theta_{0}$ accounts for any constant observed effect on the value of $y$ (similar to the intercept value in the equation for a simple straight line on a 2D plane). For a given training set, weights must be learned or found which minimise the error (defined by a cost function) between the predicted and actual values for the target variable $y$. One commonly used ML technique for solving linear regression problems is to apply gradient descent to minimise a squared error cost function."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"This notebook demonstrates how a simple linear regression problem may be solved (by simple linear regression, we mean that there is only one indepenent variable).\n",
|
|
"Simple linear regression creates a model that is equivalent to the representation of a straight line on a 2D plane, and you may remember the following well-known equation from your school days (Eqn. 2):\n",
|
|
"\n",
|
|
"\\begin{equation}\n",
|
|
" y = m x + c\n",
|
|
"\\end{equation}\n",
|
|
"\n",
|
|
"where $m$ is the slope of the line, and $c$ is the intercept value. The intercept accounts for constant factors that affect the value of $y$. Eqn. 1 above generalises Eqn. 2 to the case where there is more than one independent variable that influences the value of $y$.\n",
|
|
"\n",
|
|
"Where there is just one independent variable, Eqn. 1 and Eqn. 2 are equivalent: in Eqn. 1 $\\theta_{0}$ is the intercept ($c$ in Eqn. 2), while $\\theta_{1}$ in Eqn. 1 is equivalent to $m$ in Eqn. 2.\n",
|
|
"\n",
|
|
"Let's begin by setting up a simple training dataset in a matrix (numpy array), and then extracting the column vectors for the independent variable $x$ and the dependent variable $y$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Data in pandas dataframe: \n",
|
|
" independent variable dependent variable\n",
|
|
"0 30 70\n",
|
|
"1 40 90\n",
|
|
"2 40 100\n",
|
|
"3 50 120\n",
|
|
"4 50 130\n",
|
|
"5 50 150\n",
|
|
"6 60 160\n",
|
|
"7 70 190\n",
|
|
"8 70 200\n",
|
|
"9 80 200\n",
|
|
"10 80 220\n",
|
|
"11 80 230\n",
|
|
"Numpy 2D array:\n",
|
|
" [[ 30 70]\n",
|
|
" [ 40 90]\n",
|
|
" [ 40 100]\n",
|
|
" [ 50 120]\n",
|
|
" [ 50 130]\n",
|
|
" [ 50 150]\n",
|
|
" [ 60 160]\n",
|
|
" [ 70 190]\n",
|
|
" [ 70 200]\n",
|
|
" [ 80 200]\n",
|
|
" [ 80 220]\n",
|
|
" [ 80 230]]\n",
|
|
"x:\n",
|
|
" [[30]\n",
|
|
" [40]\n",
|
|
" [40]\n",
|
|
" [50]\n",
|
|
" [50]\n",
|
|
" [50]\n",
|
|
" [60]\n",
|
|
" [70]\n",
|
|
" [70]\n",
|
|
" [80]\n",
|
|
" [80]\n",
|
|
" [80]]\n",
|
|
"y:\n",
|
|
" [[ 70]\n",
|
|
" [ 90]\n",
|
|
" [100]\n",
|
|
" [120]\n",
|
|
" [130]\n",
|
|
" [150]\n",
|
|
" [160]\n",
|
|
" [190]\n",
|
|
" [200]\n",
|
|
" [200]\n",
|
|
" [220]\n",
|
|
" [230]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import pandas as pd\n",
|
|
"\n",
|
|
"# dataset 1. Here we write the data directly into a 2D numpy array\n",
|
|
"\"\"\"data = np.array([\n",
|
|
" [80, 28.3],\n",
|
|
" [110, 51.5],\n",
|
|
" [110, 47.3],\n",
|
|
" [130, 67.4]\n",
|
|
"])\n",
|
|
"\n",
|
|
"\"\"\"\n",
|
|
"# dataset 2. Here we load a different dataset in from the file external_data.csv using the pandas library\n",
|
|
"df = pd.read_csv(\"external_data.csv\")\n",
|
|
"print(\"Data in pandas dataframe: \\n\", df)\n",
|
|
"data = df.to_numpy()\n",
|
|
"\n",
|
|
"print(\"Numpy 2D array:\\n\", data)\n",
|
|
"\n",
|
|
"# Create and print x - the column vector of measured values for the independent variable\n",
|
|
"x = data[:,0].reshape((-1, 1)) # the -1 here means the length of the column vector will be inferred\n",
|
|
"print(\"x:\\n\", x)\n",
|
|
"\n",
|
|
"# Create and print y - the column vector of observed values for the dependent variable\n",
|
|
"y = np.array([data[:,1]]).reshape(-1,1)\n",
|
|
"print(\"y:\\n\", y)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Now we'll plot the values of $x$ and $y$ using the scatter plot from the matplotlib library (https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.scatter.html)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAG0CAYAAADJpthQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGZElEQVR4nO3de3hU1b3/8c+EJJMbSQiBXJBw0wIhEiUoiQi2hQNYpV44PYGChrvlB4hgrHhDRYQePQ0tfSoYBLGKQC3qOWKLKAKSnhCQABrUQBANmHCRkITcb/v3B4cxMwmWzEwymeH9ep55ns7aO3u+K9tmPuy91tomwzAMAQAAwMLL1QUAAAC0NwQkAAAAGwQkAAAAGwQkAAAAGwQkAAAAGwQkAAAAGwQkAAAAG96uLqA9amhoUEFBgTp27CiTyeTqcgAAwBUwDEMXLlxQdHS0vLwcuwZEQGpGQUGBunfv7uoyAACAHU6cOKFrrrnGoWMQkJrRsWNHSRd/wcHBwS6uBgAAXInS0lJ1797d8j3uCAJSMy7dVgsODiYgAQDgZpwxPIZB2gAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAeLjSqloVllQ2u62wpFKlVbVtXFH7R0ACAMCDlVbVKmXtXiW/vEcFxdYhqaC4Uskv71HK2r2EJBsEJAAAPFh5dZ3OldUov6hC49N/CEkFxZUan75H+UUVOldWo/LqOhdX2r4QkAAA8GBRIf7aODNRMWEBlpC0/9siSziKCQvQxpmJigrxd3Wp7YrJMAzD1UW0N6WlpQoJCVFJSQkPqwUAeITGV4wuuRSOokM9Ixw58/ubK0gAAFwFokP9tTw53qpteXK8x4QjZyMgAQBwFSgortT8TYes2uZvOtRk4DYuIiABAODhGt9eiwkL0OZZSVZjkghJTRGQAADwYIUllU0GZCf0CGsycPty6yRdrQhIAAB4sECztzoH+TYZkB0d+sPsts5Bvgo0e7u40vaFWWzNYBYbAMCTlFbVqry6rtmp/IUllQo0eyvYz8cFlTmXM7+/iYsAAHi4YD+fywYg1j9qHrfYAAAAbBCQAAAAbBCQAAAAbBCQAADwcKVVtZedxl9YUqnSqto2rqj9IyABAODBSqtqlbJ2r5JfbrogZEFxpZJf3qOUtXsJSTYISAAAeLDy6jqdK6tpsmp249W1z5XVqLy6zsWVti8EJAAAPFhUiH+TVbP3f1vUZHVtpvtbY6HIZrBQJADA0zS+YnSJ7era7s6Z399cQQIA4CoQHeqv5cnxVm3Lk+M9Jhw5GwEJAICrQEFxpeZvOmTVNn/ToSYDt3ERAQkAAA/X+PZaTFiANs9KshqTREhqioAEAIAHKyypbDIgO6FHWJOB25dbJ+lqRUACAMCDBZq91TnIt8mA7OjQH2a3dQ7yVaCZ59c3xiy2ZjCLDQDgSUqralVeXdfsVP7CkkoFmr0V7Ofjgsqcy5nf38RFAAA8XLCfz2UDEOsfNY9bbAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISAAAADYISACAq15pVa0KSyqb3VZYUqnSqto2rsi5PL1/raFdBaRly5bppptuUseOHdW1a1fdfffdys3NtdqnqqpKs2fPVufOnRUUFKRx48bp9OnTVvvk5+frjjvuUEBAgLp27apHHnlEdXV1bdkVAICbKK2qVcravUp+eY8Kiq1DREFxpZJf3qOUtXvdNkR4ev9aS7sKSLt27dLs2bO1Z88effjhh6qtrdWoUaNUXl5u2Wf+/Pl677339NZbb2nXrl0qKCjQvffea9leX1+vO+64QzU1Nfrf//1fvfbaa1q3bp0WLVrkii4BANq58uo6nSurUX5Rhcan/xAiCoorNT59j/KLKnSurEbl1e75D21P719rMRmGYbi6iMs5e/asunbtql27dmn48OEqKSlRly5d9Oabb+rf//3fJUlfffWV+vfvr8zMTCUmJuof//iH7rzzThUUFCgiIkKStGrVKj366KM6e/asfH19/+XnlpaWKiQkRCUlJQoODm7VPgIAXK9xWIgJC9Dy5HjN33TI8n7jzERFh/q7uky7eXr/LnHm93e7uoJkq6SkRJIUFhYmSdq/f79qa2s1cuRIyz79+vVTTEyMMjMzJUmZmZm6/vrrLeFIkkaPHq3S0lIdPny42c+prq5WaWmp1QsAcPWIDvXXxpmJigkLUH5RhcatzPSo8ODp/WsN7TYgNTQ06KGHHtLQoUMVFxcnSTp16pR8fX0VGhpqtW9ERIROnTpl2adxOLq0/dK25ixbtkwhISGWV/fu3Z3cGwBAexcd6q/lyfFWbcuT4z0mPHh6/5yt3Qak2bNnKycnRxs3bmz1z3rsscdUUlJieZ04caLVPxMA0L4UFFdq/qZDVm3zNx1qMrDZXXl6/5ytXQakOXPmaMuWLdqxY4euueYaS3tkZKRqampUXFxstf/p06cVGRlp2cd2Vtul95f2sWU2mxUcHGz1AgBcPWzH6GyelWS5HdV4YLO78vT+tYZ2FZAMw9CcOXP0zjvv6OOPP1avXr2stickJMjHx0fbt2+3tOXm5io/P19JSUmSpKSkJH3++ec6c+aMZZ8PP/xQwcHBio2NbZuOAADcRmGJdXjYODNRCT3CrMbsjE/fc9l1hNo7T+9fa2lXAWn27Nl644039Oabb6pjx446deqUTp06pcrKiyctJCRE06ZN04IFC7Rjxw7t379fU6ZMUVJSkhITEyVJo0aNUmxsrO677z4dOnRIH3zwgZ588knNnj1bZrPZld0DALRDgWZvdQ7ybTJgufHA5s5Bvgo0e7u4Uvt4ev9aS7ua5m8ymZptf/XVVzV58mRJFxeKfPjhh7VhwwZVV1dr9OjReumll6xun3377beaNWuWdu7cqcDAQKWkpOh3v/udvL2v7OQzzR8Ari6lVbUqr65TVEjTAcuFJZUKNHsr2M/HBZU5h6f37xJnfn+3q4DUXhCQAABwP1fNOkgAAACuQEACAACwQUACAACwQUACAFz1SqtqLzvNvbCkkifdX4UcCki1tbU6ceKEcnNzVVRU5KyaAABoM6VVtUpZu1fJLzddMLGguFLJL+9Rytq9hKSrTIsD0oULF7Ry5UrddtttCg4OVs+ePdW/f3916dJFPXr00IwZM7Rv377WqBUAAKcrr67TubKaJqtKN159+lxZjcqr61xcKdpSiwJSWlqaevbsqVdffVUjR47Uu+++q4MHD+rIkSPKzMzU008/rbq6Oo0aNUpjxozR0aNHW6tuAACcIirEv8mq0vu/LWqy+nRzawjBc7VoHaQJEyboySef1IABA350v+rqar366qvy9fXV1KlTHS6yrbEOEgBcfRpfMbrEdvVptG8sFNnKCEgAcHXa/22Rxq3MtLzfPCtJCT3CXFgRWqLdLBS5e/duTZo0SUlJSfruu+8kSa+//royMjIcKgoAgLZWUFyp+ZsOWbXN33SIJ91fpewOSJs3b9bo0aPl7++vAwcOqLq6WpJUUlKipUuXOq1AAABaW+PbazFhAdo8K8lqTBIh6epjd0BasmSJVq1apdWrV8vH54cH3A0dOlTZ2dlOKQ4AgNZWWFLZZEB2Qo+wJgO3L7dOEjyT3QEpNzdXw4cPb9IeEhKi4uJiR2oCAKDNBJq91TnIt8mA7OjQH2a3dQ7yVaDZ28WVoi3ZfbYjIyOVl5ennj17WrVnZGSod+/ejtYFAECbCPbz0WtTb1Z5dV2TqfzRof7a9ECiAs3eCvbzucwR4InsvoI0Y8YMzZs3T1lZWTKZTCooKND69euVmpqqWbNmObNGAABaVbCfz2XXOYoK8SccXYXsvoK0cOFCNTQ0aMSIEaqoqNDw4cNlNpuVmpqquXPnOrNGAACANuXwOkg1NTXKy8tTWVmZYmNjFRQU5KzaXIZ1kAAAcD/O/P52eMSZr6+vYmNjHT0MAABAu9GigLRgwYIr3jctLa3FxQAAALQHLQpIBw4cuKL9TCaTXcUAAAC0By0KSDt27GitOgAAANoNh57FdolhGOKZtwAAwFM4FJDWrFmjuLg4+fn5yc/PT3FxcXrllVecVRsAAIBL2D2LbdGiRUpLS9PcuXOVlJQkScrMzNT8+fOVn5+vxYsXO61IAACAtmT3OkhdunTRihUrNGHCBKv2DRs2aO7cufr++++dUqArsA4SAADux5nf33bfYqutrdXgwYObtCckJKiurs6hogAAAFzJ7oB03333aeXKlU3a09PTNXHiRIeKAgAAcCW7F4o0mUx65ZVXtG3bNiUmJkqSsrKylJ+fr/vvv9+5VQIAALQhhxaKTEhIkCQdO3ZMkhQeHq7w8HAdPnzYSeUBAAC0PRaKBAAAsOGUhSIBAAA8id3rIF3yxRdfKD8/XzU1NVbtv/zlLx09NAAAgEvYHZC+/vpr3XPPPfr8889lMpksjxq59KDa+vp651QIAADQxuy+xTZv3jz16tVLZ86cUUBAgA4fPqxPPvlEgwcP1s6dO51YIgAAQNuy+wpSZmamPv74Y4WHh8vLy0teXl669dZbtWzZMj344INNZrwBAAC4C7uvINXX16tjx46SLk7vLygokCT16NFDubm5zqkOAADABey+ghQXF6dDhw6pV69eGjJkiF544QX5+voqPT1dvXv3dmaNAAAAbcrugPTkk0+qvLxckrR48WLdeeedGjZsmDp37qxNmzY5rUAAAIC2ZjIuTT9zgqKiInXq1Mkyk81dOfNpwAAAoG048/vb4XWQGgsLC3Pm4QAAAFyixQ+rfe655xQYGGj14NrmpKWlOVQYAACAq7T4YbW1tbWW/3057n6LDQAAXN3selhtbW2tvLy8tGrVKl133XWtUhgAAICr2LUOko+Pjz777DNn1wIAANAu2L1Q5KRJk7RmzRpn1gIAANAu2D2Lra6uTmvXrtVHH32khIQEBQYGWm1nkDYAAHBXdgeknJwcDRo0SJJ05MgRq20M0gYAAO7M7oB0acA2AACAp7F7DBIAAICncngl7S+++EL5+fmqqamxav/lL3/p6KEBAO1IaVWtyqvrFBXi32RbYUmlAs3eCvbzcUFlgPPZHZC+/vpr3XPPPfr8889lMpl06ZFul8Yf1dfXO6dCAIDLlVbVKmXtXp0rq9HGmYmKDv0hJBUUV2p8+h51DvLVa1NvJiTBI9h9i23evHnq1auXzpw5o4CAAB0+fFiffPKJBg8erJ07dzqxRACAq5VX1+lcWY3yiyo0Pn2PCoorJf0QjvKLKnSurEbl1XUurhRwDrsDUmZmphYvXqzw8HB5eXnJy8tLt956q5YtW6YHH3zQmTUCAFwsKsRfG2cmKiYswBKS9n9bZAlHMWEB2jgzsdnbb4A7sjsg1dfXq2PHjpKk8PBwFRQUSJJ69Oih3Nxc51QHAGg3okOtQ9K4lZlW4ajxbTfA3dkdkOLi4nTo0CFJ0pAhQ/TCCy/on//8pxYvXqzevXs7rUAAQPsRHeqv5cnxVm3Lk+MJR/A4dgekJ598Ug0NDZKkxYsX6/jx4xo2bJj+/ve/a8WKFU4rEADQfhQUV2r+pkNWbfM3HbKMSQI8hcm4NP3MCYqKitSpUye3X0m7tLRUISEhKikpUXBwsKvLAYB2ofGA7JiwAC1Pjtf8TYe4zYZ2w5nf33ZfQZo+fXqT2WphYWFuH44AAE0VllQ2GZCd0COsycDtwhKuJMEz2B2Qzp49qzFjxqh79+565JFHdPDgQSeWBQBoTwLN3uoc5NvkSlHjgdudg3wVaHZ4/WGgXXDoFtv58+f11ltv6c0339Tu3bvVr18/TZw4Ub/+9a/Vs2dPJ5bZtrjFBgBNsZI22jtnfn87bQzSyZMntWHDBq1du1ZHjx5VXZ37LhZGQAIAwP20izFIjdXW1urTTz9VVlaWvvnmG0VERDjjsAAAAC7hUEDasWOHZsyYoYiICE2ePFnBwcHasmWLTp486az6AAAA2pzdo+m6deumoqIijRkzRunp6Ro7dqzMZrMzawMAAHAJuwPSM888o1/96lcKDQ11YjkAAACuZ3dAmjFjhjPrAACPwEwvwDM4ZZC2s3zyyScaO3asoqOjZTKZ9O6771ptnzx5skwmk9VrzJgxVvsUFRVp4sSJCg4OVmhoqKZNm6aysrI27AWAq1VpVa1S1u5V8st7mjx6o6C4Uskv71HK2r0qrap1UYUArlS7Ckjl5eWKj4/Xn//858vuM2bMGBUWFlpeGzZssNo+ceJEHT58WB9++KG2bNmiTz75RDNnzmzt0gFA5dV1OldWY1lV+lJIavyIjnNlNSqvdt9lUICrhVOfxeZMJpNJ77zzju6++25L2+TJk1VcXNzkytIlX375pWJjY7Vv3z4NHjxYkrR161b94he/0MmTJxUdHX1Fn806SADsxfPKANdpF+sg5efnq7lsZRiG8vPzHSrqx+zcuVNdu3ZV3759NWvWLJ07d86yLTMzU6GhoZZwJEkjR46Ul5eXsrKyLnvM6upqlZaWWr0AwB6NH72RX1ShcSszCUeAG7I7IPXq1Utnz55t0l5UVKRevXo5VNTljBkzRn/5y1+0fft2/ed//qd27dql22+/XfX19ZKkU6dOqWvXrlY/4+3trbCwMJ06deqyx122bJlCQkIsr+7du7dK/QCuDtGh/lqeHG/Vtjw5nnAEuBG7A5JhGDKZTE3ay8rK5Ofn51BRlzN+/Hj98pe/1PXXX6+7775bW7Zs0b59+7Rz506HjvvYY4+ppKTE8jpx4oRzCgZwVSoortT8TYes2uZvOtRk4DaA9qvF0/wXLFgg6eIYoaeeekoBAQGWbfX19crKytINN9zgtAJ/TO/evRUeHq68vDyNGDFCkZGROnPmjNU+dXV1KioqUmRk5GWPYzabWeQSgFP82Bik8el7uM0GuIkWB6QDBw5IungF6fPPP5evr69lm6+vr+Lj45Wamuq8Cn/EyZMnde7cOUVFRUmSkpKSVFxcrP379yshIUGS9PHHH6uhoUFDhgxpk5oAXL0KS6zD0aUwtHFmoqV9fPoebXogsdl1kgC0Hy0OSDt27JAkTZkyRX/84x+dOsurrKxMeXl5lvfHjx/XwYMHFRYWprCwMD377LMaN26cIiMjdezYMf32t7/Vtddeq9GjR0uS+vfvrzFjxmjGjBlatWqVamtrNWfOHI0fP/6KZ7ABgL0Czd7qHHTxH42NrxQ1Dkmdg3wVaLZ7jV4AbaRdTfPfuXOnfvaznzVpT0lJ0cqVK3X33XfrwIEDKi4uVnR0tEaNGqXnnntOERERln2Lioo0Z84cvffee/Ly8tK4ceO0YsUKBQUFXXEdTPMHYC9W0gZcx5nf3w4FpO3bt2v79u06c+aMGhoarLatXbvWocJciYAEAID7ceb3t93XeZ999lktXrxYgwcPVlRUVLMz2gAAANyR3QFp1apVWrdune677z5n1gMAAOBydq+DVFNTo1tuucWZtQAAALQLdgek6dOn680333RmLQAAAO2C3bfYqqqqlJ6ero8++kgDBw6Uj4/1rIy0tDSHiwMAAHAFuwPSZ599ZlkxOycnx2obA7YBAIA7szsgXVowEgAAwNPYPQZJknbv3q1Jkybplltu0XfffSdJev3115WRkeGU4gAAAFzB7oC0efNmjR49Wv7+/srOzlZ1dbUkqaSkREuXLnVagQAAAG3N7oC0ZMkSrVq1SqtXr7YaoD106FBlZ2c7pTgAAABXsDsg5ebmavjw4U3aQ0JCVFxc7EhNAAAALmV3QIqMjFReXl6T9oyMDPXu3duhogAAAFzJ7oA0Y8YMzZs3T1lZWTKZTCooKND69euVmpqqWbNmObNGAACANmX3NP+FCxeqoaFBI0aMUEVFhYYPHy6z2azU1FTNnTvXmTUCAAC0KZNhGIYjB6ipqVFeXp7KysoUGxuroKAgZ9XmMqWlpQoJCVFJSYmCg4NdXQ4AALgCzvz+tvsK0iW+vr6KjY119DAAAADtRosC0oIFC654X57FBgAA3FWLAtKBAwes3mdnZ6uurk59+/aVJB05ckQdOnRQQkKC8yoEAABoYy0KSI2fv5aWlqaOHTvqtddeU6dOnSRJ58+f15QpUzRs2DDnVgkAANCG7B6k3a1bN23btk0DBgywas/JydGoUaNUUFDglAJdgUHaAOxVWlWr8uo6RYX4N9lWWFKpQLO3gv18mvlJAI5y5ve33esglZaW6uzZs03az549qwsXLjhUFAC4o9KqWqWs3avkl/eooLjSaltBcaWSX96jlLV7VVpV66IKAVwpuwPSPffcoylTpujtt9/WyZMndfLkSW3evFnTpk3Tvffe68waAcAtlFfX6VxZjfKLKjQ+/YeQVFBcqfHpe5RfVKFzZTUqr65zcaUA/hW7b7FVVFQoNTVVa9euVW3txX8NeXt7a9q0aXrxxRcVGBjo1ELbErfYANircRiKCQvQ8uR4zd90yPJ+48xERYc2vf0GwHHO/P52eKHI8vJyHTt2TJLUp08ftw5GlxCQADiicUi6hHAEtL52MQbpksDAQA0cOFADBw70iHAEAI6KDvXX8uR4q7blyfGEI8CNOLSS9vbt27V9+3adOXNGDQ0NVtvWrl3rUGEA4K4Kiis1f9Mhq7b5mw5xBQlwI3ZfQXr22Wc1atQobd++Xd9//73Onz9v9QKAq5HtGKTNs5IUExbQZOA2gPbN7jFIUVFReuGFF3Tfffc5uyaXYwwSAHsUllycym87INs2NG16ILHZdZIAOKZdjEGqqanRLbfc4tCHA4AnCTR7q3OQb5MB2dGh/to4M1ExYQHqHOSrQLPDzwkH0MrsvoL06KOPKigoSE899ZSza3I5riABsBcraQOu48zvb7v/GVNVVaX09HR99NFHGjhwoHx8rP8Pn5aW5lBhAOCOgv18LhuAuK0GuA+7A9Jnn32mG264QdLF5681ZjKZHCoKAADAlewOSDt27HBmHQAAAO2GQwtF7t69W5MmTdItt9yi7777TpL0+uuvKyMjwynFAQAAuILdAWnz5s0aPXq0/P39lZ2drerqaklSSUmJli5d6rQCAcCdlFbVqrCk+bWOCksqVVpV28YVAbCH3QFpyZIlWrVqlVavXm01QHvo0KHKzs52SnEA4E5Kq2qVsnavkl9uuiBkQfHFNZJS1u4lJAFuwO6AlJubq+HDhzdpDwkJUXFxsSM1AYBbKq+u07mymiarZjdeKPJcWY3Kq+tcXCmAf8XugBQZGam8vLwm7RkZGerdu7dDRQGAO4oK+WFByEshaf+3RVaraG+cySragDuwOyDNmDFD8+bNU1ZWlkwmkwoKCrR+/XqlpqZq1qxZzqwRANxG41Wz84sqNG5lZpNHjwBo/+ye5r9w4UI1NDRoxIgRqqio0PDhw2U2m5Wamqq5c+c6s0YAcCvRof5anhyvcSszLW3Lk+MJR4AbsftRI5fU1NQoLy9PZWVlio2NVVBQkLNqcxkeNQLAEY3HHF3CFSSg9bWLh9Ve4uvrq/79++umm27yiHAEAI5oHI5iwgK0eVaS1Zgk29ltANonhwLSmjVrFBcXJz8/P/n5+SkuLk6vvPKKs2oDALdSWFLZZEB2Qo+wJgO3L7dOEoD2w+4xSIsWLVJaWprmzp2rpKQkSVJmZqbmz5+v/Px8LV682GlFAoA7CDR7q3OQryRZ3U67NHB7fPoedQ7yVaDZ7j+9ANqI3WOQunTpohUrVmjChAlW7Rs2bNDcuXP1/fffO6VAV2AMEgB7lVbVqry6rtmp/IUllQo0eyvYz6eZnwTgKGd+f9v9z5ja2loNHjy4SXtCQoLq6lgEDcDVKdjP57IBiPWPAPdh9xik++67TytXrmzSnp6erokTJzpUFAAAgCs5dCN8zZo12rZtmxITEyVJWVlZys/P1/33368FCxZY9ktLS3OsSgAAgDZkd0DKycnRoEGDJEnHjh2TJIWHhys8PFw5OTmW/Uwmk4MlAgAAtC27A9KOHTucWQcAAEC74fBCkQAAAJ7GoYC0e/duTZo0SUlJSfruu+8kSa+//royMjKcUhwAAIAr2B2QNm/erNGjR8vf318HDhxQdXW1JKmkpERLly51WoEAAABtze6AtGTJEq1atUqrV6+Wj88Pa34MHTpU2dnZTikOAADAFewOSLm5uRo+fHiT9pCQEBUXFztSEwAAgEvZHZAiIyOVl5fXpD0jI0O9e/d2qCgAAABXsjsgzZgxQ/PmzVNWVpZMJpMKCgq0fv16paamatasWc6sEQAAoE3ZvQ7SwoUL1dDQoBEjRqiiokLDhw+X2WxWamqq5s6d68waAQAA2pTJMAzDkQPU1NQoLy9PZWVlio2NVVBQkLNqcxlnPg0YAAC0DWd+f7foClLj56v9Kzx/DQAAuKsWBaQDBw5Yvc/OzlZdXZ369u0rSTpy5Ig6dOighIQE51UIAADQxloUkBo/fy0tLU0dO3bUa6+9pk6dOkmSzp8/rylTpmjYsGHOrRIAAKAN2T0GqVu3btq2bZsGDBhg1Z6Tk6NRo0apoKDAKQW6AmOQAABwP878/rZ7mn9paanOnj3bpP3s2bO6cOGCQ0UBAAC4kt0B6Z577tGUKVP09ttv6+TJkzp58qQ2b96sadOm6d5773VmjQAAAG3K7nWQVq1apdTUVP36179WbW3txYN5e2vatGl68cUXnVYgAM9TWlWr8uo6RYX4N9lWWFKpQLO3gv18mvlJAGgbdl9BCggI0EsvvaRz587pwIEDOnDggIqKivTSSy8pMDDQrmN+8sknGjt2rKKjo2UymfTuu+9abTcMQ4sWLVJUVJT8/f01cuRIHT161GqfoqIiTZw4UcHBwQoNDdW0adNUVlZmbzcBOFlpVa1S1u5V8st7VFBcabWtoLhSyS/vUcravSqtqnVRhQDgQEC6JDAwUAMHDtTAgQPtDkaXlJeXKz4+Xn/+85+b3f7CCy9oxYoVWrVqlbKyshQYGKjRo0erqqrKss/EiRN1+PBhffjhh9qyZYs++eQTzZw506G6ADhPeXWdzpXVKL+oQuPTfwhJBcWVGp++R/lFFTpXVqPy6joXVwrgaubwStqtxWQy6Z133tHdd98t6eLVo+joaD388MNKTU2VJJWUlCgiIkLr1q3T+PHj9eWXXyo2Nlb79u3T4MGDJUlbt27VL37xC508eVLR0dFX9NnMYgNaV+MwFBMWoOXJ8Zq/6ZDl/caZiYoObXr7DQB+TLuYxdbWjh8/rlOnTmnkyJGWtpCQEA0ZMkSZmZmSpMzMTIWGhlrCkSSNHDlSXl5eysrKuuyxq6urVVpaavUC0HqiQ/21cWaiYsIClF9UoXErMwlHANoVtwlIp06dkiRFRERYtUdERFi2nTp1Sl27drXa7u3trbCwMMs+zVm2bJlCQkIsr+7duzu5egC2okP9tTw53qpteXI84QhAu+A2Aak1PfbYYyopKbG8Tpw44eqSAI9XUFyp+ZsOWbXN33SoycBtAHAFtwlIkZGRkqTTp09btZ8+fdqyLTIyUmfOnLHaXldXp6KiIss+zTGbzQoODrZ6AWg9tmOQNs9KstxuazxwGwBcpUXrIC1YsOCK901LS2txMT+mV69eioyM1Pbt23XDDTdIujgYKysrS7NmzZIkJSUlqbi4WPv377c8MPfjjz9WQ0ODhgwZ4tR6ANinsMQ6HF0ac7RxZqKlfXz6Hm16ILHZdZIAoC20KCAdOHDA6n12drbq6urUt29fSdKRI0fUoUMHSzhpqbKyMuXl5VneHz9+XAcPHlRYWJhiYmL00EMPacmSJbruuuvUq1cvPfXUU4qOjrbMdOvfv7/GjBmjGTNmaNWqVaqtrdWcOXM0fvz4K57BBqB1BZq91TnIV5KsBmQ3Dkmdg3wVaLZ7HVsAcJjd0/zT0tK0c+dOvfbaa+rUqZMk6fz585oyZYqGDRumhx9+uMXH3Llzp372s581aU9JSdG6detkGIaefvpppaenq7i4WLfeeqteeukl/eQnP7HsW1RUpDlz5ui9996Tl5eXxo0bpxUrVigoKOiK62CaP9C6WEkbQGtw5ve33QGpW7du2rZtmwYMGGDVnpOTo1GjRqmgoMChwlyJgAQAgPtpF+sglZaW6uzZs03az549qwsXLjhUFAAAgCvZHZDuueceTZkyRW+//bZOnjypkydPavPmzZo2bZruvfdeZ9YIAADQpuweBblq1Sqlpqbq17/+tWprLz5U0tvbW9OmTdOLL77otAIBAADamsPPYisvL9exY8ckSX369HH4gbXtAWOQgNbFIG0AraFdjEGSpN27d+uBBx7Qb37zG3Xu3FmBgYF6/fXXlZGR4VBRADxXaVWtUtbuVfLLTReELCiuVPLLe5Sydq9Kq2pdVCEAOBCQNm/erNGjR8vf31/Z2dmqrq6WJJWUlGjp0qVOKxCAZymvrtO5spomq2Y3Xl37XFmNyqvrXFwpgKuZ3QFpyZIlWrVqlVavXi0fnx8uhQ8dOlTZ2dlOKQ6A54kKubggZONHi+z/tqjJ6tqsog3AlewOSLm5uRo+fHiT9pCQEBUXFztSEwAPd2nV7EshadzKzCaPHgEAV7I7IEVGRlo9FuSSjIwM9e7d26GiAHi+6FB/LU+Ot2pbnhxPOALQLtgdkGbMmKF58+YpKytLJpNJBQUFWr9+vVJTUy0PjwWAyykortT8TYes2uZvOtRk4DYAuILd6yAtXLhQDQ0NGjFihCoqKjR8+HCZzWalpqZq7ty5zqwRgIdpPCA7JixAy5PjNX/TIcuYJG6zAXA1h9dBqqmpUV5ensrKyhQbG9uih8K2V6yDBLSewpKLU/ltxxzZhqZNDzBQG0DLtIt1kKZPn66dO3fK19dXsbGxuvnmmz0iHAFoXYFmb3UO8m0yILvxwO3OQb4KNNt9gRsAHGb3X6CzZ89qzJgx6tKli8aPH69JkyYpPj7+X/8ggKtasJ+PXpt6c7MraUeH+mvTA4mspA3A5ey+gvTf//3fKiws1FNPPaV9+/Zp0KBBGjBggJYuXapvvvnGiSUC8DTBfj6XvX0WFeJPOALgcg6PQbrk5MmT2rBhg9auXaujR4+qrs59V8FlDBIAAO6nXYxBaqy2tlaffvqpsrKy9M033ygiIsIZhwUAAHAJhwLSjh07NGPGDEVERGjy5MkKDg7Wli1bdPLkSWfVBwAA0ObsHqTdrVs3FRUVacyYMUpPT9fYsWNlNpudWRsAAIBL2B2QnnnmGf3qV79SaGioE8sBAABwPbtusdXW1mrjxo06e/ass+sBAABwObsCko+Pjz777DNn1wIAANAu2D1Ie9KkSVqzZo0zawEAAGgX7B6DVFdXp7Vr1+qjjz5SQkKCAgMDrbanpaU5XBwAAIAr2B2QcnJyNGjQIEnSkSNHrLaZTCbHqgIAAHAhuwPSjh07nFkHAABAu+GUlbQBAAA8iUMBaffu3Zo0aZKSkpL03XffSZJef/11ZWRkOKU4AAAAV7A7IG3evFmjR4+Wv7+/Dhw4oOrqaklSSUmJli5d6rQCAQAA2prdAWnJkiVatWqVVq9eLR8fH0v70KFDlZ2d7ZTiAAAAXMHugJSbm6vhw4c3aQ8JCVFxcbEjNQEAALiU3QEpMjJSeXl5TdozMjLUu3dvh4oCAABwJbsD0owZMzRv3jxlZWXJZDKpoKBA69evV2pqqmbNmuXMGgEAANqU3esgLVy4UA0NDRoxYoQqKio0fPhwmc1mpaamau7cuc6sEQAAoE2ZDMMwHDlATU2N8vLyVFZWptjYWAUFBTmrNpcpLS1VSEiISkpKFBwc7OpycJUqrapVeXWdokL8m2wrLKlUoNlbwX4+zfwkAFydnPn9bfcttsrKSlVUVMjX11exsbGKiIjQK6+8om3btjlUEICL4Shl7V4lv7xHBcWVVtsKiiuV/PIepazdq9KqWhdVCACeze6AdNddd+kvf/mLJKm4uFhDhgzR73//e911111auXKl0woErkbl1XU6V1aj/KIKjU//ISQVFFdqfPoe5RdV6FxZjcqr61xcKQB4JrsDUnZ2toYNGyZJ+tvf/qaIiAh9++23+stf/qIVK1Y4rUDgahQV4q+NMxMVExZgCUn7vy2yhKOYsABtnJnY7O03AIDj7A5IFRUV6tixoyRp27Ztuvfee+Xl5aXExER9++23TisQuFpFh1qHpHErM63CUXQo4QgAWovdAenaa6/Vu+++qxMnTuiDDz7QqFGjJElnzpxhYDPgJNGh/lqeHG/Vtjw5nnAEAK3M7oC0aNEipaamqmfPnrr55puVlJQk6eLVpBtvvNFpBQJXs4LiSs3fdMiqbf6mQ00GbgMAnMuhaf6nTp1SYWGh4uPj5eV1MWvt3btXwcHB6tevn9OKbGtM80d70HhAdkxYgJYnx2v+pkPcZgOAy3Dm97fD6yBJ0qVDmEwmRw/VLhCQ4GqFJRen8tuGIdvQtOkBBmoDwCXtYh0kSVqzZo3i4uLk5+cnPz8/xcXF6ZVXXnGoIABSoNlbnYN8m1wpajxwu3OQrwLNdi+GDwD4EXb/dV20aJHS0tI0d+5cy/ijzMxMzZ8/X/n5+Vq8eLHTigSuNsF+Pnpt6s3NrqQdHeqvTQ8kspI2ALQiu2+xdenSRStWrNCECROs2jds2KC5c+fq+++/d0qBrsAtNgAA3E+7uMVWW1urwYMHN2lPSEhQXR2r+wIAAPdld0C67777mn2kSHp6uiZOnOhQUQAAAK7UojFICxYssPxvk8lkeThtYmKiJCkrK0v5+fm6//77nVslAABAG2pRQDpw4IDV+4SEBEnSsWPHJEnh4eEKDw/X4cOHnVQeAABA22tRQNqxY0dr1QEAANBuOLSISnFxsdasWaMvv/xSkjRgwABNnTpVISEhTikOAADAFewepP3pp5+qT58+Wr58uYqKilRUVKS0tDT16dNH2dnZzqwRAACgTdm9DtKwYcN07bXXavXq1fL2vnghqq6uTtOnT9fXX3+tTz75xKmFtiXWQQIAwP20i2ex+fv768CBA00eSvvFF19o8ODBqqiocKgwVyIgAQDgftrFQpHBwcHKz89v0n7ixAl17NjRoaIAAABcye6AlJycrGnTpmnTpk06ceKETpw4oY0bN2r69OlNHj8CAADgTuyexfZf//VfMplMuv/++y2PFvHx8dGsWbP0u9/9zmkFAgAAtDW7xyBdUlFRYVkosk+fPgoICHBKYa7EGCQAANyPM7+/HVoHSZICAgJ0/fXXO3oYAACAdqNFY5CaG5T9Y7777rsW7Q8AANAetCgg3XTTTXrggQe0b9++y+5TUlKi1atXKy4uTps3b3a4QAAAgLbWoltsX3zxhZ5//nn927/9m/z8/JSQkKDo6Gj5+fnp/Pnz+uKLL3T48GENGjRIL7zwgn7xi1+0Vt0AAACtxq5B2pWVlXr//feVkZGhb7/9VpWVlQoPD9eNN96o0aNHKy4urjVqbTMM0gYAwP20i5W0PRkBCQAA99MuVtIGAADwVG4VkJ555hmZTCarV+NnwVVVVWn27Nnq3LmzgoKCNG7cOJ0+fdqFFQMAAHfkVgFJkgYMGKDCwkLLKyMjw7Jt/vz5eu+99/TWW29p165dKigo0L333uvCagEAgDuye6HIEydOqHv37s6s5Yp4e3srMjKySXtJSYnWrFmjN998Uz//+c8lSa+++qr69++vPXv2KDExsa1LBQAAbsruK0j9+vXTokWLVFFR4cx6/qWjR48qOjpavXv31sSJEy2LV+7fv1+1tbUaOXKkVY0xMTHKzMz80WNWV1ertLTU6gUAAK5edgekDz/8UB988IGuu+46rVu3zoklXd6QIUO0bt06bd26VStXrtTx48c1bNgwXbhwQadOnZKvr69CQ0OtfiYiIkKnTp360eMuW7ZMISEhlpcrrowBAID2w+Fp/n/5y1/0xBNPqGvXrvrDH/6gYcOGOau2f6m4uFg9evRQWlqa/P39NWXKFFVXV1vtc/PNN+tnP/uZ/vM///Oyx6murrb6udLSUnXv3p1p/gAAuJF2Nc3//vvvV25uru644w7dfvvt+vd//3cdP37c0cNekdDQUP3kJz9RXl6eIiMjVVNTo+LiYqt9Tp8+3eyYpcbMZrOCg4OtXgAA4OrltFlso0aN0vTp0/XOO+8oNjZWv/3tb1VWVuaswzerrKxMx44dU1RUlBISEuTj46Pt27dbtufm5io/P19JSUmtWgcAAPAsds9iW7Vqlfbt26d9+/bpyy+/lJeXl+Li4vSb3/xG8fHx2rhxo2JjY/X2229r8ODBTik2NTVVY8eOVY8ePVRQUKCnn35aHTp00IQJExQSEqJp06ZpwYIFCgsLU3BwsObOnaukpCRmsAEAgBaxOyA9//zzGjJkiO6//34lJiYqISFB/v7+lu0zZ87U0qVLNXnyZOXk5Dil2JMnT2rChAk6d+6cunTpoltvvVV79uxRly5dJEnLly+Xl5eXxo0bp+rqao0ePVovvfSSUz4bAABcPVr1WWynT59WdHS06uvrW+sjWgXPYgMAwP20q0HaP6Zr1676+OOPW/MjAAAAnK5VA5LJZNJtt93Wmh8BAADgdG73LDYAAIDWRkACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACwQUACAACw4bEB6c9//rN69uwpPz8/DRkyRHv37nV1SQAAwE14ZEDatGmTFixYoKefflrZ2dmKj4/X6NGjdebMGVeXBgAA3IBHBqS0tDTNmDFDU6ZMUWxsrFatWqWAgACtXbvW1aUBAAA34O3qApytpqZG+/fv12OPPWZp8/Ly0siRI5WZmdnsz1RXV6u6utryvqSkRJJUWlrausUCAACnufS9bRiGw8fyuID0/fffq76+XhEREVbtERER+uqrr5r9mWXLlunZZ59t0t69e/dWqREAALSec+fOKSQkxKFjeFxAssdjjz2mBQsWWN4XFxerR48eys/Pd/gX3N6Ulpaqe/fuOnHihIKDg11djlPRN/dE39wTfXNfnty/kpISxcTEKCwszOFjeVxACg8PV4cOHXT69Gmr9tOnTysyMrLZnzGbzTKbzU3aQ0JCPO4/nkuCg4Ppmxuib+6JvrknT+6b5Nn98/JyfIi1xw3S9vX1VUJCgrZv325pa2ho0Pbt25WUlOTCygAAgLvwuCtIkrRgwQKlpKRo8ODBuvnmm/WHP/xB5eXlmjJliqtLAwAAbsAjA1JycrLOnj2rRYsW6dSpU7rhhhu0devWJgO3L8dsNuvpp59u9rabu6Nv7om+uSf65p48uW+SZ/fPmX0zGc6YCwcAAOBBPG4MEgAAgKMISAAAADYISAAAADYISAAAADYISDb+/Oc/q2fPnvLz89OQIUO0d+9eV5dkl08++URjx45VdHS0TCaT3n33XavthmFo0aJFioqKkr+/v0aOHKmjR4+6ptgWWLZsmW666SZ17NhRXbt21d13363c3FyrfaqqqjR79mx17txZQUFBGjduXJOFQ9ujlStXauDAgZbF25KSkvSPf/zDst1d+9Wc3/3udzKZTHrooYcsbe7cv2eeeUYmk8nq1a9fP8t2d+6bJH333XeaNGmSOnfuLH9/f11//fX69NNPLdvd9e9Jz549m5w3k8mk2bNnS3Lv81ZfX6+nnnpKvXr1kr+/v/r06aPnnnvO6hll7nreJOnChQt66KGH1KNHD/n7++uWW27Rvn37LNud0jcDFhs3bjR8fX2NtWvXGocPHzZmzJhhhIaGGqdPn3Z1aS3297//3XjiiSeMt99+25BkvPPOO1bbf/e73xkhISHGu+++axw6dMj45S9/afTq1cuorKx0TcFXaPTo0carr75q5OTkGAcPHjR+8YtfGDExMUZZWZlln9/85jdG9+7dje3btxuffvqpkZiYaNxyyy0urPrK/M///I/x/vvvG0eOHDFyc3ONxx9/3PDx8TFycnIMw3Dfftnau3ev0bNnT2PgwIHGvHnzLO3u3L+nn37aGDBggFFYWGh5nT171rLdnftWVFRk9OjRw5g8ebKRlZVlfP3118YHH3xg5OXlWfZx178nZ86csTpnH374oSHJ2LFjh2EY7n3enn/+eaNz587Gli1bjOPHjxtvvfWWERQUZPzxj3+07OOu580wDOM//uM/jNjYWGPXrl3G0aNHjaefftoIDg42Tp48aRiGc/pGQGrk5ptvNmbPnm15X19fb0RHRxvLli1zYVWOsw1IDQ0NRmRkpPHiiy9a2oqLiw2z2Wxs2LDBBRXa78yZM4YkY9euXYZhXOyHj4+P8dZbb1n2+fLLLw1JRmZmpqvKtFunTp2MV155xWP6deHCBeO6664zPvzwQ+O2226zBCR379/TTz9txMfHN7vN3fv26KOPGrfeeutlt3vS35N58+YZffr0MRoaGtz+vN1xxx3G1KlTrdruvfdeY+LEiYZhuPd5q6ioMDp06GBs2bLFqn3QoEHGE0884bS+cYvt/9TU1Gj//v0aOXKkpc3Ly0sjR45UZmamCytzvuPHj+vUqVNWfQ0JCdGQIUPcrq8lJSWSZHkw4f79+1VbW2vVt379+ikmJsat+lZfX6+NGzeqvLxcSUlJHtOv2bNn64477rDqh+QZ5+3o0aOKjo5W7969NXHiROXn50ty/779z//8jwYPHqxf/epX6tq1q2688UatXr3ast1T/p7U1NTojTfe0NSpU2Uymdz+vN1yyy3avn27jhw5Ikk6dOiQMjIydPvtt0ty7/NWV1en+vp6+fn5WbX7+/srIyPDaX3zyJW07fH999+rvr6+yWrbERER+uqrr1xUVes4deqUJDXb10vb3EFDQ4MeeughDR06VHFxcZIu9s3X11ehoaFW+7pL3z7//HMlJSWpqqpKQUFBeueddxQbG6uDBw+6db8kaePGjcrOzrYaJ3CJu5+3IUOGaN26derbt68KCwv17LPPatiwYcrJyXH7vn399ddauXKlFixYoMcff1z79u3Tgw8+KF9fX6WkpHjM35N3331XxcXFmjx5siT3/29y4cKFKi0tVb9+/dShQwfV19fr+eef18SJEyW59/dAx44dlZSUpOeee079+/dXRESENmzYoMzMTF177bVO6xsBCW5r9uzZysnJUUZGhqtLcZq+ffvq4MGDKikp0d/+9jelpKRo165dri7LYSdOnNC8efP04YcfNvlXnye49K9ySRo4cKCGDBmiHj166K9//av8/f1dWJnjGhoaNHjwYC1dulSSdOONNyonJ0erVq1SSkqKi6tznjVr1uj2229XdHS0q0txir/+9a9av3693nzzTQ0YMEAHDx7UQw89pOjoaI84b6+//rqmTp2qbt26qUOHDho0aJAmTJig/fv3O+0zuMX2f8LDw9WhQ4cmMxROnz6tyMhIF1XVOi71x537OmfOHG3ZskU7duzQNddcY2mPjIxUTU2NiouLrfZ3l775+vrq2muvVUJCgpYtW6b4+Hj98Y9/dPt+7d+/X2fOnNGgQYPk7e0tb29v7dq1SytWrJC3t7ciIiLcun+2QkND9ZOf/ER5eXluf+6ioqIUGxtr1da/f3/LLURP+Hvy7bff6qOPPtL06dMtbe5+3h555BEtXLhQ48eP1/XXX6/77rtP8+fP17JlyyS5/3nr06ePdu3apbKyMp04cUJ79+5VbW2tevfu7bS+EZD+j6+vrxISErR9+3ZLW0NDg7Zv366kpCQXVuZ8vXr1UmRkpFVfS0tLlZWV1e77ahiG5syZo3feeUcff/yxevXqZbU9ISFBPj4+Vn3Lzc1Vfn5+u+9bcxoaGlRdXe32/RoxYoQ+//xzHTx40PIaPHiwJk6caPnf7tw/W2VlZTp27JiioqLc/twNHTq0yVIaR44cUY8ePSS599+TS1599VV17dpVd9xxh6XN3c9bRUWFvLysv+I7dOighoYGSZ5x3iQpMDBQUVFROn/+vD744APdddddzuubs0aVe4KNGzcaZrPZWLdunfHFF18YM2fONEJDQ41Tp065urQWu3DhgnHgwAHjwIEDhiQjLS3NOHDggPHtt98ahnFxCmRoaKjx3//938Znn31m3HXXXW4xvXPWrFlGSEiIsXPnTqvpuRUVFZZ9fvOb3xgxMTHGxx9/bHz66adGUlKSkZSU5MKqr8zChQuNXbt2GcePHzc+++wzY+HChYbJZDK2bdtmGIb79utyGs9iMwz37t/DDz9s7Ny50zh+/Ljxz3/+0xg5cqQRHh5unDlzxjAM9+7b3r17DW9vb+P55583jh49aqxfv94ICAgw3njjDcs+7vr3xDAuzlaOiYkxHn300Sbb3Pm8paSkGN26dbNM83/77beN8PBw47e//a1lH3c+b1u3bjX+8Y9/GF9//bWxbds2Iz4+3hgyZIhRU1NjGIZz+kZAsvGnP/3JiImJMXx9fY2bb77Z2LNnj6tLssuOHTsMSU1eKSkphmFcnOL51FNPGREREYbZbDZGjBhh5ObmurboK9BcnyQZr776qmWfyspK4//9v/9ndOrUyQgICDDuueceo7Cw0HVFX6GpU6caPXr0MHx9fY0uXboYI0aMsIQjw3Dffl2ObUBy5/4lJycbUVFRhq+vr9GtWzcjOTnZap0gd+6bYRjGe++9Z8TFxRlms9no16+fkZ6ebrXdXf+eGIZhfPDBB4akZut15/NWWlpqzJs3z4iJiTH8/PyM3r17G0888YRRXV1t2cedz9umTZuM3r17G76+vkZkZKQxe/Zso7i42LLdGX0zGUajZTUBAADAGCQAAABbBCQAAAAbBCQAAAAbBCQAAAAbBCQAAAAbBCQAAAAbBCQAAAAbBCQAAAAbBCQAAAAbBCQAAAAbBCTAA507d05du3bVN998Y2n76U9/qoceesjhYzvrOG3BnWq9Uvb06Up+pqXHHT9+vH7/+9+3qA7AnXi7ugAAzvf888/rrrvuUs+ePS1tb7/9tnx8fFxXlIf46U9/qhtuuEF/+MMfXPL57eU8Pvnkkxo+fLimT5+ukJAQV5cDOB1XkAAPU1FRoTVr1mjatGlW7WFhYerYsaOLqoIz1NTUtJvzGBcXpz59+uiNN95wdSlAqyAgAe3chg0b5O/vr8LCQkvblClTNHDgQJWUlDTZ/+9//7vMZrMSExOt2m1vofz0pz/Vgw8+qN/+9rcKCwtTZGSknnnmGaufKS8v1/3336+goCBFRUU1uaXS0NCgZcuWqVevXvL391d8fLz+9re/NfncOXPmaM6cOQoJCVF4eLieeuopGYZxxcdxRq3O+JzJkydr165d+uMf/yiTySSTyWR1G/OS9PR0RUdHq6Ghwar9rrvu0tSpUyVJW7du1a233qrQ0FB17txZd955p44dO9bs7+6hhx5SeHi4Ro8e3eQ8XslxJKmuru5Hz0FLfk+XjB07Vhs3bmz2GIDbMwC0aw0NDcbAgQONOXPmGIZhGIsWLTKuueYa4+TJk83u/+CDDxpjxoxp0n7bbbcZ8+bNs3ofHBxsPPPMM8aRI0eM1157zTCZTMa2bdss+8yaNcuIiYkxPvroI+Ozzz4z7rzzTqNjx46W4yxZssTo16+fsXXrVuPYsWPGq6++apjNZmPnzp1WnxMUFGTMmzfP+Oqrr4w33njDCAgIMNLT0y37/KvjOKNWZ3xOcXGxkZSUZMyYMcMoLCw0CgsLjbq6uia/66KiIsPX19f46KOPLG3nzp2zavvb3/5mbN682Th69Khx4MABY+zYscb1119v1NfXN/ndPfLII8ZXX31lfPXVV03OY0uO82PnoPFxr+S8GoZh/OMf/zB8fX2NqqqqJr8DwN0RkAA38N577xlms9lYsmSJ0alTJyMnJ+ey+951113G1KlTm7Q3F5BuvfVWq31uuukm49FHHzUMwzAuXLhg+Pr6Gn/9618t28+dO2f4+/sb8+bNM6qqqoyAgADjf//3f62OMW3aNGPChAlWn9O/f3+joaHB0vboo48a/fv3NwzDuKLjOFqrsz6nud/j5dieh5dfftmIjo62Ci6NnT171pBkfP7551afdeONN1rt968+/3LH+bFz0Pi4V3peDcMwDh06ZEgyvvnmm8vWA7grBmkDbuDOO+9UbGysFi9erG3btmnAgAGX3beyslJ+fn5XdNyBAwdavY+KitKZM2ckSceOHVNNTY2GDBli2R4WFqa+fftKkvLy8lRRUaF/+7d/szpGTU2NbrzxRqu2xMREmUwmy/ukpCT9/ve/V319/RUfx5FaW1Lvj31OS0ycOFEzZszQSy+9JLPZrPXr12v8+PHy8ro4suHo0aNatGiRsrKy9P3331tux+Xn5ysuLs5ynISEhB/9nCs9zo+dgw4dOljaW3Je/f39JV0c9wZ4GgIS4Aa2bt2qr776SvX19YqIiPjRfcPDw3X+/PkrOq7tbCiTydRk3MzllJWVSZLef/99devWzWqb2Wy+omO05DiO1NqWn3PJ2LFjZRiG3n//fd10003avXu3li9fbrW9R48eWr16tWW8UlxcnGpqaqyOExgY+C8/50qOc6Vacl6LiookSV26dLHrs4D2jIAEtHPZ2dn6j//4D61Zs0br1q3TU089pbfeeuuy+994441OmVnUp08f+fj4KCsrSzExMZKk8+fP68iRI7rtttsUGxsrs9ms/Px83XbbbT96rKysLKv3e/bs0XXXXacOHTq06Dj21irJKZ8jSb6+vqqvr/+X+/n5+enee+/V+vXrlZeXp759+2rQoEGSLq5TlZubq9WrV2vYsGGSpIyMjBbX0pLj/Ng5aKwlv6ecnBxdc801Cg8Pb3HtQHtHQALasW+++UZ33HGHHn/8cU2YMEG9e/dWUlKSsrOzLV+2tkaPHq3HHntM58+fV6dOnez+7KCgIE2bNk2PPPKIOnfurK5du+qJJ56w3CLq2LGjUlNTNX/+fDU0NOjWW29VSUmJ/vnPfyo4OFgpKSmWY+Xn52vBggV64IEHlJ2drT/96U+WWWYtOY69tTrrcySpZ8+eysrK0jfffKOgoCCFhYVZfU5jEydO1J133qnDhw9r0qRJlvZOnTqpc+fOSk9PV1RUlPLz87Vw4cIr+vzGWnKcHzsHjbXk97R7926NGjWqxXUD7oCABLRTRUVFGjNmjO666y7Ll96QIUN0++236/HHH9fWrVub/bnrr79egwYN0l//+lc98MADDtXw4osvqqysTGPHjlXHjh318MMPWy0t8Nxzz6lLly5atmyZvv76a4WGhmrQoEF6/PHHrY5z//33q7KyUjfffLM6dOigefPmaebMmS0+jiO1OutzUlNTlZKSotjYWFVWVur48eNWC3I29vOf/1xhYWHKzc3Vr3/9a0u7l5eXNm7cqAcffFBxcXHq27evVqxYoZ/+9KdXXEdLj/OvzkFjV/J7qqqq0rvvvnvZ/w4Bd2cyjMsshAHAbb3//vt65JFHlJOTc9mrG23F1StPo3WsXLlS77zzjrZt2+bqUoBWwRUkwAPdcccdOnr0qL777jt1797d1eXAA/n4+OhPf/qTq8sAWg0BCfBQnvaQVrQv06dPd3UJQKviFhsAAIANnsUGAABgg4AEAABgg4AEAABgg4AEAABgg4AEAABgg4AEAABgg4AEAABgg4AEAABgg4AEAABgg4AEAABg4/8D42LYm+NTbOAAAAAASUVORK5CYII=\n",
|
|
"text/plain": [
|
|
"<Figure size 640x480 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"import matplotlib.pyplot as plt\n",
|
|
"plt.scatter(x, y, marker=\"x\")\n",
|
|
"plt.xlim([0, max(x)+10])\n",
|
|
"plt.ylim([0, max(y)+10])\n",
|
|
"plt.xlabel(\"$x$ (independent variable)\")\n",
|
|
"plt.ylabel(\"$y$ (observed dependent variable)\")\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Now that we have an idea what our data looks like, let's try fitting a simple linear regression model to it. We're going to use the LinearRegression model from scikit-learn (https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Model intercept: -20.000000000000057\n",
|
|
"Model coefficient: 3.000000000000001\n",
|
|
"Model predictions:\n",
|
|
" [[ 70.]\n",
|
|
" [100.]\n",
|
|
" [100.]\n",
|
|
" [130.]\n",
|
|
" [130.]\n",
|
|
" [130.]\n",
|
|
" [160.]\n",
|
|
" [190.]\n",
|
|
" [190.]\n",
|
|
" [220.]\n",
|
|
" [220.]\n",
|
|
" [220.]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"from sklearn.linear_model import LinearRegression\n",
|
|
"model = LinearRegression() # create a new LinearRegression object that will hold our model\n",
|
|
"model.fit(x, y) # the .fit() method calculates the model parameters\n",
|
|
"theta_0 = model.intercept_[0] # the intercept accounts for constant effects, \\theta_0 in Eqn. 1 above\n",
|
|
"theta_1 = model.coef_[0][0] # this is the \"slope\", \\theta_1 in Eqn. 2 above\n",
|
|
"predictions = model.predict(x)\n",
|
|
"print(\"Model intercept:\", theta_0) \n",
|
|
"print(\"Model coefficient:\", theta_1) \n",
|
|
"print(\"Model predictions:\\n\", predictions) "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Next let's plot the linear regression model, along with the original values and the values predicted by our model."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 37,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+40lEQVR4nO3deZxN9f/A8debCFkrS9/IqKTse7JlnRnbZIkUhSyliIooRXzT8lNfUdZEFAnZIkwYDTMtliQhWoaIjH1rGDPv3x/3zDSY4d6ZuXPnzryfj8d9zD333nPO+56Zue/7OZ/PeX9EVTHGGGOuJYevAzDGGOMfLGEYY4xxiyUMY4wxbrGEYYwxxi2WMIwxxrjlOl8HkBo333yzBgQE+DoMY4zxK5s3bz6iqkVTu75fJoyAgAA2bdrk6zCMMcaviMjetKxvp6SMMca4xRKGMcYYt1jCMMYY4xa/7MNITmxsLPv37ycmJsbXoZh0lCdPHkqWLEmuXLl8HYox2V6WSRj79++nQIECBAQEICK+DsekA1Xl6NGj7N+/nzJlyvg6HGOyvSxzSiomJoabbrrJkkUWIiLcdNNN1mo0JpPIMgkDsGSRBdnv1JjMI0slDGOMMcnbvn17mrdhCSMd5cyZk6pVq1KxYkXatGnDiRMnfB1SouHDh7N69eo0b2fdunW0bt0agKVLl/Lmm2+meZvGGO+aPXs21apVS/N2LGGko7x587J161a2b9/OjTfeyIQJE9K8zYsXL6ZDZDBq1CiaNWuWLttKEBISwtChQ9N1m8aY9DV27Fi6du1K/fr107wtSxhect9993HgwAEAfvvtN4KDg6lRowYNGjRg165diY/XqVOHSpUq8fLLL5M/f37A9S2+QYMGhISEUL58eeLi4hg8eDC1atWicuXKTJkyBYCDBw/SsGHDxFbN+vXriYuLo3v37lSsWJFKlSoxduxYALp3786CBQsAWLNmDdWqVaNSpUo8/vjjnD9/HnCVXBkxYgTVq1enUqVKiXGm5KOPPqJfv36J23/mmWeoW7cut99+e+K+AMaMGZMY+4gRI9LrEBtjrkJVGTJkCM899xwPPvggK1asSPM2s8yw2qQGDhzI1q1b03WbVatW5d1333XrtXFxcaxZs4aePXsC0KdPHyZPnkzZsmX57rvveOqpp1i7di0DBgxgwIABPPzww0yePPmSbWzZsoXt27dTpkwZpk6dSqFChdi4cSPnz5+nXr16BAYGsnDhQoKCghg2bBhxcXGcO3eOrVu3cuDAgcTzlZefFouJiaF79+6sWbOGu+66i8cee4xJkyYxcOBAAG6++Wa2bNnCxIkTefvtt5k2bZrbx+jgwYNs2LCBXbt2ERISwoMPPkhoaCh79uzh+++/R1UJCQkhPDychg0bur1dY4xnYmNj6d27NzNnzqRv376899575MyZM83btRZGOvrnn3+oWrUqJUqU4O+//6Z58+acOXOGyMhIOnbsSNWqVXniiSc4ePAgAN988w0dO3YE4JFHHrlkW7Vr10689iA0NJRZs2ZRtWpV7r33Xo4ePcqePXuoVasWM2bM4NVXX+Wnn36iQIEC3H777fz+++/079+flStXUrBgwUu2+8svv1CmTBnuuusuALp160Z4eHji8+3btwegRo0aREVFefT+27ZtS44cOShfvjx///13YuyhoaFUq1aN6tWrs2vXLvbs2ePRdo0x7jt37hzt2rVj5syZjBw5kgkTJqRLsoAs2sJwtyWQ3hL6MM6dO0dQUBATJkyge/fuFC5c2OMWzw033JB4X1V57733CAoKuuJ14eHhLF++nO7du/Pcc8/x2GOP8eOPP7Jq1SomT57MvHnzmD59utv7vf766wFXB76n/ScJ6ybEnPDzxRdf5IknnvBoW8YYzx07dozWrVvz3XffMXny5HT/v7MWhhfky5eP8ePH884775AvXz7KlCnD/PnzAdcH6I8//ghAnTp1+PzzzwGYO3duitsLCgpi0qRJxMbGArB7927Onj3L3r17KV68OL1796ZXr15s2bKFI0eOEB8fT4cOHXjttdfYsmXLJdsqV64cUVFR/PrrrwB8/PHH3H///el+DJLGPn36dM6cOQPAgQMHOHz4sNf2Z0x29eeff1K/fn02b97M/PnzvfIlLUu2MDKDatWqUblyZT799FNmz55N3759ee2114iNjaVz585UqVKFd999l65duzJ69GiCg4MpVKhQstvq1asXUVFRVK9eHVWlaNGiLF68mHXr1jFmzBhy5cpF/vz5mTVrFgcOHKBHjx7Ex8cD8MYbb1yyrTx58jBjxgw6duzIxYsXqVWrFk8++aTXjkNgYCA7d+7kvvvuAyB//vx88sknFCtWzGv7NCa72blzJ4GBgZw6dYpVq1bRqFEjr+xHEk4deJuIlAM+S/LQ7cBwYJbzeAAQBXRS1eNX21bNmjX18gmUdu7cyT333JOOEXvfuXPnyJs3LyLC3Llz+fTTT1myZImvw8p0/PF3a0xG+eabb2jdujW5c+dm5cqVVKlSJcXXishmVa2Z2n1lWAtDVX8BqgKISE7gALAIGAqsUdU3RWSoszwko+Lypc2bN9OvXz9UlcKFC3vU12CMMcuXL6djx47ceuutrFq1ittvv92r+/PVKammwG+quldEHgAaOY/PBNaRTRJGgwYNEvszjDHGE7NmzeLxxx+natWqfPnllxlymtdXnd6dgU+d+8VV9aBz/xBQPLkVRKSPiGwSkU3R0dEZEaMxxmRKY8aMoVu3bjRq1IiwsLAM6xPM8IQhIrmBEGD+5c+pq0Ml2U4VVZ2qqjVVtWbRokW9HKUxxmQ+8fHxDBo0iBdeeIGHHnqI5cuXU6BAgQzbvy9aGC2ALar6t7P8t4jcAuD8tDGXxhhzmdjYWLp168Y777xD//79mTNnziXXPmUEXySMh/n3dBTAUqCbc78bYMOEjDHZyuWjVS9fPnv2LCEhIXzyySeMHj2acePGkSNHxn98Z+geReQGoDmwMMnDbwLNRWQP0MxZ9kv79+/ngQceoGzZstxxxx0MGDCACxcuAJcW6stMEgoeZpbtGJPdjP1qN6OW7bikOsKoZTsY+9VuAI4cOUKTJk0IDQ3lgw8+4KWXXvLZxGIZmjBU9ayq3qSqJ5M8dlRVm6pqWVVtpqrHMiiWqy6nZnvt27enbdu27Nmzh927d3PmzBmGDRuWpu1eTXqVPjfG+IaqciomlhkRUYlJY9SyHcyIiOJUTCxRUVHUr1+fbdu2sXDhQnr16uXTeLNlaZBrZfTUWLt2LXny5KFHjx6AqxbT2LFjmT59OufOnQNcl+43atSIsmXLMnLkSMDV1GzVqhVVqlShYsWKfPaZ69rGzZs3c//991OjRg2CgoISCxY2atSIgQMHUrNmTUaPHk3p0qUTr+o+e/YspUqVIjY2NsWS6n/88Qf33XdfYkn15AwdOvSSuTxeffVV3n77bc6cOUPTpk0Ty58nd5Fh0gmWAPr168dHH3101fc0fvx4ypcvT+XKlencuXPqfgHG+CERYXjr8vSoF8CMiCjKvPglMyKi6FEvgPYB8dSvX59Dhw4RGhrKAw884OtwXR+W/narUaOGXm7Hjh1XPJac+Ph4fXXpdi09ZJm+unR7ssupMW7cOB04cOAVj1etWlV//PFHnTFjhpYoUUKPHDmi586d0woVKujGjRt1wYIF2qtXr8TXnzhxQi9cuKD33XefHj58WFVV586dqz169FBV1fvvv1/79u2b+PqQkBBdu3Zt4ut69uypqqpNmjTR3bt3q6rqt99+q40bN1ZV1TZt2ujMmTNVVfX999/XG2644YqYt2zZog0bNkxcvueee3Tfvn0aGxurJ0+eVFXV6OhoveOOOxKPV8J2wsLCtFWrVonrPv300zpjxoyrvqdbbrlFY2JiVFX1+PHjV8Tj7u/WGH8VHx+vpYcsS7yFh4dr4cKF9ZZbbtFt27al236ATZqGz16PL9xz+iFiVDUu/dOX9yVkdIAZEVHMiIgCoEe9AIa3Lu/Vc4PNmzfnpptuAlxlxDds2EDLli15/vnnGTJkCK1bt6ZBgwZs376d7du307x5c8A1v8Ytt9ySuJ2HHnrokvufffYZjRs3Zu7cuTz11FOXlFRPkDBJUkRERGLBw0cffZQhQ668RrJatWocPnyYv/76i+joaIoUKZLYcnnppZcIDw8nR44cHDhwgL///psSJUpc873/8ssvKb6nypUr06VLF9q2bUvbtm09OaTG+D1VZeGg/2PDtHf4z6kjzMxbgGZvnyPgjttZtWoVAQEBvg4x0TUThojkwHWhXRegFnAeuF5EjgDLgSmq+qtXo0xnCUkjIVkAaU4W5cuXv2SWOYBTp06xb98+7rzzTrZs2XLF9kWEu+66iy1btvDll1/y8ssv07RpU9q1a0eFChX45ptvkt1X0tLnISEhvPTSSxw7dozNmzfTpEkTzp49e9WS6u68z44dO7JgwQIOHTqUmKBmz55NdHQ0mzdvJleuXAQEBBATE3PJetddd13iKTIg8XlVTfE9LV++nPDwcL744gtGjx7NTz/9xHXXWV1Mk/UlJIsW40eQ7+J5pgN9/jlFVYSnmnSidOnSvg7xEu70YYQBdwAvAiVUtZSqFgPqA98Cb4lIVy/GmO7U6bNIKmmfRmo0bdqUc+fOMWvWLMD1Dfr555+ne/fu5MuXD4CvvvqKY8eO8c8//7B48WLq1avHX3/9Rb58+ejatSuDBw9my5YtlCtXjujo6MQP19jYWH7++edk95s/f35q1arFgAEDaN26NTlz5qRgwYIpllSvV69eYin12bNnp/h+HnroIebOncuCBQsSWyonT56kWLFi5MqVi7CwMPbu3XvFeqVLl2bHjh2cP3+eEydOsGbNGoAU31N8fDx//vknjRs35q233uLkyZOJpdCNyepEhOafjCfvxfO8AfTEVTdpHcqDC6f7bDRUStxJGM1U9b+quk1VE786quoxVf1cVTtwaRXaTC0hWSR0LP3xRsvEDqe0JA0RYdGiRcyfP5+yZcty1113kSdPHl5//fXE19SuXZsOHTpQuXJlOnToQM2aNfnpp5+oXbs2VatWZeTIkbz88svkzp2bBQsWMGTIEKpUqULVqlWJjIxMcd8PPfQQn3zyySWnqmbPns2HH35IlSpVqFChQmIH9bhx45gwYQKVKlVKnHM8ORUqVOD06dPceuutiaeOunTpwqZNm6hUqRKzZs3i7rvvvmK9UqVK0alTJypWrEinTp2oVq0aQIrvKS4ujq5du1KpUiWqVavGM888Q+HChT069sb4s/yH/+JZ4CXgEeALID9QMPrgVdfzBbfLm4sr1XUBblfVUSJyG64Wx/feDDA5aS1vPvar3ZyKiU08DZWQRArmycWzze/yRsgmDay8ucmqLly4QPciRfj03DkGAu+Q5Ft86dLg4TTJ15KR5c0nAvFAE2AUcBr4HFe/hl95tvldqGpicy+hTyOzNf+MMVnX6dOn6dChA1+dO8ebuXLxQmwsiZ9A+fLB6NG+DC9ZnlyHca+qPg3EAKhrkqPcXokqAyTXAW2MMRkhOjqaJk2asHbtWqZPn86QGTOQ0qVBxNWymDoVunTxdZhX8KSFEetMfKQAIlIUV4vDGGOMm6KioggMDGT//v0sXrz43wtdM2GCuJwnCWM8rhnyionIaOBBIPlLhY0xxlxh27ZtBAcHExMTw+rVq6lbt66vQ/KI2wlDVWeLyGZco74EaKuqO70WmTHGZCHh4eGEhISQP39+1q9fT4UKFXwdksc8ujpKVXcBu7wUizHGZEmLFy+mc+fOlClThlWrVnHbbbf5OqRUuWant4icFpFTzu2K+xkRpL/ImTMnVatWpWLFinTs2DGx6GBqdO/ePfHK8V69erFjx44UX7tu3bqrXqeRkoCAAI4cOZLqGNN7O8ZkRR988AEdOnSgatWqbNiwwW+TBbiRMFS1gKoWdG5X3M+IIP1F3rx52bp1K9u3byd37txMnjz5kudTW4582rRplC9fPsXnU5swjDHeo6q89tpr9OnTh6CgINasWZNYS85fuT2sVkTyiMhzIrJQRD4XkYEiksebwXnV7NkQEAA5crh+XqVMRmo0aNCAX3/9lXXr1tGgQQNCQkIoX748cXFxDB48mFq1alG5cmWmTJkCuP64+vXrR7ly5WjWrBmHD/87U22jRo1IuFBx5cqVVK9enSpVqtC0aVOioqKYPHkyY8eOpWrVqqxfv57o6Gg6dOhArVq1qFWrFhEREQAcPXqUwMBAKlSoQK9evZK9qn3y5MkMHjw4cTnpxE9t27alRo0aVKhQgalTp16xblRUFBUrVkxcfvvtt3n11VcBUiy3Pn/+fCpWrEiVKlVo2LBhWg65MZlGXFwc/fv355VXXuHRRx9lyZIll9SA81vulrUF5gEfAo2d2wfA/LSUyk3tLS3lzVVV9ZNPVPPlU4V/b/nyuR5Pg4QS37GxsRoSEqITJ07UsLAwzZcvn/7++++qqjplyhT973//q6qqMTExWqNGDf3999/1888/12bNmunFixf1wIEDWqhQIZ0/f76qukqab9y4UQ8fPqwlS5ZM3NbRo0dVVXXEiBE6ZsyYxDgefvhhXb9+vaqq7t27V++++25VVe3fv7+OHDlSVVWXLVumgEZHR1/yHg4fPqx33HFH4nJwcHDithL2l1Ce/ciRI6qqWrp0aY2OjtY//vhDK1SokLjumDFjdMSIEaqacrn1ihUr6v79+1U1+dLmqlbe3PiXmJgY7dSpkwI6aNAgjYuL83VIicjA8uYVVTXpeZEwEUn5xHpmNmwYXN6/cO6c6/E0jIX+559/qFq1KuBqYfTs2ZPIyEhq165NmTJlAAgNDWXbtm2J/RMnT55kz549hIeH8/DDD5MzZ07+85//0KRJkyu2/+2339KwYcPEbd14443JxrF69epL+jxOnTrFmTNnCA8PZ+FC1+y4rVq1okiRIlesW7RoUW6//Xa+/fZbypYty65du6hXrx7gmuho0aJFgGsyqD179rjVxL5aufV69erRvXt3OnXqRPv27a+5LWMys1OnTtGuXTvWrl3LmDFjGDRokK9DSleeJIwtIlJHVb8FEJF7gU3XWCdz2rfPs8fdlNCHcbmkTVFV5b333iMoKOiS13z55Zdp2ndS8fHxfPvtt+TJk7ozhp07d2bevHncfffdtGvXDhFh3bp1rF69mm+++YZ8+fLRqFEjt0ubx8fHp1huffLkyXz33XcsX76cGjVqsHnzZr8/z2uyp7///puWLVuybds2Zs2axaOPPurrkNKdO6OkfhKRbUANIFJEokTkD+AbwKMiViJSWEQWiMguEdkpIveJyI0i8pWI7HF+Xvm1N72lNEohA0YvBAUFMWnSJGJjYwHYvXs3Z8+epWHDhnz22WfExcVx8OBBwsLCrli3Tp06hIeH88cffwBw7Jhr+vMCBQpw+vTpxNcFBgby3nvvJS4nfFA3bNiQOXPmALBixQqOHz+ebIzt2rVjyZIlfPrpp4lTpp48eZIiRYqQL18+du3axbfffnvFesWLF+fw4cMcPXqU8+fPs2zZMoCrllv/7bffuPfeexk1ahRFixblzz//dPNIGpN5/P7779SrV49du3axdOnSLJkswL1O79ZAGyAYKAPcDzRy7rfwcH/jgJWqejdQBdgJDAXWqGpZYI2z7F2jR7uKeyWVQcW+evXqRfny5alevToVK1bkiSee4OLFi7Rr146yZctSvnx5HnvsMe67774r1i1atChTp06lffv2VKlSJbGceZs2bVi0aFFip/f48ePZtGkTlStXpnz58omjtUaMGEF4eDgVKlRg4cKFKQ7vK1KkCPfccw979+6ldu3aAAQHB3Px4kXuuecehg4dSp06da5YL1euXAwfPpzatWvTvHnzS8qfp1RuffDgwVSqVImKFStSt25dqlSpkrYDbEwG27p1K3Xr1uX48eOsWbOGFi08/Vj0I2npAPHkBhQC/sApqZ7k8V+AW5z7twC/XGtbae70VnV1cJcurSri+pnGDm/jPdbpbTKrtWvXaoECBbRUqVJ+8XdKRs7p7ZwuKgsknhxX1XA3Vy8DRAMzRKQKsBkYABRX1YSZQg4BxVPYdx+gD5A+F7506eIXxb6MMZnTggUL6NKlC3feeSerVq2iZMmSvg7J6zy5DqMXEA6sAkY6P1/1YF/XAdWBSapaDTjLZaefnAyY7IxOqjpVVWuqas2iRYt6sFtjjElfkydPplOnTtSsWZP169dni2QBns2HMQDXZEl7VbUxUA044cH6+4H9qvqds7wAVwL5W0RuAXB+Hk5h/WvSNMzJbTIn+52azERVefXVV+nbty+tWrXiq6++SnF4e1bkScKIUdUYABG5Xl2FCMu5u7KqHgL+FJGEdZoCO4ClQDfnsW7AEg9iSpQnTx6OHj1qHzBZiKpy9OjRVA8PNiY9xcXF8dRTTzFy5Eh69OjBokWLyHf54JkszpM+jP0iUhhYDKwWkWPAXg/31x+YLSK5gd+BHriS1jwR6elsr5OH2wSgZMmS7N+/n+jo6NSsbjKpPHnyZJvmvsm8YmJi6NKlCwsXLmTo0KG8/vrr2XKWTknNN3IRuR8oiGuIbGy6R3UNNWvW1ITaSsYY400nT56kbdu2rFu3jrFjxzJw4EBfh5RqIrJZVT26fi6pa7YwRGSDqtYXkdNc2iEtzrJVrDXGZEkHDx6kRYsW/Pzzz8yePZtHHnnE1yH51DUThpMsBKigqmmrnWGMMX7i119/JTAwkMOHD7N8+XICAwN9HZLPudXp7Qx3Xe7lWIwxJlPYvHkzdevW5fTp06xdu9aShcOTUVJbRKSW1yIxxphMYPXq1TRq1Ih8+fKxYcOGxPI4xrOEcS/wjYj8JiLbkhQlNMaYLGHevHm0bNmSMmXKEBkZSblybl85kC14Mqw26NovMcYY//T+++/zzDPPUL9+fZYuXUrhwoV9HVKm43bCUFVPr7kwxphMT1UZPnw4r732Gg888ACffvopefPm9XVYmVJGFh80xphM5eLFi/Tt25dp06bRq1cvJk2axHXXefSxmK24fWSc4oMDgJLAVqAOrkmUrpxL1BhjMrl//vmHRx55hMWLF/Pyyy8zatSobHn1ticysvigMcZkCidOnCAoKIglS5Ywfvx4/vvf/1qycIMnba8YVY0RkcTig0kKCRpjjF/466+/CA4OZteuXXz66aeJM1eaa0tt8cGvROQ4nhcfNMYYn9m9ezeBgYEcPXqUL7/8kmbNmvk6JL/iScLYANygqq+KSBiuKVdXeicsY4xJXxs3bqRly5aICOvWraNGjRq+DsnveNKHkR8IFZH1QCXgO1W94J2wjDEm/YSGhtK4cWMKFChARESEJYtUcjthqOpIVa0APA3cAnwtIqu9Fpkxxrjp8mkaki7PmTOHVq1aceeddxIREUHZsmUzOrwsIzUDjg8Dh4CjQLH0DccYYzwz9qvdnIqJZXjr8ogIqsqoZTsomCcXOXasYODAgdx///0sWbKEQoUK+Tpcv+bJdRhP4ZoNrygwH+itqju8FZgxxlyLqnIqJpYZEVEADG9dnlHLdjB9wx/c9ttiNnz+Ie3bt2f27Nk21W868KSFUQoYqKpbvRSLMcZ4REQY3ro8ADMiopgREYXGx3HjlhlsWLOYJ554ggkTJpAzZ04fR5o1eFJL6kVvBmKMMakhIgxfMYFhk6dwXuN5CNfkPSNGjGDEiBF2QV46ytCiKSISBZwG4oCLqlpTRG4EPgMCgCigk6oez8i4jDH+S596CiZN4hQQAkQCE4C+hw9bskhnngyrTS+NVbVqkonIhwJrVLUssMZZNsaYa1JV4qdM4QDQENgIzAOeAuKnTLli9JRJG18kjMs9AMx07s8E2vouFGOMPxERfomPpy6wD9eVxA86z+WIj7cWRjpzO2GIyFvuPHYNiuviv80i0sd5rLiqHnTuHwKKp7D/PiKySUQ2RUdHe7hbY0xW9N1339EAuAB8DTRO8pxYR3e686SF0TyZx1p4uL/6qlrdWe9pEWmY9El1tR+TbUOq6lRVramqNYsWLerhbo0xWc2KFSto0qQJRQoWJAJX+exL9OmTzFomLa6ZMESkr4j8BJRz5vJOuP0B/OTJzlT1gPPzMLAIqA38LSK3OPu6BdeFgcYYk6KPP/6YkJAQypUrR8Tu3dzRty8ktChy5oS+fWHiRN8GmQXJtTqFRKQQUAR4g0s7pE+r6jG3dyRyA5BDVU87978CRgFNgaOq+qaIDAVuVNUXrratmjVr6qZNm9zdtTEmC3nnnXcYNGgQTZo0YdGiRRQsWNDXIfkNEdmcZMCRx645rFZVTwInRaQH0B7X8NfrnJ2jqqPc3FdxYJHTCXUdMEdVV4rIRmCeiPTEVS69k8fvwhiT5cXHxzNkyBDefvttOnbsyMcff8z1118PuEZLJe3gvnzZpA9PrsNYDJwENgPnPd2Rqv4OVEnm8aO4WhnGGJOs2NhYevXqxaxZs3j66acZN25c4tXbV6sl9Wzzu3wcedbiScIoqarBXovEGGOScfbsWTp27MiKFSv473//y7BhwxJbDynVkpoREUWPegHW0khnniSMSBGppKoedXQbY0xqHT16lFatWrFx40amTp1K7969L3k+uVpSAD3qBSS2OEz68WRYbX1gi4j84oyS+klEtnkrMGNM9rZv3z7q16/P1q1b+fzzz69IFgmSJo0Eliy8w5MWhqfXXBhjTKr8/PPPBAUFcebMGUJDQ2nYsGGKr03os0hq1LIdljS8wJMWxj6gAdBNVffiusAu2auyjTEmtSIiIqhfvz7x8fGEh4e7lSwS+iz+eKMlPeoFMCMiilHLdlgtqXTmSQtjIhAPNMF1/cRp4HOglhfiMsZkQ1988QWdOnXitttuY9WqVQQEBFz19SJCwTy5LumzSDg9VTBPLmthpDNPEsa9qlpdRH4AUNXjIpLbS3EZY7KZGTNm0Lt3b6pVq8aXX36JuyWAnm1+1yWjoRKShiWL9OfJKalYEcmJU+tJRIrianEYY0yqqSpvvvkmjz/+OE2bNiUsLMztZJHg8uRgycI7PEkY43HVfyouIqOBDcDrXonKGJMtxMfH89xzz/Hiiy/y8MMP88UXX5A/f35fh2VS4MkUrbNFZDP/XpXdVlV3eicsY0xWd+HCBXr06MGcOXMYMGAA//vf/8iRIzNM0WNScs2EISLPpfBUCxFpoar/S+eYjDFZ3JkzZ+jQoQOhoaG88cYbDBkyxE4j+QF3WhgFnJ/lcI2IWuostwG+90ZQxpisKzo6mlatWrFlyxamT59Ojx49fB2ScZM71WpHAohIOFBdVU87y68Cy70anTEmS4mKiiIoKIh9+/axaNEi2rRp4+uQjAc8GVZbHNdMiAkuYBfuGWPctG3bNoKDg/nnn39YvXo19erV83VIxkOeJIxZwPcisshZbgt8lN4BGWOynvXr19OmTRvy58/Phg0bqFChgq9DMqng9pAEVR0N9ACOO7ceqvqGtwIzxmQNixcvpnnz5pQoUYLIyEhLFn7MkxYGqroF2OKlWIwxWcwHH3zAk08+Sa1atVi2bBk333yzr0MyaeB2whCR64EOJJmiFfBkilZjTDahqowePZpXXnmFFi1aMH/+fG644QZfh2XSyJMWxhLSMEWrMSZ7iIuLY8CAAUyYMIGuXbsyffp0cuXK5euwTDrI8ClanXpUm4ADqtpaRMoAc4GbcCWjR1X1wtW2YYzJnM6fP89jjz3GvHnzeP755/m///s/u3o7C/HkNxkpIpXSYZ8DgKQlRd4Cxqrqnbg603umwz6MMRns1KlTtGzZknnz5jFmzBjefvttSxZZTIZO0SoiJYFWwDRnWXDNr7HAeclMXMN1jTF+5O+//6Zx48Z8/fXXzJw5k0GDBvk6JOMFGT1F67vAC/xbbuQm4ISqXnSW9wO3JreiiPQB+gDcdttt6RCKMSY9/P777wQGBvLXX3+xdOlSWrZs6euQjJdk2BStItIaOKyqmz0L0UVVp6pqTVWt6WmtfGOMd2zdupW6dety/Phx1q5da8kii/MkYUwE7gMedpZPAxM8WL8eECIiUbg6uZsA44DCIpLQ0ikJHPBgm8YYHwkLC6Nhw4bkzp2bDRs2UKdOHV+HZLzMk4Rxr6o+DcSAa4pWwO0pWlX1RVUtqaoBQGdgrap2AcKAB52XdcM1fNcYk4ktWLCA4OBgSpUqRUREBPfcc4+vQzIZIDNM0ToEeE5EfsXVp/FhOmzTGOMlkydPplOnTtSsWZP169dTqlQpX4dkMognnd6XT9H6IPByanaqquuAdc7934HaqdmOMSbjqCojR45k5MiRtGrVinnz5pEvXz5fh2UykE3Raoy5pri4OPr168fkyZPp3r07U6dOtau3syFPaknlAVriGikVD+QWkT9UNcZbwRljfC8mJoYuXbqwcOFCXnjhBd58802bTjWb8nQ+jNO4Tk0BPAJ8DHRM76CMMZnDyZMnadu2LevWreN///sfzz77rK9DMj7kScKoqKrlkyyHiciO9A7IGJM5HDx4kBYtWvDzzz/zySef0KVLF1+HZHzMk4SxRUTqqOq3ACJyL64igsaYLObXX38lMDCQw4cPs2zZMoKCgnwdkskEPEkYNXAVINznLN8G/CIiPwGqqpXTPTpjTIbbvHkzLVq0ID4+nrVr11K7tg1iNC6eJIw0lzY3xmRuq1evpl27dtx4442EhoZSrlw5X4dkMpE015JS1b3OsjHGj82bN4+WLVsSEBBAZGSkJQtzhYysJWWMyaTef/99OnfuzL333kt4eDi33pps0WiTzWVYLSljTOajqrzyyiv079+fNm3aEBoaSpEiRXwdlsmkPOnD8FYtKWOMD1y8eJG+ffsybdo0evbsyeTJk7nuOk8+Ekx240kLI6GWVDGnltQG4HWvRGWM8ap//vmHjh07Mm3aNIYNG8YHH3xgycJcU2prSQlWS8oYv3TixAlCQkLYsGED48ePp3///r4OyfgJj75SqOouYJeXYjHGeNlff/1FcHAwu3btYs6cOXTu3NnXIRk/cs2EISLPXe15Vf1f+oVjjPGW3bt3ExgYyNGjR/nyyy9p1qyZr0MyfsadFkYB52c5oBaw1FluA3zvjaCMMelr48aNtGzZEhEhLCyMmjVr+jok44eumTBUdSSAiIQD1VX1tLP8KrDcq9EZY9IsNDSU9u3bU7RoUUJDQylbtqyvQzJ+ypNRUsWBC0mWLziPGWMyqTlz5tCqVSvuuOMOIiMjLVmYNPF0PozvRWSRs9wW+Ci9AzLGpI9x48YxcOBA7r//fpYsWUKhQoV8HZLxc263MFR1NNADOO7ceqjqG+6uLyJ5ROR7EflRRH4WkYRTXWVE5DsR+VVEPhMRu3rcmDRQVV588UUGDhxI+/btWblypSULky48HVa7BdiSyn2dB5qo6hkRyQVsEJEVwHPAWFWdKyKTgZ7ApFTuw5hs7eLFi/Tp04cZM2bQp08fJk6cSM6cOX0dlskiPOnDSBN1OeMs5nJuCjQBFjiPz8R1qssY46Fz587Rrl07ZsyYwfDhw5k8ebIlC5OuMrQWgFOLajNwJ65Kt78BJ1T1ovOS/YCVyTTGQ8eOHSMkJITIyEgmTJjAU0895euQTBaUoRfuqWocUFVECuOqS3W3u+uKSB+gD8Btt93m7mrGZAmqiogku7x//36Cg4PZs2cP8+bN48EHH/RVmCaL88mFe6p6QkTCcM2vUVhErnNaGSWBAymsMxWYClCzZk1NzX6N8Udjv9rNqZhYhrcuj4igqoxatoOCeXIRXDKOoKAgTpw4wYoVK2jSpImvwzVZ2DX7MFR1pHPxXklcF+49r6rP45rj2+2v+iJS1GlZICJ5gebATiAMSPhK1A1Y4tE7MCYLU1VOxcQyIyKKUct2JCaLGRFR7PhxM/Xr1+fChQt8/fXXliyM13nSh5HWC/duAWY6/Rg5gHmqukxEdgBzReQ14AfgQw+2aUyWJiIMb10egBkRUcyIiAKgfp4/mTPieUqUKEFoaCh33HGHD6M02UVaL9yb6e7KqroNqJbM478DtT2Iw5hsRUQYfmorPSc9y39OHeG9PPkZFHuOSpUrs2LFCooXt4ILJmN4Mh/GaOe6iQbOQz1U9QfvhGWMSaCzZxPbsxclz8fwDjAo5jSNJAeL+/alkCULk4Hcvg5DXEMyygOFVHUccFRErGVgjBepKiefHcx152MYDAzC1eG3UuPRYSNQtfEfJuN4cuHeRFyjmh52lk/jupbCGOMlIkK+6IP0AN4GngLmAtcDhY4cumSorTHe5knCuFdVnwZiAFT1OGB1n4zxorNnz9Iub15mAaOA94GEa7fFrkcyGcyTTu9YZ4STgmuYLBDvlaiMMRw9epTWrVvzfUwMU3Lnps+FJIMU8+WD0aN9F5zJljxpYYzHdXV2MREZDWwAXvdKVMZkc/v27aN+/fr88MMPzF+wgD7Tp0Pp0iDi+jl1KnTp4uswTTbjySip2SKyGWgKCNBWVXd6LTJjsqmff/6Z4OBgTp06xapVq7j//vtdT1iCMD7mdsJwakp9pqrW0W2Ml0RGRtK6dWuuv/56wsPDqVKliq9DMiaRJ6ekCgChIrJeRPqJiA0ANyYdLVu2jGbNmnHTTTcRGRlpycJkOp7MuDdSVSsAT+Mq8/G1iKz2WmTGZCMfffQRbdu2pXz58kRERFCmTBlfh2TMFVIzgdJh4BBwFCiWvuEYk72oKm+99RY9evSgcePGhIWFUayY/VuZzMmTK72fEpF1wBrgJqC3qlb2VmDGZHXx8fE8//zzDB06lM6dO7N8+XIKFChw7RWN8RG3Or2dsiA1gIGqutWrERmTDVy4cIEePXowZ84c+vfvz7vvvkuOHBk2Y7IxqeLWX6i6Ctbca8nCmLQ7c+YMISEhzJkzh9dff51x48ZZsjB+wZMrvTeLSC1V3ei1aIzJ4qKjo2nVqhWbN29m2rRp9OzZ09chGeM2TxLGvUAXEdkLnMV18Z5aP4Yx7tm7dy+BgYHs27ePRYsWERIS4uuQjPGIJwkjyGtRGJPF/fTTTwQHB3Pu3DlCQ0Np0KDBtVcyJpPxpDTIXm8GYkxWtX79etq0acMNN9xAeHg4lSpV8nVIxqSKRxMoiUhXERnuLN9mEygZc3VLliwhMDCQEiVKEBkZacnC+LUMm0BJREqJSJiI7BCRn0VkgPP4jSLylYjscX4W8SAmYzKtadOm0b59eypXrsyGDRsoXbq0r0MyJk0ycgKli8DzqloeqAM8LSLlgaHAGlUti+uiwKEebNOYTEdVef311+nduzfNmzdnzZo13Hzzzb4Oy5g08yRhpGkCJVU9qKpbnPungZ3ArcADwEznZTOBth7EZEymEh8fz4ABAxg2bBhdunRh6dKl5M+f39dhGZMufDKBkogEANWA74DiqnrQeeoQYFVwjV86f/48jzzyCO+99x7PPvsss2bNIndum8XYZB2pnUAJ4AFV3eXpDkUkP/A5rjIjp5JOYq+qKiKawnp9gD4At9lcxiaTOX36NO3bt2f16tW89dZbDB48mKR/28ZkBZ6MkuoIHHAmULoReF1EqnuyMxHJhStZzFbVhc7Df4vILc7zt+CqhnsFVZ2qqjVVtWbRokU92a0xXnX48OHESrMzZszghRdesGRhsiRPTkm9oqqnRaQ+0AT4EJjk7spOAcMPgZ2q+r8kTy0Fujn3uwFLPIjJGJ/6448/qFevHjt27GDx4sV0797d1yEZ4zWeJIw452cr4ANVXY5no6TqAY8CTURkq3NrCbwJNBeRPUAzZ9mYTO/HH3+kbt26HD16lNWrV9O6dWtfh2SMV3lSGuSAiEwBAoG3ROR6PJuxbwOu+lPJaZrC48ZkSl9//TUhISEULFiQNWvWUL58eV+HZIzXedLC6ASsApqr6glc/RiDvRGUMZnZwoULCQoK4tZbbyUyMtKShck2PEkY8UAZ4P9E5HOgMxDulaiMyaSmTJlCx44dqVatGuvXr6dUqVK+DsmYDONJwpgF3AO8B7wPlAc+9kZQxmQ2qsqoUaN48sknCQ4OZvXq1dx0002+DsuYDOVJH0ZFp6xHgjAR2ZHeARmT2cTFxfHMM88wceJEHnvsMaZNm0auXLl8HZYxGc6TFsYWEamTsCAi9wKb0j8kYzKP8+fP07lzZyZOnMjgwYP56KOPLFmYbOuaLQwR+QlX/ahcQKSI7HOeug3w+EpvY/zFqVOnaNu2LWFhYbz99ts8//zzvg7JGJ9y55SUDS432c6hQ4do0aIF27dv5+OPP6Zr166+DskYn7tmwkg6056IVAES5pZcr6o/eiswY3zlt99+IzAwkEOHDrF06VJatGjh65CMyRQ8qSU1AJgNFHNun4hIf28FZowv/PDDD9StW5eTJ0+ydu1aSxbGJOHJKKmeuCZROgsgIm8B3+AaZmuM31u7di1t27alcOHChIaGcvfdd/s6JGMyFU9GSQn/1pPCuW8lOU2WMH/+fFq0aMFtt91GZGSkJQtjkuFJC2MG8J2ILHKW2+KqPmuMX5s4cSL9+vWjbt26fPHFFxQpYtPKG5McT4oH/g/oARxzbj1U9V0vxWWM16kqw4cP5+mnn6Z169aEhoZasjDmKty5DkNUVQGcObm3XO01xviDuLg4nnrqKaZOncrjjz/OlClTuO46TxrcxmQ/7rQwwkSkv4hcMi+qiOQWkSYiMpN/J0AyJtOLiYmhY8eOTJ06lRdffJFp06ZZsjDGDe78lwQDjwOfikgZ4ASQF1eyCQXeVdUfvBahMenoxIkTPPDAA4SHh/Puu+8yYMAAX4dkjN9w58K9GGAiMNGZk/tm4B9nTgxj/MbBgwcJDg5m586dzJkzh4cfftjXIRnjVzxqh6tqLHDQS7EY4zV79uwhMDCQ6Oholi1bRmBgoK9DMsbv2Ilbk+Vt2rSJli1boqqEhYVRq1YtX4dkjF/y5MI9Y/zOV199RePGjcmXLx8RERGWLIxJA09qSX3lFB9MFRGZLiKHRWR7ksdudLa7x/lpg+BNupk7dy6tWrWiTJkyREZGctddd/k6JGP8mictjCHAuyIyQ0RuScW+PsI14iqpocAaVS0LrHGWjUmz8ePH8/DDD1OnTh3Cw8P5z3/+4+uQjPF7nlzpvUVVGwPLgJUiMkJE8nqwfjiuK8STegCY6dyfiavciDGppqoMGzaMAQMG0LZtW1atWkXhwoV9HZYxWYJHfRgiIsAvwCSgP7BHRB5Nw/6Lq2rCqKtDQPGr7LuPiGwSkU3R0dFp2KXJqi5evEjv3r15/fXX6d27N/PnzydvXre/0xhjrsGTPowI4AAwFrgV6A40AmqLyNS0BuKUFkmxvIiqTlXVmqpas2jRomndncli/vnnHzp06MCHH37IK6+8YqU+jPECT/6j+gA7kqkZ1V9EdqZy/3+LyC2qetDpFzmcyu2YbOz48eOEhIQQERHBe++9R79+/XwdkjFZkid9GD9fpcBgq1Tufyn/1qHqBixJ5XZMNnXgwAEaNmzId999x9y5cy1ZGONF6dJmV9Xfr/UaEfkU1ymsm0VkPzACeBOYJyI9gb1Ap/SIx2QPu3btIigoiGPHjrFixQqaNm3q65CMydIy7CSvqqZUuMf+y43Hvv/+e1q2bEnOnDn5+uuvqV69uq9DMibLsyu9jd9ZuXIljRs3plChQkRERFiyMCaDWMIwfmX27Nm0adOGsmXLEhERwZ133unrkIzJNixhGL8xduxYunbtSv369fn6668pUaKEr0MyJluxhGEyPVVlyJAhPPfcc3To0IEVK1ZQqFAhX4dlTLZjVzaZTC02NpbevXszc+ZMnnzySd5//31y5szp67CMyZashWEyrXPnztGuXTtmzpzJq6++ysSJEy1ZGOND1sIwmdKxY8do3bo13377LZMmTeLJJ5/0dUjGZHuWMEym8+effxIUFMRvv/3G/Pnz6dChg69DMsZgCcNkMjt37iQwMJBTp06xatUqGjVq5OuQjDEO68MwmcY333xD/fr1iY2N5euvv7ZkYUwmYwnDZArLly+nadOmFClShMjISKpWrerrkIwxl7GEYXxu1qxZPPDAA9xzzz1ERERw++23+zokY0wyLGEYnxozZgzdunWjUaNGhIWFUbx4ipMuGmN8zBKG8Yn4+HgGDRrECy+8QKdOnVi+fDkFCxb0dVjGmKuwUVImw8XGxvL444/zySef0K9fP8aNG0eOHPbdxZjMzhKGyVBnz57lwQcfZOXKlbz22mu89NJLiIivwzLGuMEShskwR44coVWrVmzatImpU6fSu3dvX4dkjPGAJQyTIfbu3UtQUBBRUVF8/vnntG3b1tchGWM8ZAnDeN327dsJDg7mzJkzhIaG0rBhQ1+HZIxJhUzR0ygiwSLyi4j8KiJDfR2PST8bNmygQYMGxMfHEx4ebsnCGD/m84QhIjmBCUALoDzwsIiU921UJj188cUXNG/enGLFihEZGUnlypV9HZIxJg18njCA2sCvqvq7ql4A5gIP+Dgmk0bTp0+nXbt2VKxYkQ0bNhAQEODrkIwxaZQZ+jBuBf5MsrwfuPfyF4lIH6CPs3heRLZnQGzecjNwxNdBpJJHsW/atIlixYp5MRyP+fOxB4vf1/w9/nJpWTkzJAy3qOpUYCqAiGxS1Zo+DinV/Dl+f44dLH5fs/h9S0Q2pWX9zHBK6gBQKslySecxY4wxmUhmSBgbgbIiUkZEcgOdgaU+jskYY8xlfH5KSlUvikg/YBWQE5iuqj9fY7Wp3o/Mq/w5fn+OHSx+X7P4fStN8YuqplcgxhhjsrDMcErKGGOMH7CEYYwxxi1+lTD8rYSIiJQSkTAR2SEiP4vIAOfxG0XkKxHZ4/ws4utYr0ZEcorIDyKyzFkuIyLfOb+Hz5zBCpmSiBQWkQUisktEdorIff50/EXkWedvZ7uIfCoieTLz8ReR6SJyOOl1Uikdb3EZ77yPbSJS3XeRpxj7GOdvZ5uILBKRwkmee9GJ/RcRCfJJ0EkkF3+S554XERWRm53lVB17v0kYflpC5CLwvKqWB+oATzsxDwXWqGpZYI2znJkNAHYmWX4LGKuqdwLHgZ4+ico944CVqno3UAXX+/CL4y8itwLPADVVtSKuQSGdydzH/yMg+LLHUjreLYCyzq0PMCmDYkzJR1wZ+1dARVWtDOwGXgRw/o87AxWcdSY6n1G+9BFXxo+IlAICgX1JHk7VsfebhIEflhBR1YOqusW5fxrXh9WtuOKe6bxsJtDWJwG6QURKAq2Aac6yAE2ABc5LMm38IlIIaAh8CKCqF1T1BH50/HGNZMwrItcB+YCDZOLjr6rhwLHLHk7peD8AzFKXb4HCInJLhgSajORiV9VQVb3oLH6L6zoxcMU+V1XPq+ofwK+4PqN8JoVjDzAWeAFIOsIpVcfenxJGciVEbvVRLB4TkQCgGvAdUFxVDzpPHQKK+youN7yL648t3lm+CTiR5J8oM/8eygDRwAznlNo0EbkBPzn+qnoAeBvXN8ODwElgM/5z/BOkdLz97X/6cWCFc98vYheRB4ADqvrjZU+lKn5/Shh+S0TyA58DA1X1VNLn1DWuOVOObRaR1sBhVd3s61hS6TqgOjBJVasBZ7ns9FMmP/5FcH0TLAP8B7iBZE45+JPMfLyvRkSG4TrFPNvXsbhLRPIBLwHD02ub/pQw/LKEiIjkwpUsZqvqQufhvxOaf87Pw76K7xrqASEiEoXrFGATXH0ChZ1TJJC5fw/7gf2q+p2zvABXAvGX498M+ENVo1U1FliI63fiL8c/QUrH2y/+p0WkO9Aa6KL/XrjmD7HfgevLxo/O/3BJYIuIlCCV8ftTwvC7EiLO+f4PgZ2q+r8kTy0Fujn3uwFLMjo2d6jqi6paUlUDcB3vtaraBQgDHnRelpnjPwT8KSIJFTqbAjvwk+OP61RUHRHJ5/wtJcTvF8c/iZSO91LgMWfETh3gZJJTV5mCiATjOiUboqrnkjy1FOgsIteLSBlcncff+yLGlKjqT6paTFUDnP/h/UB15/8idcdeVf3mBrTENVLhN2CYr+NxI976uJrf24Ctzq0lrn6ANcAeYDVwo69jdeO9NAKWOfdvx/XP8SswH7je1/FdJe6qwCbnd7AYKOJPxx8YCewCtgMfA9dn5uMPfIqrvyXW+YDqmdLxBgTXyMffgJ9wjQbLbLH/iutcf8L/7+Qkrx/mxP4L0CIzHvvLno8Cbk7LsbfSIMYYY9ziT6ekjDHG+JAlDGOMMW6xhGGMMcYtljCMMca4xRKGMcYYt1jCMMYY4xZLGMYYY9xiCcP4nIjkFZGvE8pDi0ikh+u/KiKDvBOd9/Ytrrk6nkrvmJJs/6rHUUQCkps7IcnzZ67yXG4RCU9SosRkA5YwTGbwOLBQVeMAVLWuj+PJKIWBdE8YTrmHHN48juqaYmAN8JC39mEyH0sYxqvENeNgc+f+ayLyXjIv60KSekgJ32ydb8A7ReQDcc06FyoieZ3nhonIbhHZAJRLsm5XEfleRLaKyBRxzRYY4MyaNtvZ3gKnkue11vFo36nY1pvAHc5rx1y2nTdF5Okky4ktGRFZLCKbnW31SXKsfhGRWbjKiJRK2kJIbh3HdSkdl6u9J+epxc7vzmQXvq5/YresfcM1gdE6XB8sy4Gclz2fGzh02WNnnJ8BuEpKV3WW5wFdgRq46t/kAwriqvczCLgH+ALI5bx+IvCYsx0F6jmPTwcGOfevto7b+07ltgKA7Skct2rA10mWdwClnPsJtZjy4koONznbigfqXH4cr7FOssclYf2U3pNzPycQ7eu/Mbtl3M3OPxqvUtVwp9Lqc0AjdU47JXEzcOIqm/hDVbc69zfj+pC7GVikTvVQEUmoWtwU1wf6RtcuyYurlHY48KeqRjiv+wTX1KdvX2MdT/Z9rf0nt60NKb1pVf1BRIqJyH+AosBxVU2Y8OYZEWnn3C+Fq1LqIWCvumZPS05K66R0XK71nlDVOBG5ICIF1DWjpMniLGEYrxKRSsAtwNEUPlT+AfJcZRPnk9yPw/WBleLugJmq+uJlMQRw5aQ96sY6nuw7vbcFrkq0DwIlgM+cbTXCNU/Gfap6TkTW8e/xO5tsUFdfJ6XjctX3lMT1QIwb78VkAdaHYbxGXJPlzMY1a9wZcc0tcAlVPQ7kFJGrJY3LhQNtndFVBYA2zuNrgAdFpJiz/xtFpLTz3G0icp9z/xH+/XZ/tXU82XdqtnUaKHCV5z/DNQ/Jg7iSB0AhXK2NcyJyN1DnKusnuNo6KR2Xa74nEbkJOKKuyZ1MNmAJw3iF03m6EHheVXcC/wVGpPDyUFxzh7hFVbfg+jD9Edccyxudx3cALwOhIrIN+ApX6wZccxY8LSI7cc2JMcmNddzedyq3dRSIEJHtl3d6O8//jCuhHNB/J7dZiaujeieuTvOUTkEldbV1kj0ubr6nxrj6pUw2YfNhGJ8TkerAs6r6qJe2H4Br8qeK3th+diUiC4Ghqrrb17GYjGEtDONzzrf2sCTDNU0mJ65pkhdbssherIVhjDHGLdbCMMYY4xZLGMYYY9xiCcMYY4xbLGEYY4xxiyUMY4wxbrGEYYwxxi2WMIwxxrjl/wGn93LM5A0n8AAAAABJRU5ErkJggg==\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"plt.scatter(x, y, marker=\"x\", label=\"Observed values\")\n",
|
|
"plt.scatter(x, predictions, marker=\"o\", color=\"red\", label=\"Predicted values\")\n",
|
|
"reg_x_vals = np.linspace(0, max(x)+10, 10)\n",
|
|
"reg_y_vals = np.array([(theta_1 * x_i) + theta_0 for x_i in reg_x_vals])\n",
|
|
"plt.plot(reg_x_vals, reg_y_vals, color=\"black\", label=\"Regression line\")\n",
|
|
"plt.xlim([0, max(x)+10])\n",
|
|
"plt.ylim([0, max(y)+10])\n",
|
|
"plt.xlabel(\"$x$ (independent variable)\")\n",
|
|
"plt.ylabel(\"$y$ (observed dependent variable)\")\n",
|
|
"plt.legend()\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Now let's quantify how accurately our model fits the training dataset. A commonly used metric for regression models is root-mean-square-error (RMSE), which may be calculated as:\n",
|
|
"\n",
|
|
"\\begin{equation}\n",
|
|
" \\label{eqn:rmse}\n",
|
|
" RMSE = \\sqrt{\\frac{1}{m}\\Sigma_{i=1}^{m}{\\Big(\\hat{y}_i - \\vec{y}_i\\Big)^2}}\n",
|
|
"\\end{equation}\n",
|
|
"\n",
|
|
"where $\\hat{y}$ is the vector of predictions from our model for each $x$ value in our training dataset, $\\vec{y}$ is a vector containing the actual observed values for each x value in our training dataset.\n",
|
|
"\n",
|
|
"The coefficent of determination $R^2$ is another commonly used metric for regression models.\n",
|
|
"\n",
|
|
"N.B. for the sake of simplicity, the RMSE and $R^2$ metrics are calculated using training data only. To evaluate the model properly we should also claculate these metrics separately on on independent test data (more on this in later lectures!) "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 38,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"RMSE: 1.8894703065186529\n",
|
|
"Coefficient of determination: 0.9815885948385498\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# this function calculates RMSE as per the equation above\n",
|
|
"def rmse(predictions, targets):\n",
|
|
" return np.sqrt(((predictions - targets) ** 2).mean())\n",
|
|
"\n",
|
|
"# calculate and print the RMSE\n",
|
|
"calculated_rmse1 = rmse(predictions, y)\n",
|
|
"print(\"RMSE:\", calculated_rmse1)\n",
|
|
"\n",
|
|
"# calculate and print the coefficient of determination\n",
|
|
"r_sq = model.score(x, y)\n",
|
|
"print(\"Coefficient of determination:\", r_sq)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Developing a linear regression model using numpy and the singular value decomposition\n",
|
|
"This section will demonstrate an alternative method to develop a least squares linear regression model using the Singular Value Decomposition (SVD).\n",
|
|
"\n",
|
|
"Matrix decomposition techniques allow a given matrix to be factorised as a product of matrices. Eigendecomposition is one of the most widely used matrix decomposition methods, which allows a matrix to be decomposed into a set of eigenvectors and eigenvalues. The SVD may be used to factorise a real or complex matrix into singular values and singular vectors. All real matrices have a SVD, which makes it more generally applicable than other matrix factorisation methods (e.g. eigendecompositions are not defined for $m \\times n$ matrices, i.e. matrices which are not square). The SVD may be used to factorise a matrix $A$ as follows:\n",
|
|
"\n",
|
|
"\\begin{equation}\n",
|
|
"\\label{eqn:svd}\n",
|
|
"A = U\\Sigma V^{T}\n",
|
|
"\\end{equation}\n",
|
|
"\n",
|
|
"where: \n",
|
|
"* $A$ is a $m \\times n$ matrix.\n",
|
|
"* $U$ is a $m \\times m$ orthogonal matrix (i.e. a matrix with rows and columns comprised of orthogonal unit vectors) containing the **left singular vectors**. The left singular vectors are a set of orthonormal (i.e. both orthogonal and normal) eigenvectors of $AA^{T}$.\n",
|
|
"* $\\Sigma$ is a $m \\times n$ diagonal matrix, where each of the non-negative diagonal entries $\\{\\sigma_{1},\\sigma_{2},...\\}$ is a **singular value**. The singular values are are the square roots of the eigenvalues of $A^{T}A$ and of $AA^{T}$. Note that there are $\\mathrm{min}(m,n)$ singular values.\n",
|
|
"* $V$ is an $n \\times n$ orthogonal matrix containing the **right singular vectors**. The right singular vectors are a set of orthonormal eigenvectors of $A^{T}A$.\n",
|
|
"\n",
|
|
"\\begin{equation}\n",
|
|
"\\label{eqn:svd_k}\n",
|
|
"A_k = U_k \\Sigma_k V^{T}_k\n",
|
|
"\\end{equation}\n",
|
|
"\n",
|
|
"Taking a linear algebra perspective, a linear regression problem may be modelled as follows (Eqn. 3):\n",
|
|
"\n",
|
|
"\\begin{equation}\n",
|
|
" X \\vec{\\theta} = \\vec{y}\n",
|
|
"\\end{equation}\n",
|
|
"\n",
|
|
"where $X$ is a matrix containing the measured values for the independent variables $\\{x_{1},x_{2},...,x_{j}\\}$, $\\vec{\\theta}$ is a column vector containing the set of model weights $\\{\\theta_{0},\\theta_{1},\\theta_{2},...,\\theta_{j}\\}$ and $\\vec{y}$ is a column vector containing the observations of the dependent variable."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Note that when constructing $X$, to account for constant observed effects on the value of the target variable $y$, an all-ones column vector should be appended to the left of the columns containing the observations of the dependent variables, so that the format of $X$ is as follows (Eqn. 4):\n",
|
|
"\n",
|
|
"\\begin{equation}\n",
|
|
"X=\n",
|
|
" \\begin{bmatrix}\n",
|
|
" 1 & x_{1,1} & x_{1,2} & ... & x_{1,j} \\\\\n",
|
|
" 1 & x_{2,1} & x_{2,2} & ... & x_{2,j} \\\\\n",
|
|
" ... & ... & ... & ... & ... \\\\\n",
|
|
" 1 & x_{m,1} & x_{m,2} & ... & x_{m,j} \n",
|
|
" \\end{bmatrix}\n",
|
|
"\\end{equation}\n",
|
|
"\n",
|
|
"where $m$ is the number of data points in the training set, and the notation $x_{m,j}$ refers to the observation in the $m$th data point of the $j$th independent variable."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"If $X$ is invertible, the set of model weights could easily computed by premultiplying each side of Eqn. 3 by $X^{-1}$ to obtain the following result (Eqn. 5):\n",
|
|
"\n",
|
|
"\\begin{align}\n",
|
|
" \\label{eqn:linregmatrixnaive}\n",
|
|
" X^{-1} X \\vec{\\theta} = X^{-1} \\vec{y} \\nonumber\\\\\n",
|
|
" \\vec{\\theta} = X^{-1} \\vec{y}\n",
|
|
"\\end{align}"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"However, a matrix is only invertible if it is square (i.e. $m=n$) and if all of its columns are linearly independent. Therefore, while Eqn. 5 could be used to calculate $\\vec{\\theta}$ exactly for a given $X$ and $\\vec{y}$, assuming that $X$ is invertible greatly restricts its applicability, requiring all observed values of the target variable to fit the model exactly. If $X$ has more rows than columns (i.e. $m$ > $n$) it is possible for Eqn. 5 to have no solution, whereas if $X$ has more columns than rows (i.e. $m$ < $n$) it is possible for Eqn. 5 to have multiple different solutions."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"A more general result to compute an approximation of $\\vec{\\theta}$ which minimises the Euclidian norm $|| X \\vec{\\theta} - {\\vec{y}} ||_2$ may be developed using the Moore-Penrose psuedoinverse. The psuedoinverse $A^{+}$ of a given matrix $A$ may be calculated as (Eqn. 6):\n",
|
|
"\n",
|
|
"\\begin{equation}\n",
|
|
" \\label{eqn:psuedoinverse}\n",
|
|
" A^{+} = V \\Sigma^{+} U^{T}\n",
|
|
"\\end{equation}\n",
|
|
"\n",
|
|
"where $U$, $\\Sigma$ and $V$ are the SVD of A. $\\Sigma^{+}$ is obtained by first taking the reciprocal of the nonzero entries of the diagonal matrix $\\Sigma$ containing the singular values and then taking the transpose of the resulting matrix."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"To compute the approximation of the weights $\\vec{\\theta}$, $X^{+}$ may be used as a substitute for $X^{-1}$ in Eqn. 5 to obtain the following result (Eqn. 7):\n",
|
|
"\n",
|
|
"\\begin{align}\n",
|
|
" \\label{eqn:linregmatrixpsuedo}\n",
|
|
" \\vec{\\theta} = X^{+} \\vec{y}\n",
|
|
"\\end{align}"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 39,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"X:\n",
|
|
" [[ 1. 80.]\n",
|
|
" [ 1. 110.]\n",
|
|
" [ 1. 110.]\n",
|
|
" [ 1. 130.]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# the dimensions of the matrix X\n",
|
|
"m = len(data)\n",
|
|
"n = len(data[0]) # add an additional column for the all-ones vector\n",
|
|
"\n",
|
|
"# Create and print X - the matrix of observed values for the independent variables\n",
|
|
"X = np.c_[np.ones(m), data[:,0]]\n",
|
|
"print(\"X:\\n\", X)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 40,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Psuedoinverse of X:\n",
|
|
" [[ 2.56862745e+00 3.92156863e-02 3.92156863e-02 -1.64705882e+00]\n",
|
|
" [-2.15686275e-02 1.96078431e-03 1.96078431e-03 1.76470588e-02]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# calculate and print the Moore-Penrose psuedoinverse of X\n",
|
|
"X_pi = np.linalg.pinv(X)\n",
|
|
"print(\"Psuedoinverse of X:\\n\", X_pi)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 41,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Vector of caclulated model weights:\n",
|
|
" [[-34.44509804]\n",
|
|
" [ 0.7727451 ]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# calculate and print the values in the model weights vector theta using Eqn. 7\n",
|
|
"theta = np.dot(X_pi, y)\n",
|
|
"print(\"Vector of caclulated model weights:\\n\",theta)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 42,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[[27.3745098 ]\n",
|
|
" [50.55686275]\n",
|
|
" [50.55686275]\n",
|
|
" [66.01176471]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# calculate and print the predicted y values when using the values from X\n",
|
|
"y_hat = np.dot(X, theta)\n",
|
|
"print(y_hat)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Finally, let's quantify how accurately our model fits the training dataset as before."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 43,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"1.889470306518652\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# calculate and print the RMSE\n",
|
|
"calculated_rmse2 = rmse(y_hat, y)\n",
|
|
"print(calculated_rmse2)"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "ct4101",
|
|
"language": "python",
|
|
"name": "ct4101"
|
|
},
|
|
"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.10.7"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|