预测模型的建立（一）

一.灰色预测GM11模型的原理

        在本讲我们首先将介绍灰色预测模型背后的原理，同时也将在原始的灰色预测模型的基础上作出改进与创新，然后将简要介绍神经网络在数据预测中的应用，在这里，我们只简单的介绍神经网络的思想，以及如何操作matlab中的神经网络工具箱以对数据进行预测，在本讲最后，将谈谈个人对于数据预测的一些看法。

        灰色预测模型可以理解为是给我们的灰色系统建模。那么，我们就要区分什么是灰色系统。

        灰色系统是系统的部分信息已知而部分信息未知，灰色预测是对既含有已知信息又含有不确定信息的系统进行预测。灰色预测通过对原始数据进行生成处理来寻找系统变动的规律，并生成有较强规律性的数据序列，然后建立相应的微分方程模型，可见，灰色预测模型的原理其实是个微分方程模型。

        我们使用的灰色模型叫GM(1,1)模型，其中，第一个'1'表示建立的微分方程是一阶的，后面的‘1’表示只有一个变量。这里我们用随机数模拟生成了12个时期对应的12个原始数据。

        这里的X(0)表示我们的原始数据，注意我们的原始数据一定要得是非负数据列，经过一次累加得到新的数据列X(1)，X(1)我们也称为X(0)的1-AGO序列，即累加序列。通过解释我们可看出，X(1)数据列的第m个元素等于X(0)的前m个元素累加之和得到的结果。又定义了一个Z(1)数据列称为数列X(1)的紧邻均值生成数列，通过式子可看出Z(1)数列的第m个元素等于X(1)序列的第m个元素和第m-1个元素的加权平均，其中，权重是0.5。由于涉及X(1)的第m-1个元素，所以，Z(1)数据列是从m=2开始的。通过观察图片，原始数据是完全随机的，但经过一次累加之和的累加数据我们削减了数据的随机性，数据大致上是在一条直线上，而紧邻均值生成序列的直线拟合效果上是更好的。

        这里，我们引入一个新的定义，叫GM(1,1)模型的基本形式，是一个方程形式，其中，b表示灰作用量，-a表示发展系数。其实就是将原始数据序列视为因变量y，Z(1)序列视为自变量x，进行一元线性回归。回归方程中的k就是-a，b就是灰作用量。还介绍了该方程的矩阵表达形式。既然，我们将方程转化成了一元线性回归的表达式，那么我们就可以用最小二乘法得到a，b参数的估计值了。

        在这里，我们也模拟建立了一个线性回归，最后得到的灰作用量等于13.819，发展系数等于0.0203。

        同一元线性回归，这里我们通过最小二乘法来估计发展系数和灰作用量的值。基本思路与一元线性回归的一样，同样，我们令最后估计出的k和b的值等于残差平方和取最小时对应的k和b，令残差平方和为L，矩阵Y等于原始数据的列向量，矩阵X的第一列均为1，第二列为紧邻均值生成序列的列向量，矩阵β为要估计的b和k，通过矩阵的乘法和减法，我们可发现，L就等于(Y-Xβ)矩阵乘(Y-Xβ)矩阵的转置。

        这里，我们求β矩阵的估计值是通过L矩阵对β矩阵求偏导得出来的。这里涉及到了矩阵的求导，可参考下述网址来理解矩阵求导的过程。总之，最后我们得出β矩阵的估计值等于X的转置乘X的逆乘上X的转置再乘上Y。

        在之前，我们曾提到过完全多重共线性的问题。这里最后我们求解出的β矩阵的估计值的表达式要求X的转置乘X一定要可逆，事实上，如果我们不满足这个条件的话，就会引发完全多重共线性问题。具体的解释可参考ppt。

        利用OLS估计我们最后可以得到a和b的估计值。这里，我们再次对GM(1,1)模型的基本形式进行变形。方程左边X(0)的第k项可以写成是X(1)的第k项减去第k-1项，然后我们可以通过牛顿莱布尼兹公式将等式转换乘一个定积分的形式，即将一个离散的数值转换成一个连续型数值，总之，我们将方程左边转换成定积分函数。同时，方程右边的Z(1)也可以转换成定积分函数的形式，依靠它的定义，我们使用定积分的几何意义可以近似相等。总之，最后，我们可将GM(1,1)模型的基本形式转换成微分方程的形式。该微分方程称为GM(1,1)模型的白化方程，而之前的GM(1,1)模型的基本形式则被称为灰色微分方程。

        可能有疑问是为什么我们要经过那么多次变换将灰色微分方程转换成白化方程的形式呢。其实，都是为了便于求解，在高数中学到过微分方程的求解方式，给出一个微分方程，只需给出一个初始值，我们便可以求出对应的解。前面提到过，这里X(1)t的定义范围是从2-n，我们不妨直接令t=m+1，这样，我们m的定义范围就是1-n-1了。

        这里给出了一些微分方程的求解方式。

        通过白化方程以及初始值原始数据的第一项，最终我们可以求出对应的解。即可以求出累加数据的第m+1项，然后我们再将累加数据转换为原始数据，便可得到原始数据的第m+1项达到一个预测的效果。在这里，如果我们要对原始数据进行预测，只需要在上式最后的求解结果中令m大于等于我们原始数据的样本个数即可。

        这里说GM(1,1)模型的本质是有条件的指数拟合，其实，不难理解，因为我们最后求出的对应的解跟指数拟合求出的解十分相似。在这里，我们忽略最后求出解的常数项。

        既然我们得出了GM(1,1)模型的本质是有条件的指数拟合，那么在建模前，我们就需要对数据进行分析。数据具有准指数规律是使用灰色系统建模的理论基础，这里我们称数据检验为准指数规律的检验。

        在这里，我们定义了累加序列的级比和原始数据的光滑比。对于GM(1,1)模型，我们需要判断累加一次后的序列X(1)是否具有准指数规律。ppt中给出了级比的定义式，等于累加序列的第k项与第k-1项的比值，要求我们的级比收敛于一个区间内，并且该区间长度小于0.5，就称累加一次之后的序列具有准指数规律。

        通过对级比定义式的变形，我们可得到原始数据的光滑比的定义，随着我们数据的增加，光滑比最后会逐渐接近0，同样的，要使累加序列具有准指数规律，就需要保证光滑比收敛在0-0.5之内。

        在实际建模中，我们可得到光滑比的向量，需要计算出光滑比在0-0.5之内的占比，占比越高越好，一般光滑比的前两项可能不符合要求，我们重点关注后面光滑比的项数。

        这里，我们得到估计出的发展系数的估计值，具体参考的是一本教材中说GM(1,1)适用情况和发展系数的大小有很大关系，针对不同的发展系数，GM(1,1)适用预测的情况也不一样，有适合于中期和长期数据的预测，也有适合于短期预测的。但这个结论并不严谨，参考下表的预测误差，可发现，发展系数越小的时候预测的越准确，根本体现不出对短期数据和中长期数据预测的适应性如何。

        在我们使用GM(1,1)模型对未来数据进行预测时，我们需要先检验GM(1,1)模型对原数据的拟合程度。一般有两种检验方法，第一种是残差检验。

        我们通过计算平均相对残差的值来判断GM(1,1)对原数据的拟合效果。图中，给了我们一种评价的方式，即认为平均相对残差的值小于10%时，我们认为GM(1,1)模型对原数据的拟合效果非常不错。但这个标准的值得结合实际预测的情景，比如，我们在进行物理领域内一些要求特别高的预测情景时，可能我们就要规定认为平均相对残差的值小于5%甚至更小的时候，GM(1,1)模型对原数据的拟合效果不错。

        第二种检验方式是级比偏差检验。

        最后是对GM(1,1)模型的拓展。通常情况下，我们建立的GM(1,1)模型被称为全数据GM(1,1)，在这种模型下，我们使用了所有的原始数据；我们只使用部分的原始数据建立的GM(1,1)模型称为部分数据GM(1,1)；第三种情况是当我们预测出新的原始数据时，在下一次使用GM(1,1)模型时，我们会将使用预测出的新的原始数据建立GM(1,1)模型，这种建立的模型为新信息GM(1,1)；在第四种情况，我们在置入最新信息时，去掉最老信息建立的模型为新陈代谢GM(1,1)。

        从模拟精度来看，新陈代谢的模型高，从预测角度看，新陈代谢模型是最理想的模型。其实，这也不难理解，随着系统的发展，老数据的信息意义将逐步降低，在不断补充新信息的同时，及时地去掉老信息，会使建模序列更能反映出系统在目前地特征。

        最后，我们再强调一下使用灰色预测的前提条件，这也为我们代码的书写提供了基础。首先，数据一定要是以年份度量的非负数据，这里也指出了，如果数据是月份或者季度数据要使用上一讲学过的时间序列模型。时间序列模型我已经学过了，我也会尽快在csdn中更新出来。

        第二点是数据要经过准指数规律的检验，其次数据的期数较短且和其他数据之间的关联性不强，但也不能太短了，比如只有三期数据的话，我们自己带入计算预测出一个数值就好了，同时，要是数据期数较长，一般使用传统的时间序列模型比较合适。

