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

Mlearner

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

zqyzixin

牛啊,想不到的强帖

Mlearner

求高手解答。

Mlearner

luoshichao123 发表于 2012-5-16 07:32
: Y) ]4 z* m- e3 p5 l) v9 Q9 ^这个程序自己编啊，原理很简单的

网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x09 N; o1 z3 Y# x* o) U
   
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点
   8 g# |4 P. @5 i\" z! f
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);
   4 N$ P2 B3 I2 H, J& G+ O
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);
   2 R* D1 H. P1 D. k+ b5 m
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);
   5 I* s$ {9 B& a6 ~: t5 S
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);
   / `+ r. N+ B1 u, P$ K5 r4 y
7. for i=1:length(x0)0 J# D% d8 E- q1 f3 x0 o
   
8.     for j=1:i
   + F3 D6 O. r/ E& A* e
9.         x1(i)=x1(i)+x0(j);5 H% u& R: _7 B: d, T. q# }
   
10.     end
    ; Z7 H4 ]6 K# ]8 I# C, S\" S
11. end4 D  {4 o& x: f% G! k
    
12. for i=1:length(x0)-1
    6 x- x) u! S9 H9 ~\" Z
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));\" E0 a9 @% g# b4 Z2 D& o
    
14.     B(i,2)=1;
    + {# x$ J4 o\" _- B+ r
15.     yn(i)=x0(i+1);0 U- z; u2 c/ _; Q; w$ N$ e
    
16. end
    5 T2 V3 x4 S8 e7 Q& w4 \- ]
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计4 P$ i& P9 l+ I; N# j
    
18. for k=1:length(x0)+T
    ; }* d5 F\" A\" K- J) c
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);1 T6 e# G, Y) I1 L
    
20. end
    & i# D& \. u; ?. U
21. Hatx0(1)=Hatx1(1);) }9 i1 A- [2 h( N; I* ~9 G& i
    
22. for k=2:length(x0)+T
    5 J/ {: L7 u) L1 V$ V$ Z& @1 ?
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值9 {* I2 p2 ]\" r* j. T+ A
    
24. end
    & R- ^% `; D* }1 }6 S! o# \: U
25. for i=1:length(x0) %开始模型检验+ J% k, C1 u3 |) G2 e
    
26.     epsilon(i)=x0(i)-Hatx0(i);( M0 }0 O: {$ y\" |
    
27.     omega(i)=(epsilon(i)/x0(i))*100;
    * f: Z# _' U9 g8 O3 u+ f0 y
28. end
    6 t# ]1 T0 n  h+ x+ x5 B! j  L
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据
    \" w. [- M7 v$ @8 ^\" ^4 d
30. c=std(epsilon)/std(x0);p=0;\" ^: t, q1 m) {4 [5 C- o
    
31. for i=1:length(x0)
    0 u/ u3 ^( y, [7 v+ L5 s. p
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)
    6 g* Z& w1 R2 z4 z\" Z& \
33.         p=p+1;
    ( p# \9 f1 G  H3 h
34.     end: H) F\" T, d6 F6 ]4 b9 u, D0 T
    
35. end! v2 q5 `; _' f3 D
    
36. p=p/length(x0)
      U! ^0 e7 Y8 S3 X5 G  \2 n
37. if p>0.95 & c<0.35\" ]7 ], T. }# e
    
38.     disp('The model is good,and the forecast is:'),; `2 H; T6 [, U; Q/ N8 C, ]4 M
    
39.     disp(Hatx0(length(x0)+T)): G0 `! y& H' \! c& _4 q/ V
    
40. elseif p>0.85 & c<0.5
    ; E$ d& V  U3 _0 n! ]! r/ z. D# J
41.     disp('The model is eligibility,and the forecast is:'),8 W4 Y$ }\" b: D. a3 x! b' z0 w\" L  }
    
42.     disp(Hatx0(length(x0)+T))0 d1 g$ Y3 u7 t* A+ v
    
43. elseif p>0.7 & c>0.65
    8 q5 g5 a/ L6 B& O4 g
44.     disp('The model is not good,and the forecast is:'),
    # J0 j4 M) P8 t% D7 ?
45.     disp(Hatx0(length(x0)+T))
    . r& e& n! Z$ f, ]
46. else p<=0.7 & c>0.655 M/ j' L( ~\" B\" t' {4 E
    
47.     disp('The model is bad and try again')0 _\" D. ^8 j7 ~7 Q, V/ E
    
48. end( v: c6 f( E0 U0 y
    
49. for i=1:length(x0)
    - M6 I. a0 @$ h  M+ r' G
50.     Hatx00(i)=Hatx0(i);* [$ e% q; f/ ^7 |$ T/ r+ }, Y( q8 F
    
51. end
      l& X7 I) S' y$ Y1 t  O# b
52. z=1:length(x0);
    0 [1 Y0 u' n3 e; G( L- D
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察
    4 {1 P% w' r: H
54. text(2,x0(2),'History data: real line')( R, D% |  l0 w
    
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')& G. s0 n\" c. X
    
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2]
   9 @( r4 q1 m  X\" J8 Y* ?
2. Warning: Input arguments must be scalar.5 u3 Y& N8 `8 G: N4 V# R8 a4 _
   
3. > In fungry1 at 4  u7 K: s4 U' r% e2 f5 x
   
4. Warning: Input arguments must be scalar.7 |  p- p; F, ^# U/ p& l6 A
   
5. > In fungry1 at 5; ~& u: _, ]2 V- b) Q
   
6. Warning: Matrix is singular to working precision.3 w& p3 `; L+ r  D- C/ p+ A
   
7. > In fungry1 at 17
   3 _5 G* e- N; w. E: {
8. + S) X) J3 W* i' N5 F: L) P
   
9. HatA =
   1 j/ B4 Y/ d2 I2 O7 e, a1 h; P
10. 9 F$ N0 H# Q\" H  Z  _# g* a& |
    
11.      0  e: v+ b7 o1 c6 W# T% G* S7 M$ k
    
12.      0  |+ r2 k+ \$ u2 t7 V$ L1 h
    
13. 3 ^5 q  r: _/ O& Y  J6 D
    
14. 
    $ K. X0 q( @. U
15. p =. `9 F# _\" ]( Q& n
    
16. 8 d' i7 Y; H% q, E% v\" ^, T. k
    
17.      0! k2 Q$ l8 C5 h5 c2 t+ A9 A/ x
    
18. ' ~& s, X1 b3 k; J1 O- h' s9 D
    
19. 
    ! I3 J. ]\" @9 @( g6 D
20. ans =
    ' g7 S# Q  q. ~4 w3 q
21. 
    8 F8 S( P1 Q; A* y
22.      0
    8 `\" d0 r# m* W2 Y6 L$ |3 {
23. 9 [, J6 z! O. S$ @4 Q
    
24. The model is bad and try again% v% _* c) q9 Z% Z
    
25. ??? Index exceeds matrix dimensions.8 l\" ~\" b0 a7 l0 w$ d
    
26. 
    8 d6 T( e9 v- p' N& G5 {
27. Error in ==> fungry1 at 54/ V* |+ n: P+ r) b( d
    
28. text(2,x0(2),'History data: real line')3 ~5 ]5 \/ n% \/ u/ u2 |
    
29. 
    , B, T: U\" m0 J
30. >>

复制代码

是怎么回事。

luoshichao123

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