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

Mlearner

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

luoshichao123

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

Mlearner

luoshichao123 发表于 2012-5-16 07:32
5 f5 m$ s. z0 N5 Y9 V0 ~这个程序自己编啊，原理很简单的

网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x0: \% o) L8 G: q: g6 J
   
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点, A2 a) o' |9 B! k# q1 Y
   
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);
   5 b6 W& ?$ j2 {, t* ~
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);
   ( A\" w% m* J\" M# M0 u# A
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);$ ~* V, d' a0 h
   
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);
   $ Y2 e+ @: D8 ]2 C( t/ K
7. for i=1:length(x0)
   3 K$ v8 }# i, d# M% Z
8.     for j=1:i$ ?5 j- E- G9 j\" G% l( H% G% J
   
9.         x1(i)=x1(i)+x0(j);( x- e- h, W) f$ t; U+ [# `
   
10.     end
    \" W$ A2 z2 B4 H+ k) g
11. end
    - d8 p$ i2 w: m( R! y; ~7 X
12. for i=1:length(x0)-15 y3 W$ F9 b/ C8 z
    
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));
    , W2 _8 c# u6 e1 p
14.     B(i,2)=1;
    , g& p; ~( Z! S, ~
15.     yn(i)=x0(i+1);& K. ?- ~) j  {2 Y! J0 Y& ^
    
16. end6 Q4 \; }7 A$ g( H( c* U0 A
    
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计
    . k4 g\" |* j% t, `2 z  \, \+ S0 j
18. for k=1:length(x0)+T1 P7 n2 ]# L- ~) ?6 G$ o. L
    
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);2 n# r* y/ r. d0 q
    
20. end
    3 R$ ?+ |, w- o4 m% V
21. Hatx0(1)=Hatx1(1);% ~- S\" o2 p/ L+ D; n7 N9 r9 h3 T
    
22. for k=2:length(x0)+T
    - l8 c  I* O\" w: w* `% f; c$ c. Q
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值: |0 O; i2 g1 P7 y+ a- }
    
24. end
    8 ^  ]! O+ u0 I
25. for i=1:length(x0) %开始模型检验
    ' z; o# p- J1 N2 v- c) f2 u
26.     epsilon(i)=x0(i)-Hatx0(i);
    & T, e( r, g, m% b& ?4 u
27.     omega(i)=(epsilon(i)/x0(i))*100;7 m: ]' S8 w  R9 K5 b! D
    
28. end6 M: _5 {7 p2 A8 ~! v4 m1 q
    
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据
    - G\" Z5 q4 f, W3 Y* |
30. c=std(epsilon)/std(x0);p=0;
    - b0 J5 @+ ^' W6 `2 s
31. for i=1:length(x0)
    ( b& K: \' X: n6 c7 C8 C) X/ ~
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)% c' Q+ m9 ^5 L. J6 T3 M
    
33.         p=p+1;
    * M7 _7 x1 z+ O# E! e
34.     end6 r: N. r1 G7 y& x
    
35. end. |/ v, y) A7 F
    
36. p=p/length(x0)/ V+ o/ c' j) G
    
37. if p>0.95 & c<0.35
      x% e5 K6 B: J& i; ~$ X( k/ b
38.     disp('The model is good,and the forecast is:'),
    , r! P; R: `$ b6 ~3 x' }
39.     disp(Hatx0(length(x0)+T))\" c8 a- J/ C, y9 M3 u  w% H
    
40. elseif p>0.85 & c<0.5
    4 H  A9 d, O. R* H  q
41.     disp('The model is eligibility,and the forecast is:'),. x4 z: [% p0 q. @
    
42.     disp(Hatx0(length(x0)+T))\" a) D6 c( X0 r  P, M
    
43. elseif p>0.7 & c>0.659 S- C1 ]: H' E% e; f
    
44.     disp('The model is not good,and the forecast is:'),# a: }5 I- e  S+ k) G
    
45.     disp(Hatx0(length(x0)+T))7 j\" c8 v. ^* p# `0 C$ o
    
46. else p<=0.7 & c>0.65
    3 `! s' X4 P$ x7 H+ J0 P\" m
47.     disp('The model is bad and try again')
    : p5 M( g4 m. Y5 }$ ]9 j
48. end. s8 @5 L, `: q# W  s9 b
    
49. for i=1:length(x0)% r, ~( K% ?& u! T& a# e
    
50.     Hatx00(i)=Hatx0(i);8 J% O3 a( K; v* A7 A* @
    
51. end
    . }$ C: j- R& B( K& z! V* f1 }
52. z=1:length(x0);
    + }) [2 u; @2 B9 G- Z' g9 v7 H
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察5 [# c& H- O! g# u  ?+ _% m+ J, y
    
54. text(2,x0(2),'History data: real line')
    + K- x7 m# Z- w! J( Z9 u/ y
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')
    6 c3 A* i4 _. D5 H
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2]2 N9 Y: I! f: h+ M2 L
   
2. Warning: Input arguments must be scalar.0 t5 M( J8 |' @6 ]& k& C- R
   
3. > In fungry1 at 4
   # N1 i$ f& X\" K9 t9 q
4. Warning: Input arguments must be scalar.
   0 `\" D  q  n% l2 [- W& d
5. > In fungry1 at 5
   9 i/ T) j, V$ {7 P0 o' g9 w
6. Warning: Matrix is singular to working precision.
   8 j1 c/ D9 Q, o) W! K  u) B
7. > In fungry1 at 171 ^* k. C2 J. H  O
   
8. ) x0 t- b9 I' \0 J) J
   
9. HatA =
   ) l4 v2 X) [3 \9 e& h2 N
10. 5 \( j) L) W1 r\" O. t
    
11.      0) S6 m, P0 z5 ^  R) Z
    
12.      0
    & i' R  U1 \- P/ ?, A7 |8 \
13. 
    ! w\" b8 Z; s( [- B/ x1 c3 j( c4 Q5 `* Y
14. 
    # p  u$ Q, i! ]! v\" ~: H3 V
15. p =! {- Z8 |2 [( y$ N* ^! ?, ?& V
    
16. 
    \" B! _; R\" t3 C- F8 T
17.      0: T3 V* r  c7 g\" R2 P# @+ z
    
18. 
    2 q) L4 P- Y/ Q( M' H) @
19. 
    8 B8 z7 ?! z  E, @
20. ans =
    ; i2 B4 b! ~; |4 G0 I( {
21. 9 ^9 b. r+ J\" B( \5 _; t
    
22.      0
    9 h/ G9 X0 s& m4 m
23. 
    , R, }; h' k& c; `- _8 O
24. The model is bad and try again$ S; g/ u% ~  t# V. y
    
25. ??? Index exceeds matrix dimensions.
    5 h  v; S5 h. c5 V
26. 8 {# E$ g( Z0 g\" T% w% G( P
    
27. Error in ==> fungry1 at 548 k/ m# J\" U& _/ Y! C5 x
    
28. text(2,x0(2),'History data: real line')& e) ^3 D, I7 ?9 f4 F% S/ C2 w
    
29. 
    * a9 F4 U' B$ D% N
30. >>

复制代码

是怎么回事。

Mlearner

求高手解答。

zqyzixin

牛啊,想不到的强帖