灰色預(yù)測(cè)算法及相關(guān)程序1引言 2算法的基本原理 2.1GM(1,1)模型： 22.2生成數(shù) 2.2.1累加生成 2.2.2累減生成 3算法的具體實(shí)現(xiàn)流程 3.1算法流程圖 3.2實(shí)現(xiàn)步驟 3.3數(shù)據(jù)準(zhǔn)備與預(yù)處理 4算法程序?qū)崿F(xiàn) 4.1程序使用說明 4.2程序源代碼 64.3程序運(yùn)行 4.3.1程序運(yùn)行及運(yùn)行環(huán)境說明 4.3.2輸入數(shù)據(jù) 4.3.3輸出數(shù)據(jù) 5參考文獻(xiàn) 灰色預(yù)測(cè)算法1引言灰色預(yù)測(cè)(greyprediction)是利用灰色系統(tǒng)理論就灰色系統(tǒng)所作的預(yù)測(cè).灰色系統(tǒng)理論認(rèn)為，盡管系統(tǒng)表象復(fù)雜，數(shù)據(jù)散亂，信息不充分，但作為系統(tǒng)，它必然有整體功能和內(nèi)在規(guī)律，必然是有序的.現(xiàn)有的分析方法大多依據(jù)過去的大量數(shù)據(jù)，按照統(tǒng)計(jì)方法分析其規(guī)律，這樣不僅受數(shù)據(jù)量的限制，而且準(zhǔn)確程度不高.而灰色系統(tǒng)理論把隨機(jī)量看作是在一定范圍內(nèi)變化的灰色量，對(duì)灰色量的處理不是尋求它的統(tǒng)計(jì)規(guī)律和概率分布，而是對(duì)原始數(shù)據(jù)加以處理，將雜亂無章的原始數(shù)據(jù)變?yōu)橐?guī)律性較強(qiáng)的生成數(shù)據(jù)，通過對(duì)生成數(shù)據(jù)建立動(dòng)態(tài)模型，來挖掘系統(tǒng)內(nèi)部信息并充分利用信息進(jìn)行分析預(yù)測(cè).目前，灰色系統(tǒng)理論用于預(yù)測(cè)主要通過GM(m，n)模型，該模型是灰色系統(tǒng)理論的量化體現(xiàn)，可用于以下幾個(gè)方面的預(yù)測(cè)：數(shù)列預(yù)測(cè)：對(duì)某個(gè)事物發(fā)展變化的大小與時(shí)間進(jìn)行預(yù)測(cè).災(zāi)變預(yù)測(cè)：預(yù)測(cè)災(zāi)變發(fā)生的時(shí)間或者說是異常值出現(xiàn)時(shí)區(qū)的分布.如人體的血壓過高或過低的時(shí)間預(yù)測(cè).季節(jié)性災(zāi)變預(yù)測(cè)：對(duì)發(fā)生在每年特定時(shí)區(qū)的事件和命題作預(yù)測(cè).拓?fù)漕A(yù)測(cè)：即事物整體的預(yù)測(cè)，亦稱波形預(yù)測(cè).其特點(diǎn)是對(duì)于預(yù)先給定的多組數(shù)值建立GM(1，1)模型群，根據(jù)預(yù)測(cè)結(jié)果構(gòu)造出整個(gè)波形.系統(tǒng)預(yù)測(cè)：對(duì)系統(tǒng)中眾多變量間相互協(xié)調(diào)關(guān)系的發(fā)展變化所進(jìn)行的預(yù)測(cè).2算法的基本原理2.1GM(1，1)模型：灰色模型GM(1，1)GM(1，1)的含義為1階，1個(gè)變量的灰色模型，它是在數(shù)據(jù)生成的基礎(chǔ)上建立如下灰微分方程：式中濾。)(左)為原始序列，^⑴=AGOx(0)，z⑴(k)=0.5x(1)(k)+0.5x(D(k-1).a稱為發(fā)展系數(shù)，它反映X⑴和X(。)的發(fā)展態(tài)勢(shì)；b稱為灰作用量，它的大小反映數(shù)據(jù)變化的關(guān)系.對(duì)序列乙⑴=｛乙⑴(2)，乙⑴⑶，…，z⑴(〃)｝，因?yàn)閦⑴(k)為X(1)(k)與m)(kT)的平均值，故記z⑴為MEANX⑴，艮口z⑴=MEANX⑴x(0)(k)+az⑴(k)=b的白化型為：初始值用X(1)⑴=x(0)⑴，則其解為：該式用于預(yù)測(cè)時(shí)稱為時(shí)間響應(yīng)函數(shù)，表示為累減還原：其中(a,b)可通過最小二乘求解。2.2生成數(shù)灰色模型是將隨機(jī)數(shù)經(jīng)生成后變?yōu)橛行虻纳蓴?shù)據(jù)，然后建立微分方程，尋找生成數(shù)據(jù)的規(guī)律，再將運(yùn)算結(jié)果還原的一種方法，其基礎(chǔ)是數(shù)據(jù)的生成.常用的生成方式有累加生成和累減生成.2.2.1累加生成累加生成是將原始數(shù)據(jù)通過累加以生成新的數(shù)列.記原始數(shù)列為X(0):記尤(0)的生成數(shù)列AGOX(0)為X⑴：其中x(1)(1)=x(0)(1),x⑴(k)=x⑴(k-1)+x(0)(k),k=2,3,...,n稱尤⑴為X(0)的一次累加生成，記為1-AGO定義X(0)的2—AGO為X⑵：X⑵二AGOX⑴.一般定義X(0)的r次AGO^X°):X°)=AGO^(r-1).原始數(shù)據(jù)經(jīng)累加生成后，其隨機(jī)性明顯減小，而規(guī)律性將增加.對(duì)于非負(fù)的數(shù)據(jù)序列，累加生成的是一個(gè)遞增的序列.