二.灰色预测代码的讲解

        这里，我们通过一个例题来讲解灰色预测的代码。

        这便是预测题目中的一些固定的模板套路了，我们重点关注今天这节GM(1,1)模型代码的讲解。

        以上便是我们使用GM(1,1)模型解决预测问题时的代码讲解了。以文字的表述描述出每一段代码的含义，代码也是依靠以上步骤来编写的。

%%  输入原始数据并做出时间序列图
clear;clc
year =[1995:1:2004]';  % 横坐标表示年份，写成列向量的形式（加'就表示转置）
x0 = [174,179,183,189,207,234,220.5,256,270,285]';  %原始数据序列，写成列向量的形式（加'就表示转置）
% year = [2009:2015]; % 其实本程序写成了行向量也可以，因为我怕你们真的这么写了，所以在后面会有判断。
% x0 = [730, 679, 632, 599, 589, 532, 511];
% year = [2010:2017]';   % 该数据很特殊，可以通过准指数规律检验，但是预测效果却很差
% x0 = [1.321,0.387,0.651,0.985,1.235,0.987,0.854,1.021]';
% year = [2014:2017]';
% x0 = [2.874,3.278,3.337,3.390]';

% 画出原始数据的时间序列图
figure(1); % 因为我们的图形不止一个，因此要设置编号
plot(year,x0,'o-'); grid on;  % 原式数据的时间序列图
set(gca,'xtick',year(1:1:end))  % 设置x轴横坐标的间隔为1
xlabel('年份');  ylabel('排污总量');  % 给坐标轴加上标签


