# -*- coding: utf-8 -*-
# BP神经网络信用风险实证分析 - 使用平衡模拟数据

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import time
from datetime import datetime
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_split, cross_val_score, GridSearchCV
from sklearn.metrics import confusion_matrix, roc_curve, auc, accuracy_score, precision_score, recall_score, f1_score, roc_auc_score
from sklearn.neural_network import MLPClassifier
import warnings
warnings.filterwarnings('ignore')

# 记录开始时间
start_time = datetime.now()
print("="*70)
print(f"【运行日志】BP神经网络信用风险实证分析")
print(f"开始时间：{start_time.strftime('%Y-%m-%d %H:%M:%S')}")
print("="*70)

# ============================================================
# Step 1: 生成平衡的模拟数据
# ============================================================
print("\n[Step 1] 生成模拟数据...")
step_start = time.time()

np.random.seed(42)
n_samples = 5000  # 总样本量
n_features = 14

# 生成特征
X = pd.DataFrame(np.random.randn(n_samples, n_features),
                 columns=['loan_amnt', 'term', 'int_rate', 'installment', 'annual_inc',
                          'dti', 'delinq_2yrs', 'inq_last_6mths', 'open_acc', 'pub_rec',
                          'revol_bal', 'revol_util', 'total_acc', 'mort_acc'])

# 构造违约概率（与特征1、2、3相关）
prob = 1 / (1 + np.exp(-(X.iloc[:,0] * 0.8 + X.iloc[:,1] * 0.5 + X.iloc[:,2] * 0.3)))
# 设定违约率为40%左右
y = (prob > np.percentile(prob, 60)).astype(int)

print(f"生成{n_samples}条模拟数据")
print(f"违约样本：{sum(y)}，正常样本：{len(y)-sum(y)}")
print(f"违约率：{y.mean():.2%}")

step_time = time.time() - step_start
print(f"[Step 1] 完成，耗时：{step_time:.2f}秒")

# ============================================================
# Step 2: 标准化
# ============================================================
print("\n[Step 2] 数据标准化...")
step_start = time.time()

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

print(f"✓ 标准化完成")
step_time = time.time() - step_start
print(f"[Step 2] 完成，耗时：{step_time:.2f}秒")

# ============================================================
# Step 3: PCA降维
# ============================================================
print("\n[Step 3] PCA降维...")
step_start = time.time()

pca = PCA()
pca.fit(X_scaled)
cumsum = np.cumsum(pca.explained_variance_ratio_)
n_comp = np.argmax(cumsum >= 0.85) + 1
print(f"累积方差贡献率：{cumsum[n_comp-1]:.4f}")
print(f"选择主成分数：{n_comp}")

pca = PCA(n_components=n_comp)
X_pca = pca.fit_transform(X_scaled)

step_time = time.time() - step_start
print(f"[Step 3] 完成，耗时：{step_time:.2f}秒")

# ============================================================
# Step 4: 划分数据集
# ============================================================
print("\n[Step 4] 划分训练集和测试集...")
step_start = time.time()

X_train, X_test, y_train, y_test = train_test_split(X_pca, y, test_size=0.2, random_state=42, stratify=y)

print(f"训练集：{len(X_train)}条，测试集：{len(X_test)}条")
step_time = time.time() - step_start
print(f"[Step 4] 完成，耗时：{step_time:.2f}秒")

# ============================================================
# Step 5: BP神经网络
# ============================================================
print("\n[Step 5] BP神经网络训练...")
step_start = time.time()

# 参数搜索
params = {
    'hidden_layer_sizes': [(10,), (15,), (10,5)],
    'activation': ['relu', 'tanh'],
    'learning_rate_init': [0.001, 0.01],
    'alpha': [0.0001, 0.001]
}

mlp = MLPClassifier(max_iter=500, random_state=42, early_stopping=True)
grid = GridSearchCV(mlp, params, cv=5, scoring='roc_auc', n_jobs=-1)
grid.fit(X_train, y_train)

print(f"最佳参数：{grid.best_params_}")
print(f"最佳交叉验证AUC：{grid.best_score_:.4f}")

best_model = grid.best_estimator_

step_time = time.time() - step_start
print(f"[Step 5] 完成，耗时：{step_time:.2f}秒")

# ============================================================
# Step 6: 模型评估
# ============================================================
print("\n[Step 6] 模型评估...")
step_start = time.time()

# 预测
y_pred = best_model.predict(X_test)
y_proba = best_model.predict_proba(X_test)[:, 1]

# 计算指标
acc = accuracy_score(y_test, y_pred)
prec = precision_score(y_test, y_pred)
rec = recall_score(y_test, y_pred)
f1 = f1_score(y_test, y_pred)
auc_score = roc_auc_score(y_test, y_proba)

# 交叉验证
cv_scores = cross_val_score(best_model, X_train, y_train, cv=5, scoring='accuracy')

# 混淆矩阵
cm = confusion_matrix(y_test, y_pred)
tn, fp, fn, tp = cm.ravel()

print(f"\n【模型性能指标】")
print(f"准确率：{acc:.4f}")
print(f"精确率：{prec:.4f}")
print(f"召回率：{rec:.4f}")
print(f"F1分数：{f1:.4f}")
print(f"AUC值：{auc_score:.4f}")
print(f"5折CV准确率：{cv_scores.mean():.4f} ± {cv_scores.std():.4f}")

print(f"\n【混淆矩阵】")
print(f"           预测正常  预测违约")
print(f"实际正常      {tn:>6}    {fp:>6}")
print(f"实际违约      {fn:>6}    {tp:>6}")

