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

Mlearner

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

luoshichao123

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

Mlearner

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

网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x0
   : [. Q! M3 L6 a! ~2 m
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点, |( k\" {+ M  R5 r4 @
   
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);
   , }4 x0 t% y! p- U$ y4 Z
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);( Y- ]: _, m5 p! t
   
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);8 a6 l5 V  v\" s; e9 d$ g/ F
   
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);
   : u9 @- C8 w7 a
7. for i=1:length(x0)7 N% L' I& |& C2 z1 _$ a
   
8.     for j=1:i
     t$ ~2 U4 [; s  W4 ]$ P
9.         x1(i)=x1(i)+x0(j);7 S$ O  P2 N& v* I! l; Y
   
10.     end
    ; j0 Z$ @\" Z3 T
11. end
    0 R2 X* M( m$ @+ |; z( U! S8 E
12. for i=1:length(x0)-1
    / w8 s. s5 p4 o6 `8 \: l
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));# T5 C$ E* q5 _- ~. D; c
    
14.     B(i,2)=1;( p\" f7 V6 o3 M) `$ H  R7 O
    
15.     yn(i)=x0(i+1);  {& Z3 f. F+ l1 X8 d  i4 m8 F. O
    
16. end
    ! b( s; v& ]) Q5 M; M1 C
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计8 s, u7 O7 H) t( G
    
18. for k=1:length(x0)+T+ H$ i; ^5 Q3 l- u$ C
    
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);' k7 E\" e# r( F9 n
    
20. end3 R6 X1 p* X3 h8 p\" v! q2 J
    
21. Hatx0(1)=Hatx1(1);& I8 P, U7 l$ V' q5 V( o% h; E
    
22. for k=2:length(x0)+T
    & {4 [/ h7 ?( c$ P
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值. T/ t, O' l1 ^) I1 i, w\" @
    
24. end
    9 `  u1 w: |8 X3 |4 \
25. for i=1:length(x0) %开始模型检验
    : x) Z8 K* a5 u+ Z
26.     epsilon(i)=x0(i)-Hatx0(i);5 g: ?* b2 q& V6 B
    
27.     omega(i)=(epsilon(i)/x0(i))*100;3 t% g. j- o. H- n
    
28. end, I8 h\" _' D( e* Z% @
    
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据
    * c9 c7 g\" u! n0 C% N
30. c=std(epsilon)/std(x0);p=0;
    3 N1 a9 Z3 W5 f. l0 d( w
31. for i=1:length(x0)
    8 b- v8 F% v5 T
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)
    , e- S* k: a! V# h1 d
33.         p=p+1;; `7 N0 K; c5 }0 k2 P4 k+ j) ~. i- t. T
    
34.     end
    - t# a, ]- b1 e6 q) f' N
35. end
    & L2 O* c; a( ]
36. p=p/length(x0)
    $ L% u2 q4 g2 ~- y1 \7 j3 j7 O
37. if p>0.95 & c<0.35  Y( q$ B3 d/ U3 M
    
38.     disp('The model is good,and the forecast is:'),
    ' M# w  h  i' j3 c# D
39.     disp(Hatx0(length(x0)+T))
    : j/ @7 U: |4 r4 Y2 f; d! P
40. elseif p>0.85 & c<0.5! A. U9 O3 K* k( k& e
    
41.     disp('The model is eligibility,and the forecast is:'),, @% c) d( L: M, ?  C
    
42.     disp(Hatx0(length(x0)+T))& B6 v' E* t' Q$ w
    
43. elseif p>0.7 & c>0.65
    ( s3 e! _3 I\" S
44.     disp('The model is not good,and the forecast is:'),( i% j2 s5 y5 }
    
45.     disp(Hatx0(length(x0)+T))6 |( n. S) r6 w$ a0 b
    
46. else p<=0.7 & c>0.653 g$ ^+ w  e- w\" K* [
    
47.     disp('The model is bad and try again')
    4 d2 E! i1 |7 u- N# d7 ]
48. end
    ) ~/ q$ ^! t: I; C$ q
49. for i=1:length(x0)
    \" X, W* T1 \3 x3 I$ L/ E\" ?$ w4 T
50.     Hatx00(i)=Hatx0(i);
    , x3 c9 I( O% `8 `1 V
51. end9 `; b- e! C! H5 @+ ?
    
52. z=1:length(x0);7 h4 ?9 t% ?, G) f/ l
    
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察5 y  u$ F  H4 {
    
54. text(2,x0(2),'History data: real line')
    0 \- r1 X* c1 W2 Q4 U
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')( Z! R+ `7 \* R1 j' W, \3 R
    
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2]
   # |) J8 z  a\" Q. l) S$ m' i
2. Warning: Input arguments must be scalar.
   . i4 W, i% ?\" W  y+ [
3. > In fungry1 at 4
   / z/ O) M2 y4 m
4. Warning: Input arguments must be scalar.( J3 }: v, d6 c
   
5. > In fungry1 at 5
   3 n& r# r0 C- e
6. Warning: Matrix is singular to working precision., \2 M( N! Y4 d* W7 @4 @4 v( K2 J
   
7. > In fungry1 at 172 _' R4 `; ~/ P& M4 @
   
8. 
   ( _) u9 b0 X: q! _$ d1 _  c0 G
9. HatA =
   / x7 F$ b; v& `+ C' _, X
10. 6 L0 T# g& v1 b) b6 ^
    
11.      0
    7 R9 x% M0 m8 @2 c8 U1 W, i; d
12.      0
    / v! f7 n( q/ U! ?' |$ N6 c5 f! l6 o; M& V
13. 
    ! I6 V+ ?4 N( [1 @4 u0 L, T
14. ( E' X1 ^$ J  y9 M, d\" i  i
    
15. p =
    4 {; y+ r# s! U, x
16. - ^5 o\" ^( E2 A. ^7 u
    
17.      0
    + |' Y6 P9 _+ I) Q
18. 
    ! G/ k2 y4 Q4 A( T9 k4 I8 t: f8 y
19. 
    1 B! O7 f7 v8 f7 W# k! p! m. s/ b
20. ans =
    ! B6 z; d% [; }: l\" p+ l  v
21. $ n4 k  l- c6 A. m9 W
    
22.      05 f( F, I5 Z7 g6 T8 S$ b1 J8 c# Z
    
23. ( i# p& Y9 H! F: z+ M. [+ ]+ g
    
24. The model is bad and try again
    $ [2 g8 ?  G7 P5 h
25. ??? Index exceeds matrix dimensions.: D' p1 A. @! _. n- M7 u, V% S
    
26. 
    6 I3 k; B! t2 M1 P
27. Error in ==> fungry1 at 54
    # I  F4 r3 Q1 `8 K) o1 F% W) `! g
28. text(2,x0(2),'History data: real line')7 l, v$ F. q4 F7 B
    
29. 
    : v( t' [- h! f, A8 w\" x, F
30. >>

复制代码

是怎么回事。

Mlearner

求高手解答。

zqyzixin

牛啊,想不到的强帖