Linear Regression Explained — Math, Code, and Real-World Intuition

import numpy as np

# --- Dataset: house size (sq ft) vs sale price ($1000s) ---
house_sizes = np.array([750, 900, 1100, 1300, 1500, 1700, 1900, 2100, 2300, 2500], dtype=float)
sale_prices  = np.array([150, 175, 200, 230, 260, 290, 310, 340, 375, 400], dtype=float)

n = len(house_sizes)  # number of training examples

# --- Ordinary Least Squares: calculate slope (m) and intercept (b) ---
# Formula: m = (n*Σxy - Σx*Σy) / (n*Σx² - (Σx)²)

sum_x  = np.sum(house_sizes)           # Σx
sum_y  = np.sum(sale_prices)           # Σy
sum_xy = np.sum(house_sizes * sale_prices)  # Σxy
sum_x2 = np.sum(house_sizes ** 2)     # Σx²

# Numerator and denominator for slope
slope_numerator   = (n * sum_xy) - (sum_x * sum_y)
slope_denominator = (n * sum_x2) - (sum_x ** 2)

slope     = slope_numerator / slope_denominator   # m
intercept = (sum_y - slope * sum_x) / n           # b

print(f"Slope     (m): {slope:.4f}")
print(f"Intercept (b): {intercept:.4f}")
print(f"Model equation: price = {slope:.4f} * size + ({intercept:.4f})")

# --- Make a prediction ---
new_house_size = 1800  # sq ft
predicted_price = slope * new_house_size + intercept
print(f"\nPredicted price for {new_house_size} sq ft house: ${predicted_price:.1f}k")

# --- Calculate R² (coefficient of determination) ---
# R² tells you how much variance in y your model explains
predictions   = slope * house_sizes + intercept
ss_residual   = np.sum((sale_prices - predictions) ** 2)  # unexplained variance
ss_total      = np.sum((sale_prices - np.mean(sale_prices)) ** 2)  # total variance
r_squared     = 1 - (ss_residual / ss_total)
print(f"R² score: {r_squared:.4f}  (1.0 = perfect fit)")

import numpy as np
import matplotlib
matplotlib.use('Agg')  # use non-interactive backend for script execution
import matplotlib.pyplot as plt

# --- Synthetic dataset: study hours vs exam score ---
np.random.seed(42)
study_hours = np.linspace(1, 10, 50)                         # 50 points from 1 to 10 hours
exam_scores = 6.5 * study_hours + 20 + np.random.randn(50) * 5  # true slope=6.5, intercept=20, plus noise

# --- Normalise features to help gradient descent converge faster ---
# Without this, large feature values cause the gradients to explode
hours_mean = np.mean(study_hours)
hours_std  = np.std(study_hours)
study_hours_norm = (study_hours - hours_mean) / hours_std

# --- Initialise parameters ---
learning_rate  = 0.05
num_iterations = 300
n_samples      = len(study_hours_norm)

weight    = 0.0   # slope  — start at zero, gradient descent will adjust
bias      = 0.0   # intercept
loss_history = []  # track loss so we can plot convergence

# --- Gradient Descent Loop ---
for iteration in range(num_iterations):
    # Forward pass: compute predictions with current weight and bias
    predictions = weight * study_hours_norm + bias

    # Compute residuals (errors)
    residuals = predictions - exam_scores

    # Compute Mean Squared Error loss
    mse_loss = np.mean(residuals ** 2)
    loss_history.append(mse_loss)

    # Compute gradients (partial derivatives of MSE w.r.t weight and bias)
    # dL/dw = (2/n) * Σ(residual * x)  — how much loss changes as we nudge weight
    # dL/db = (2/n) * Σ(residual)      — how much loss changes as we nudge bias
    grad_weight = (2 / n_samples) * np.dot(residuals, study_hours_norm)
    grad_bias   = (2 / n_samples) * np.sum(residuals)

    # Update parameters: move opposite to the gradient (downhill)
    weight -= learning_rate * grad_weight
    bias   -= learning_rate * grad_bias

    if iteration % 50 == 0:
        print(f"Iteration {iteration:>3}: MSE Loss = {mse_loss:.4f} | w = {weight:.4f} | b = {bias:.4f}")

print(f"\nFinal weight (normalised): {weight:.4f}")
print(f"Final bias:               {bias:.4f}")

# --- De-normalise weight to interpret it in original units ---
original_scale_slope = weight / hours_std
original_scale_intercept = bias - (weight * hours_mean / hours_std)
print(f"\nReal-world slope:     {original_scale_slope:.4f} (exam points per study hour)")
print(f"Real-world intercept: {original_scale_intercept:.4f} (score with 0 hours studied)")

# --- Predict score for a student who studied 7 hours ---
hours_studied = 7.0
hours_norm    = (hours_studied - hours_mean) / hours_std
predicted_score = weight * hours_norm + bias
print(f"\nPredicted score for {hours_studied} hours of study: {predicted_score:.1f}/100")

import numpy as np
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

# --- Synthetic dataset: predict employee salary ---
# Features: years_experience, num_certifications, weekly_hours_worked, team_size
np.random.seed(7)
n_employees = 200

