import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime

# 設置中文字體
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False

# 生成更豐富的模擬數據
np.random.seed(42)
n_samples = 1000

raw_loan_data = pd.DataFrame({
    'applicant_id': range(1, n_samples + 1),
    'age': np.random.randint(20, 60, n_samples),
    'annual_income': np.random.normal(80000, 30000, n_samples).clip(30000, 200000),
    'debt_amount': np.random.normal(150000, 80000, n_samples).clip(0, 500000),
    'work_years': np.random.exponential(5, n_samples).clip(0, 30),
    'phone_tenure_months': np.random.exponential(24, n_samples).clip(1, 120),
    'late_payments_3y': np.random.poisson(0.5, n_samples),
    'credit_card_utilization': np.random.beta(2, 5, n_samples),  # 信用卡使用率
    'number_of_credit_cards': np.random.poisson(2, n_samples),
    'has_mortgage': np.random.binomial(1, 0.3, n_samples),
    'application_hour': np.random.randint(0, 24, n_samples)
})

# 生成目標變量（違約標籤）
def generate_default_label(df):
    """基於業務規則生成違約標籤"""
    default_proba = (
        (df['age'] < 30) * 0.3 +
        (df['debt_amount'] / df['annual_income'] > 3) * 0.4 +
        (df['late_payments_3y'] >= 3) * 0.6 +
        (df['work_years'] < 2) * 0.2 +
        (df['credit_card_utilization'] > 0.8) * 0.3 +
        np.random.normal(0, 0.1, len(df))
    )
    return (default_proba > 0.5).astype(int)

raw_loan_data['default'] = generate_default_label(raw_loan_data)

def create_risk_features(df):
    """創建金融風控特徵"""
    df_fe = df.copy()
    
    print("開始特徵工程...")
    
    # 1. 單變量分析驅動的特徵
    print("1. 創建單變量特徵...")
    df_fe['has_serious_delinquency'] = (df_fe['late_payments_3y'] >= 3).astype(int)
    df_fe['high_credit_utilization'] = (df_fe['credit_card_utilization'] > 0.8).astype(int)
    df_fe['is_young_applicant'] = (df_fe['age'] < 25).astype(int)
    
    # 2. 多變量分析驅動的特徵
    print("2. 創建多變量交互特徵...")
    df_fe['debt_to_income_ratio'] = df_fe['debt_amount'] / df_fe['annual_income']
    df_fe['young_high_debt_risk'] = ((df_fe['age'] < 30) & (df_fe['debt_to_income_ratio'] > 2.5)).astype(int)
    df_fe['delinquency_income_risk'] = ((df_fe['late_payments_3y'] > 1) & (df_fe['annual_income'] < 50000)).astype(int)
    
    # 3. 業務分析驅動的特徵
    print("3. 創建業務特徵...")
    # 穩定性特徵
    df_fe['job_stability_score'] = np.log1p(df_fe['work_years']) / np.log1p(30)
    df_fe['life_stability_score'] = np.minimum(df_fe['phone_tenure_months'] / 60, 1.0)
    
    # 申請行為特徵
    df_fe['is_off_hours_application'] = ((df_fe['application_hour'] < 6) | (df_fe['application_hour'] > 22)).astype(int)
    
    # 信用行為特徵
    df_fe['credit_behavior_score'] = (
        (1 - df_fe['credit_card_utilization']) * 0.4 +
        (df_fe['number_of_credit_cards'] <= 5).astype(int) * 0.3 +
        (df_fe['late_payments_3y'] == 0).astype(int) * 0.3
    )
    
    # 4. 綜合風險評分
    print("4. 計算綜合風險評分...")
    df_fe['comprehensive_risk_score'] = (
        df_fe['has_serious_delinquency'] * 0.25 +
        np.minimum(df_fe['debt_to_income_ratio'], 5) * 0.15 +
        df_fe['young_high_debt_risk'] * 0.15 +
        (1 - df_fe['job_stability_score']) * 0.15 +
        (1 - df_fe['life_stability_score']) * 0.1 +
        df_fe['is_off_hours_application'] * 0.05 +
        (1 - df_fe['credit_behavior_score']) * 0.15
    )
    
    # 風險等級劃分
    df_fe['risk_level'] = pd.cut(df_fe['comprehensive_risk_score'], 
                               bins=[0, 0.3, 0.6, 1.0],
                               labels=['低風險', '中風險', '高風險'])
    
    return df_fe

# 執行特徵工程
engineered_data = create_risk_features(raw_loan_data)

print("\n特徵工程完成！")
print(f"原始特徵數量: {len(raw_loan_data.columns)}")
print(f"工程後特徵數量: {len(engineered_data.columns)}")
print(f"違約率: {engineered_data['default'].mean():.2%}")

