Linear regression from scratch, then with scikit-learn

The problem

Ten students. Hours studied per day, and the marks they scored.

import numpy as np

hours = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=float)
marks = np.array([35, 42, 48, 51, 60, 63, 70, 74, 80, 86], dtype=float)

We want a line marks = m * hours + b that sits as close to these points as possible.

What "close" means

For each student, the error is the real mark minus the predicted mark. Some errors are positive and some negative, so adding them lets them cancel. Square them first, then take the mean. That is mean squared error:

def mse(y_true, y_pred):
    return np.mean((y_true - y_pred) ** 2)

Squaring also punishes a big miss much more than a small one. An error of 10 contributes 100; an error of 2 contributes 4.

Training means finding the m and b that make this number smallest.

The closed-form solution

For a straight line you do not need to search. Calculus gives the answer directly:

m = sum((x - x_mean) * (y - y_mean)) / sum((x - x_mean) ** 2)
b = y_mean - m * x_mean

In code:

import numpy as np

hours = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=float)
marks = np.array([35, 42, 48, 51, 60, 63, 70, 74, 80, 86], dtype=float)

x_mean = hours.mean()
y_mean = marks.mean()

m = np.sum((hours - x_mean) * (marks - y_mean)) / np.sum((hours - x_mean) ** 2)
b = y_mean - m * x_mean

print("slope     :", round(m, 4))      # 5.5455
print("intercept :", round(b, 4))      # 30.4

pred = m * hours + b
print("MSE       :", round(np.mean((marks - pred) ** 2), 4))

You get a slope of about 5.5455 and an intercept of 30.4.

Read the parameters, do not just print them

• Slope 5.55: in this data, one extra hour of study per day goes with

about 5.5 more marks.

• Intercept 30.4: the line's value at zero hours.

The intercept deserves care. Nobody in this dataset studied zero hours. The smallest value is 1. So 30.4 is the line extended past the data, and the data says nothing about whether it is true. Predicting outside your observed range is called extrapolation, and it is where models embarrass you.

Feed in 40 hours a day and the formula happily returns 252 marks out of 100. The maths is fine. The prediction is nonsense.

Measuring the fit with R squared

ss_res = np.sum((marks - pred) ** 2)
ss_tot = np.sum((marks - marks.mean()) ** 2)
r2 = 1 - ss_res / ss_tot
print("R squared:", round(r2, 4))      # about 0.996

R squared is the fraction of the variation in marks that the line accounts for. 0.996 is very high, and honestly, that is because this small dataset was chosen to be nearly a straight line. Real student data gives something like 0.3 to 0.6, and that is normal. Marks depend on many things, and hours studied is only one of them.

The same thing with scikit-learn

import numpy as np
from sklearn.linear_model import LinearRegression

hours = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=float).reshape(-1, 1)
marks = np.array([35, 42, 48, 51, 60, 63, 70, 74, 80, 86], dtype=float)

model = LinearRegression()
model.fit(hours, marks)

print("slope     :", round(model.coef_[0], 4))       # 5.5455
print("intercept :", round(model.intercept_, 4))     # 30.4
print("R squared :", round(model.score(hours, marks), 4))
print("predicted for 5.5 hours:", model.predict([[5.5]]).round(2))

Identical numbers. scikit-learn is not doing anything you did not just do by hand; it is doing it with better numerical stability and with more features.

More than one feature

import numpy as np
from sklearn.linear_model import LinearRegression

X = np.array([
    [1, 55], [2, 60], [3, 65], [4, 70], [5, 75],
    [6, 80], [7, 85], [8, 88], [9, 92], [10, 95],
], dtype=float)                      # hours, attendance
y = np.array([35, 42, 48, 51, 60, 63, 70, 74, 80, 86], dtype=float)

model = LinearRegression().fit(X, y)
print(dict(zip(["hours", "attendance"], model.coef_.round(3))))
print("intercept:", round(model.intercept_, 3))

Here the two features rise together, which is called multicollinearity. The coefficients become unstable: change one data point and they can swing wildly, or one can even turn negative. The predictions stay reasonable, but you must not read the individual coefficients as "the effect of attendance".

What linear regression cannot do

• It fits straight lines only. A curved relationship needs transformed

features or a different model.

• One extreme outlier pulls the whole line, because errors are squared.
• It finds association, not cause.