Linear Regression &
Gradient Descent

import numpy as np
import matplotlib.pyplot as plt

class LinearRegressionScratch:
    """
    Linear Regression from scratch.
    Supports: Normal Equation, Batch GD, SGD, Mini-batch GD.
    """

    def __init__(self, method='normal', lr=0.01, n_iters=1000,
                 batch_size=32, tol=1e-6, verbose=False):
        """
        Parameters:
        -----------
        method : str, one of 'normal', 'batch_gd', 'sgd', 'mini_batch_gd'
        lr     : float, learning rate (alpha)
        n_iters: int, max iterations for GD methods
        batch_size: int, for mini-batch GD
        tol    : float, convergence tolerance
        verbose: bool, print progress
        """
        self.method = method
        self.lr = lr
        self.n_iters = n_iters
        self.batch_size = batch_size
        self.tol = tol
        self.verbose = verbose
        self.weights = None       # includes bias
        self.cost_history = []

    def _add_bias(self, X):
        """Add a column of 1s for the bias term."""
        return np.c_[np.ones(X.shape[0]), X]

    def _compute_cost(self, X, y):
        """Compute MSE cost: J = (1/n)||y - Xw||^2"""
        n = len(y)
        predictions = X @ self.weights
        cost = (1 / n) * np.sum((y - predictions) ** 2)
        return cost

    def _compute_gradient(self, X, y):
        """Compute gradient: (2/n) * X^T(Xw - y)"""
        n = len(y)
        predictions = X @ self.weights
        gradient = (2 / n) * X.T @ (predictions - y)
        return gradient

    def fit_normal_equation(self, X, y):
        """Closed-form solution: w = (X^TX)^{-1} X^Ty"""
        X_b = self._add_bias(X)
        # Use np.linalg.pinv for numerical stability
        self.weights = np.linalg.pinv(X_b.T @ X_b) @ X_b.T @ y
        if self.verbose:
            cost = self._compute_cost(X_b, y)
            print(f"Normal Equation — MSE: {cost:.6f}")

    def fit_batch_gd(self, X, y):
        """Full-batch Gradient Descent."""
        X_b = self._add_bias(X)
        n, p = X_b.shape
        self.weights = np.zeros(p)   # initialize to zeros
        self.cost_history = []

        for i in range(self.n_iters):
            gradient = self._compute_gradient(X_b, y)
            self.weights -= self.lr * gradient
            cost = self._compute_cost(X_b, y)
            self.cost_history.append(cost)

            if self.verbose and i % 100 == 0:
                print(f"Iter {i:4d} | MSE: {cost:.6f}")

            # Convergence check
            if len(self.cost_history) > 1:
                if abs(self.cost_history[-1] - self.cost_history[-2]) < self.tol:
                    if self.verbose:
                        print(f"Converged at iteration {i}")
                    break

    def fit_sgd(self, X, y):
        """Stochastic Gradient Descent (one sample at a time)."""
        X_b = self._add_bias(X)
        n, p = X_b.shape
        self.weights = np.zeros(p)
        self.cost_history = []

        for epoch in range(self.n_iters):
            indices = np.random.permutation(n)
            for idx in indices:
                xi = X_b[idx:idx+1]
                yi = y[idx:idx+1]
                gradient = (2) * xi.T @ (xi @ self.weights - yi)
                self.weights -= self.lr * gradient

            cost = self._compute_cost(X_b, y)
            self.cost_history.append(cost)
            if self.verbose and epoch % 100 == 0:
                print(f"Epoch {epoch:4d} | MSE: {cost:.6f}")

    def fit_mini_batch_gd(self, X, y):
        """Mini-batch Gradient Descent."""
        X_b = self._add_bias(X)
        n, p = X_b.shape
        self.weights = np.zeros(p)
        self.cost_history = []

        for epoch in range(self.n_iters):
            indices = np.random.permutation(n)
            for start in range(0, n, self.batch_size):
                end = min(start + self.batch_size, n)
                batch_idx = indices[start:end]
                X_batch = X_b[batch_idx]
                y_batch = y[batch_idx]
                gradient = (2 / len(y_batch)) * X_batch.T @ (X_batch @ self.weights - y_batch)
                self.weights -= self.lr * gradient

            cost = self._compute_cost(X_b, y)
            self.cost_history.append(cost)

    def fit(self, X, y):
        """Fit model using the chosen method."""
        X = np.array(X, dtype=np.float64)
        y = np.array(y, dtype=np.float64)

        if self.method == 'normal':
            self.fit_normal_equation(X, y)
        elif self.method == 'batch_gd':
            self.fit_batch_gd(X, y)
        elif self.method == 'sgd':
            self.fit_sgd(X, y)
        elif self.method == 'mini_batch_gd':
            self.fit_mini_batch_gd(X, y)
        else:
            raise ValueError(f"Unknown method: {self.method}")
        return self

    def predict(self, X):
        """Make predictions."""
        X = np.array(X, dtype=np.float64)
        X_b = self._add_bias(X)
        return X_b @ self.weights

    def score(self, X, y):
        """Compute R² score."""
        y = np.array(y, dtype=np.float64)
        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)

    def get_metrics(self, X, y):
        """Return dict of all metrics."""
        y = np.array(y, dtype=np.float64)
        y_pred = self.predict(X)
        n = len(y)
        p = len(self.weights) - 1  # exclude bias

        mse = np.mean((y - y_pred) ** 2)
        rmse = np.sqrt(mse)
        mae = np.mean(np.abs(y - y_pred))
        r2 = self.score(X, y)
        adj_r2 = 1 - ((1 - r2) * (n - 1) / (n - p - 1))

        return {
            'MSE': mse, 'RMSE': rmse, 'MAE': mae,
            'R²': r2, 'Adj R²': adj_r2
        }

    def plot_cost_history(self):
        """Plot cost vs iterations (for GD methods)."""
        if not self.cost_history:
            print("No cost history (use a GD method).")
            return
        plt.figure(figsize=(10, 5))
        plt.plot(self.cost_history, color='#059669', linewidth=2)
        plt.xlabel('Iteration')
        plt.ylabel('MSE Cost')
        plt.title('Cost Function Convergence')
        plt.grid(True, alpha=0.3)
        plt.show()

    def plot_fit(self, X, y, feature_idx=0):
        """Plot data + regression line (for 1D visualization)."""
        X = np.array(X, dtype=np.float64)
        y = np.array(y, dtype=np.float64)
        y_pred = self.predict(X)

        plt.figure(figsize=(10, 6))
        plt.scatter(X[:, feature_idx], y, color='#0891b2',
                    s=80, alpha=0.7, label='Actual', zorder=5)

        # Sort for smooth line
        sort_idx = X[:, feature_idx].argsort()
        plt.plot(X[sort_idx, feature_idx], y_pred[sort_idx],
                 color='#059669', linewidth=2.5, label='Predicted')

        plt.xlabel(f'Feature {feature_idx}')
        plt.ylabel('Target')
        plt.title('Linear Regression Fit')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.show()

