MATLAB实现灰色预测

        久违了，前段时间由于学习压力大，就没怎么更新MATLAB相关的内容，今天实在学不进去了，换个内容更新一下~

        本贴介绍灰色预测模型，这也是数学建模竞赛常见算法中的一员，和许多预测模型一样——底层原理是根据已知数据对未知进行预测~

一.理论部分

二.例题讲解

如下是有关棉花生产的数据：

1.传统GM模型的代码

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;  
    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;  
    x0_hat=zeros(n,1);  x0_hat(1)=x0(1);   
    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) ;  
    eta = abs(1-(1-0.5*a)/(1+0.5*a)*(1./class_ratio));  
end

2.原始数据检验

此处选用【单产】作为示例~

year =[1995:1:2001]';  % 横坐标表示年份，写成列向量的形式（加'就表示转置）
yield= [1017
1036.5
792
861
901.5
922.5
916.5
976.5
1024.5
1003.5
1069.5
1168.5
1228.5
1023
1144.5
]';  %原始数据序列，写成列向量的形式（加'就表示转置）

ERROR = 0;  % 建立一个错误指标，一旦出错就指定为1
% 判断是否有负数元素
if sum(yield<0) > 0 
    disp('灰色预测的时间序列中不能有负数！')
    ERROR = 1;
end

% 判断数据量是否太少
n = length(yield);  % 计算原始数据的长度
disp(strcat('原始数据的长度为',num2str(n)))    
if n<=3
    disp('数据量太小')
    ERROR = 1;
end

% 数据太多时提示可考虑使用其他方法
if n>10
    disp('考虑使用其他的方法')
end

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

3.准指数规律检验

if ERROR == 0    
    disp('准指数规律检验')
    x1 = cumsum(yield);  
    rho = yield(2:end) ./ x1(1:end-1) ;   
    rho
   
    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,'临界线')   
    set(gca,'xtick',year(2:1:end)) 
    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('灰色预测模型不适合你的数据！')
        ERROR = 1;
    end
end

4.传统的GM（1,1）预测

        if n > 7
            test_num = 3;
        else
            test_num = 2;
        end
        train_yield = yield(1:end-test_num);  
        disp('训练数据是: ')
        disp(mat2str(train_yield')) 
        test_yield = yield(end-test_num+1:end); 
        disp('试验数据是: ')
        disp(mat2str(test_yield')) 

disp(' ')
   disp('***下面是传统的GM(1,1)模型预测的详细过程***')
   result1 = gm11(train_yield, test_num);

5.评估误差精度

    %% 残差检验
    average_relative_residuals = mean(relative_residuals); 
    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(' ')

答案如下，大家自己尝试（每个单独预测一遍，因为GM（1,1）只针对一维数据~）

三.实战案例

1.2022年美赛C题

        根据现有数据，通过灰色预测模型，预测比特币和黄金两种波动资产的走向，拟合优度较高，残差与级比偏差均很低，预测模型的可信度较高~

2.2022亚太赛C题

上图是准指数检验~

 通过BP神经网络、多元线性回归以及灰色预测3种方式预测气温变化~结果可信度较高

代码007普通

打赏 收藏 海报 链接