Chapter 5: Multiple & Polynomial Regression

import numpy as np
import matplotlib.pyplot as plt

class MultipleLinearRegression:
    """Multiple Linear Regression using the Normal Equation.
    
    Solves: W = (XᵀX)⁻¹XᵀY exactly — no iterations needed.
    """
    
    def __init__(self):
        self.weights = None   # includes bias as w[0]
        self.r_squared = None
        self.adj_r_squared = None
    
    def fit(self, X, y):
        """Fit the model using the Normal Equation.
        
        Parameters:
            X: numpy array of shape (m, n) — feature matrix
            y: numpy array of shape (m,) — target values
        """
        m, n = X.shape
        
        # Add bias column (column of 1s)
        X_b = np.column_stack([np.ones(m), X])
        
        # Normal Equation: W = (XᵀX)⁻¹XᵀY
        XtX = X_b.T @ X_b            # (n+1) x (n+1)
        XtY = X_b.T @ y              # (n+1) x 1
        self.weights = np.linalg.solve(XtX, XtY)  # more stable than inv()
        
        # Compute R² and Adjusted R²
        y_pred = X_b @ self.weights
        ss_res = np.sum((y - y_pred) ** 2)
        ss_tot = np.sum((y - np.mean(y)) ** 2)
        self.r_squared = 1 - ss_res / ss_tot
        self.adj_r_squared = 1 - ((1 - self.r_squared) * (m - 1)) / (m - n - 1)
        
        return self
    
    def predict(self, X):
        """Predict target values for new data."""
        m = X.shape[0]
        X_b = np.column_stack([np.ones(m), X])
        return X_b @ self.weights
    
    @property
    def intercept(self):
        return self.weights[0]
    
    @property
    def coefficients(self):
        return self.weights[1:]


# ===== DEMO: House Price Prediction =====
np.random.seed(42)
m = 100  # samples

# Generate synthetic data: price = 5 + 3*area + 8*rooms + noise
area = np.random.uniform(5, 20, m)       # 100 sq ft units
rooms = np.random.randint(1, 6, m).astype(float)
noise = np.random.normal(0, 3, m)
price = 5 + 3 * area + 8 * rooms + noise

X = np.column_stack([area, rooms])
y = price

# Fit the model
model = MultipleLinearRegression()
model.fit(X, y)

print("=" * 50)
print("MULTIPLE LINEAR REGRESSION RESULTS")
print("=" * 50)
print(f"Intercept (w₀):         {model.intercept:.4f}")
print(f"Area coefficient (w₁):  {model.coefficients[0]:.4f}")
print(f"Rooms coefficient (w₂): {model.coefficients[1]:.4f}")
print(f"R²:                     {model.r_squared:.6f}")
print(f"Adjusted R²:            {model.adj_r_squared:.6f}")
print(f"\nTrue coefficients: w₀=5, w₁=3, w₂=8")

# Predictions
y_pred = model.predict(X)
print(f"\nFirst 5 predictions vs actual:")
for i in range(5):
    print(f"  Predicted: {y_pred[i]:.2f}, Actual: {y[i]:.2f}, "
          f"Error: {abs(y[i]-y_pred[i]):.2f}")

import numpy as np
import matplotlib.pyplot as plt

class PolynomialRegression:
    """Polynomial Regression from scratch.
    
    Transforms x → [1, x, x², ..., xᵈ] then applies Normal Equation.
    """
    
    def __init__(self, degree=2):
        self.degree = degree
        self.weights = None
    
    def _create_features(self, x):
        """Create polynomial feature matrix."""
        m = len(x)
        X = np.ones((m, self.degree + 1))
        for d in range(1, self.degree + 1):
            X[:, d] = x ** d
        return X
    
    def fit(self, x, y):
        """Fit polynomial to data."""
        X = self._create_features(x)
        # Normal equation: W = (XᵀX)⁻¹XᵀY
        self.weights = np.linalg.solve(X.T @ X, X.T @ y)
        return self
    
    def predict(self, x):
        """Predict using fitted polynomial."""
        X = self._create_features(x)
        return X @ self.weights
    
    def r_squared(self, x, y):
        """Compute R² score."""
        y_pred = self.predict(x)
        ss_res = np.sum((y - y_pred) ** 2)
        ss_tot = np.sum((y - np.mean(y)) ** 2)
        return 1 - ss_res / ss_tot


