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

Mlearner

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

luoshichao123

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

Mlearner

luoshichao123 发表于 2012-5-16 07:32
% t' D/ z4 \  H; S0 V' s$ a这个程序自己编啊，原理很简单的

4 l" l& |) J9 [0 H! h- k网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x08 N- U* W9 e+ t: U( _; f
   
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点4 {# s5 T0 m1 d\" N. \' C  K9 G
   
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);. B\" V6 e9 \# U- p
   
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);/ J3 I' X! Z# A3 f
   
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);
   ! s) ~) b+ N2 E. g4 }1 C
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);) n8 z2 n# o% @6 {/ L# d
   
7. for i=1:length(x0)  E! a+ y) A( \; _( m
   
8.     for j=1:i
   & Z$ `% e# \9 Y4 g+ R- C
9.         x1(i)=x1(i)+x0(j);# _\" x7 H! v, h' A7 b9 h
   
10.     end
    % d+ u1 D\" _; E4 ^1 |+ p
11. end
      ?( F7 M0 ~1 U: \  `2 y\" C8 Q
12. for i=1:length(x0)-1
    ! ?* h* R* |! M
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));
    $ {8 B+ F, E; ?1 D$ L, W8 ~
14.     B(i,2)=1;
    * h+ o. K/ C& V+ W9 {
15.     yn(i)=x0(i+1);9 {) J. O( L& I$ V! S6 I
    
16. end6 I; j1 K0 \$ D
    
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计
    + @2 D5 r7 C( w; Q  C) c
18. for k=1:length(x0)+T- ]5 S2 k& u\" J5 @# b4 Z! N
    
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);$ M5 \- L5 z, P# ?
    
20. end
    + L' s: Q) o! F$ X+ `( R
21. Hatx0(1)=Hatx1(1);
    ; j4 q1 _9 O/ `/ J8 l2 r& K3 e3 z9 `
22. for k=2:length(x0)+T7 M8 s% Q' j( i$ o1 `& V+ V
    
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值0 k( G9 F  }' t9 H
    
24. end
    , C& `7 P+ P( q! ?) w7 H( X& W
25. for i=1:length(x0) %开始模型检验
    \" @3 Y# G$ A; \  r& r, w
26.     epsilon(i)=x0(i)-Hatx0(i);
    : y1 O6 e! t+ v5 ?\" z- Z0 [( g
27.     omega(i)=(epsilon(i)/x0(i))*100;
    , e- f/ Q  u8 w; H3 X
28. end) q# ?0 q4 O. `. n4 D' k3 |& C
    
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据
    ! N# R2 k  _6 V/ M; `
30. c=std(epsilon)/std(x0);p=0;9 \  r\" z8 l0 {9 b( s, ~
    
31. for i=1:length(x0)
    ' W5 N- y/ l& a; ]! n
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)
    7 m9 p\" |( x- A) J( P2 F
33.         p=p+1;
    3 t, t/ S5 H  C4 e' h9 P+ k\" G
34.     end. I0 i\" z; h7 B\" q# \* W+ V
    
35. end
    4 H, b3 {9 \. |: k% Z7 k
36. p=p/length(x0)
    * S6 ?8 b6 H0 ?, v
37. if p>0.95 & c<0.35
    - B+ K  ^( j% z* k$ Z7 h5 O
38.     disp('The model is good,and the forecast is:'),7 E! v4 w2 H- h) Z
    
39.     disp(Hatx0(length(x0)+T))
    9 c7 Z+ Y$ B& W
40. elseif p>0.85 & c<0.5* W5 w/ u: u) E* S/ A7 q3 w\" |- E
    
41.     disp('The model is eligibility,and the forecast is:'),7 K4 z9 `\" {1 R3 k$ [% T\" L5 o
    
42.     disp(Hatx0(length(x0)+T))
    ' o8 Z7 t; h+ w0 k8 V
43. elseif p>0.7 & c>0.656 ~9 z2 _7 C2 t; O7 y) `* Y
    
44.     disp('The model is not good,and the forecast is:'),
    7 m. U\" U' k9 T) g, Y) t6 U+ i0 A
45.     disp(Hatx0(length(x0)+T))
    7 H, o0 m+ o8 C) X* l' ]\" k
46. else p<=0.7 & c>0.65
    ! ?% T8 ?4 D  N+ d( ]\" Y/ b6 O\" e
47.     disp('The model is bad and try again')
    ! |0 S; m( W; u- u) d5 V
48. end9 _5 s+ G: X/ Z. V. i
    
49. for i=1:length(x0)
    ( P% a) l) u3 u6 B
50.     Hatx00(i)=Hatx0(i);
    2 b. C& ^' ~( H0 ?! c$ g
51. end
    , J3 D: c* I# v7 C3 M7 z* D
52. z=1:length(x0);6 G: I7 q- g3 }! f
    
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察
    , ^: C\" p+ w( F, i) G. j
54. text(2,x0(2),'History data: real line')* C* \9 l  x\" ?3 g9 p& H
    
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')
    6 s- v& T& j: M, I' F  Y9 J6 F
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2]# l8 d) x5 K9 x) _3 H& Q/ u( N: s
   
2. Warning: Input arguments must be scalar.
   : {9 w! A\" }/ Z! T6 A5 S
3. > In fungry1 at 4$ Y6 K; `, L) t* Y$ G4 ^
   
4. Warning: Input arguments must be scalar.
   . G0 R) D( A$ y) \& r, O1 |+ O
5. > In fungry1 at 5
   % C* M\" j4 k7 @4 B% c: j
6. Warning: Matrix is singular to working precision.! v& c; I0 _% o# q8 c: x
   
7. > In fungry1 at 171 i- S) A* M+ F3 U9 A
   
8. 
   ' |0 o9 F4 f0 S
9. HatA =\" Q6 ?* L9 J$ n  J
   
10. 
    ' e: v6 Y3 e' ~( H$ {2 C
11.      0
    1 T\" c6 P/ o0 P# E9 ?
12.      0
    * U0 d/ I7 H5 _! v
13. 
    ) k+ ^1 A& T4 n: ~. t
14. 
    3 @1 I4 m  @' K- b
15. p =! l/ o/ M4 {# N& {( E0 A
    
16. 
    : s% k; P/ g+ w6 {0 Z/ s, r' U
17.      02 o) ~' H\" K) p  g2 J. c* |
    
18. 2 Z% P2 `4 S- Q7 l. b# `
    
19. / f6 o6 I$ T& }0 e* E) ?/ O
    
20. ans =
    * g6 L7 }( C\" c$ u+ a
21. ' e/ s1 {: V\" R2 b3 L1 I
    
22.      05 f- t: ]$ `' p* J
    
23. % M$ f/ d& ^5 m7 t3 ?( I8 V# A
    
24. The model is bad and try again: b3 h' V( Z0 p, ?4 k5 T2 E. t$ U
    
25. ??? Index exceeds matrix dimensions.1 o( T  b5 F2 A# [7 @8 x
    
26. 
    1 a* j( c& K. ^9 o9 E& \1 p: H
27. Error in ==> fungry1 at 54
    ) U% z, I) I( X4 X, d
28. text(2,x0(2),'History data: real line')
    . J) x% P\" w# B2 E. R% e9 U
29. 
    \" ^% R$ U0 z/ n. B! ^! E/ U
30. >>

复制代码

是怎么回事。

Mlearner

求高手解答。

zqyzixin

牛啊,想不到的强帖