{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"import sys\n",
"import os\n",
"if not any(path.endswith('textbook') for path in sys.path):\n",
" sys.path.append(os.path.abspath('../../..'))\n",
"from textbook_utils import *\n",
"\n",
"csv_file = 'data/Full24hrdataset.csv'\n",
"usecols = ['Date', 'ID', 'region', 'PM25FM', 'PM25cf1']\n",
"full = (pd.read_csv(csv_file, usecols=usecols, parse_dates=['Date'])\n",
" .dropna())\n",
"full.columns = ['date', 'id', 'region', 'pm25aqs', 'pm25pa']\n",
"\n",
"GA = full.loc[(full['id'] == 'GA1') , :]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(sec:linear_simple_fit)=\n",
"# Fitting the Simple Linear Model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To fit a model, we use the notion of loss minimization that was first introduced in\n",
"{numref}`Chapter %s `).\n",
"The first step is to choose the loss function. Then, we find which values of\n",
"$\\theta_0$ and $\\theta_1$ give the smallest loss for our data.\n",
"We call $\\theta_0$ and $\\theta_1$ *model parameters*, and we often represent them as the vector ${\\boldsymbol{\\theta}} = [\\theta_0, \\theta_1]$.\n",
"\n",
"To fit a linear model with a numeric outcome variable, we typically choose squared loss, where for a given data point $(x,y)$ and model parameters ${\\boldsymbol{\\theta}}$, the squared loss is:\n",
"\n",
"$$\n",
"{\\cal l}({\\boldsymbol{\\theta}}, x, y) = [y - (\\theta_0 + \\theta_1 x)]^2\n",
"$$\n",
"\n",
"Notice that this loss function squares the error between the outcome $y$ and the prediction, $\\theta_0 + \\theta_1 x$, given $x$. That's why squared loss is also called *squared error*. We saw in the previous section that these errors are in the vertical direction on the scatter plot, meaning for a specific $x$, the error is the vertical difference between the points $(x, y)$ and $(x, \\theta_0 + \\theta_1 x)$. (See {numref}`Figure %s ` for a diagram that includes two sample errors). "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For a data set with $n$ points: $(x_1, y_1), \\ldots, (x_n, y_n) $, the average squared loss (ASE) is:\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"L({\\boldsymbol{\\theta}}, {\\mathbf{x}}, {\\mathbf{y}})\n",
" &= \\frac{1}{n} \\sum_{i} {\\cal l}({\\boldsymbol{\\theta}}, x_i, y_i) \\\\\n",
" &= \\frac{1}{n} \\sum_{i}[y_i - (\\theta_0 + \\theta_1 x_i)]^2.\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"Again, ${\\mathbf{x}}$ is the vector $[x_1, \\ldots, x_n]$ and ${\\mathbf{y}}$ is similarly defined."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have chosen the loss function, we want to find $\\hat{\\boldsymbol{\\theta}}$ that minimizes the loss. We use calculus to do this."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Minimizing the Loss"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the simple linear model, the average squared error is a function of two model parameters. This means that if we use calculus to find the minimizing parameter values, we need to find the partial derivatives of the ASE with respect to $\\theta_0$ and $\\theta_1$. We can also find these minimizing values by other techniques:\n",
"\n",
"+ *Gradient descent*. We can use numerical optimization techniques, such as gradient descent, when the loss function is more complex and calculus is not an option (see {numref}`Chapter %s `.\n",
"+ *Quadratic formula*. Since $L({\\boldsymbol{\\theta}},{\\mathbf{x}},{\\mathbf{y})}$ is a quadratic function of ${\\boldsymbol{\\theta}}$ we can use the quadratic formula (along with a lot of algebra) to solve for the minimizing model parameter values. We guide you through the quadratic approach in the Exercises. \n",
"+ *Geometric argument*. Later in this chapter, we use a geometric interpretation of least squares to fit multiple linear models. This approach relates to the Pythagorean theorem and has several intuitive benefits."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We choose calculus to optimize the simple linear model since it is quick and straightforward. To begin, we take the partial derivatives of the average squared error with respect to each parameter:\n",
" \n",
"$$\n",
"\\begin{aligned}\n",
"\\frac{\\partial}{\\partial \\theta_0} L({\\boldsymbol{\\theta}}, {\\mathbf{x}}, {\\mathbf{y}})\n",
" &= \\frac{1}{n} \\sum_{i} 2 (y_i - \\theta_0 - \\theta_1 x_i ) (-1)\\\\\n",
" & \\\\ \n",
"\\frac{\\partial}{\\partial \\theta_1} L({\\boldsymbol{\\theta}}, {\\mathbf{x}}, {\\mathbf{y}})\n",
" &= \\frac{1}{n} \\sum_{i} 2 (y_i - \\theta_0 - \\theta_1 x_i) (-x_i) \n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then, we set the partial derivatives equal to 0, and simplify a bit by multiplying both sides of the equations by $-n/2$ to get:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{aligned}\n",
" 0 &= \\sum_{i} (y_i - \\hat{\\theta}_0 - \\hat{\\theta}_1 x_i) \\\\\n",
" 0 &= \\sum_{i} (y_i - \\hat{\\theta}_0 - \\hat{\\theta}_1 x_i)x_i \\\\\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These equations are called the *normal equations*. \n",
"In the first equation, we see that $\\hat{\\theta}_0$ can be represented as a function of $\\hat{\\theta}_1$.\n",
"\n",
"$$\n",
"\\hat{\\theta}_0 = \\bar{y} - \\hat{\\theta}_1 \\bar{x}\n",
"$$\n",
"\n",
"Plugging this value into the second equation, gives us:\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" 0 &= \\sum_{i} (y_i - \\bar y + \\hat{\\theta}_1 \\bar x - \\hat{\\theta}_1 x_i ) x_i \\\\\n",
" &= \\sum_{i} [(y_i - \\bar y) - \\hat{\\theta}_1 ( x_i - \\bar x)]x_i \\\\ \n",
"\\hat{\\theta}_1 &= \\frac{\\sum_{i} (y_i - \\bar y)x_i} {\\sum_{i}( x_i - \\bar x)x_i} \\\\\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After some algebra, we can represent $\\hat{\\theta}_1$ in terms of quantities that we are familiar with:\n",
"\n",
"$$\n",
"\\hat{\\theta}_1 = r({\\mathbf{x}},{\\mathbf{y}}) \\frac{SD({\\mathbf{y}})}{SD({\\mathbf{x}})}\n",
"$$\n",
"\n",
"This representation says that, a point on the fitted line at $x$ can be written as\n",
"\n",
"$$ \n",
"\\begin{aligned}\n",
"\\hat{y} &= \\hat{\\theta}_0 + \\hat{\\theta}_1 x \\\\\n",
"&= \\bar{y} - \\hat{\\theta}_1 \\bar{x} + \\hat{\\theta}_1 x \\\\\n",
"&= \\bar{y} + r({\\mathbf{x}},{\\mathbf{y}}) SD({\\mathbf{y}}) \\frac{(x - \\bar{x})}{SD({\\mathbf{x}})}\\\\ \n",
"\\end{aligned}\n",
"$$\n",
"\n",
"This equation for the line has a nice interpretation: for a given $x$ value, we find how many standard deviations above (or below) average it is, and then predict $y$ to be $r$ times as many standard deviations above (or below) its average. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since `pandas.Series` objects have built-in methods to compute\n",
"$SD(\\mathbf{x})$, $SD(\\mathbf{y})$, and $r(\\mathbf{x}, \\mathbf{y})$,\n",
"we can quickly define functions that fit the simple linear model to our data:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def theta_1(x, y):\n",
" r = x.corr(y)\n",
" return r * y.std() / x.std()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def theta_0(x, y):\n",
" return y.mean() - theta_1(x, y) * x.mean()"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"Now, we can fit the model by computing\n",
"$\\hat{\\theta}_0$ and $\\hat{\\theta}_1$ on the data set."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Model: f(x) = -3.36 + 2.10x\n"
]
}
],
"source": [
"t1 = theta_1(GA['pm25aqs'], GA['pm25pa'])\n",
"t0 = theta_0(GA['pm25aqs'], GA['pm25pa'])\n",
"print(f'Model: f(x) = {t0:.2f} + {t1:.2f}x')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The model we fitted matches the line that we used in the previous section. This is not by accident. We purposely chose the least squares line as our fitted model. Of course, can use the functionality provided in scikit learn for finding the model parameters."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"from sklearn.linear_model import LinearRegression\n",
"\n",
"y = GA[['pm25pa']]\n",
"X = GA[['pm25aqs']]\n",
"\n",
"reg = LinearRegression().fit(X, y)"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Intercept: [-3.36] Slope: [2.1]\n"
]
}
],
"source": [
"print(f\"Intercept: {reg.intercept_} Slope: {reg.coef_[0]}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Rather than computing the model parameters with summary statistics, we will want to use the `LinearRegression` method in scikit learn because it offers numerically stable algorithms. This is especially important when fitting multiple variables. We introduce the multiple linear model in the next section."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}