# === Example: Bangalore flat price prediction ===
import numpy as np

# Data: [area_sqft, bedrooms, age_years]
X = np.array([
    [600,  1, 10],
    [800,  2,  7],
    [1000, 2,  5],
    [1200, 3,  3],
    [1500, 3,  1],
    [1800, 4,  0],
    [900,  2,  8],
    [1100, 3,  4],
])
y = np.array([30, 45, 55, 68, 82, 105, 40, 62])  # Price in ₹ Lakhs

# --- Method 1: Normal Equation ---
model_ne = LinearRegressionScratch(method='normal', verbose=True)
model_ne.fit(X, y)
print("Weights (Normal Eq):", model_ne.weights)
print("Metrics:", model_ne.get_metrics(X, y))

# --- Method 2: Batch Gradient Descent ---
# Feature scaling is CRITICAL for GD!
X_scaled = (X - X.mean(axis=0)) / X.std(axis=0)

model_gd = LinearRegressionScratch(
    method='batch_gd', lr=0.01, n_iters=1000, verbose=True
)
model_gd.fit(X_scaled, y)
print("R² (Batch GD):", model_gd.score(X_scaled, y))
model_gd.plot_cost_history()

# --- Predict new flat ---
new_flat = np.array([[1300, 3, 2]])
price = model_ne.predict(new_flat)
print(f"Predicted price for 1300 sqft, 3BHK, 2yr old: ₹{price[0]:.1f} Lakhs")

class StandardScaler:
    """Standardize features: z = (x - mean) / std"""

    def fit(self, X):
        self.mean_ = np.mean(X, axis=0)
        self.std_ = np.std(X, axis=0)
        self.std_[self.std_ == 0] = 1  # avoid division by zero
        return self

    def transform(self, X):
        return (X - self.mean_) / self.std_

    def fit_transform(self, X):
        return self.fit(X).transform(X)

# Why is scaling critical for GD?
# Without scaling: area~1000, bedrooms~3, age~5
# The gradient for 'area' is ~300x larger than for 'bedrooms'
# GD oscillates wildly on the area axis while barely moving on bedrooms
# With scaling: all features have similar ranges → smooth convergence