# ===== OVERFITTING DEMONSTRATION =====
np.random.seed(42)

# True function: y = 0.5x² - 2x + 1 + noise
x_train = np.sort(np.random.uniform(-3, 3, 15))
y_true = 0.5 * x_train**2 - 2 * x_train + 1
y_train = y_true + np.random.normal(0, 1.5, len(x_train))

x_test = np.linspace(-3.5, 3.5, 200)

degrees = [1, 2, 3, 5, 12]
print("=" * 55)
print(f"{'Degree':<10} {'Train R²':<15} {'Overfit?'}")
print("=" * 55)

fig, axes = plt.subplots(1, len(degrees), figsize=(20, 4))
for ax, deg in zip(axes, degrees):
    model = PolynomialRegression(degree=deg)
    model.fit(x_train, y_train)
    r2 = model.r_squared(x_train, y_train)
    
    overfit = "✓ YES" if deg > 4 else "  No"
    print(f"  {deg:<10} {r2:<15.6f} {overfit}")
    
    # Plot
    y_fit = model.predict(x_test)
    ax.scatter(x_train, y_train, c='#059669', s=40, zorder=5)
    ax.plot(x_test, y_fit, color='#0891b2', linewidth=2)
    ax.set_title(f'Degree {deg} (R²={r2:.3f})')
    ax.set_ylim(-15, 25)
    ax.set_xlabel('x')

plt.tight_layout()
plt.savefig('overfitting_demo.png', dpi=150, bbox_inches='tight')
plt.show()
print("\n⚠ Notice: Degree 12 has R²≈1.0 on training but will fail on test data!")

import numpy as np

def compute_vif(X, feature_names=None):
    """Compute Variance Inflation Factor for each feature.
    
    VIF_j = 1 / (1 - R²_j)
    where R²_j is from regressing feature j on all other features.
    """
    n_features = X.shape[1]
    vif_values = []
    
    if feature_names is None:
        feature_names = [f"x{i+1}" for i in range(n_features)]
    
    for j in range(n_features):
        # Target: feature j
        y_j = X[:, j]
        # Predictors: all other features
        X_others = np.delete(X, j, axis=1)
        
        # Add bias column
        m = X_others.shape[0]
        X_b = np.column_stack([np.ones(m), X_others])
        
        # Fit regression: feature j ~ all others
        w = np.linalg.solve(X_b.T @ X_b, X_b.T @ y_j)
        y_pred = X_b @ w
        
        # R² of this regression
        ss_res = np.sum((y_j - y_pred) ** 2)
        ss_tot = np.sum((y_j - np.mean(y_j)) ** 2)
        r2_j = 1 - ss_res / ss_tot
        
        # VIF
        vif_j = 1.0 / (1.0 - r2_j) if r2_j < 1 else float('inf')
        vif_values.append(vif_j)
    
    # Print results
    print("=" * 55)
    print("VARIANCE INFLATION FACTOR (VIF) ANALYSIS")
    print("=" * 55)
    print(f"{'Feature':<15} {'VIF':<10} {'R²_j':<10} {'Status'}")
    print("-" * 55)
    for name, vif in zip(feature_names, vif_values):
        r2_j = 1 - 1/vif
        if vif > 10:
            status = "❌ SEVERE"
        elif vif > 5:
            status = "⚠️  HIGH"
        elif vif > 2:
            status = "🔶 MODERATE"
        else:
            status = "✅ OK"
        print(f"{name:<15} {vif:<10.2f} {r2_j:<10.4f} {status}")
    
    return vif_values


# ===== DEMO =====
np.random.seed(42)
m = 200

