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

Mlearner

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

luoshichao123

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

Mlearner

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

网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x05 L5 k- }6 W) H+ _; s
   
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点
   % B$ Z, e9 Q5 L: |7 Z
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);: f  i; O3 {3 n' @1 X) Q$ R5 F
   
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);& @* t8 J. O$ z5 t! H* F
   
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);
   4 r. Z9 M4 t3 e
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);5 P8 \; N2 ~% n$ M6 m  A
   
7. for i=1:length(x0)9 n% |# l/ y4 p# L, L
   
8.     for j=1:i
   3 o9 k% F. l* W  _3 w! ~
9.         x1(i)=x1(i)+x0(j);
   8 X0 _8 ~$ U, F: P6 E6 r( E9 @
10.     end
    + Q8 Q1 H1 \7 J$ S2 U3 H. [* a
11. end. I3 f; k% b1 q3 ]; ]
    
12. for i=1:length(x0)-1& P3 x# l\" ?1 h$ E( c5 |, Q+ C
    
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));
    5 L6 y, r4 y1 a% ?5 h  C  x$ P
14.     B(i,2)=1;, ]; s( [' M% Y/ q
    
15.     yn(i)=x0(i+1);
    + R; ~- V: ~0 g8 t0 A
16. end
    \" E, A) P- [( B3 ?
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计
    ( q! c( I; G' K7 L
18. for k=1:length(x0)+T9 [( N* ^# Q# R0 w( `) w* M' x4 i
    
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);
    8 G# [+ ~8 q* N, B
20. end! R5 B+ L% a) S  p5 v# T( F4 |; G. [
    
21. Hatx0(1)=Hatx1(1);
    7 U2 I! L9 M4 p, |! P& P0 {1 m
22. for k=2:length(x0)+T, k\" J0 _( h  U8 `$ r+ j
    
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值% X4 B/ B4 m6 M
    
24. end1 Q& T6 P$ z\" @/ R1 R- s  s
    
25. for i=1:length(x0) %开始模型检验
    7 u\" f7 a% d! {: Y/ o0 Y
26.     epsilon(i)=x0(i)-Hatx0(i);
      e5 S8 U% ?: R( S; ~  V
27.     omega(i)=(epsilon(i)/x0(i))*100;
    0 p# |3 I- V0 u* g% X; b6 }3 W
28. end
    2 b) X+ S- }& w
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据: p2 g. S* `7 E6 C9 [- |7 l
    
30. c=std(epsilon)/std(x0);p=0;
    , ~4 p- n, o- @* h
31. for i=1:length(x0)
    9 t\" h+ U( h: E/ ^\" R$ \, l
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)* ?6 b& a; B# q# h2 X
    
33.         p=p+1;
    ' @5 i- b0 s. J9 _* Y: q6 h
34.     end! d' u; \- `\" M1 |* I) [/ N
    
35. end
    0 H( y6 y\" c' P; R; i+ G
36. p=p/length(x0); ~) n& \6 Q) o+ j' U0 w
    
37. if p>0.95 & c<0.35
    - y: t7 R; i% Z! W
38.     disp('The model is good,and the forecast is:'),
    : \5 X2 @1 y. U; w, a6 W\" {
39.     disp(Hatx0(length(x0)+T))\" @- L& f9 p) L
    
40. elseif p>0.85 & c<0.5  a- A' l6 c9 S6 D9 q) O
    
41.     disp('The model is eligibility,and the forecast is:'),
    : \, @8 ^. F* `  s( ?
42.     disp(Hatx0(length(x0)+T))+ C) A7 a! n9 x' K% U) g3 v  c  C) |
    
43. elseif p>0.7 & c>0.65
    8 L/ G5 Y) m* \\" M
44.     disp('The model is not good,and the forecast is:'),
    7 q% l/ j; Y# k$ b1 l  v
45.     disp(Hatx0(length(x0)+T))
    7 B* O. L+ }2 ^4 ]+ t4 O! ]2 V
46. else p<=0.7 & c>0.65
      R* m* h: N& Y\" ?9 e
47.     disp('The model is bad and try again')! s9 ^4 A7 b9 B
    
48. end
    ( I) c, r9 @\" a/ u: t
49. for i=1:length(x0)
    ' h  B\" C6 u& ]- g1 p- o
50.     Hatx00(i)=Hatx0(i);6 }. k, j- i2 [# Y' x. {- m
    
51. end
    ( }. u  ^' @9 M8 T+ F
52. z=1:length(x0);
    # Y+ i1 |/ C/ n$ S1 }# A
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察
    0 R* l, b6 w# S5 Y* Z5 Q1 r! M) a3 ?! |
54. text(2,x0(2),'History data: real line')' u* q1 Y4 t' v) {8 m9 J2 e8 K
    
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')( I7 f7 r+ Q+ n5 `+ @
    
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2]
   / m\" d3 |' @8 L, @- p: s3 I
2. Warning: Input arguments must be scalar.9 S5 D- I$ e3 a% X
   
3. > In fungry1 at 4
   , B9 r! V. h! j! I) h. U
4. Warning: Input arguments must be scalar.
   * `! F/ R( m/ J; r' H
5. > In fungry1 at 5% b, M3 o! L% D
   
6. Warning: Matrix is singular to working precision.4 V7 E8 j7 g- Z. P0 h
   
7. > In fungry1 at 17
   ) c! p2 m9 x& d
8. 
   , d- k  }  R5 w/ _+ ]\" V5 w
9. HatA =
   : Y* {/ J) X  J
10. 
    8 C' F9 n' s/ V7 I* E
11.      0, ~3 z6 c& v- u! T
    
12.      0
    8 B6 [0 c  ^/ u# t
13. 
    ' i% Q; v& m/ o* ~: x% i
14. % c) Y' Y& J, n% s! K, J  p4 C
    
15. p =
    $ m0 ~: F4 J1 P, c
16. 
    ! c# [7 z& @4 D. S
17.      0
    . S- S2 O+ Z$ q. Q\" ?
18. 
    - Q) F0 S; @4 ]3 R
19. 
    \" s2 q* G. `9 L* B- J: B( I. x+ s
20. ans =\" T8 o* P! [9 _7 ~\" L+ N
    
21. 
    * `. T( S2 u2 \
22.      0- c, R0 O; E/ t7 g* }5 H$ c* G/ i
    
23. 
      W# q; B  ]1 S9 x- m, U. s
24. The model is bad and try again, c9 t9 ^- \& Q. J$ _( B
    
25. ??? Index exceeds matrix dimensions.
    8 P+ N8 |2 r1 ]1 n3 A2 u6 @2 W
26. 7 [* J0 W; Z6 X5 T7 ~& @+ a9 o
    
27. Error in ==> fungry1 at 54
    \" k  `* D+ \, [\" p: s# j: S# Y
28. text(2,x0(2),'History data: real line')
    ' p, O! t' I$ r& X* D) P
29. 
    ! m7 ~! R( a7 O\" B) s- z; W
30. >>

复制代码

是怎么回事。

Mlearner

求高手解答。

zqyzixin

牛啊,想不到的强帖