预测方法

预测类算法就是一种通过学习历史数据来预测未来事件或趋势的方法。它会根据过去的经验和规律，根据已知的数据来做出预测，帮助我们预测未来可能会发生的事情。这里记录几种在数学建模里较为常见且泛用性较高的算法。

先列出一些常见数据的类型：

1. 横截面数据（Cross-Sectional Data）：是在某一特定时间点上收集的数据。它包括不同个体或单位在同一时间点上的观察结果。例如，在一项调查中，收集了不同家庭在2019年的收入数据，这些数据就构成了横截面数据。通常用于描述群体特征、进行比较分析等。可以采用多元回归分析。

2. 时间序列数据（Time Series Data）：是按照时间顺序记录的数据。它包括同一变量在不同时间点上的观察结果，可以用来分析随时间变化的趋势和模式。例如，某公司每个季度的销售额数据就构成了时间序列数据。可以采用时间序列分析。

3. 面板数据（Panel Data）：结合了横截面数据和时间序列数据的特点，包括了多个个体（横截面）在多个时间点上的观察结果。面板数据可用于研究个体差异以及随时间的变化。例如，对不同地区居民的收入情况进行跟踪调查，即包含了不同个体在不同时间点上的数据。很少见，遇到了再说。

回归分析方法

回归分析是一种用来预测或解释变量之间关系的统计方法。它的主要思想是通过观察已知的数据来建立一个数学模型，然后利用这个模型来预测未来的结果或者解释变量之间的关系，根据现有的数据来预测未来的结果。对于决策和规划都非常有帮助。

回归分析应该识别重要变量，判断相关性的方向，并且估计权重（回归系数）。

首先是简单线性回归。假设我们有两个变量，一个是用来帮助预测的变量X（称为自变量）,另一个是我们想要预测的变量Y（称为因变量）。比如，我们想要预测一个人的体重，可能会用身高作为帮助预测的变量。简单线性回归就是要找到一个直线方程，通过自变量（比如身高）来预测因变量（比如体重）。这个直线方程可以告诉我们，随着自变量的变化，因变量会怎么变化。

简单线性回归虽然可以用来描述变量之间的相关性，但却不能证明因果性。两个变量之间的相关性并不一定意味着它们之间存在因果关系，也就是说，一个变量的变化是否会导致另一个变量的变化，需要更进一步的研究来确定。举例来说，假设我们发现一个城市的人均收入和房价呈现正相关，即人均收入越高，该城市的房价也越高。这种相关性看起来很有道理，但我们不能因此就断定高收入导致了高房价，或者是高房价影响了人们的收入水平。实际上还可能存在其他因素，比如地理位置、人口密度等，也可能影响了这两个变量之间的关系。

让我们来举一个更极端的例子。假设我们观察到某地区消防队员人数与火灾数量之间存在正相关关系，也就是说，我们得出的结论是随着消防队员人数的增加，火灾的数量也在增加。因此我们应该减少消防人员的数量，这样火灾就会越来越少了！很显然，这个结论存在相关性，但没有因果性，我们不能因此得出消防队员人数的增加导致了火灾数量的增加的结论。实际上，更可能的解释是，这种相关性是由于第三个隐藏变量引起的，比如人口密度、城市规模等因素。另外，也可能是由于数据的选择或者测量方法等问题造成了这种看似荒谬的相关性。

因此，我们更应该进行多元回归分析，在简单线性回归的基础上考虑了更多的自变量。比如，我们想预测一个人的收入，不仅仅会考虑教育水平，还会考虑工作经验、所在地区等因素。多元回归可以帮助我们建立一个更复杂的模型，来解释所有这些自变量对因变量的影响。

多元回归分析中，存在内生性的问题。具体来说，当解释变量与误差项之间存在内生性时，解释变量不再是外生的，即它们不再是独立于误差项的。这意味着解释变量的变化可能是由于未观测到的因素或其他内部因素引起的，而不仅仅是外部因素的结果。内生性问题可能导致统计模型中的参数估计不一致，使得我们无法准确地推断解释变量对因变量的真实影响。它还可能模糊或混淆解释变量与因变量之间的因果关系，使得我们难以确定哪个是因果变量，哪个是结果变量。

核心解释变量通常是研究中我们最感兴趣的变量，也是对因果关系进行推断的变量。在一个统计模型中，核心解释变量通常是被认为对因变量有直接影响的变量。控制变量则是指在研究中被引入并控制在统计模型中的其他变量。它们被用来减少内生性问题的影响，以确保核心解释变量与因变量之间的关系是被准确估计的。通过引入控制变量，可以尽可能地排除其他可能对因变量产生影响的变量，从而更准确地衡量核心解释变量的影响。应保证核心解释变量不存在内生性。

进行回归分析,代码如下：

% 读取 Excel 文件数据
data = xlsread('高数成绩数据.xlsx', 'Sheet1', 'A2:E');

% 提取自变量和因变量
X = data(:, 2:end); % 自变量：高考数学、高考总分、班干与否、平时成绩
y = data(:, 1); % 因变量：期末成绩

% 进行多元线性回归分析
mdl = fitlm(X, y);