# 創建特徵工程可視化
def plot_feature_engineering_analysis(df):
    """繪製特徵工程分析圖"""
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))
    fig.suptitle('金融風控特徵工程深度分析', fontsize=16, fontweight='bold')
    
    # 1. 單變量分析 - 債務收入比分佈
    axes[0,0].hist(df['debt_to_income_ratio'].clip(0, 5), bins=30, alpha=0.7, color='skyblue', edgecolor='black')
    axes[0,0].axvline(x=3, color='red', linestyle='--', label='風險閾值 (3.0)')
    axes[0,0].set_xlabel('債務收入比')
    axes[0,0].set_ylabel('頻數')
    axes[0,0].set_title('單變量分析（債務收入比分佈）')
    axes[0,0].legend()
    axes[0,0].grid(True, alpha=0.3)
    
    # 2. 多變量分析 - 年齡 vs 債務收入比
    scatter = axes[0,1].scatter(df['age'], df['debt_to_income_ratio'].clip(0, 5), 
                               c=df['default'], alpha=0.6, cmap='coolwarm', s=20)
    axes[0,1].axhline(y=3, color='red', linestyle='--', alpha=0.7, label='高負債風險線')
    axes[0,1].axvline(x=30, color='orange', linestyle='--', alpha=0.7, label='年輕風險線')
    axes[0,1].set_xlabel('年齡')
    axes[0,1].set_ylabel('債務收入比')
    axes[0,1].set_title('多變量分析（年齡 vs 債務收入比）')
    axes[0,1].legend()
    plt.colorbar(scatter, ax=axes[0,1], label='是否違約')
    
    # 3. 業務分析 - 穩定性特徵分析
    stability_risk = df.groupby(pd.cut(df['job_stability_score'], bins=5))['default'].mean()
    axes[0,2].bar(range(len(stability_risk)), stability_risk.values, color=['green', 'lightgreen', 'yellow', 'orange', 'red'])
    axes[0,2].set_xlabel('工作穩定性分數分組')
    axes[0,2].set_ylabel('違約率')
    axes[0,2].set_title('業務分析（工作穩定性 vs 違約率）')
    axes[0,2].set_xticks(range(len(stability_risk)))
    axes[0,2].set_xticklabels([f'組{i+1}' for i in range(len(stability_risk))], rotation=45)
    
    # 4. 特徵重要性 - 風險評分對比
    risk_features = ['has_serious_delinquency', 'debt_to_income_ratio', 'young_high_debt_risk', 
                    'job_stability_score', 'is_off_hours_application']
    feature_importance = [0.25, 0.15, 0.15, 0.15, 0.05]
    feature_labels = ['嚴重逾期', '債務收入比', '年輕高負債', '工作穩定性', '非工作時間申請']
    
    axes[1,0].barh(feature_labels, feature_importance, color=['red', 'orange', 'yellow', 'lightgreen', 'lightblue'])
    axes[1,0].set_xlabel('特徵權重')
    axes[1,0].set_title('特徵權重分配（在綜合風險評分中的重要性）')
    for i, v in enumerate(feature_importance):
        axes[1,0].text(v + 0.01, i, f'{v:.2f}', va='center')
    
    # 5. 風險等級分佈
    risk_level_counts = df['risk_level'].value_counts()
    colors = ['lightgreen', 'gold', 'lightcoral']
    wedges, texts, autotexts = axes[1,1].pie(risk_level_counts.values, labels=risk_level_counts.index, 
                                            autopct='%1.1f%%', colors=colors, startangle=90)
    axes[1,1].set_title('風險等級分佈（最終客户風險分層）')
    
    # 6. 模型效果驗證 - 風險評分 vs 違約率
    risk_bins = pd.cut(df['comprehensive_risk_score'], bins=10)
    performance_data = df.groupby(risk_bins)['default'].mean()
    axes[1,2].plot(range(len(performance_data)), performance_data.values, marker='o', linewidth=2, markersize=6)
    axes[1,2].set_xlabel('風險評分分組')
    axes[1,2].set_ylabel('實際違約率')
    axes[1,2].set_title('特徵工程效果驗證（風險評分預測能力）')
    axes[1,2].set_xticks(range(len(performance_data)))
    axes[1,2].set_xticklabels([f'組{i+1}' for i in range(len(performance_data))], rotation=45)
    axes[1,2].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()

# 生成可視化
plot_feature_engineering_analysis(engineered_data)

# 輸出關鍵特徵統計
print("\n關鍵特徵統計:")
feature_stats = engineered_data[['debt_to_income_ratio', 'comprehensive_risk_score', 'default']].describe()
print(feature_stats.round(4))

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime, timedelta
import warnings
warnings.filterwarnings('ignore')

# 設置中文字體
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 60)
print("電商推薦系統特徵工程完整示例")
print("=" * 60)

# 1. 創建電商用户行為數據（簡化版本避免複雜數據類型）
np.random.seed(42)
n_samples = 2000