def diagnostic_plots(y_true, y_pred, feature_name="Feature"):
    """Generate diagnostic plots for regression analysis."""
    residuals = y_true - y_pred
    fig, axes = plt.subplots(1, 3, figsize=(18, 5))

    # 1. Residual Plot
    axes[0].scatter(y_pred, residuals, color='#0891b2', alpha=0.6, s=50)
    axes[0].axhline(y=0, color='#dc2626', linestyle='--', linewidth=1.5)
    axes[0].set_xlabel('Predicted Values')
    axes[0].set_ylabel('Residuals')
    axes[0].set_title('Residual Plot')
    axes[0].grid(True, alpha=0.3)

    # 2. Q-Q Plot (manual implementation)
    sorted_residuals = np.sort(residuals)
    n = len(sorted_residuals)
    theoretical = np.array([
        np.sqrt(2) * 0.4769 *
        ((2 * (i + 1) - 1) / (2 * n) - 0.5)
        for i in range(n)
    ])  # Simplified; in practice, use scipy.stats.norm.ppf
    axes[1].scatter(theoretical, sorted_residuals,
                    color='#059669', alpha=0.7, s=50)
    # Reference line
    lims = [min(theoretical.min(), sorted_residuals.min()),
            max(theoretical.max(), sorted_residuals.max())]
    axes[1].plot(lims, lims, 'r--', linewidth=1.5)
    axes[1].set_xlabel('Theoretical Quantiles')
    axes[1].set_ylabel('Sample Quantiles')
    axes[1].set_title('Q-Q Plot')
    axes[1].grid(True, alpha=0.3)

    # 3. Actual vs Predicted
    axes[2].scatter(y_true, y_pred, color='#8b5cf6', alpha=0.6, s=50)
    lims = [min(y_true.min(), y_pred.min()),
            max(y_true.max(), y_pred.max())]
    axes[2].plot(lims, lims, 'r--', linewidth=1.5)
    axes[2].set_xlabel('Actual Values')
    axes[2].set_ylabel('Predicted Values')
    axes[2].set_title('Actual vs Predicted')
    axes[2].grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

# Usage:
# y_pred = model.predict(X)
# diagnostic_plots(y, y_pred)

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

# === Generate sample Indian housing data ===
np.random.seed(42)
n_samples = 500

area = np.random.uniform(500, 2500, n_samples)
bedrooms = np.random.randint(1, 5, n_samples)
age = np.random.uniform(0, 20, n_samples)
floor = np.random.randint(1, 20, n_samples)

# True relationship + noise
price = 0.05 * area + 8 * bedrooms - 1.5 * age + 0.3 * floor + 5 + \
        np.random.normal(0, 5, n_samples)

X = np.column_stack([area, bedrooms, age, floor])
y = price

# Split & 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)

# === Build TensorFlow Model ===
model = tf.keras.Sequential([
    tf.keras.layers.Dense(1, input_shape=(4,),
                          kernel_initializer='normal',
                          name='linear_layer')
])

# Linear regression = 1 Dense layer with no activation!
# This is equivalent to ŷ = Xw + b

model.compile(
    optimizer=tf.keras.optimizers.SGD(learning_rate=0.01),
    loss='mse',               # Mean Squared Error
    metrics=['mae']            # Also track MAE
)

model.summary()

# === Train ===
history = model.fit(
    X_train_s, y_train,
    epochs=100,
    batch_size=32,
    validation_split=0.15,
    verbose=1
)

# === Evaluate ===
test_loss, test_mae = model.evaluate(X_test_s, y_test)
print(f"\nTest MSE: {test_loss:.4f}")
print(f"Test MAE: {test_mae:.4f}")

# === Extract learned weights ===
weights, bias = model.get_layer('linear_layer').get_weights()
print(f"\nLearned weights: {weights.flatten()}")
print(f"Learned bias: {bias[0]:.4f}")

# === Plot training history ===
import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

ax1.plot(history.history['loss'], label='Train MSE', color='#059669')
ax1.plot(history.history['val_loss'], label='Val MSE', color='#0891b2')
ax1.set_xlabel('Epoch'); ax1.set_ylabel('MSE')
ax1.set_title('Loss Convergence'); ax1.legend(); ax1.grid(True, alpha=0.3)

ax2.plot(history.history['mae'], label='Train MAE', color='#059669')
ax2.plot(history.history['val_mae'], label='Val MAE', color='#0891b2')
ax2.set_xlabel('Epoch'); ax2.set_ylabel('MAE')
ax2.set_title('MAE Over Training'); ax2.legend(); ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# === Predict new flat ===
new_flat = np.array([[1300, 3, 2, 10]])
new_flat_s = scaler.transform(new_flat)
pred_price = model.predict(new_flat_s)
print(f"Predicted price: ₹{pred_price[0][0]:.1f} Lakhs")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, SGDRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.pipeline import Pipeline
from sklearn.metrics import (
    mean_squared_error, mean_absolute_error, r2_score
)

