Ordinary Least Squares (OLS)

Ordinary Least Squares (OLS)#

Ordinary Least Squares (OLS) From Scratch#

In this notebook, we’ll implement Ordinary Least Squares regression from scratch using NumPy, and visualize how it behaves under different levels of noise in the data.

Overview#

We’ll follow these steps:

  1. Generate synthetic data: \( y = 3x + 2 + \varepsilon \)

  2. Implement OLS manually via matrix operations.

  3. Visualize the linear fit for varying noise levels.

  4. Discuss limitations of OLS with high noise, outliers, and model misspecification.

import numpy as np
import matplotlib.pyplot as plt

# Set seed for reproducibility
np.random.seed(42)

def generate_data(n_samples=100, noise_std=1.0):
    x = np.linspace(0, 10, n_samples)
    ######## Your code here ###########
    y =
    ##################################
    return x.reshape(-1, 1), y


# Stack 1s in first column of Beta matrix - learn the intercept
def add_intercept(X):
    return np.hstack([np.ones((X.shape[0], 1)), X])

def ols_fit(X, y):
    # OLS: beta = (X^T X)^(-1) X^T y
    ######## Your code here ###########
    beta =
    ##################################
    return beta

def predict(X, beta):
    ######## Your code here ###########
    y_hat =
    ##################################
    return y_hat

def r2_score(y_true, y_pred):
    # R²: 1 - (SS_res / SS_tot)
    # SS_res = ∑(yᵢ - ŷᵢ)²
    # SS_tot = ∑(yᵢ - ȳ)²
    ######## Your code here ###########
    ss_res =
    ss_tot =
    ##################################
    return 1 - ss_res / ss_tot

def plot_fit(x, y, y_pred, title):
    plt.scatter(x, y, label='Data', alpha=0.6)
    plt.plot(x, y_pred, color='red', label='OLS fit', linewidth=2)
    plt.title(title)
    plt.xlabel("x")
    plt.ylabel("y")
    plt.legend()
    plt.grid(True)
    plt.show()

# Try different noise levels
for noise in [0.1, 1.0, 5.0, 15.0]:
    x, y = generate_data(noise_std=noise)
    X_design = add_intercept(x)
    beta_hat = ols_fit(X_design, y)
    y_pred = predict(X_design, beta_hat)

    r2 = r2_score(y, y_pred)
    title = f"Noise Std = {noise} | Coeffs: {beta_hat.round(2)} | R² = {r2:.4f}"
    plot_fit(x, y, y_pred, title)

Additional Depth#

You can find additional depth (if interested) here:

https://colab.research.google.com/github/fbeilstein/machine_learning/blob/master/lecture_09_linear_regression.ipynb#scrollTo=28nzh_ZhkQh1

Contents