# 生成基礎數據
raw_behavior_data = pd.DataFrame({
    'user_id': np.random.randint(1, 501, n_samples),
    'item_id': np.random.randint(1, 101, n_samples),
    'behavior_type': np.random.choice(['view', 'click', 'add_cart', 'purchase'], n_samples, 
                                    p=[0.5, 0.3, 0.15, 0.05]),
    'category': np.random.choice(['電子產品', '服裝', '家居', '美妝', '食品'], n_samples),
    'price': np.random.lognormal(5, 1, n_samples).clip(10, 2000),
    'rating': np.random.normal(4.2, 0.5, n_samples).clip(3.0, 5.0),
    'user_value_score': np.random.choice([1, 2, 3, 4], n_samples, p=[0.5, 0.3, 0.15, 0.05]),
    'hour': np.random.randint(0, 24, n_samples)
})

# 生成點擊標籤
def generate_click_label(df):
    click_proba = (
        (df['behavior_type'] == 'click').astype(int) * 0.3 +
        (df['behavior_type'] == 'add_cart').astype(int) * 0.5 +
        (df['price'] < 100).astype(int) * 0.1 +
        (df['user_value_score'] >= 3).astype(int) * 0.1 +
        np.random.normal(0, 0.1, len(df))
    )
    return (click_proba > 0.3).astype(int)

raw_behavior_data['is_click'] = generate_click_label(raw_behavior_data)

def create_recommendation_features(df):
    """創建電商推薦特徵"""
    df_fe = df.copy()
    
    print(">>> 開始電商推薦特徵工程處理...")
    
    # 單變量特徵
    df_fe['is_high_value_user'] = (df_fe['user_value_score'] >= 3).astype(int)
    df_fe['is_premium_product'] = (df_fe['price'] > df_fe['price'].quantile(0.8)).astype(int)
    df_fe['is_high_rating'] = (df_fe['rating'] > 4.5).astype(int)
    df_fe['is_weekend'] = np.random.choice([0, 1], len(df_fe), p=[0.7, 0.3])  # 簡化
    df_fe['is_evening'] = ((df_fe['hour'] >= 18) & (df_fe['hour'] <= 23)).astype(int)
    
    # 多變量特徵
    user_avg_price = df_fe.groupby('user_id')['price'].mean()
    df_fe['user_item_price_affinity'] = df_fe.apply(
        lambda x: x['price'] / user_avg_price.get(x['user_id'], x['price']), axis=1
    ).fillna(1)
    
    # 業務特徵
    df_fe['product_competitiveness'] = (
        (df_fe['rating'] / 5) * 0.4 +
        (1 - (df_fe['price'] / df_fe['price'].max())) * 0.3 +
        df_fe['is_high_rating'] * 0.3
    )
    
    # 綜合評分
    df_fe['recommendation_score'] = (
        df_fe['product_competitiveness'] * 0.4 +
        df_fe['user_item_price_affinity'].clip(0, 2) * 0.3 +
        df_fe['user_value_score'] / 4 * 0.2 +
        df_fe['is_evening'] * 0.1
    )
    
    df_fe['recommendation_level'] = pd.cut(df_fe['recommendation_score'],
                                         bins=[0, 0.3, 0.6, 1.0],
                                         labels=['低推薦度', '中推薦度', '高推薦度'])
    
    print(">>> 電商推薦特徵工程完成！")
    return df_fe

# 執行特徵工程
engineered_data = create_recommendation_features(raw_behavior_data)

