1. 项目背景与核心价值为什么用混合模型预测COVID数量最近几年无论是做研究还是做项目但凡涉及到时间序列预测尤其是像COVID-19新增病例、住院人数这类数据大家都会发现一个头疼的问题数据太“妖”了。它既有明显的趋势性比如一波波的疫情高峰又有很强的季节性比如周末检测量下降还时不时给你来个突变新毒株出现、防控政策调整。传统的单一模型比如ARIMA对付平稳序列还行遇到这种复杂情况预测曲线跟真实数据比起来经常是“差之毫厘谬以千里”。所以混合模型Hybrid Model的思路就火起来了。核心思想很简单一个模型搞不定那就找几个各有所长的模型“组团”让它们优势互补。这次我们要聊的这个项目标题是“基于BILSTM, GMDH和遗传神经网络实现COVID数量的回归预测”就是一个非常典型的、思路清晰的混合模型实战案例。它把深度学习的序列建模能力BILSTM、数据驱动的自组织网络GMDH和进化算法的优化能力遗传算法优化的神经网络拧成了一股绳。这个组合拳的价值在哪里我结合自己处理类似公共卫生数据的经验来说说。首先BILSTM双向长短期记忆网络是处理时间序列的利器它能同时考虑过去和未来的信息来理解当前状态对于捕捉病例数据前后期的依赖关系非常有效。但LSTM类模型有个老毛病容易过拟合特别是在数据量不是特别巨大的时候很多地区的COVID数据也就几百到几千条而且超参数层数、神经元数、学习率调起来很费劲效果不稳定。这时候GMDHGroup Method of Data Handling数据处理分组方法就派上用场了。它本质上是一种迭代的自组织算法能够从一堆可能的模型结构中自动筛选出复杂度适中、预测性能最好的那个。你可以把它理解为一个“模型结构自动搜索器”。把它和神经网络结合就是GMDH神经网络它能帮我们找到一个相对更优的网络架构避免人工设计网络的盲目性。那“遗传神经网络”又是什么角色你可以把它看作是GMDH神经网络的“高级教练”。GMDH确定了网络的大致结构有几层层之间怎么连但每层有多少个神经元、神经元之间的连接权重是多少这些参数还是需要优化的。用标准的梯度下降法训练容易陷入局部最优解。遗传算法Genetic Algorithm模仿生物进化通过选择、交叉、变异来全局搜索最优参数组合正好能弥补这个不足让神经网络的“大脑”变得更聪明。所以这个项目的核心逻辑链条是这样的用BILSTM从原始时间序列中提取高层次的特征和时序依赖关系将这些特征作为输入送入一个由GMDH方法初步构建、再由遗传算法精细优化参数的前馈神经网络也就是遗传神经网络中进行最终的非线性回归预测。BILSTM负责“理解”时间模式GMDH负责“设计”预测器结构遗传算法负责“打磨”预测器参数三者各司其职最终目标就是得到一个对COVID数量变化更稳健、更准确的预测模型。对于公共卫生决策、医疗资源调度来说哪怕预测精度提升几个百分点都可能带来显著的社会效益。接下来我们就一步步拆解如何在MATLAB里把这个想法实现出来。2. 核心组件深度解析BILSTM、GMDH与遗传算法如何协同在动手写代码之前我们必须把这三个核心组件的原理、在项目中的具体职责以及它们之间如何传递数据搞清楚。一知半解地调包出了问题根本无从下手。2.1 BILSTM双向的时序特征提取器LSTM大家都很熟了通过门控机制遗忘门、输入门、输出门来解决长期依赖问题。而BILSTMBidirectional LSTM则是在此基础上增加了一个反向层。正向LSTM层按时间顺序t1, 2, ..., T处理序列反向LSTM层按逆序tT, ..., 2, 1处理序列。在每一个时间步t模型会将正向层和反向层在该时刻的隐藏状态连接起来作为最终的输出。为什么在疫情预测中需要双向因为病例数往往不是孤立的。今天的病例数不仅受过去几天的影响传播周期也可能隐含着对未来几天趋势的“预示”例如检测阳性率上升可能预示着未来几日报告病例数将增加。BILSTM通过结合过去和未来的上下文信息能更全面地“理解”当前时间点在整个疫情曲线中的位置。在这个项目中BILSTM扮演的角色是“特征提取器”或“编码器”。它的任务不是直接输出预测值而是将原始的、可能包含多个特征如每日新增病例、检测数、移动平均等的时序数据转换为一组固定长度的、富含时序信息的特征向量。这个特征向量将作为下游GMDH-遗传神经网络的输入。一个关键细节是BILSTM的输出维度。假设我们设置BILSTM层有NumHiddenUnits个神经元由于是双向每个时间步的输出维度将是2 * NumHiddenUnits。我们通常取最后一个时间步对应序列末尾的输出或者对所有时间步的输出做全局平均池化来得到一个1 x (2*NumHiddenUnits)的特征向量。这个向量就封装了整个输入序列的时序模式。2.2 GMDH自组织的网络结构设计师GMDH是一种归纳式的建模方法它的核心思想是“迭代筛选”和“外推准则”。它从一个非常简单的模型比如两个变量的二次多项式开始通过不断组合产生更复杂的候选模型并使用一个独立的验证集来评估这些候选模型的性能只保留性能最好的那些进入下一轮迭代直到模型复杂度增加不再带来预测精度的提升。将其应用于神经网络时GMDH用于确定网络的最佳结构应该有几层隐藏层每层应该有多少个神经元神经元之间如何连接传统的做法是手动试错或者凭经验而GMDH将这个流程自动化了。在项目中GMDH的作用是“结构搜索”。它的输入是BILSTM提取出的特征向量假设维度为M以及对应的目标值未来N天的病例数。GMDH算法会尝试多种不同的神经元组合方式例如第一层尝试5、10、15个神经元第二层尝试3、5个神经元等在验证集上评估每种结构对应的网络预测效果最终选择一个在验证集上误差最小的网络结构。这个结构就是一个未经训练或仅初步训练的前馈神经网络。注意纯GMDH算法本身可以构建一个完整的多项式网络。但在与神经网络结合的语境下更常见的做法是用GMDH思想来优化一个普通前馈神经网络的结构。MATLAB中没有直接的“GMDH神经网络”工具箱我们需要自己实现这个结构搜索过程或者利用其迭代选择的思想。2.3 遗传算法神经网络的终极调参教练遗传算法是受生物进化论启发的一类优化算法。它维护一个“种群”种群中的每个“个体”代表问题的一个可能解在这里就是神经网络的权重和偏置参数集。算法通过模拟自然选择选择适应度高的个体、交叉交换两个个体的部分基因产生后代和变异随机改变个体的部分基因来迭代进化种群最终找到适应度最高预测误差最小的个体即最优的网络参数。在项目中遗传算法扮演“参数优化器”的角色。当GMDH为我们选定了一个神经网络结构比如输入层-M个节点隐藏层1-10个节点隐藏层2-5个节点输出层-1个节点后这个网络有大量的权重和偏置参数需要确定。如果用反向传播BP算法容易陷入局部最小值且对学习率等超参数敏感。遗传算法则从全局出发能更有效地找到一组更优的参数从而让这个既定结构的网络发挥出最大潜力。三者的协作流程可以概括为数据准备将COVID时间序列数据整理成监督学习格式例如用过去14天的数据预测未来7天的数据。BILSTM特征提取用训练集训练一个BILSTM网络然后用它处理所有数据训练、验证、测试将每个样本序列转换为一个特征向量。GMDH结构搜索在BILSTM特征向量和对应目标值构成的新的数据集上运行GMDH算法或模拟其思想确定一个最优的前馈神经网络结构。遗传算法参数优化固定上一步确定的结构使用遗传算法在这个结构上搜索最优的权重和偏置参数。预测与评估使用优化好的“BILSTM特征提取 GMDH结构遗传算法参数”的混合模型对测试集进行预测并评估性能。理解了这套逻辑代码实现就有了清晰的路线图。下面我们进入实战环节。3. MATLAB环境准备与数据预处理实战工欲善其事必先利其器。在写核心算法之前我们必须把环境和数据准备好。这里我会分享一些在MATLAB中做时间序列预测时容易踩坑的细节。3.1 工具箱与依赖检查这个项目需要用到几个关键的MATLAB工具箱Deep Learning Toolbox用于构建和训练BILSTM网络。这是核心依赖。Global Optimization Toolbox或自定义遗传算法代码用于实现遗传算法优化。MATLAB的Global Optimization Toolbox提供了ga函数可以直接调用非常方便。Statistics and Machine Learning Toolbox用于数据标准化、划分数据集等常规操作。在脚本开头可以通过以下代码检查并给出提示% 检查必要工具箱 requiredToolboxes {Deep Learning Toolbox, Global Optimization Toolbox, Statistics and Machine Learning Toolbox}; for i 1:length(requiredToolboxes) if ~license(test, requiredToolboxes{i}) || isempty(ver(requiredToolboxes{i})) error(请确保安装并激活 %s 工具箱。, requiredToolboxes{i}); end end disp(必要工具箱检查通过。);3.2 COVID-19数据获取与清洗数据源可以选择公开数据集如约翰斯·霍普金斯大学JHU的GitHub仓库、Our World in Data等。这里假设我们已经从一个CSV文件中加载了数据包含date和new_cases两列。关键步骤1处理缺失值与异常值疫情数据常因报告延迟等原因出现负值或异常大值。简单的处理方法是% 假设 data 是一个timetable或table包含Date和Cases变量 % 1. 将负值替换为NaN或前后插值 data.Cases(data.Cases 0) NaN; % 2. 使用移动中位数或Hampel滤波器去除脉冲噪声 window 7; % 7天窗口 data.Cases hampel(data.Cases, window); % hampel滤波能识别并替换离群点 % 3. 对剩余的NaN进行线性插值 data.Cases fillmissing(data.Cases, linear);实操心得对于病例数据我更喜欢用hampel滤波器而不是简单的移动平均。因为它能更鲁棒地识别并处理那些由于单日数据录入错误产生的极端离群点同时更好地保留真实的波动趋势。关键步骤2构建监督学习数据集我们要用过去lookback天的数据预测未来horizon天的数据。这是一个多步预测问题。lookback 14; % 用过去14天预测 horizon 7; % 预测未来7天 cases data.Cases; % 一维病例数列 numSteps length(cases); X []; % 特征序列 Y []; % 目标序列 for i 1:(numSteps - lookback - horizon 1) X(end1, :) cases(i:ilookback-1); % 输入过去lookback天 Y(end1, :) cases(ilookback:ilookbackhorizon-1); % 输出未来horizon天 end现在X是一个m x lookback的矩阵Y是m x horizon的矩阵其中m是样本数。关键步骤3数据标准化与数据集划分时序数据必须按时间顺序划分不能随机打乱% 按时间顺序划分70%训练15%验证15%测试 mTotal size(X, 1); idxTrain 1:floor(0.7 * mTotal); idxVal floor(0.7*mTotal)1 : floor(0.85*mTotal); idxTest floor(0.85*mTotal)1 : mTotal; XTrain X(idxTrain, :); YTrain Y(idxTrain, :); XVal X(idxVal, :); YVal Y(idxVal, :); XTest X(idxTest, :); YTest Y(idxTest, :); % 标准化仅使用训练集的均值和标准差 muX mean(XTrain, 1); sigmaX std(XTrain, 0, 1); muY mean(YTrain, 1); sigmaY std(YTrain, 0, 1); XTrainNorm (XTrain - muX) ./ sigmaX; XValNorm (XVal - muX) ./ sigmaX; XTestNorm (XTest - muX) ./ sigmaX; YTrainNorm (YTrain - muY) ./ sigmaY; % 注意验证集和测试集的目标值标准化仅用于训练过程中的比较最终评估需反标准化踩坑提醒标准化时一定要用训练集的统计量这是一个非常容易犯的错误。用整个数据集包含未来信息的均值和标准差来标准化会导致数据泄露严重高估模型性能。验证集和测试集必须使用从训练集计算得到的mu和sigma。数据准备好了下一步就是搭建并训练我们的BILSTM特征提取器。4. BILSTM特征提取器的构建、训练与特征生成在这一步我们的目标是训练一个BILSTM网络让它学会从一段病例序列中提取出有代表性的特征。这个网络本身不直接输出预测值而是输出一个特征向量。4.1 网络架构设计与层配置在MATLAB的Deep Learning Toolbox中构建序列网络非常直观。对于回归任务一个典型的BILSTM网络结构如下inputSize 1; % 输入特征维度我们只有病例数这一个特征但可以扩展为多变量 numHiddenUnits 64; % LSTM隐藏单元数这是一个可调超参数 numFeatures 2 * numHiddenUnits; % BILSTM的输出维度 layers [ sequenceInputLayer(inputSize, Name, input) % 序列输入层 bilstmLayer(numHiddenUnits, OutputMode, last, Name, bilstm) % 关键OutputMode设为last只取最后一个时间步的输出 fullyConnectedLayer(numFeatures, Name, fc_feature) % 全连接层将BILSTM输出映射到特征空间 regressionLayer(Name, output) % 回归层用于训练 ];这里有几个关键设计点bilstmLayer的OutputMode参数我们设置为last。这意味着对于每个输入序列BILSTM层只输出最后一个时间步的隐藏状态连接了正向和反向。这个1 x (2*numHiddenUnits)的向量就是我们想要的特征向量。如果设置为sequence则会输出每个时间步的状态后续处理会更复杂。全连接层fullyConnectedLayer这一层不是必须的但加上它有两个好处。一是可以增加非线性让特征表达更丰富二是可以强制将特征向量的维度固定为我们想要的numFeatures。numFeatures是一个超参数可以小于或等于2*numHiddenUnits起到降维或保持的作用。训练目标注意这个网络是用未来horizon天的真实值标准化后的YTrainNorm作为目标来训练的。虽然我们最终不直接用它的输出做预测但通过这个回归任务训练可以迫使BILSTM学会提取与未来趋势高度相关的特征。4.2 训练选项与技巧训练循环神经网络需要仔细设置选项。options trainingOptions(adam, ... % 优化器 MaxEpochs, 150, ... % 最大迭代次数 MiniBatchSize, 32, ... % 批大小 InitialLearnRate, 0.001, ... % 初始学习率 GradientThreshold, 1, ... % 梯度阈值防止梯度爆炸 Shuffle, every-epoch, ... % 每个epoch打乱数据 ValidationData, {XValNorm, YValNorm}, ... % 验证集 ValidationFrequency, 30, ... % 每30次迭代验证一次 Verbose, false, ... % 不显示详细训练过程 Plots, training-progress, ...% 绘制训练进度图 LearnRateSchedule, piecewise, ... LearnRateDropFactor, 0.5, ... LearnRateDropPeriod, 80); % 学习率衰减Shuffle, every-epoch对于时间序列通常不建议打乱因为会破坏时序依赖。但在这里我们的每个样本一段历史序列已经是独立的打乱有助于提升训练效果防止模型记忆顺序。学习率衰减这是一个非常实用的技巧。训练后期降低学习率有助于模型收敛到更优的点。4.3 训练网络与提取特征训练完成后我们需要“拆解”这个网络取出特征提取部分。% 训练网络 netBILSTM trainNetwork(XTrainNorm, YTrainNorm, layers, options); % 提取特征提取器去掉最后的回归层 featureLayers netBILSTM.Layers(1:end-1); % 取除最后regressionLayer外的所有层 netFeature assembleNetwork(featureLayers); % 重新组装成网络 % 使用特征提取器处理所有数据 XTrainFeature predict(netFeature, XTrainNorm); XValFeature predict(netFeature, XValNorm); XTestFeature predict(netFeature, XTestNorm);现在XTrainFeature是一个size(XTrainNorm, 1) x numFeatures的矩阵。原始的时序数据被转换成了高维特征空间中的点。这个特征矩阵就是后续GMDH和遗传神经网络的输入。经验之谈训练BILSTM时如果验证集损失很早就开始上升而训练集损失持续下降这是典型的过拟合。可以尝试1) 增加Dropout层在bilstmLayer后加dropoutLayer(0.5)2) 减少numHiddenUnits3) 使用更早停止ValidationPatience, 20。我们的目标是让BILSTM学到通用的时序模式而不是记住训练集噪声。5. 模拟GMDH思想进行神经网络结构搜索如前所述MATLAB没有现成的GMDH神经网络工具箱。因此我们需要模拟GMDH的核心思想——通过迭代评估来自动选择网络结构。这里我提供一个简化但实用的实现方案使用网格搜索Grid Search结合验证集性能来寻找一个较优的前馈网络结构。5.1 定义结构搜索空间我们假设GMDH要为我们确定一个包含至多两个隐藏层的网络结构。我们定义一些候选的神经元数量。% 结构搜索空间 hidden1_candidates [8, 16, 32, 64]; % 第一隐藏层神经元数候选 hidden2_candidates [0, 4, 8, 16]; % 第二隐藏层神经元数候选0表示没有该层 activation_candidates {relu, tanh}; % 激活函数候选 numFeatures size(XTrainFeature, 2); % BILSTM提取的特征维度 outputSize horizon; % 预测未来horizon天 bestValRMSE inf; bestArchitecture struct();我们用一个双重循环来遍历所有可能的组合。5.2 网格搜索与验证评估对于每一种结构我们快速训练一个浅层网络使用标准的trainlm或trainscg算法并在验证集上评估其性能。for h1 hidden1_candidates for h2 hidden2_candidates if h2 0 % 单隐藏层结构 layers [featureInputLayer(numFeatures) fullyConnectedLayer(h1) reluLayer() % 或根据候选选择 fullyConnectedLayer(outputSize) regressionLayer]; else % 双隐藏层结构 layers [featureInputLayer(numFeatures) fullyConnectedLayer(h1) reluLayer() fullyConnectedLayer(h2) reluLayer() fullyConnectedLayer(outputSize) regressionLayer]; end % 转换为layerGraph以兼容trainNetwork如果使用Deep Learning Toolbox % 但这里为了快速评估我们使用传统的fitrnet需要Statistics and Machine Learning Toolbox % 它更轻量适合快速结构搜索。 % 方案A使用fitrnet (R2021a以上) % net fitrnet(XTrainFeature, YTrainNorm, LayerSizes, [h1, h2(h2~0)], Activations, relu); % YValPred predict(net, XValFeature); % 方案B使用trainNetwork但设置少量epochs快速训练 optionsFast trainingOptions(sgdm, ... MaxEpochs, 50, ... MiniBatchSize, 32, ... InitialLearnRate, 0.01, ... Verbose, false, ... Plots, none); % 注意需要将数据转换为适合序列输入的格式这里我们使用featureInputLayer % 因此需要将数据以n*1的cell数组形式传入每个cell是一个样本的特征向量 XTrainCell num2cell(XTrainFeature, 1); % 转置并转换使每列为一个样本 XValCell num2cell(XValFeature, 1); % 但trainNetwork对featureInputLayer的支持在回归任务上可能不如fitrnet直接。 % 鉴于结构搜索是独立步骤我推荐使用更简单的交叉验证配合浅层网络。 % 这里给出一个更通用的实现使用MATLAB的feedforwardnet并手动训练 net feedforwardnet([h1, h2(h2~0)]); % 创建网络 net.trainFcn trainlm; % Levenberg-Marquardt训练快 net.performFcn mse; % 均方误差 net.divideFcn divideind; % 按索引划分 net.divideParam.trainInd 1:size(XTrainFeature,1); net.divideParam.valInd []; net.divideParam.testInd []; net.trainParam.epochs 50; % 快速训练 net.trainParam.showWindow false; % 训练并预测 [net, tr] train(net, XTrainFeature, YTrainNorm); YValPred net(XValFeature); % 计算验证集RMSE需反标准化到原始量纲进行比较 YValPredOrig YValPred .* sigmaY muY; YValOrig YVal .* sigmaY muY; % YVal是未标准化的原始目标值 currentRMSE sqrt(mean((YValPredOrig - YValOrig).^2, all)); fprintf(结构 [%d, %d] - Val RMSE: %.2f\n, h1, h2, currentRMSE); if currentRMSE bestValRMSE bestValRMSE currentRMSE; bestArchitecture.hidden1 h1; bestArchitecture.hidden2 h2; bestArchitecture.net net; % 保存当前最佳网络 end end end fprintf(\n最佳结构: 隐藏层1 %d, 隐藏层2 %d, 验证集RMSE %.2f\n, ... bestArchitecture.hidden1, bestArchitecture.hidden2, bestValRMSE);这个网格搜索过程模拟了GMDH的“评估-选择”思想。我们遍历了多种可能的结构用验证集性能作为选择依据最终找到了一个在验证集上表现最好的网络结构bestArchitecture。重要提示这个搜索过程计算量可能较大。在实际项目中可以根据经验先缩小搜索范围或者使用更高效的超参数优化方法如贝叶斯优化。这里为了清晰展示GMDH思想采用了最直观的网格搜索。6. 遗传算法优化神经网络权重现在我们有了一个确定的网络结构例如输入层 - 64个神经元 - 16个神经元 - 输出层。接下来我们要用遗传算法来优化这个网络的权重和偏置而不是用传统的反向传播。6.1 将网络参数编码为染色体遗传算法处理的是“染色体”即一维向量。我们需要将神经网络的所有权重W和偏置b展平并拼接成一个长向量。% 假设 bestArchitecture.net 是我们找到的最佳结构网络 net bestArchitecture.net; % 获取当前网络的初始权重和偏置 IW net.IW{1,1}; % 输入层到第一隐藏层的权重 b1 net.b{1}; % 第一隐藏层的偏置 LW1 net.LW{2,1}; % 第一隐藏层到第二隐藏层的权重如果存在 b2 net.b{2}; % 第二隐藏层的偏置 LW2 net.LW{3,2}; % 第二隐藏层到输出层的权重如果存在 b3 net.b{3}; % 输出层的偏置 % 将它们展平并拼接成一个向量 initialParams [IW(:); b1(:); LW1(:); b2(:); LW2(:); b3(:)];这个initialParams向量就是我们的“染色体”。遗传算法的任务就是找到使网络预测误差最小的染色体。6.2 定义适应度函数适应度函数是遗传算法的核心它评估一条染色体即一组网络参数的好坏。对于回归问题我们通常用均方误差MSE或均方根误差RMSE的负数作为适应度因为遗传算法默认是求最大值我们需要最小化误差。function fitness gaFitness(params, net, XTrain, YTrain, XVal, YVal) % params: 染色体向量 % net: 网络模板用于获取结构信息 % XTrain, YTrain: 训练数据用于计算训练误差可加入正则化 % XVal, YVal: 验证数据用于计算适应度防止过拟合 % 1. 将染色体解码回网络参数 [decodedNet, paramInfo] decodeChromosome(params, net); % 2. 使用解码后的网络进行预测 YValPred decodedNet(XVal); % 3. 计算验证集上的均方根误差RMSE mse mean((YValPred - YVal).^2, all); rmse sqrt(mse); % 4. 适应度 -RMSE 遗传算法最大化适应度所以我们最小化RMSE fitness -rmse; % 可选加入正则化项L2正则化到适应度函数中惩罚大的权重防止过拟合 % lambda 0.001; % weightPenalty lambda * sum(params.^2); % fitness -rmse - weightPenalty; end其中decodeChromosome是一个自定义函数负责根据网络结构信息将一维参数向量重新赋值给网络的各个权重矩阵和偏置向量。这是遗传算法优化神经网络最关键也最繁琐的一步需要仔细处理矩阵的维度。6.3 配置与运行遗传算法MATLAB的Global Optimization Toolbox提供了ga函数。我们需要配置其选项。% 定义优化问题维度染色体长度 nvars length(initialParams); % 定义变量的上下界权重和偏置的范围 lb -10 * ones(1, nvars); % 下界 ub 10 * ones(1, nvars); % 上界 % 设置遗传算法选项 optionsGA optimoptions(ga, ... PopulationSize, 100, ... % 种群大小 MaxGenerations, 200, ... % 最大进化代数 FunctionTolerance, 1e-6, ... % 函数值容忍度 PlotFcn, gaplotbestf, ... % 绘制最佳适应度曲线 Display, iter, ... % 显示迭代信息 UseParallel, true); % 如果可用使用并行计算加速 % 运行遗传算法 % 注意适应度函数需要能接收单个参数向量。我们通过匿名函数将其他参数固定。 fitnessFunc (params) gaFitness(params, net, XTrainFeature, YTrainNorm, XValFeature, YValNorm); [bestParams, bestFitness] ga(fitnessFunc, nvars, [], [], [], [], lb, ub, [], optionsGA); fprintf(遗传算法优化完成。最佳适应度负RMSE: %.4f\n, bestFitness);遗传算法运行时间可能较长取决于种群大小、代数和网络规模。UseParallel, true可以显著加速计算。6.4 解码最优参数并构建最终模型得到最优染色体bestParams后我们将其解码回网络得到最终优化好的GMDH-遗传神经网络。% 解码最优参数 optimizedNet decodeChromosome(bestParams, net); % 现在 optimizedNet 就是我们的最终预测模型第二部分 % 完整的混合模型是BILSTM特征提取器 (netFeature) 优化后的前馈网络 (optimizedNet)至此我们完成了混合模型的所有组件构建和优化。接下来就是将它们串联起来进行最终的预测和评估。7. 模型集成、预测与性能评估现在我们有了两个训练好的组件netFeatureBILSTM特征提取器和optimizedNetGMDH结构遗传算法优化的前馈网络。预测流程是新数据 - 标准化 - BILSTM提取特征 - 优化网络预测 - 反标准化。7.1 构建集成预测函数function YPred hybridPredict(netFeature, optimizedNet, X, muX, sigmaX, muY, sigmaY) % X: 原始输入数据维度 m x lookback % netFeature: BILSTM特征提取网络 % optimizedNet: 优化后的前馈网络 % muX, sigmaX: 用于标准化的训练集统计量 % muY, sigmaY: 用于反标准化的训练集统计量 % 1. 标准化输入 XNorm (X - muX) ./ sigmaX; % 2. BILSTM提取特征 (注意predict函数可能要求特定输入格式) % 如果netFeature是SeriesNetwork或DAGNetworkpredict接受数值数组 XFeatures predict(netFeature, XNorm); % 3. 优化网络预测 YPredNorm optimizedNet(XFeatures); % 注意网络输入是特征按列排列 % 4. 反标准化预测结果 YPred YPredNorm .* sigmaY muY; end7.2 在测试集上进行最终评估使用测试集模型从未见过的数据来公平地评估混合模型的泛化能力。% 使用集成模型进行测试集预测 YTestPred hybridPredict(netFeature, optimizedNet, XTest, muX, sigmaX, muY, sigmaY); % 计算关键评估指标 testRMSE sqrt(mean((YTestPred - YTest).^2, all)); testMAE mean(abs(YTestPred - YTest), all); testMAPE mean(abs((YTestPred - YTest) ./ YTest), all) * 100; % 平均绝对百分比误差注意YTest可能为0需处理 fprintf(测试集性能评估:\n); fprintf(RMSE: %.2f\n, testRMSE); fprintf(MAE: %.2f\n, testMAE); fprintf(MAPE: %.2f%%\n, testMAPE); % 为了对比我们也可以评估单独使用BILSTM直接回归的性能 % YTestPred_BILSTM predict(netBILSTM, XTestNorm); % YTestPred_BILSTM_Orig YTestPred_BILSTM .* sigmaY muY; % testRMSE_BILSTM sqrt(mean((YTestPred_BILSTM_Orig - YTest).^2, all)); % fprintf(单独BILSTM测试集RMSE: %.2f\n, testRMSE_BILSTM);7.3 结果可视化与分析可视化是理解模型表现的关键。% 绘制预测值与真实值对比图以第一个预测步长为例或绘制多步预测的序列 figure(Position, [100, 100, 1200, 500]); subplot(1,2,1); plotIdx 1:min(100, size(YTest,1)); % 绘制前100个测试样本 plot(YTest(plotIdx, 1), b-o, LineWidth, 1.5, DisplayName, 真实值); hold on; plot(YTestPred(plotIdx, 1), r--s, LineWidth, 1.5, DisplayName, 预测值); xlabel(测试样本序号); ylabel(病例数); title(测试集预测对比 (未来第1天)); legend(Location, best); grid on; % 绘制预测误差分布直方图 subplot(1,2,2); errors YTestPred - YTest; histogram(errors(:,1), 50); % 查看未来第1天预测误差的分布 xlabel(预测误差); ylabel(频数); title(预测误差分布); grid on; % 计算并绘制整个预测视界horizon的平均误差 meanErrorOverHorizon mean(abs(YTestPred - YTest), 1); figure; plot(1:horizon, meanErrorOverHorizon, m-^, LineWidth, 2); xlabel(预测天数 (未来第n天)); ylabel(平均绝对误差 (MAE)); title(预测误差随预测天数的变化); grid on;通过分析这些图表我们可以判断预测趋势是否一致对比图应显示预测曲线能大致跟随真实曲线的波动。误差分布误差直方图应近似正态分布且均值接近0说明无系统性偏差。预测不确定性误差通常随着预测天数增加而增大这符合常识。7.4 与基线模型对比为了体现混合模型的优势务必与一些基线模型对比例如朴素预测用最近一天的值预测未来所有天持久化预测。线性回归直接用历史数据对每个未来天做线性回归。单一BILSTM仅使用我们训练的第一个BILSTM网络带回归层进行预测。标准前馈网络用相同的BILSTM特征但网络结构随机设定并用BP算法训练。制作一个对比表格会非常直观模型测试集RMSE测试集MAE备注朴素预测计算值计算值基准线线性回归计算值计算值线性模型单一BILSTM计算值计算值深度时序模型标准前馈网络计算值计算值人工设计结构BPBILSTMGMDHGA本模型计算值计算值混合模型如果混合模型的RMSE和MAE显著低于其他模型就证明了这种集成策略的有效性。8. 项目总结、避坑指南与扩展思考走完整个流程你会发现这个项目虽然思路清晰但实现起来细节颇多很容易在某个环节卡住。这里我总结几个最容易出问题的地方和对应的解决方案。避坑指南1数据标准化与信息泄露这是最隐蔽也最致命的错误。一定要牢记任何从数据中计算的统计量均值、标准差、最小值、最大值都必须且只能从训练集计算。验证集和测试集必须使用训练集的统计量进行转换。在代码中要把muX、sigmaX、muY、sigmaY作为从训练集“学来”的参数保存好并在预测新数据时严格使用它们。避坑指南2BILSTM的输入输出格式MATLAB的Deep Learning Toolbox对输入格式要求严格。对于序列数据训练时XTrainNorm通常应该是numFeatures x numTimeSteps x numSamples的三维数组对于单变量numFeatures1或者是numTimeSteps x numSamples的二维数组此时numFeatures在sequenceInputLayer中指定。而在使用predict函数提取特征时要确保输入格式与训练时一致。如果遇到维度错误仔细检查size并使用permute、reshape等函数调整。避坑指南3遗传算法的适应度函数设计适应度函数是遗传算法的导航仪。除了验证集误差强烈建议加入正则化项如L2正则化来惩罚大的权重这能有效防止过拟合提升模型泛化能力。另外适应度函数的计算开销很大每次都要前向传播整个验证集如果数据量大可以考虑在适应度函数中对验证集进行采样或者使用更小的验证子集来加速。避坑指南4GMDH结构搜索的耗时我们实现的网格搜索是“暴力”方法。如果搜索空间大如隐藏层数、每层神经元数、激活函数类型都作为候选计算量会爆炸。在实际应用中可以考虑分阶段搜索先粗搜如神经元数以2的幂次增长再在最优点附近细搜。使用更高效的优化器MATLAB的bayesopt函数贝叶斯优化非常适合这类黑盒函数优化问题能更快地找到较优结构。借鉴经典网络结构根据经验对于这类回归问题1-2个隐藏层每层神经元数在输入维度和输出维度之间通常效果就不错。可以大大缩小搜索范围。扩展思考1引入更多特征本项目只使用了历史病例数。实际上COVID-19的传播受多种因素影响。可以考虑加入外部特征疫苗接种率、移动指数、政府干预指数如牛津字符串ency指数、天气数据温度、湿度。衍生特征病例数的7天移动平均、环比增长率、加速度等。 这些特征可以作为多变量输入给BILSTM调整inputSize让模型捕捉更复杂的关系。扩展思考2概率预测与不确定性量化点预测给出一个具体值往往不够。决策者更关心预测的置信区间。可以对模型进行扩展例如让神经网络的输出层输出均值和方差使用负对数似然损失。使用蒙特卡洛DropoutMC Dropout在预测时多次推理用输出的分布来估计不确定性。采用分位数回归直接预测不同分位数如10%50%90%对应的值。扩展思考3模型部署与更新模型不是一劳永逸的。疫情数据分布会随时间变化概念漂移。在实际部署中需要建立模型监控和定期更新机制。例如可以每周用最新数据对模型进行微调Fine-tuning或者设定一个阈值当预测误差连续超过阈值时触发模型重训练。这个项目从构思到实现完整地展示了一个混合机器学习模型解决实际时间序列预测问题的 pipeline。它不仅仅是代码的堆砌更是对问题本质的思考为什么单一模型不够、对工具特性的理解每个组件擅长什么以及对工程细节的把握数据怎么处理、参数怎么调。希望这份详细的拆解能让你在复现或借鉴这个思路时少走一些弯路更深入地理解模型背后的“为什么”和“怎么做”。 混合模型预测COVID-19病例：BILSTM、GMDH与遗传算法实战
1. 项目背景与核心价值为什么用混合模型预测COVID数量最近几年无论是做研究还是做项目但凡涉及到时间序列预测尤其是像COVID-19新增病例、住院人数这类数据大家都会发现一个头疼的问题数据太“妖”了。它既有明显的趋势性比如一波波的疫情高峰又有很强的季节性比如周末检测量下降还时不时给你来个突变新毒株出现、防控政策调整。传统的单一模型比如ARIMA对付平稳序列还行遇到这种复杂情况预测曲线跟真实数据比起来经常是“差之毫厘谬以千里”。所以混合模型Hybrid Model的思路就火起来了。核心思想很简单一个模型搞不定那就找几个各有所长的模型“组团”让它们优势互补。这次我们要聊的这个项目标题是“基于BILSTM, GMDH和遗传神经网络实现COVID数量的回归预测”就是一个非常典型的、思路清晰的混合模型实战案例。它把深度学习的序列建模能力BILSTM、数据驱动的自组织网络GMDH和进化算法的优化能力遗传算法优化的神经网络拧成了一股绳。这个组合拳的价值在哪里我结合自己处理类似公共卫生数据的经验来说说。首先BILSTM双向长短期记忆网络是处理时间序列的利器它能同时考虑过去和未来的信息来理解当前状态对于捕捉病例数据前后期的依赖关系非常有效。但LSTM类模型有个老毛病容易过拟合特别是在数据量不是特别巨大的时候很多地区的COVID数据也就几百到几千条而且超参数层数、神经元数、学习率调起来很费劲效果不稳定。这时候GMDHGroup Method of Data Handling数据处理分组方法就派上用场了。它本质上是一种迭代的自组织算法能够从一堆可能的模型结构中自动筛选出复杂度适中、预测性能最好的那个。你可以把它理解为一个“模型结构自动搜索器”。把它和神经网络结合就是GMDH神经网络它能帮我们找到一个相对更优的网络架构避免人工设计网络的盲目性。那“遗传神经网络”又是什么角色你可以把它看作是GMDH神经网络的“高级教练”。GMDH确定了网络的大致结构有几层层之间怎么连但每层有多少个神经元、神经元之间的连接权重是多少这些参数还是需要优化的。用标准的梯度下降法训练容易陷入局部最优解。遗传算法Genetic Algorithm模仿生物进化通过选择、交叉、变异来全局搜索最优参数组合正好能弥补这个不足让神经网络的“大脑”变得更聪明。所以这个项目的核心逻辑链条是这样的用BILSTM从原始时间序列中提取高层次的特征和时序依赖关系将这些特征作为输入送入一个由GMDH方法初步构建、再由遗传算法精细优化参数的前馈神经网络也就是遗传神经网络中进行最终的非线性回归预测。BILSTM负责“理解”时间模式GMDH负责“设计”预测器结构遗传算法负责“打磨”预测器参数三者各司其职最终目标就是得到一个对COVID数量变化更稳健、更准确的预测模型。对于公共卫生决策、医疗资源调度来说哪怕预测精度提升几个百分点都可能带来显著的社会效益。接下来我们就一步步拆解如何在MATLAB里把这个想法实现出来。2. 核心组件深度解析BILSTM、GMDH与遗传算法如何协同在动手写代码之前我们必须把这三个核心组件的原理、在项目中的具体职责以及它们之间如何传递数据搞清楚。一知半解地调包出了问题根本无从下手。2.1 BILSTM双向的时序特征提取器LSTM大家都很熟了通过门控机制遗忘门、输入门、输出门来解决长期依赖问题。而BILSTMBidirectional LSTM则是在此基础上增加了一个反向层。正向LSTM层按时间顺序t1, 2, ..., T处理序列反向LSTM层按逆序tT, ..., 2, 1处理序列。在每一个时间步t模型会将正向层和反向层在该时刻的隐藏状态连接起来作为最终的输出。为什么在疫情预测中需要双向因为病例数往往不是孤立的。今天的病例数不仅受过去几天的影响传播周期也可能隐含着对未来几天趋势的“预示”例如检测阳性率上升可能预示着未来几日报告病例数将增加。BILSTM通过结合过去和未来的上下文信息能更全面地“理解”当前时间点在整个疫情曲线中的位置。在这个项目中BILSTM扮演的角色是“特征提取器”或“编码器”。它的任务不是直接输出预测值而是将原始的、可能包含多个特征如每日新增病例、检测数、移动平均等的时序数据转换为一组固定长度的、富含时序信息的特征向量。这个特征向量将作为下游GMDH-遗传神经网络的输入。一个关键细节是BILSTM的输出维度。假设我们设置BILSTM层有NumHiddenUnits个神经元由于是双向每个时间步的输出维度将是2 * NumHiddenUnits。我们通常取最后一个时间步对应序列末尾的输出或者对所有时间步的输出做全局平均池化来得到一个1 x (2*NumHiddenUnits)的特征向量。这个向量就封装了整个输入序列的时序模式。2.2 GMDH自组织的网络结构设计师GMDH是一种归纳式的建模方法它的核心思想是“迭代筛选”和“外推准则”。它从一个非常简单的模型比如两个变量的二次多项式开始通过不断组合产生更复杂的候选模型并使用一个独立的验证集来评估这些候选模型的性能只保留性能最好的那些进入下一轮迭代直到模型复杂度增加不再带来预测精度的提升。将其应用于神经网络时GMDH用于确定网络的最佳结构应该有几层隐藏层每层应该有多少个神经元神经元之间如何连接传统的做法是手动试错或者凭经验而GMDH将这个流程自动化了。在项目中GMDH的作用是“结构搜索”。它的输入是BILSTM提取出的特征向量假设维度为M以及对应的目标值未来N天的病例数。GMDH算法会尝试多种不同的神经元组合方式例如第一层尝试5、10、15个神经元第二层尝试3、5个神经元等在验证集上评估每种结构对应的网络预测效果最终选择一个在验证集上误差最小的网络结构。这个结构就是一个未经训练或仅初步训练的前馈神经网络。注意纯GMDH算法本身可以构建一个完整的多项式网络。但在与神经网络结合的语境下更常见的做法是用GMDH思想来优化一个普通前馈神经网络的结构。MATLAB中没有直接的“GMDH神经网络”工具箱我们需要自己实现这个结构搜索过程或者利用其迭代选择的思想。2.3 遗传算法神经网络的终极调参教练遗传算法是受生物进化论启发的一类优化算法。它维护一个“种群”种群中的每个“个体”代表问题的一个可能解在这里就是神经网络的权重和偏置参数集。算法通过模拟自然选择选择适应度高的个体、交叉交换两个个体的部分基因产生后代和变异随机改变个体的部分基因来迭代进化种群最终找到适应度最高预测误差最小的个体即最优的网络参数。在项目中遗传算法扮演“参数优化器”的角色。当GMDH为我们选定了一个神经网络结构比如输入层-M个节点隐藏层1-10个节点隐藏层2-5个节点输出层-1个节点后这个网络有大量的权重和偏置参数需要确定。如果用反向传播BP算法容易陷入局部最小值且对学习率等超参数敏感。遗传算法则从全局出发能更有效地找到一组更优的参数从而让这个既定结构的网络发挥出最大潜力。三者的协作流程可以概括为数据准备将COVID时间序列数据整理成监督学习格式例如用过去14天的数据预测未来7天的数据。BILSTM特征提取用训练集训练一个BILSTM网络然后用它处理所有数据训练、验证、测试将每个样本序列转换为一个特征向量。GMDH结构搜索在BILSTM特征向量和对应目标值构成的新的数据集上运行GMDH算法或模拟其思想确定一个最优的前馈神经网络结构。遗传算法参数优化固定上一步确定的结构使用遗传算法在这个结构上搜索最优的权重和偏置参数。预测与评估使用优化好的“BILSTM特征提取 GMDH结构遗传算法参数”的混合模型对测试集进行预测并评估性能。理解了这套逻辑代码实现就有了清晰的路线图。下面我们进入实战环节。3. MATLAB环境准备与数据预处理实战工欲善其事必先利其器。在写核心算法之前我们必须把环境和数据准备好。这里我会分享一些在MATLAB中做时间序列预测时容易踩坑的细节。3.1 工具箱与依赖检查这个项目需要用到几个关键的MATLAB工具箱Deep Learning Toolbox用于构建和训练BILSTM网络。这是核心依赖。Global Optimization Toolbox或自定义遗传算法代码用于实现遗传算法优化。MATLAB的Global Optimization Toolbox提供了ga函数可以直接调用非常方便。Statistics and Machine Learning Toolbox用于数据标准化、划分数据集等常规操作。在脚本开头可以通过以下代码检查并给出提示% 检查必要工具箱 requiredToolboxes {Deep Learning Toolbox, Global Optimization Toolbox, Statistics and Machine Learning Toolbox}; for i 1:length(requiredToolboxes) if ~license(test, requiredToolboxes{i}) || isempty(ver(requiredToolboxes{i})) error(请确保安装并激活 %s 工具箱。, requiredToolboxes{i}); end end disp(必要工具箱检查通过。);3.2 COVID-19数据获取与清洗数据源可以选择公开数据集如约翰斯·霍普金斯大学JHU的GitHub仓库、Our World in Data等。这里假设我们已经从一个CSV文件中加载了数据包含date和new_cases两列。关键步骤1处理缺失值与异常值疫情数据常因报告延迟等原因出现负值或异常大值。简单的处理方法是% 假设 data 是一个timetable或table包含Date和Cases变量 % 1. 将负值替换为NaN或前后插值 data.Cases(data.Cases 0) NaN; % 2. 使用移动中位数或Hampel滤波器去除脉冲噪声 window 7; % 7天窗口 data.Cases hampel(data.Cases, window); % hampel滤波能识别并替换离群点 % 3. 对剩余的NaN进行线性插值 data.Cases fillmissing(data.Cases, linear);实操心得对于病例数据我更喜欢用hampel滤波器而不是简单的移动平均。因为它能更鲁棒地识别并处理那些由于单日数据录入错误产生的极端离群点同时更好地保留真实的波动趋势。关键步骤2构建监督学习数据集我们要用过去lookback天的数据预测未来horizon天的数据。这是一个多步预测问题。lookback 14; % 用过去14天预测 horizon 7; % 预测未来7天 cases data.Cases; % 一维病例数列 numSteps length(cases); X []; % 特征序列 Y []; % 目标序列 for i 1:(numSteps - lookback - horizon 1) X(end1, :) cases(i:ilookback-1); % 输入过去lookback天 Y(end1, :) cases(ilookback:ilookbackhorizon-1); % 输出未来horizon天 end现在X是一个m x lookback的矩阵Y是m x horizon的矩阵其中m是样本数。关键步骤3数据标准化与数据集划分时序数据必须按时间顺序划分不能随机打乱% 按时间顺序划分70%训练15%验证15%测试 mTotal size(X, 1); idxTrain 1:floor(0.7 * mTotal); idxVal floor(0.7*mTotal)1 : floor(0.85*mTotal); idxTest floor(0.85*mTotal)1 : mTotal; XTrain X(idxTrain, :); YTrain Y(idxTrain, :); XVal X(idxVal, :); YVal Y(idxVal, :); XTest X(idxTest, :); YTest Y(idxTest, :); % 标准化仅使用训练集的均值和标准差 muX mean(XTrain, 1); sigmaX std(XTrain, 0, 1); muY mean(YTrain, 1); sigmaY std(YTrain, 0, 1); XTrainNorm (XTrain - muX) ./ sigmaX; XValNorm (XVal - muX) ./ sigmaX; XTestNorm (XTest - muX) ./ sigmaX; YTrainNorm (YTrain - muY) ./ sigmaY; % 注意验证集和测试集的目标值标准化仅用于训练过程中的比较最终评估需反标准化踩坑提醒标准化时一定要用训练集的统计量这是一个非常容易犯的错误。用整个数据集包含未来信息的均值和标准差来标准化会导致数据泄露严重高估模型性能。验证集和测试集必须使用从训练集计算得到的mu和sigma。数据准备好了下一步就是搭建并训练我们的BILSTM特征提取器。4. BILSTM特征提取器的构建、训练与特征生成在这一步我们的目标是训练一个BILSTM网络让它学会从一段病例序列中提取出有代表性的特征。这个网络本身不直接输出预测值而是输出一个特征向量。4.1 网络架构设计与层配置在MATLAB的Deep Learning Toolbox中构建序列网络非常直观。对于回归任务一个典型的BILSTM网络结构如下inputSize 1; % 输入特征维度我们只有病例数这一个特征但可以扩展为多变量 numHiddenUnits 64; % LSTM隐藏单元数这是一个可调超参数 numFeatures 2 * numHiddenUnits; % BILSTM的输出维度 layers [ sequenceInputLayer(inputSize, Name, input) % 序列输入层 bilstmLayer(numHiddenUnits, OutputMode, last, Name, bilstm) % 关键OutputMode设为last只取最后一个时间步的输出 fullyConnectedLayer(numFeatures, Name, fc_feature) % 全连接层将BILSTM输出映射到特征空间 regressionLayer(Name, output) % 回归层用于训练 ];这里有几个关键设计点bilstmLayer的OutputMode参数我们设置为last。这意味着对于每个输入序列BILSTM层只输出最后一个时间步的隐藏状态连接了正向和反向。这个1 x (2*numHiddenUnits)的向量就是我们想要的特征向量。如果设置为sequence则会输出每个时间步的状态后续处理会更复杂。全连接层fullyConnectedLayer这一层不是必须的但加上它有两个好处。一是可以增加非线性让特征表达更丰富二是可以强制将特征向量的维度固定为我们想要的numFeatures。numFeatures是一个超参数可以小于或等于2*numHiddenUnits起到降维或保持的作用。训练目标注意这个网络是用未来horizon天的真实值标准化后的YTrainNorm作为目标来训练的。虽然我们最终不直接用它的输出做预测但通过这个回归任务训练可以迫使BILSTM学会提取与未来趋势高度相关的特征。4.2 训练选项与技巧训练循环神经网络需要仔细设置选项。options trainingOptions(adam, ... % 优化器 MaxEpochs, 150, ... % 最大迭代次数 MiniBatchSize, 32, ... % 批大小 InitialLearnRate, 0.001, ... % 初始学习率 GradientThreshold, 1, ... % 梯度阈值防止梯度爆炸 Shuffle, every-epoch, ... % 每个epoch打乱数据 ValidationData, {XValNorm, YValNorm}, ... % 验证集 ValidationFrequency, 30, ... % 每30次迭代验证一次 Verbose, false, ... % 不显示详细训练过程 Plots, training-progress, ...% 绘制训练进度图 LearnRateSchedule, piecewise, ... LearnRateDropFactor, 0.5, ... LearnRateDropPeriod, 80); % 学习率衰减Shuffle, every-epoch对于时间序列通常不建议打乱因为会破坏时序依赖。但在这里我们的每个样本一段历史序列已经是独立的打乱有助于提升训练效果防止模型记忆顺序。学习率衰减这是一个非常实用的技巧。训练后期降低学习率有助于模型收敛到更优的点。4.3 训练网络与提取特征训练完成后我们需要“拆解”这个网络取出特征提取部分。% 训练网络 netBILSTM trainNetwork(XTrainNorm, YTrainNorm, layers, options); % 提取特征提取器去掉最后的回归层 featureLayers netBILSTM.Layers(1:end-1); % 取除最后regressionLayer外的所有层 netFeature assembleNetwork(featureLayers); % 重新组装成网络 % 使用特征提取器处理所有数据 XTrainFeature predict(netFeature, XTrainNorm); XValFeature predict(netFeature, XValNorm); XTestFeature predict(netFeature, XTestNorm);现在XTrainFeature是一个size(XTrainNorm, 1) x numFeatures的矩阵。原始的时序数据被转换成了高维特征空间中的点。这个特征矩阵就是后续GMDH和遗传神经网络的输入。经验之谈训练BILSTM时如果验证集损失很早就开始上升而训练集损失持续下降这是典型的过拟合。可以尝试1) 增加Dropout层在bilstmLayer后加dropoutLayer(0.5)2) 减少numHiddenUnits3) 使用更早停止ValidationPatience, 20。我们的目标是让BILSTM学到通用的时序模式而不是记住训练集噪声。5. 模拟GMDH思想进行神经网络结构搜索如前所述MATLAB没有现成的GMDH神经网络工具箱。因此我们需要模拟GMDH的核心思想——通过迭代评估来自动选择网络结构。这里我提供一个简化但实用的实现方案使用网格搜索Grid Search结合验证集性能来寻找一个较优的前馈网络结构。5.1 定义结构搜索空间我们假设GMDH要为我们确定一个包含至多两个隐藏层的网络结构。我们定义一些候选的神经元数量。% 结构搜索空间 hidden1_candidates [8, 16, 32, 64]; % 第一隐藏层神经元数候选 hidden2_candidates [0, 4, 8, 16]; % 第二隐藏层神经元数候选0表示没有该层 activation_candidates {relu, tanh}; % 激活函数候选 numFeatures size(XTrainFeature, 2); % BILSTM提取的特征维度 outputSize horizon; % 预测未来horizon天 bestValRMSE inf; bestArchitecture struct();我们用一个双重循环来遍历所有可能的组合。5.2 网格搜索与验证评估对于每一种结构我们快速训练一个浅层网络使用标准的trainlm或trainscg算法并在验证集上评估其性能。for h1 hidden1_candidates for h2 hidden2_candidates if h2 0 % 单隐藏层结构 layers [featureInputLayer(numFeatures) fullyConnectedLayer(h1) reluLayer() % 或根据候选选择 fullyConnectedLayer(outputSize) regressionLayer]; else % 双隐藏层结构 layers [featureInputLayer(numFeatures) fullyConnectedLayer(h1) reluLayer() fullyConnectedLayer(h2) reluLayer() fullyConnectedLayer(outputSize) regressionLayer]; end % 转换为layerGraph以兼容trainNetwork如果使用Deep Learning Toolbox % 但这里为了快速评估我们使用传统的fitrnet需要Statistics and Machine Learning Toolbox % 它更轻量适合快速结构搜索。 % 方案A使用fitrnet (R2021a以上) % net fitrnet(XTrainFeature, YTrainNorm, LayerSizes, [h1, h2(h2~0)], Activations, relu); % YValPred predict(net, XValFeature); % 方案B使用trainNetwork但设置少量epochs快速训练 optionsFast trainingOptions(sgdm, ... MaxEpochs, 50, ... MiniBatchSize, 32, ... InitialLearnRate, 0.01, ... Verbose, false, ... Plots, none); % 注意需要将数据转换为适合序列输入的格式这里我们使用featureInputLayer % 因此需要将数据以n*1的cell数组形式传入每个cell是一个样本的特征向量 XTrainCell num2cell(XTrainFeature, 1); % 转置并转换使每列为一个样本 XValCell num2cell(XValFeature, 1); % 但trainNetwork对featureInputLayer的支持在回归任务上可能不如fitrnet直接。 % 鉴于结构搜索是独立步骤我推荐使用更简单的交叉验证配合浅层网络。 % 这里给出一个更通用的实现使用MATLAB的feedforwardnet并手动训练 net feedforwardnet([h1, h2(h2~0)]); % 创建网络 net.trainFcn trainlm; % Levenberg-Marquardt训练快 net.performFcn mse; % 均方误差 net.divideFcn divideind; % 按索引划分 net.divideParam.trainInd 1:size(XTrainFeature,1); net.divideParam.valInd []; net.divideParam.testInd []; net.trainParam.epochs 50; % 快速训练 net.trainParam.showWindow false; % 训练并预测 [net, tr] train(net, XTrainFeature, YTrainNorm); YValPred net(XValFeature); % 计算验证集RMSE需反标准化到原始量纲进行比较 YValPredOrig YValPred .* sigmaY muY; YValOrig YVal .* sigmaY muY; % YVal是未标准化的原始目标值 currentRMSE sqrt(mean((YValPredOrig - YValOrig).^2, all)); fprintf(结构 [%d, %d] - Val RMSE: %.2f\n, h1, h2, currentRMSE); if currentRMSE bestValRMSE bestValRMSE currentRMSE; bestArchitecture.hidden1 h1; bestArchitecture.hidden2 h2; bestArchitecture.net net; % 保存当前最佳网络 end end end fprintf(\n最佳结构: 隐藏层1 %d, 隐藏层2 %d, 验证集RMSE %.2f\n, ... bestArchitecture.hidden1, bestArchitecture.hidden2, bestValRMSE);这个网格搜索过程模拟了GMDH的“评估-选择”思想。我们遍历了多种可能的结构用验证集性能作为选择依据最终找到了一个在验证集上表现最好的网络结构bestArchitecture。重要提示这个搜索过程计算量可能较大。在实际项目中可以根据经验先缩小搜索范围或者使用更高效的超参数优化方法如贝叶斯优化。这里为了清晰展示GMDH思想采用了最直观的网格搜索。6. 遗传算法优化神经网络权重现在我们有了一个确定的网络结构例如输入层 - 64个神经元 - 16个神经元 - 输出层。接下来我们要用遗传算法来优化这个网络的权重和偏置而不是用传统的反向传播。6.1 将网络参数编码为染色体遗传算法处理的是“染色体”即一维向量。我们需要将神经网络的所有权重W和偏置b展平并拼接成一个长向量。% 假设 bestArchitecture.net 是我们找到的最佳结构网络 net bestArchitecture.net; % 获取当前网络的初始权重和偏置 IW net.IW{1,1}; % 输入层到第一隐藏层的权重 b1 net.b{1}; % 第一隐藏层的偏置 LW1 net.LW{2,1}; % 第一隐藏层到第二隐藏层的权重如果存在 b2 net.b{2}; % 第二隐藏层的偏置 LW2 net.LW{3,2}; % 第二隐藏层到输出层的权重如果存在 b3 net.b{3}; % 输出层的偏置 % 将它们展平并拼接成一个向量 initialParams [IW(:); b1(:); LW1(:); b2(:); LW2(:); b3(:)];这个initialParams向量就是我们的“染色体”。遗传算法的任务就是找到使网络预测误差最小的染色体。6.2 定义适应度函数适应度函数是遗传算法的核心它评估一条染色体即一组网络参数的好坏。对于回归问题我们通常用均方误差MSE或均方根误差RMSE的负数作为适应度因为遗传算法默认是求最大值我们需要最小化误差。function fitness gaFitness(params, net, XTrain, YTrain, XVal, YVal) % params: 染色体向量 % net: 网络模板用于获取结构信息 % XTrain, YTrain: 训练数据用于计算训练误差可加入正则化 % XVal, YVal: 验证数据用于计算适应度防止过拟合 % 1. 将染色体解码回网络参数 [decodedNet, paramInfo] decodeChromosome(params, net); % 2. 使用解码后的网络进行预测 YValPred decodedNet(XVal); % 3. 计算验证集上的均方根误差RMSE mse mean((YValPred - YVal).^2, all); rmse sqrt(mse); % 4. 适应度 -RMSE 遗传算法最大化适应度所以我们最小化RMSE fitness -rmse; % 可选加入正则化项L2正则化到适应度函数中惩罚大的权重防止过拟合 % lambda 0.001; % weightPenalty lambda * sum(params.^2); % fitness -rmse - weightPenalty; end其中decodeChromosome是一个自定义函数负责根据网络结构信息将一维参数向量重新赋值给网络的各个权重矩阵和偏置向量。这是遗传算法优化神经网络最关键也最繁琐的一步需要仔细处理矩阵的维度。6.3 配置与运行遗传算法MATLAB的Global Optimization Toolbox提供了ga函数。我们需要配置其选项。% 定义优化问题维度染色体长度 nvars length(initialParams); % 定义变量的上下界权重和偏置的范围 lb -10 * ones(1, nvars); % 下界 ub 10 * ones(1, nvars); % 上界 % 设置遗传算法选项 optionsGA optimoptions(ga, ... PopulationSize, 100, ... % 种群大小 MaxGenerations, 200, ... % 最大进化代数 FunctionTolerance, 1e-6, ... % 函数值容忍度 PlotFcn, gaplotbestf, ... % 绘制最佳适应度曲线 Display, iter, ... % 显示迭代信息 UseParallel, true); % 如果可用使用并行计算加速 % 运行遗传算法 % 注意适应度函数需要能接收单个参数向量。我们通过匿名函数将其他参数固定。 fitnessFunc (params) gaFitness(params, net, XTrainFeature, YTrainNorm, XValFeature, YValNorm); [bestParams, bestFitness] ga(fitnessFunc, nvars, [], [], [], [], lb, ub, [], optionsGA); fprintf(遗传算法优化完成。最佳适应度负RMSE: %.4f\n, bestFitness);遗传算法运行时间可能较长取决于种群大小、代数和网络规模。UseParallel, true可以显著加速计算。6.4 解码最优参数并构建最终模型得到最优染色体bestParams后我们将其解码回网络得到最终优化好的GMDH-遗传神经网络。% 解码最优参数 optimizedNet decodeChromosome(bestParams, net); % 现在 optimizedNet 就是我们的最终预测模型第二部分 % 完整的混合模型是BILSTM特征提取器 (netFeature) 优化后的前馈网络 (optimizedNet)至此我们完成了混合模型的所有组件构建和优化。接下来就是将它们串联起来进行最终的预测和评估。7. 模型集成、预测与性能评估现在我们有了两个训练好的组件netFeatureBILSTM特征提取器和optimizedNetGMDH结构遗传算法优化的前馈网络。预测流程是新数据 - 标准化 - BILSTM提取特征 - 优化网络预测 - 反标准化。7.1 构建集成预测函数function YPred hybridPredict(netFeature, optimizedNet, X, muX, sigmaX, muY, sigmaY) % X: 原始输入数据维度 m x lookback % netFeature: BILSTM特征提取网络 % optimizedNet: 优化后的前馈网络 % muX, sigmaX: 用于标准化的训练集统计量 % muY, sigmaY: 用于反标准化的训练集统计量 % 1. 标准化输入 XNorm (X - muX) ./ sigmaX; % 2. BILSTM提取特征 (注意predict函数可能要求特定输入格式) % 如果netFeature是SeriesNetwork或DAGNetworkpredict接受数值数组 XFeatures predict(netFeature, XNorm); % 3. 优化网络预测 YPredNorm optimizedNet(XFeatures); % 注意网络输入是特征按列排列 % 4. 反标准化预测结果 YPred YPredNorm .* sigmaY muY; end7.2 在测试集上进行最终评估使用测试集模型从未见过的数据来公平地评估混合模型的泛化能力。% 使用集成模型进行测试集预测 YTestPred hybridPredict(netFeature, optimizedNet, XTest, muX, sigmaX, muY, sigmaY); % 计算关键评估指标 testRMSE sqrt(mean((YTestPred - YTest).^2, all)); testMAE mean(abs(YTestPred - YTest), all); testMAPE mean(abs((YTestPred - YTest) ./ YTest), all) * 100; % 平均绝对百分比误差注意YTest可能为0需处理 fprintf(测试集性能评估:\n); fprintf(RMSE: %.2f\n, testRMSE); fprintf(MAE: %.2f\n, testMAE); fprintf(MAPE: %.2f%%\n, testMAPE); % 为了对比我们也可以评估单独使用BILSTM直接回归的性能 % YTestPred_BILSTM predict(netBILSTM, XTestNorm); % YTestPred_BILSTM_Orig YTestPred_BILSTM .* sigmaY muY; % testRMSE_BILSTM sqrt(mean((YTestPred_BILSTM_Orig - YTest).^2, all)); % fprintf(单独BILSTM测试集RMSE: %.2f\n, testRMSE_BILSTM);7.3 结果可视化与分析可视化是理解模型表现的关键。% 绘制预测值与真实值对比图以第一个预测步长为例或绘制多步预测的序列 figure(Position, [100, 100, 1200, 500]); subplot(1,2,1); plotIdx 1:min(100, size(YTest,1)); % 绘制前100个测试样本 plot(YTest(plotIdx, 1), b-o, LineWidth, 1.5, DisplayName, 真实值); hold on; plot(YTestPred(plotIdx, 1), r--s, LineWidth, 1.5, DisplayName, 预测值); xlabel(测试样本序号); ylabel(病例数); title(测试集预测对比 (未来第1天)); legend(Location, best); grid on; % 绘制预测误差分布直方图 subplot(1,2,2); errors YTestPred - YTest; histogram(errors(:,1), 50); % 查看未来第1天预测误差的分布 xlabel(预测误差); ylabel(频数); title(预测误差分布); grid on; % 计算并绘制整个预测视界horizon的平均误差 meanErrorOverHorizon mean(abs(YTestPred - YTest), 1); figure; plot(1:horizon, meanErrorOverHorizon, m-^, LineWidth, 2); xlabel(预测天数 (未来第n天)); ylabel(平均绝对误差 (MAE)); title(预测误差随预测天数的变化); grid on;通过分析这些图表我们可以判断预测趋势是否一致对比图应显示预测曲线能大致跟随真实曲线的波动。误差分布误差直方图应近似正态分布且均值接近0说明无系统性偏差。预测不确定性误差通常随着预测天数增加而增大这符合常识。7.4 与基线模型对比为了体现混合模型的优势务必与一些基线模型对比例如朴素预测用最近一天的值预测未来所有天持久化预测。线性回归直接用历史数据对每个未来天做线性回归。单一BILSTM仅使用我们训练的第一个BILSTM网络带回归层进行预测。标准前馈网络用相同的BILSTM特征但网络结构随机设定并用BP算法训练。制作一个对比表格会非常直观模型测试集RMSE测试集MAE备注朴素预测计算值计算值基准线线性回归计算值计算值线性模型单一BILSTM计算值计算值深度时序模型标准前馈网络计算值计算值人工设计结构BPBILSTMGMDHGA本模型计算值计算值混合模型如果混合模型的RMSE和MAE显著低于其他模型就证明了这种集成策略的有效性。8. 项目总结、避坑指南与扩展思考走完整个流程你会发现这个项目虽然思路清晰但实现起来细节颇多很容易在某个环节卡住。这里我总结几个最容易出问题的地方和对应的解决方案。避坑指南1数据标准化与信息泄露这是最隐蔽也最致命的错误。一定要牢记任何从数据中计算的统计量均值、标准差、最小值、最大值都必须且只能从训练集计算。验证集和测试集必须使用训练集的统计量进行转换。在代码中要把muX、sigmaX、muY、sigmaY作为从训练集“学来”的参数保存好并在预测新数据时严格使用它们。避坑指南2BILSTM的输入输出格式MATLAB的Deep Learning Toolbox对输入格式要求严格。对于序列数据训练时XTrainNorm通常应该是numFeatures x numTimeSteps x numSamples的三维数组对于单变量numFeatures1或者是numTimeSteps x numSamples的二维数组此时numFeatures在sequenceInputLayer中指定。而在使用predict函数提取特征时要确保输入格式与训练时一致。如果遇到维度错误仔细检查size并使用permute、reshape等函数调整。避坑指南3遗传算法的适应度函数设计适应度函数是遗传算法的导航仪。除了验证集误差强烈建议加入正则化项如L2正则化来惩罚大的权重这能有效防止过拟合提升模型泛化能力。另外适应度函数的计算开销很大每次都要前向传播整个验证集如果数据量大可以考虑在适应度函数中对验证集进行采样或者使用更小的验证子集来加速。避坑指南4GMDH结构搜索的耗时我们实现的网格搜索是“暴力”方法。如果搜索空间大如隐藏层数、每层神经元数、激活函数类型都作为候选计算量会爆炸。在实际应用中可以考虑分阶段搜索先粗搜如神经元数以2的幂次增长再在最优点附近细搜。使用更高效的优化器MATLAB的bayesopt函数贝叶斯优化非常适合这类黑盒函数优化问题能更快地找到较优结构。借鉴经典网络结构根据经验对于这类回归问题1-2个隐藏层每层神经元数在输入维度和输出维度之间通常效果就不错。可以大大缩小搜索范围。扩展思考1引入更多特征本项目只使用了历史病例数。实际上COVID-19的传播受多种因素影响。可以考虑加入外部特征疫苗接种率、移动指数、政府干预指数如牛津字符串ency指数、天气数据温度、湿度。衍生特征病例数的7天移动平均、环比增长率、加速度等。 这些特征可以作为多变量输入给BILSTM调整inputSize让模型捕捉更复杂的关系。扩展思考2概率预测与不确定性量化点预测给出一个具体值往往不够。决策者更关心预测的置信区间。可以对模型进行扩展例如让神经网络的输出层输出均值和方差使用负对数似然损失。使用蒙特卡洛DropoutMC Dropout在预测时多次推理用输出的分布来估计不确定性。采用分位数回归直接预测不同分位数如10%50%90%对应的值。扩展思考3模型部署与更新模型不是一劳永逸的。疫情数据分布会随时间变化概念漂移。在实际部署中需要建立模型监控和定期更新机制。例如可以每周用最新数据对模型进行微调Fine-tuning或者设定一个阈值当预测误差连续超过阈值时触发模型重训练。这个项目从构思到实现完整地展示了一个混合机器学习模型解决实际时间序列预测问题的 pipeline。它不仅仅是代码的堆砌更是对问题本质的思考为什么单一模型不够、对工具特性的理解每个组件擅长什么以及对工程细节的把握数据怎么处理、参数怎么调。希望这份详细的拆解能让你在复现或借鉴这个思路时少走一些弯路更深入地理解模型背后的“为什么”和“怎么做”。