Regression Analysis Modeling

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Core Regression Analysis Engine
Comprehensive regression modeling with multiple algorithms and automated feature engineering
"""

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split, cross_val_score, GridSearchCV
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.preprocessing import StandardScaler, LabelEncoder, PolynomialFeatures
from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error, mean_absolute_percentage_error
import warnings
warnings.filterwarnings('ignore')

class RegressionAnalyzer:
    """Comprehensive regression analysis engine with multiple algorithms"""

    def __init__(self, random_state=42, chinese_font='SimHei'):
        """Initialize the regression analyzer"""
        self.random_state = random_state
        self.chinese_font = chinese_font
        self.models = {}
        self.scalers = {}
        self.feature_names = None
        self.target_name = None
        self.best_model = None
        self.best_model_name = None

    def load_and_validate_data(self, file_path, target_column, **kwargs):
        """
        Load and validate data for regression analysis

        Args:
            file_path (str): Path to the CSV file
            target_column (str): Name of the target variable column
            **kwargs: Additional parameters for pd.read_csv()

        Returns:
            tuple: (X, y) features and target
        """
        print("=== 数据加载与验证 ===")

        # Load data with different encodings
        try:
            df = pd.read_csv(file_path, encoding='utf-8', **kwargs)
        except UnicodeDecodeError:
            try:
                df = pd.read_csv(file_path, encoding='gbk', **kwargs)
            except:
                df = pd.read_csv(file_path, encoding='latin-1', **kwargs)

        print(f"数据集形状: {df.shape}")
        print(f"列名: {list(df.columns)}")

        # Check if target column exists
        if target_column not in df.columns:
            # Try to find similar column names
            possible_targets = [col for col in df.columns
                              if target_column.lower() in col.lower() or
                                 col.lower() in target_column.lower()]
            if possible_targets:
                target_column = possible_targets[0]
                print(f"自动识别目标列: {target_column}")
            else:
                raise ValueError(f"找不到目标列: {target_column}")

        # Separate features and target
        X = df.drop(columns=[target_column])
        y = df[target_column]

        self.feature_names = list(X.columns)
        self.target_name = target_column

        # Data validation
        print(f"\n数据质量检查:")
        print(f"- 缺失值: {X.isnull().sum().sum()}")
        print(f"- 目标变量缺失值: {y.isnull().sum()}")
        print(f"- 特征数量: {X.shape[1]}")
        print(f"- 样本数量: {X.shape[0]}")

        return X, y

    def preprocess_data(self, X, y, handle_missing='auto', handle_outliers='iqr'):
        """
        Preprocess data with missing value handling and outlier detection

        Args:
            X (pd.DataFrame): Features
            y (pd.Series): Target
            handle_missing (str): How to handle missing values
            handle_outliers (str): How to handle outliers

        Returns:
            tuple: (X_processed, y_processed)
        """
        print("\n=== 数据预处理 ===")

        X_processed = X.copy()
        y_processed = y.copy()

        # Handle missing values
        if handle_missing == 'auto':
            # Numerical columns: median imputation
            numerical_cols = X_processed.select_dtypes(include=[np.number]).columns
            for col in numerical_cols:
                if X_processed[col].isnull().sum() > 0:
                    median_val = X_processed[col].median()
                    X_processed[col].fillna(median_val, inplace=True)
                    print(f"- {col}: 用中位数 {median_val:.2f} 填充 {X_processed[col].isnull().sum()} 个缺失值")

            # Categorical columns: mode imputation
            categorical_cols = X_processed.select_dtypes(include=['object']).columns
            for col in categorical_cols:
                if X_processed[col].isnull().sum() > 0:
                    mode_val = X_processed[col].mode()[0] if len(X_processed[col].mode()) > 0 else 'Unknown'
                    X_processed[col].fillna(mode_val, inplace=True)
                    print(f"- {col}: 用众数 '{mode_val}' 填充 {X_processed[col].isnull().sum()} 个缺失值")

        # Remove rows with missing target
        missing_target = y_processed.isnull().sum()
        if missing_target > 0:
            mask = ~y_processed.isnull()
            X_processed = X_processed[mask]
            y_processed = y_processed[mask]
            print(f"- 移除 {missing_target} 行目标变量缺失值")

        # Handle outliers in target variable
        if handle_outliers == 'iqr' and y_processed.dtype in [np.number, 'int64', 'float64']:
            Q1 = y_processed.quantile(0.25)
            Q3 = y_processed.quantile(0.75)
            IQR = Q3 - Q1
            lower_bound = Q1 - 1.5 * IQR
            upper_bound = Q3 + 1.5 * IQR

            outlier_mask = (y_processed >= lower_bound) & (y_processed <= upper_bound)
            outliers_removed = len(y_processed) - outlier_mask.sum()

            if outliers_removed > 0:
                X_processed = X_processed[outlier_mask]
                y_processed = y_processed[outlier_mask]
                print(f"- 移除 {outliers_removed} 个异常值 (IQR方法)")

        print(f"预处理后数据形状: {X_processed.shape}")
        return X_processed, y_processed

    def encode_categorical_features(self, X):
        """
        Encode categorical features to numerical

        Args:
            X (pd.DataFrame): Features with categorical columns

        Returns:
            pd.DataFrame: Features with encoded categorical variables
        """
        print("\n=== 分类特征编码 ===")

        X_encoded = X.copy()
        categorical_cols = X_encoded.select_dtypes(include=['object']).columns

        if len(categorical_cols) > 0:
            print(f"发现 {len(categorical_cols)} 个分类特征:")

            label_encoders = {}

            for col in categorical_cols:
                # Use label encoding for high-cardinality features
                if X_encoded[col].nunique() > 10:
                    le = LabelEncoder()
                    X_encoded[col] = le.fit_transform(X_encoded[col].astype(str))
                    label_encoders[col] = le
                    print(f"- {col}: 标签编码 ({X_encoded[col].nunique()} 个类别)")
                else:
                    # Use one-hot encoding for low-cardinality features
                    dummies = pd.get_dummies(X_encoded[col], prefix=col, drop_first=True)
                    X_encoded = pd.concat([X_encoded.drop(columns=[col]), dummies], axis=1)
                    print(f"- {col}: 独热编码 ({X_encoded[col].nunique()} 个类别)")

            self.label_encoders = label_encoders

        self.feature_names = list(X_encoded.columns)
        print(f"编码后特征数量: {len(self.feature_names)}")

        return X_encoded

    def create_interaction_features(self, X):
        """
        Create interaction features for improved model performance

        Args:
            X (pd.DataFrame): Original features

        Returns:
            pd.DataFrame: Features with interactions
        """
        print("\n=== 交互特征生成 ===")

        X_interactions = X.copy()
        numerical_cols = X.select_dtypes(include=[np.number]).columns

        # Create pairwise interactions for top correlated features
        if len(numerical_cols) >= 2:
            # Calculate correlation matrix
            corr_matrix = X[numerical_cols].corr().abs()

            # Find top correlated pairs (excluding self-correlation)
            top_pairs = []
            for i in range(len(corr_matrix.columns)):
                for j in range(i+1, len(corr_matrix.columns)):
                    corr_val = corr_matrix.iloc[i, j]
                    if corr_val > 0.3:  # Threshold for creating interaction
                        col1, col2 = corr_matrix.columns[i], corr_matrix.columns[j]
                        top_pairs.append((col1, col2, corr_val))

            # Sort by correlation and take top 5
            top_pairs.sort(key=lambda x: x[2], reverse=True)
            top_pairs = top_pairs[:5]

            for col1, col2, corr in top_pairs:
                interaction_name = f"{col1}_x_{col2}"
                X_interactions[interaction_name] = X[col1] * X[col2]
                print(f"- 创建交互特征: {interaction_name} (相关系数: {corr:.3f})")

        # Create polynomial features for important numerical features
        if len(numerical_cols) > 0:
            # Select features with highest variance
            variances = X[numerical_cols].var()
            top_variance_features = variances.nlargest(3).index

            for col in top_variance_features:
                squared_name = f"{col}_squared"
                X_interactions[squared_name] = X[col] ** 2
                print(f"- 创建平方特征: {squared_name}")

        print(f"交互特征生成完成，总特征数: {X_interactions.shape[1]}")
        return X_interactions

    def train_models(self, X, y, test_size=0.2, cv_folds=5):
        """
        Train multiple regression models and compare performance

        Args:
            X (pd.DataFrame): Features
            y (pd.Series): Target
            test_size (float): Test set proportion
            cv_folds (int): Cross-validation folds

        Returns:
            dict: Model performance results
        """
        print("\n=== 模型训练与评估 ===")

        # Split data
        X_train, X_test, y_train, y_test = train_test_split(
            X, y, test_size=test_size, random_state=self.random_state
        )

        print(f"训练集大小: {X_train.shape[0]}")
        print(f"测试集大小: {X_test.shape[0]}")

        # Define models
        models = {
            'Linear Regression': LinearRegression(),
            'Ridge Regression': Ridge(alpha=1.0, random_state=self.random_state),
            'Lasso Regression': Lasso(alpha=1.0, random_state=self.random_state),
            'Decision Tree': DecisionTreeRegressor(random_state=self.random_state),
            'Random Forest': RandomForestRegressor(n_estimators=100, random_state=self.random_state),
            'Gradient Boosting': GradientBoostingRegressor(random_state=self.random_state)
        }

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

        results = {}

        for name, model in models.items():
            print(f"\n--- {name} ---")

            # Choose appropriate data scaling
            if 'Linear' in name or 'Ridge' in name or 'Lasso' in name:
                X_tr, X_te = X_train_scaled, X_test_scaled
            else:
                X_tr, X_te = X_train, X_test

            # Train model
            model.fit(X_tr, y_train)

            # Make predictions
            y_train_pred = model.predict(X_tr)
            y_test_pred = model.predict(X_te)

            # Calculate metrics
            train_r2 = r2_score(y_train, y_train_pred)
            test_r2 = r2_score(y_test, y_test_pred)
            train_mae = mean_absolute_error(y_train, y_train_pred)
            test_mae = mean_absolute_error(y_test, y_test_pred)
            train_rmse = np.sqrt(mean_squared_error(y_train, y_train_pred))
            test_rmse = np.sqrt(mean_squared_error(y_test, y_test_pred))

            # Cross-validation
            cv_scores = cross_val_score(model, X_tr, y_train, cv=cv_folds, scoring='r2')

            # Store results
            results[name] = {
                'model': model,
                'scaler': scaler if 'Linear' in name or 'Ridge' in name or 'Lasso' in name else None,
                'train_r2': train_r2,
                'test_r2': test_r2,
                'train_mae': train_mae,
                'test_mae': test_mae,
                'train_rmse': train_rmse,
                'test_rmse': test_rmse,
                'cv_mean': cv_scores.mean(),
                'cv_std': cv_scores.std(),
                'predictions': y_test_pred,
                'y_test': y_test
            }

            # Store model and scaler
            self.models[name] = model
            if 'Linear' in name or 'Ridge' in name or 'Lasso' in name:
                self.scalers[name] = scaler

            # Print results
            print(f"训练集 R²: {train_r2:.4f}")
            print(f"测试集 R²: {test_r2:.4f}")
            print(f"测试集 MAE: {test_mae:.4f}")
            print(f"测试集 RMSE: {test_rmse:.4f}")
            print(f"交叉验证 R²: {cv_scores.mean():.4f} (±{cv_scores.std():.4f})")

        # Select best model based on test R²
        best_model_name = max(results.keys(), key=lambda k: results[k]['test_r2'])
        self.best_model = results[best_model_name]['model']
        self.best_model_name = best_model_name

        print(f"\n🏆 最佳模型: {best_model_name} (测试集 R²: {results[best_model_name]['test_r2']:.4f})")

        self.results = results
        return results

    def get_feature_importance(self, model_name=None, top_n=10):
        """
        Get feature importance from trained models

        Args:
            model_name (str): Name of the model (use best model if None)
            top_n (int): Number of top features to return

        Returns:
            pd.DataFrame: Feature importance with rankings
        """
        if model_name is None:
            model_name = self.best_model_name

        if model_name not in self.models:
            raise ValueError(f"Model '{model_name}' not found in trained models")

        model = self.models[model_name]

        # Get feature importance based on model type
        if hasattr(model, 'feature_importances_'):
            # Tree-based models
            importance = model.feature_importances_
        elif hasattr(model, 'coef_'):
            # Linear models
            importance = np.abs(model.coef_)
        else:
            raise ValueError(f"Model '{model_name}' does not support feature importance")

        # Create importance DataFrame
        feature_importance = pd.DataFrame({
            'feature': self.feature_names,
            'importance': importance
        })

        # Sort and rank
        feature_importance = feature_importance.sort_values('importance', ascending=False)
        feature_importance['rank'] = range(1, len(feature_importance) + 1)
        feature_importance['importance_pct'] = (feature_importance['importance'] /
                                               feature_importance['importance'].sum() * 100)

        return feature_importance.head(top_n)

    def predict(self, X, model_name=None):
        """
        Make predictions using trained model

        Args:
            X (pd.DataFrame): Features for prediction
            model_name (str): Model to use (best model if None)

        Returns:
            np.array: Predictions
        """
        if model_name is None:
            model_name = self.best_model_name

        if model_name not in self.models:
            raise ValueError(f"Model '{model_name}' not found in trained models")

        model = self.models[model_name]

        # Apply scaling if needed
        if model_name in self.scalers:
            X_scaled = self.scalers[model_name].transform(X)
            return model.predict(X_scaled)
        else:
            return model.predict(X)

    def run_complete_analysis(self, file_path, target_column,
                            create_interactions=False, **kwargs):
        """
        Run complete regression analysis pipeline

        Args:
            file_path (str): Path to data file
            target_column (str): Target variable column
            create_interactions (bool): Whether to create interaction features
            **kwargs: Additional parameters

        Returns:
            dict: Complete analysis results
        """
        print("🚀 开始完整回归分析")
        print("=" * 50)

        # 1. Load and validate data
        X, y = self.load_and_validate_data(file_path, target_column, **kwargs)

        # 2. Preprocess data
        X_processed, y_processed = self.preprocess_data(X, y)

        # 3. Handle categorical features
        X_encoded = self.encode_categorical_features(X_processed)

        # 4. Create interaction features if requested
        if create_interactions:
            X_final = self.create_interaction_features(X_encoded)
        else:
            X_final = X_encoded

        # 6. Train models
        results = self.train_models(X_final, y_processed)

        # 7. Get feature importance
        feature_importance = self.get_feature_importance()

        # 8. Create summary
        summary = {
            'data_shape': X_final.shape,
            'feature_count': X_final.shape[1],
            'sample_count': len(y_processed),
            'best_model': self.best_model_name,
            'best_r2': results[self.best_model_name]['test_r2'],
            'best_mae': results[self.best_model_name]['test_mae'],
            'feature_importance': feature_importance
        }

        print(f"\n✅ 分析完成！")
        print(f"最佳模型: {self.best_model_name}")
        print(f"测试集 R²: {summary['best_r2']:.4f}")
        print(f"测试集 MAE: {summary['best_mae']:.4f}")

        return {
            'results': results,
            'summary': summary,
            'X_final': X_final,
            'y_final': y_processed,
            'feature_importance': feature_importance
        }

def main():
    """Example usage"""
    analyzer = RegressionAnalyzer()

    # Example with housing price data
    file_path = "房价预测数据.csv"
    target_column = "房价"

    if pd.io.common.file_exists(file_path):
        analysis_results = analyzer.run_complete_analysis(
            file_path,
            target_column,
            create_interactions=True
        )

        print(f"\n特征重要性排名:")
        print(analysis_results['feature_importance'])

    else:
        print(f"数据文件未找到: {file_path}")
        print("请提供正确的CSV数据文件路径")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Housing Price Prediction Example
Complete example of housing price prediction using multiple regression models
"""

