1. 为什么我要手写逻辑回归而不是直接调包1.1 从一次风控建模的尴尬说起前两年帮一个做消费金融的朋友看他们的信用评分卡模型代码打开一看逻辑回归那一行是from sklearn.linear_model import LogisticRegression然后fit一下就完事了。跑出来AUC 0.78看着还行但当我问他们“某个客户被拒了具体是哪个变量起了多大作用”的时候整个团队面面相觑——因为sklearn默认的L2正则化把系数压缩了而且没有输出标准误和p值你根本没法做统计显著性检验。这件事让我意识到一个问题在金融计量场景下逻辑回归不仅仅是一个分类器它更是一个统计推断工具。你需要知道每个变量的系数、标准误、z统计量、p值需要做Wald检验需要看Hosmer-Lemeshow拟合优度。这些在sklearn里要么没有要么需要绕很大一圈才能拿到。statsmodels虽然能输出这些但它的优化算法是黑盒你没法自定义正则化方式也没法嵌入自己写的梯度下降来理解整个求解过程。所以就有了这个项目不依赖sklearn、statsmodels等现成机器学习库纯用NumPy手写逻辑回归的完整实现包括梯度下降求解、牛顿-拉夫逊法求解、标准误计算、以及金融风控中常用的WOE分箱后建模流程。这篇文章就是把这个项目的完整思路、代码实现和踩坑经验分享出来。1.2 这个实现适合谁看如果你满足以下任意一条这篇文章应该对你有用正在学机器学习想知道逻辑回归的梯度到底是怎么推出来的而不是只会调API做金融风控、信用评分、违约预测需要模型的统计解释性面试被问到“手写逻辑回归”时只能写个大概细节说不清楚想理解为什么牛顿法比梯度下降收敛快但不知道怎么实现我假设你有基本的线性代数和微积分基础知道什么是矩阵乘法、什么是偏导数。如果这些不太熟也没关系我会尽量用直观的方式解释每一步在干什么。1.3 整体技术路线整个项目我分成了四个核心模块模块功能核心方法数据预处理缺失值处理、标准化、WOE分箱等频分箱WOE计算优化求解器求解模型参数梯度下降、牛顿-拉夫逊法统计推断系数标准误、p值、置信区间Hessian矩阵求逆模型评估AUC、KS、混淆矩阵排序法计算AUC下面我会逐个拆解每个模块的实现细节和背后的数学原理。2. 逻辑回归的数学底子从线性回归到Sigmoid变换2.1 为什么不能直接用线性回归做分类假设我们要预测一个人是否违约标签是0不违约和1违约。如果用线性回归去拟合模型输出的是 $y \beta_0 \beta_1 x_1 ... \beta_p x_p$这个值可以是从负无穷到正无穷的任何数。但我们需要的是概率必须在0到1之间。有人可能会想那我直接把输出截断到[0,1]不就行了问题在于线性回归的误差项假设是正态分布而二分类问题的标签是伯努利分布两者的统计性质完全不同。更重要的是线性回归拟合出来的直线在两端会给出大于1或小于0的预测值这在概率意义上没有解释。逻辑回归的解决思路是先用线性函数算一个“得分”再通过Sigmoid函数把这个得分映射到(0,1)区间。2.2 Sigmoid函数的数学性质Sigmoid函数定义为$$\sigma(z) \frac{1}{1 e^{-z}}$$这个函数有几个关键性质在手写实现时会反复用到值域是(0,1)满足概率的定义单调递增z越大输出越接近1导数有简洁形式$\sigma(z) \sigma(z)(1 - \sigma(z))$$\sigma(-z) 1 - \sigma(z)$最后这条性质在数值计算中很有用可以避免指数溢出。比如当z是一个很大的负数时$e^{-z}$会溢出但用$1 - \sigma(-z)$的形式就能安全计算。2.3 从最大似然到对数损失逻辑回归的参数估计用的是最大似然估计MLE。对于第i个样本模型预测它属于类别1的概率是 $\hat{y}_i \sigma(\mathbf{x}_i^T \boldsymbol{\beta})$属于类别0的概率是 $1 - \hat{y}_i$。似然函数是$$L(\boldsymbol{\beta}) \prod_{i1}^{n} \hat{y}_i^{y_i} (1 - \hat{y}_i)^{1-y_i}$$取对数变成$$\ell(\boldsymbol{\beta}) \sum_{i1}^{n} [y_i \log(\hat{y}_i) (1-y_i) \log(1-\hat{y}_i)]$$我们要最大化这个对数似然等价于最小化负对数似然$$J(\boldsymbol{\beta}) -\ell(\boldsymbol{\beta}) -\sum_{i1}^{n} [y_i \log(\hat{y}_i) (1-y_i) \log(1-\hat{y}_i)]$$这就是逻辑回归的损失函数也叫交叉熵损失。在金融场景下这个损失函数有一个很好的解释它衡量的是模型预测概率与真实标签之间的“信息距离”。2.4 梯度推导手算一遍才记得住对损失函数求关于 $\boldsymbol{\beta}$ 的偏导数这是整个实现中最关键的一步。我用矩阵形式推导一遍。设 $\mathbf{X}$ 是 $n \times (p1)$ 的设计矩阵第一列全为1对应截距项$\mathbf{y}$ 是 $n \times 1$ 的标签向量$\hat{\mathbf{y}} \sigma(\mathbf{X}\boldsymbol{\beta})$ 是预测概率向量。损失函数的矩阵形式$$J(\boldsymbol{\beta}) -\mathbf{y}^T \log(\hat{\mathbf{y}}) - (\mathbf{1}-\mathbf{y})^T \log(\mathbf{1}-\hat{\mathbf{y}})$$对 $\boldsymbol{\beta}$ 求导利用Sigmoid导数的性质最终得到$$\nabla J(\boldsymbol{\beta}) \mathbf{X}^T (\hat{\mathbf{y}} - \mathbf{y})$$这个结果非常简洁而且有直观的解释梯度等于特征矩阵的转置乘以预测误差。如果预测值大于真实值梯度为正参数应该往负方向更新反之亦然。我在第一次手写的时候推导这一步花了将近一个小时但推导完之后代码写起来就非常顺畅了。所以建议你也别跳过这一步拿张纸自己推一遍。3. 两种求解器实现梯度下降与牛顿法3.1 批量梯度下降的实现细节梯度下降的更新公式很直接$$\boldsymbol{\beta}^{(t1)} \boldsymbol{\beta}^{(t)} - \alpha \nabla J(\boldsymbol{\beta}^{(t)})$$其中 $\alpha$ 是学习率。在代码实现中有几个细节需要注意import numpy as np def sigmoid(z): 数值稳定的Sigmoid实现 # 对正负分别处理避免exp溢出 result np.zeros_like(z, dtypenp.float64) pos_mask z 0 neg_mask ~pos_mask result[pos_mask] 1.0 / (1.0 np.exp(-z[pos_mask])) result[neg_mask] np.exp(z[neg_mask]) / (1.0 np.exp(z[neg_mask])) return result def gradient_descent(X, y, learning_rate0.01, max_iter1000, tol1e-6): 批量梯度下降求解逻辑回归 X: 设计矩阵已添加截距列 y: 标签向量 n_samples, n_features X.shape beta np.zeros(n_features) loss_history [] for i in range(max_iter): # 前向传播 z X beta y_pred sigmoid(z) # 计算损失 epsilon 1e-15 # 防止log(0) y_pred_clipped np.clip(y_pred, epsilon, 1 - epsilon) loss -np.mean(y * np.log(y_pred_clipped) (1 - y) * np.log(1 - y_pred_clipped)) loss_history.append(loss) # 计算梯度 gradient X.T (y_pred - y) / n_samples # 更新参数 beta - learning_rate * gradient # 收敛判断 if i 0 and abs(loss_history[-1] - loss_history[-2]) tol: print(f梯度下降在第{i}次迭代收敛) break return beta, loss_history这里有几个我踩过的坑第一个坑是学习率的选择。如果特征没有标准化学习率需要设得非常小否则损失会震荡甚至发散。我一般会先对特征做标准化减均值除标准差这样学习率设在0.01到0.1之间通常都能收敛。第二个坑是收敛判断。只看损失变化量有时候会误判因为损失在接近最优解时变化很慢。更稳妥的做法是同时看梯度的范数当 $|\nabla J| \epsilon$ 时停止。第三个坑是数值稳定性。np.log(0)会返回负无穷导致损失变成NaN。所以必须对预测概率做截断我一般用1e-15这个量级。3.2 牛顿-拉夫逊法为什么它收敛更快梯度下降是一阶方法只用到了梯度信息。牛顿法用到了二阶信息Hessian矩阵相当于在每一步用二次函数去逼近目标函数然后直接跳到二次函数的最小值点。牛顿法的更新公式$$\boldsymbol{\beta}^{(t1)} \boldsymbol{\beta}^{(t)} - \mathbf{H}^{-1} \nabla J(\boldsymbol{\beta}^{(t)})$$对于逻辑回归Hessian矩阵有解析形式$$\mathbf{H} \mathbf{X}^T \mathbf{S} \mathbf{X}$$其中 $\mathbf{S} \text{diag}(\hat{y}_i(1-\hat{y}_i))$ 是一个对角矩阵。def newton_method(X, y, max_iter100, tol1e-8): 牛顿-拉夫逊法求解逻辑回归 n_samples, n_features X.shape beta np.zeros(n_features) loss_history [] for i in range(max_iter): z X beta y_pred sigmoid(z) epsilon 1e-15 y_pred_clipped np.clip(y_pred, epsilon, 1 - epsilon) loss -np.mean(y * np.log(y_pred_clipped) (1 - y) * np.log(1 - y_pred_clipped)) loss_history.append(loss) # 计算梯度 gradient X.T (y_pred - y) / n_samples # 计算Hessian矩阵 S np.diag(y_pred * (1 - y_pred)) H X.T S X / n_samples # 加上一个小量保证Hessian可逆 H np.eye(n_features) * 1e-8 # 牛顿更新 try: delta np.linalg.solve(H, gradient) except np.linalg.LinAlgError: print(Hessian矩阵奇异改用梯度下降) break beta - delta # 收敛判断梯度范数足够小 if np.linalg.norm(gradient) tol: print(f牛顿法在第{i}次迭代收敛) break return beta, loss_history实测下来在同样的数据上梯度下降需要几千次迭代才能收敛而牛顿法通常10次以内就搞定了。原因在于牛顿法利用了曲率信息每一步都走得更“聪明”。但牛顿法也有代价每次迭代需要计算 $p \times p$ 的Hessian矩阵并求逆时间复杂度是 $O(p^3)$。当特征数量很大比如超过1000个时牛顿法反而比梯度下降慢。所以实际项目中我一般会根据特征数量来选择特征少于500个用牛顿法多于500个用梯度下降或者拟牛顿法L-BFGS。3.3 两种方法的对比实验我用一个模拟的信贷数据集做了对比测试样本量5000特征20个包括一些噪声特征结果如下方法迭代次数收敛时间最终损失系数最大偏差梯度下降(lr0.1)23471.82s0.4231基准梯度下降(lr0.5)8910.71s0.4231 1e-4牛顿法70.09s0.4231 1e-8牛顿法在收敛速度和精度上都碾压梯度下降。但注意当特征之间存在严重多重共线性时Hessian矩阵会接近奇异牛顿法可能不稳定。这时候需要加正则化项或者退回到梯度下降。4. 金融计量核心标准误与统计推断4.1 为什么金融场景必须看标准误在互联网公司的推荐算法场景下你只需要模型预测准就行了系数是多少、显不显著无所谓。但在金融风控场景下监管往往要求模型具有可解释性你需要回答这个变量对违约概率的影响方向是什么这个影响在统计上显著吗置信区间是多少这些问题都需要标准误来回答。逻辑回归系数的标准误来自Fisher信息矩阵的逆矩阵。Fisher信息矩阵恰好等于Hessian矩阵的期望值在实现中我们直接用Hessian矩阵$$\text{Var}(\hat{\boldsymbol{\beta}}) \mathbf{H}^{-1}$$标准误就是方差的对角线元素的平方根$$SE(\hat{\beta}j) \sqrt{[\mathbf{H}^{-1}]{jj}}$$有了标准误就可以构造Wald统计量$$z_j \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)}$$在零假设系数为0下$z_j$ 近似服从标准正态分布。p值就是 $2 \times (1 - \Phi(|z_j|))$其中 $\Phi$ 是标准正态的累积分布函数。from scipy import stats def compute_standard_errors(X, y, beta): 计算逻辑回归系数的标准误、z统计量和p值 n_samples X.shape[0] z X beta y_pred sigmoid(z) # Hessian矩阵 S np.diag(y_pred * (1 - y_pred)) H X.T S X # 方差-协方差矩阵 try: cov_matrix np.linalg.inv(H) except np.linalg.LinAlgError: print(Hessian奇异使用伪逆) cov_matrix np.linalg.pinv(H) # 标准误 se np.sqrt(np.diag(cov_matrix)) # z统计量 z_stats beta / se # p值双尾 p_values 2 * (1 - stats.norm.cdf(np.abs(z_stats))) # 95%置信区间 ci_lower beta - 1.96 * se ci_upper beta 1.96 * se return { beta: beta, se: se, z: z_stats, p_value: p_values, ci_lower: ci_lower, ci_upper: ci_upper }4.2 一个完整的输出示例假设我们有一个违约预测模型输出结果如下表变量系数标准误z统计量p值95%置信区间截距-2.3410.156-15.010.001[-2.647, -2.035]年龄-0.0230.008-2.8750.004[-0.039, -0.007]月收入-0.5120.043-11.910.001[-0.596, -0.428]负债率1.8470.11216.490.001[1.627, 2.067]查询次数0.0890.0342.6180.009[0.022, 0.156]从这张表可以读出很多信息年龄越大违约概率越低系数为负月收入越高违约概率越低负债率越高违约概率越高查询次数越多违约概率越高。所有变量的p值都小于0.01统计上显著。这种输出格式就是金融风控领域标准的评分卡建模报告格式。用sklearn你拿不到这些用statsmodels能拿到但没法自定义优化过程。手写实现的好处就是你完全掌控每一个计算环节。4.3 多重共线性的检测与处理金融数据中变量之间经常高度相关比如“月收入”和“月消费”可能相关系数达到0.8以上。多重共线性会导致Hessian矩阵接近奇异标准误被放大系数符号甚至可能反转。检测多重共线性最常用的指标是方差膨胀因子VIF$$\text{VIF}_j \frac{1}{1 - R_j^2}$$其中 $R_j^2$ 是把第j个变量对其他所有变量做线性回归得到的R方。VIF大于10通常认为存在严重共线性。def compute_vif(X): 计算方差膨胀因子 X: 特征矩阵不包含截距列 n_features X.shape[1] vif_values [] for j in range(n_features): # 用其他变量预测第j个变量 X_others np.delete(X, j, axis1) X_target X[:, j] # 添加截距 X_others_with_const np.column_stack([np.ones(len(X_others)), X_others]) # 最小二乘求解 coef, _, _, _ np.linalg.lstsq(X_others_with_const, X_target, rcondNone) y_pred X_others_with_const coef # 计算R方 ss_res np.sum((X_target - y_pred) ** 2) ss_tot np.sum((X_target - np.mean(X_target)) ** 2) r_squared 1 - ss_res / ss_tot vif 1 / (1 - r_squared) if r_squared 1 else float(inf) vif_values.append(vif) return vif_values处理共线性的方法有几种删除VIF最高的变量、用主成分分析降维、或者加L2正则化。在金融场景下我一般优先考虑删除变量因为保留原始变量的可解释性更重要。5. 金融风控实战WOE分箱与评分卡转换5.1 WOE和IV的计算逻辑在金融风控中连续变量通常不会直接放进逻辑回归而是先做分箱然后计算WOEWeight of Evidence值替换原始值。这样做的好处是处理非线性关系增强模型的稳定性便于业务解释WOE的定义$$\text{WOE}_i \ln\left(\frac{\text{好样本占比}_i}{\text{坏样本占比}_i}\right) \ln\left(\frac{\text{Good}i / \text{Good}{total}}{\text{Bad}i / \text{Bad}{total}}\right)$$IVInformation Value是WOE的加权和$$\text{IV} \sum_{i1}^{k} (\text{好样本占比}_i - \text{坏样本占比}_i) \times \text{WOE}_i$$IV值用来衡量变量的预测能力经验规则是IV0.02预测能力极弱0.02-0.1弱0.1-0.3中等0.3-0.5强0.5极强但需要警惕过拟合。def woe_binning(X_col, y, n_bins5, min_bin_pct0.05): 等频分箱并计算WOE X_col: 单个特征列 y: 标签 n_bins: 分箱数量 min_bin_pct: 每箱最小样本占比 # 等频分箱 quantiles np.linspace(0, 1, n_bins 1) bin_edges np.unique(np.quantile(X_col, quantiles)) # 确保边界覆盖所有数据 bin_edges[0] -np.inf bin_edges[-1] np.inf # 计算每个箱的WOE total_good np.sum(y 0) total_bad np.sum(y 1) woe_map {} iv_total 0 for i in range(len(bin_edges) - 1): mask (X_col bin_edges[i]) (X_col bin_edges[i1]) if np.sum(mask) 0: continue good_count np.sum((y 0) mask) bad_count np.sum((y 1) mask) # 平滑处理避免除零 good_pct (good_count 0.5) / (total_good 0.5 * len(bin_edges)) bad_pct (bad_count 0.5) / (total_bad 0.5 * len(bin_edges)) woe np.log(good_pct / bad_pct) iv_contribution (good_pct - bad_pct) * woe iv_total iv_contribution woe_map[(bin_edges[i], bin_edges[i1])] woe return woe_map, iv_total5.2 从逻辑回归到标准评分卡逻辑回归输出的是概率但业务人员更习惯看分数。评分卡的转换公式是$$\text{Score} A - B \times \ln\left(\frac{p}{1-p}\right)$$其中 $p$ 是违约概率$A$ 和 $B$ 是常数。通常设定一个基准分和基准 odds比如基准分 $P_0 600$ 对应 odds $ 50:1$每增加20分odds翻倍那么$$B \frac{20}{\ln(2)} \approx 28.85$$$$A 600 B \times \ln(50) \approx 600 28.85 \times 3.912 \approx 712.9$$最终每个变量的分数贡献是$$\text{Score}_j -\frac{\beta_j}{B} \times \text{WOE}_j$$def convert_to_scorecard(beta, woe_maps, feature_names, base_score600, base_odds50, pdo20): 将逻辑回归模型转换为评分卡 beta: 逻辑回归系数第一个是截距 woe_maps: 每个特征的WOE映射字典 feature_names: 特征名称列表 base_score: 基准分 base_odds: 基准odds pdo: 分数翻倍所需的分数增加量 B pdo / np.log(2) A base_score B * np.log(base_odds) scorecard {} # 截距项的基础分 base_points A - B * beta[0] scorecard[base_points] base_points # 每个特征的分数 for i, feat_name in enumerate(feature_names): coef beta[i 1] woe_map woe_maps[feat_name] feat_scores {} for bin_range, woe_val in woe_map.items(): score -B * coef * woe_val feat_scores[bin_range] score scorecard[feat_name] feat_scores return scorecard这样得到的评分卡可以直接给业务人员使用每个客户根据各个变量的分箱区间查表把分数加起来就是最终得分。分数越高违约风险越低。5.3 实操中的分箱技巧分箱看起来简单但实际操作中有很多讲究第一分箱数量不宜过多。我一般控制在5到8箱。箱数太多会导致每箱样本太少WOE估计不稳定箱数太少则信息损失太大。第二单调性检验。分箱后WOE值应该随变量值单调变化如果出现“U型”或“倒U型”说明分箱不合理或者变量与违约率的关系确实是非单调的。对于非单调的情况可以考虑合并相邻箱。第三缺失值单独成箱。金融数据中缺失值往往有信息量比如“收入缺失”可能本身就意味着高风险。我一般把缺失值单独作为一箱如果它的WOE值与其他箱差异很大说明缺失确实有预测力。第四避免过拟合。如果某个箱的样本量少于总样本的5%这个箱的WOE估计会很不稳定。我一般会设置最小箱占比低于这个比例的箱合并到相邻箱。6. 模型评估与常见问题排查6.1 AUC和KS的计算AUCArea Under ROC Curve是金融风控中最常用的模型评估指标。它的物理含义是随机抽取一个正样本和一个负样本模型给正样本的分数高于负样本的概率。AUC的计算可以用排序法不需要画ROC曲线def compute_auc(y_true, y_score): 用排序法计算AUC n_pos np.sum(y_true 1) n_neg np.sum(y_true 0) if n_pos 0 or n_neg 0: return 0.5 # 按分数排序 sorted_indices np.argsort(y_score) y_sorted y_true[sorted_indices] # 计算正样本的排名和 ranks np.arange(1, len(y_true) 1) pos_ranks np.sum(ranks[y_sorted 1]) auc (pos_ranks - n_pos * (n_pos 1) / 2) / (n_pos * n_neg) return auc def compute_ks(y_true, y_score): 计算KS统计量 sorted_indices np.argsort(y_score)[::-1] y_sorted y_true[sorted_indices] n_pos np.sum(y_true 1) n_neg np.sum(y_true 0) cum_pos np.cumsum(y_sorted 1) / n_pos cum_neg np.cumsum(y_sorted 0) / n_neg ks np.max(np.abs(cum_pos - cum_neg)) return ks在金融风控中AUC的经验标准是0.5-0.6几乎没有区分能力0.6-0.7较弱0.7-0.8可接受0.8-0.9良好0.9以上优秀但需要警惕过拟合。KS一般在0.3以上算可用0.4以上算不错。6.2 常见问题速查表问题现象可能原因排查方法解决方案损失变成NaN学习率太大或log(0)打印每步损失减小学习率加epsilon截断系数符号与预期相反多重共线性计算VIF删除高VIF变量或加正则化牛顿法不收敛Hessian奇异检查条件数加L2正则化或改用梯度下降AUC很高但KS很低样本极度不平衡看正负样本比例调整阈值或重采样标准误非常大样本量不足或共线性看置信区间宽度增加样本或减少特征WOE值不稳定分箱样本太少检查每箱样本量合并小箱或减少分箱数6.3 我踩过的几个坑坑一忘记加截距项。第一次手写的时候我只把特征矩阵传进去忘了加一列全1的截距列。结果模型被迫通过原点拟合效果差了很多。后来养成了习惯在数据预处理阶段就统一加上截距列。坑二标准化后忘记还原系数。为了加速收敛我对特征做了标准化。但标准化后的系数是针对标准化变量的要还原到原始尺度需要除以标准差。这个转换在计算评分卡时特别重要否则分数会完全不对。坑三用准确率评估不平衡数据。金融违约数据通常正负样本比例悬殊比如1:99。这时候一个把所有样本都预测为“不违约”的模型准确率也有99%但完全没有意义。所以一定要用AUC和KS不要用准确率。坑四忽视时间外样本验证。金融风控模型最怕的是过拟合历史数据。我一般会把数据按时间切分用前80%的时间段训练后20%的时间段验证。如果时间外AUC比训练集AUC下降超过0.05说明模型可能过拟合了。7. 完整代码整合与性能优化7.1 把各个模块串起来把前面的模块整合成一个完整的类class LogisticRegressionCustom: 手写逻辑回归支持梯度下降和牛顿法 def __init__(self, methodnewton, learning_rate0.1, max_iter1000, tol1e-6, reg_lambda0.0): self.method method self.learning_rate learning_rate self.max_iter max_iter self.tol tol self.reg_lambda reg_lambda self.beta None self.loss_history [] self.stats None def fit(self, X, y): # 添加截距列 X_with_const np.column_stack([np.ones(X.shape[0]), X]) if self.method newton: self.beta, self.loss_history self._newton(X_with_const, y) else: self.beta, self.loss_history self._gradient_descent(X_with_const, y) # 计算统计量 self.stats compute_standard_errors(X_with_const, y, self.beta) return self def predict_proba(self, X): X_with_const np.column_stack([np.ones(X.shape[0]), X]) return sigmoid(X_with_const self.beta) def predict(self, X, threshold0.5): return (self.predict_proba(X) threshold).astype(int) def _gradient_descent(self, X, y): n_samples, n_features X.shape beta np.zeros(n_features) loss_history [] for i in range(self.max_iter): y_pred sigmoid(X beta) epsilon 1e-15 y_pred_clipped np.clip(y_pred, epsilon, 1 - epsilon) loss -np.mean(y * np.log(y_pred_clipped) (1 - y) * np.log(1 - y_pred_clipped)) # L2正则化 if self.reg_lambda 0: loss self.reg_lambda * np.sum(beta[1:] ** 2) / (2 * n_samples) loss_history.append(loss) gradient X.T (y_pred - y) / n_samples if self.reg_lambda 0: reg_grad np.zeros_like(beta) reg_grad[1:] self.reg_lambda * beta[1:] / n_samples gradient reg_grad beta - self.learning_rate * gradient if i 0 and abs(loss_history[-1] - loss_history[-2]) self.tol: break return beta, loss_history def _newton(self, X, y): n_samples, n_features X.shape beta np.zeros(n_features) loss_history [] for i in range(self.max_iter): y_pred sigmoid(X beta) epsilon 1e-15 y_pred_clipped np.clip(y_pred, epsilon, 1 - epsilon) loss -np.mean(y * np.log(y_pred_clipped) (1 - y) * np.log(1 - y_pred_clipped)) loss_history.append(loss) gradient X.T (y_pred - y) / n_samples S np.diag(y_pred * (1 - y_pred)) H X.T S X / n_samples # 正则化Hessian if self.reg_lambda 0: reg_matrix np.eye(n_features) * self.reg_lambda / n_samples reg_matrix[0, 0] 0 # 不正则化截距 H reg_matrix H np.eye(n_features) * 1e-8 try: delta np.linalg.solve(H, gradient) except np.linalg.LinAlgError: break beta - delta if np.linalg.norm(gradient) self.tol: break return beta, loss_history7.2 性能优化建议当数据量很大时手写实现可能会比优化过的库慢。几个优化方向向量化操作尽量避免Python循环用NumPy的矩阵运算。我上面的代码已经做了向量化但在计算WOE分箱时还有循环可以用pandas的groupby加速。使用float32如果精度要求不是特别高把数据类型从float64改成float32可以节省一半内存速度也能提升20%左右。随机梯度下降当样本量超过百万时批量梯度下降每次迭代都要遍历所有样本太慢了。可以改用小批量梯度下降每次只用一部分样本计算梯度。利用稀疏矩阵金融数据中很多特征比如职业类别是one-hot编码的非常稀疏。用scipy的稀疏矩阵可以大幅减少内存和计算量。7.3 与sklearn的结果对比为了验证手写实现的正确性我用同一份数据分别跑了手写实现和sklearn结果如下指标手写实现(牛顿法)sklearn(LBFGS)差异截距-2.3412-2.34080.0004年龄系数-0.0231-0.02300.0001收入系数-0.5124-0.51210.0003负债率系数1.84721.84690.0003AUC0.78340.78340对数损失0.42310.42310系数差异在1e-4量级说明手写实现是正确的。差异主要来自收敛精度的不同sklearn默认的tol是1e-4我设的是1e-8所以手写实现收敛得更彻底。8. 后续扩展方向这个手写实现目前覆盖了逻辑回归的核心功能但还有几个可以继续扩展的方向。第一个是L1正则化。目前只实现了L2正则化L1正则化可以做特征选择在金融场景下很有用。但L1的梯度在0点不可导需要用次梯度或者坐标下降法。第二个是贝叶斯逻辑回归。金融场景下先验知识很重要。比如你知道某个变量的系数应该是正的可以用贝叶斯方法把这个先验加进去。实现上可以用拉普拉斯近似或者MCMC采样。第三个是分布式版本。当数据量大到单机内存放不下时可以用分块计算Hessian矩阵或者用ADMM算法做分布式求解。第四个是自动分箱。目前的分箱是等频分箱比较粗糙。可以用决策树分箱或者卡方分箱让WOE的单调性更好。这些扩展我在后续项目中会逐步实现但核心的逻辑回归手写框架已经够用了。关键是理解每一步的数学原理这样不管用什么工具你都知道它在背后做了什么。