# ============================================
# 1. CREATE SYNTHETIC INDIAN HOUSING DATASET
# ============================================
np.random.seed(42)
n = 1000

data = pd.DataFrame({
    'area_sqft':  np.random.uniform(400, 3000, n),
    'bedrooms':   np.random.randint(1, 6, n),
    'bathrooms':  np.random.randint(1, 4, n),
    'age_years':  np.random.uniform(0, 25, n),
    'floor':      np.random.randint(1, 25, n),
    'city_tier':  np.random.choice([1, 2, 3], n, p=[0.3, 0.4, 0.3]),
})

# Generate realistic prices (₹ Lakhs)
city_multiplier = {1: 1.5, 2: 1.0, 3: 0.6}
data['price_lakhs'] = (
    0.045 * data['area_sqft'] +
    7 * data['bedrooms'] +
    5 * data['bathrooms'] -
    1.2 * data['age_years'] +
    0.3 * data['floor'] +
    10 * data['city_tier'].map(city_multiplier) +
    np.random.normal(0, 8, n)
)

print(data.head())
print(data.describe())

# ============================================
# 2. EXPLORATORY DATA ANALYSIS
# ============================================
print("\nCorrelation with price:")
print(data.corr()['price_lakhs'].sort_values(ascending=False))

# ============================================
# 3. PREPARE FEATURES & TARGET
# ============================================
features = ['area_sqft', 'bedrooms', 'bathrooms', 'age_years', 'floor']
X = data[features].values
y = data['price_lakhs'].values

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

# ============================================
# 4. MODEL 1: BASIC LINEAR REGRESSION
# ============================================
model = LinearRegression()
model.fit(X_train, y_train)

print("\n=== Linear Regression Results ===")
print(f"Intercept: {model.intercept_:.4f}")
for feat, coef in zip(features, model.coef_):
    print(f"  {feat:15s}: {coef:+.6f}")

y_pred = model.predict(X_test)
print(f"\nTest R²:   {r2_score(y_test, y_pred):.4f}")
print(f"Test RMSE: {np.sqrt(mean_squared_error(y_test, y_pred)):.4f}")
print(f"Test MAE:  {mean_absolute_error(y_test, y_pred):.4f}")

# ============================================
# 5. CROSS-VALIDATION
# ============================================
cv_scores = cross_val_score(model, X, y, cv=5, scoring='r2')
print(f"\n5-Fold CV R² scores: {cv_scores}")
print(f"Mean CV R²: {cv_scores.mean():.4f} ± {cv_scores.std():.4f}")

# ============================================
# 6. MODEL 2: PIPELINE WITH SCALING + SGD
# ============================================
sgd_pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('sgd', SGDRegressor(max_iter=1000, tol=1e-4,
                         learning_rate='optimal',
                         random_state=42))
])
sgd_pipeline.fit(X_train, y_train)
y_pred_sgd = sgd_pipeline.predict(X_test)
print(f"\nSGD R²:   {r2_score(y_test, y_pred_sgd):.4f}")

# ============================================
# 7. MODEL 3: POLYNOMIAL REGRESSION
# ============================================
poly_pipeline = Pipeline([
    ('poly', PolynomialFeatures(degree=2, include_bias=False)),
    ('scaler', StandardScaler()),
    ('reg', LinearRegression())
])
poly_pipeline.fit(X_train, y_train)
y_pred_poly = poly_pipeline.predict(X_test)
print(f"Poly(2) R²: {r2_score(y_test, y_pred_poly):.4f}")

# ============================================
# 8. VISUALIZATION
# ============================================
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Actual vs Predicted
axes[0,0].scatter(y_test, y_pred, alpha=0.5, color='#059669', s=30)
axes[0,0].plot([y_test.min(), y_test.max()],
              [y_test.min(), y_test.max()], 'r--')
axes[0,0].set_title('Actual vs Predicted')

# Residuals
residuals = y_test - y_pred
axes[0,1].scatter(y_pred, residuals, alpha=0.5, color='#0891b2', s=30)
axes[0,1].axhline(y=0, color='red', linestyle='--')
axes[0,1].set_title('Residual Plot')

# Residual distribution
axes[1,0].hist(residuals, bins=30, color='#059669', alpha=0.7, edgecolor='white')
axes[1,0].set_title('Residual Distribution')

# Feature importance (coefficients)
coefs = pd.Series(model.coef_, index=features).sort_values()
coefs.plot.barh(ax=axes[1,1], color='#0891b2')
axes[1,1].set_title('Feature Coefficients')

plt.tight_layout()
plt.show()