clc
close all
clear all
%训练数据
XTrain = xlsread('GFD1.xlsx',1,'B2:G2976');
YTrain = xlsread('GFD1.xlsx',1,'H2:H2976');
mu = mean(XTrain,'ALL');
sig = std(XTrain,0,'ALL');
XTrain = (XTrain - mu)/sig;
YTrain = (YTrain - mu)/sig;
XTrain = XTrain';
YTrain = YTrain';
XTest = xlsread('GFD1.xlsx',1,'B2977:G3072');
YTest = xlsread('GFD1.xlsx',1,'H2977:H3072');
XTest = (XTest - mu)/sig;
YTest = (YTest - mu)/sig;
XTest = XTest';
YTest = YTest';
%% define the Deeper LSTM networks
%创建LSTM回归网络,指定LSTM层的隐含单元个数125
%序列预测,因此,输入一维,输出一维
numFeatures= 6;%输入特征的维度
numResponses = 1;%响应特征的维度
numHiddenUnits = 100;%隐藏单元个数
layers = [sequenceInputLayer(numFeatures)
lstmLayer(numHiddenUnits)
dropoutLayer(0.5)%防止过拟合
fullyConnectedLayer(numResponses)
regressionLayer];
MaxEpochs=300;%最大迭代次数
InitialLearnRate=0.005;%初始学习率
%% back up to LSTM model
options = trainingOptions('adam', ...
'MaxEpochs',MaxEpochs,...%30用于训练的最大轮次数,迭代是梯度下降算法中使用小批处理来最小化损失函数的一个步骤。一个轮次是训练算法在整个训练集上的一次全部遍历。
'MiniBatchSize',128, ... %128小批量大小,用于每个训练迭代的小批处理的大小,小批处理是用于评估损失函数梯度和更新权重的训练集的子集。
'GradientThreshold',1, ...%梯度下降阈值
'InitialLearnRate',InitialLearnRate, ...%全局学习率。默认率是0.01,但是如果网络训练不收敛,你可能希望选择一个更小的值。默认情况下,trainNetwork在整个训练过程中都会使用这个值,除非选择在每个特定的时间段内通过乘以一个因子来更改这个值。而不是使用一个小的固定学习速率在整个训练, 在训练开始时选择较大的学习率,并在优化过程中逐步降低学习率,可以帮助缩短训练时间,同时随着训练的进行,可以使更小的步骤达到最优值,从而在训练结束时进行更精细的搜索。
'LearnRateSchedule','piecewise', ...%在训练期间降低学习率的选项,默认学习率不变。'piecewise’表示在学习过程中降低学习率
'LearnRateDropPeriod',100, ...% 10降低学习率的轮次数量,每次通过指定数量的epoch时,将全局学习率与学习率降低因子相乘
'LearnRateDropFactor',0.5, ...%学习下降因子
'ValidationData',{XTrain,YTrain}, ...% ValidationData用于验证网络性能的数据,即验证集
'ValidationFrequency',1, ...%以迭代次数表示的网络验证频率,“ValidationFrequency”值是验证度量值之间的迭代次数。
'Verbose',1, ...%1指示符,用于在命令窗口中显示训练进度信息,显示的信息包括轮次数、迭代数、时间、小批量的损失、小批量的准确性和基本学习率。训练一个回归网络时,显示的是均方根误差(RMSE),而不是精度。如果在训练期间验证网络。
'Plots','training-progress');
%% Train LSTM Network
net = trainNetwork(XTrain,YTrain,layers,options);
% net = predictAndUpdateState(net,XTrain);
numTimeStepsTrain = numel(XTrain(1,:));
for i = 1:numTimeStepsTrain
[net,YPred_Train(i)] = predictAndUpdateState(net,XTrain(:,i),'ExecutionEnvironment','cpu');
end
YPred_Train = sig*YPred_Train + mu;
YTrain = sig*YTrain + mu;
numTimeStepsTest = numel(XTest(1,:));
for i = 1:numTimeStepsTest
[net,YPred(i)] = predictAndUpdateState(net,XTest(:,i),'ExecutionEnvironment','cpu');
end
YPred = sig*YPred + mu;
YTest = sig*YTest + mu;
% % 评价
ae= abs(YPred - YTest);
RMSE = (mean(ae.^2)).^0.5;
MSE = mean(ae.^2);
MAE = mean(ae);
MAPE = mean(ae./YPred);
disp('预测结果评价指标:')
disp(['RMSE = ', num2str(RMSE)])
disp(['MSE = ', num2str(MSE)])
disp(['MAE = ', num2str(MAE)])
disp(['MAPE = ', num2str(MAPE)])
%
% MSE = double(mse(YTest,YPred))
% RMSE = sqrt(MSE)
%%
figure
hold on
plot(1:length(YPred_Train),YPred_Train,'r-','LineWidth',1.5);
plot(1:length(YTrain),YTrain,'b-','LineWidth',1.5);
legend('预测值','实际值')
xlabel('训练样本序号')
%%
confidence_interval = [0.05 0.1 0.2 0.3 0.4]; % 置信区间的大小,可以根据需要调整
% 创建上下置信区间的数据
for i=1:1:length(confidence_interval)
y_upper(:,i) = YPred + confidence_interval(i)*YPred;
y_lower(:,i) = YPred - confidence_interval(i)*YPred;
end
figure
plot(1:length(YPred),YPred,'b-','LineWidth',1.5);
hold on
for i=1:1:length(confidence_interval)
% 使用 fill 函数绘制置信区间
fill([1:length(YPred), fliplr(1:length(YPred))], [y_upper(:,i)', fliplr(y_lower(:,i)')], 'c', 'FaceAlpha', 0.1);
end
plot(1:length(YTest),YTest,'r--','LineWidth',1.5);
legend('预测值','95%置信区间','90%置信区间','80%置信区间','70%置信区间','60%置信区间','真实值')
xlabel('测试样本序号')
%%
% 假设已有一组数据
data = [XTest(:,10:18);YPred(10:18)]
% 计算均值和方差
mu = mean(data);
sigma = std(data);
rows = 3;
cols = 3;
% 生成正态分布概率密度曲线的 x 值
for i=1:1:rows*cols
x = linspace(min(data(:,i))-2, max(data(:,i)), 1000);
% 计算概率密度函数值
pdf_values = normpdf(x, mu(i), sigma(i));
% 绘制概率密度曲线
% 设置子图的行数和列数
% 创建新的图形窗口
% 使用 for 循环绘制子图
subplot(rows, cols, i);
% 在每个子图中绘制相应的数据
switch i
case 1
area(x, pdf_values, 'FaceColor',[0.9, 0.9, 0.9], 'EdgeColor', [0.3010 0.1450 0.9330]);
% 添加竖直的红线
x_position = data(9+i); % 你希望添加红线的位置
line([x_position, x_position], ylim, 'Color', 'red', 'LineStyle', '-', 'LineWidth', 1);
legend('概率密度','观测值')
% 设置图标题和轴标签
title('正态分布概率密度曲线');
xlabel('X轴');
ylabel('概率密度');
case 2
area(x, pdf_values, 'FaceColor',[0.9, 0.9, 0.9], 'EdgeColor', [0.3010 0.1450 0.9330]);
% 添加竖直的红线
x_position = data(9+i); % 你希望添加红线的位置
line([x_position, x_position], ylim, 'Color', 'red', 'LineStyle', '-', 'LineWidth', 1);
legend('概率密度','观测值')
% 设置图标题和轴标签
title('正态分布概率密度曲线');
xlabel('X轴');
ylabel('概率密度');
case 3
area(x, pdf_values, 'FaceColor',[0.9, 0.9, 0.9], 'EdgeColor', [0.3010 0.1450 0.9330]);
% 添加竖直的红线
x_position = data(9+i); % 你希望添加红线的位置
line([x_position, x_position], ylim, 'Color', 'red', 'LineStyle', '-', 'LineWidth', 1);
legend('概率密度','观测值')
% 设置图标题和轴标签
title('正态分布概率密度曲线');
xlabel('X轴');
ylabel('概率密度');
case 4
area(x, pdf_values, 'FaceColor',[0.9, 0.9, 0.9], 'EdgeColor', [0.3010 0.1450 0.9330]);
% 添加竖直的红线
x_position = data(9+i); % 你希望添加红线的位置
line([x_position, x_position], ylim, 'Color', 'red', 'LineStyle', '-', 'LineWidth', 1);
legend('概率密度','观测值')
% 设置图标题和轴标签
title('正态分布概率密度曲线');
xlabel('X轴');
ylabel('概率密度');
case 5
area(x, pdf_values, 'FaceColor',[0.9, 0.9, 0.9], 'EdgeColor', [0.3010 0.1450 0.9330]);
% 添加竖直