AI可解释性：黑盒不够用了之后

import shap
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import train_test_split

# 加载数据和模型
df = pd.read_csv('customer_data.csv')
X = df.drop(['churn'], axis=1)
y = df['churn']

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

# 训练模型
model = GradientBoostingClassifier(n_estimators=100, learning_rate=0.1, max_depth=3, random_state=42)
model.fit(X_train, y_train)

# 计算全局特征重要性
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)

# 可视化特征重要性
plt.figure(figsize=(10, 8))
shap.summary_plot(shap_values, X_test, plot_type="bar")
plt.savefig('global_feature_importance.png', bbox_inches='tight')
plt.close()

# 选取具体客户进行解释
customer_idx = 42  # 那个引起争议的老客户
customer_data = X_test.iloc[[customer_idx]]

# 计算局部 SHAP 值
local_shap_values = explainer.shap_values(customer_data)

# 瀑布图展示
plt.figure(figsize=(12, 8))
shap.waterfall_plot(
    shap.Explanation(
        values=local_shap_values[0],
        base_values=explainer.expected_value[0],
        data=customer_data.iloc[0],
        feature_names=customer_data.columns.tolist()
    )
)
plt.savefig('local_explanation_customer_42.png', bbox_inches='tight')
plt.close()

# 也可以用决策图
plt.figure(figsize=(12, 8))
shap.decision_plot(
    explainer.expected_value[0],
    local_shap_values[0],
    feature_names=customer_data.columns.tolist()
)
plt.savefig('decision_plot_customer_42.png', bbox_inches='tight')
plt.close()

# 生成反事实解释
def generate_counterfactual(explainer, model, instance, target_class, max_iter=100, tolerance=0.5):
    """
    生成反事实样本
    """
    current_pred = model.predict_proba(instance)[0, 1]  # 预测为正类的概率
    original_values = instance.copy().values[0]
    modified_values = original_values.copy()

    # 对每个特征进行尝试调整
    feature_names = instance.columns.tolist()
    n_features = len(feature_names)

    for feature_idx in range(n_features):
        feature_name = feature_names[feature_idx]
        current_value = modified_values[feature_idx]

        # 尝试调整这个特征
        step_size = 0.1 * current_value if current_value != 0 else 0.1

        for direction in [-1, 1]:
            test_instance = modified_values.copy()
            test_instance[feature_idx] = current_value + direction * step_size

            # 预测新概率
            test_df = pd.DataFrame([test_instance], columns=feature_names)
            new_pred = model.predict_proba(test_df)[0, 1]

            # 如果改变了预测结果
            if (target_class == 0 and new_pred < current_pred - tolerance) or \
               (target_class == 1 and new_pred > current_pred + tolerance):
                modified_values[feature_idx] = test_instance[feature_idx]
                current_pred = new_pred
                break

    # 创建结果 DataFrame
    result_df = pd.DataFrame({
        'feature': feature_names,
        'original': original_values,
        'counterfactual': modified_values,
        'change': modified_values - original_values
    })
    # 只返回有变化的特征
    changed_features = result_df[result_df['change'] != 0]
    changed_features = changed_features.sort_values('change', key=abs, ascending=False)

    return changed_features

# 为客户 42 生成反事实解释
counterfactuals = generate_counterfactual(explainer, model, customer_data, target_class=0)
print("降低风险的建议调整：")
print(counterfactuals.head(10))

# 方案 1：采样计算，不要对所有样本都算
sample_size = min(100, len(X_test))
sample_indices = np.random.choice(len(X_test), sample_size, replace=False)
shap_values = explainer.shap_values(X_test.iloc[sample_indices])

# 方案 2：使用快速近似算法
explainer = shap.TreeExplainer(model, approximate=True)
shap_values = explainer.shap_values(X_test)

# 方案 3：对特定样本单独计算，按需解释
def explain_single_customer(customer_idx):
    customer_data = X_test.iloc[[customer_idx]]
    return explainer.shap_values(customer_data)

# 创建高度相关的特征
X_corr = X_test.copy()
X_corr['login_count_double'] = X_corr['login_count'] * 2
X_corr['login_count_triple'] = X_corr['login_count'] * 3

# 重新计算 SHAP 值
explainer_corr = shap.TreeExplainer(model)
shap_values_corr = explainer_corr.shap_values(X_corr)

# 看看特征重要性
plt.figure(figsize=(10, 8))
shap.summary_plot(shap_values_corr, X_corr, plot_type="bar")
plt.savefig('correlated_features_shap.png', bbox_inches='tight')
plt.close()

# 计算相关性矩阵
correlation_matrix = X_train.corr()

# 找出高相关性的特征对
high_corr_pairs = []
for i in range(len(correlation_matrix.columns)):
    for j in range(i+1, len(correlation_matrix.columns)):
        if abs(correlation_matrix.iloc[i, j]) > 0.8:
            high_corr_pairs.append(
                (correlation_matrix.columns[i], correlation_matrix.columns[j], correlation_matrix.iloc[i, j])
            )

print("高相关性特征对（阈值 0.8）：")
for feat1, feat2, corr in high_corr_pairs:
    print(f"{feat1} - {feat2}: {corr:.3f}")

# 根据相关性结果调整特征集

# 对同一个客户多次计算 SHAP 值
stability_test = []
for _ in range(10):
    shap_val = explainer.shap_values(customer_data)[0]
    stability_test.append(shap_val)

# 计算标准差
stability_df = pd.DataFrame(stability_test, columns=customer_data.columns.tolist())
stability_std = stability_df.std()

print("SHAP 值的标准差：")
print(stability_std.sort_values(ascending=False).head(10))