企业官网建设流程全解析

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 整体技术路线

整个项目我分成了四个核心模块：

下面我会逐个拆解每个模块的实现细节和背后的数学原理。

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_{i=1}^{n} \hat{y}_i^{y_i} (1 - \hat{y}_i)^{1-y_i}$$

取对数变成：

$$\ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} [y_i \log(\hat{y}_i) + (1-y_i) \log(1-\hat{y}_i)]$$

我们要最大化这个对数似然，等价于最小化负对数似然：

$$J(\boldsymbol{\beta}) = -\ell(\boldsymbol{\beta}) = -\sum_{i=1}^{n} [y_i \log(\hat{y}_i) + (1-y_i) \log(1-\hat{y}_i)]$$

这就是逻辑回归的损失函数，也叫交叉熵损失。在金融场景下，这个损失函数有一个很好的解释：它衡量的是模型预测概率与真实标签之间的“信息距离”。

2.4 梯度推导：手算一遍才记得住

对损失函数求关于 $\boldsymbol{\beta}$ 的偏导数，这是整个实现中最关键的一步。我用矩阵形式推导一遍。

设 $\mathbf{X}$ 是 $n \times (p+1)$ 的设计矩阵（第一列全为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}^{(t+1)} = \boldsymbol{\beta}^{(t)} - \alpha \nabla J(\boldsymbol{\beta}^{(t)})$$

其中 $\alpha$ 是学习率。在代码实现中，有几个细节需要注意：

import numpy as np def sigmoid(z): """数值稳定的Sigmoid实现""" # 对正负分别处理，避免exp溢出 result = np.zeros_like(z, dtype=np.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_rate=0.01, max_iter=1000, tol=1e-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}^{(t+1)} = \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_iter=100, tol=1e-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个（包括一些噪声特征），结果如下：

牛顿法在收敛速度和精度上都碾压梯度下降。但注意，当特征之间存在严重多重共线性时，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 一个完整的输出示例

假设我们有一个违约预测模型，输出结果如下表：

从这张表可以读出很多信息：年龄越大违约概率越低（系数为负），月收入越高违约概率越低，负债率越高违约概率越高，查询次数越多违约概率越高。所有变量的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, axis=1) 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, rcond=None) 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的计算逻辑

在金融风控中，连续变量通常不会直接放进逻辑回归，而是先做分箱，然后计算WOE（Weight 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)$$

IV（Information Value）是WOE的加权和：

$$\text{IV} = \sum_{i=1}^{k} (\text{好样本占比}_i - \text{坏样本占比}_i) \times \text{WOE}_i$$

IV值用来衡量变量的预测能力，经验规则是：IV<0.02预测能力极弱，0.02-0.1弱，0.1-0.3中等，0.3-0.5强，>0.5极强但需要警惕过拟合。

def woe_binning(X_col, y, n_bins=5, min_bin_pct=0.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[i+1]) 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[i+1])] = woe return woe_map, iv_total

5.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_score=600, base_odds=50, pdo=20): """ 将逻辑回归模型转换为评分卡 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的计算

AUC（Area 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 常见问题速查表

6.3 我踩过的几个坑

坑一：忘记加截距项。第一次手写的时候，我只把特征矩阵传进去，忘了加一列全1的截距列。结果模型被迫通过原点，拟合效果差了很多。后来养成了习惯，在数据预处理阶段就统一加上截距列。

坑二：标准化后忘记还原系数。为了加速收敛，我对特征做了标准化。但标准化后的系数是针对标准化变量的，要还原到原始尺度需要除以标准差。这个转换在计算评分卡时特别重要，否则分数会完全不对。

坑三：用准确率评估不平衡数据。金融违约数据通常正负样本比例悬殊，比如1:99。这时候一个把所有样本都预测为“不违约”的模型准确率也有99%，但完全没有意义。所以一定要用AUC和KS，不要用准确率。

坑四：忽视时间外样本验证。金融风控模型最怕的是过拟合历史数据。我一般会把数据按时间切分，用前80%的时间段训练，后20%的时间段验证。如果时间外AUC比训练集AUC下降超过0.05，说明模型可能过拟合了。

7. 完整代码整合与性能优化

7.1 把各个模块串起来

把前面的模块整合成一个完整的类：

class LogisticRegressionCustom: """ 手写逻辑回归，支持梯度下降和牛顿法 """ def __init__(self, method='newton', learning_rate=0.1, max_iter=1000, tol=1e-6, reg_lambda=0.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, threshold=0.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_history

7.2 性能优化建议

当数据量很大时，手写实现可能会比优化过的库慢。几个优化方向：

向量化操作：尽量避免Python循环，用NumPy的矩阵运算。我上面的代码已经做了向量化，但在计算WOE分箱时还有循环，可以用pandas的groupby加速。

使用float32：如果精度要求不是特别高，把数据类型从float64改成float32可以节省一半内存，速度也能提升20%左右。

随机梯度下降：当样本量超过百万时，批量梯度下降每次迭代都要遍历所有样本，太慢了。可以改用小批量梯度下降，每次只用一部分样本计算梯度。

利用稀疏矩阵：金融数据中很多特征（比如职业类别）是one-hot编码的，非常稀疏。用scipy的稀疏矩阵可以大幅减少内存和计算量。

7.3 与sklearn的结果对比

为了验证手写实现的正确性，我用同一份数据分别跑了手写实现和sklearn，结果如下：

系数差异在1e-4量级，说明手写实现是正确的。差异主要来自收敛精度的不同，sklearn默认的tol是1e-4，我设的是1e-8，所以手写实现收敛得更彻底。

8. 后续扩展方向

这个手写实现目前覆盖了逻辑回归的核心功能，但还有几个可以继续扩展的方向。

第一个是L1正则化。目前只实现了L2正则化，L1正则化可以做特征选择，在金融场景下很有用。但L1的梯度在0点不可导，需要用次梯度或者坐标下降法。

第二个是贝叶斯逻辑回归。金融场景下，先验知识很重要。比如你知道某个变量的系数应该是正的，可以用贝叶斯方法把这个先验加进去。实现上可以用拉普拉斯近似或者MCMC采样。

第三个是分布式版本。当数据量大到单机内存放不下时，可以用分块计算Hessian矩阵，或者用ADMM算法做分布式求解。

第四个是自动分箱。目前的分箱是等频分箱，比较粗糙。可以用决策树分箱或者卡方分箱，让WOE的单调性更好。

这些扩展我在后续项目中会逐步实现，但核心的逻辑回归手写框架已经够用了。关键是理解每一步的数学原理，这样不管用什么工具，你都知道它在背后做了什么。