# 修復後的可視化函數
def create_recommendation_feature_visualization(df):
    """創建電商推薦特徵工程可視化圖表"""
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))
    fig.suptitle('電商推薦系統特徵工程分析', fontsize=16, fontweight='bold')
    
    # 創建臨時數據副本用於分組
    df_temp = df.copy()
    
    # 1. 用户行為類型分佈
    ax1 = axes[0, 0]
    behavior_counts = df['behavior_type'].value_counts()
    colors = ['lightblue', 'lightgreen', 'orange', 'red']
    bars = ax1.bar(behavior_counts.index, behavior_counts.values, color=colors, alpha=0.7)
    ax1.set_title('用户行為類型分佈 (單變量分析)')
    ax1.set_xlabel('行為類型')
    ax1.set_ylabel('頻次')
    for i, v in enumerate(behavior_counts.values):
        ax1.text(i, v + 20, str(v), ha='center')
    
    # 2. 商品類別點擊率分析
    ax2 = axes[0, 1]
    category_clicks = df.groupby('category')['is_click'].mean().sort_values(ascending=False)
    bars = ax2.bar(range(len(category_clicks)), category_clicks.values, color='skyblue', alpha=0.7)
    ax2.set_title('商品類別點擊率分析 (單變量分析)')
    ax2.set_xlabel('商品類別')
    ax2.set_ylabel('點擊率')
    ax2.set_xticks(range(len(category_clicks)))
    ax2.set_xticklabels(category_clicks.index, rotation=45)
    for i, v in enumerate(category_clicks.values):
        ax2.text(i, v + 0.01, f'{v:.3f}', ha='center')
    
    # 3. 時間段點擊率分析
    ax3 = axes[0, 2]
    hourly_click_rate = df.groupby('hour')['is_click'].mean()
    ax3.plot(hourly_click_rate.index, hourly_click_rate.values, marker='o', color='purple', linewidth=2)
    ax3.set_title('時間段點擊率分析 (單變量分析)')
    ax3.set_xlabel('小時')
    ax3.set_ylabel('點擊率')
    ax3.grid(True, alpha=0.3)
    ax3.axvspan(18, 23, alpha=0.2, color='orange', label='高峯時段')
    ax3.legend()
    
    # 4. 價格區間點擊率分析
    ax4 = axes[1, 0]
    # 修復：使用數值索引避免類型轉換錯誤
    price_bins = [0, 100, 300, 1000, 2000]
    price_labels = ['0-100', '100-300', '300-1000', '1000+']
    df_temp['price_group'] = pd.cut(df_temp['price'], bins=price_bins, labels=price_labels)
    price_click_rate = df_temp.groupby('price_group')['is_click'].mean()
    bars = ax4.bar(range(len(price_click_rate)), price_click_rate.values, 
                  color=['green', 'lightgreen', 'orange', 'red'])
    ax4.set_title('價格區間點擊率分析 (多變量分析)')
    ax4.set_xlabel('價格區間')
    ax4.set_ylabel('點擊率')
    ax4.set_xticks(range(len(price_click_rate)))
    ax4.set_xticklabels(price_click_rate.index, rotation=45)
    
    # 5. 用户價值 vs 商品價格
    ax5 = axes[1, 1]
    scatter = ax5.scatter(df['user_value_score'], df['price'], 
                         c=df['is_click'], alpha=0.6, cmap='coolwarm', s=20)
    ax5.set_title('用户價值 vs 商品價格 (多變量分析)')
    ax5.set_xlabel('用户價值評分')
    ax5.set_ylabel('商品價格')
    plt.colorbar(scatter, ax=ax5, label='是否點擊')
    
    # 6. 推薦等級分佈
    ax6 = axes[1, 2]
    level_counts = df['recommendation_level'].value_counts()
    colors = ['lightgreen', 'gold', 'lightcoral']
    wedges, texts, autotexts = ax6.pie(level_counts.values, 
                                       labels=level_counts.index, 
                                       autopct='%1.1f%%', colors=colors, startangle=90)
    ax6.set_title('推薦等級分佈 (業務應用)')
    
    plt.tight_layout()
    plt.savefig('電商推薦特徵工程.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    # 清理臨時列
    if 'price_group' in df_temp.columns:
        df_temp.drop('price_group', axis=1, inplace=True)
    
    return fig

# 生成第二張圖：特徵組合與效果驗證
def create_feature_validation_plot(df):
    """創建特徵組合與效果驗證圖表"""
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))
    fig.suptitle('電商推薦特徵組合與效果驗證', fontsize=16, fontweight='bold')
    
    df_temp = df.copy()
    
    # 1. 用户價格偏好匹配度
    ax1 = axes[0, 0]
    # 修復：使用數值索引
    affinity_bins = 5
    df_temp['affinity_group'] = pd.cut(df_temp['user_item_price_affinity'], bins=affinity_bins)
    affinity_click_rate = df_temp.groupby('affinity_group')['is_click'].mean()
    ax1.plot(range(len(affinity_click_rate)), affinity_click_rate.values, 
            marker='o', color='darkblue', linewidth=2)
    ax1.set_title('用户價格偏好匹配度 (特徵組合)')
    ax1.set_xlabel('價格偏好匹配度分組')
    ax1.set_ylabel('點擊率')
    ax1.grid(True, alpha=0.3)
    
    # 2. 商品競爭力分析
    ax2 = axes[0, 1]
    competitiveness_bins = 5
    df_temp['competitiveness_group'] = pd.cut(df_temp['product_competitiveness'], bins=competitiveness_bins)
    competitiveness_click_rate = df_temp.groupby('competitiveness_group')['is_click'].mean()
    bars = ax2.bar(range(len(competitiveness_click_rate)), competitiveness_click_rate.values, 
                  color='green', alpha=0.7)
    ax2.set_title('商品競爭力 vs 點擊率 (特徵組合)')
    ax2.set_xlabel('商品競爭力分組')
    ax2.set_ylabel('點擊率')
    
    # 3. 用户價值分層效果
    ax3 = axes[0, 2]
    user_value_click_rate = df.groupby('user_value_score')['is_click'].mean()
    bars = ax3.bar(user_value_click_rate.index, user_value_click_rate.values, 
                  color=['lightblue', 'lightgreen', 'orange', 'red'])
    ax3.set_title('用户價值分層效果驗證 (業務特徵)')
    ax3.set_xlabel('用户價值評分')
    ax3.set_ylabel('點擊率')
    for i, v in enumerate(user_value_click_rate.values):
        ax3.text(i+1, v + 0.01, f'{v:.3f}', ha='center')
    
    # 4. 時間場景分析
    ax4 = axes[1, 0]
    time_scenario = df.groupby('is_evening')['is_click'].mean()
    bars = ax4.bar(['白天', '晚上'], time_scenario.values, color=['lightblue', 'lightgreen'])
    ax4.set_title('時間場景點擊率分析 (業務特徵)')
    ax4.set_xlabel('時間場景')
    ax4.set_ylabel('點擊率')
    
    # 5. 綜合評分效果驗證
    ax5 = axes[1, 1]
    score_bins = 6
    df_temp['score_group'] = pd.cut(df_temp['recommendation_score'], bins=score_bins)
    score_performance = df_temp.groupby('score_group')['is_click'].mean()
    ax5.plot(range(len(score_performance)), score_performance.values, 
            marker='s', linewidth=2, markersize=6, color='red')
    ax5.set_title('推薦評分預測能力驗證 (效果驗證)')
    ax5.set_xlabel('推薦評分分組')
    ax5.set_ylabel('實際點擊率')
    ax5.grid(True, alpha=0.3)
    
    # 6. 特徵重要性總結
    ax6 = axes[1, 2]
    features_importance = {
        '商品競爭力': 0.4,
        '價格匹配度': 0.3,
        '用户價值': 0.2,
        '時間場景': 0.1
    }
    ax6.barh(list(features_importance.keys()), list(features_importance.values()), 
             color=['red', 'orange', 'yellow', 'lightgreen'])
    ax6.set_title('特徵權重分配 (模型可解釋性)')
    ax6.set_xlabel('特徵權重')
    ax6.set_xlim(0, 0.5)
    for i, (k, v) in enumerate(features_importance.items()):
        ax6.text(v + 0.01, i, f'{v:.2f}', va='center')
    
    plt.tight_layout()
    plt.savefig('電商推薦特徵驗證.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    # 清理臨時列
    temp_cols = ['affinity_group', 'competitiveness_group', 'score_group']
    for col in temp_cols:
        if col in df_temp.columns:
            df_temp.drop(col, axis=1, inplace=True)
    
    return fig

# 生成兩張圖表
print("\n>>> 生成第一張圖：特徵發現與分析過程...")
create_recommendation_feature_visualization(engineered_data)

print("\n>>> 生成第二張圖：特徵組合與效果驗證...")
create_feature_validation_plot(engineered_data)

# 輸出統計信息
print("\n" + "="*50)
print("電商推薦特徵工程統計摘要")
print("="*50)

print(f"用户行為記錄數量: {len(engineered_data):,}")
print(f"整體點擊率: {engineered_data['is_click'].mean():.2%}")
print(f"原始特徵數量: {len(raw_behavior_data.columns)}")
print(f"衍生特徵數量: {len(engineered_data.columns) - len(raw_behavior_data.columns)}")

print(f"\n推薦等級分佈:")
level_dist = engineered_data['recommendation_level'].value_counts()
for level, count in level_dist.items():
    percentage = count / len(engineered_data) * 100
    print(f"  {level}: {count:,}條 ({percentage:.1f}%)")

print(f"\n關鍵指標統計:")
print(f"  高價值用户比例: {engineered_data['is_high_value_user'].mean():.2%}")
print(f"  高評分商品比例: {engineered_data['is_high_rating'].mean():.2%}")
print(f"  平均推薦評分: {engineered_data['recommendation_score'].mean():.3f}")

print(f"\n>>> 電商推薦特徵工程分析完成！")

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')

# 設置中文字體
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 60)
print("糖尿病視網膜病變特徵工程完整示例")
print("=" * 60)

