灰色系統

灰色預測是對既含有已知信息又含有不確定信息的系統進行預測，就是對在一定范圍內變化的、與時間有關的灰色過程進行預測。灰色預測對原始數據進行生成處理來尋找系統變動的規律，并生成有較強規律性的數據序列，然后建立相應的微分方程模型，從而預測事物未來發展趨勢的狀況。

GM(1,1)模型：Grey(Gray)Model

GM(1,1)是使用原始的離散非負數據列，通過一次累加生成削弱隨機性的
較有規律的新的離散數據列，然后通過建立微分方程模型，得到在離散點處的解經過累減生成的原始數據的近似估計值，從而預測原始數據的后續發展。
注：第一個‘1’表示微分方程是一階的，后面的‘1’表示只有一個變量。

GM(1,1)原理介紹

其中累加數據為相鄰兩項相加的和，緊鄰均值生成序列為相鄰兩項相加并求均值。

我們利用OLS就可以得到a和b，然后就可以將灰色微分方程：
轉化為白化方程：

擬合值的計算：

準指數規律的檢驗

GM(1,1)模型的檢驗

1.殘差檢驗

2.級比偏差檢驗

GM(1,1)模型的拓展

GM(1,1)模型代碼的步驟

畫出原始數據的時間序列圖，并判斷原始數據中是否有負數或期數是否低于4期，如果是的話則報錯，否則執行下一步；

對一次累加后的數據進行準指數規律檢驗，返回兩個指標：
指標1：光滑比小于0.5的數據占比（一般要大于60%）
指標2：除去前兩個時期外，光滑比小于0.5的數據占比（一般大于90%）
并讓用戶決定數據是否滿足準指數規律，滿足則輸入1，不滿足則輸入0

如果上一步用戶輸入0，則程序停止；如果輸入1，則繼續下面的步驟。

讓用戶輸入需要預測的后續期數，并判斷原始數據的期數：
4.1 數據期數為4：
分別計算出傳統的GM(1,1)模型、新信息GM(1,1)模型和新陳代謝GM(1,1)模型對于未來期數的預測結果，為了保證結果的穩健性，對三個結果求平均值作為預測值。
4.2 數據期數為5,6或7：
取最后兩期為試驗組，前面的n‐2期為訓練組；用訓練組的數據分別訓練三種GM模型，并將訓練出來的模型分別用于預測試驗組的兩期數據；利用試驗組兩期的真實數據和預測出來的兩期數據，可分別計算出三個模型的SSE；選擇SSE最小的模型作為我們建模的模型。
4.3 數據期數大于7：
取最后三期為試驗組，其他的過程和4.2類似。

輸出并繪制圖形顯示預測結果，并進行殘差檢驗和級比偏差檢驗

例題

主函數

