基于Matlab的LSTM神经网络时序预测(完整代码+范例数据文件)
本文部分理論引用了:4 Strategies for Multi-Step Time Series Forecasting
Multivariate Time Series Forecasting with LSTMs in Keras (machinelearningmastery.com)
(更新于2022年10月,最近需要做相關(guān)的大型實驗,對這個LSTM神經(jīng)網(wǎng)絡(luò)預(yù)測有了更深的認(rèn)識,也發(fā)現(xiàn)以往有些概念是錯誤的,在這里記錄一下。有問題或者有其他預(yù)測需求可以直接私聊我一起交流,我有空的時候會看的。)
(1.首先題目的“時序多步預(yù)測”改為了“時序預(yù)測”,其實這個代碼算是一個多次的單步預(yù)測,而且并不需要一次性去預(yù)測兩步(如果需要改的話,需要改神經(jīng)網(wǎng)絡(luò)預(yù)測中的交替時間步,可以用一個XTrain對應(yīng)兩個YTrain中的兩個數(shù),那就可以多步預(yù)測了。原理和結(jié)果詳情可見置頂?shù)牡诙恼?#xff0c;里面基于Python提到了一次性預(yù)測兩步的問題),因為LSTM不斷在更新網(wǎng)絡(luò),所以這樣反而精度下降了。但是那個在前文提到,交替一個時間步的預(yù)測效果是最好的,因為很顯然,只預(yù)測一步當(dāng)然比預(yù)測兩步的效果好。)
(2.在文章中區(qū)分出lag和多步預(yù)測(預(yù)測步數(shù)為a)的概念。這個概念在之前有些混淆。lag指的是利用當(dāng)前數(shù)據(jù)去預(yù)測下一個數(shù)據(jù)時中間隔了多少秒。預(yù)測步數(shù)a指的是利用當(dāng)前數(shù)據(jù)去預(yù)測下一幾個時間步的數(shù)據(jù)。)
(網(wǎng)上資料繁雜眾多,而且到底如何預(yù)測也沒有一個準(zhǔn)確答案,基本都是教如何訓(xùn)練神經(jīng)網(wǎng)絡(luò)以及驗證神經(jīng)網(wǎng)絡(luò)的準(zhǔn)確度的,但是都沒有一個神經(jīng)網(wǎng)絡(luò)的實戰(zhàn),本文將好好記錄一下我這些年的所思所感以及LSTM神經(jīng)網(wǎng)絡(luò)在實驗中的應(yīng)用,所用源代碼+數(shù)據(jù)(非實驗所用,經(jīng)過簡化處理)將會以文件的方式供大家下載,有需要的自己去下載使用。希望能給大家啟發(fā)。)
LSTM時序預(yù)測神經(jīng)網(wǎng)絡(luò)與NAR、NARx等時序預(yù)測神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)有著本質(zhì)上的不同,像其他幾篇博客所述,通過觀察工作區(qū)發(fā)現(xiàn),NAR和NARx的輸入實際上是“延遲向量”。LSTM神經(jīng)網(wǎng)絡(luò)因為“門”這種特殊的概念引入,使其有著更長、更好的時間預(yù)測效果。?本文將根據(jù)LSTM有兩種預(yù)測方式,一種是根據(jù)預(yù)測值去更新網(wǎng)絡(luò)狀態(tài),另外一種是根據(jù)觀測值去更新網(wǎng)絡(luò)狀態(tài)。
目錄
一、根據(jù)預(yù)測值去更新網(wǎng)絡(luò)狀態(tài)(單序列數(shù)據(jù)用)
Step1:加載數(shù)據(jù),并進(jìn)行數(shù)據(jù)預(yù)處理
Step2:交替一個時間步
Step3:定義LSTM網(wǎng)絡(luò)架構(gòu)并訓(xùn)練
Step4:初始化網(wǎng)絡(luò)狀態(tài)
Step5:利用LSTM網(wǎng)絡(luò)進(jìn)行單步預(yù)測
Step6:將預(yù)測結(jié)果與實際數(shù)據(jù)對比
Step7:重置網(wǎng)絡(luò)狀態(tài)
二、根據(jù)觀測值去更新網(wǎng)絡(luò)狀態(tài)(實驗用)
Step1:加載數(shù)據(jù),并進(jìn)行數(shù)據(jù)預(yù)處理
Step2:構(gòu)造交替5秒時間間隔的時間序列作為神經(jīng)網(wǎng)絡(luò)的輸入和輸出。
Step3:神經(jīng)網(wǎng)絡(luò)訓(xùn)練
Step4:神經(jīng)網(wǎng)絡(luò)預(yù)測+精度反映
Step5:神經(jīng)網(wǎng)絡(luò)實際應(yīng)用
一、根據(jù)預(yù)測值去更新網(wǎng)絡(luò)狀態(tài)(單序列數(shù)據(jù)用)
數(shù)據(jù)集、源碼地址:(直接運行就能用)
放在評論區(qū)了,放在文字里面沒法通過審核。
實現(xiàn)的是:已知前面幾年的數(shù)據(jù),往后預(yù)測2023年的100天的數(shù)據(jù)。
Step1:加載數(shù)據(jù),并進(jìn)行數(shù)據(jù)預(yù)處理
clc clear load dataset.mat %加載數(shù)據(jù)(double型,剔除了2023年的的第一個數(shù)據(jù),總共為2191個。時序預(yù)測沒有實際時間,只有事情發(fā)生的順序) data=data(:,2)'; %不轉(zhuǎn)置的話,無法訓(xùn)練lstm網(wǎng)絡(luò),顯示維度不對。 %% 序列的前2000個用于訓(xùn)練,后191個用于驗證神經(jīng)網(wǎng)絡(luò),然后往后預(yù)測200個數(shù)據(jù)。 dataTrain = data(1:2000); %定義訓(xùn)練集 dataTest = data(2001:2191); %該數(shù)據(jù)是用來在最后與預(yù)測值進(jìn)行對比的 %% 數(shù)據(jù)預(yù)處理 mu = mean(dataTrain); %求均值 sig = std(dataTrain); %求均差 dataTrainStandardized = (dataTrain - mu) / sig;?LSTM對于數(shù)據(jù)標(biāo)準(zhǔn)化要求很高。且這里只對訓(xùn)練集進(jìn)行標(biāo)準(zhǔn)化的原因是經(jīng)過神經(jīng)網(wǎng)絡(luò)中的值只有訓(xùn)練集,因此無須對測試集進(jìn)行標(biāo)準(zhǔn)化。
Step2:交替一個時間步
%% 輸入的每個時間步,LSTM網(wǎng)絡(luò)學(xué)習(xí)預(yù)測下一個時間步,這里交錯一個時間步效果最好。 XTrain = dataTrainStandardized(1:end-1); YTrain = dataTrainStandardized(2:end);大多數(shù)LSTM相關(guān)的資料還是使用默認(rèn)的a=1來將時間序列拆分為具有輸入和輸出分量的樣本。我在自己的其他試驗中使用不同的a來觀察(a=1,2,5,10)。這個a的意思是規(guī)定了一次性往后預(yù)測多少步。設(shè)置不同的預(yù)測步數(shù)a的目的是觀察是否能夠改善預(yù)測模型性能。但據(jù)某些驗證表明,增加滯后并不能改善預(yù)測模型的性能。其實也很好理解,單次預(yù)測精度肯定是比預(yù)測兩步的精度要高的。
a=1即輸入[[1],[2],[3],...,[10]]對應(yīng)著輸出[[2],[3],[4],...[11]],a為2的時候即[[1],[2],[3],...,[10]]對應(yīng)著[[2,3],[3,4],[4,5],...,[10,11]。
這里的交替時間步的作用是什么呢?
Step3:定義LSTM網(wǎng)絡(luò)架構(gòu)并訓(xùn)練
%% 一維特征lstm網(wǎng)絡(luò)訓(xùn)練 numFeatures = 1; %特征為一維 numResponses = 1; %輸出也是一維 numHiddenUnits = 200; %創(chuàng)建LSTM回歸網(wǎng)絡(luò),指定LSTM層的隱含單元個數(shù)200。可調(diào)layers = [ ...sequenceInputLayer(numFeatures) %輸入層lstmLayer(numHiddenUnits) % lstm層,如果是構(gòu)建多層的LSTM模型,可以修改。fullyConnectedLayer(numResponses) %為全連接層,是輸出的維數(shù)。regressionLayer]; %其計算回歸問題的半均方誤差模塊 。即說明這不是在進(jìn)行分類問題。%指定訓(xùn)練選項,求解器設(shè)置為adam, 1000輪訓(xùn)練。 %梯度閾值設(shè)置為 1。指定初始學(xué)習(xí)率 0.01,在 125 輪訓(xùn)練后通過乘以因子 0.2 來降低學(xué)習(xí)率。 options = trainingOptions('adam', ...'MaxEpochs',1000, ...'GradientThreshold',1, ...'InitialLearnRate',0.01, ... 'LearnRateSchedule','piecewise', ...%每當(dāng)經(jīng)過一定數(shù)量的時期時,學(xué)習(xí)率就會乘以一個系數(shù)。'LearnRateDropPeriod',400, ... %乘法之間的紀(jì)元數(shù)由“ LearnRateDropPeriod”控制。可調(diào)'LearnRateDropFactor',0.15, ... %乘法因子由參“ LearnRateDropFactor”控制,可調(diào)'Verbose',0, ... %如果將其設(shè)置為true,則有關(guān)訓(xùn)練進(jìn)度的信息將被打印到命令窗口中。默認(rèn)值為true。'Plots','training-progress'); %構(gòu)建曲線圖 將'training-progress'替換為none net = trainNetwork(XTrain,YTrain,layers,options);需要提醒的是:在MATLAB中,其已經(jīng)將LSTM網(wǎng)絡(luò)進(jìn)行封裝,類似工具箱的形式,因此不再涉及底層的“門”概念。如果利用Python去編程,將會涉及到。
這里有三個參數(shù)很重要,第一個是numHiddenUnits,第二是LearnRateDropPeriod,第三個是LearnRateDropFactor,這三個需要根據(jù)你的數(shù)據(jù)集去改。如果RMSE曲線下降太慢,那么就可能是以下原因:(在本代碼中,需要將LearnRateDropPeriod改大。)
Step4:初始化網(wǎng)絡(luò)狀態(tài)
net = predictAndUpdateState(net,XTrain); %將新的XTrain數(shù)據(jù)用在網(wǎng)絡(luò)上進(jìn)行初始化網(wǎng)絡(luò)狀態(tài) [net,YPred] = predictAndUpdateState(net,YTrain(end)); %用訓(xùn)練的最后一步來進(jìn)行預(yù)測第一個預(yù)測值,給定一個初始值。這是用預(yù)測值更新網(wǎng)絡(luò)狀態(tài)特有的。訓(xùn)練的最后一個,即預(yù)測的第一個。
Step5:利用LSTM網(wǎng)絡(luò)進(jìn)行單步預(yù)測
%% 進(jìn)行用于驗證神經(jīng)網(wǎng)絡(luò)的數(shù)據(jù)預(yù)測(用預(yù)測值更新網(wǎng)絡(luò)狀態(tài)) for i = 2:291 %從第二步開始,這里進(jìn)行191次單步預(yù)測(191為用于驗證的預(yù)測值,100為往后預(yù)測的值。一共291個)[net,YPred(:,i)] = predictAndUpdateState(net,YPred(:,i-1),'ExecutionEnvironment','cpu'); %predictAndUpdateState函數(shù)是一次預(yù)測一個值并更新網(wǎng)絡(luò)狀態(tài) end在其他博客提到的NAR和NARx中,遞歸預(yù)測會一直累積誤差,以至于使性能迅速下降,但LSTM對比其他網(wǎng)絡(luò)獨有的是其有這個特殊的函數(shù)predictAndUpdateState能夠在每次預(yù)測時更新網(wǎng)絡(luò)狀態(tài),因此其能夠預(yù)測更長的序列。
Step6:將預(yù)測結(jié)果與實際數(shù)據(jù)對比
%% 驗證神經(jīng)網(wǎng)絡(luò) YPred = sig*YPred + mu; %使用先前計算的參數(shù)對預(yù)測去標(biāo)準(zhǔn)化。 rmse = sqrt(mean((YPred(1:191)-dataTest).^2)) ; %計算均方根誤差 (RMSE)。 subplot(2,1,1) plot(dataTrain(1:end)) %先畫出前面2000個數(shù)據(jù),是訓(xùn)練數(shù)據(jù)。 hold on idx = 2001:(2000+191); %為橫坐標(biāo) plot(idx,YPred(1:191),'.-') %顯示預(yù)測值 hold off xlabel("Time") ylabel("Case") title("Forecast") legend(["Observed" "Forecast"]) subplot(2,1,2) plot(data) xlabel("Time") ylabel("Case") title("Dataset") %% 繼續(xù)往后預(yù)測2023年的數(shù)據(jù) figure(2) idx = 2001:(2000+291); %為橫坐標(biāo) plot(idx,YPred(1:291),'.-') %顯示預(yù)測值 hold offStep7:重置網(wǎng)絡(luò)狀態(tài)
net = resetState(net);如果需要進(jìn)行其他的預(yù)測,如進(jìn)行下一輪300個的預(yù)測。需要重置一下網(wǎng)絡(luò)狀態(tài)。
訓(xùn)練結(jié)果:
?
左圖表示的是,從2191個開始,往后預(yù)測了100個。右圖表示的是利用2001:2191個數(shù)據(jù)來驗證神經(jīng)網(wǎng)絡(luò)的訓(xùn)練結(jié)果,這里還是比較好的。
但在預(yù)測的時候,不能單從訓(xùn)練集和預(yù)測集的數(shù)量以及比例去評價這個數(shù)據(jù)集的好壞,從理論上來說,訓(xùn)練集應(yīng)包含多個周期的變化,不然肯定無法通過此去進(jìn)行預(yù)測。而且100個數(shù)據(jù)去預(yù)測30個數(shù)據(jù)這并不意味著用1000個數(shù)據(jù)去預(yù)測后面300個是一樣的表現(xiàn),兩者是不相同的。
使預(yù)測模型更加精確的方法有:增加預(yù)測輸入的維度,即加入不同的特征(如下文所提到的多個波高儀)。同時需要注意的是,該模型僅適用于單序列預(yù)測(即僅有波浪力這個數(shù)據(jù),去預(yù)測后面的波浪力,因此其實它的預(yù)測精度也是不高的。)
二、根據(jù)觀測值去更新網(wǎng)絡(luò)狀態(tài)(實驗用、多序列)
在實驗中的需求是:利用5個波高儀的波高數(shù)據(jù)去預(yù)測去一段時間間隔(5秒)后的波浪力。因此神經(jīng)網(wǎng)絡(luò)訓(xùn)練所需要的是足夠數(shù)量的波高、波浪力數(shù)據(jù)樣本,每一步將5個波高儀的數(shù)據(jù)去跟5秒后的波浪力數(shù)據(jù)進(jìn)行對應(yīng)即可。由于實驗數(shù)據(jù)不能泄露,因此大家可參透代碼根據(jù)自己的實際情況進(jìn)行使用。
Step1:加載數(shù)據(jù),并進(jìn)行數(shù)據(jù)預(yù)處理
%% 將所有測得的數(shù)據(jù)切分為7500個數(shù)據(jù)和對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化,消除變量之間的量綱關(guān)系,使數(shù)據(jù)具有可比性。 load testdata.mat %其中第一列為時間,第二列至第六列為波高數(shù)據(jù),第七列為波浪力數(shù)據(jù)。 dt=0.01 %這里的每一步的間隔為0.01秒,可根據(jù)實際實驗采樣頻率情況來定。 time_lag=500; %這里指的是根據(jù)波高力數(shù)據(jù)預(yù)測5秒后的波浪力,這個間距可以自己根據(jù)實際來調(diào)整。 N_total=7500; %這里指的是7500個時間步 t=testdata(:,1); %時間,共有7500個。 std_scale=std(testdata(:,2:7)); %對波高數(shù)據(jù)和波浪力數(shù)據(jù)求標(biāo)準(zhǔn)差 testdata_norm(:,2:7) =normalize(testdata(:,2:7)); %將數(shù)據(jù)歸一化 lag=500;Step2:構(gòu)造交替5秒時間間隔的時間序列作為神經(jīng)網(wǎng)絡(luò)的輸入和輸出。
%% 首先先將時間步分割好 t_input=t(1:end-lag); %7000=7500-500。1:7000個 t_output=t(lag+1:end); %501:7500個%% 分割波高數(shù)據(jù)和波浪力數(shù)據(jù) height_input(1:5,:)=testdata_norm(1:end-lag,2:6)'; %因此這里height_input中的每一行是一個波高儀的數(shù)據(jù),一共有5個波高儀,每個波高儀收集了1-7500個時間步的波高信息。然后這里波高作為輸入,是只截取第1-7000個時間步的數(shù)據(jù)作為輸入。 height_output(1,:)=testdata_norm(lag+1:end,7)'; %這里指的是第501~7500個時間步的波浪力作為輸出。%% 決定網(wǎng)絡(luò)的輸入、輸出、總數(shù) net_input=height_input; %1:7000,指代的是波高信息 net_output=height_output; %501:7500,指代的是波浪力信息 sample_size=length(height_output); %樣本總數(shù)7000Step3:神經(jīng)網(wǎng)絡(luò)訓(xùn)練
%% 訓(xùn)練神經(jīng)網(wǎng)絡(luò)參數(shù)設(shè)定 numHiddenUnits =5; %指定LSTM層的隱含單元個數(shù)為5 train_ratio= 0.8; %劃分用于神經(jīng)網(wǎng)絡(luò)訓(xùn)練的數(shù)據(jù)比例為總數(shù)的80% LearnRateDropPeriod=50; %乘法之間的紀(jì)元數(shù)由“ LearnRateDropPeriod”控制 LearnRateDropFactor=0.5; %乘法因子由參“ LearnRateDropFactor”控制,%% 定義訓(xùn)練時的時間步。 numTimeStepsTrain = floor(train_ratio*numel(net_input(1,:))); %一共為7000個*0.8=5600個%% 交替一個時間步,可以交替多個時間步,但這里交替一個時間步的效果其實是最好的,詳見那篇開頭第二篇建立在python之上的文章。(真正想要改變往后預(yù)測時間的長短為lag這個變量。) XTrain = net_input(:, 1: numTimeStepsTrain+1); %1~5601,XTrain---input,一共為5601個 YTrain = net_output(:, 2: numTimeStepsTrain+2); %2~5602,YTrain---expected output,為5601個%% 輸入有5個特征,輸出有1個特征。 numFeatures = numel(net_input(:,1)); %5 numResponses = numel(net_output(:,1)); %1 layers = [ ...sequenceInputLayer(numFeatures) %輸入層為5lstmLayer(numHiddenUnits) %lstm層,構(gòu)建5層的LSTM模型,fullyConnectedLayer(numResponses) %為全連接層,是輸出的維數(shù)。regressionLayer]; %其計算回歸問題的半均方誤差模塊 。即說明這不是在進(jìn)行分類問題。options = trainingOptions('adam', ... %指定訓(xùn)練選項,求解器設(shè)置為adam, 1000輪訓(xùn)練。'MaxEpochs',150, ... %最大訓(xùn)練周期為150'GradientThreshold',1, ... %梯度閾值設(shè)置為 1'InitialLearnRate',0.01, ... %指定初始學(xué)習(xí)率 0.01'LearnRateSchedule','piecewise', ... %每當(dāng)經(jīng)過一定數(shù)量的時期時,學(xué)習(xí)率就會乘以一個系數(shù)。'LearnRateDropPeriod', LearnRateDropPeriod, ... 'LearnRateDropFactor',LearnRateDropFactor, ... %在50輪訓(xùn)練后通過乘以因子 0.5 來降低學(xué)習(xí)率。'Verbose',0, ... %如果將其設(shè)置為true,則有關(guān)訓(xùn)練進(jìn)度的信息將被打印到命令窗口中,0即是不打印 。'Plots','training-progress'); %構(gòu)建曲線圖 ,不想構(gòu)造就將'training-progress'替換為nonenet = trainNetwork(XTrain,YTrain,layers,options); %訓(xùn)練神經(jīng)網(wǎng)絡(luò) save('LSTM_net', 'net'); %將net保存為LSTM_netendStep4:神經(jīng)網(wǎng)絡(luò)預(yù)測+精度反映
%% 定義Xtest和Ytest,其目的是仿照神經(jīng)網(wǎng)絡(luò)訓(xùn)練中XTrain和YTrain的樣式去構(gòu)造一個測試集 %% input_Test指的是實驗中得到的波高數(shù)據(jù)。 %% output_Test指的就是參與跟預(yù)測值進(jìn)行對比的測量值。即expected output numTimeStepsTrain2 = floor(0.1*sample_size); %一共為7000個*0.1=700個。可以自己再隨意取一個比例,去嘗試一下。當(dāng)然也可以跟上面的numTimeStepsTrain保持一致。 input_Test = net_input(:, numTimeStepsTrain2+1: end-1); %7000個中取701個~6999,一共是6299個 output_Test = net_output(:, numTimeStepsTrain2+2: end); %7000個取702:7000,為6299個,因為 numTimeStepsTrain是700.%% 這里首先用input_Train來初始化神經(jīng)網(wǎng)絡(luò)。經(jīng)過測試,不是將整個input_Train放進(jìn)去最好,放其中一點比例即可。 input_Train = net_input(:, floor(numTimeStepsTrain*0.9): numTimeStepsTrain); %630~700 net = predictAndUpdateState(net, input_Train); %% 預(yù)測量預(yù)定義 rec_step=numel(output_Test); %滾動預(yù)測6299個。跟后面的j循環(huán)有關(guān)。 YPred=zeros(rec_step,1); %6299個預(yù)測,1個特征。這個預(yù)測是和Test來比對,看是否正確的。%% 神經(jīng)網(wǎng)絡(luò)預(yù)測(這個也是之后實際預(yù)測需要用到的)for j=1: rec_step[net,YPred0] = predictAndUpdateState(net, input_Test(:, j)); %用input_Test代入剛剛用input_Train來更新的網(wǎng)絡(luò)得到第一個輸出并得到對應(yīng)的預(yù)測值。YPred(j, :) = YPred0; %記錄輸出的預(yù)測值。end%% 神經(jīng)網(wǎng)絡(luò)精度反映。這里的Y都是波浪力 YTest = output_Test(:,1:rec_step)'; %這一步可以簡化。其實只是一個轉(zhuǎn)置。 NRMSE=sqrt(sum((YPred-YTest).^2)/rec_step/(max(YTest)-min(YTest)));這里需要重點說明的是用input_Train來初始化神經(jīng)網(wǎng)絡(luò)的問題。 因為在LSTM中,state cell是實時在更新的,但對于整個網(wǎng)絡(luò)需要給它一個初始化的值,但不一定是利用整個input_Train的效果最好,因此在實際使用中input_Train可以考慮給予一個比例,來讓其初始化。
Step5:神經(jīng)網(wǎng)絡(luò)實際應(yīng)用
rec_step=2; %這里只預(yù)測兩步 input_Test=[1;2;3;4;5,6;4;8;9;7]; %假設(shè)第一個時間步波高數(shù)據(jù)分別為1,2,3,4,5,第二個時間步波高數(shù)據(jù)分別為6;4;8;9;7。%% 神經(jīng)網(wǎng)絡(luò)預(yù)測for j=1: rec_step[net,YPred0] = predictAndUpdateState(net, input_Test(:, j)); %用input_Test代入剛剛用input_Train來更新的網(wǎng)絡(luò)得到第一個輸出并得到對應(yīng)的預(yù)測值。YPred(j, :) = YPred0; %記錄輸出2個的預(yù)測值。end假設(shè)已經(jīng)同時得到了5個波高儀數(shù)據(jù),需要預(yù)測5秒(在之前的神經(jīng)網(wǎng)絡(luò)訓(xùn)練中已經(jīng)預(yù)設(shè)好了)后的波浪力。每次得到的波高數(shù)據(jù)應(yīng)該是放在for循環(huán)里面,去不斷的進(jìn)行單步預(yù)測。
?左邊是數(shù)據(jù)集,右邊是預(yù)測,我在預(yù)測中的輸入是2秒前的波高,預(yù)測了2秒后的波浪力。這個相當(dāng)于多次的單步預(yù)測,還是相當(dāng)?shù)臏?zhǔn)確的。
總結(jié)
以上是生活随笔為你收集整理的基于Matlab的LSTM神经网络时序预测(完整代码+范例数据文件)的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: log4j自定义配置文件(SpringM
- 下一篇: 郸城二高2021年高考成绩查询时间,河南