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>>> 医疗特征工程分析完成！")