step_time = time.time() - step_start
print(f"[Step 6] 完成，耗时：{step_time:.2f}秒")

# ============================================================
# Step 7: 绘图
# ============================================================
print("\n[Step 7] 生成图表...")
step_start = time.time()

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 混淆矩阵
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', ax=axes[0,0],
            xticklabels=['正常','违约'], yticklabels=['正常','违约'])
axes[0,0].set_title('混淆矩阵')

# ROC曲线
fpr, tpr, _ = roc_curve(y_test, y_proba)
axes[0,1].plot(fpr, tpr, 'b-', lw=2, label=f'AUC={auc_score:.4f}')
axes[0,1].plot([0,1], [0,1], 'r--')
axes[0,1].set_xlabel('假正率')
axes[0,1].set_ylabel('真正率')
axes[0,1].set_title('ROC曲线')
axes[0,1].legend()

# 损失曲线
if hasattr(best_model, 'loss_curve_') and len(best_model.loss_curve_) > 0:
    axes[1,0].plot(best_model.loss_curve_, 'g-', lw=2)
    axes[1,0].set_xlabel('迭代次数')
    axes[1,0].set_ylabel('损失')
    axes[1,0].set_title('训练损失曲线')
    axes[1,0].grid(True, alpha=0.3)

# PCA方差
x = np.arange(1, len(cumsum)+1)
axes[1,1].bar(x[:n_comp], pca.explained_variance_ratio_, alpha=0.7)
axes[1,1].plot(x, cumsum, 'ro-', label='累积方差')
axes[1,1].axhline(y=0.85, color='g', linestyle='--', label='85%阈值')
axes[1,1].set_xlabel('主成分')
axes[1,1].set_ylabel('方差贡献率')
axes[1,1].set_title('PCA方差贡献率')
axes[1,1].legend()

plt.tight_layout()
plt.savefig('模型评估结果.png', dpi=200)
print("✓ 图表已保存")

step_time = time.time() - step_start
print(f"[Step 7] 完成，耗时：{step_time:.2f}秒")

# ============================================================
# 运行日志汇总
# ============================================================
end_time = datetime.now()
total_time = (end_time - start_time).total_seconds()

print("\n" + "="*70)
print("【运行完成】")
print("="*70)
print(f"""
结束时间：{end_time.strftime('%Y-%m-%d %H:%M:%S')}
总运行时间：{total_time:.2f}秒 ({total_time/60:.2f}分钟)

数据规模：{X.shape[0]}条，{X.shape[1]}个特征
PCA维度：{n_comp}个主成分

模型性能：
- 准确率：{acc:.4f}
- 精确率：{prec:.4f}
- 召回率：{rec:.4f}
- F1分数：{f1:.4f}
- AUC值：{auc_score:.4f}
""")

# 保存日志
log = f"""BP神经网络信用风险分析运行日志
生成时间：{end_time.strftime('%Y-%m-%d %H:%M:%S')}

运行时间：{total_time:.2f}秒
样本量：{X.shape[0]}
特征数：{X.shape[1]}
PCA维度：{n_comp}

准确率：{acc:.4f}
精确率：{prec:.4f}
召回率：{rec:.4f}
F1分数：{f1:.4f}
AUC值：{auc_score:.4f}

混淆矩阵：
           预测正常  预测违约
实际正常      {tn}      {fp}
实际违约      {fn}      {tp}
"""

with open('运行日志.txt', 'w', encoding='utf-8') as f:
    f.write(log)
print("✓ 运行日志已保存")

print("\n【全部完成！】")

======================================================================
【运行日志】BP神经网络信用风险实证分析
开始时间：2026-06-27 12:10:52
======================================================================

[Step 1] 生成模拟数据...
生成5000条模拟数据
违约样本：2000，正常样本：3000
违约率：40.00%
[Step 1] 完成，耗时：0.01秒

[Step 2] 数据标准化...
✓ 标准化完成
[Step 2] 完成，耗时：0.02秒

[Step 3] PCA降维...
累积方差贡献率：0.8678
选择主成分数：12
[Step 3] 完成，耗时：0.02秒

[Step 4] 划分训练集和测试集...
训练集：4000条，测试集：1000条
[Step 4] 完成，耗时：0.01秒

[Step 5] BP神经网络训练...
最佳参数：{'activation': 'tanh', 'alpha': 0.001, 'hidden_layer_sizes': (15,), 'learning_rate_init': 0.01}
最佳交叉验证AUC：0.9605
[Step 5] 完成，耗时：14.76秒

[Step 6] 模型评估...

【模型性能指标】
准确率：0.8970
精确率：0.8722
召回率：0.8700
F1分数：0.8711
AUC值：0.9625
5折CV准确率：0.8875 ± 0.0108

【混淆矩阵】
           预测正常  预测违约
实际正常         549        51
实际违约          52       348
[Step 6] 完成，耗时：0.69秒

[Step 7] 生成图表...
✓ 图表已保存
[Step 7] 完成，耗时：1.08秒

======================================================================
【运行完成】
======================================================================

结束时间：2026-06-27 12:11:08
总运行时间：16.61秒 (0.28分钟)

数据规模：5000条，14个特征
PCA维度：12个主成分

模型性能：
- 准确率：0.8970
- 精确率：0.8722
- 召回率：0.8700
- F1分数：0.8711
- AUC值：0.9625

✓ 运行日志已保存

【全部完成！】