Regression: Simple, Multiple & Polynomial

Simple Linear Regression

Simple Linear Regression models the relationship between one feature $x$ and a continuous target $y$ as a straight line. It's the foundation of everything in supervised learning.

Simple Linear Model$$\hat{y} = w_1 x + b$$$w_1$: slope (weight). $b$: intercept (bias). Goal: find values that minimise prediction error.

MSE Loss & Closed-Form Solution$$J(w,b)=\frac{1}{m}\sum_{i=1}^m(\hat{y}^{(i)}-y^{(i)})^2 \qquad w^*=\frac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sum(x_i-\bar{x})^2},\quad b^*=\bar{y}-w^*\bar{x}$$

Python

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score

# Synthetic data: house size vs price
np.random.seed(42)
X_simple = np.random.uniform(500, 3000, 100).reshape(-1, 1)   # sq ft
y = 150 * X_simple.squeeze() + 50000 + np.random.randn(100) * 30000

model = LinearRegression()
model.fit(X_simple, y)

print(f"Slope ($/sqft):  {model.coef_[0]:.2f}")
print(f"Intercept:       {model.intercept_:.2f}")
print(f"R² score:        {r2_score(y, model.predict(X_simple)):.4f}")

plt.figure(figsize=(8, 4))
plt.scatter(X_simple, y, alpha=0.5, color="#4a8fa8", label="data")
x_line = np.linspace(500, 3000, 100).reshape(-1, 1)
plt.plot(x_line, model.predict(x_line), color="#b85c2a", lw=2, label="fit")
plt.xlabel("Size (sq ft)"); plt.ylabel("Price ($)")
plt.title("Simple Linear Regression"); plt.legend(); plt.tight_layout(); plt.show()

Multiple Linear Regression

Extend to $n$ features. Each feature gets its own weight. The prediction is a dot product between the weight vector and the feature vector:

Multiple Linear Model$$\hat{y} = \mathbf{w}^\top\mathbf{x} + b = w_1x_1 + w_2x_2 + \cdots + w_nx_n + b$$

Normal Equation (Closed Form)$$\mathbf{w}^* = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}$$Exact solution, but $O(n^3)$ — use gradient descent for $n > 10^4$ features.

Python

import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error, r2_score
import numpy as np

# Simulate a housing dataset
np.random.seed(0)
n = 500
df = pd.DataFrame({
    "size":       np.random.uniform(500, 3500, n),
    "bedrooms":   np.random.randint(1, 6, n),
    "age":        np.random.randint(0, 50, n),
    "dist_center":np.random.uniform(1, 30, n),
})
df["price"] = (120*df["size"] + 15000*df["bedrooms"]
               - 800*df["age"] - 3000*df["dist_center"]
               + 40000 + np.random.randn(n)*20000)

X = df.drop("price", axis=1)
y = df["price"]

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

scaler = StandardScaler()
X_train_s = scaler.fit_transform(X_train)
X_test_s  = scaler.transform(X_test)

model = LinearRegression()
model.fit(X_train_s, y_train)
y_pred = model.predict(X_test_s)

print(f"RMSE: ${np.sqrt(mean_squared_error(y_test, y_pred)):,.0f}")
print(f"R²:   {r2_score(y_test, y_pred):.4f}")

