ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

Matlab手写逻辑回归:从数学原理到多变量概率预测模型实现

Matlab手写逻辑回归:从数学原理到多变量概率预测模型实现 很多朋友第一次看到逻辑回归这四个字第一反应就是——这玩意儿是个回归模型吧我当年也是在Matlab里跑完一段代码看着输出的0.73、0.86这种概率值才回过神来这家伙其实是披着回归外衣的分类神器顺带还能兼职干概率预测的活。今天咱们就把这个事儿聊透用Matlab从零搭一个基于逻辑回归的多变量预测模型不调工具箱手写核心算法把每一步的数学原理和代码对应起来。这篇文章适合谁刚入门机器学习、被各种术语绕晕的Matlab用户以及那些想搞明白模型到底在算什么而非单纯调包的同学。1. 先纠正一个误解逻辑回归怕不是个回归模型很多教材上来就给公式不说人话导致一大半人把逻辑回归和线性回归混在一起。这里我用自己的理解给你捋一遍。1.1 名字里的回归到底从哪来逻辑回归的英文是Logistic Regression名字确实是回归但它干的事是分类。之所以叫回归是因为它最早是从线性回归的底子上长出来的——先算一个连续的线性组合再通过一个函数把它压成0到1之间的概率值最后根据概率做分类决策。换句话说回归是它的计算手段分类是它的最终目的。这个先算连续值再映射成概率的思路非常关键也正是它能用来做概率预测的原因。比如你想预测一个用户会不会流失逻辑回归输出的不是简单的会/不会而是一个流失概率0.87。这个0.87放在业务里就可以直接当风险分数用而不是只拿它跟0.5比大小。1.2 sigmoid函数一切的核心逻辑回归的灵魂就是sigmoid函数公式长这样[ h_\theta(x) \frac{1}{1 e^{-\theta^T x}} ]其中 ( \theta^T x ) 就是多个特征的线性加权求和。sigmoid的作用是把任意实数值压缩到(0,1)区间。你想想如果特征加权求和算出来是10或者是-20直接拿这个数去做分类肯定不行因为它的范围没有边界但经过sigmoid一压10对应接近1的概率-20对应接近0的概率这才符合概率的定义。用生活化的方式理解线性回归像是给一个螺丝拧螺母拧多少圈都有对应位置逻辑回归则是在这个基础上加了一个限位器不管你怎么用力拧最后的结果都被限制在一个固定区间里。这个限位器就是sigmoid。1.3 多变量逻辑回归的数学表达所谓多变量就是输入特征不止一个。举例来说假设你想预测某个用户会不会点击广告特征可能有历史点击次数、页面停留时间、设备类型编码、时段特征等等。假设总共有n个特征那么[ z \theta_0 \theta_1 x_1 \theta_2 x_2 \cdots \theta_n x_n ]在Matlab里这个求和用矩阵乘法一句话就能写完z X * theta;这里X是样本矩阵每一行是一个样本每一列是一个特征theta是参数向量。X乘theta的结果就是一个列向量每个元素是该样本的线性加权值。然后再整体套sigmoidh 1 ./ (1 exp(-z));这行代码就是逻辑回归的前向计算简单到让人怀疑。1.4 它到底适合解决什么问题逻辑回归适合解决以下几类问题二分类问题是/否、流失/不流失、点击/不点击需要输出概率值做后续决策的场景特征和标签之间大致呈单调关系或者经过变换后呈单调关系对模型可解释性要求高的业务比如金融风控、医疗诊断如果你手里的问题恰好是这些类型逻辑回归就是一个性价比极高的起步模型。哪怕后面要上更复杂的模型逻辑回归也值得先跑一遍作为baseline。2. 多变量数据在Matlab里怎么组织才不出乱子写模型之前数据组织是第一步。很多人模型跑不通不是算法问题而是数据格式一开始就没搞对。2.1 特征矩阵和标签向量的组织规范在Matlab里逻辑回归的数据输入一般是两个变量X样本特征矩阵维度为 m × nm是样本数n是特征数y标签向量维度为 m × 1取值一般是0和1一个容易踩的坑是特征矩阵的维度方向搞反。Excel表格导进来的数据通常是行是样本、列是特征这个和Matlab的习惯是一致的。但有些数据是从Python转过来的或者用reshape时没注意转置之后维度就错了。建议每次加载完数据先看一眼sizesize(X) size(y)确认X是 m×n、y是 m×1 之后再往下走别偷懒。2.2 特征归一化为什么不能直接拿原始数据训练多变量数据最大的问题就是特征量纲不一致。举个特别典型的例子预测房价特征里面积可能是80到200的数值房龄是1到30的数值而周边学校排名是1到500的数值。如果不做任何处理梯度下降的时候数值大的特征会主导梯度导致模型训练极其不稳定甚至不收敛。我在Matlab里最常用的归一化方式是z-score标准化公式是[ x \frac{x - \mu}{\sigma} ]即每个特征减去均值再除以标准差。处理后的特征均值为0标准差为1。Matlab里用zscore函数一行搞定[X_norm, mu, sigma] zscore(X);注意这里返回了mu和sigma它们是训练集上的均值和标准差。等到模型上线预测新数据时必须用训练集保存下来的mu和sigma对新数据做同样的变换不能重新计算新数据的均值和标准差否则数据分布就变了模型输出就乱套了。还有另一种常用方法是min-max归一化把数据压到[0,1]区间公式是[ x \frac{x - min}{max - min} ]Matlab本身有个mapminmax函数但它是按行操作的和特征按列的常规用法不太一致。我实测下来如果X是 m×n 矩阵用mapminmax得先转置用完再转回来稍不注意就把维度搞乱了。所以我的习惯是大多数情况直接用zscore省心且可控。2.3 训练集和测试集划分的原则模型训练之前一定要把数据分成训练集和测试集。这个原则看起来简单但实际做的时候有讲究。第一划分要在归一化之前还是之后答案是先划分再归一化。如果先整体归一化再划分测试集的信息就已经混进了训练集这叫数据泄露。具体表现是测试集上表现虚高到真实场景就拉胯。第二划分要保持类别比例。假设你的正样本占总体的20%那么训练集和测试集里也应该大致是20%的正样本。在Matlab里最稳妥的方式是按类别分层抽样% 假设y是0/1标签向量 pos_idx find(y 1); neg_idx find(y 0); % 每类各取80%做训练20%做测试 pos_perm randperm(length(pos_idx)); neg_perm randperm(length(neg_idx)); pos_train_idx pos_idx(pos_perm(1:round(0.8*length(pos_idx)))); pos_test_idx pos_idx(pos_perm(round(0.8*length(pos_idx))1:end)); % 负样本同理这样做的好处是保证训练集和测试集的正负样本比例接近模型评估才靠谱。如果直接randperm全体样本索引运气不好的时候测试集里可能全是负样本准确率看着90%实际上啥也没学到。3. 手写梯度下降从数学公式到Matlab代码核心算法部分来了。逻辑回归的训练过程本质上就是调整参数theta让模型的预测结果尽量接近真实标签。这个尽量接近怎么量化需要定义一个损失函数。3.1 交叉熵损失函数是怎么来的逻辑回归用的损失函数是交叉熵形式是[ J(\theta) -\frac{1}{m} \sum_{i1}^{m} \left[ y^{(i)} \log(h^{(i)}) (1 - y^{(i)}) \log(1 - h^{(i)}) \right] ]为什么不直接用均方误差MSE呢原因很简单因为sigmoid函数是S型的如果沿用线性回归的MSE损失整个损失函数会是一个非凸函数里面有大量的局部最小值梯度下降很容易陷进去出不来。而交叉熵损失和sigmoid组合之后整个损失函数是凸函数理论上只有一个全局最小值梯度下降稳稳收敛。这个数学性质是逻辑回归好用的重要基础。从直觉上理解交叉熵当真实标签y1时损失函数只剩下 (- \log(h)) 这一项h越接近1损失越小越接近0损失无穷大。也就是说模型信誓旦旦地预测0但真实是1受到的惩罚极其严重。这种错误越离谱惩罚越大的设计正是分类问题需要的。3.2 梯度计算的向量化写法有了损失函数接下来求梯度。对参数theta_j求偏导最后整理成矩阵形式[ \frac{\partial J}{\partial \theta} \frac{1}{m} X^T (h - y) ]这个结果特别优雅。h是模型输出概率向量y是真实标签向量两者相减得到的是一个误差向量维度 m×1X^T 是 n×m乘出来的结果是 n×1和theta维度一致。每一轮迭代就是theta theta - alpha * (1/m) * X * (h - y);alpha是学习率控制每一步迈多大。为什么梯度是这个形式你可以简单理解为误差(h-y)相当于我错得有多离谱然后通过X^T把这份误差按照每个特征的方向分摊回去告诉每个theta应该往哪个方向调整、调多少。这里不需要手动为每个特征单独写梯度公式矩阵乘法一次性完成了所有参数的更新。3.3 完整的训练函数代码整合上面的内容一个完整的逻辑回归训练函数如下function [theta, loss_history] train_logistic(X, y, alpha, num_iters) % 训练逻辑回归模型 % X: m x n 特征矩阵已归一化 % y: m x 1 标签向量0/1 % alpha: 学习率 % num_iters: 迭代次数 m length(y); % 在X前面加一列1对应截距项theta_0 X [ones(m, 1), X]; theta zeros(size(X, 2), 1); loss_history zeros(num_iters, 1); for iter 1:num_iters % 前向计算 z X * theta; h 1 ./ (1 exp(-z)); % 计算损失交叉熵 loss -1/m * sum(y .* log(h 1e-5) (1 - y) .* log(1 - h 1e-5)); loss_history(iter) loss; % 梯度下降更新 gradient (1/m) * X * (h - y); theta theta - alpha * gradient; end end注意代码里计算对数时加了1e-5这是一个小保护防止h恰好是0或1导致log(0)算出无穷大。这种数值稳定性处理是实战代码和教材代码的重要区别教材里不管这个实际跑起来就会遇到NaN莫名其妙的出现。3.4 学习率和迭代次数的经验取值学习率alpha怎么选我的经验是从0.01开始试。训练完看一眼loss_history是不是在稳定下降figure; plot(loss_history); xlabel(迭代次数); ylabel(损失值); title(训练损失曲线);如果损失曲线震荡或者往上走说明学习率大了减小到0.001再试。如果损失下降太慢5000次迭代还没到底说明学习率小了适当加大到0.05。迭代次数没有一个万能值。有的人固定1000次有的人固定10000次我建议不要死记数字而是看损失曲线变平了就停。更优雅的方案是设置一个容忍度当连续若干次迭代损失下降量小于某个阈值比如1e-5就提前终止训练这样既省时间又避免过拟合。4. 模型效果怎么评估准确率远不够看模型训练完了先别急着开心评估环节比训练更考验功夫。很多人跑完模型只看一个准确率这是远远不够的。4.1 混淆矩阵和相关指标二分类问题的预测结果有四种情况真实正样本被预测为正真正例TP真实正样本被预测为负假负例FN真实负样本被预测为正假正例FP真实负样本被预测为负真负例TN由此衍生出几个核心指标准确率Accuracy (TPTN) / (TPTNFPFN)看整体预测对的占比精确率Precision TP / (TPFP)预测为正的里面有多少是真正正确的召回率Recall TP / (TPFN)真实为正的里面有多少被成功找出来F1分数 2×Precision×Recall / (PrecisionRecall)精确率和召回率的调和平均在Matlab里用confusionmat函数可以直接得到混淆矩阵% 假设pred是预测标签y_test是真实标签 pred double(h_test 0.5); % 默认阈值0.5 C confusionmat(y_test, pred); % C(1,1)是TNC(1,2)是FPC(2,1)是FNC(2,2)是TP为什么说不只看准确率因为类别不平衡的时候准确率会骗人。举个例子100个样本里只有5个正样本模型全部预测成负样本准确率95%但模型实际上什么也没学会。这时候召回率是0一下就暴露问题了。4.2 ROC曲线和AUC值怎么看ROC曲线和AUC值是评估分类模型概率输出质量的金标准。ROC曲线的横轴是假正例率FPR纵轴是真正例率TPR曲线上每一个点对应一个概率阈值。把阈值从1往下扫描到0记录下来每个阈值下的FPR和TPR连线就得到ROC曲线。Matlab里可以用perfcurve函数% h_test是模型输出的概率值 [roc_x, roc_y, ~, auc] perfcurve(y_test, h_test, 1); figure; plot(roc_x, roc_y, b-, LineWidth, 2); xlabel(假正例率 (FPR)); ylabel(真正例率 (TPR)); title([ROC Curve (AUC , num2str(auc), )]);AUC值的含义是随机抽一个正样本和一个负样本模型给正样本打分高于负样本的概率。AUC0.5说明模型和抛硬币没区别AUC0.8以上说明模型有较好的区分能力AUC逼近1说明模型太强但也要警惕过拟合。4.3 概率阈值怎么选才合理逻辑回归输出的概率本身不同业务场景对阈值的要求完全不同。默认0.5是欠考虑的我吃过大亏。如果你做的是医疗筛查类的任务宁可多报假阳性也不愿意漏掉真病人那阈值就应该往下调比如0.3把更多样本预测为正。反过来如果你是做广告点击率预估预算有限希望投放给最有把握的用户那阈值应该上调到0.7左右宁可错过一些潜在点击也要把预算花在刀刃上。阈值调整的实际操作很简单先用perfcurve算出不同阈值下的TPR和FPR再根据业务诉求画一条收益曲线找到收益最大化的点。这个选阈值的过程往往是逻辑回归项目里最体现业务理解的部分。5. 实战踩坑那些让我白加班的Matlab细节做逻辑回归这几个月踩过不少坑挑几个典型的分享出来希望能帮你省下一些本来没必要熬的夜。5.1 不归一化的典型症状有一回我直接把原始数据丢进去训练损失曲线一开始看起来在下降但降到某个点之后开始锯齿状震荡怎么也压不下去。我当时以为是学习率的问题把alpha从0.01调到0.001没用又调到0.0001下降倒是稳定了但速度慢得像蜗牛。后来才反应过来是数据里有个特征取值范围到几千其他特征只有个位数梯度被这个主导特征带着走其他特征几乎学不到东西。zscore归一化之后同样用0.01的学习率几十轮迭代就收敛了。教训如果损失曲线出现诡异的锯齿或者收敛极慢先别急着调学习率检查一下特征归一化做了没有。5.2 数据泄露归一化顺序反了这个坑更隐蔽。我之前按先归一化再划分训练测试集的顺序处理数据测试集上AUC高达0.95当时还挺高兴。一上线模型表现直接崩到0.6。排查了半天终于发现是归一化环节出了问题——整个数据集做zscore时测试集的均值和标准差已经被混进训练过程里了。模型在偷看测试集的信息自然虚高。正确做法是先划分训练集和测试集然后在训练集上计算mu和sigma再用同一组参数去变换测试集。我在前面第2章写代码时特意把这个点拎出来强调因为这个错误特别隐蔽报错也不会有只有上线后才露出马脚。5.3 类别不平衡的时候损失函数直接躺平另一个翻车轮是正负样本比例严重失调。当时手里的数据大概99%是负样本、1%是正样本直接用原始数据做梯度下降训练出来的模型把所有样本都判成负样本AUC只有0.5准确率却高达99%。解决办法有两种思路。一是对少数类做加权在损失函数里给正样本的误差乘以一个放大系数等于告诉模型分错正样本的代价更大。比如正样本权重设为99负样本权重设为1。二是试试对少数类过采样或对多数类欠采样让比例没那么悬殊。二选一或者组合用都可以我的经验是加权损失在Matlab里最好实现改一行代码就够% 假设pos_weight是正样本的惩罚系数 gradient (1/m) * X * (h - y) .* [1; ones(n,1)]; % 具体加权方式需要按样本乘权重更标准的写法是 weight_vec ones(m, 1); weight_vec(y 1) pos_weight; gradient (1/m) * X * ((h - y) .* weight_vec);5.4 手写实现训不出来对比一下内置优化器如果手写的梯度下降怎么调学习率都不收敛还有一个备用方案直接用Matlab内置的fminunc函数让Matlab自己选择最优的下降方式。它的用法需要把损失函数和梯度封装到同一个函数里function [J, grad] costFunction(theta, X, y) m length(y); h 1 ./ (1 exp(-X * theta)); J -1/m * sum(y .* log(h 1e-5) (1 - y) .* log(1 - h 1e-5)); grad (1/m) * X * (h - y); end options optimset(GradObj, on, MaxIter, 400); initial_theta zeros(size(X, 2), 1); [theta_opt] fminunc((t) costFunction(t, X, y), initial_theta, options);fminunc用的是BFGS这类拟牛顿方法比朴素梯度下降收敛快得多还不用手动调学习率。不过它对你的损失函数有个要求必须是光滑可导的。交叉熵加sigmoid完全满足这个条件所以放心用。我现在的习惯是调试阶段用手写梯度下降理解原理真正跑实验用fminunc省时间。两条腿走路既学到了东西又高效。6. 选型思考为什么用逻辑回归而不是BP神经网络做预测同一个预测任务Matlab里也能轻松调用BP神经网络工具箱。所以有一段时间我一直在纠结到底什么时候用逻辑回归什么时候上神经网络。跑了很多对比实验之后总结出了一些经验判断。6.1 两者的本质差异逻辑回归可以看成是一个没有隐藏层的神经网络它的输出层激活函数恰好是sigmoid。BP神经网络相比之下多了一层或者多层隐藏层意味着它能拟合更复杂的非线性关系。这个差异决定了它们的适用边界。如果特征和标签之间的关系本身是近似线性的或者经过特征工程后可以变成近似线性的逻辑回归完全够用而且训练快、参数少、不会过拟合得太离谱。只有当你确信特征之间存在着复杂的非线性交互线性边界实在切不开数据分布时才需要上神经网络。6.2 什么时候逻辑回归更合适我现在的经验是以下几类情况优先选逻辑回归样本量不大几千到几万神经网络很容易过拟合业务方要求模型可解释得能说清楚每个特征对结果的影响方向和幅度模型要频繁上线做实时预测计算资源有限需要一个稳定可靠的baseline做对照逻辑回归的可解释性是它最大的杀手锏。训练完成后theta参数的正负直接告诉你这个特征和预测目标正相关还是负相关大小告诉你影响程度。这在金融风控、医疗辅助诊断这类场景是硬需求哪怕神经网络AUC高上0.02业务方也不敢用一个说不清逻辑的模型做决策。6.3 二分类扩展成多分类预测如果你手里的任务不只是是/否而是要预测多个类别逻辑回归也不是不能干。常见方案是OvR一对多或者Softmax回归。Softmax可以理解为逻辑回归在多分类上的推广输出不再是单个概率而是一个概率向量每个分量对应一个类别的概率。Matlab里实现Softmax回归也不复杂核心是把sigmoid换成softmax函数function h softmax(z) % z是 m x k 矩阵每行是一个样本在所有类别上的得分 exp_z exp(z - max(z, [], 2)); % 减去最大值防止数值溢出 h exp_z ./ sum(exp_z, 2); end这里的技巧是每行减去该行最大值防止exp算出来太大导致Inf。这一步叫数值稳定化很多教程不讲但实际跑起来特别重要。6.4 用交叉验证选模型时的一个细节不管最终用逻辑回归还是神经网络模型选择阶段我强烈建议用K折交叉验证来评估稳定性。Matlab里可以用cvpartitioncv cvpartition(y, KFold, 5); for fold 1:cv.NumTestSets train_idx cv.training(fold); test_idx cv.test(fold); % 在train_idx上训练在test_idx上评估 end这样做的好处是能看出模型在不同数据子集上的表现波动。如果逻辑回归在五折上的AUC分别是0.82、0.81、0.83、0.79、0.84而神经网络是0.79、0.88、0.75、0.86、0.72说明神经网络的稳定性差很多这时候我会更倾向于选择逻辑回归。毕竟在实际部署中模型的不稳定比略低的精度更让人头疼。最后再说一个实用的小细节运行完模型后把关键参数和评估指标存下来。我习惯把theta向量、mu、sigma、测试集上的ROC数据统一保存成一个.mat文件这样下次复现实验结果、调整阈值、或者写技术报告的时候分分钟就能取出来用。模型上线后一旦业务反馈有问题这套可追溯的存档能帮你节省大量排障时间。这就是逻辑回归在实战中比黑盒模型可爱的地方——它所有的逻辑都摊在桌面上你随时能打开来看。
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进