%% 因为我们要使用GM(1,1)模型，其适用于数据期数较短的非负时间序列
ERROR = 0;  % 建立一个错误指标，一旦出错就指定为1
% 判断是否有负数元素
if sum(x0<0) > 0  % x0<0返回一个逻辑数组(0-1组成)，如果有数据小于0，则所在位置为1，如果原始数据均为非负数，那么这个逻辑数组中全为0，求和后也是0~
    disp('亲，灰色预测的时间序列中不能有负数哦')
    ERROR = 1;
end

% 判断数据量是否太少
n = length(x0);  % 计算原始数据的长度
disp(strcat('原始数据的长度为',num2str(n)))    % strcat()是连接字符串的函数，第一讲学了，可别忘了哦
if n<=3
    disp('亲，数据量太小，我无能为力哦')
    ERROR = 1;
end

% 数据太多时提示可考虑使用其他方法（不报错）
if n>10
    disp('亲，这么多数据量，一定要考虑使用其他的方法哦，例如ARIMA，指数平滑等')
end

% 判断数据是否为列向量，如果输入的是行向量则转置为列向量
if size(x0,1) == 1
    x0 = x0';
end
if size(year,1) == 1
    year = year';
end


%% 对一次累加后的数据进行准指数规律的检验(注意，这个检验有时候即使能通过，也不一定能保证预测结果非常好，例如上面的第三组数据)
if ERROR == 0   % 如果上述错误均没有发生时，才能执行下面的操作步骤
    disp('------------------------------------------------------------')
    disp('准指数规律检验')
    x1 = cumsum(x0);   % 生成1-AGO序列，cumsum是累加函数哦~    注意：1.0e+03 *0.1740的意思是科学计数法,即10^3*0.1740 = 174
    rho = x0(2:end) ./ x1(1:end-1) ;   % 计算光滑度rho(k) = x0(k)/x1(k-1)
    
    % 画出光滑度的图形，并画上0.5的直线，表示临界值
    figure(2)
    plot(year(2:end),rho,'o-',[year(2),year(end)],[0.5,0.5],'-'); grid on;
    text(year(end-1)+0.2,0.55,'临界线')   % 在坐标(year(end-1)+0.2,0.55)上添加文本
    set(gca,'xtick',year(2:1:end))  % 设置x轴横坐标的间隔为1
    xlabel('年份');  ylabel('原始数据的光滑度');  % 给坐标轴加上标签
    
    
    disp(strcat('指标1：光滑比小于0.5的数据占比为',num2str(100*sum(rho<0.5)/(n-1)),'%'))
    disp(strcat('指标2：除去前两个时期外，光滑比小于0.5的数据占比为',num2str(100*sum(rho(3:end)<0.5)/(n-3)),'%'))
    disp('参考标准：指标1一般要大于60%, 指标2要大于90%，你认为本例数据可以通过检验吗？')
    
    Judge = input('你认为可以通过准指数规律的检验吗？可以通过请输入1，不能请输入0：');
    if Judge == 0
        disp('亲，灰色预测模型不适合你的数据哦~ 请考虑其他方法吧 例如ARIMA，指数平滑等')
        ERROR = 1;
    end
    disp('------------------------------------------------------------')
