Implementing Linear Regression Algorithm from Scratch in Python

Linear regression is virtually the simplest ML algorithm but is widely used from predictive analytics to hypothesis testing. While its assumptions may not always hold in practice, understanding its fundamentals provides a strong foundation for tackling more complex models.

The linear regression model is based on the linear equation

Equation for linear regression

where y is the dependent variable, X is the independent variable aka feature vector, w is the slope aka weight and b is the intercept or bias. It assumes that the dependent variable y is linearly dependent on X.

To build the model (that is, find optimal w and b), the data is fit to this equation by optimizing a goodness of fit metric. Popular metrics for evaluating the model are:

• Sum of Squared Residuals (SSR)
• Mean Squared Error (MSE)
• R² (square of the Pearson correlation coefficient — between 0 and 1 — the higher, the better)
• p-value for parameter significance

Let’s jump right in to coding it!

Generate Some Linearly Correlated Data

I am using a weight of 2 and a randomized bias for each datapoint.

import numpy as np
import matplotlib.pyplot as plt

std = 10
X = np.linspace(0, 100, 50, endpoint=True)
y = 2*X + np.random.normal(0, std, size=len(X)) 

The scatter plot for the data

Initialize Weight and Bias

I am initializing them based on the slope from highest and lowest given data points (supposedly the slope), and y for X ~ 0 (supposedly the intercept), respectively. They could just as well be initialized as zeros.

weight = (np.max(y) - np.min(y)) / (X[np.argmax(y)] - X[np.argmin(y)])
bias = y[np.min(np.where(np.round(X) == 0)[0])]

Define Function to Calculate Mean Squared Error

For goodness of fit, I am choosing the Mean Squared Error (MSE). Other metrics could also have been chosen.

The following function takes the weight w and the bias b, and calculates the MSE as

Equation for mean squared error

where n is the number of training points.

def MSE(weight, bias):
    pred = np.dot(X, weight) + bias
    MSE = np.average((pred - y)**2)
    return MSE

Define Function for Gradient Descent

To optimize the above goodness of fit metric, we need to find the minima of its gradient with respect to the variables we’re trying to find the optimal value for (that is, w and b). With some calculus, you would get the following equations

The gradient of MSE with respect to the weight

The gradient of MSE with respect to the bias

where y-hat is the prediction based on the linear model.

Get Sunia Tanweer’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

The following function finds the derivatives, and modifies the weights and bias based on the derivatives’ values, while using a learning rate for damping the speed of the descent.

def grad_descent(w, b, lr=0.0001):
    pred = np.dot(X, weight) + bias
    
    ddw = -2*np.average(X*(y - pred))
    ddb = -2*np.average(y - pred)

    w = w - lr*ddw
    b = b - lr*ddb

    err = MSE(w, b)

    return w, b, err

Run All Epochs to Find Optimal Parameter Values

Now, we just need to run the grad_descent function above for a number of training rounds (called epochs) to find the optimal value for the weight and bias. Technically, this step is what’s called grad_descent.

for epoch in range(10000):
    weight, bias, error = grad_descent(weight, bias)
    if epoch % 1000 == 0:
        print(f'Iteration {epoch}, MSE {error}.')

Iteration 0, MSE 81.09709273798977.
Iteration 1000, MSE 69.7282969059099.
Iteration 2000, MSE 69.72566892765921.
Iteration 3000, MSE 69.72329818187762.
Iteration 4000, MSE 69.7211594900678.
Iteration 5000, MSE 69.71923013826103.
Iteration 6000, MSE 69.71748963578318.
Iteration 7000, MSE 69.71591949763376.
Iteration 8000, MSE 69.71450304816604.
Iteration 9000, MSE 69.71322524398347.

The MSE has stabilized showing that the parameters are now optimized and we can plot the best fit line on top of the data used for training.

Make Predictions

pred = np.dot(X, weight) + bias

Training data and the fitted regression model

Bringing it all together in a class

You might want to keep using your linear regression code for different datasets. So a neater way of doing it is by creating a class.

class LinearRegressionEstimator():
    def __init__(self, lr=0.0001, epochs=10000):
        self.lr = lr
        self.epochs = epochs
        self.w = None
        self.b = None
        self.error = None
    
    def initialize_parameters(self, X, y):

        self.w = (np.max(y) - np.min(y)) / (X[np.argmax(y)] - X[np.argmin(y)])

        indices = np.where(np.round(X) == 0)[0]
        if indices.size > 0:
            self.b = y[np.min(indices)]
        else:
            self.b = np.average(y)

        self.error = 0

    def update_MSE(self, X, y):
        pred = self.predict(X)
        self.error = np.average((pred - y)**2)
    
    def grad_descent(self, X, y):
        pred = self.predict(X)
        
        ddw = -2*np.average(X*(y - pred))
        ddb = -2*np.average(y - pred)

        self.w = self.w - self.lr*ddw
        self.b = self.b - self.lr*ddb

        self.update_MSE(X, y)
    
    def fit(self, X, y):

        self.initialize_parameters(X, y)

        for _ in range(self.epochs):
            self.grad_descent(X, y)

    def predict(self, X):
        return np.dot(X, self.w) + self.b
    
    def return_parameters(self):
        print(f'Weight is {self.w}, bias is {self.b}, MSE is {self.error}.')
        return self.w, self.b

Running the same dataset through this class —

• Initialize the class
• Use the fit method to optimize parameters for the specific data
• Use predict method to get predictions for the train or test data

linreg = LinearRegressionEstimator(lr=0.0001, epochs=10000)
linreg.fit(X, y)

weight, bias = linreg.return_parameters()

Weight is 1.9844443976822137, bias is 1.477535428810894, MSE is 21.71950809693673.

pred = linreg.predict(X)

Get Some New Data for Testing

(or you could have just split the data into test and train initially, instead of doing this — for a real dataset, split the data first and use the train portion to fit the model and the test portion for testing)

X_test = np.linspace(105, 130, 7, endpoint=True)
y_test = 2*X_test + np.random.normal(0, std, size=len(X_test))

predictions = linreg.predict(X_test)
print(f"MSE: {np.average((predictions - y_test)**2)}")

MSE: 71.94673813442024

Training data, test data, linear regression fitted line and the test predictions

For this toy-dataset, the model fits great since the data is very linearly correlated. In practice, this isn’t always true but linear regression is still widely used as a baseline for advanced ML algorithms — making it an important tool to master.

You can find the notebook on my github.