# 1. 創建醫療數據（使用之前的代碼）
np.random.seed(42)
n_samples = 800

raw_medical_data = pd.DataFrame({
    'patient_id': range(1, n_samples + 1),
    'age': np.random.randint(30, 80, n_samples),
    'diabetes_years': np.random.exponential(8, n_samples).clip(1, 30),
    'hbA1c': np.random.normal(7.5, 1.5, n_samples).clip(5.0, 12.0),
    'systolic_bp': np.random.normal(140, 20, n_samples).clip(90, 200),
    'diastolic_bp': np.random.normal(85, 15, n_samples).clip(60, 120),
    'bmi': np.random.normal(28, 5, n_samples).clip(18, 45),
    'cholesterol': np.random.normal(200, 40, n_samples).clip(150, 300),
    'smoking_status': np.random.choice(['不吸煙', '已戒煙', '吸煙'], n_samples, p=[0.4, 0.3, 0.3]),
    'family_history': np.random.binomial(1, 0.3, n_samples),
    'microaneurysms_count': np.random.poisson(3, n_samples),
    'hemorrhage_area': np.random.exponential(0.5, n_samples).clip(0, 5),
    'exudates_area': np.random.exponential(0.3, n_samples).clip(0, 4),
    'vessel_tortuosity': np.random.normal(1.5, 0.4, n_samples).clip(1.0, 3.0),
    'macular_thickness': np.random.normal(280, 30, n_samples).clip(200, 400),
    'optic_disc_ratio': np.random.normal(0.4, 0.1, n_samples).clip(0.2, 0.6)
})

