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

Mlearner

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

zqyzixin

牛啊,想不到的强帖

Mlearner

求高手解答。

Mlearner

luoshichao123 发表于 2012-5-16 07:32
$ l& I% T+ J# f7 q% x# C- }这个程序自己编啊，原理很简单的

网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x0
   ; S- p3 {; m/ x8 J
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点7 l2 S, g7 U0 O0 E' c
   
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);
   # \( A& j0 |# K  t, }; M/ P7 [
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);
     F3 L\" P) C3 F+ Z
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);3 M- O4 m  n8 L' ~! q) p
   
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);$ M/ @+ c+ t4 c0 U: \) g
   
7. for i=1:length(x0)% T; K; p8 D1 [/ S
   
8.     for j=1:i$ c' t6 ^8 A# ]! {6 L+ z. M$ ]
   
9.         x1(i)=x1(i)+x0(j);  b: ^# L6 |9 ?! q- k
   
10.     end$ ]! d( R; t2 h& n% T6 N
    
11. end
    : Z4 g/ h. M7 M, M2 C7 e0 k
12. for i=1:length(x0)-1' ?' |7 A5 G, c; f! T( A2 E
    
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));
    1 |2 V: q6 w2 O$ n
14.     B(i,2)=1;$ V. H3 C) J/ Z5 d
    
15.     yn(i)=x0(i+1);\" l  ^: v; X: L) I! s5 F; B8 a9 N& t
    
16. end  x' b$ U4 ~( o! {) b1 o1 K. x
    
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计. V\" n6 z: W2 i1 o! F2 e9 f
    
18. for k=1:length(x0)+T4 T. f\" Y1 T# f. b* E$ s
    
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);
    * X& _4 b4 q3 M) k( }( _1 f
20. end, J7 ?* F( b# }' ?
    
21. Hatx0(1)=Hatx1(1);
    8 [) p( V8 T, w$ u$ I2 T
22. for k=2:length(x0)+T/ L& {# D# {! y/ ]6 g
    
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值4 k0 ]! K( m& n# Y0 n8 M  V# z
    
24. end8 X) `1 l9 [& Z; c
    
25. for i=1:length(x0) %开始模型检验6 N6 L. D7 m! T
    
26.     epsilon(i)=x0(i)-Hatx0(i);
    2 n: \3 m\" v, S! d) _3 _
27.     omega(i)=(epsilon(i)/x0(i))*100;
    $ W) j' D# {& H% J, r. V
28. end4 m9 N7 ~- h/ e7 X1 V
    
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据
    : f8 l+ L; J( e+ h3 C
30. c=std(epsilon)/std(x0);p=0;
    $ ~2 A8 z7 u8 f* |% I6 g\" Q
31. for i=1:length(x0)
    # ^( J. q8 L/ K1 p
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)
    3 R9 O\" g/ y3 p5 W/ S
33.         p=p+1;
    \" o- S/ o/ k) u) w. G
34.     end5 g9 x  ^. e' b% `. J
    
35. end# W! V% Q0 }; O6 s+ [- j$ ^+ c
    
36. p=p/length(x0)( K9 E8 B- D+ y1 R5 ^+ M  D* A
    
37. if p>0.95 & c<0.35
    9 h1 @% V- t# x2 ^+ m8 D. ]
38.     disp('The model is good,and the forecast is:'),$ A0 s  R0 s. ^' }# R8 |' ~1 U
    
39.     disp(Hatx0(length(x0)+T))
    \" d- D+ M1 s\" I6 F, X
40. elseif p>0.85 & c<0.5
    ( W4 U; e7 ~7 L0 d7 ~6 _6 r4 ^4 k
41.     disp('The model is eligibility,and the forecast is:'),
    4 U, j! l; i, W5 \
42.     disp(Hatx0(length(x0)+T))7 S) x6 \) o\" r3 _, S
    
43. elseif p>0.7 & c>0.652 k5 r8 _: {; m
    
44.     disp('The model is not good,and the forecast is:'),1 s% \1 W9 A6 r5 ?: P1 ]- \+ {, g
    
45.     disp(Hatx0(length(x0)+T))/ ]9 Z8 ^( {2 P& L5 Z4 ~2 }
    
46. else p<=0.7 & c>0.65
    , ?# D1 D' O- @4 l1 e7 q
47.     disp('The model is bad and try again')
    6 }2 x8 [( J8 M' p, @
48. end
    ; j+ q2 e- r) L% ~
49. for i=1:length(x0)
    , u! ^- J# j% h) C6 i+ ]! H. t
50.     Hatx00(i)=Hatx0(i);
    % X( q6 w- b. A; m6 ], K
51. end
    7 j9 o5 @  h4 C! R, W& I
52. z=1:length(x0);3 {\" p3 Y$ _8 [% \
    
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察6 K9 `+ j! b: B/ @
    
54. text(2,x0(2),'History data: real line')# M! l- S4 S7 n9 q2 @
    
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')
      D\" z9 [& t5 p( U
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2]/ f' b8 U2 F: m# Q
   
2. Warning: Input arguments must be scalar.
   % i+ M  u\" K  M3 C0 u) G9 Q
3. > In fungry1 at 4
   ' M9 F9 q. |* T* g
4. Warning: Input arguments must be scalar.$ I' u6 w  K- D- t+ `9 d
   
5. > In fungry1 at 5
   6 t, P- d4 ~. ^# Z1 \  B
6. Warning: Matrix is singular to working precision.
   6 V  ?! Q: [2 ^5 |# Y\" v2 G2 l
7. > In fungry1 at 17
   9 }. N; d0 y* B
8. 
   / c& N2 R7 S, b' ~
9. HatA =
   1 G' W) z- n7 ]\" w% g3 {5 z
10. % i9 p* i) v4 L
    
11.      0- {. q8 M4 y4 i. v. M
    
12.      0
    : Y1 y7 I  j  ~: ]3 E/ ?
13. 
    $ s& j, F# @6 r9 |\" K+ j6 {\" G0 \
14. ( T# Y2 E2 ~% `7 t
    
15. p =% j3 X) R' j; ]7 W
    
16. 
    * T\" b* F0 O3 _. j8 h& e  N+ C
17.      0$ `- c# x/ c/ `( |' L! A4 M
    
18. 
    ; u6 @* Z( T. r
19. 1 f. C' g) T6 b1 Z\" G- y7 \: M2 q
    
20. ans =
    2 p4 j$ B6 d3 ^\" y  ~- }3 y/ y
21. 
    7 b* W; ^2 X+ i1 r* o
22.      0
    # [% ]& \! L5 S  |
23. / `0 d- i' i\" O+ J( k! Y+ [
    
24. The model is bad and try again8 R2 k8 M0 Z0 `9 k9 x6 \
    
25. ??? Index exceeds matrix dimensions.$ Z$ _  I* `+ g. Z3 O\" `% A
    
26. 8 U/ Y' N\" ~5 N9 e& Q4 Y& x
    
27. Error in ==> fungry1 at 543 b: o( _1 m; I
    
28. text(2,x0(2),'History data: real line')  L/ R& ?5 d4 @, K
    
29. 
    , y5 p' ?2 j9 q4 z3 _. G  S
30. >>

复制代码

是怎么回事。

luoshichao123

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