# Create features with known multicollinearity
x1 = np.random.normal(0, 1, m)
x2 = 0.9 * x1 + 0.1 * np.random.normal(0, 1, m)  # highly correlated with x1
x3 = np.random.normal(0, 1, m)                      # independent
x4 = 0.5 * x1 + 0.5 * x3 + np.random.normal(0, 0.3, m)  # moderate correlation

X = np.column_stack([x1, x2, x3, x4])
feature_names = ["x1 (original)", "x2 (≈0.9·x1)", "x3 (independent)", "x4 (mixed)"]

vif = compute_vif(X, feature_names)

import numpy as np

def compute_aic_bic(y, y_pred, k):
    """Compute AIC and BIC for a regression model.
    
    Args:
        y: actual values
        y_pred: predicted values
        k: number of parameters (including intercept)
    
    Returns:
        (AIC, BIC)
    """
    m = len(y)
    rss = np.sum((y - y_pred) ** 2)
    
    # AIC = m * ln(RSS/m) + 2k
    aic = m * np.log(rss / m) + 2 * k
    
    # BIC = m * ln(RSS/m) + k * ln(m)
    bic = m * np.log(rss / m) + k * np.log(m)
    
    return aic, bic


def forward_selection(X, y, feature_names=None, criterion='aic'):
    """Forward feature selection using AIC or BIC.
    
    Starts with no features, greedily adds the best feature at each step.
    """
    m, n = X.shape
    if feature_names is None:
        feature_names = [f"x{i}" for i in range(n)]
    
    selected = []
    remaining = list(range(n))
    
    print("FORWARD FEATURE SELECTION")
    print("=" * 60)
    
    # Baseline: intercept-only model
    y_mean = np.full(m, np.mean(y))
    best_aic, best_bic = compute_aic_bic(y, y_mean, k=1)
    best_score = best_aic if criterion == 'aic' else best_bic
    
    print(f"Step 0: Intercept only → {criterion.upper()} = {best_score:.2f}")
    
    step = 1
    while remaining:
        scores = {}
        for j in remaining:
            # Try adding feature j
            trial = selected + [j]
            X_trial = X[:, trial]
            X_b = np.column_stack([np.ones(m), X_trial])
            w = np.linalg.solve(X_b.T @ X_b, X_b.T @ y)
            y_pred = X_b @ w
            aic, bic = compute_aic_bic(y, y_pred, k=len(trial)+1)
            scores[j] = aic if criterion == 'aic' else bic
        
        # Best candidate
        best_j = min(scores, key=scores.get)
        if scores[best_j] < best_score:
            best_score = scores[best_j]
            selected.append(best_j)
            remaining.remove(best_j)
            print(f"Step {step}: Add '{feature_names[best_j]}' "
                  f"→ {criterion.upper()} = {best_score:.2f}")
            step += 1
        else:
            print(f"\nStopping: no feature improves {criterion.upper()}")
            break
    
    print(f"\nSelected features: {[feature_names[i] for i in selected]}")
    return selected


# ===== DEMO =====
np.random.seed(42)
m = 150
x1 = np.random.normal(0, 1, m)   # useful
x2 = np.random.normal(0, 1, m)   # useful
x3 = np.random.normal(0, 1, m)   # noise
x4 = np.random.normal(0, 1, m)   # noise
X = np.column_stack([x1, x2, x3, x4])
y = 3*x1 + 5*x2 + np.random.normal(0, 2, m)  # only x1, x2 matter

names = ["area", "rooms", "noise_1", "noise_2"]
selected = forward_selection(X, y, names, criterion='bic')

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

# ===== Generate Data =====
np.random.seed(42)
m = 500
X = np.column_stack([
    np.random.uniform(5, 25, m),        # area (100 sq ft)
    np.random.randint(1, 6, m),          # rooms
    np.random.uniform(0.5, 15, m),       # distance to metro (km)
    np.random.uniform(1, 30, m),         # age of building (years)
    np.random.uniform(50, 100, m),       # safety score
])
feature_names = ['area', 'rooms', 'metro_dist', 'building_age', 'safety']