# Feature importance via coefficients (after scaling)
coef_df = pd.Series(model.coef_, index=X.columns).sort_values(key=abs, ascending=False)
print("
Feature coefficients (scaled):")
print(coef_df)

⚠️

Assumptions of linear regression: linearity between features and target, homoscedasticity (constant error variance), no multicollinearity (features not highly correlated with each other), and normally distributed residuals. Always check these with residual plots before trusting your model.

Polynomial Regression

Linear regression can only fit straight lines (or hyperplanes). For curved relationships, polynomial regression adds powers of the original features as new columns — then fits a standard linear model on top.

Polynomial Features (degree $d$)$$x \to [1,\, x,\, x^2,\, x^3,\, \ldots,\, x^d]$$Still a linear model — linear in the parameters $w_0,w_1,\ldots,w_d$. But can fit nonlinear curves.

Python

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score

# Non-linear data
np.random.seed(7)
X = np.sort(np.random.uniform(-3, 3, 80)).reshape(-1, 1)
y = 0.5*X.squeeze()**3 - 2*X.squeeze() + np.random.randn(80) * 1.5

fig, axes = plt.subplots(1, 3, figsize=(15, 4), sharey=True)
x_plot = np.linspace(-3, 3, 200).reshape(-1, 1)

for ax, degree, color in zip(axes, [1, 3, 10], ["#4a8fa8","#5c8a58","#b85c2a"]):
    pipe = Pipeline([
        ("poly",  PolynomialFeatures(degree=degree, include_bias=False)),
        ("scaler", __import__("sklearn.preprocessing", fromlist=["StandardScaler"]).StandardScaler()),
        ("model", LinearRegression()),
    ])
    pipe.fit(X, y)
    cv = cross_val_score(pipe, X, y, cv=5, scoring="neg_mean_squared_error")
    rmse = np.sqrt(-cv.mean())

    ax.scatter(X, y, alpha=0.5, color=color, s=20)
    ax.plot(x_plot, pipe.predict(x_plot), color=color, lw=2)
    ax.set_title(f"Degree {degree}
CV RMSE: {rmse:.2f}")
    ax.set_xlabel("x")

axes[0].set_ylabel("y")
fig.suptitle("Polynomial Regression — Underfitting vs. Overfitting", fontsize=12)
plt.tight_layout(); plt.show()
# Degree 1 underfits, Degree 3 fits well, Degree 10 overfits

Ridge, Lasso and ElasticNet

Polynomial features increase model complexity fast — regularisation is essential to prevent overfitting.

Python

from sklearn.linear_model import Ridge, Lasso, ElasticNet, RidgeCV
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score
import numpy as np

# High-degree polynomial + Ridge regularisation
pipe_ridge = Pipeline([
    ("poly",   PolynomialFeatures(degree=8, include_bias=False)),
    ("scaler", StandardScaler()),
    ("model",  Ridge(alpha=1.0)),          # alpha = lambda (regularisation strength)
])
pipe_lasso = Pipeline([
    ("poly",   PolynomialFeatures(degree=8, include_bias=False)),
    ("scaler", StandardScaler()),
    ("model",  Lasso(alpha=0.01)),         # Lasso drives some coefficients to exactly 0
])

for name, pipe in [("Ridge", pipe_ridge), ("Lasso", pipe_lasso)]:
    cv = cross_val_score(pipe, X, y, cv=5, scoring="neg_mean_squared_error")
    print(f"{name}: CV RMSE = {np.sqrt(-cv.mean()):.4f}")

Residual Analysis

Python

import matplotlib.pyplot as plt
import scipy.stats as stats

model.fit(X_train_s, y_train)
residuals = y_test.values - model.predict(X_test_s)

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
# Residuals vs. Fitted (should look random — no pattern)
axes[0].scatter(model.predict(X_test_s), residuals, alpha=0.5, color="#4a8fa8")
axes[0].axhline(0, color="red", ls="--"); axes[0].set_title("Residuals vs. Fitted")
axes[0].set_xlabel("Fitted values"); axes[0].set_ylabel("Residuals")
# Histogram (should look Gaussian)
axes[1].hist(residuals, bins=25, color="#5c8a58", edgecolor="white")
axes[1].set_title("Residual Distribution"); axes[1].set_xlabel("Residual")
# Q-Q plot (points should lie on the diagonal line)
stats.probplot(residuals, plot=axes[2]); axes[2].set_title("Normal Q-Q Plot")
plt.tight_layout(); plt.show()

Summary

Simple LR: $\hat{y}=wx+b$. Minimise MSE. Closed-form solution via Normal Equation.
Multiple LR: extends to $n$ features. Check for multicollinearity and heteroscedasticity.
Polynomial regression: add $x^2, x^3, \ldots$ as features then fit linear model. Always regularise.
Ridge shrinks all weights; Lasso zeroes some out. Use Ridge as default, Lasso for feature selection.
Residual plots are non-negotiable — always verify assumptions after fitting.