import pandas as pd
import numpy as np
from core_regression import RegressionAnalyzer
from feature_engineering import FeatureEngineering
from model_evaluation import ModelEvaluator
from prediction_visualizer import PredictionVisualizer

def create_sample_housing_data():
    """Create realistic sample housing data"""
    print("创建示例房价数据...")

    np.random.seed(42)
    n_samples = 1000

    # Generate realistic housing features
    data = []

    for i in range(n_samples):
        # Basic features
        area = np.random.lognormal(4.5, 0.3)  # Log-normal distribution for area
        area = np.clip(area, 40, 300)  # Clip to reasonable range

        # Number of rooms depends on area
        max_rooms = int(area / 25)
        min_rooms = max(1, max_rooms - 2)
        rooms = np.random.randint(min_rooms, max_rooms + 1)

        # Bathrooms
        bathrooms = np.random.randint(1, min(rooms + 1, 4))

        # Floor information
        total_floors = np.random.choice([6, 12, 18, 30, 40], p=[0.3, 0.3, 0.2, 0.1, 0.1])
        floor = np.random.randint(1, total_floors + 1)

        # Building age
        build_year = np.random.choice(np.arange(1990, 2024))
        age = 2024 - build_year

        # Location features (distance in meters)
        subway_distance = np.random.exponential(800) * (1 + 0.3 * np.random.randn())
        subway_distance = np.clip(subway_distance, 50, 5000)

        school_distance = np.random.exponential(600) * (1 + 0.3 * np.random.randn())
        school_distance = np.clip(school_distance, 100, 4000)

        mall_distance = np.random.exponential(1000) * (1 + 0.3 * np.random.randn())
        mall_distance = np.clip(mall_distance, 200, 6000)

        # Categorical features with price impact
        decoration = np.random.choice(['毛坯', '简装修', '精装修', '豪华装修'],
                                   p=[0.1, 0.3, 0.4, 0.2])

        direction = np.random.choice(['南', '东南', '东', '西南', '西', '北'],
                                   p=[0.3, 0.2, 0.15, 0.15, 0.1, 0.1])

        estate_type = np.random.choice(['普通小区', '高档小区', '豪华小区'],
                                     p=[0.4, 0.4, 0.2])

        # Calculate base price using realistic factors
        # Base price per square meter
        base_price_per_sqm = 8000

        # Location adjustments
        location_multiplier = 1.0
        location_multiplier *= (1 - subway_distance / 10000)  # Closer to subway = higher price
        location_multiplier *= (1 - school_distance / 15000)  # Closer to school = higher price
        location_multiplier *= (1 - mall_distance / 20000)   # Closer to mall = higher price

        # Building quality adjustments
        decoration_multiplier = {'毛坯': 0.7, '简装修': 0.85, '精装修': 1.0, '豪华装修': 1.3}[decoration]
        direction_multiplier = {'南': 1.1, '东南': 1.05, '东': 1.0, '西南': 0.95, '西': 0.9, '北': 0.85}[direction]
        estate_multiplier = {'普通小区': 0.8, '高档小区': 1.0, '豪华小区': 1.25}[estate_type]

        # Age adjustment (newer buildings are more expensive)
        age_multiplier = 1.0 - (age / 100) * 0.3  # Max 30% reduction for very old buildings

        # Floor adjustment (higher floors in tall buildings are more expensive)
        floor_ratio = floor / total_floors
        floor_multiplier = 1.0 + (floor_ratio - 0.5) * 0.1 * (total_floors / 30)

        # Calculate final price
        price_per_sqm = (base_price_per_sqm *
                        location_multiplier *
                        decoration_multiplier *
                        direction_multiplier *
                        estate_multiplier *
                        age_multiplier *
                        floor_multiplier *
                        (1 + 0.1 * np.random.randn()))  # Add some noise

        total_price = price_per_sqm * area

        # Ensure reasonable price ranges
        total_price = np.clip(total_price, 300000, 15000000)

        data.append({
            '房屋ID': i + 1,
            '面积': round(area, 1),
            '房间数': rooms,
            '卫生间数': bathrooms,
            '楼层': floor,
            '总楼层': total_floors,
            '建造年份': build_year,
            '地铁距离': round(subway_distance),
            '学校距离': round(school_distance),
            '商场距离': round(mall_distance),
            '装修等级': decoration,
            '朝向': direction,
            '小区类型': estate_type,
            '房价': round(total_price, 2)
        })

    df = pd.DataFrame(data)
    df.to_csv('sample_housing_data.csv', index=False, encoding='utf-8-sig')

    print(f"生成了 {len(df)} 条房屋数据")
    print(f"价格范围: ￥{df['房价'].min():,.0f} - ￥{df['房价'].max():,.0f}")
    print(f"平均价格: ￥{df['房价'].mean():,.0f}")

    return df