# True model: price = 5 + 3*area + 7*rooms - 2*metro + noise
y = 5 + 3*X[:,0] + 7*X[:,1] - 2*X[:,2] - 0.5*X[:,3] + 0.3*X[:,4]
y += np.random.normal(0, 5, m)

# Train-test split
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Scale features
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Fit model
model = LinearRegression()
model.fit(X_train_scaled, y_train)

# Results
y_pred = model.predict(X_test_scaled)
print("SKLEARN MULTIPLE REGRESSION RESULTS")
print("=" * 50)
print(f"Intercept: {model.intercept_:.4f}")
for name, coef in zip(feature_names, model.coef_):
    print(f"  {name:>15}: {coef:>8.4f}")
print(f"\nTrain R²: {model.score(X_train_scaled, y_train):.6f}")
print(f"Test R²:  {r2_score(y_test, y_pred):.6f}")
print(f"RMSE:     {np.sqrt(mean_squared_error(y_test, y_pred)):.4f}")
print(f"MAE:      {mean_absolute_error(y_test, y_pred):.4f}")

# Cross-validation
cv_scores = cross_val_score(model, X_train_scaled, y_train, cv=5, scoring='r2')
print(f"\n5-Fold CV R²: {cv_scores.mean():.6f} ± {cv_scores.std():.6f}")

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

# ===== Generate Nonlinear Data =====
np.random.seed(42)
x = np.sort(np.random.uniform(-3, 3, 80))
y = 0.5 * x**3 - 2 * x**2 + x + 3 + np.random.normal(0, 3, 80)
X = x.reshape(-1, 1)

# ===== Compare Different Polynomial Degrees =====
print("POLYNOMIAL DEGREE COMPARISON (5-Fold CV)")
print("=" * 60)
print(f"{'Degree':<10} {'Num Features':<15} {'CV R² Mean':<15} {'CV R² Std'}")
print("-" * 60)

results = {}
for degree in range(1, 11):
    pipeline = Pipeline([
        ('poly', PolynomialFeatures(degree=degree, include_bias=False)),
        ('scaler', StandardScaler()),
        ('reg', LinearRegression())
    ])
    
    cv_scores = cross_val_score(pipeline, X, y, cv=5, scoring='r2')
    n_features = PolynomialFeatures(degree=degree, include_bias=False).fit_transform(X).shape[1]
    
    results[degree] = {
        'mean': cv_scores.mean(),
        'std': cv_scores.std(),
        'n_features': n_features
    }
    
    marker = " ◄── BEST" if degree == 3 else ""
    print(f"  {degree:<10} {n_features:<15} {cv_scores.mean():<15.6f} {cv_scores.std():.6f}{marker}")

# ===== Best Model: Degree 3 =====
best_pipe = Pipeline([
    ('poly', PolynomialFeatures(degree=3, include_bias=False)),
    ('scaler', StandardScaler()),
    ('reg', LinearRegression())
])
best_pipe.fit(X, y)

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

for ax, deg, title in zip(axes, [1, 3, 9], 
                           ['Underfit', 'Good Fit', 'Overfit']):
    pipe = Pipeline([
        ('poly', PolynomialFeatures(degree=deg, include_bias=False)),
        ('scaler', StandardScaler()),
        ('reg', LinearRegression())
    ])
    pipe.fit(X, y)
    y_plot = pipe.predict(x_plot)
    
    ax.scatter(x, y, c='#059669', alpha=0.6, s=25, label='Data')
    ax.plot(x_plot, y_plot, c='#0891b2', linewidth=2, label=f'Degree {deg}')
    ax.set_title(f'{title} (degree={deg})')
    ax.set_ylim(-30, 40)
    ax.legend()

plt.tight_layout()
plt.savefig('poly_comparison.png', dpi=150)
plt.show()


# ===== Feature Interactions =====
print("\n\nFEATURE INTERACTION EXAMPLE")
print("=" * 60)
# 2D features with interaction
X_2d = np.column_stack([
    np.random.uniform(0, 10, 200),
    np.random.uniform(0, 10, 200)
])
y_2d = 2*X_2d[:,0] + 3*X_2d[:,1] + 0.5*X_2d[:,0]*X_2d[:,1] + np.random.normal(0,2,200)