def generate_dr_label(df):
    dr_proba = (
        (df['diabetes_years'] > 10) * 0.3 +
        (df['hbA1c'] > 8.0) * 0.25 +
        (df['microaneurysms_count'] > 5) * 0.15 +
        (df['hemorrhage_area'] > 1.0) * 0.15 +
        (df['exudates_area'] > 0.5) * 0.1 +
        (df['systolic_bp'] > 160) * 0.05 +
        np.random.normal(0, 0.05, len(df))
    )
    return (dr_proba > 0.4).astype(int)

raw_medical_data['has_dr'] = generate_dr_label(raw_medical_data)

def create_medical_features(df):
    df_fe = df.copy()
    
    print(">>> 開始醫療特徵工程處理...")
    
    # 單變量特徵
    df_fe['high_hba1c_risk'] = (df_fe['hbA1c'] > 7.0).astype(int)
    df_fe['long_diabetes_risk'] = (df_fe['diabetes_years'] > 10).astype(int)
    df_fe['hypertension_risk'] = ((df_fe['systolic_bp'] > 140) | (df_fe['diastolic_bp'] > 90)).astype(int)
    df_fe['significant_microaneurysms'] = (df_fe['microaneurysms_count'] > 5).astype(int)
    df_fe['significant_hemorrhage'] = (df_fe['hemorrhage_area'] > 1.0).astype(int)
    
    # 多變量特徵
    df_fe['metabolic_syndrome_score'] = (
        (df_fe['hbA1c'] > 7.0).astype(int) * 0.3 +
        (df_fe['bmi'] > 30).astype(int) * 0.2 +
        (df_fe['systolic_bp'] > 140).astype(int) * 0.2 +
        (df_fe['cholesterol'] > 240).astype(int) * 0.2 +
        (df_fe['smoking_status'] == '吸煙').astype(int) * 0.1
    )
    
    df_fe['lesion_severity_score'] = (
        np.minimum(df_fe['microaneurysms_count'] / 10, 1) * 0.3 +
        np.minimum(df_fe['hemorrhage_area'] / 3, 1) * 0.4 +
        np.minimum(df_fe['exudates_area'] / 2, 1) * 0.3
    )
    
    # 業務特徵
    df_fe['dr_mild'] = ((df_fe['microaneurysms_count'] >= 1) & (df_fe['microaneurysms_count'] <= 5)).astype(int)
    df_fe['dr_moderate'] = ((df_fe['microaneurysms_count'] > 5) | (df_fe['hemorrhage_area'] >= 0.5)).astype(int)
    df_fe['dr_severe'] = ((df_fe['microaneurysms_count'] > 15) | (df_fe['hemorrhage_area'] >= 2.0)).astype(int)
    
    df_fe['macular_edema_risk'] = (
        (df_fe['macular_thickness'] > 300) * 0.6 +
        (df_fe['exudates_area'] > 0.3) * 0.4
    )
    
    # 綜合評分
    df_fe['comprehensive_dr_risk_score'] = (
        df_fe['metabolic_syndrome_score'] * 0.3 +
        df_fe['lesion_severity_score'] * 0.4 +
        df_fe['macular_edema_risk'] * 0.3
    )
    
    df_fe['recommended_referral'] = (df_fe['comprehensive_dr_risk_score'] > 0.5).astype(int)
    df_fe['risk_level'] = pd.cut(df_fe['comprehensive_dr_risk_score'],
                               bins=[0, 0.3, 0.6, 1.0],
                               labels=['低風險', '中風險', '高風險'])
    
    print(">>> 醫療特徵工程完成！")
    return df_fe

# 執行特徵工程
engineered_medical_data = create_medical_features(raw_medical_data)

