求助，10个数据的灰色模型代码

Mlearner

问题同上，需要可执行，在网上下了几个，各种错误，难以执行，希望程序中有详细的说明。

luoshichao123

这个程序自己编啊，原理很简单的

Mlearner

luoshichao123 发表于 2012-5-16 07:32
这个程序自己编啊，原理很简单的

网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x0
   4 k# X8 |) J* O\" M
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点
     z; U& ^' k' V6 m* }/ I
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);
   9 \# l5 k2 {6 U$ u% N2 W
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);
   $ G3 J: s  D8 l
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);1 r, _7 o# j9 N% h9 c
   
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);
   9 v/ K$ y8 A  U. N7 b8 @8 ]
7. for i=1:length(x0)+ u! N+ j! C* i: M+ N. T5 q4 y
   
8.     for j=1:i
   ) l$ L* X- B\" x! L: R
9.         x1(i)=x1(i)+x0(j);
   / Q% k- B1 W0 D5 j
10.     end5 z' |) D$ |2 U% q6 ?
    
11. end
    # v/ l! E. `: f$ R  u# r
12. for i=1:length(x0)-1
    9 V* f; @/ u5 B9 A
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));' V8 D5 @, w  K9 m) b7 j
    
14.     B(i,2)=1;$ _+ i; A: H( j9 y+ [# A
    
15.     yn(i)=x0(i+1);
    ! x' t0 t/ N3 ?0 e6 P/ N+ G
16. end2 ]6 h9 b7 B* H6 C4 G
    
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计$ d0 {+ k8 r& W
    
18. for k=1:length(x0)+T5 x- `& ^0 p; Q9 T4 p8 \
    
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);# b9 A0 q9 ?' H2 \
    
20. end
    ( j! S( `6 ]4 g8 K  A: S
21. Hatx0(1)=Hatx1(1);
    9 h2 z4 M  R, K# A7 K9 r2 E
22. for k=2:length(x0)+T
    % `! e# x' W8 d$ H8 ^& s
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值
    - {# ^* |3 b$ ~- Q* s, U0 e) b
24. end$ ~! ?% d0 T\" T, L
    
25. for i=1:length(x0) %开始模型检验% m: c9 y4 P) o\" U' j. z4 j7 J
    
26.     epsilon(i)=x0(i)-Hatx0(i);
    : x+ K( Q  Z7 P6 v$ o: T
27.     omega(i)=(epsilon(i)/x0(i))*100;1 `1 E+ N9 b. P. q9 \- g0 P% y3 D
    
28. end
    7 f; D* M+ I, i9 _- ^
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据: V. S% A$ g+ [% Y
    
30. c=std(epsilon)/std(x0);p=0;' ^. L# r4 v5 _% Y
    
31. for i=1:length(x0)
    7 o8 P: ?' q; b. j
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)
    4 P' I7 K6 t+ f# ]& z; G4 k
33.         p=p+1;& _8 ^: l: a, g( k  K
    
34.     end: i% @; j\" {- j* q; Q7 ?0 e
    
35. end2 F. j# b1 T4 \+ x' v7 h) E
    
36. p=p/length(x0)% r\" Y& K# a\" J, d+ O$ F9 K& R: \
    
37. if p>0.95 & c<0.35
    8 Z9 ]# V3 A/ G7 Y- J
38.     disp('The model is good,and the forecast is:'),+ @7 s' @8 s! P7 G
    
39.     disp(Hatx0(length(x0)+T))- I# S) E' M; ?& j
    
40. elseif p>0.85 & c<0.5- ~4 O' x5 H. U1 x8 f$ W
    
41.     disp('The model is eligibility,and the forecast is:'),\" y& A1 m/ ~  }0 s$ |$ V
    
42.     disp(Hatx0(length(x0)+T))
    8 y' U7 L5 F6 A) g. L8 g
43. elseif p>0.7 & c>0.65
    8 ~; `: V8 |+ g! q1 w, \
44.     disp('The model is not good,and the forecast is:'),+ L\" D% e$ D& b\" }4 h
    
45.     disp(Hatx0(length(x0)+T)); r& _( r; @: Q0 S# w: |. i
    
46. else p<=0.7 & c>0.65
    ! e4 h5 O\" c) s+ E1 r/ v+ y- J
47.     disp('The model is bad and try again')
    & f7 G; L  w, U6 @5 J
48. end
    , z& l' ~5 Q3 |: u
49. for i=1:length(x0)
    ; t' ]; j0 f1 l\" _. {( m5 |
50.     Hatx00(i)=Hatx0(i);4 B0 G- y. x/ \0 Q- D
    
51. end
    # k/ ~( ^( c* s0 |) [7 `
52. z=1:length(x0);
    4 W/ ]2 N1 |6 c7 |% w  ^/ y: r
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察
    ! a% N) t. g0 U
54. text(2,x0(2),'History data: real line')0 k% R* x0 a$ g3 K\" f
    
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line'); s  Y: h+ L. o- {( K
    
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2]
   # m+ z& s+ |# J  n8 a, c  b& l
2. Warning: Input arguments must be scalar.\" B) m\" g. v7 {2 a/ Y6 M7 a# v; S2 x$ k
   
3. > In fungry1 at 48 K5 A2 F0 l' L2 \/ e\" [3 i
   
4. Warning: Input arguments must be scalar.9 W! g0 V* x+ D- L4 {  M( E
   
5. > In fungry1 at 5! p, w! e) G+ f, ?8 h\" _9 U
   
6. Warning: Matrix is singular to working precision.
   5 I# T& r, U+ N/ e' d0 d7 I4 D
7. > In fungry1 at 17
   * T+ C5 ~# S. A1 `
8. 
   5 v; g. p& Y# V% e1 o0 J& o
9. HatA =9 n& H/ b. s# ~2 K7 E
   
10. 
    9 s! Q0 R  `/ r% Z! c) G: D
11.      0; p7 p8 A6 K. L) w* z! [
    
12.      0) J- Y$ u4 E7 [' y) r
    
13. , x6 L  K! `8 ?. }
    
14.   |+ L0 r4 C$ h- p/ g3 }
    
15. p =\" _! }( d. M3 \7 V' r
    
16. 7 z2 L0 N5 |/ }8 U4 U
    
17.      0
    9 ~- ^5 N1 W4 J+ s# b2 r
18.   ?& M6 M- K8 p- f0 t
    
19. 6 o/ R( u7 }; ^2 W% A
    
20. ans =
    ( D- w; j5 Q! _. L6 V3 F6 K+ i& a9 n7 F! p8 O
21. 7 _) l2 ?& R- u
    
22.      0) J6 x8 ^) `# w/ J7 R! L4 @
    
23. 
    - l1 }9 ?3 P8 B  b
24. The model is bad and try again) G6 r3 l% g8 d3 ]3 W
    
25. ??? Index exceeds matrix dimensions.
    4 F- n7 l1 ^* ^\" [% n* ^3 j
26. 
      C8 o. R) U7 N& q\" m2 W
27. Error in ==> fungry1 at 545 t, }4 y- d% F& e
    
28. text(2,x0(2),'History data: real line')
    5 h3 W0 v+ g8 y* J  a9 y
29. 5 `3 D& H4 o2 h: h; y6 b
    
30. >>

复制代码

是怎么回事。

Mlearner

求高手解答。

zqyzixin

牛啊,想不到的强帖