QIAN数据:Lasso风控模型Python应用

import numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom sklearn.linear_model import LogisticRegression, LassoCVfrom sklearn.model_selection import train_test_splitfrom sklearn.metrics import roc_auc_score, roc_curvefrom sklearn.preprocessing import StandardScaler, PolynomialFeaturesimport warningswarnings.filterwarnings('ignore')plt.rcParams['font.sans-serif'] = ['Microsoft YaHei', 'DejaVu Sans']plt.rcParams['axes.unicode_minus'] = Falsenp.random.seed(42)

# ---------- 加载 German Credit ----------url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data'cols = [    'checking_status', 'duration', 'credit_history', 'purpose',    'credit_amount', 'savings', 'employment', 'installment_rate',    'personal_status', 'other_debtors', 'residence_since',    'property', 'age', 'other_installment', 'housing',    'existing_credits', 'job', 'dependents', 'telephone', 'foreign_worker',    'target']df = pd.read_csv(url, sep=' ', header=None, names=cols)df['target'] = df['target'].map({1: 0, 2: 1})print(f'原始数据: {df.shape[0]} 条, {df.shape[1]-1} 个原始变量')print(f'坏账率: {df.target.mean():.1%}')# 预期输出：原始数据: 1000 条, 20 个原始变量#         坏账率: 30.0%

# ---------- 特征工程 ----------# 1. 类别变量 → 独热编码cat_cols = ['checking_status', 'credit_history', 'purpose', 'savings',            'employment', 'personal_status', 'other_debtors', 'property',            'other_installment', 'housing', 'job', 'telephone', 'foreign_worker']df_encoded = pd.get_dummies(df[cat_cols], drop_first=False)n_cat = df_encoded.shape[1]print(f'独热编码后类别变量: {n_cat} 个')# 2. 数值变量 → 多项式（degree=2，含交互项）num_cols = ['duration', 'credit_amount', 'installment_rate',            'residence_since', 'age', 'existing_credits', 'dependents']poly = PolynomialFeatures(degree=2, include_bias=False)num_poly = poly.fit_transform(df[num_cols])n_num = num_poly.shape[1]print(f'多项式扩展后数值变量: {n_num} 个')# 3. 合并X_engineered = np.hstack([df_encoded.values, num_poly])feature_names = list(df_encoded.columns) + \    [f'poly_{i}' for i in range(n_num)]y = df['target'].valuesn_total = X_engineered.shape[1]print(f'最终特征总数: {n_total} 个（{n_cat} 类别独热 + {n_num} 多项式）')# 预期输出：独热编码后类别变量: 54 个#         多项式扩展后数值变量: 35 个#         最终特征总数: 89 个（54 类别独热 + 35 多项式）

# 标准化scaler = StandardScaler()X_scaled = scaler.fit_transform(X_engineered)

# 划分X_train, X_test, y_train, y_test = train_test_split(    X_scaled, y, test_size=0.3, random_state=42)# LassoCV自动选λlasso_cv = LassoCV(    alphas=np.logspace(-3, 1, 100),    cv=5, max_iter=10000, random_state=42)lasso_cv.fit(X_train, y_train)selected_idx = np.abs(lasso_cv.coef_) > 1e-4n_selected = selected_idx.sum()selected_features = [feature_names[i] for i in range(n_total)                     if np.abs(lasso_cv.coef_[i]) > 1e-4]print(f'最优 λ: {lasso_cv.alpha_:.4f}')print(f'Lasso 保留变量: {n_selected} / {n_total}')print(f'保留的工程变量: {n_selected} 个')for f in selected_features[:15]:    idx = feature_names.index(f)    print(f'  + {f}: {lasso_cv.coef_[idx]:.4f}')print(f'...（共 {n_selected} 个）')

alphas = np.logspace(-0.00005, -3, 50)coef_path = []for alpha in alphas:    lasso = LogisticRegression(        penalty='l1', C=1/alpha, solver='saga',        max_iter=10000, random_state=42    )    lasso.fit(X_train, y_train)    coef_path.append(lasso.coef_[0])coef_path = np.array(coef_path)# 只画非零系数的变量路径（前15条）nonzero_at_end = np.abs(coef_path[-1, :]) > 0.01top_indices = np.argsort(np.max(np.abs(coef_path), axis=0))[::-1][:15]plt.figure(figsize=(10, 6))for idx in top_indices:    plt.plot(alphas, coef_path[:, idx], linewidth=1.2,             label=feature_names[idx] if len(feature_names[idx]) < 20             else feature_names[idx][:18] + '...')plt.xscale('log')plt.xlabel('λ (log scale)', fontsize=12)plt.ylabel('回归系数', fontsize=12)plt.title('Lasso正则化路径（German Credit, 89→25个变量）', fontsize=14)plt.axvline(x=lasso_cv.alpha_, color='red', linestyle='--',            alpha=0.7, label=f'最优 λ = {lasso_cv.alpha_:.4f}')plt.legend(loc='best', fontsize=7)plt.tight_layout()plt.savefig('fig_lasso_path_german.png', dpi=150, bbox_inches='tight')plt.close()print('正则化路径图已保存')

