matlab GM(1,1)画出曲线图

function y=gm()
%??????a(????)?b(????)
syms a b;
c = [a b]';

%???? A
A = [41462, 28969, 32251, 47694, 32265, 27819, 33539, 36994];%VA
% A2=[70999	71282	85415	93747	101423	109150	115276	119349];%OH
% A3=[89981	86793	78577	72096	77318	75351	72376	68751];%PA
% A4=[29588	28285	27502	26820	27077	25811	26530	28870];%KY
% A5=[8668	9310	9429	9062	6926	5345	5405	3672];%WV
n=length(A)
% n1 = length(A1);
% n2 = length(A2);
% n3 = length(A3);
% n4 = length(A4);
% n5 = length(A5);
B=cumsum(A)
%????? A ??????? B
% B = cumsum(A1);
% B = cumsum(A2);
% B = cumsum(A3);
% B = cumsum(A4);
% B = cumsum(A5);

%??? 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';

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:8);
%??????
epsilon = A - H
%???????Q??
%????????
delta = abs(epsilon./A)
%??????Q
disp('????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 = 2010:2017;
t2 = 2010:2027;

plot(t1, A,'ro'); hold on;
plot(t2, G, 'g-');

xlabel('Year');
ylabel('TDRS');
legend('VA','Predict');
title('????');
grid on;