灰色系统模型G(1,1)&G(1,N)

G(1,1):

    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）

    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