ARTICLE DETAIL

资讯详情

深耕编程入门与网站建设的一线实战洞察。

Logistic回归用于时间序列预测:MATLAB手写实现全流程

Logistic回归用于时间序列预测:MATLAB手写实现全流程 前两天一个朋友甩了个论文味很重的标题给我问我“基于Logistic线性二项分布回归模型的时间序列预测算法这到底是个什么鬼”我第一反应是这不就是换了层学术外衣的Logistic回归嘛。拆开看核心思路其实很简单把时间序列先“翻译”成上升/下降、超阈值/不超阈值这类二项分布事件然后用Logistic回归去估计事件发生的概率最后再把概率反推成预测值。这套做法在金融涨跌预测、工业设备异常预警、用户流失预测里相当常见。它的优势是数据量要求比LSTM低得多输出的是概率而不是一个冷冰冰的点估计而且每个特征对预测结果的影响都可以摊开来看。这篇文章我会把整个流程用MATLAB从零手写一遍——包括滑窗特征构造、数值稳定的Sigmoid函数、梯度下降训练、AUC评估、多步滚动预测不依赖额外工具箱拿到就能跑。1. 先把概念拆清楚线性二项分布回归到底是个什么东西1.1 它的正式名字叫Logistic回归“Logistic线性二项分布回归模型”这名字看着唬人其实就是广义线性模型GLM里的Logistic回归。在统计教材里它经常被描述成“响应变量服从二项分布用logit链接函数建立线性预测项与概率之间的关系”。什么意思呢。假设我们关心某个事件是否发生比如明天的股价比今天高、下一台设备会不会故障、下一位用户会不会流失。这类事件的结果只有0或1单次试验服从的是伯努利分布把它重复多次、统计成功次数那就是标准的二项分布。Logistic回归要做的就是估计这个二项分布中“成功概率”P(y1)到底是多少。它用一个S型曲线——也就是Sigmoid函数——把线性组合压缩到(0,1)区间[ p \frac{1}{1 e^{-z}}, \quad z \beta_0 \beta_1 x_1 \beta_2 x_2 \dots ]换个写法就是logit变换[ \log\left(\frac{p}{1-p}\right) \beta_0 \beta_1 x_1 \beta_2 x_2 \dots ]左边那个对数几率是“事件发生与不发生的比值取对数”。Logistic曲线在这个模型里做的是“闸门”的工作不管z是负一百还是正一百输出永远落在0到1之间。这个性质太重要了因为概率天然就得在这个范围内。1.2 二项分布和时间序列是怎么扯上关系的很多人一听到“时间序列预测”脑子里想的都是ARIMA、指数平滑、LSTM预测出来的都是连续数值。但Logistic回归预测的是概率怎么跟时间序列挂钩答案是把连续序列离散成“状态”。最常见的做法有两种涨跌二值化y(t) 1表示序列在t1时刻比t时刻高0表示低。阈值二值化y(t) 1表示序列在t时刻超过设定阈值0表示没超过。这样一来原本连续的回归问题就变成了一个二分类/二项分布问题。模型学习的不是“下一个值具体是多少”而是“下一个值会不会涨”“会不会超阈值”。这在很多业务场景里其实比精确预测更有用——很多时候你只需要知道方向或者只需要知道风险概率。1.3 与ARIMA、LSTM相比它强在哪里弱在哪里我用一张表说明不同模型在时间序列预测里的定位模型输出形式数据量需求可解释性典型适用场景ARIMA连续值中等较强平稳或差分平稳的线性序列LSTM连续值大较差复杂非线性、长依赖序列Logistic回归概率小很强状态/方向/超阈值预测中小样本LSTM确实能拟合很复杂的时序模式但调参成本高几十条数据根本喂不饱它。ARIMA对非线性趋势和突变处理能力弱。Logistic回归夹在中间适合那些数据量不大、业务上需要解释预测依据、输出概率可以直接对接决策的场景。很多工业预警项目里一条流水线的故障样本可能一年才几百条这种数据上LSTM就是灾难Logistic回归加上人工特征反而非常稳。2. 特征工程时间序列怎么变成回归问题的输入2.1 响应变量Y的构造方式用我之前说的涨跌二值化来构造响应变量。给定一个序列x(1), x(2), …, x(T)定义[ y(t) I\big(x(t1) x(t)\big) ]也就是把“下一时刻是否高于当前时刻”作为预测目标。这里要特别提醒一个新手常犯的错构造y(t)时不能用未来数据偷看。比如你要预测的是t时刻之后的方向特征里只能出现t时刻及之前的信息而y(t)可以用t1时刻的真实值去判断这个顺序不能乱。MATLAB里直接用diff可以快速构造y zeros(T, 1); y(1:end-1) double(diff(x) 0); y(end) y(end-1); % 最后一个点没有未来值沿用上一个2.2 输入特征滑窗滞后项Logistic回归本身不是一个序列模型它不会自动“记住”历史所以我们必须把历史信息加工成特征喂给它。最基础的特征是滑窗滞后项。假设窗口长度是L在第k个时间点取最近L个观测值作为原始特征滞后1到L期的序列值窗口内净变化量窗口最后一个值减去第一个值相当于动量窗口均值起平滑作用如果数据有波动聚集现象可以再加窗口标准差。这些特征组合起来让模型同时看到“当前水平”“变化幅度”和“短期趋势”。如果你有领域知识还可以继续加特征比如同环比、星期几、上一周期的极值这些都不影响整体框架。特征矩阵的构造逻辑是第k行样本用x(k)到x(kL-1)这一段预测第kL时刻的方向。写成代码function [X, Y] buildFeatures(x, L) n length(x); m n - L; X zeros(m, L 2); Y zeros(m, 1); for k 1:m X(k, 1:L) x(k:kL-1); % 滞后L期 X(k, L1) x(kL-1) - x(k); % 窗口净变化 X(k, L2) mean(x(k:kL-1)); % 窗口均值 Y(k) double(x(kL) x(kL-1)); end end2.3 训练测试切分时间序列不能随机打乱这是老生常谈但必须强调普通分类任务习惯用randperm随机打乱数据做交叉验证但时间序列绝对不行。随机打乱等于把未来信息混进训练集让你的验证结果虚高得离谱。正确做法是按时序切分。比如前70%训练、后30%测试split floor(0.7 * size(X, 1)); Xtrain X(1:split, :); Ytrain Y(1:split, :); Xtest X(split1:end, :); Ytest Y(split1:end, :);如果你有耐心可以做滚动窗口验证固定训练窗口长度每预测完一段就往前滚把真实值纳入训练集。这是时序任务里最靠谱的验证方式比一次性切分更能反映模型在真实环境中的表现。3. MATLAB源码级实现从零手写五步流程3.1 生成一份带趋势和季节性的模拟数据为了让大家能把整个流程跑通我先造一份仿真数据。它包含线性趋势、正弦季节成分和随机噪声形态上很像日常业务里的销售、用电负荷序列。rng(42); % 固定随机种子结果可复现 T 600; t (1:T); trend 0.03 .* t; % 线性趋势 season 5 .* sin(0.02 .* t 0.5); % 季节波动 x trend season randn(T, 1) .* 0.6; y zeros(T, 1); y(1:end-1) double(diff(x) 0); y(end) y(end-1);这里T600训练集大约有390个样本测试集接近200个样本足够演示Logistic回归的收敛行为。如果你只有一两百条数据也能跑只是AUC波动会大一些。3.2 标准化为什么这步不能省Sigmoid函数对输入特征的尺度非常敏感。如果特征数值在几百上千的量级而权重初始化在0附近那么z Xβ的值可能极其大或极其小Sigmoid输出直接怼到0或1梯度消失训练直接废掉。正确的预处理是做z-score标准化mu mean(Xtrain); sd std(Xtrain); XtrainS (Xtrain - mu) ./ (sd eps); XtestS (Xtest - mu) ./ (sd eps);注意一个坑测试集标准化必须用训练集的均值和标准差不能用测试集自己算。因为真实预测场景里你拿到的只有当前历史数据未来的均值标准差你是算不出来的。用测试集自己的统计量相当于偷看了未来的分布。3.3 数值稳定的Sigmoid与梯度下降训练Sigmoid函数的标准写法是1 ./ (1 exp(-z))看起来很简洁但z很小比如-800时exp(-z)会溢出成Inf结果变成0。这不是理论问题是实实在在的数值问题。我在工程里习惯拆成两段function p sigmoidStable(z) p zeros(size(z)); pos z 0; p(pos) 1 ./ (1 exp(-z(pos))); % z0 时不会溢出 p(~pos) exp(z(~pos)) ./ (1 exp(z(~pos))); % z0 时 exp(z)不会溢出 end逻辑是z为正时用常规写法z为负时对exp(-z)做等价变换避免指数爆炸。损失函数用交叉熵。因为Logistic回归没有解析解我用批量梯度下降去优化。梯度推导的关键一步是对数似然对β求导后梯度恰好等于X转置乘以“预测概率减真实标签”的平均值形式非常干净function [beta, Jhist] trainLogistic(X, Y, alpha, maxIter, lambda) X [ones(size(X, 1), 1), X]; % 加一列截距项 m size(X, 1); beta zeros(size(X, 2), 1); Jhist zeros(maxIter, 1); for iter 1:maxIter p sigmoidStable(X * beta); grad (X * (p - Y)) / m; if lambda 0 % 可选L2正则防止过拟合 grad(2:end) grad(2:end) lambda * beta(2:end) / m; end beta beta - alpha * grad; J -mean(Y .* log(max(p, eps)) (1-Y) .* log(max(1-p, eps))); if lambda 0 J J lambda * sum(beta(2:end).^2) / (2*m); end Jhist(iter) J; end end调用时我用alpha0.05、迭代1500次、lambda1e-3。标准化之后这个学习率收敛稳定损失曲线在500次迭代之后基本平缓。作为对比如果你装了Statistics and Machine Learning Toolbox可以用fitglm直接验算mdl fitglm(XtrainS, Ytrain, Distribution, binomial);内置函数和手写版得到的权重方向基本一致。但自己写一遍你对“梯度从哪里来”“为什么损失不降了”的感受会完全不同。3.4 输出概率与连续值多步重构模型训练好之后测试集上的预测就是一个概率p表示“下一时刻上涨的概率”。方向判断很直接p大于0.5判上涨小于0.5判下跌。但如果你要的是连续预测值那Logistic只给方向还不够需要做一个幅度重构。我在项目里常用一个启发式方案用历史平均步长作为基准幅度再根据概率偏离0.5的程度缩放。p接近0.5说明模型没把握步长小p越接近0或1说明信号越确定步长越接近历史均值。H 20; stepMean mean(diff(x)); predSeq zeros(H, 1); win x(end-L1:end); % 最新L个值 xPrev x(end); for h 1:H feat win; featRow [feat, feat(end)-feat(1), mean(feat)]; featRowS (featRow - mu) ./ (sd eps); p sigmoidStable([1, featRowS] * beta); dirSign sign(p - 0.5); stepLen abs(stepMean) * (0.5 abs(p - 0.5)); cur xPrev dirSign * stepLen; predSeq(h) cur; win [win(2:end); cur]; xPrev cur; end这段代码里有个关键细节预测出的每个新值会被拼回窗口作为下一次预测的输入。这就是滚动预测也叫递归预测。缺点是误差会累积但优点是代码简单、适合短步数预测。4. 调参与验证关键参数的量级与选择原则4.1 学习率、迭代次数与正则化Logistic回归的损失函数是凸函数理论上梯度下降一定能收敛到全局最优前提是学习率别太大也别太小。在标准化后的数据上我习惯从alpha0.05开始试。判断发散的信号很简单看Jhist曲线。如果损失震荡上升、或者直接变成NaN说明学习率过大。如果损失下降极其缓慢几百次迭代还在稳步下降说明学习率偏小可以加倍。实际项目里我一般看500次迭代后的曲线形态平缓了基本就行。L2正则项lambda取多大我给的1e-3是保守值。特征维度只有7个过拟合风险不大正则影响很小。如果特征维度扩大到几十上百个建议用交叉验证从1e-4到1e-1里选。正则化的本质是给权重加先验——别长得太离谱换来的往往是测试集上的稳定性提升。4.2 滑窗长度L怎么定滑窗L是特征工程里最敏感的参数。L太小模型看不到足够的趋势信息L太大特征维度上升、样本信息被重复利用反而过拟合。我的经验是先用偏自相关函数PACF看序列的自相关截断点。观察原始模拟数据的话滞后5期以内的自相关已经覆盖了大多数信息。所以示例里用L5是合理的。如果数据有明显的周期比如日频业务数据有周周期可以把L加到7或14让窗口覆盖完整周期。更稳妥的做法是设一个候选集合比如{3,5,7,10}每个L跑一遍滚动验证选AUC最高的。4.3 用AUC而不是准确率评估在涨跌二值化任务里准确率是很有欺骗性的指标。假如测试集55%的样本是上涨那么模型只要永远预测“上涨”准确率就有55%看起来不错实际毫无价值。AUCROC曲线下面积对类别不平衡不敏感能更客观地反映排序能力。手写一个快速AUC很方便function auc fastAUC(p, y) [~, I] sort(p, descend); r zeros(size(p)); r(I) 1:numel(p); n1 sum(y 1); n0 numel(y) - n1; if n1 0 || n0 0 auc NaN; return; end auc (sum(r(y 1)) - n1*(n11)/2) / (n1 * n0); end在模拟数据上这个简单特征组合的测试集AUC能到0.62到0.68之间。别小看这个数字——很多日频金融涨跌预测能稳定做到0.55以上就已经有实际价值了。对于纯随机数据AUC会在0.5上下波动跑出来的结果能帮你判断特征到底有没有信息量。5. 常见问题排查与工具链经验长时间用MATLAB做这类模型的训练踩坑主要集中在几个地方我列几个高频问题。问题1训练时损失变成NaN。这是最常见的故障多半是Sigmoid溢出导致的。检查一下代码里有没有直接用1/(1exp(-z))如果有换成稳定版。另外检查特征是否标准化特征尺度差异太大也会导致梯度爆掉。问题2预测结果偏向一边全都预测为0。这是类别不平衡的典型症状。涨跌样本比例如果极端模型会退化成“永远预测多数类”。处理办法有两个。一是用加权交叉熵给少数类样本更高的损失权重二是在预测阶段调整阈值比如用训练集正样本比例作为概率阈值而不是死板的0.5。问题3多步预测越来越平或越来越飘。这是滚动递归预测的误差累积效应。前几步误差小但误差会反馈到窗口里越往后越不可靠。针对这个问题的实用技巧是不要一口气预测几十步而是每预测一步就等待真实值落地把真实值重新放入窗口更新特征。如果必须长步数预测建议每次滚动窗口里保留最近的真实值而不是全用预测值填充。问题4MATLAB版本差异与工具箱缺失。我这套手写代码只用了基本的矩阵运算、sort、mean、std这些基础函数所以没有Statistics Toolbox的机器也能跑。如果你装了工具箱可以用fitglm和perfcurve对照。R2023a之后MATLAB默认字符编码切换到了UTF-8早期用GBK保存的中文注释脚本打开会乱码。解决办法是用编辑器打开后“另存为”选UTF-8编码。这个坑我在老项目脚本上碰过不止一次尤其从旧电脑迁移代码时最容易踩。问题5老是有“为什么不用现成函数”的疑问。fitglm一行确实能出结果但我仍然建议至少手写一遍训练过程。Logistic回归的代码实现量并不大手写能让你彻底搞懂Sigmoid的数值陷阱、梯度怎么算、正则怎么加。一旦上手后面看任何复杂模型的损失函数和优化脚本都不会心虚。6. 这个模型的适用边界与可扩展方向6.1 适合与不适合的场景从我的实际项目经验来看Logistic回归做时间序列预测最适合三种场景第一样本量不大的二值状态预测比如设备故障、用户流失第二业务上需要概率输出的场合比如风控评分里预测逾期概率第三需要向非技术背景领导解释特征影响力的场合Logistic的系数可以直观换算成风险评分。不适合的场景也很明确如果序列存在强非线性衰减、突变模式或者你需要高精度的连续数值预测Logistic回归只能给你方向和粗糙幅度这时候要么换回归模型要么把Logistic作为上游状态判定、再叠加一个误差纠正模块。6.2 从二分类扩展到有序多分类如果你觉得“涨跌”太粗糙可以把序列按照涨跌幅分成多个档位大跌、小跌、持平、小涨、大涨。这就是有序Logistic回归Ordinal Logistic Regression。MATLAB里可以分别构造多个二分类模型做one-vs-rest也可以用mnrfit配合累积logit。思路跟二分类完全一致只是在链接函数上加了一个累积概率的约束。6.3 后续还可以怎么玩我觉得这个模型最有意思的扩展方向是“概率校准加决策优化”。Logistic回归输出的概率虽然落在0到1但在真实业务里往往不够准确常用的校准方法是Platt缩放或Isotonic回归。校准之后再设置决策阈值可以从单纯的预测变成“制定行动策略”——比如预测概率超过0.7才触发报警低于0.3才允许自动通过。最后再分享一个我个人的使用习惯任何时序预测模型上线之前我都要先跑一遍滚动验证而不是用前70%训练、后30%测试这种一次性切分。时间序列最怕跨时间段的分布漂移一次性切分只验证了“同一个历史区间内的泛化能力”滚动验证才能模拟真实部署时“每走一步就更新一次”的状态。这个模型虽然代码量不大但它把Logistic回归从评分卡工具搬到了时序预测场景核心价值在于概率输出和可解释性。跑通之后建议你把滑窗长度、特征组合、正则化系数做一轮网格搜索大概率还能再往上提几个点的AUC。
返回列表