G(1,1):

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
clc;clear;
%建立符号变量a(发展系数)和b(灰作用量)
syms a b;
c = [a b]';


%原始数列 A
A = [174, 179, 183, 189, 207, 234, 220.5, 256, 270, 285];
n = length(A);

%对原始数列 A 做累加得到数列 B
B = 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('预测数据为:');
G

%模型检验

H = G(1:10);
%计算残差序列
epsilon = A - H;

%法一:相对残差Q检验
%计算相对误差序列
delta = abs(epsilon./A);
%计算相对误差Q
disp('相对残差Q检验:')
Q = mean(delta)

%法二:方差比C检验
disp('方差比C检验:')
C = std(epsilon, 1)/std(A, 1)
q
%法三:小误差概率P检验
S1 = std(A, 1);
tmp = find(abs(epsilon - mean(epsilon))< 0.6745 * S1);
disp('小误差概率P检验:')
P = length(tmp)/n

%绘制曲线图
t1 = 1995:2004;
t2 = 1995:2014;

plot(t1, A,'ro'); hold on;
plot(t2, G, 'g-');
xlabel('x'); ylabel('y');
legend('实际','预测');
title('量增长曲线');
grid on;

G(1,n)

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
A=[560823,542386,604834,591248,583031,640636,575688,689637,570790,519574,614677];
x0=[104,101.8,105.8,111.5,115.97,120.03,113.3,116.4,105.1,83.4,73.3;
135.6,140.2,140.1,146.9,144,143,133.3,135.7,125.8,98.5,99.8;
131.6,135.5,142.6,143.2,142.2,138.4,138.4,135,122.5,87.2,96.5;
54.2,54.9,54.8,56.3,54.5,54.6,54.9,54.8,49.3,41.5,48.9];
[n,m]=size(x0);
AGO=cumsum(A);
T=1;
x1=zeros(n,m+T);

for k=1:(m-1)
Z(k)=(AGO(k)+AGO(k+1))/2; %Z(i)为xi(1)的紧邻均值生成序列
end
for i=1:n
for j=1:m
for k=1:j
x1(i,j)=x1(i,j)+x0(i,k);%原始数据一次累加,得到xi(1)
end
end
end
x11=x1(:,1:m);
X=x1(:,2:m)';%截取矩阵
Yn =A;%Yn为常数项向量
Yn(1)=[]; %从第二个数开始,即x(2),x(3)...
Yn=Yn';
%Yn=A(:,2:m)';
B=[-Z',X];
C=((B'*B)\(B'*Yn))';%由公式建立GM(1,n)模型
a=C(1);
b=C(:,2:n+1);
F=[];
F(1)=A(1);
u=zeros(1,m);
for i=1:m
for j=1:n
u(i)=u(i)+(b(j)*x11(j,i));
end
end
for k=2:m
F(k)=(A(1)-u(k-1)/a)/exp(a*(k-1))+u(k-1)/a;
end
G=[];
G(1)=A(1);
for k=2:m
G(k)=F(k)-F(k-1);%两者做差还原原序列,得到预测数据
end
t1=1:m;
t2=1:m;
plot(t1,A,'bo--');
hold on;
plot(t2,G,'r*-');
title('销售预测结果');
legend('真实值','预测值');

https://blog.csdn.net/wuli_dear_wang/article/details/80587650