# 全变量Logistic回归（无惩罚）lr_full = LogisticRegression(penalty=None, max_iter=10000)lr_full.fit(X_train, y_train)# Lasso稀疏模型（只保留被选中的变量）X_train_sparse = X_train[:, selected_idx]X_test_sparse = X_test[:, selected_idx]lr_sparse = LogisticRegression(penalty=None, max_iter=10000)lr_sparse.fit(X_train_sparse, y_train)# AUC对比auc_full = roc_auc_score(y_test, lr_full.decision_function(X_test))auc_sparse = roc_auc_score(y_test, lr_sparse.decision_function(X_test_sparse))print('=== 模型对比 ===')print(f'{"模型":<25}{"变量数":<10}{"AUC":<10}')print('-' * 45)print(f'{"全变量模型(89)":<25}{n_total:<10}{auc_full:.4f}')print(f'{"Lasso稀疏模型(25)":<25}{n_selected:<10}{auc_sparse:.4f}')print(f'{"AUC变化":<25}{"":<10}{auc_sparse - auc_full:+.4f}')

# ============ WOE编码函数 ============def calc_woe(series, target):    """对单个变量做WOE编码，返回(WOE值序列, 分组→WOE映射)"""    raw = series.to_numpy()    s = pd.Series(raw)    if not pd.api.types.is_numeric_dtype(s):        groups = s  # 类别变量：直接用原始分组    else:        try:            groups = pd.qcut(s, q=min(10, s.nunique()),                             duplicates='drop')        except:            groups = pd.cut(s, bins=min(10, s.nunique()))    tbl = pd.crosstab(groups, target, margins=False)    for c in [0, 1]:        if c not in tbl.columns: tbl[c] = 0    total_good = max(tbl[0].sum(), 1)    total_bad = max(tbl[1].sum(), 1)    p_good = tbl[0] / total_good    p_bad = tbl[1] / total_bad    p_good = p_good.replace(0, 0.0001)    p_bad = p_bad.replace(0, 0.0001)    woe_series = np.log(p_good / p_bad)    woe_map = dict(zip(tbl.index, woe_series))    return series.map(woe_map), woe_mapdef build_woe_model(train_df, test_df, var_list):    """对变量列表做WOE编码 → 无惩罚LR → 返回AUC和分数"""    Xtr_parts, Xte_parts = [], []    for var in var_list:        woe_tr, woe_map = calc_woe(train_df[var], train_df['target'])        woe_te = test_df[var].map(woe_map).fillna(0)        Xtr_parts.append(woe_tr.values.reshape(-1, 1))        Xte_parts.append(woe_te.values.reshape(-1, 1))    Xtr_woe = np.hstack(Xtr_parts)    Xte_woe = np.hstack(Xte_parts)    lr = LogisticRegression(penalty=None, max_iter=10000, random_state=42)    lr.fit(Xtr_woe, train_df['target'])    s_tr = lr.decision_function(Xtr_woe)    s_te = lr.decision_function(Xte_woe)    return (roc_auc_score(train_df['target'], s_tr),            roc_auc_score(test_df['target'], s_te),            s_tr, s_te)# 准备原始DataFrameoriginal_cols = cat_cols + num_cols  # 20个原始变量名# （cat_cols和num_cols在特征工程中已定义过，此处复用）# ============ 模型1：全变量(89) — 89个标准化特征 ============lr_all = LogisticRegression(penalty=None, max_iter=10000, random_state=42)lr_all.fit(X_train, y_train)s_all_tr = lr_all.decision_function(X_train)s_all_te = lr_all.decision_function(X_test)auc_all_tr = roc_auc_score(y_train, s_all_tr)auc_all_te = roc_auc_score(y_test, s_all_te)# ============ 模型2：全变量WOE(20) ============auc_woe20_tr, auc_woe20_te, s_woe20_tr, s_woe20_te = build_woe_model(    df_train, df_test, original_cols)# ============ 模型3：Lasso直出(25) ============X_train_sparse = X_train[:, selected_idx]X_test_sparse = X_test[:, selected_idx]lr_sparse = LogisticRegression(penalty=None, max_iter=10000, random_state=42)lr_sparse.fit(X_train_sparse, y_train)s_sp_tr = lr_sparse.decision_function(X_train_sparse)s_sp_te = lr_sparse.decision_function(X_test_sparse)auc_sp_tr = roc_auc_score(y_train, s_sp_tr)auc_sp_te = roc_auc_score(y_test, s_sp_te)# ============ 模型4：Lasso+WOE(14) ============# 将Lasso选中的展平特征映射回原始变量def map_back(var_name):    for cat in cat_cols:        if var_name == cat or var_name.startswith(cat + '_'):            return cat    for num in num_cols:        if num in var_name:            return num    return Nonelasso_orig_vars = set()for f in selected_features:    v = map_back(f)    if v: lasso_orig_vars.add(v)lasso_orig_vars = sorted(lasso_orig_vars)print(f'Lasso选中的展平特征 → {len(lasso_orig_vars)}个原始变量:')print(', '.join(lasso_orig_vars))auc_lw_tr, auc_lw_te, s_lw_tr, s_lw_te = build_woe_model(    df_train, df_test, lasso_orig_vars)# ============ KS + Lift ============def calc_ks(y_true, score):    order = np.argsort(score)    pos = np.cumsum(y_true[order])    neg = np.cumsum(1 - y_true[order])    return np.max(np.abs(pos / pos[-1] - neg / neg[-1]))def cum_capture_d3(y_true, score):    order = np.argsort(score)[::-1]    n = len(order)    return y_true[order][:n//10*3].sum() / max(y_true.sum(), 1) * 100# 打印对比表print(f'\n{"="*60}')print(f'{"四模型总对比":^60}')print(f'{"="*60}')print(f'{"模型":<20s}{"变量":<6s}{"AUC训练":<10s}{"AUC测试":<10s}'      f' {"KS训练":<10s}{"KS测试":<10s}')print('-'*60)for name, nv, s_tr, s_te in [    ('全变量(89)', 89, s_all_tr, s_all_te),    ('全变量WOE(20)', 20, s_woe20_tr, s_woe20_te),    (f'Lasso直出({n_selected})', n_selected, s_sp_tr, s_sp_te),    (f'Lasso+WOE({len(lasso_orig_vars)})', len(lasso_orig_vars), s_lw_tr, s_lw_te),]:    auc_tr = roc_auc_score(y_train, s_tr)    auc_te = roc_auc_score(y_test, s_te)    ks_tr = calc_ks(y_train, s_tr)    ks_te = calc_ks(y_test, s_te)    print(f'{name:<20s}{nv:<6d}{auc_tr:<10.4f}{auc_te:<10.4f} '          f'{ks_tr:<10.4f}{ks_te:<10.4f}')# ============ 10等分Lift对比 ============def calc_decile_lift(y_true, score, n_deciles=10):    """计算10等分Lift表，返回每档的(坏账数, 坏账率, Lift系数)"""    order = np.argsort(score)[::-1]  # 风险从高到低    n = len(order)    decile_size = n // n_deciles    results = []    base_rate = y_true.mean()    for d in range(n_deciles):        start = d * decile_size        end = start + decile_size if d < n_deciles - 1 else n        idx = order[start:end]        n_bad = y_true[idx].sum()        rate = n_bad / len(idx)        lift = rate / base_rate if base_rate > 0 else 0        results.append((int(n_bad), rate, lift))    return resultsmodel_scores = [    ('全变量(89)', 89, s_all_te),    ('全变量WOE(20)', 20, s_woe20_te),    (f'Lasso直出({n_selected})', n_selected, s_sp_te),    (f'Lasso+WOE({len(lasso_orig_vars)})', len(lasso_orig_vars), s_lw_te),]print()print('完整的10等分Lift对比（测试集）：')print()header = f'{"分档":<10}'for name, _, _ in model_scores:    header += f'{name:<30}'print(header)for d in range(10):    row = f'D{d+1:<9}'    for _, _, score in model_scores:        results = calc_decile_lift(y_test, score)        n_bad, rate, lift = results[d]        row += f'{n_bad:<4}{rate*100:>5.1f}%  {lift:>5.2f}x{"":8}'    print(row)# D1-D3累计捕获率print()print('D1-D3累计捕获率（最高风险三档捕捉到的坏账占比）：')print()print(f'{"模型":<25}{"训练":<12}{"测试":<12}')print('-' * 49)for name, _, s_tr in [    ('全变量(89)', 89, s_all_tr),    ('全变量WOE(20)', 20, s_woe20_tr),    (f'Lasso直出({n_selected})', n_selected, s_sp_tr),    (f'Lasso+WOE({len(lasso_orig_vars)})', len(lasso_orig_vars), s_lw_tr),]:    # Find matching test score    for nm, _, s_te in model_scores:        if nm == name:            c3_tr = cum_capture_d3(y_train, s_tr)            c3_te = cum_capture_d3(y_test, s_te)            print(f'{name:<25}{c3_tr:>5.1f}%{"":>8}{c3_te:>5.1f}%')            break# 单调性检查print()print('单调性检查（各档坏账率是否严格递减）：')print()checks = ['全变量(89)', '全变量WOE(20)', f'Lasso直出({n_selected})', f'Lasso+WOE({len(lasso_orig_vars)})']for name, _, score in model_scores:    results = calc_decile_lift(y_test, score)    violations = 0    out = f'{name:<20s}'    for d in range(1, 10):        if results[d][1] <= results[d-1][1]:            out += f' D{d}→D{d+1}✓'        else:            out += f' D{d}→D{d+1}×({results[d-1][2]:.2f}→{results[d][2]:.2f})'            violations += 1    out += f'  违反{ violations}次'    print(out)