%% 輸入原始數據并做出時間序列圖 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 % x0<0返回一個邏輯數組(0-1組成)，如果有數據小于0，則所在位置為1，如果原始數據均為非負數，那么這個邏輯數組中全為0，求和后也是0~disp('親，灰色預測的時間序列中不能有負數哦')ERROR = 1; end% 判斷數據量是否太少 n = length(x0); % 計算原始數據的長度 disp(strcat('原始數據的長度為',num2str(n))) % strcat()是連接字符串的函數，第一講學了，可別忘了哦 if n<=3disp('親，數據量太小，我無能為力哦')ERROR = 1; end% 數據太多時提示可考慮使用其他方法（不報錯） if n>10disp('親，這么多數據量，一定要考慮使用其他的方法哦，例如ARIMA，指數平滑等') end% 判斷數據是否為列向量，如果輸入的是行向量則轉置為列向量 if size(x0,1) == 1x0 = x0'; end if size(year,1) == 1year = year'; end%% 對一次累加后的數據進行準指數規律的檢驗(注意，這個檢驗有時候即使能通過，也不一定能保證預測結果非常好，例如上面的第三組數據) if ERROR == 0 % 如果上述錯誤均沒有發生時，才能執行下面的操作步驟disp('------------------------------------------------------------')disp('準指數規律檢驗')x1 = cumsum(x0); % 生成1-AGO序列，cumsum是累加函數哦~ 注意：1.0e+03 *0.1740的意思是科學計數法,即10^3*0.1740 = 174rho = 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軸橫坐標的間隔為1xlabel('年份'); 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('親，灰色預測模型不適合你的數據哦~ 請考慮其他方法吧 例如ARIMA，指數平滑等')ERROR = 1;enddisp('------------------------------------------------------------') 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 > 7test_num = 3;elsetest_num = 2;endtrain_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軸橫坐標的間隔為1legend('試驗組的真實數據','傳統GM(1,1)預測結果','新信息GM(1,1)預測結果','新陳代謝GM(1,1)預測結果') % 注意：如果lengend擋著了圖形中的直線，那么lengend的位置可以自己手動拖動xlabel('年份'); ylabel('排污總量'); % 給坐標軸加上標簽% 計算誤差平方和SSESSE1 = 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<SSE2if SSE1<SSE3choose = 1; % SSE1最小，選擇傳統GM(1,1)模型elsechoose = 3; % SSE3最小，選擇新陳代謝GM(1,1)模型endelseif SSE2<SSE3choose = 2; % SSE2最小，選擇新信息GM(1,1)模型elsechoose = 3; % SSE3最小，選擇新陳代謝GM(1,1)模型endModel = {'傳統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 == 2result = new_gm11(x0, predict_num);endif choose == 3result = metabolism_gm11(x0, predict_num);end%% 輸出使用最佳的模型預測出來的結果disp('------------------------------------------------------------')disp('對原始數據的擬合結果：')for i = 1:ndisp(strcat(num2str(year(i)), ' ： ',num2str(x0_hat(i))))enddisp(strcat('往后預測',num2str(predict_num),'期的得到的結果：'))for i = 1:predict_numdisp(strcat(num2str(year(end)+i), ' ： ',num2str(result(i))))end%% 如果只有四期數據，那么我們就沒必要選擇何種模型進行預測，直接對三種模型預測的結果求一個平均值~elsedisp('因為數據只有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:ndisp(strcat(num2str(year(i)), ' ： ',num2str(x0_hat(i))))enddisp(strcat('傳統GM(1,1)往后預測',num2str(predict_num),'期的得到的結果：'))for i = 1:predict_numdisp(strcat(num2str(year(end)+i), ' ： ',num2str(result1(i))))enddisp(strcat('新信息GM(1,1)往后預測',num2str(predict_num),'期的得到的結果：'))for i = 1:predict_numdisp(strcat(num2str(year(end)+i), ' ： ',num2str(result2(i))))enddisp(strcat('新陳代謝GM(1,1)往后預測',num2str(predict_num),'期的得到的結果：'))for i = 1:predict_numdisp(strcat(num2str(year(end)+i), ' ： ',num2str(result3(i))))enddisp(strcat('三種方法求平均得到的往后預測',num2str(predict_num),'期的得到的結果：'))for i = 1:predict_numdisp(strcat(num2str(year(end)+i), ' ： ',num2str(result(i))))endend%% 繪制相對殘差和級比偏差的圖形(注意：因為是對原始數據的擬合效果評估，所以三個模型都是一樣的哦~~~)figure(4)subplot(2,1,1) % 繪制子圖（將圖分塊）plot(year(2:end), relative_residuals,'*-'); grid on; % 原數據中的各時期和相對殘差legend('相對殘差'); xlabel('年份');set(gca,'xtick',year(2:1:end)) % 設置x軸橫坐標的間隔為1subplot(2,1,2)plot(year(2:end), eta,'o-'); grid on; % 原數據中的各時期和級比偏差legend('級比偏差'); xlabel('年份');set(gca,'xtick',year(2:1:end)) % 設置x軸橫坐標的間隔為1disp(' ')disp('****下面將輸出對原數據擬合的評價結果***')%% 殘差檢驗average_relative_residuals = mean(relative_residuals); % 計算平均相對殘差 mean函數用來均值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(' ')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軸橫坐標的間隔為1xlabel('年份'); ylabel('排污總量'); % 給坐標軸加上標簽 end

傳統GM(1,1): gm11

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 +by = 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-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)); % 帶入公式直接計算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

新信息的GM(1,1): new_gm11

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

新陳代謝的GM(1,1): metabolism_gm11

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

總結

以上是生活随笔為你收集整理的预测模型（数学建模）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。