# Experiment 3: Multiple Linear Regression (House Price Prediction)
# Requirements: pip install numpy pandas matplotlib seaborn scikit-learn statsmodels

# Step 1: Environment Setup & Synthetic Data Generation
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.api as sm

# Set random seed for reproducibility
np.random.seed(42)
n_samples = 500

# Generate synthetic features
sqft = np.random.normal(2100, 500, n_samples).clip(800, 4500)
bedrooms = np.random.choice([1, 2, 3, 4, 5], size=n_samples, p=[0.1, 0.25, 0.4, 0.2, 0.05])
dist_transit = np.random.exponential(scale=4.0, size=n_samples).clip(0.5, 20.0)
age = np.random.uniform(1, 40, n_samples)
school_rating = np.random.uniform(3.0, 9.8, n_samples)

# True linear relationship + Gaussian noise
noise = np.random.normal(0, 25000, n_samples)
price = (150 * sqft) + (12000 * bedrooms) - (8500 * dist_transit) - (1500 * age) + (9000 * school_rating) + 50000 + noise

# Construct DataFrame
df = pd.DataFrame({
    'Square_Footage': sqft,
    'Number_of_Bedrooms': bedrooms,
    'Distance_to_Transit_km': dist_transit,
    'Age_of_Property': age,
    'School_Rating': school_rating,
    'House_Price': price
})

print("Dataset Head:")
print(df.head())

# Step 2: Train-Test Split (80:20)
# Define Independent Features (X) and Target Variable (Y)
X = df[['Square_Footage', 'Number_of_Bedrooms', 'Distance_to_Transit_km', 'Age_of_Property', 'School_Rating']]
Y = df['House_Price']

# Split data into 80% training and 20% testing sets
X_train, X_test, Y_train, Y_test = train_test_split(
    X, Y, test_size=0.20, random_state=42
)

print(f"Training Set Shape: {X_train.shape}")
print(f"Testing Set Shape:  {X_test.shape}")

# Step 3: Model Fitting & Parameter Estimation
# Instantiate and fit Linear Regression model
model = LinearRegression()
model.fit(X_train, Y_train)

# Extract intercept and slope coefficients
intercept = model.intercept_
coefficients = pd.DataFrame({
    'Feature': X.columns,
    'Coefficient (Beta)': model.coef_
})

print(f"Intercept (Beta_0): ${intercept:,.2f}\n")
print("Learned Feature Coefficients:")
print(coefficients.to_string(index=False))

# Step 4: Predictions & Quantitative Performance Evaluation
# Model Predictions on Test Data
Y_pred = model.predict(X_test)

# Calculate Evaluation Metrics
mae = mean_absolute_error(Y_test, Y_pred)
mse = mean_squared_error(Y_test, Y_pred)
rmse = np.sqrt(mse)
r2 = r2_score(Y_test, Y_pred)

# Adjusted R-Squared
n = len(Y_test)
p = X_test.shape[1]
adj_r2 = 1 - (1 - r2) * (n - 1) / (n - p - 1)

print("--- MODEL PERFORMANCE METRICS ---")
print(f"Mean Absolute Error (MAE)      : ${mae:,.2f}")
print(f"Mean Squared Error (MSE)        : ${mse:,.2f}")
print(f"Root Mean Squared Error (RMSE)  : ${rmse:,.2f}")
print(f"R-squared Score (R²)            : {r2:.4f}")
print(f"Adjusted R-squared Score       : {adj_r2:.4f}")

# Step 5: Residual Diagnostics & VIF Check
# Calculate Residuals
residuals = Y_test - Y_pred

# Calculate Variance Inflation Factor (VIF) for Multicollinearity
vif_data = pd.DataFrame()
vif_data["Feature"] = X.columns
vif_data["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]

print("--- VARIANCE INFLATION FACTOR (VIF) ---")
print(vif_data.to_string(index=False))

# Statsmodels OLS detailed summary
X_train_sm = sm.add_constant(X_train)
ols_model = sm.OLS(Y_train, X_train_sm).fit()
print(ols_model.summary())