poly2d = PolynomialFeatures(degree=2, include_bias=False, interaction_only=False)
X_poly = poly2d.fit_transform(X_2d)
print(f"Original features: {X_2d.shape[1]}")
print(f"Poly features (degree=2): {X_poly.shape[1]}")
print(f"Feature names: {poly2d.get_feature_names_out(['x1','x2'])}")

model_2d = LinearRegression().fit(X_poly, y_2d)
print(f"\nCoefficients:")
for name, coef in zip(poly2d.get_feature_names_out(['x1','x2']), model_2d.coef_):
    print(f"  {name:>10}: {coef:.4f}")
print(f"  Intercept: {model_2d.intercept_:.4f}")

import numpy as np
import pandas as pd
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.stattools import durbin_watson
import statsmodels.api as sm
from scipy import stats

# ===== Create Dataset =====
np.random.seed(42)
m = 300

data = pd.DataFrame({
    'area': np.random.uniform(5, 25, m),
    'rooms': np.random.randint(1, 6, m),
    'bathrooms': np.random.randint(1, 4, m),
    'metro_km': np.random.uniform(0.5, 15, m),
    'age_years': np.random.uniform(1, 30, m),
})
# Add correlated feature (total_rooms ≈ rooms + bathrooms)
data['total_rooms'] = data['rooms'] + data['bathrooms'] + np.random.normal(0, 0.5, m)

# Target
y = (5 + 3*data['area'] + 7*data['rooms'] + 5*data['bathrooms']
     - 2*data['metro_km'] - 0.5*data['age_years']
     + np.random.normal(0, 5, m))

# ===== VIF Calculation =====
X = data.values
print("VIF ANALYSIS")
print("=" * 50)
for i, name in enumerate(data.columns):
    vif = variance_inflation_factor(X, i)
    flag = " ← ⚠️ HIGH!" if vif > 5 else ""
    print(f"  {name:>15}: VIF = {vif:.2f}{flag}")

# ===== Full OLS Report with statsmodels =====
X_ols = sm.add_constant(data)
ols_model = sm.OLS(y, X_ols).fit()
print("\n" + "=" * 50)
print("OLS REGRESSION SUMMARY (via statsmodels)")
print("=" * 50)
print(ols_model.summary())

# ===== Durbin-Watson Test =====
dw = durbin_watson(ols_model.resid)
print(f"\nDurbin-Watson statistic: {dw:.4f}")
if 1.5 < dw < 2.5:
    print("  → No significant autocorrelation ✅")
else:
    print("  → Autocorrelation detected ⚠️")

# ===== Q-Q Plot of Residuals =====
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# Residual vs Fitted
axes[0].scatter(ols_model.fittedvalues, ols_model.resid, alpha=0.5, c='#059669', s=15)
axes[0].axhline(y=0, color='red', linestyle='--')
axes[0].set_xlabel('Fitted Values')
axes[0].set_ylabel('Residuals')
axes[0].set_title('Residuals vs Fitted')

# Q-Q plot
stats.probplot(ols_model.resid, dist="norm", plot=axes[1])
axes[1].set_title('Q-Q Plot of Residuals')

# Histogram of residuals
axes[2].hist(ols_model.resid, bins=25, color='#0891b2', alpha=0.7, edgecolor='white')
axes[2].set_title('Residual Distribution')
axes[2].set_xlabel('Residual')

plt.tight_layout()
plt.savefig('diagnostics.png', dpi=150)
plt.show()

import numpy as np
import tensorflow as tf
from tensorflow import keras
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

# ===== Generate Synthetic Data (5 features) =====
np.random.seed(42)
tf.random.set_seed(42)
m = 1000