# 第一張圖：特徵發現與分析過程
def create_feature_discovery_plot(df):
    """創建特徵發現與分析過程圖表"""
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))
    fig.suptitle('糖尿病視網膜病變 - 特徵發現與分析過程', fontsize=16, fontweight='bold')
    
    # 1. 單變量分析 - HbA1c分佈
    ax1 = axes[0, 0]
    sns.histplot(data=df, x='hbA1c', hue='has_dr', bins=30, ax=ax1, alpha=0.7)
    ax1.axvline(x=7.0, color='red', linestyle='--', linewidth=2, label='臨牀閾值 (7.0%)')
    ax1.set_title('HbA1c分佈與風險閾值 (單變量分析)', fontsize=12)
    ax1.set_xlabel('糖化血紅蛋白 (%)')
    ax1.set_ylabel('患者數量')
    ax1.legend()
    
    # 2. 單變量分析 - 糖尿病病程風險
    ax2 = axes[0, 1]
    years_bins = [0, 5, 10, 15, 30]
    years_labels = ['0-5年', '5-10年', '10-15年', '15+年']
    df_temp = df.copy()
    df_temp['diabetes_years_group'] = pd.cut(df_temp['diabetes_years'], bins=years_bins, labels=years_labels)
    years_risk = df_temp.groupby('diabetes_years_group')['has_dr'].mean()
    years_risk.plot(kind='bar', ax=ax2, color=['lightgreen', 'yellow', 'orange', 'red'])
    ax2.set_title('糖尿病病程 vs 患病率 (單變量分析)', fontsize=12)
    ax2.set_xlabel('糖尿病病程')
    ax2.set_ylabel('視網膜病變患病率')
    
    # 3. 多變量分析 - 代謝指標交互
    ax3 = axes[0, 2]
    scatter = ax3.scatter(df['hbA1c'], df['diabetes_years'], 
                         c=df['has_dr'], alpha=0.6, cmap='coolwarm', s=30)
    ax3.set_title('HbA1c vs 糖尿病病程 (多變量分析)', fontsize=12)
    ax3.set_xlabel('糖化血紅蛋白 (%)')
    ax3.set_ylabel('糖尿病病程 (年)')
    ax3.axvline(x=7.0, color='red', linestyle='--', alpha=0.5)
    ax3.axhline(y=10, color='red', linestyle='--', alpha=0.5)
    plt.colorbar(scatter, ax=ax3, label='是否患病')
    
    # 4. 多變量分析 - 病變特徵相關性
    ax4 = axes[1, 0]
    lesion_features = ['microaneurysms_count', 'hemorrhage_area', 'exudates_area', 'vessel_tortuosity']
    lesion_corr = df[lesion_features + ['has_dr']].corr()
    sns.heatmap(lesion_corr, annot=True, cmap='RdYlBu_r', center=0, ax=ax4, fmt='.2f')
    ax4.set_title('影像病變特徵相關性 (多變量分析)', fontsize=12)
    
    # 5. 業務特徵 - 臨牀分期分佈
    ax5 = axes[1, 1]
    dr_stages = ['dr_mild', 'dr_moderate', 'dr_severe']
    dr_counts = [df[stage].sum() for stage in dr_stages]
    colors = ['lightblue', 'orange', 'red']
    ax5.bar(['輕度', '中度', '重度'], dr_counts, color=colors, alpha=0.7)
    ax5.set_title('臨牀分期分佈 (業務特徵)', fontsize=12)
    ax5.set_ylabel('患者數量')
    for i, v in enumerate(dr_counts):
        ax5.text(i, v + 5, str(v), ha='center')
    
    # 6. 業務特徵 - 黃斑水腫風險
    ax6 = axes[1, 2]
    edema_bins = [0, 0.3, 0.6, 1.0]
    edema_labels = ['低風險', '中風險', '高風險']
    df_temp['edema_risk_group'] = pd.cut(df_temp['macular_edema_risk'], bins=edema_bins, labels=edema_labels)
    edema_risk_analysis = df_temp.groupby('edema_risk_group')['has_dr'].mean()
    edema_risk_analysis.plot(kind='bar', ax=ax6, color=['green', 'yellow', 'red'])
    ax6.set_title('黃斑水腫風險 vs 患病率 (業務特徵)', fontsize=12)
    ax6.set_xlabel('黃斑水腫風險等級')
    ax6.set_ylabel('視網膜病變患病率')
    
    plt.tight_layout()
    plt.savefig('糖尿病視網膜病變_特徵發現.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    return fig

# 第二張圖：特徵組合與效果驗證
def create_feature_validation_plot(df):
    """創建特徵組合與效果驗證圖表"""
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))
    fig.suptitle('糖尿病視網膜病變 - 特徵組合與效果驗證', fontsize=16, fontweight='bold')
    
    # 1. 特徵組合 - 代謝綜合徵評分
    ax1 = axes[0, 0]
    metabolic_bins = 5
    df_temp = df.copy()
    df_temp['metabolic_group'] = pd.cut(df_temp['metabolic_syndrome_score'], bins=metabolic_bins)
    metabolic_risk = df_temp.groupby('metabolic_group')['has_dr'].mean()
    ax1.plot(range(len(metabolic_risk)), metabolic_risk.values, marker='o', 
             color='purple', linewidth=2, markersize=6)
    ax1.set_title('代謝綜合徵評分 vs 患病率 (特徵組合)', fontsize=12)
    ax1.set_xlabel('代謝評分分組')
    ax1.set_ylabel('視網膜病變患病率')
    ax1.grid(True, alpha=0.3)
    
    # 2. 特徵組合 - 病變嚴重度評分
    ax2 = axes[0, 1]
    lesion_bins = 5
    df_temp['lesion_group'] = pd.cut(df_temp['lesion_severity_score'], bins=lesion_bins)
    lesion_risk = df_temp.groupby('lesion_group')['has_dr'].mean()
    ax2.plot(range(len(lesion_risk)), lesion_risk.values, marker='s', 
             color='brown', linewidth=2, markersize=6)
    ax2.set_title('病變嚴重度評分 vs 患病率 (特徵組合)', fontsize=12)
    ax2.set_xlabel('病變嚴重度分組')
    ax2.set_ylabel('視網膜病變患病率')
    ax2.grid(True, alpha=0.3)
    
    # 3. 綜合評分效果驗證
    ax3 = axes[0, 2]
    risk_bins = 8
    df_temp['risk_score_group'] = pd.cut(df_temp['comprehensive_dr_risk_score'], bins=risk_bins)
    score_performance = df_temp.groupby('risk_score_group')['has_dr'].mean()
    ax3.plot(range(len(score_performance)), score_performance.values, 
            marker='o', linewidth=2, markersize=6, color='darkred')
    ax3.set_title('綜合風險評分預測能力 (效果驗證)', fontsize=12)
    ax3.set_xlabel('風險評分分組')
    ax3.set_ylabel('實際患病率')
    ax3.grid(True, alpha=0.3)
    
    # 4. 最終風險分層
    ax4 = axes[1, 0]
    risk_level_counts = df['risk_level'].value_counts()
    colors = ['lightgreen', 'gold', 'lightcoral']
    wedges, texts, autotexts = ax4.pie(risk_level_counts.values, 
                                       labels=risk_level_counts.index, 
                                       autopct='%1.1f%%', colors=colors, startangle=90)
    ax4.set_title('最終風險等級分佈 (臨牀決策支持)', fontsize=12)
    
    # 5. 轉診建議分析
    ax5 = axes[1, 1]
    referral_analysis = df.groupby('recommended_referral')['has_dr'].mean()
    referral_analysis.plot(kind='bar', ax=ax5, color=['lightblue', 'red'], alpha=0.7)
    ax5.set_title('轉診建議效果驗證 (業務應用)', fontsize=12)
    ax5.set_xlabel('是否建議轉診')
    ax5.set_ylabel('實際患病率')
    ax5.set_xticklabels(['否', '是'], rotation=0)
    
    # 6. 特徵重要性總結
    ax6 = axes[1, 2]
    features_importance = {
        '代謝綜合徵': 0.3,
        '病變嚴重度': 0.4,
        '黃斑水腫': 0.3
    }
    ax6.barh(list(features_importance.keys()), list(features_importance.values()), 
             color=['red', 'orange', 'lightgreen'])
    ax6.set_title('特徵權重分配 (模型可解釋性)', fontsize=12)
    ax6.set_xlabel('特徵權重')
    ax6.set_xlim(0, 0.5)
    for i, (k, v) in enumerate(features_importance.items()):
        ax6.text(v + 0.01, i, f'{v:.2f}', va='center')
    
    plt.tight_layout()
    plt.savefig('糖尿病視網膜病變_特徵驗證.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    return fig

# 生成兩張圖表
print("\n>>> 生成第一張圖：特徵發現與分析過程...")
create_feature_discovery_plot(engineered_medical_data)

print("\n>>> 生成第二張圖：特徵組合與效果驗證...")
create_feature_validation_plot(engineered_medical_data)

# 輸出統計信息
print("\n" + "="*50)
print("醫療特徵工程統計摘要")
print("="*50)

print(f"患者樣本數量: {len(engineered_medical_data):,}")
print(f"糖尿病視網膜病變患病率: {engineered_medical_data['has_dr'].mean():.2%}")
print(f"原始特徵數量: {len(raw_medical_data.columns)}")
print(f"衍生特徵數量: {len(engineered_medical_data.columns) - len(raw_medical_data.columns)}")

print(f"\n風險等級分佈:")
risk_dist = engineered_medical_data['risk_level'].value_counts()
for level, count in risk_dist.items():
    percentage = count / len(engineered_medical_data) * 100
    print(f"  {level}: {count:,}人 ({percentage:.1f}%)")

print(f"\n>>> 醫療特徵工程分析完成！")