Linear Regression from Scratch

Difficulty: Advanced

Linear regression is one of the most fundamental algorithms in machine learning and statistics. It models the relationship between a dependent variable (y) and one or more independent variables (x) by fitting a straight line through the data. The equation of a simple linear regression line is y = mx + b, where m is the slope (how much y changes for each unit change in x) and b is the y-intercept (the value of y when x is 0).

The Ordinary Least Squares (OLS) method finds the line that minimizes the sum of squared differences between the actual y values and the predicted y values. The formulas for slope and intercept are derived from calculus: m = sum((xi - x_mean) * (yi - y_mean)) / sum((xi - x_mean)^2), and b = y_mean - m * x_mean. These closed-form solutions make linear regression fast to compute even on large datasets.

R-squared (coefficient of determination) measures how well the regression line fits the data. It is computed as R^2 = 1 - (SS_res / SS_tot), where SS_res is the sum of squared residuals (errors between predictions and actuals) and SS_tot is the total sum of squares (variance in the actual values around their mean). An R^2 of 1.0 means the model perfectly explains all variance; an R^2 of 0.0 means the model is no better than simply predicting the mean.

The cost function (also called loss function) in linear regression is typically the Mean Squared Error (MSE), which is the average of the squared differences between predictions and actual values. Gradient descent is an iterative optimization algorithm that can also be used to find the optimal slope and intercept by repeatedly adjusting them in the direction that reduces the MSE. While OLS gives an exact solution, gradient descent is more general and scales to problems where closed-form solutions do not exist.

Linear regression assumes a linear relationship between features and target, independence of observations, homoscedasticity (constant variance of residuals), and normally distributed residuals. Violations of these assumptions can lead to unreliable predictions. In practice, you should always visualize residuals and check these assumptions before trusting the model.

Code examples

Simple Linear Regression with OLS

def linear_regression(x, y):
    n = len(x)
    x_mean = sum(x) / n
    y_mean = sum(y) / n
    
    numerator = sum((xi - x_mean) * (yi - y_mean) for xi, yi in zip(x, y))
    denominator = sum((xi - x_mean) ** 2 for xi in x)
    
    slope = numerator / denominator
    intercept = y_mean - slope * x_mean
    return slope, intercept

# Dataset: hours studied vs exam score
hours = [1, 2, 3, 4, 5]
scores = [52, 56, 61, 64, 69]

slope, intercept = linear_regression(hours, scores)
print(f"Slope: {slope:.2f}")
print(f"Intercept: {intercept:.2f}")
print(f"Equation: y = {slope:.2f}x + {intercept:.2f}")

# Predict for 6 hours
prediction = slope * 6 + intercept
print(f"Predicted score for 6 hours: {prediction:.2f}")

We compute the slope and intercept using the OLS formulas. x_mean = 3, y_mean = (52+56+61+64+69)/5 = 302/5 = 60.4. Numerator = (-2)(-8.4)+(-1)(-4.4)+(0)(0.6)+(1)(3.6)+(2)(8.6) = 16.8+4.4+0+3.6+17.2 = 42.0. Denominator = 4+1+0+1+4 = 10. Slope = 42/10 = 4.2. Intercept = 60.4 - 4.2*3 = 60.4 - 12.6 = 47.8. Prediction for 6 hours: 4.2*6 + 47.8 = 73.0.

Computing R-Squared

def r_squared(x, y, slope, intercept):
    y_mean = sum(y) / len(y)
    ss_tot = sum((yi - y_mean) ** 2 for yi in y)
    ss_res = sum((yi - (slope * xi + intercept)) ** 2 for xi, yi in zip(x, y))
    return 1 - (ss_res / ss_tot)

hours = [1, 2, 3, 4, 5]
scores = [52, 56, 61, 64, 69]

# Using previously computed values
slope = 4.2
intercept = 47.8

r2 = r_squared(hours, scores, slope, intercept)
print(f"R-squared: {r2:.4f}")
print(f"The model explains {r2 * 100:.1f}% of the variance in scores.")

R-squared of 0.9955 means the linear model explains 99.5% of the variation in exam scores based on hours studied. SS_tot = (52-60.4)^2 + (56-60.4)^2 + (61-60.4)^2 + (64-60.4)^2 + (69-60.4)^2 = 70.56+19.36+0.36+12.96+73.96 = 177.2. The exact OLS values slope=4.2, intercept=47.8 give predictions: 52.0, 56.2, 60.4, 64.6, 68.8. SS_res = 0+0.04+0.36+0.36+0.04 = 0.8. R^2 = 1 - 0.8/177.2 = 0.9955.

Gradient Descent for Linear Regression

def gradient_descent(x, y, lr=0.01, epochs=1000):
    m, b = 0.0, 0.0  # Initialize slope and intercept
    n = len(x)
    
    for epoch in range(epochs):
        y_pred = [m * xi + b for xi in x]
        dm = (-2/n) * sum((yi - yp) * xi for xi, yi, yp in zip(x, y, y_pred))
        db = (-2/n) * sum(yi - yp for yi, yp in zip(y, y_pred))
        m -= lr * dm
        b -= lr * db
    
    return m, b

x = [1, 2, 3, 4, 5]
y = [2, 4, 5, 4, 5]

slope, intercept = gradient_descent(x, y, lr=0.01, epochs=1000)
print(f"Slope: {slope:.2f}")
print(f"Intercept: {intercept:.2f}")
predictions = [round(slope * xi + intercept, 2) for xi in x]
print(f"Predictions: {predictions}")

Gradient descent iteratively adjusts the slope and intercept to minimize MSE. Starting from (0, 0), it converges to approximately slope=0.62 and intercept=2.14 after 1000 epochs. The OLS solution for this data is m=0.6 and b=2.2, so gradient descent is very close but not exact due to the finite number of iterations and learning rate.

Making Predictions and Computing Residuals

# Given a fitted model
slope = 2.5
intercept = 1.0

# Data points
x_values = [1, 2, 3, 4, 5]
y_actual = [3.5, 6.2, 8.0, 10.5, 13.1]

# Predictions and residuals
print("x | actual | predicted | residual")
print("-" * 40)
residuals = []
for x, y in zip(x_values, y_actual):
    y_pred = slope * x + intercept
    residual = y - y_pred
    residuals.append(residual)
    print(f"{x} | {y:6.1f}  | {y_pred:9.1f} | {residual:+.1f}")

mse = sum(r ** 2 for r in residuals) / len(residuals)
print(f"\nMSE: {mse:.2f}")

For each data point, we compute the predicted value using y = 2.5x + 1.0, then the residual (actual - predicted). The residuals show the model slightly underpredicts for x=2 and overpredicts for x=3,4,5. MSE = (0^2 + 0.2^2 + 0.5^2 + 0.5^2 + 0.4^2)/5 = (0 + 0.04 + 0.25 + 0.25 + 0.16)/5 = 0.70/5 = 0.14.

Key points

Concepts covered

Linear regression, Slope and intercept, Least squares method, R-squared, Cost function