end

%% 当数据量大于4时，我们利用试验组来选择使用传统的GM(1,1)模型、新信息GM(1,1)模型还是新陈代谢GM(1,1)模型； 如果数据量等于4，那么我们直接对三种方法求一个平均来进行预测
if ERROR == 0   % 如果上述错误均没有发生时，才能执行下面的操作步骤
    if  n > 4  % 数据量大于4时，将数据分为训练组和试验组(根据原数据量大小n来取，n为5-7个则取最后两年为试验组，n大于7则取最后三年为试验组)
        disp('因为原数据的期数大于4，所以我们可以将数据组分为训练组和试验组')   % 注意，如果试验组的个数只有1个，那么三种模型的结果完全相同，因此至少要取2个试验组
        if n > 7
            test_num = 3;
        else
            test_num = 2;
        end
        train_x0 = x0(1:end-test_num);  % 训练数据
        disp('训练数据是: ')
        disp(mat2str(train_x0'))  % mat2str可以将矩阵或者向量转换为字符串显示, 这里加一撇表示转置，把列向量变成行向量方便观看
        test_x0 =  x0(end-test_num+1:end); % 试验数据
        disp('试验数据是: ')
        disp(mat2str(test_x0'))  % mat2str可以将矩阵或者向量转换为字符串显示
        disp('------------------------------------------------------------')
        
        % 使用三种模型对训练数据进行训练，返回的result就是往后预测test_num期的数据
        disp(' ')
        disp('***下面是传统的GM(1,1)模型预测的详细过程***')
        result1 = gm11(train_x0, test_num); %使用传统的GM(1,1)模型对训练数据，并预测后test_num期的结果
        disp(' ')
        disp('***下面是进行新信息的GM(1,1)模型预测的详细过程***')
        result2 = new_gm11(train_x0, test_num); %使用新信息GM(1,1)模型对训练数据，并预测后test_num期的结果
        disp(' ')
        disp('***下面是进行新陈代谢的GM(1,1)模型预测的详细过程***')
        result3 = metabolism_gm11(train_x0, test_num); %使用新陈代谢GM(1,1)模型对训练数据，并预测后test_num期的结果
        
        % 现在比较三种模型对于试验数据的预测结果
        disp(' ')
        disp('------------------------------------------------------------')
        % 绘制对试验数据进行预测的图形（对于部分数据，可能三条直线预测的结果非常接近）
        test_year = year(end-test_num+1:end);  % 试验组对应的年份
        figure(3)
        plot(test_year,test_x0,'o-',test_year,result1,'*-',test_year,result2,'+-',test_year,result3,'x-'); grid on;
        set(gca,'xtick',year(end-test_num+1): 1 :year(end))  % 设置x轴横坐标的间隔为1
        legend('试验组的真实数据','传统GM(1,1)预测结果','新信息GM(1,1)预测结果','新陈代谢GM(1,1)预测结果')  % 注意：如果lengend挡着了图形中的直线，那么lengend的位置可以自己手动拖动
        xlabel('年份');  ylabel('排污总量');  % 给坐标轴加上标签
        % 计算误差平方和SSE
        SSE1 = sum((test_x0-result1).^2);
        SSE2 = sum((test_x0-result2).^2);
        SSE3 = sum((test_x0-result3).^2);
        disp(strcat('传统GM(1,1)对于试验组预测的误差平方和为',num2str(SSE1)))
        disp(strcat('新信息GM(1,1)对于试验组预测的误差平方和为',num2str(SSE2)))
        disp(strcat('新陈代谢GM(1,1)对于试验组预测的误差平方和为',num2str(SSE3)))
        if SSE1<SSE2
            if SSE1<SSE3
                choose = 1;  % SSE1最小，选择传统GM(1,1)模型
            else
                choose = 3;  % SSE3最小，选择新陈代谢GM(1,1)模型
            end
        elseif SSE2<SSE3
            choose = 2;  % SSE2最小，选择新信息GM(1,1)模型
        else
            choose = 3;  % SSE3最小，选择新陈代谢GM(1,1)模型
        end
        Model = {'传统GM(1,1)模型','新信息GM(1,1)模型','新陈代谢GM(1,1)模型'};
        disp(strcat('因为',Model(choose),'的误差平方和最小，所以我们应该选择其进行预测'))
        disp('------------------------------------------------------------')
        
        %% 选用误差最小的那个模型进行预测
        predict_num = input('请输入你要往后面预测的期数： ');
        % 计算使用传统GM模型的结果，用来得到另外的返回变量：x0_hat, 相对残差relative_residuals和级比偏差eta
        [result, x0_hat, relative_residuals, eta] = gm11(x0, predict_num);  % 先利用gm11函数得到对原数据拟合的详细结果
        
        % % 判断我们选择的是哪个模型，如果是2或3，则更新刚刚由模型1计算出来的预测结果
        if choose == 2
            result = new_gm11(x0, predict_num);
        end
        if choose == 3
            result = metabolism_gm11(x0, predict_num);
        end
        
        %% 输出使用最佳的模型预测出来的结果
        disp('------------------------------------------------------------')
        disp('对原始数据的拟合结果：')
        for i = 1:n
            disp(strcat(num2str(year(i)), ' ： ',num2str(x0_hat(i))))
        end
        disp(strcat('往后预测',num2str(predict_num),'期的得到的结果：'))
        for i = 1:predict_num
            disp(strcat(num2str(year(end)+i), ' ： ',num2str(result(i))))
        end
        
        %% 如果只有四期数据，那么我们就没必要选择何种模型进行预测，直接对三种模型预测的结果求一个平均值~
    else
        disp('因为数据只有4期，因此我们直接将三种方法的结果求平均即可~')
        predict_num = input('请输入你要往后面预测的期数： ');
        disp(' ')
        disp('***下面是传统的GM(1,1)模型预测的详细过程***')
        [result1, x0_hat, relative_residuals, eta] = gm11(x0, predict_num);
        disp(' ')
        disp('***下面是进行新信息的GM(1,1)模型预测的详细过程***')
        result2 = new_gm11(x0, predict_num);
        disp(' ')
        disp('***下面是进行新陈代谢的GM(1,1)模型预测的详细过程***')
        result3 = metabolism_gm11(x0, predict_num);
        result = (result1+result2+result3)/3;
        
        disp('对原始数据的拟合结果：')
        for i = 1:n
            disp(strcat(num2str(year(i)), ' ： ',num2str(x0_hat(i))))
        end
        disp(strcat('传统GM(1,1)往后预测',num2str(predict_num),'期的得到的结果：'))
        for i = 1:predict_num
            disp(strcat(num2str(year(end)+i), ' ： ',num2str(result1(i))))
        end
        disp(strcat('新信息GM(1,1)往后预测',num2str(predict_num),'期的得到的结果：'))
        for i = 1:predict_num
            disp(strcat(num2str(year(end)+i), ' ： ',num2str(result2(i))))
        end
        disp(strcat('新陈代谢GM(1,1)往后预测',num2str(predict_num),'期的得到的结果：'))
        for i = 1:predict_num
            disp(strcat(num2str(year(end)+i), ' ： ',num2str(result3(i))))
        end
        disp(strcat('三种方法求平均得到的往后预测',num2str(predict_num),'期的得到的结果：'))
        for i = 1:predict_num
            disp(strcat(num2str(year(end)+i), ' ： ',num2str(result(i))))
        end
    end
    
    %% 绘制相对残差和级比偏差的图形(注意：因为是对原始数据的拟合效果评估，所以三个模型都是一样的哦~~~)
    figure(4)
    subplot(2,1,1)  % 绘制子图（将图分块）
    plot(year(2:end), relative_residuals,'*-'); grid on;   % 原数据中的各时期和相对残差
    legend('相对残差'); xlabel('年份');
    set(gca,'xtick',year(2:1:end))  % 设置x轴横坐标的间隔为1
    subplot(2,1,2)
    plot(year(2:end), eta,'o-'); grid on;   % 原数据中的各时期和级比偏差
    legend('级比偏差'); xlabel('年份');
    set(gca,'xtick',year(2:1:end))  % 设置x轴横坐标的间隔为1
    disp(' ')
    disp('****下面将输出对原数据拟合的评价结果***')
    
    %% 残差检验
    average_relative_residuals = mean(relative_residuals);  % 计算平均相对残差 mean函数用来均值
    disp(strcat('平均相对残差为',num2str(average_relative_residuals)))
    if average_relative_residuals<0.1
        disp('残差检验的结果表明：该模型对原数据的拟合程度非常不错')
    elseif average_relative_residuals<0.2
        disp('残差检验的结果表明：该模型对原数据的拟合程度达到一般要求')
    else
        disp('残差检验的结果表明：该模型对原数据的拟合程度不太好，建议使用其他模型预测')
    end
    
    %% 级比偏差检验
    average_eta = mean(eta);   % 计算平均级比偏差
    disp(strcat('平均级比偏差为',num2str(average_eta)))
    if average_eta<0.1
        disp('级比偏差检验的结果表明：该模型对原数据的拟合程度非常不错')
    elseif average_eta<0.2
        disp('级比偏差检验的结果表明：该模型对原数据的拟合程度达到一般要求')
    else
        disp('级比偏差检验的结果表明：该模型对原数据的拟合程度不太好，建议使用其他模型预测')
    end
    disp(' ')
    disp('------------------------------------------------------------')
    
    %% 绘制最终的预测效果图
    figure(5)  % 下面绘图中的符号m:洋红色 b:蓝色
    plot(year,x0,'-o',  year,x0_hat,'-*m',  year(end)+1:year(end)+predict_num,result,'-*b' );   grid on;
    hold on;
    plot([year(end),year(end)+1],[x0(end),result(1)],'-*b')
    legend('原始数据','拟合数据','预测数据')  % 注意：如果lengend挡着了图形中的直线，那么lengend的位置可以自己手动拖动
    set(gca,'xtick',[year(1):1:year(end)+predict_num])  % 设置x轴横坐标的间隔为1
    xlabel('年份');  ylabel('排污总量');  % 给坐标轴加上标签

function [result, x0_hat, relative_residuals, eta] = gm11(x0, predict_num)
    % 函数作用：使用传统的GM(1,1)模型对数据进行预测
    %     x0：要预测的原始数据
    %     predict_num： 向后预测的期数
    % 输出变量 （注意，实际调用时该函数时不一定输出全部结果，就像corrcoef函数一样~，可以只输出相关系数矩阵，也可以附带输出p值矩阵）
    %     result：预测值
    %     x0_hat：对原始数据的拟合值
    %     relative_residuals： 对模型进行评价时计算得到的相对残差
    %     eta： 对模型进行评价时计算得到的级比偏差

    n = length(x0); % 数据的长度
    x1=cumsum(x0); % 计算一次累加值
    z1 = (x1(1:end-1) + x1(2:end)) / 2;  % 计算紧邻均值生成数列（长度为n-1）
    % 将从第二项开始的x0当成y，z1当成x，来进行一元回归  y = kx +b
    y = x0(2:end); x = z1;
    % 下面的表达式就是第四讲拟合里面的哦~ 但是要注意，此时的样本数应该是n-1，少了一项哦
    k = ((n-1)*sum(x.*y)-sum(x)*sum(y))/((n-1)*sum(x.*x)-sum(x)*sum(x));
    b = (sum(x.*x)*sum(y)-sum(x)*sum(x.*y))/((n-1)*sum(x.*x)-sum(x)*sum(x));
    a = -k;  %注意：k = -a哦
    % 注意： -a就是发展系数,  b就是灰作用量
    
    disp('现在进行GM(1,1)预测的原始数据是: ')
    disp(mat2str(x0'))  % mat2str可以将矩阵或者向量转换为字符串显示
    disp(strcat('最小二乘法拟合得到的发展系数为',num2str(-a),'，灰作用量是',num2str(b)))
    disp('***************分割线***************')
    x0_hat=zeros(n,1);  x0_hat(1)=x0(1);   % x0_hat向量用来存储对x0序列的拟合值，这里先进行初始化
    for m = 1: n-1
        x0_hat(m+1) = (1-exp(a))*(x0(1)-b/a)*exp(-a*m);
    end
    result = zeros(predict_num,1);  % 初始化用来保存预测值的向量
    for i = 1: predict_num
        result(i) = (1-exp(a))*(x0(1)-b/a)*exp(-a*(n+i-1)); % 带入公式直接计算
    end

    % 计算绝对残差和相对残差
    absolute_residuals = x0(2:end) - x0_hat(2:end);   % 从第二项开始计算绝对残差，因为第一项是相同的
    relative_residuals = abs(absolute_residuals) ./ x0(2:end);  % 计算相对残差，注意分子要加绝对值，而且要使用点除
    % 计算级比和级比偏差
    class_ratio = x0(2:end) ./ x0(1:end-1) ;  % 计算级比 sigma(k) = x0(k)/x0(k-1)
    eta = abs(1-(1-0.5*a)/(1+0.5*a)*(1./class_ratio));  % 计算级比偏差
end

function [result] = new_gm11(x0, predict_num)
% 函数作用：使用新信息的GM(1,1)模型对数据进行预测
% 输入变量
%     x0：要预测的原始数据
%     predict_num： 向后预测的期数
% 输出变量
%     result：预测值
    result = zeros(predict_num,1);  % 初始化用来保存预测值的向量
    for i = 1 : predict_num  
        result(i) = gm11(x0, 1);  % 将预测一期的结果保存到result中
        x0 = [x0; result(i)];  % 更新x0向量，此时x0多了新的预测信息
    end
end

function [result] = metabolism_gm11(x0, predict_num)
% 函数作用：使用新陈代谢的GM(1,1)模型对数据进行预测
% 输入变量
%     x0：要预测的原始数据
%     predict_num： 向后预测的期数
% 输出变量
%     result：预测值
    result = zeros(predict_num,1);  % 初始化用来保存预测值的向量
    for i = 1 : predict_num  
        result(i) = gm11(x0, 1);  % 将预测一期的结果保存到result中
        x0 = [x0(2:end); result(i)];  % 更新x0向量，此时x0多了新的预测信息，并且删除了最开始的那个向量
    end
end