X = np.column_stack([
    np.random.uniform(5, 25, m),         # area
    np.random.randint(1, 6, m).astype(float),  # rooms
    np.random.uniform(0.5, 15, m),       # metro distance
    np.random.uniform(1, 30, m),         # building age
    np.random.uniform(50, 100, m),       # safety score
])

# Non-linear relationship with interactions
y = (5 + 3*X[:,0] + 7*X[:,1] - 2*X[:,2] - 0.5*X[:,3] + 0.3*X[:,4]
     + 0.1*X[:,0]*X[:,1]    # interaction term
     - 0.05*X[:,0]**2       # quadratic term
     + np.random.normal(0, 5, m))

# Split and scale
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 1: Simple Linear (no hidden layers) =====
model_linear = keras.Sequential([
    keras.layers.Dense(1, input_shape=(5,), activation='linear',
                       kernel_initializer='normal')
], name='Linear_Model')

model_linear.compile(optimizer='adam', loss='mse', metrics=['mae'])
model_linear.fit(X_train_s, y_train, epochs=100, batch_size=32, 
                 verbose=0, validation_split=0.1)

# ===== Model 2: Deep Regression (captures nonlinearity) =====
model_deep = keras.Sequential([
    keras.layers.Dense(64, input_shape=(5,), activation='relu'),
    keras.layers.BatchNormalization(),
    keras.layers.Dropout(0.2),
    keras.layers.Dense(32, activation='relu'),
    keras.layers.BatchNormalization(),
    keras.layers.Dropout(0.1),
    keras.layers.Dense(16, activation='relu'),
    keras.layers.Dense(1, activation='linear')
], name='Deep_Regression')

model_deep.compile(
    optimizer=keras.optimizers.Adam(learning_rate=0.001),
    loss='mse',
    metrics=['mae']
)

# Train with early stopping
early_stop = keras.callbacks.EarlyStopping(
    monitor='val_loss', patience=15, restore_best_weights=True
)

history = model_deep.fit(
    X_train_s, y_train,
    epochs=200, batch_size=32,
    validation_split=0.15,
    callbacks=[early_stop],
    verbose=0
)

# ===== Evaluate Both =====
from sklearn.metrics import r2_score, mean_squared_error

for name, model in [('Linear', model_linear), ('Deep', model_deep)]:
    y_pred = model.predict(X_test_s, verbose=0).flatten()
    r2 = r2_score(y_test, y_pred)
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))
    print(f"{name:>8} Model → R²: {r2:.6f}, RMSE: {rmse:.4f}")

# ===== Model 3: Functional API (Multi-Input) =====
# For more complex architectures
input_continuous = keras.Input(shape=(3,), name='continuous')  # area, metro, age
input_discrete = keras.Input(shape=(2,), name='discrete')     # rooms, safety

# Process branches separately
cont_branch = keras.layers.Dense(32, activation='relu')(input_continuous)
cont_branch = keras.layers.Dense(16, activation='relu')(cont_branch)

disc_branch = keras.layers.Dense(16, activation='relu')(input_discrete)

# Merge
merged = keras.layers.Concatenate()([cont_branch, disc_branch])
merged = keras.layers.Dense(16, activation='relu')(merged)
output = keras.layers.Dense(1, activation='linear')(merged)

model_multi = keras.Model(
    inputs=[input_continuous, input_discrete],
    outputs=output,
    name='Multi_Input_Model'
)
model_multi.compile(optimizer='adam', loss='mse')

# Prepare multi-input data
X_train_cont = X_train_s[:, [0, 2, 3]]  # area, metro, age
X_train_disc = X_train_s[:, [1, 4]]     # rooms, safety
X_test_cont = X_test_s[:, [0, 2, 3]]
X_test_disc = X_test_s[:, [1, 4]]

model_multi.fit(
    [X_train_cont, X_train_disc], y_train,
    epochs=100, batch_size=32, verbose=0, validation_split=0.1
)

y_pred_multi = model_multi.predict([X_test_cont, X_test_disc], verbose=0).flatten()
print(f"{'Multi':>8} Model → R²: {r2_score(y_test, y_pred_multi):.6f}")

model_multi.summary()