years_experience       = np.random.uniform(0, 20, n_employees)
num_certifications     = np.random.randint(0, 6, n_employees).astype(float)
weekly_hours_worked    = np.random.uniform(35, 60, n_employees)
team_size              = np.random.randint(2, 20, n_employees).astype(float)

# Ground truth: salary is influenced by all features (with some noise)
# True relationship: salary = 3500*exp + 4000*certs + 500*hours + 200*team + 30000 + noise
true_salary = (
    3500  * years_experience +
    4000  * num_certifications +
    500   * weekly_hours_worked +
    200   * team_size +
    30000 +
    np.random.randn(n_employees) * 5000  # realistic noise
)

# --- Build a tidy DataFrame ---
df = pd.DataFrame({
    'years_experience'    : years_experience,
    'num_certifications'  : num_certifications,
    'weekly_hours_worked' : weekly_hours_worked,
    'team_size'           : team_size,
    'annual_salary'       : true_salary
})

print("Dataset shape:", df.shape)
print(df.describe().round(1))

# --- Split into features and target ---
feature_columns = ['years_experience', 'num_certifications', 'weekly_hours_worked', 'team_size']
X = df[feature_columns].values
y = df['annual_salary'].values

# --- Train/test split: 80% train, 20% test ---
# shuffle=True is default — always use it for regression
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# --- Feature scaling: critical for gradient-based solvers ---
# StandardScaler transforms each feature to mean=0, std=1
# FIT only on training data — never peek at test data during fitting
scaler  = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)   # fit AND transform training set
X_test_scaled  = scaler.transform(X_test)        # only TRANSFORM test set (use training stats)

# --- Train the model ---
model = LinearRegression()
model.fit(X_train_scaled, y_train)

# --- Inspect learned parameters ---
print("\n--- Learned Model Parameters ---")
for feature_name, weight in zip(feature_columns, model.coef_):
    print(f"  {feature_name:<25}: {weight:>10,.2f}")
print(f"  {'intercept':<25}: {model.intercept_:>10,.2f}")

# --- Evaluate on test set ---
y_predictions = model.predict(X_test_scaled)

rmse      = np.sqrt(mean_squared_error(y_test, y_predictions))
r_squared = r2_score(y_test, y_predictions)

print(f"\n--- Model Performance on Held-Out Test Set ---")
print(f"  RMSE : ${rmse:>10,.0f}  (average prediction error in dollars)")
print(f"  R²   : {r_squared:>10.4f}  (proportion of salary variance explained)")

# --- Predict salary for a new employee ---
new_employee = pd.DataFrame([{
    'years_experience'    : 8.0,
    'num_certifications'  : 3.0,
    'weekly_hours_worked' : 45.0,
    'team_size'           : 10.0
}])
new_employee_scaled = scaler.transform(new_employee[feature_columns].values)
predicted_salary = model.predict(new_employee_scaled)[0]
print(f"\nPredicted salary for new employee profile: ${predicted_salary:,.0f}/year")

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.metrics import r2_score

np.random.seed(0)

# --- Scenario: electricity consumption vs temperature ---
# True relationship is quadratic (U-shaped) — think AC in summer AND heating in winter
temperature_c    = np.linspace(-10, 40, 100)
consumption_kwh  = 0.8 * temperature_c**2 - 20 * temperature_c + 300 + np.random.randn(100) * 15

X = temperature_c.reshape(-1, 1)  # sklearn expects 2D input
y = consumption_kwh

# --- Model 1: Naive linear regression ---
linear_model = LinearRegression()
linear_model.fit(X, y)
y_pred_linear = linear_model.predict(X)
r2_linear     = r2_score(y, y_pred_linear)

residuals_linear = y - y_pred_linear

print("=== Linear Model ===")
print(f"R² score: {r2_linear:.4f}")
print(f"Max residual: {np.max(np.abs(residuals_linear)):.1f} kWh  ← suspiciously large")
print(f"Residual pattern: {'U-shaped — model is missing the curve' if r2_linear < 0.7 else 'OK'}")

# --- Model 2: Polynomial regression (degree 2) ---
# make_pipeline chains preprocessing + model cleanly
poly_model = make_pipeline(
    PolynomialFeatures(degree=2, include_bias=False),  # adds x² as a feature
    LinearRegression()
)
poly_model.fit(X, y)
y_pred_poly  = poly_model.predict(X)
r2_poly      = r2_score(y, y_pred_poly)

residuals_poly = y - y_pred_poly

print("\n=== Polynomial Model (degree 2) ===")
print(f"R² score: {r2_poly:.4f}  ← much better")
print(f"Max residual: {np.max(np.abs(residuals_poly)):.1f} kWh")
print(f"Residual pattern: {'Random — model captured the curve correctly' if r2_poly > 0.95 else 'Still biased'}")

# --- Check residual distribution ---
print("\n=== Residual Statistics ===")
print(f"Linear  residuals — mean: {np.mean(residuals_linear):>7.2f}, std: {np.std(residuals_linear):.2f}")
print(f"Poly    residuals — mean: {np.mean(residuals_poly):>7.2f}, std: {np.std(residuals_poly):.2f}")
print("\nPro tip: residuals should have mean ≈ 0 and show no pattern when plotted against predictions.")
print("If they fan out or curve, your model has a systematic error your metrics are hiding.")