GM(1-n)(灰色模型代码)

%灰色预测模型GM（1，n）模型的matlab源代码，包括预测模型的建立，以及模型的精度检验指标c,p的计算%假设预测3步，N=3%如在命令窗口键入：%gm=ycgm1n([1.6,1.7,2,1.8,1.9],[2,2.4,3,3.2,3.1],[3,3.1,3.2,3.5,2.8],3)functionGM=ycgm1n(data1,data2,data3,N)%data1:纵摇，data2:升沉，data3:波浪T=length(data1);PYX1=data1;PYX2=data2;PYX3=data3;%进行数据预处理，这里用初值化X0_1=PYX1./PYX1(1);X0_2=PYX2./PYX2(1);X0_3=PYX3./PYX3(1);%用AGO生成一阶累加生成模块X1_1(1)=X0_1(1);X1_2(1)=X0_2(1);X1_3(1)=X0_3(1);fori=2:TX1_1(i)=X1_1(i-1)+X0_1(i);X1_2(i)=X1_2(i-1)+X0_2(i);X1_3(i)=X1_3(i-1)+X0_3(i);end%构造累加矩阵Bfori=1:T-1M1(i)=(0.5*(X1_1(i)+X1_1(i+1)));M2(i)=(0.5*(X1_2(i)+X1_2(i+1)));M3(i)=(0.5*(X1_3(i)+X1_3(i+1)));endB1=zeros(T-1,3);fori=1:(T-1)B1(i,1)=-M1(i);%-(X1_1(i)+X1_1(i+1)))/2;B1(i,2)=X1_2(i+1);B1(i,3)=X1_3(i+1);endB2=zeros(T-1,2);fori=1:(T-1)B2(i,1)=-M2(i);%-(X1_2(i)+X1_2(i+1)))/2;B2(i,2)=X1_3(i+1);endB3=zeros(T-1,2);fori=1:(T-1)B3(i,1)=-M3(i);%-(X1_3(i)+X1_3(i+1)))/2;B3(i,2)=1;endsaveB1B1;saveB2B2;saveB3B3;%构造常数项向量Yfori=2:TY1(i-1)=X0_1(i);Y2(i-1)=X0_2(i);Y3(i-1)=X0_3(i);endHCS1=inv(B1'*B1)*B1'*Y1';%用最小二乘法求灰参数HCS1H1=HCS1';%H1=[a,b2,b3]HCS2=inv(B2'*B2)*B2'*Y2';%用最小二乘法求灰参数HCS2H2=HCS2';%H2=[a,b3]HCS3=inv(B3'*B3)*B3'*Y3';%用最小二乘法求灰参数HCS3H3=HCS3';%H3=[b,a]%计算出X3的累加序列fori=1:T+NYCX13(i)=(X0_3(1)-H3(2)/H3(1))*exp(-1*H3(1)*(i-1))+H3(2)/H3(1);endfori=2:T+N%K(i)=XR1(i)-XR1(i-1);YCX0_3(i)=YCX13(i)-YCX13(i-1);endYCX0_3(1)=X0_3(1);%对参数作alpha，beta变换H2=H2./(1+0.5*H2(1));%还原计算出X2的预测值YCX0_2(1)=X0_2(1);fori=2:TYCX0_2(i)=H2(2).*X1_3(i)-H2(1).*X1_2(i-1);endYCX12(T)=X1_2(T);fori=T+1:T+NYCX0_2(i)=H2(2).*YCX13(i)-H2(1).*YCX12(i-1);YCX12(i)=YCX0_2(i)+YCX12(i-1);end%对参数作alpha，beta变换H1=H1./(1+0.5*H1(1));%还原计算出X1的预测值YCX0_1(1)=X0_1(1);fori=2:TYCX0_1(i)=H1(2).*X1_2(i)+H1(3).*X1_3(i)-H1(1).*X1_1(i-1);endYCX11(T)=X1_1(T);fori=T+1:T+NYCX0_1(i)=H1(2).*YCX12(i)+H1(3).*YCX13(i)-H1(1).*YCX11(i-1);YCX11(i)=YCX0_1(i)+YCX11(i-1);end%数据还原GM=YCX0_1;%.*PYX1(1);saveGMGM;e0(1,T-1)=zeros;fori=1:T-1%求X1到X5的残差值e0e0(i)=(X0_1(i+1)-YCX0_1(i+1))/X0_1(i+1);%1-YCX0_1(i+1)/X0_1(i+1);endsavee0e0;e0_average=sum(abs(e0))/length(e0)p=1-e0_average;X_average=mean(X0_1)%求原始数据x0均值s1=std(PYX1)%求原始数据的标准差s2=std(e0)c=s2/s1%计算方差比c,c<0.35为好end