02. Mastering Linear Regression

ŷ = w₀ + w₁x₁ + w₂x₂ + ... + wₚxₚ

L(w) = (1/n) Σ(yᵢ - ŷᵢ)² = (1/n)||y - Xw||²

w = (XᵀX)⁻¹Xᵀy

class LinearRegressionScratch:
    def fit(self, X, y):
        # Add bias term (add column of ones to X)
        X_b = np.c_[np.ones((len(X), 1)), X]
        # Normal equation
        w = np.linalg.inv(X_b.T @ X_b) @ X_b.T @ y
        self.bias = w[0]
        self.weights = w[1:]
        return self
 
    def predict(self, X):
        return X @ self.weights + self.bias

w_{t+1} = w_t - α * ∇L(w_t)

∇L(w) = -(2/n) Xᵀ(y - Xw)

class LinearRegressionGD:
    def fit(self, X, y, lr=0.01, n_iter=1000):
        self.weights = np.zeros(X.shape[1])
        self.bias = 0
        self.loss_history = []
 
        for _ in range(n_iter):
            y_pred = X @ self.weights + self.bias
            error = y_pred - y
 
            # Compute gradients
            dw = (1/len(y)) * (X.T @ error)
            db = (1/len(y)) * np.sum(error)
 
            # Update weights
            self.weights -= lr * dw
            self.bias -= lr * db
 
            # Record loss
            loss = (1/(2*len(y))) * np.sum((y - y_pred)**2)
            self.loss_history.append(loss)
 
        return self

L1 (Lasso): Diamond-shaped constraint
  → Contacts loss function contour at vertex → Some coefficients exactly 0

L2 (Ridge): Circular constraint
  → Contacts loss function contour on curve → All coefficients shrink but not 0

from sklearn.linear_model import RidgeCV, LassoCV
 
# Auto-search optimal alpha with Cross-Validation
ridge_cv = RidgeCV(alphas=np.logspace(-4, 4, 50), cv=5)
ridge_cv.fit(X_train_scaled, y_train)
print(f"Best alpha: {ridge_cv.alpha_}")
 
lasso_cv = LassoCV(alphas=np.logspace(-4, 1, 50), cv=5, max_iter=10000)
lasso_cv.fit(X_train_scaled, y_train)
print(f"Best alpha: {lasso_cv.alpha_}")

VIF = 1 / (1 - R²)

def calculate_vif(X):
    vif_data = []
    for col in X.columns:
        X_temp = X.drop(columns=[col])
        r2 = LinearRegression().fit(X_temp, X[col]).score(X_temp, X[col])
        vif = 1 / (1 - r2) if r2 < 1 else float('inf')
        vif_data.append({'Feature': col, 'VIF': vif})
    return pd.DataFrame(vif_data)

from scipy import stats
 
residuals = y_test - y_pred
 
# Basic statistics
print(f"Mean: {residuals.mean():.6f}")  # Should be close to 0
print(f"Skewness: {stats.skew(residuals):.4f}")  # Should be close to 0
print(f"Kurtosis: {stats.kurtosis(residuals):.4f}")  # Should be close to 0
 
# Normality test
_, p_value = stats.shapiro(residuals[:5000])
print(f"Shapiro-Wilk p-value: {p_value:.6f}")

from sklearn.linear_model import LinearRegression, Ridge, Lasso, RidgeCV
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import r2_score, mean_squared_error
 
# Scaling (required for GD and regularization!)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
 
# Train and compare models
models = {
    'OLS': LinearRegression(),
    'Ridge': Ridge(alpha=1.0),
    'Lasso': Lasso(alpha=0.1)
}
 
for name, model in models.items():
    model.fit(X_train_scaled, y_train)
    y_pred = model.predict(X_test_scaled)
    print(f"{name}: R²={r2_score(y_test, y_pred):.4f}, "
          f"RMSE={np.sqrt(mean_squared_error(y_test, y_pred)):.4f}")