# Complex demo: Load California Housing, clean, compute descriptive stats, feature engineering, model with LinearRegression, evaluate, visualize distributions/correlations
import pandas as pd
import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
import seaborn as sns

def advanced_housing_exploration_demo():
    # Load data
    housing = fetch_california_housing()
    df = pd.DataFrame(housing.data, columns=housing.feature_names)
    df['Target'] = housing.target
    
    # Data cleaning
    df = df.dropna()
    df = df[df.duplicated() == False]
    
    # Descriptive stats
    stats = df.describe()
    correlation = df.corr()
    print("Descriptive Statistics:\n", stats)
    print("Correlation Matrix:\n", correlation)
    
    # Feature engineering
    df['RoomsPerPop'] = df['AveRooms'] / (df['Population'] + 1e-5)
    
    # Preprocessing
    scaler = StandardScaler()
    X = df.drop('Target', axis=1)
    y = df['Target']
    X_scaled = scaler.fit_transform(X)
    X_train, X_test, y_train, y_test = 
        train_test_split(X_scaled, y, test_size=0.25, random_state=42)
    
    # Model and evaluate
    model = LinearRegression()
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    mse = mean_squared_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    print(f"MSE: {mse:.2f}, R2: {r2:.2f}")
    
    # Visualize distributions
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    axes = axes.ravel()
    
    for i, col in enumerate(df.columns[:6]):
        sns.histplot(df[col], kde=True, ax=axes[i])
        axes[i].set_title(f'Distribution of {col}')
        axes[i].annotate(f'Mean: {df[col].mean():.2f}', 
            xy=(0.05, 0.95), xycoords='axes fraction', fontsize=10)
    plt.tight_layout()
    plt.show()
    
    # Visualize correlations
    plt.figure(figsize=(10, 8))
    sns.heatmap(correlation, annot=True, cmap='coolwarm')
    plt.title('Correlation Heatmap for Housing Data')
    plt.annotate('High Corr', xy=(0, 0), xytext=(1, 1), arrowprops=dict(facecolor='black', shrink=0.05))
    plt.show()
    
    # Visualize predictions
    plt.figure(figsize=(8, 6))
    plt.scatter(y_test, y_pred, c='blue', alpha=0.5)
    plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'k--', lw=2)
    plt.title('Housing Regression Predictions')
    plt.xlabel('True Values')
    plt.ylabel('Predictions')
    plt.annotate('Best Fit', xy=(y_test.mean(), y_pred.mean()), xytext=(y_test.mean() + 0.5, y_pred.mean() + 0.5), arrowprops=dict(facecolor='black', shrink=0.05))
    plt.show()

advanced_housing_exploration_demo()