% 打印回归系数
disp(mdl.Coefficients);

% 打印回归统计信息
disp(mdl.Rsquared);
disp(sqrt(mdl.MSE));

或者使用stata软件，导入后使用:

clear
cls 
// 导入数据
// import excel "奶粉数据.xlsx", sheet("Sheet1") firstrow
import excel "奶粉数据.xlsx", sheet("Sheet1") firstrow
// 定量变量的描述性统计

summarize 团购价元 评价量 商品毛重kg
// 定性变量的频数分布，并得到相应字母开头的虚拟变量
tabulate 配方,gen(A)
tabulate 奶源产地 ,gen(B)
tabulate 国产或进口 ,gen(C)
tabulate 适用年龄岁 ,gen(D)
tabulate 包装单位 ,gen(E)
tabulate 分类 ,gen(F)
tabulate 段位 ,gen(G)

// Stata会自动剔除多重共线性的变量
regress 评价量 团购价元 商品毛重kg A1 A2 A3 B1 B2 B3 B4 B5 B6 B7 B8 B9 C1 C2 D1 D2 D3 D4 D5 E1 E2 E3 E4 F1 F2 G1 G2 G3 G4
est store m2
reg2docx m2 using m2.docx, replace

// 得到标准化回归系数
regress 评价量 团购价元 商品毛重kg, b 

// 画出残差图
regress 评价量 团购价元 商品毛重kg A1 A2 A3 B1 B2 B3 B4 B5 B6 B7 B8 B9 C1 C2 D1 D2 D3 D4 D5 E1 E2 E3 E4 F1 F2 G1 G2 G3 G4
rvfplot 
// 残差与拟合值的散点图
graph export a1.png ,replace
// 残差与自变量团购价的散点图
rvpplot  团购价元
graph export a2.png ,replace

// 描述性统计并给出分位数对应的数值
summarize 评价量,d

// 作评价量的概率密度估计图
kdensity 评价量 
graph export a3.png ,replace

// 异方差BP检验
estat hettest ,rhs iid

// 异方差怀特检验
estat imtest,white

// 使用OLS + 稳健的标准误
regress 评价量 团购价元 商品毛重kg A1 A2 A3 B1 B2 B3 B4 B5 B6 B7 B8 B9 C1 C2 D1 D2 D3 D4 D5 E1 E2 E3 E4 F1 F2 G1 G2 G3 G4, r
est store m3
reg2docx m3 using m3.docx, replace

// 计算VIF
estat  vif

// 逐步回归（一定要注意完全多重共线性的影响）
// 向前逐步回归（后面的r表示稳健的标准误）
stepwise reg 评价量 团购价元 商品毛重kg A1 A3 B1 B2 B3 B4 B5 B6 B7 B9 C1 D1 D2 D3 D4 E1 E2 E3 F1 G1 G2 G3,  r pe(0.05)
// 向后逐步回归（后面的r表示稳健的标准误）
stepwise reg 评价量 团购价元 商品毛重kg A1 A3 B1 B2 B3 B4 B5 B6 B7 B9 C1 D1 D2 D3 D4 E1 E2 E3 F1 G1 G2 G3,  r pr(0.05)
// 向后逐步回归的同时使用标准化回归系数（在r后面跟上一个b即可）
stepwise reg 评价量 团购价元 商品毛重kg A1 A3 B1 B2 B3 B4 B5 B6 B7 B9 C1 D1 D2 D3 D4 E1 E2 E3 F1 G1 G2 G3,  r b pr(0.05)

时间序列分析

时间序列分析是一种用来研究随着时间推移而变化的数据的方法。这种数据可能是按照时间顺序记录下来的，比如每天、每月或每年的数据。时间序列分析的目标是揭示数据中存在的模式、趋势和周期性，并且利用这些信息做出预测或者推断。时间序列分析根据时间和数值性质的不同，可以分为时期序列和时点序列。时期序列反映现象在一定时期内发展的结果，而时点序列反映现象在一定时点上的瞬间水平。时期序列可加，时点序列不可加。通常情况下，时间序列分析可以帮助我们回答以下问题：

• 趋势：数据是否存在逐渐增加或减少的趋势？这种趋势是线性的还是非线性的？

• 季节性：数据是否存在重复出现的季节性模式？比如一年中特定月份或季度的销售数据是否有规律性的波动？

• 周期性：数据是否存在长期的循环或者周期性变化？比如经济周期、商业周期等。

• 随机性：数据中是否存在随机的、不规律的波动？

为了回答这些问题，时间序列分析会使用各种统计方法和模型，最常用的是季节分解，指数平滑方法和ARIMA模型。

• 季节分解（Seasonal Decomposition）：季节分解是指将时间序列数据分解为趋势、季节性和残差三个部分的过程。趋势表示长期变化趋势，季节性表示周期性变化，残差则是除去趋势和季节性之后的随机波动部分。这种分解有助于我们更好地理解数据的结构和规律，以便进行预测和分析。