def run_housing_price_example():
    """Run complete housing price prediction example"""
    print("🏠 开始房价预测示例")
    print("=" * 50)

    # 1. Create sample data
    housing_df = create_sample_housing_data()

    # 2. Initialize analyzers
    analyzer = RegressionAnalyzer()
    fe = FeatureEngineering()
    evaluator = ModelEvaluator()
    visualizer = PredictionVisualizer()

    # 3. Advanced feature engineering
    print("\n=== 高级特征工程 ===")

    # Create derived features
    housing_df['房龄'] = 2024 - housing_df['建造年份']
    housing_df['楼层比例'] = housing_df['楼层'] / housing_df['总楼层']
    housing_df['房间密度'] = housing_df['面积'] / housing_df['房间数']

    # Distance score (lower is better, so we invert)
    housing_df['交通便利性'] = (housing_df['地铁距离'].max() - housing_df['地铁距离']) / housing_df['地铁距离'].max()
    housing_df['学区便利性'] = (housing_df['学校距离'].max() - housing_df['学校距离']) / housing_df['学校距离'].max()
    housing_df['购物便利性'] = (housing_df['商场距离'].max() - housing_df['商场距离']) / housing_df['商场距离'].max()

    # Combined convenience score
    housing_df['综合便利性'] = (housing_df['交通便利性'] +
                               housing_df['学区便利性'] +
                               housing_df['购物便利性']) / 3

    # Price per square meter for analysis
    housing_df['单价'] = housing_df['房价'] / housing_df['面积']

    print(f"新增特征数量: {housing_df.shape[1] - 13}")
    print("新增特征包括: 房龄、楼层比例、房间密度、交通便利性、学区便利性、购物便利性、综合便利性、单价")

    # 4. Run regression analysis with interaction features
    print("\n=== 回归模型训练 ===")
    analysis_results = analyzer.run_complete_analysis(
        housing_df,
        '房价',
        create_interactions=True
    )

    # 5. Detailed model evaluation
    print("\n=== 模型详细评估 ===")

    # Residual analysis for best model
    best_model_name = analyzer.best_model_name
    best_results = analysis_results['results'][best_model_name]

    residual_analysis = evaluator.perform_residual_analysis(
        best_results['y_test'],
        best_results['predictions'],
        best_model_name
    )

    # Learning curve analysis
    learning_analysis = evaluator.analyze_learning_curves(
        analyzer.best_model,
        analysis_results['X_final'],
        analysis_results['y_final']
    )

    # 6. Feature importance analysis
    print("\n=== 特征重要性分析 ===")
    feature_importance = analysis_results['feature_importance']

    print("Top 10 影响房价的关键因素:")
    for idx, row in feature_importance.head(10).iterrows():
        print(f"  {idx + 1}. {row['feature']}: {row['importance']:.4f} ({row['importance_pct']:.1f}%)")

    # 7. Create visualizations
    print("\n=== 生成可视化分析 ===")

    # Comprehensive dashboard
    visualizer.create_comprehensive_dashboard(
        analysis_results['results'],
        feature_importance,
        save_path='housing_price_dashboard.png'
    )

    # Individual analysis plots
    visualizer.create_individual_analysis_plots(
        analysis_results['results'],
        output_dir='housing_analysis_plots'
    )

    # 8. Generate comprehensive report
    print("\n=== 生成分析报告 ===")
    evaluation_report = evaluator.generate_evaluation_report(
        analysis_results['results'],
        save_to_file=True
    )

    # 9. Business insights and price predictions
    print("\n=== 房价预测业务洞察 ===")

    # Analyze price predictions by different segments
    y_pred = best_results['predictions']
    y_true = best_results['y_test']

    # Calculate prediction accuracy by price ranges
    price_analysis = pd.DataFrame({
        '实际价格': y_true,
        '预测价格': y_pred,
        '绝对误差': np.abs(y_true - y_pred),
        '相对误差': np.abs(y_true - y_pred) / y_true * 100
    })

    # Create price segments
    price_analysis['价格区间'] = pd.cut(
        price_analysis['实际价格'],
        bins=[0, 500000, 1000000, 2000000, np.inf],
        labels=['经济型(≤50万)', '中档型(50-100万)', '高档型(100-200万)', '豪华型(>200万)']
    )

    print("\n不同价格区间的预测准确性:")
    for segment in price_analysis['价格区间'].cat.categories:
        segment_data = price_analysis[price_analysis['价格区间'] == segment]
        if len(segment_data) > 0:
            avg_error = segment_data['相对误差'].mean()
            print(f"  {segment}: 平均预测误差 {avg_error:.1f}%")

    # 10. Feature impact analysis
    print(f"\n=== 关键特征对房价的影响 ===")

    # Analyze categorical features impact
    categorical_analysis = {}

    # Decoration level impact
    decoration_impact = housing_df.groupby('装修等级')['单价'].mean().sort_values(ascending=False)
    print("\n装修等级对单价的影响:")
    for level, price in decoration_impact.items():
        print(f"  {level}: ￥{price:,.0f}/m²")

    # Direction impact
    direction_impact = housing_df.groupby('朝向')['单价'].mean().sort_values(ascending=False)
    print("\n朝向对单价的影响:")
    for direction, price in direction_impact.items():
        print(f"  {direction}: ￥{price:,.0f}/m²")

    # Estate type impact
    estate_impact = housing_df.groupby('小区类型')['单价'].mean().sort_values(ascending=False)
    print("\n小区类型对单价的影响:")
    for estate, price in estate_impact.items():
        print(f"  {estate}: ￥{price:,.0f}/m²")

    # 11. Model performance summary
    print(f"\n=== 模型性能总结 ===")
    best_metrics = analysis_results['results'][best_model_name]['metrics']
    print(f"最佳模型: {best_model_name}")
    print(f"R² 分数: {best_metrics['test_r2']:.4f}")
    print(f"平均绝对误差: ￥{best_metrics['test_mae']:,.0f}")
    print(f"均方根误差: ￥{best_metrics['test_rmse']:,.0f}")

    # Business interpretation
    mean_price = housing_df['房价'].mean()
    mae_percentage = (best_metrics['test_mae'] / mean_price) * 100
    print(f"平均预测误差百分比: {mae_percentage:.2f}%")

    print(f"\n=== 投资建议 ===")
    print("基于模型分析的投资建议:")
    print("1. 关注靠近地铁、学校、商圈的地段")
    print("2. 优先选择精装修或豪华装修的房产")
    print("3. 南向和东南朝向的房产具有更高价值")
    print("4. 高档小区和豪华小区有更好的增值潜力")
    print("5. 中高楼层（特别是电梯房）价格优势明显")

    print(f"\n✅ 房价预测分析完成！")
    print(f"生成的文件:")
    print(f"- sample_housing_data.csv: 示例房价数据")
    print(f"- housing_price_dashboard.png: 综合分析仪表板")
    print(f"- housing_analysis_plots/: 详细分析图表")
    print(f"- model_evaluation_report.md: 评估报告")

