灰色预测
clc,clear;
close all;
syms a u;
c=[a,u]';%构成矩阵
A=[2006 1138 993 1027 1155 1175 1299 1786 1781 1570 1387 1635 2795 1679 1329 1585 1406 1260 1814 1556 2143 1963 1555 1217 1738 1261 1141 1099 1274 1325 1197 1596 1923 1789 1606 2174 3089 1746 1571 1545 1593 1480 1925 1290 2637 1907 1922 1660 2355 1104 1093 1085 1125 1257 1055 1547 1927 1881 1490 2156 3147 1706 1599 1544 1738 1747 2380 2049 2713 2543 2139 2105 2378 1461 1295 1313 1360 1582 1425 1656 1822 1602 1404 1837 2589 1052 956 895 972 1124 1620 1280 1492 1587 1330 1219 2064 1347 1325 1273 1461 1481 1511 2000 2137 1643 1464 1706 2236 1065 994 907 948 997 1788 1394 1787 1435 1225 1221 1900 1259 1357 509 557 1375 1548 1802 1995 1720 1468 1625 2515 1371 1439 1642 1344 1247 1737 1434 1782 1595 1348 1083 1864 1367 1299 1301 1308 1382 1401 1837 2045 1790 1409 1179 1852 796 780 602 584 558 831 639 943 960 859 775 976 581 563 547 635 664 684 828 981 889 1303 1490 2246 1335 1270 1281 1210 1078 1664 1337 1632 1454 1333 1146 1721 1059 963 1000 1107 1132 1086 1399 1575 1507 1469 1650 2305 1328 1323 1314 1317 1443 1963 2201 2577 1900 1606 1390 1859 1009 1001 940 981 1239 1023 1284 1510 1497 1448 1677 2188 1499 1313 1351 1237 1454 2031 1615 2230 1899 1574 1666];
%输入数据,可以修改
Ago=cumsum(A);%原始数据一次累加,得到1-AGO序列xi(1)。
n=length(A);%原始数据个数
for k=1:(n-1)
Z(k)=(Ago(k)+Ago(k+1))/2; %Z(i)为xi(1)的紧邻均值生成序列
end
Yn =A;%Yn为常数项向量
Yn(1)=[]; %从第二个数开始,即x(2),x(3)...
Yn=Yn';
E=[-Z;ones(1,n-1)]';%累加生成数据做均值
c=(E'*E)\(E'*Yn);%利用公式求出a,u
c= c';
a=c(1);%得到a的值
u=c(2);%得到u的值
F=[];
F(1)=A(1);
for k=2:(n)
F(k)=(A(1)-u/a)/exp(a*(k-1))+u/a;%求出GM(1,1)模型公式
end
G=[];
G(1)=A(1);
for k=2:(n)
G(k)=F(k)-F(k-1);%两者作差还原原序列,得到预测数据
end
t1=1:n;
t2=1:n;
plot(t1,A,'bo--');
hold on;
plot(t2,G,'r*-');
title('预测结果');
legend('真实值','预测值');
%后验差检验
e=A-G;
q=e/A;%相对误差
s1=var(A);
s2=var(e);
c=s2/s1;%方差比
len=length(e);
p=0; %小误差概率
for i=1:len
if(abs(e(i))<0.6745*s1)
p=p+1;
end
end
p=p/len;