-
Notifications
You must be signed in to change notification settings - Fork 0
/
HepatitisB-Region
100 lines (98 loc) · 3.28 KB
/
HepatitisB-Region
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
clc;clear;
data1=xlsread('data6.xlsx');
data1=data1(1:84,:);
num=72;
data=data1(1:num,:);
test1=data1(num+1:end,:);
y = log(data);
T = length(y);
Mdl = arima('Constant',0,'D',1,'Seasonality',12,...
'MALags',1,'SMALags',12);
EstMdl = estimate(Mdl,y);
res = infer(EstMdl,y);
stres = res/sqrt(EstMdl.Variance);
figure(1);
qqplot(stres)
subplot(2,1,1)
autocorr(stres)
subplot(2,1,2)
parcorr(stres)
[h,p] = lbqtest(stres,'lags',[5,10,15],'dof',[3,8,13]);
yfit=y-res;
yfit2008=exp(yfit);
y2008=exp(y);
resinput=y2008-yfit2008;
resinput=resinput(14:end);
a=13;
xtrianres=resinput(1:end-a);
ytrainres=resinput(a+1:end);
mu = mean(resinput);
sig = std(resinput);
XTrain = (xtrianres - mu) / sig;
XTrain=XTrain';
YTrain = (ytrainres - mu) / sig;
YTrain=YTrain';
numFeatures = 1;
numResponses = 1;
numHiddenUnits = 60*3;
layers = [ ...
sequenceInputLayer(numFeatures)
lstmLayer(numHiddenUnits)
fullyConnectedLayer(numResponses)
regressionLayer];
options = trainingOptions('adam', ...
'MaxEpochs',250, ...
'GradientThreshold',1, ...
'InitialLearnRate',0.005, ...
'LearnRateSchedule','piecewise', ...
'LearnRateDropPeriod',125, ...
'LearnRateDropFactor',0.2, ...
'Verbose',0, ...
'Plots','training-progress');
net = trainNetwork(XTrain,YTrain,layers,options);
dataTrainStandardized = (resinput - mu) / sig;
dataTrainStandardized=dataTrainStandardized';
Ytrainres=predict(net,dataTrainStandardized);
Ytrainres=Ytrainres';
trainresoutput=Ytrainres(1:end-a);
testresoutput=Ytrainres(end-a+1:end-a+12);
Mdl1 = estimate(Mdl,y);
yF1 = forecast(Mdl1,12,'Y0',y);
res2=test1-yF1;
yF11=exp(yF1);
yfit2008_2018=yfit2008(14+a:end)+trainresoutput;
y2008_2018=y2008(14+a:end);
yF1=yF11+testresoutput;
figure(10);
plot(yfit2008_2018,'r--','LineWidth',1.5);
hold on;
plot(y2008_2018,'b','LineWidth',2);
title('2008-2018年LSTM_SARIMA组合模型拟合丙肝发病率曲线')
legend('拟合值','真实值','Location','NorthWest')
num2=length(y2008_2018);
mse2 = sqrt(sum((y2008_2018-yfit2008_2018).^2)) ./ num2;
mae2 = mean(abs(y2008_2018-yfit2008_2018));
rmse2= sqrt(mean((yfit2008_2018-y2008_2018).^2));
R2 = 1 - (sum((yfit2008_2018-y2008_2018).^2) / sum((y2008_2018- mean(y2008_2018)).^2));
disp(['2013-2019年LSTM_sarima组合模型预测发病率均方误差MSE为:',num2str(mse2)])
disp(['2013-2019年预测发病率平均绝对误差MAE为:',num2str(mae2)])
disp(['2013-2019年预测发病率均方根误差RMSE为:',num2str(rmse2)])
disp(['2013-2019年预测发病率决定系数为:',num2str(R2)])
figure(11);
plot(yF1,'r--','LineWidth',1.5)
hold on
plot(test1,'b','LineWidth',2)
axis tight
title('2019年LSTM_sarima组合模型预测丙肝发病率曲线')
legend('预测值','真实值','Location','NorthWest')
hold off
num1=12;
mse1 = sqrt(sum((test1-yF1).^2)) ./ num1;
mae1 = mean(abs(test1- yF1));
rmse1 = sqrt(mean((yF1-test1).^2));
R1 = 1 - (sum((yF1-test1).^2) / sum((test1- mean(test1)).^2));
disp(['2020年LSTM_sarima组合模型预测发病率均方误差MSE为:',num2str(mse1)])
disp(['2020年LSTM_sarima组合模型预测发病率平均绝对误差MAE为:',num2str(mae1)])
disp(['2020年LSTM_sarima组合模型预测发病率均方根误差RMSE为:',num2str(rmse1)])
disp(['2020年LSTM_sarima组合模型预测发病率决定系数为:',num2str(R1)])
%save filename14.mat;