if __name__ == "__main__":
    run_housing_price_example()

房屋ID,面积,房间数,卫生间数,楼层,总楼层,建造年份,地铁距离,学校距离,商场距离,装修等级,朝向,小区类型,房价
1,120,3,2,15,30,2010,500,800,300,精装修,南,高档小区,850000
2,85,2,1,8,18,2015,300,600,200,简装修,东,普通小区,520000
3,150,4,2,25,35,2008,800,1200,1500,豪华装修,南,豪华小区,1280000
4,95,3,1,12,20,2012,200,500,100,精装修,南,高档小区,680000
5,200,5,3,1,40,2020,1000,1500,800,豪华装修,东南,豪华小区,1880000
6,75,2,1,5,15,2018,150,300,50,简装修,西,普通小区,450000
7,110,3,2,18,25,2011,600,900,400,精装修,南,高档小区,750000
8,160,4,2,20,30,2014,400,700,200,精装修,南,高档小区,980000
9,90,2,1,10,20,2016,250,450,80,简装修,东,普通小区,580000
10,130,3,2,22,32,2013,700,1100,500,精装修,南,高档小区,920000
11,100,3,1,14,25,2017,350,650,150,简装修,北,普通小区,620000
12,180,5,3,8,45,2009,900,1400,1200,豪华装修,东南,豪华小区,1450000
13,80,2,1,6,18,2019,180,320,70,简装修,西,普通小区,480000
14,140,4,2,28,35,2010,650,1000,600,精装修,南,高档小区,880000
15,105,3,2,16,28,2015,450,750,300,精装修,南,高档小区,720000
16,125,3,2,11,22,2012,550,850,250,精装修,东,高档小区,780000
17,95,2,1,9,16,2018,280,480,90,简装修,南,普通小区,590000
18,170,4,3,19,38,2011,750,1200,700,豪华装修,南,豪华小区,1380000
19,85,2,1,13,20,2016,220,400,110,简装修,东,普通小区,560000
20,155,4,2,26,33,2013,680,1050,550,精装修,南,高档小区,960000

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Advanced Feature Engineering Tools
Comprehensive feature engineering for regression analysis
"""

import pandas as pd
import numpy as np
from datetime import datetime, timedelta
from sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler
from sklearn.feature_selection import SelectKBest, f_regression, mutual_info_regression
import warnings
warnings.filterwarnings('ignore')

class FeatureEngineering:
    """Advanced feature engineering toolkit for regression analysis"""

    def __init__(self):
        """Initialize feature engineering toolkit"""
        self.scalers = {}
        self.encoders = {}
        self.feature_transformers = {}
        self.feature_selection_results = {}

    def extract_temporal_features(self, df, date_columns, reference_date=None):
        """
        Extract comprehensive temporal features from date columns

        Args:
            df (pd.DataFrame): Input dataframe
            date_columns (list): List of date column names
            reference_date (datetime): Reference date for calculations

        Returns:
            pd.DataFrame: DataFrame with temporal features
        """
        print("\n=== 时间特征工程 ===")

        df_temporal = df.copy()

        if reference_date is None:
            reference_date = datetime.now()

        for col in date_columns:
            if col in df_temporal.columns:
                # Ensure datetime format
                df_temporal[col] = pd.to_datetime(df_temporal[col], errors='coerce')

                # Basic temporal features
                df_temporal[f'{col}_year'] = df_temporal[col].dt.year
                df_temporal[f'{col}_month'] = df_temporal[col].dt.month
                df_temporal[f'{col}_day'] = df_temporal[col].dt.day
                df_temporal[f'{col}_dayofweek'] = df_temporal[col].dt.dayofweek
                df_temporal[f'{col}_dayofyear'] = df_temporal[col].dt.dayofyear
                df_temporal[f'{col}_quarter'] = df_temporal[col].dt.quarter
                df_temporal[f'{col}_week'] = df_temporal[col].dt.isocalendar().week

                # Cyclical features for better seasonality capture
                df_temporal[f'{col}_month_sin'] = np.sin(2 * np.pi * df_temporal[f'{col}_month'] / 12)
                df_temporal[f'{col}_month_cos'] = np.cos(2 * np.pi * df_temporal[f'{col}_month'] / 12)
                df_temporal[f'{col}_day_sin'] = np.sin(2 * np.pi * df_temporal[f'{col}_day'] / 31)
                df_temporal[f'{col}_day_cos'] = np.cos(2 * np.pi * df_temporal[f'{col}_day'] / 31)

                # Time since reference date
                df_temporal[f'{col}_days_since_ref'] = (reference_date - df_temporal[col]).dt.days
                df_temporal[f'{col}_weeks_since_ref'] = df_temporal[f'{col}_days_since_ref'] / 7
                df_temporal[f'{col}_months_since_ref'] = df_temporal[f'{col}_days_since_ref'] / 30.44

                # Weekend and holiday indicators
                df_temporal[f'{col}_is_weekend'] = (df_temporal[f'{col}_dayofweek'] >= 5).astype(int)

                print(f"- 为 {col} 创建 {14} 个时间特征")

        return df_temporal


    def create_polynomial_features(self, X, degree=2, interaction_only=False):
        """
        Create polynomial features for non-linear relationships

        Args:
            X (pd.DataFrame): Input features
            degree (int): Polynomial degree
            interaction_only (bool): Create only interaction terms

        Returns:
            pd.DataFrame: DataFrame with polynomial features
        """
        print(f"\n=== 多项式特征工程 (degree={degree}) ===")

        numerical_cols = X.select_dtypes(include=[np.number]).columns

        if len(numerical_cols) == 0:
            print("- 没有数值型特征可用于多项式变换")
            return X.copy()

        # Select features for polynomial transformation (avoid too many features)
        if len(numerical_cols) > 8:
            # Select features with highest variance
            variances = X[numerical_cols].var()
            selected_features = variances.nlargest(8).index
            print(f"- 从 {len(numerical_cols)} 个数值特征中选择方差最高的 8 个")
        else:
            selected_features = numerical_cols

        X_numerical = X[selected_features]

        # Create polynomial features
        poly = PolynomialFeatures(degree=degree, interaction_only=interaction_only, include_bias=False)
        X_poly = poly.fit_transform(X_numerical)

        # Get feature names
        poly_feature_names = poly.get_feature_names_out(selected_features)

        # Create DataFrame
        poly_df = pd.DataFrame(X_poly, columns=poly_feature_names, index=X.index)

        # Combine with original features
        X_combined = pd.concat([X.drop(columns=selected_features), poly_df], axis=1)

        print(f"- 创建了 {len(poly_feature_names)} 个多项式特征")
        print(f"- 总特征数: {X_combined.shape[1]}")

        return X_combined

    def create_aggregation_features(self, df, group_columns, agg_columns,
                                  agg_functions=['mean', 'std', 'min', 'max', 'count']):
        """
        Create aggregation features for grouped data

        Args:
            df (pd.DataFrame): Input dataframe
            group_columns (list): Columns to group by
            agg_columns (list): Columns to aggregate
            agg_functions (list): Aggregation functions

        Returns:
            pd.DataFrame: DataFrame with aggregation features
        """
        print(f"\n=== 聚合特征工程 ===")

        if not all(col in df.columns for col in group_columns):
            print("- 分组列不存在，跳过聚合特征")
            return df.copy()

        agg_df = df.copy()

        # Create aggregation features
        for group_col in group_columns:
            if group_col not in df.columns:
                continue

            for agg_col in agg_columns:
                if agg_col not in df.columns:
                    continue

                for func in agg_functions:
                    if func == 'count':
                        # Count is special case
                        feature_name = f"{group_col}_{agg_col}_count"
                        agg_df[feature_name] = df.groupby(group_col)[agg_col].transform('count')
                    else:
                        feature_name = f"{group_col}_{agg_col}_{func}"
                        agg_df[feature_name] = df.groupby(group_col)[agg_col].transform(func)

        print(f"- 为分组列 {group_columns} 创建聚合特征")
        print(f"- 新增特征数: {agg_df.shape[1] - df.shape[1]}")

        return agg_df

    def create_ratio_features(self, X, ratio_pairs):
        """
        Create ratio features from pairs of columns

        Args:
            X (pd.DataFrame): Input features
            ratio_pairs (list): List of tuples (numerator, denominator, new_name)

        Returns:
            pd.DataFrame: DataFrame with ratio features
        """
        print(f"\n=== 比例特征工程 ===")

        X_ratios = X.copy()
        created_features = 0

        for numerator, denominator, new_name in ratio_pairs:
            if numerator in X.columns and denominator in X.columns:
                # Handle division by zero
                with np.errstate(divide='ignore', invalid='ignore'):
                    X_ratios[new_name] = np.where(
                        X[denominator] != 0,
                        X[numerator] / X[denominator],
                        0  # or np.nan, depending on preference
                    )
                created_features += 1
                print(f"- 创建比例特征: {new_name} = {numerator} / {denominator}")
            else:
                print(f"- 跳过 {new_name}: 缺少必要的列")

        print(f"- 总共创建了 {created_features} 个比例特征")
        return X_ratios

    def create_binning_features(self, X, binning_config):
        """
        Create binned/discretized features

        Args:
            X (pd.DataFrame): Input features
            binning_config (dict): Configuration for binning

        Returns:
            pd.DataFrame: DataFrame with binned features
        """
        print(f"\n=== 分箱特征工程 ===")

        X_binned = X.copy()

        for col, config in binning_config.items():
            if col not in X.columns:
                print(f"- 跳过 {col}: 列不存在")
                continue

            if isinstance(config, dict):
                bin_type = config.get('type', 'quantile')
                bins = config.get('bins', 5)
                labels = config.get('labels', None)
                suffix = config.get('suffix', '_binned')
            else:
                # Simple configuration
                bin_type = 'quantile'
                bins = config
                labels = None
                suffix = '_binned'

            try:
                if bin_type == 'quantile':
                    X_binned[f"{col}{suffix}"] = pd.qcut(X[col], q=bins, labels=labels, duplicates='drop')
                elif bin_type == 'uniform':
                    X_binned[f"{col}{suffix}"] = pd.cut(X[col], bins=bins, labels=labels)
                else:
                    print(f"- 不支持的分箱类型: {bin_type}")
                    continue

                print(f"- 为 {col} 创建分箱特征 (类型: {bin_type}, 分箱数: {bins})")

            except Exception as e:
                print(f"- 为 {col} 创建分箱特征失败: {str(e)}")

        return X_binned

    def select_features(self, X, y, method='f_regression', k=10):
        """
        Feature selection using statistical methods

        Args:
            X (pd.DataFrame): Features
            y (pd.Series): Target
            method (str): Selection method
            k (int): Number of features to select

        Returns:
            tuple: (selected_features, selection_scores)
        """
        print(f"\n=== 特征选择 (method={method}, k={k}) ===")

        # Ensure X is numeric
        X_numeric = X.select_dtypes(include=[np.number])

        if X_numeric.shape[1] == 0:
            print("- 没有数值型特征可供选择")
            return X, {}

        # Feature selection
        if method == 'f_regression':
            selector = SelectKBest(score_func=f_regression, k=min(k, X_numeric.shape[1]))
        elif method == 'mutual_info':
            selector = SelectKBest(score_func=mutual_info_regression, k=min(k, X_numeric.shape[1]))
        else:
            raise ValueError(f"Unsupported selection method: {method}")

        # Fit selector
        X_selected = selector.fit_transform(X_numeric, y)
        selected_features = X_numeric.columns[selector.get_support()]
        selection_scores = selector.scores_[selector.get_support()]

        # Create results DataFrame
        results_df = pd.DataFrame({
            'feature': selected_features,
            'score': selection_scores
        }).sort_values('score', ascending=False)

        # Create selected features DataFrame
        X_selected_df = pd.DataFrame(X_selected, columns=selected_features, index=X.index)

        # Add non-numeric columns back
        non_numeric_cols = X.select_dtypes(exclude=[np.number]).columns
        for col in non_numeric_cols:
            X_selected_df[col] = X[col]

        print(f"- 从 {X_numeric.shape[1]} 个特征中选择了 {len(selected_features)} 个")
        print("- Top 10 特征:")
        for _, row in results_df.head(10).iterrows():
            print(f"  {row['feature']}: {row['score']:.4f}")

        self.feature_selection_results[method] = results_df

        return X_selected_df, results_df

    def scale_features(self, X, method='standard', columns=None):
        """
        Scale numerical features

        Args:
            X (pd.DataFrame): Input features
            method (str): Scaling method
            columns (list): Columns to scale (all numeric if None)

        Returns:
            pd.DataFrame: Scaled features
        """
        print(f"\n=== 特征缩放 (method={method}) ===")

        X_scaled = X.copy()

        if columns is None:
            columns = X.select_dtypes(include=[np.number]).columns

        if method == 'standard':
            scaler = StandardScaler()
        elif method == 'minmax':
            scaler = MinMaxScaler()
        elif method == 'robust':
            scaler = RobustScaler()
        else:
            raise ValueError(f"Unsupported scaling method: {method}")

        # Fit and transform
        X_scaled[columns] = scaler.fit_transform(X[columns])

        # Store scaler
        self.scalers[method] = scaler

        print(f"- 使用 {method} 方法缩放 {len(columns)} 个数值特征")

        return X_scaled

def main():
    """Example usage"""
    fe = FeatureEngineering()

    # Create sample data
    sample_data = pd.DataFrame({
        'feature1': [1, 2, 3, 4, 5],
        'feature2': [10, 20, 30, 40, 50],
        'feature3': ['A', 'B', 'A', 'C', 'B'],
        'date': pd.date_range('2023-01-01', periods=5),
        'target': [100, 150, 200, 250, 300]
    })

    # Test temporal features
    temporal_result = fe.create_temporal_features(
        sample_data,
        date_columns=['date']
    )
    print("时间特征结果:")
    print(temporal_result)

    # Test polynomial features
    X = sample_data[['feature1', 'feature2']]
    poly_result = fe.create_polynomial_features(X, degree=2)
    print("\n多项式特征结果:")
    print(poly_result)

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Advanced Model Evaluation and Comparison
Comprehensive model evaluation with diagnostics and visualization
"""

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.metrics import (r2_score, mean_absolute_error, mean_squared_error,
                           mean_absolute_percentage_error, explained_variance_score)
