????????久違了，前段時間由于學習壓力大，就沒怎么更新MATLAB相關的內容，今天實在學不進去了，換個內容更新一下~

????????本貼介紹灰色預測模型，這也是數學建模競賽常見算法中的一員，和許多預測模型一樣——底層原理是根據已知數據對未知進行預測~

---

一.理論部分

二.例題講解

如下是有關棉花生產的數據：

年份	單產	種子費	化肥費	農藥費	機械費	灌溉費
	kg/公頃	元/公頃	元/公頃	元/公頃	元/公頃	元/公頃
1990	1017	106.05	495.15	305.1	45.9	56.1
1991	1036.5	113.55	561.45	343.8	68.55	93.3
1992	792	104.55	584.85	414	73.2	104.55
1993	861	132.75	658.35	453.75	82.95	107.55
1994	901.5	174.3	904.05	625.05	114	152.1
1995	922.5	230.4	1248.75	834.45	143.85	176.4
1996	916.5	238.2	1361.55	720.75	165.15	194.25
1997	976.5	260.1	1337.4	727.65	201.9	291.75
1998	1024.5	270.6	1195.8	775.5	220.5	271.35
1999	1003.5	286.2	1171.8	610.95	195	284.55
2000	1069.5	282.9	1151.55	599.85	190.65	277.35
2001	1168.5	317.85	1105.8	553.8	211.05	290.1
2002	1228.5	319.65	1213.05	513.75	231.6	324.15
2003	1023	368.4	1274.1	567.45	239.85	331.8
2004	1144.5	466.2	1527.9	487.35	408	336.15

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-1x0_hat(m+1) = (1-exp(a))*(x0(1)-b/a)*exp(-a*m);endresult = zeros(predict_num,1);  for i = 1: predict_numresult(i) = (1-exp(a))*(x0(1)-b/a)*exp(-a*(n+i-1)); endabsolute_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<=3disp('數據量太小')ERROR = 1;
end% 數據太多時提示可考慮使用其他方法
if n>10disp('考慮使用其他的方法')
end% 判斷數據是否為列向量，如果輸入的是行向量則轉置為列向量
if size(yield,1) == 1yield = yield';
end
if size(year,1) == 1year = year';
end

3.準指數規律檢驗

if ERROR == 0    disp('準指數規律檢驗')x1 = cumsum(yield);  rho = yield(2:end) ./ x1(1:end-1) ;   rhofigure(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 == 0disp('灰色預測模型不適合你的數據！')ERROR = 1;end
end

4.傳統的GM（1,1）預測

        if n > 7test_num = 3;elsetest_num = 2;endtrain_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.1disp('該模型對原數據的擬合程度非常不錯！')elseif average_relative_residuals<0.2disp('該模型對原數據的擬合程度達到一般要求！')elsedisp('該模型對原數據的擬合程度不太好！')end%% 級比偏差檢驗average_eta = mean(eta);   % 計算平均級比偏差disp(strcat('平均級比偏差為',num2str(average_eta)))if average_eta<0.1disp('該模型對原數據的擬合程度非常不錯！')elseif average_eta<0.2disp('該模型對原數據的擬合程度達到一般要求！')elsedisp('該模型對原數據的擬合程度不太好！')enddisp(' ')

答案如下，大家自己嘗試（每個單獨預測一遍，因為GM（1,1）只針對一維數據~）

年份	單產	種子費	化肥費	農藥費	機械費	灌溉費
	kg/公頃	元/公頃	元/公頃	元/公頃	元/公頃	元/公頃
2005	1122	449.85	1703.25	555.15	402.3	358.8
2006	1276.5	537	1888.5	637.2	480.75	428.4
2007	1233	565.5	2009.85	715.65	562.05	456.9

?

三.實戰案例

1.2022年美賽C題

?

????????根據現有數據，通過灰色預測模型，預測比特幣和黃金兩種波動資產的走向，擬合優度較高，殘差與級比偏差均很低，預測模型的可信度較高~

2.2022亞太賽C題

上圖是準指數檢驗~

?通過BP神經網絡、多元線性回歸以及灰色預測3種方式預測氣溫變化~結果可信度較高