• 指数平滑方法（Exponential Smoothing Method）：指数平滑方法是一种常用的预测方法，适用于没有明显趋势和季节性的时间序列数据。该方法通过对历史数据赋予不同权重，来进行未来数值的预测。简单指数平滑、双指数平滑和三指数平滑是常见的指数平滑方法，它们在对历史数据赋予权重时有不同的考虑因素，适用于不同类型的数据。

• ARIMA模型（Autoregressive Integrated Moving Average Model）：ARIMA模型是一种经典的时间序列分析和预测模型，适用于具有趋势和季节性的时间序列数据。ARIMA模型包括自回归（AR）、差分（I）和移动平均（MA）三个部分，分别代表了数据的自相关、趋势和残差成分。通过对这三个部分的组合，ARIMA模型可以对未来的数值进行预测。

对于时间序列分析，我们将使用季节ARIMA模型（SARIMA），这是表述最全面的时间序列预测模型之一，其他的模型都可以由它简化变型后得到。其建模思想是，将预测对象随时间推移而形成的数据序列看作一个随机序列，时间序列是一组依赖时间*的随机变量，构成该时间序列的单个序列值虽然具有不确定性，但整个序列的变化却是有一定的规律性，可以用数学模型近似描述。这组随机变量所具有的依存关系或自相关性表征了预测对象发展的延续性，而一旦这种自相关性被相应的数学模型描述出来，就可以从时间序列的过去值预测其未来值，通过SARIMA模型可以消除趋势性和季节性，转化为平稳时间序列进行建模。

受限于水平，具体过程完全依赖于spss进行处理，请移步帮助文档。

此处内嵌网页，便于查阅。

这个方法在预测中短期的变化时效果较好，长期的预测可能会遇到外界环境变化的问题。

灰色预测模型

存在三种系统，白色、黑色和灰色系统。它们之间的区别主要在于数据是否具有可预测性，以及对系统内部机制的了解程度。

• 白色系统（White System）：白色系统指的是没有相关性的随机过程或时间序列，其各个时刻的观测值彼此独立且方差相等，因此无法通过历史数据来预测未来的数据。白噪声就是一种典型的白色系统，其在时间轴上呈现为完全随机、无法预测的波动。

• 黑色系统（Black System）：黑色系统指的是具有确定性规律的时间序列，其内部机制已经被充分理解，可以通过历史数据来预测未来的数据。黑色系统中的趋势和周期性成分都是明显的，由此可以建立相应的数学模型。例如，股票市场中的价格波动就属于一种黑色系统。

• 灰色系统（Grey System）：灰色系统指的是既不是白色系统也不是黑色系统的时间序列，其内部机制尚未被完全理解，但对其进行预测和分析仍然有一定的可行性。灰色系统中的趋势和周期性成分相对不太明显，其模型参数通常需要通过数据拟合来确定。

而灰色预测模型是一种对含有不确定因素的系统进行预测的方法。灰色预测通过鉴别系统因素之间发展趋势的相异程度，即进行关联分析，并对原始数据进行生成处理来寻找系统变动的规律，生成有较强规律性的数据序列，然后建立相应的微分方程模型，从而预测事物未来发展趋势的状况。在处理较少的特征值数据时，不需要数据的样本空间足够大，就可以解决历史数据少，序列完整性低，可靠性低的问题，可以生成得到规律较强的生成序列。

我们根据数据量的不同，选择使用传统的GM(1,1)模型、新信息GM(1,1)模型还是新陈代谢GM(1,1)模型进行预测，然后比较三种模型对试验数据的预测结果。

主函数

仅作为例子，代码如下：

clear;clc
year =[1995:1:2004]';  % 横坐标表示年份，写成列向量的形式（加'就表示转置）
x0 = [174,179,183,189,207,234,220.5,256,270,285]';  %原始数据序列，写成列向量的形式（加'就表示转置）


% 画出原始数据的时间序列图
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  
    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('y坐标数据');  % 给坐标轴加上标签
end

相关模型的代码实现

最后给出传统的GM(1,1)模型、新信息GM(1,1)模型还是新陈代谢GM(1,1)模型的代码。

1. 传统的GM(1,1)模型:

• x0：要预测的原始数据
• predict_num： 向后预测的期数
• result：预测值
• x0_hat：对原始数据的拟合值
• relative_residuals： 对模型进行评价时计算得到的相对残差
• eta： 对模型进行评价时计算得到的级比偏差

注意，实际调用时该函数时不一定输出全部结果，就像corrcoef函数一样~，可以只输出相关系数矩阵，也可以附带输出p值矩阵

代码如下：

function [result, x0_hat, relative_residuals, eta] = gm11(x0, predict_num)



    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;
    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

2. 新信息GM(1,1)模型，代码如下：

function [result] = new_gm11(x0, predict_num)
    result = zeros(predict_num,1);  % 初始化用来保存预测值的向量
    for i = 1 : predict_num  
        result(i) = gm11(x0, 1);  % 将预测一期的结果保存到result中
        x0 = [x0; result(i)];  % 更新x0向量，此时x0多了新的预测信息
    end
end

3. 新陈代谢GM(1,1)模型

function [result] = metabolism_gm11(x0, predict_num)

    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

以上。