from sklearn.model_selection import learning_curve, validation_curve
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

class ModelEvaluator:
    """Advanced model evaluation toolkit with comprehensive diagnostics"""

    def __init__(self, chinese_font='SimHei'):
        """Initialize model evaluator"""
        self.chinese_font = chinese_font
        plt.rcParams['font.sans-serif'] = [chinese_font]
        plt.rcParams['axes.unicode_minus'] = False

        self.evaluation_results = {}
        self.comparison_results = {}

    def calculate_comprehensive_metrics(self, y_true, y_pred, model_name="Model"):
        """
        Calculate comprehensive evaluation metrics

        Args:
            y_true (array-like): True values
            y_pred (array-like): Predicted values
            model_name (str): Name of the model

        Returns:
            dict: Comprehensive metrics
        """
        # Basic metrics
        r2 = r2_score(y_true, y_pred)
        mae = mean_absolute_error(y_true, y_pred)
        rmse = np.sqrt(mean_squared_error(y_true, y_pred))
        mse = mean_squared_error(y_true, y_pred)

        # Additional metrics
        mape = mean_absolute_percentage_error(y_true, y_pred)
        explained_var = explained_variance_score(y_true, y_pred)

        # Custom metrics
        residuals = y_true - y_pred
        mean_residual = np.mean(residuals)
        std_residual = np.std(residuals)

        # Symmetric Mean Absolute Percentage Error
        smape = np.mean(2 * np.abs(residuals) / (np.abs(y_true) + np.abs(y_pred))) * 100

        # Relative metrics
        y_mean = np.mean(y_true)
        relative_mae = mae / y_mean * 100
        relative_rmse = rmse / y_mean * 100

        metrics = {
            'model_name': model_name,
            'r2_score': r2,
            'mae': mae,
            'rmse': rmse,
            'mse': mse,
            'mape': mape,
            'smape': smape,
            'explained_variance': explained_var,
            'mean_residual': mean_residual,
            'std_residual': std_residual,
            'relative_mae': relative_mae,
            'relative_rmse': relative_rmse,
            'n_samples': len(y_true),
            'mean_target': y_mean,
            'std_target': np.std(y_true)
        }

        self.evaluation_results[model_name] = metrics
        return metrics

    def perform_residual_analysis(self, y_true, y_pred, model_name="Model", save_plots=True):
        """
        Perform comprehensive residual analysis

        Args:
            y_true (array-like): True values
            y_pred (array-like): Predicted values
            model_name (str): Name of the model
            save_plots (bool): Whether to save plots

        Returns:
            dict: Residual analysis results
        """
        print(f"\n=== {model_name} 残差分析 ===")

        residuals = y_true - y_pred
        standardized_residuals = residuals / np.std(residuals)
        predicted_values = y_pred

        # Statistical tests
        # 1. Normality test for residuals
        shapiro_stat, shapiro_p = stats.shapiro(residuals[:5000] if len(residuals) > 5000 else residuals)

        # 2. Homoscedasticity test (Breusch-Pagan)
        # Simplified version - check correlation between residuals and predicted values
        homoscedasticity_corr = np.corrcoef(np.abs(residuals), predicted_values)[0, 1]

        # 3. Independence test (Durbin-Watson approximation)
        if len(residuals) > 1:
            dw_stat = np.sum(np.diff(residuals)**2) / np.sum(residuals**2)
        else:
            dw_stat = 2.0  # Perfect independence

        # Create diagnostic plots
        if save_plots:
            self._create_residual_plots(y_true, y_pred, model_name)

        results = {
            'model_name': model_name,
            'residuals_mean': np.mean(residuals),
            'residuals_std': np.std(residuals),
            'residuals_skewness': stats.skew(residuals),
            'residuals_kurtosis': stats.kurtosis(residuals),
            'shapiro_stat': shapiro_stat,
            'shapiro_p_value': shapiro_p,
            'normality_test': 'Pass' if shapiro_p > 0.05 else 'Fail',
            'homoscedasticity_correlation': homoscedasticity_corr,
            'homoscedasticity_test': 'Pass' if abs(homoscedasticity_corr) < 0.3 else 'Fail',
            'durbin_watson_stat': dw_stat,
            'independence_test': 'Pass' if 1.5 < dw_stat < 2.5 else 'Fail'
        }

        # Print summary
        print(f"残差均值: {results['residuals_mean']:.4f}")
        print(f"残差标准差: {results['residuals_std']:.4f}")
        print(f"正态性检验: {results['normality_test']} (p-value: {shapiro_p:.4f})")
        print(f"同方差性检验: {results['homoscedasticity_test']} (相关性: {homoscedasticity_corr:.4f})")
        print(f"独立性检验: {results['independence_test']} (DW统计量: {dw_stat:.4f})")

        return results

    def _create_residual_plots(self, y_true, y_pred, model_name):
        """Create comprehensive residual diagnostic plots"""
        residuals = y_true - y_pred
        predicted_values = y_pred

        # Create subplot layout
        fig, axes = plt.subplots(2, 3, figsize=(18, 12))
        fig.suptitle(f'{model_name} - 残差诊断图', fontsize=16, fontweight='bold')

        # 1. Residuals vs Fitted
        axes[0, 0].scatter(predicted_values, residuals, alpha=0.6, s=20)
        axes[0, 0].axhline(y=0, color='r', linestyle='--')
        axes[0, 0].set_xlabel('预测值')
        axes[0, 0].set_ylabel('残差')
        axes[0, 0].set_title('残差 vs 拟合值')
        axes[0, 0].grid(True, alpha=0.3)

        # 2. Q-Q Plot
        stats.probplot(residuals, dist="norm", plot=axes[0, 1])
        axes[0, 1].set_title('Q-Q图 (正态性检验)')
        axes[0, 1].grid(True, alpha=0.3)

        # 3. Scale-Location Plot
        axes[0, 2].scatter(predicted_values, np.sqrt(np.abs(residuals)), alpha=0.6, s=20)
        axes[0, 2].set_xlabel('预测值')
        axes[0, 2].set_ylabel('√|残差|')
        axes[0, 2].set_title('Scale-Location图')
        axes[0, 2].grid(True, alpha=0.3)

        # 4. Histogram of Residuals
        axes[1, 0].hist(residuals, bins=30, alpha=0.7, color='skyblue', edgecolor='black')
        axes[1, 0].set_xlabel('残差')
        axes[1, 0].set_ylabel('频次')
        axes[1, 0].set_title('残差直方图')
        axes[1, 0].grid(True, alpha=0.3)

        # 5. Actual vs Predicted
        min_val = min(y_true.min(), y_pred.min())
        max_val = max(y_true.max(), y_pred.max())
        axes[1, 1].scatter(y_true, y_pred, alpha=0.6, s=20)
        axes[1, 1].plot([min_val, max_val], [min_val, max_val], 'r--', alpha=0.8)
        axes[1, 1].set_xlabel('实际值')
        axes[1, 1].set_ylabel('预测值')
        axes[1, 1].set_title('实际值 vs 预测值')
        axes[1, 1].grid(True, alpha=0.3)

        # 6. Residuals vs Order (if index is available)
        axes[1, 2].plot(range(len(residuals)), residuals, alpha=0.6)
        axes[1, 2].axhline(y=0, color='r', linestyle='--')
        axes[1, 2].set_xlabel('观测顺序')
        axes[1, 2].set_ylabel('残差')
        axes[1, 2].set_title('残差 vs 观测顺序')
        axes[1, 2].grid(True, alpha=0.3)

        plt.tight_layout()
        plt.savefig(f'{model_name}_residual_analysis.png', dpi=300, bbox_inches='tight')
        plt.close()

    def compare_models(self, model_results, save_comparison=True):
        """
        Compare multiple models with comprehensive metrics

        Args:
            model_results (dict): Dictionary with model results
            save_comparison (bool): Whether to save comparison plots

        Returns:
            pd.DataFrame: Model comparison results
        """
        print("\n=== 模型比较分析 ===")

        # Create comparison DataFrame
        comparison_data = []

        for model_name, results in model_results.items():
            if 'metrics' in results:
                metrics = results['metrics']
            elif hasattr(results, '__dict__'):
                # Handle model objects
                metrics = self.calculate_comprehensive_metrics(
                    results['y_test'], results['predictions'], model_name
                )
            else:
                continue

            comparison_data.append(metrics)

        comparison_df = pd.DataFrame(comparison_data)

        if len(comparison_df) == 0:
            print("没有可比较的模型结果")
            return pd.DataFrame()

        # Rank models by different metrics
        ranking_metrics = ['r2_score', 'mae', 'rmse', 'mape']
        for metric in ranking_metrics:
            if metric in comparison_df.columns:
                if metric == 'r2_score':
                    comparison_df[f'{metric}_rank'] = comparison_df[metric].rank(ascending=False)
                else:
                    comparison_df[f'{metric}_rank'] = comparison_df[metric].rank(ascending=True)

        # Calculate overall rank
        rank_columns = [col for col in comparison_df.columns if col.endswith('_rank')]
        if rank_columns:
            comparison_df['overall_rank'] = comparison_df[rank_columns].mean(axis=1)
            comparison_df = comparison_df.sort_values('overall_rank')

        self.comparison_results = comparison_df

        # Print comparison summary
        print("\n模型性能排名:")
        for idx, row in comparison_df.iterrows():
            print(f"{idx + 1}. {row['model_name']}: "
                  f"R²={row.get('r2_score', 'N/A'):.4f}, "
                  f"MAE={row.get('mae', 'N/A'):.4f}, "
                  f"RMSE={row.get('rmse', 'N/A'):.4f}")

        if save_comparison:
            self._create_comparison_plots(comparison_df)

        return comparison_df

    def _create_comparison_plots(self, comparison_df):
        """Create model comparison visualizations"""
        fig, axes = plt.subplots(2, 2, figsize=(15, 12))
        fig.suptitle('模型性能比较', fontsize=16, fontweight='bold')

        # 1. R² Score Comparison
        if 'r2_score' in comparison_df.columns:
            bars = axes[0, 0].bar(comparison_df['model_name'], comparison_df['r2_score'])
            axes[0, 0].set_title('R² 分数比较')
            axes[0, 0].set_ylabel('R² 分数')
            axes[0, 0].tick_params(axis='x', rotation=45)
            axes[0, 0].grid(True, alpha=0.3)

            # Add value labels on bars
            for bar, value in zip(bars, comparison_df['r2_score']):
                axes[0, 0].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01,
                               f'{value:.4f}', ha='center', va='bottom')

        # 2. MAE Comparison
        if 'mae' in comparison_df.columns:
            bars = axes[0, 1].bar(comparison_df['model_name'], comparison_df['mae'], color='orange')
            axes[0, 1].set_title('平均绝对误差比较')
            axes[0, 1].set_ylabel('MAE')
            axes[0, 1].tick_params(axis='x', rotation=45)
            axes[0, 1].grid(True, alpha=0.3)

            # Add value labels
            for bar, value in zip(bars, comparison_df['mae']):
                axes[0, 1].text(bar.get_x() + bar.get_width()/2, bar.get_height(),
                               f'{value:.2f}', ha='center', va='bottom')

        # 3. RMSE Comparison
        if 'rmse' in comparison_df.columns:
            bars = axes[1, 0].bar(comparison_df['model_name'], comparison_df['rmse'], color='green')
            axes[1, 0].set_title('均方根误差比较')
            axes[1, 0].set_ylabel('RMSE')
            axes[1, 0].tick_params(axis='x', rotation=45)
            axes[1, 0].grid(True, alpha=0.3)

            # Add value labels
            for bar, value in zip(bars, comparison_df['rmse']):
                axes[1, 0].text(bar.get_x() + bar.get_width()/2, bar.get_height(),
                               f'{value:.2f}', ha='center', va='bottom')

        # 4. Overall Ranking
        if 'overall_rank' in comparison_df.columns:
            bars = axes[1, 1].bar(comparison_df['model_name'], comparison_df['overall_rank'], color='purple')
            axes[1, 1].set_title('综合排名 (越小越好)')
            axes[1, 1].set_ylabel('排名分数')
            axes[1, 1].tick_params(axis='x', rotation=45)
            axes[1, 1].grid(True, alpha=0.3)
            axes[1, 1].invert_yaxis()  # Lower rank is better

            # Add value labels
            for bar, value in zip(bars, comparison_df['overall_rank']):
                axes[1, 1].text(bar.get_x() + bar.get_width()/2, bar.get_height(),
                               f'{value:.1f}', ha='center', va='bottom')

        plt.tight_layout()
        plt.savefig('model_comparison.png', dpi=300, bbox_inches='tight')
        plt.close()

    def analyze_learning_curves(self, model, X, y, cv=5, train_sizes=None, save_plot=True):
        """
        Analyze learning curves to diagnose model performance

        Args:
            model: Trained model object
            X (pd.DataFrame): Features
            y (pd.Series): Target
            cv (int): Cross-validation folds
            train_sizes (array): Training sizes to evaluate
            save_plot (bool): Whether to save the plot

        Returns:
            dict: Learning curve analysis results
        """
        print("\n=== 学习曲线分析 ===")

        if train_sizes is None:
            train_sizes = np.linspace(0.1, 1.0, 10)

        # Calculate learning curves
        train_sizes_abs, train_scores, val_scores = learning_curve(
            model, X, y, cv=cv, train_sizes=train_sizes,
            scoring='r2', random_state=42, n_jobs=-1
        )

        # Calculate statistics
        train_mean = np.mean(train_scores, axis=1)
        train_std = np.std(train_scores, axis=1)
        val_mean = np.mean(val_scores, axis=1)
        val_std = np.std(val_scores, axis=1)

        # Analyze curve characteristics
        final_train_score = train_mean[-1]
        final_val_score = val_mean[-1]
        overfitting_gap = final_train_score - final_val_score

        # Determine if model is overfitting, underfitting, or well-balanced
        if overfitting_gap > 0.1:
            model_status = "过拟合"
            recommendation = "增加正则化、获取更多数据、简化模型"
        elif final_val_score < 0.5:
            model_status = "欠拟合"
            recommendation = "增加特征复杂度、使用更强大的模型"
        else:
            model_status = "平衡良好"
            recommendation = "模型性能良好，可以考虑微调"

        results = {
            'model_status': model_status,
            'final_train_score': final_train_score,
            'final_val_score': final_val_score,
            'overfitting_gap': overfitting_gap,
            'recommendation': recommendation,
            'train_sizes': train_sizes_abs,
            'train_scores_mean': train_mean,
            'train_scores_std': train_std,
            'val_scores_mean': val_mean,
            'val_scores_std': val_std
        }

        print(f"模型状态: {model_status}")
        print(f"最终训练分数: {final_train_score:.4f}")
        print(f"最终验证分数: {final_val_score:.4f}")
        print(f"过拟合程度: {overfitting_gap:.4f}")
        print(f"建议: {recommendation}")

        if save_plot:
            self._create_learning_curve_plot(results)

        return results

    def _create_learning_curve_plot(self, results):
        """Create learning curve visualization"""
        plt.figure(figsize=(10, 6))

        # Plot training scores
        plt.plot(results['train_sizes'], results['train_scores_mean'], 'o-',
                color='blue', label='训练分数')
        plt.fill_between(results['train_sizes'],
                        results['train_scores_mean'] - results['train_scores_std'],
                        results['train_scores_mean'] + results['train_scores_std'],
                        alpha=0.1, color='blue')

        # Plot validation scores
        plt.plot(results['train_sizes'], results['val_scores_mean'], 'o-',
                color='red', label='验证分数')
        plt.fill_between(results['train_sizes'],
                        results['val_scores_mean'] - results['val_scores_std'],
                        results['val_scores_mean'] + results['val_scores_std'],
                        alpha=0.1, color='red')

        plt.xlabel('训练样本数')
        plt.ylabel('R² 分数')
        plt.title('学习曲线分析')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.ylim(0, 1)

        plt.tight_layout()
        plt.savefig('learning_curves.png', dpi=300, bbox_inches='tight')
        plt.close()

    def generate_evaluation_report(self, model_results, save_to_file=True):
        """
        Generate comprehensive evaluation report

        Args:
            model_results (dict): Dictionary with model results
            save_to_file (bool): Whether to save report to file

        Returns:
            str: Evaluation report
        """
        print("\n=== 生成评估报告 ===")

        # Compare models
        comparison_df = self.compare_models(model_results, save_comparison=False)

        if len(comparison_df) == 0:
            return "没有模型结果可用于生成报告"

        # Generate report
        report = f"""
# 回归模型评估报告
# Regression Model Evaluation Report

生成时间: {pd.Timestamp.now().strftime('%Y-%m-%d %H:%M:%S')}

## 模型性能排名 | Model Performance Ranking

| 排名 | 模型名称 | R² 分数 | MAE | RMSE | MAPE | 综合评分 |
|------|----------|---------|-----|------|------|----------|
"""

        for idx, (_, row) in enumerate(comparison_df.iterrows(), 1):
            report += f"| {idx} | {row['model_name']} | {row.get('r2_score', 'N/A'):.4f} | {row.get('mae', 'N/A'):.2f} | {row.get('rmse', 'N/A'):.2f} | {row.get('mape', 'N/A'):.4f} | {row.get('overall_rank', 'N/A'):.2f} |\n"

        # Best model analysis
        best_model = comparison_df.iloc[0]
        report += f"""

## 最佳模型分析 | Best Model Analysis

**模型名称:** {best_model['model_name']}

### 性能指标 | Performance Metrics
- **R² 分数:** {best_model.get('r2_score', 'N/A'):.4f}
- **平均绝对误差 (MAE):** {best_model.get('mae', 'N/A'):.2f}
- **均方根误差 (RMSE):** {best_model.get('rmse', 'N/A'):.2f}
- **平均绝对百分比误差 (MAPE):** {best_model.get('mape', 'N/A'):.4f}
- **解释方差比:** {best_model.get('explained_variance', 'N/A'):.4f}

### 模型特征 | Model Characteristics
- **样本数量:** {best_model.get('n_samples', 'N/A')}
- **目标变量均值:** {best_model.get('mean_target', 'N/A'):.2f}
- **相对误差 (MAE):** {best_model.get('relative_mae', 'N/A'):.2f}%

## 诊断建议 | Diagnostic Recommendations

"""

        # Add specific recommendations based on model performance
        if best_model.get('r2_score', 0) > 0.8:
            report += "✅ **模型表现优秀** - 模型具有很强的预测能力\n"
        elif best_model.get('r2_score', 0) > 0.6:
            report += "⚠️ **模型表现良好** - 模型具有较好的预测能力，仍有改进空间\n"
        else:
            report += "❌ **模型表现需要改进** - 建议尝试特征工程、更复杂的模型或数据预处理\n"

        # Overfitting analysis
        if 'train_r2' in best_model and 'test_r2' in best_model:
            overfitting = best_model['train_r2'] - best_model['test_r2']
            if overfitting > 0.2:
                report += f"\n⚠️ **存在过拟合** - 训练和测试分数差距较大 ({overfitting:.3f})\n"
                report += "   建议: 增加正则化、使用交叉验证或获取更多训练数据\n"

        # Feature engineering suggestions
        report += f"""

## 改进建议 | Improvement Recommendations

### 1. 数据层面 | Data Level
- 考虑添加更多相关特征
- 处理异常值和缺失值
- 检查数据质量和一致性

### 2. 特征工程 | Feature Engineering
- 尝试特征交互和多项式特征
- 使用特征选择方法
- 考虑特征缩放和标准化

### 3. 模型层面 | Model Level
- 尝试不同的回归算法
- 调整超参数
- 使用集成学习方法

### 4. 验证方法 | Validation
- 使用时间序列交叉验证（如果数据有时间依赖）
- 考虑分层抽样
- 增加验证的稳定性

## 可视化文件 | Generated Visualizations

以下可视化文件已生成:
- `model_comparison.png`: 模型性能比较图
- `learning_curves.png`: 学习曲线分析图
- `*_residual_analysis.png`: 各模型残差诊断图

---

*报告由回归模型评估系统自动生成*
"""

        if save_to_file:
            with open('model_evaluation_report.md', 'w', encoding='utf-8') as f:
                f.write(report)
            print("评估报告已保存到: model_evaluation_report.md")

        return report

def main():
    """Example usage"""
    evaluator = ModelEvaluator()

    # Create sample data
    np.random.seed(42)
    y_true = np.random.normal(100, 20, 1000)
    y_pred_1 = y_true + np.random.normal(0, 5, 1000)  # Good model
    y_pred_2 = y_true + np.random.normal(0, 15, 1000)  # Poor model

    # Calculate metrics
    metrics_1 = evaluator.calculate_comprehensive_metrics(y_true, y_pred_1, "Good Model")
    metrics_2 = evaluator.calculate_comprehensive_metrics(y_true, y_pred_2, "Poor Model")

    print("Good Model Metrics:")
    for key, value in metrics_1.items():
        if key != 'model_name':
            print(f"  {key}: {value:.4f}")

    # Residual analysis
    residual_results_1 = evaluator.perform_residual_analysis(y_true, y_pred_1, "Good Model")

if __name__ == "__main__":
    main()