matlab 灰色系统预测 GM(1,1) 数学建模

佛自业障

本文代码主要是基于邓聚龙教授在20实际80年代提出的灰色系统理论。
GM0.m
5 f7 u  o; i0 \8 T! n; H%该函数为GM（1，1）模型返回还原值
6 y5 g, W! U6 o1 Jfunction f=GM0(x0,t)  %数据数列
[M,N]=size(x0);        %算出数据数列的大小
x1(1)=x0(1);           %累加生成数列
for i=2:N;     
2 B7 q) x( E( e* f) M    x1(i)=x1(i-1)+x0(i);
end
x2=[];              %累加生成数列均值生成数列
for j=1N-1);     
$ v7 p( e1 v0 E% \    x2(j)=(x1(j)+x1(j+1))/2;
end
x=x0;              %数据数列镜像
; n. Z& O! l4 v1 [2 \. ~7 Ax(1)=[];           %删除第一个数据
Y=x';              %数据列向量
6 S% l, ]* p+ K! O2 N* rglobal a;
global b;
B(:,1)=-x2';
B(:,2)=1;
A=inv(B'*B)*B'*Y;   %求参量a，b组成的参数向量
) y4 A9 u7 n) i5 l# Wa=A(1,1);           %求参数a
# T/ k2 L( p7 f2 W% a7 V; w7 ob=A(2,1);           %求参数b  
f=(1-exp(a))*(x0(1)-b/a)*exp(-a*(t-1));
6 v$ I+ r/ b4 s$ [f

& [: O4 g. L$ I- Y/ W6 ^/ V1 EGM1.m:
. o9 M8 V3 n' {: o; O%该函数为GM（1，1）模型中数据数列进行光滑比检验
- j. e* v9 w, b- p4 v/ b6 |3 qfunction f=GM1(x0)   %数据数列
) C' K, @; g* b: ~# T& P; g6 ^N=max(size(x0));       %算出数据数列的大小  
x1=cumsum(x0);         %累加生成数列
global J;
1 p) A" D5 i* B% pglobal J1;
global J2;
x0(1)=[];
7 y  a, `7 q5 ~: h8 J! Zx1(N)=[];
$ i0 z0 U* \7 sglobal r;
r=x0./x1;  
for j=2N-1);           %判断数据数列是否满足准光滑条件1   
   if(r(j)>=0.5||r(j)<0)         
       J1=0;         
       break;     
   else
2 d' r& h! E0 y6 j& n       J1=1;     
2 L0 S0 U3 V; w) M& H) O   end
end
for l=1N-2);           %判断数据数列是否满足准光滑条件2     
    if((r(l+1)/r(l))>=1)         
        J2=0;         
        break;     
    else
) n  U1 o  X4 e% `" }        J2=1;     
  h! h  S" O, v3 u& v4 Q    end
end
' D) S5 d1 }6 LJ=J1+J2;  
if(J==2)                 %判断数据数列是否为准光滑数列     
1 [& M+ M# l7 @4 v& H+ c; C    disp('数据为准光滑数列')
1 w4 ~& d! A, C0 q& M" H5 Kelse
4 A2 J) b/ I4 a6 [( x) T    disp('数据不是准光滑数列')
end

' x/ d! f- K0 z. I: cGM2.m
%该函数为GM（1，1）模型还原值参数计算
function f=GM2(x0)  %数据数列
[M,N]=size(x0);      %算出数据数列的大小  
x1(1)=x0(1);         %累加生成数列
for i=2:N;      
% K  N! X4 N+ g8 }) D    x1(i)=x1(i-1)+x0(i);
# A$ [7 Q  m& s( |1 vend
x2=[];              %累加生成数列均值生成数列
6 i- ]! c8 U6 j5 cfor j=1N-1);      
    x2(j)=(x1(j)+x1(j+1))/2;
0 k, N1 S. ~  L) ^9 L3 S7 y; }end
  x=x0;              %数据数列镜像      
3 Q* u- [; m% k! k1 P8 {  x(1)=[];           %删除第一个数据
  Y=x';              %数据列向量
  global a;
  global b;
$ ?+ H7 |. |7 @' K0 v1 F# f. v- z  B(:,1)=-x2';
  B(:,2)=1;  
  A=inv(B'*B)*B'*Y;   %求参量a，b组成的参数向量
  a=A(1,1);           %求参数a
( C8 f% X- B* w, F. u! f% V8 o  disp('参数a为：')
  a
0 ~) K4 R; e, u7 G5 [* E) x  b=A(2,1);           %求参数b
  disp('参数b为：')
8 f$ L4 c7 A6 P  y2 _  b

* O. u; j& _" A+ R# @0 SGM3.m
%该程序实现G(1,1)模型的精度检验
+ m, C/ S, J9 @9 ?, }%包括平均相对误差，绝对关联度，均方差比值，小误差概率检验
function f=GM3(x0)
N=max(size(x0));
x=GM0(x0,1:N);  %利用已有程序GM2得出数据列模型估计值
% K" s2 ]  ?" m: n1 r  ?6 \. w0 P. ^x(1)=x0(1);     %更正第一个估计值
& c* I3 f8 K6 Udisp('模型模拟估计值为')
/ ^9 P3 T0 O' P) Z, \6 ?- Ux
* G! A/ @9 s7 B4 W0 mA=x-x0;         %计算绝对残差序列
" \+ R) Q3 |( b9 v7 M. e6 Ldisp('模型估计值绝对残差序列为：')
A
G=abs(A);
Amin=min(G);    %计算最小绝对值绝对残差
$ F, n( m9 E! W- m8 `3 c# t. ZAmax=max(G);    %计算最大绝对值绝对残差
B=A./x0;        %计算相对误差序列
- Q3 l* X8 ?# L- r7 bdisp('模型估计值相对误差序列为：')
6 k1 C$ B9 t6 kB  
8 E6 L* }% i9 c3 G' J2 K7 |P=sum(abs(B))/N;     %计算平均相对误差
disp('模型估计值平均相对误差为：')
P  
* L- u; _- N3 e2 }. u  J( Xfor i=1:1:N       %通过循环计算关联系数序列     
    D(i)=(Amin+0.5*Amax)/(G(i)+0.5*Amax);
$ O1 W) l( r1 G# p7 z5 gend
R=sum(D)/N;  
disp('关联度为：')
; w0 J& x1 f$ H( v, g; b% vR  
) w. _& _8 \2 S% w+ e, _; zx_=sum(x0)/N;    %计算数据的均值  
' e+ ], a0 [7 `* GS1=(sum((x0-x_).^2)/(N-1))^0.5;   %计算数据序列方均差
; ^2 X8 O4 s: A) x# @9 b1 BA_=sum(A)/N;     %计算残差平均值  
. g! V) Y0 }* kS2=(sum((A-A_).^2)/(N-1))^0.5;    %计算残差序列方均差
C=S2/S1;         %计算方均差比值
1 ?- J; I+ }$ @' W. Z  h- Vdisp('均方差比值为：')
C  
+ V3 a7 o# V: C/ y. HS0=0.6745*S1;
E=A-A_;
F=find(E<S0);
M=max(size(F)); %计算小残差个数
p=M/N;          %计算小误差概率
disp('小误差概率为：')
6 ^2 S7 S4 k" P0 J2 Qp
) t! E" c; o9 h  l" k. I8 ~2 D
4 ~$ g% y  |) z% U1 d3 n/ r( t