%建立符號變數a(發展係數)和b(灰作用量)
syms a b;
c =[a b]';%原始數列 AA= x0';
n =length(A);%對原始數列 A 做累加得到數列 BB=cumsum(A);%對數列 B 做緊鄰均值生成
for i =2:n
C(i)=(B(i)+B(i -1))/2;
end
C(1)=[];%構造資料矩陣
B=[-C;ones(1,n-1)];Y=A;Y(1)=[];Y=Y';%使用最小二乘法計算引數 a(發展係數)和b(灰作用量)
c =inv(B*B')*B*Y;
c = c';
a =c(1); b =c(2);%預測後續資料
F=[];F(1)=A(1);for i =2:(n+10)F(i)=(A(1)-b/a)/exp(a*(i-1))+ b/a;
end
%對數列 F 累減還原,得到預測出的資料
G=[];G(1)=A(1);for i =2:(n+10)G(i)=F(i)-F(i-1);%得到預測出來的資料
end
disp('預測資料為:');%模型檢驗
H=G(1:243);%計算殘差序列
epsilon =A-H;%法一:相對殘差Q檢驗
%計算相對誤差序列
delta =abs(epsilon./A);%計算相對誤差Qdisp('相對殘差Q檢驗:')Q=mean(delta)%法二:方差比C檢驗
disp('方差比C檢驗:')C=std(epsilon,1)/std(A,1)%法三:小誤差概率P檢驗
S1=std(A,1);
tmp =find(abs(epsilon -mean(epsilon))<0.6745*S1);disp('小誤差概率P檢驗:')P=length(tmp)/n
%繪製曲線圖
t1 =1:243;
t2 =1:253;plot(t1,A,'ro'); hold on;plot(t2,G,'g-');xlabel('');ylabel('');legend('實際','預測');title('');
grid on;