2.2.2累減生成累減生成是累加生成的逆運(yùn)算，它是通過將原始序列前后兩個(gè)數(shù)據(jù)相減生成新的數(shù)據(jù)序列.累減生成可將累加生成還原為非生成數(shù)列，在建模中獲得增量信息.即:IAG停⑴=尤(0)3算法的具體實(shí)現(xiàn)流程3.1算法流程圖3.2實(shí)現(xiàn)步驟(1)數(shù)據(jù)檢驗(yàn)與處理:具體算法：根據(jù)原始序列雙0)的級(jí)比a(0)(k)的大小，判斷GM(1,1)建模的可行性.a(0)(k)的定義為:a(0)(k)可揭示數(shù)列x(0)的平滑情況或指數(shù)律的符合情況，若a(0)(k)為常數(shù)，則X(0)具有白指數(shù)律(即確切的指數(shù)律)，當(dāng)然，此時(shí)無建立GM(1，1)模型的價(jià)值.GM(1,1)建模希望的是a(o)(k)的取值區(qū)間長(zhǎng)度(覆蓋的測(cè)度)接近于零而不等于零，其可容區(qū)為(0.1353,7.389)，即參數(shù)或級(jí)比的可容區(qū)是GM(1，1)建模的基本條件，但不是實(shí)用條件，要想建立滿意有效的GM想建立滿意有效的GM(11)模型，參數(shù)或級(jí)比應(yīng)落于界區(qū)內(nèi)?設(shè)n是原始數(shù)據(jù)的個(gè)數(shù)，則我們有:a的界區(qū)：ae/—2-1，若a(0)(k)e/—2-1en+1,en+1en+1,en+1kJkJa(0)(k)的界區(qū)：a(0)(k)e，則認(rèn)為X(0)是可作GM(1，1)建模的.(2)建立模型GM(1,1)具體算法：利用算法思想的介紹建立模型，則可得累加預(yù)測(cè)值和預(yù)測(cè)值:模型的白化型為：初始值用^⑴(1)f(o)(1)，則其解為：該式用于預(yù)測(cè)時(shí)稱為時(shí)間響應(yīng)函數(shù)，表示為累減還原x(0)(k+1)=x⑴(k+1)一x⑴(k)；檢驗(yàn)預(yù)測(cè)值：殘差檢驗(yàn)：設(shè)x(0)為原始數(shù)列設(shè)x(0)為通過模型(11-52)和(11-53)得到的預(yù)測(cè)數(shù)據(jù)列則稱△(k)=x(0)(k)-x(0)(k)為殘差，稱為GM(1，1)的殘差相對(duì)值，若殘差小于0.1，則稱精度達(dá)到較高要求；若殘差介于0.1與0.2之間，則稱精度達(dá)到一般要求；若殘差大于0.2,則模型不可靠。級(jí)比偏差值檢驗(yàn)：與原始數(shù)列級(jí)比一樣，我們定義模型級(jí)比推導(dǎo)可知，模型級(jí)比與k無關(guān)(因?yàn)椋鹸(0)(k)｝符合指數(shù)律)其中a為發(fā)展系數(shù).定義級(jí)比偏差一般要求Ip(k)1<20%，最好Ip(k)1<10%；由模型GM(1,1)作相應(yīng)的預(yù)測(cè)。3.3數(shù)據(jù)準(zhǔn)備與預(yù)處理輸入原始數(shù)據(jù)列。對(duì)原始數(shù)據(jù)做級(jí)比檢驗(yàn)，若不滿足要求，需要對(duì)原始數(shù)據(jù)作平移變換，直到滿足要求為止。4算法程序?qū)崿F(xiàn)4.1程序使用說明程序采用Matlab語言平移量由預(yù)測(cè)的人自行輸入，輸入值應(yīng)逐漸增大Matlab函數(shù)：GM(A)——灰色預(yù)測(cè)函數(shù)；A 初始輸入數(shù)據(jù)向量n——A的長(zhǎng)度lamda 級(jí)比向量P——需要做平移的常數(shù)B——一次累加向量alpha 為分辨率系數(shù)，一般取0.5D——B的灰導(dǎo)數(shù)a——發(fā)展系統(tǒng)b 灰作用量，且u=(a,b)Y——提取A第2列后的所有值所構(gòu)成的向量Z——B的加權(quán)均值，也稱為均值數(shù)列Z(白化背景值)T——由Y擴(kuò)展得到C——原始累加數(shù)據(jù)預(yù)測(cè)值，且滿足C(1)=A(1)C1——原始數(shù)據(jù)(A)預(yù)測(cè)值，且滿足C1(1)=A(1)C(i+1)——累加數(shù)據(jù)第i次預(yù)測(cè)值C1(j)——原始數(shù)據(jù)第j次預(yù)測(cè)值eps(i)——預(yù)測(cè)值C1(i)與原始數(shù)據(jù)A(i)的相對(duì)誤差，即殘差ruo(i) 關(guān)于lamda的級(jí)比檢驗(yàn)表達(dá)式m 輸入的值t 預(yù)測(cè)天數(shù)，且t=m-nCt——第t次累加數(shù)據(jù)預(yù)測(cè)值Ct_1——第t-1次累加數(shù)據(jù)預(yù)測(cè)值pre_Ct 第t次原始數(shù)據(jù)預(yù)測(cè)值，且pre_Ct=Ct-Ct_14.2程序源代碼functionGM(A)%數(shù)據(jù)檢驗(yàn)與處理n=length(A);%lamda為級(jí)比向量a=A(1);forj=1:n-1lamda(j)=A(j)/A(j+1);%判斷能否通過級(jí)比檢驗(yàn)ifexp(-2/(n+1))<lamda(j)&&lamda(j)<exp(2/(n+1))disp('可通過級(jí)比檢驗(yàn)');continue;elsedisp('需要做平移變換，自行決定平移系數(shù)，這里的P為任意選取的平移常數(shù)')；endwhileexp(-2/(n+1))>lamda(j)||lamda(j)>exp(2/(n+1))P二input('平移系數(shù)c=');fori=1:nA(i)=A(i)+P;endforj=1:n-1lamda(j)=A(j)/A(j+1);endendP=A(1)-a;endA;lamda;%對(duì)處理后的數(shù)據(jù)做一次累加，B為累加數(shù)列B二zeros(1,n);B(1)=A(1);fori=2:nj=i-1;B(i)=A(i)+B(j);endB;%求B的加權(quán)均值，即得均值數(shù)列Z(白化背景值)%alpha為分辨率系數(shù)，一般取0.5alpha=0.5;fori=1:n-1Z(i)二alpha*B(i+1)+(1-alpha)*B(i);end%求B的灰導(dǎo)數(shù)D，由B數(shù)列做差可得D=A;%求發(fā)展系統(tǒng)a與灰作用量b,建立方程化簡(jiǎn)可得%提取A第2列后的所有值，稱為數(shù)據(jù)向量YY=A(2:n);%將Y擴(kuò)展為數(shù)據(jù)矩陣TT=(cat(1,-Z,ones(1,n-1)))';%由最小二乘法可得u=(a,b)u=(inv(T'*T))*T'*Y';%建立GM(1,1)模型求預(yù)測(cè)值%C是原始累加模型預(yù)測(cè)值，C1是原始數(shù)據(jù)預(yù)測(cè)值(均從原始數(shù)據(jù)第二項(xiàng)開始比較)，可計(jì)算誤差C(1)=A(1);fori=1:n-1C(i+1)=(A(1)-u(2,1)/u(1,1))*exp(-u(1,1)*i)+u(2,1)/u(1,1);endC1(1)=A(1);forj=1:n-1C1(j+1)=C(j+1)-C(j);endC;C1;%殘差檢驗(yàn)disp('以下為殘差檢驗(yàn)')；fori=1:neps(i)=(A(i)-C1(i))/A(i);ifeps(i)<0.1disp('精度達(dá)到較高要求')；elseifeps(i)<0.2&&eps(i)>0.1disp('精度達(dá)到一般要求')；elseifeps(i)>0.2disp('不能通過檢驗(yàn)')；endend%級(jí)比偏差檢驗(yàn)disp('以下為級(jí)比偏差檢驗(yàn)')；fori=1:n-1ruo(i)=1-(1-0.5*u(1,1))/(1+0.5*u(1,1))*lamda(i);ifruo(i)<0.1disp('精度達(dá)到較高要求')；elseifruo(i)<0.2&&ruo(i)>0.1disp('精度達(dá)到一般要求')；elseifruo(i)>0.2disp('不能通過檢驗(yàn)')；endend%數(shù)據(jù)檢驗(yàn)到此結(jié)束%改進(jìn)預(yù)測(cè)算法 等位新信息遞補(bǔ)法，用于中長(zhǎng)期動(dòng)態(tài)預(yù)測(cè)%預(yù)測(cè)其后第t天%Ct為累加數(shù)列的預(yù)測(cè)值%pre_Ct為原始數(shù)據(jù)的預(yù)測(cè)值m=input（'請(qǐng)輸入x='）；whilem<n||m==ndisp('請(qǐng)重新輸入大于原始數(shù)據(jù)個(gè)數(shù)的值:')；m=input('x=');endP1=0;P0=A(1)-a;fort=1:m-nCt=(A(1)-u(2,1)/u(1,1))*exp(-u(1,1)*(n))+u(2,1)/u(1,1);Ct_1=(A(1)-u(2,1)/u(1,1))*exp(-u(1,1)*(n-1))+u(2,1)/u(1,1);pre_Ct=Ct-Ct_1;disp('最終平移量：’)；%最終的平移量ifP1==0P=P0elseP=P1+P0end%判斷預(yù)測(cè)最終是否應(yīng)減去平移量ifP>0disp(['第',num2str(t),天的預(yù)測(cè)值為：’]);pre_Ct=pre_Ct-PelseifP==0disp(['第',num2str(t),天的預(yù)測(cè)值為：’]);pre_Ctend%等位新信息遞補(bǔ)，即除去第一個(gè)元素，引入一個(gè)預(yù)測(cè)的新元素fori=1:n-1A(i)=A(i+1);endA(n)=pre_Ct+P;%判斷是否結(jié)束預(yù)測(cè)ift==m-nbreak;else%未結(jié)束，繼續(xù)檢驗(yàn)新得到的數(shù)據(jù)%判斷能否通過級(jí)比檢驗(yàn)a=A(1);forj=1:n-1lamda(j)=A(j)/A(j+1);ifexp(-2/(n+1))<lamda(j)&&lamda(j)<exp(2/(n+1))disp('可通過級(jí)比檢驗(yàn)');continue;elsedisp('需要做平移變換，自行決定平移系數(shù)，這里的P為任意選取的平移常數(shù)')；endwhileexp(-2/(n+1))>lamda(j)||lamda(j)>exp(2/(n+1))P1=input('平移系數(shù)c=');%整體平移數(shù)據(jù)fori=1:nA(i)=A(i)+P1;end%計(jì)算級(jí)比forj=1:n-1lamda(j)=A(j)/A(j+1);endend%此循環(huán)所得到的平移總量P1=A(1)-a;P1;endend%對(duì)處理后的數(shù)據(jù)做一次累加，B為累加數(shù)列B(1)=A(1);fori=2:nj=i-1;B(i)=A(i)+B(j);endB;%求B的加權(quán)均值，即得均值數(shù)列Z(白化背景值)%alpha為分辨率系數(shù)，一般取0.5alpha=0.5;fori=1:n-1Z(i)二alpha*B(i+1)+(1-alpha)*B(i);end%求B的灰導(dǎo)數(shù)D，由B數(shù)列做差可得D=A;%求發(fā)展系統(tǒng)a與灰作用量b,建立方程化簡(jiǎn)可得%提取A第2列后的所有值，稱為數(shù)據(jù)向量YY=A(2:n);%將Y擴(kuò)展為數(shù)據(jù)矩陣TT=(cat(1,-Z,ones(1,n-1)))';%由最小二乘法可得u=(a,b)u=(inv(T'*T))*T'*Y';end4.3程序運(yùn)行4.3.1程序運(yùn)行及運(yùn)行環(huán)境說明運(yùn)行環(huán)境：windows操作系統(tǒng)運(yùn)行軟件：MATLAB適用范圍：只能進(jìn)行短期或中期的預(yù)測(cè)4.3.2輸入數(shù)據(jù)輸入實(shí)際數(shù)據(jù)向量：a=[2.28,2.98,3.39,4.24,6.86,8.64,1