In machine learning, we learn(fit) a model to some data using only a small sample of the data called the training set. The model is then used to make predictions about similar but unseen data. Linear regression is a model that represents data using a stright line. Linear regression got its name for historical reasons. In simple language it means line fitting.
But why do we need learn or fit a model, in our case a line? Because a model lets us approximate the data using only a few quantities(parameters) and draw conclusions that might not be evident looking at the bulk of the data.
In the rest of the post I will show you why linear regression counts as machine learning, when it is OK to use linear regression and how it actually works.
Most typical introductions to linear regression use matrix algebra. In this post, however, I will use simple algebra and calculus to derive all equations. Interestingly, Python has a very powerful symbolic math library called SymPy that will let us do all calculation like we do them with pen and paper.
The data we will use here comes from CMU's Body Fat dataset. This data is stored in a csv file. And in order to work with this data I need to access it inside the notebook. For this purpose I am using the Pandas library.
First I import the Pandas library and then use the read_csv() method to load our data. After that I print the first few rows from the data.
import pandas as pn
# Read data and print some rows
dataframe = pn.read_csv('bodyfat.csv', sep="\s+")
dataframe.head()
As we can see from the data for the first 5 persons, there are a lot of columns. It'd be interesting to see how the column values are related to each other. A good approach is to use a scatter plot matrix to visualize dependencies between different columns. Data visualization is considered a crucial step while dealing with machine learning algorithms. The code below extracts a subset of the data, columns 3 to 8 and visualizes the relationship between every pair of columns from this subset using the Plotly library.
import numpy as np
df = pn.DataFrame(dataframe, columns=list(dataframe.columns.values)[3:8])
import plotly as py
import plotly.figure_factory as ff
import plotly.graph_objs as go
py.offline.init_notebook_mode()
fig = ff.create_scatterplotmatrix(df,diag='histogram',height=800, width=800)
py.offline.iplot(fig)
The plots falling on the main diagonal are actually histograms not scatter plots as it doesn't make sense to plot a column against itself. The histogram are interesting in their own self but we want to focus on relationships between pairs of columns which are all off the diagonal. There are some relationships that are almost linear. To me "Weight vs Abdomen" relationship appears to be a good choice. So let's focus on this relationship and visualize it inside a single zoomed in plot.
This relationship is ideal for fiting a line because it is 2-dimensional and all the data points seem to lie along a line. If the data points were not distributed like this trying to fit a line would not have made much sense. The goal now is to find such a line which is close enough to most of the points in the data set.
Before findinf the best line we need to 'hide' some of our data points. Once the best line is found we will use these hidden points to measure the accuracy of our predictions.
I'm going to split data set into two parts. The first part consists of 60% data and second one 40%. In ML domain these parts have special names: training data, and testing data. The training data is used to find the line while the test data is later used to find if this line provides good approximation for new unseen values.
One important thing to notice while spliting data is that both train and test sets should have almost same distribution. To ensure this, we split the data randomly.
# Split data between train and test set
nTrain = int(0.6 * weight_abdomen_data.shape[0])
train = weight_abdomen_data.iloc[0:nTrain]
test = weight_abdomen_data.iloc[nTrain :]
abdomen_train = train.iloc[:,1]
weight_train = train.iloc[:,0]
abdomen_test = test.iloc[:,1]
weight_test = test.iloc[:,0]
In the plot below we can see the how the data in two sets look like.
Now the question is how do we find a line that approximates the (training) data well. In the figure below the orange line passing through the middle of the data appears to be a good approximation. However, the green line which also passes through the center of the data does not appear to be such a good approxiation.
x_min = abdomen_train.min()
x_max = abdomen_train.max()
y_min = 120
y_max = 320
green_line = go.Scatter(x=[x_min,x_max], y=[y_min,y_max], mode='lines', name='line 1')
red_line = go.Scatter(x=[x_min,x_max], y=[80,400], mode='lines', name='line 2')
scatter_layout = go.Layout(
title='Candidate Lines',
autosize=True,
width=500,
height=500,
xaxis=dict(title='Abdomen (cm)'),
yaxis=dict(title='Weight (lbs)'))
fig = go.Figure(data=[trace_scatter_train,green_line, red_line], layout= scatter_layout)
py.offline.iplot(fig)
We need an objective measure of how good a line. Since the whole point of fitting a line is to predict weight values from abdomen sizes, one way to measure the goodness-of-fit of a line is to measure the difference between the actual and predicted weights for every abdomen size in the training set, which is the sum of lengths of the vertical red lines in the plot below.
green_slope = (y_max-y_min)/(x_max-x_min)
green_intercept = -x_min*green_slope+y_min
scatter_layout = go.Layout(
title='Prediction Errors',
autosize=True,
width=1000,
height=500,
xaxis=dict(title='Abdomen (cm)'),
yaxis=dict(title='Weight (lbs)'),
showlegend=False)
traces = [trace_scatter_train,green_line]
for x,y in zip(abdomen_train,weight_train):
y_pred = green_slope * x + green_intercept
traces.append(go.Scatter(x=[x,x],y=[y,y_pred],mode='lines',line = dict(color = ('rgb(255, 0, 0)'))))
fig = go.Figure(data=traces, layout= scatter_layout)
py.offline.iplot(fig)
If for a given abdomen mesurement $x$, the actual weight is $y$ and the weight predicted using our line is $\hat y$ the error of prediction is $\lvert y-\hat y \rvert$. Taking the absolute value ensures that errors with opposing signs add up instead of cancelling out. The total error is $\sum_1^n {\lvert y-\hat y \rvert}$ if there are $n$ (abdomen,weight) pairs in the data. The average prediction error is $\sum_1^n {\lvert y-\hat y \rvert} \over n$. This mean abolute error is computed below using sympy but first an introduction to Sympy is in order.
As promised earlier I will use algebra and calculus to find the desired line. And in order to achieve this I'm going to use a symbolic mathematics library in Python called SymPy. In Sympy we can define variables, use operators and functions to combine the variables in to expressions and perform mathematical operations on these expression to get new expression just as we do with pen and paper.
SymPy.symbols() takes a string of variable names separated by spaces or commas, and creates Symbols out of them. For example, to define a quadratic $y=a*x^2+b*x+c$, we use the code below.
import sympy as sp
# Define the RHS variables as SymPy symbols.
a,b,c,x = sp.symbols('a b c x')
y = a*x*x + b*x + c
If I print the SymPy variable $y$, an expression gets pretty-printed.
sp.pretty_print(y)
Suppose I know the value of $x$ is 80 then I can substitute this in the expression for $y$ as below.
z = y.subs([(x,80)])
print(z)
In our case we have the equation of a line $y=m*x+c$ and a sequence of data in the form of $(x_i,y_i)$, where the subscript $i$ stands for the $i$th data point. To find the absolute cost, for the $i$th data point, we substitute $x_i$ in the equation of our line. This gives us the predicted value $\hat y_i$ for that data point. We then subtract the predicted values from the actual value $y_i$ and take the absolute value of the difference to find the error of prediction.
x_i,y_i,m,c = sp.symbols('x_i,y_i,m,c')
yhat_i = m*x_i+c
abs_error_i = sp.Abs(y_i-yhat_i)
sp.pretty_print(abs_error_i)
The to find the mean absolute error(MAE) we need to sum up the absoluter error over all $(x_i,y_i)$ and then divide by the number of datpoints.
MAE = 0
count = len(abdomen_train)
for abd,wt in zip(abdomen_train,weight_train):
MAE += abs_error_i.subs([(x_i,abd),(y_i,wt)])
MAE = sp.Mul(MAE,1/count,evaluate=False)
sp.pretty_print(MAE,num_columns=100,wrap_line=False)
The mean absolute error is a long, unwieldy expression. It is a function of $m$ and $c$ but because of the absolute values, there is no way to simplify it. It is also not possible to analytically minimize this expression. This brings us to mean squared error (MSE). MSE is different from MAE in that for every $(x_i,y_i)$ the error is calculated as the square of the difference $(y-\hat y)^2$.
sqrd_error_i = (y_i-yhat_i)**2
sp.pretty_print(sqrd_error_i)
The MSE for our data is calculated below. Just like the MAE calculated above, MSE has 151 terms but becasue the individual terms can be expanded and the resulting expression simplified, we are left with a second order polynomial in $m$ and $c$, which has 5 terms.
MSE = 0
count = len(abdomen_train)
for abd,wt in zip(abdomen_train,weight_train):
MSE += sqrd_error_i.subs([(x_i,abd),(y_i,wt)])
MSE = sp.simplify(MSE)
MSE /= count
sp.pretty_print(MSE,num_columns=100,wrap_line=False)
Although it is not evident from the simplified expression, the MSE can only take positive values because in its original form it is a sum of squared terms. Further, the MSE grows larger for large positive or negative values of $m$ or $c$ and is an elliptic paraboloid. Because of these properties this MSE can have at most one minimum that can be found by differentiating the MSE with respect to $m$ and $c$. This could be done in SymPy and yields two linear expressions.
# Find differentiation w.r.t m and c
grad_m = sp.diff(MSE, m)
grad_c = sp.diff(MSE, c)
sp.pretty_print(grad_m,num_columns=100,wrap_line=False)
sp.pretty_print(grad_c,num_columns=100,wrap_line=False)
We can solve the two equations and get the exact values of $m$ and $c$ that minimizes the MSE.
# Solve the two equations
solution = sp.solvers.solve([grad_m, grad_c])
sp.pretty_print(solution, use_unicode=True,num_columns=100)
Geometrically the MSE is the sum of sqaured lengths of the red error lines. The values of $m$ and $c$ we have got minimze the sum of squared lengths.
Let us plot the best line. In the figure below the left plot shows the training data while the right plot shows the cost function, which another name for the error function. The best line can be seen in the left plot the cost corresponding to this line is the position of the tiny sphericial marker in the right surface plot. You can play with the sliders and verify that the cost is indeed minimzed for our solution $(m,c)$.
reset_output()
output_notebook()
show_plots()
show_sliders()
Understanding linear regression requires us to wrap our head around this transformation from the data space $(x,y)$ to parameter space $(m,c)$. The takeaway is that there are an infinite number of lines that can drawn in the data space. Each line has a unique slope $m$ and a unique intercept $c$ and an associated cost or error. Each line is thus a point in the parameter space, where the parameters $m$ and $c$ can be thought of as the x and y coordinates and the cost as the z coordinate.
A brute force method finding the best line would be to try all values of $m$ and $c$ and find the pair of values that minimizes the cost. However $m$ and $c$ are real valued and unbounded and so this brute force method will never terminate.
Fortunately we were able to determine a formula for the cost in terms of $m$ and $c$ and find the best pair $(m,c)$ that minimizes the cost. If we use the means squared error measure as our cost we can find a closed-form formula for finding the MSE (or cost) as a function of $m$ and $c$. The properties of this particular cost function lets us find the cost minimizing values of $m$ and $c$ analytically. We also saw that the mean absolute error MAE can not be minimized in this fashion.