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

Mlearner

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

luoshichao123

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

Mlearner

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

! e, F! v0 \0 {' W! @网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x0
   0 l& N6 `  T4 o9 |3 Q, _
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点, O. i\" N1 o6 }- e# |% E; e9 z& L
   
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);% S8 z6 U( l8 G8 K( }8 ]0 w7 P
   
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);
   # A1 N2 V; d& v\" B: d3 q
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);
   8 d6 J4 p. f3 X+ ]) f) q( R+ }
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);& j8 T; G5 A& ?, d2 r! Q8 M
   
7. for i=1:length(x0)2 d! b' G- d, b: W$ G6 G0 `* s
   
8.     for j=1:i: P8 \( z% r. c% A5 M
   
9.         x1(i)=x1(i)+x0(j);4 [* n) i% ]% S# H. L3 U
   
10.     end4 y* S7 n1 x& k+ @9 z& P
    
11. end
    / L- e/ x  L# x2 C\" r
12. for i=1:length(x0)-1. `- `  N2 w( l# F1 Q. z
    
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));
    & @: _2 v# x' d( S  G# u! @2 l
14.     B(i,2)=1;
    . n: W) i/ E, w5 [& r
15.     yn(i)=x0(i+1);; p5 N8 u: R9 [
    
16. end
    & U+ ^\" x6 d# z* g$ M. _6 b
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计# W* s4 z* O4 z1 h* W$ M
    
18. for k=1:length(x0)+T
    . _3 J, S4 x: `! S4 [& S
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);
    - z$ K% }! s2 [* H6 n% a
20. end. G9 ?& o: n* a' S* R+ K' @
    
21. Hatx0(1)=Hatx1(1);
    6 J% ^  n1 |' X. j8 y0 z
22. for k=2:length(x0)+T
    9 @- Z' V' h) E! I
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值  g! V% z, }  u
    
24. end
    ! v6 b: I# C, Q- m+ p9 x
25. for i=1:length(x0) %开始模型检验  w$ L/ S0 g& @
    
26.     epsilon(i)=x0(i)-Hatx0(i);. J: ]# `% r: s9 G
    
27.     omega(i)=(epsilon(i)/x0(i))*100;* a! V3 _$ E/ N* L1 S2 i/ s\" x
    
28. end  O6 h/ g: _; l3 D
    
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据
    \" I7 m  g8 Q6 @
30. c=std(epsilon)/std(x0);p=0;2 I* `# {5 x+ @! n
    
31. for i=1:length(x0)+ \! Y* H8 O7 u8 q4 ?# j2 h
    
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)
    ( Y* \- H& u* X! m/ |* c; }* f
33.         p=p+1;5 J3 d( C( z) Y7 k0 v- P
    
34.     end6 k9 _: r! x: _* _
    
35. end
    ; Y! o7 ]8 O+ I. ~/ ?
36. p=p/length(x0)
    2 H' @+ a- S& S+ U
37. if p>0.95 & c<0.35
    . e* a8 S2 d6 A9 W+ |3 ~) n, y
38.     disp('The model is good,and the forecast is:'),$ \3 ]; ]- j$ k2 x  C2 R, e
    
39.     disp(Hatx0(length(x0)+T))
    4 [  T0 f, `* s4 F6 X. s
40. elseif p>0.85 & c<0.5# o/ @: Z+ A* d
    
41.     disp('The model is eligibility,and the forecast is:'),* ^$ D7 x: u; a2 K3 M
    
42.     disp(Hatx0(length(x0)+T))) g- Q\" P: W# F# x- c4 i
    
43. elseif p>0.7 & c>0.65
    8 H. E1 A) N7 m\" G
44.     disp('The model is not good,and the forecast is:'),
    5 X1 \5 j% @+ }! N! h7 x
45.     disp(Hatx0(length(x0)+T))\" k3 w9 m- u- @: S5 G
    
46. else p<=0.7 & c>0.650 ]: L5 p0 c. D) K7 P/ B
    
47.     disp('The model is bad and try again')
    7 g/ G# D7 f8 S! ^8 R, Q- x: d6 z
48. end9 x6 Q% h( `- h8 C$ o
    
49. for i=1:length(x0)
    : `, Y4 z1 p5 ~
50.     Hatx00(i)=Hatx0(i);
    ( F9 A3 C0 a- S- `/ Q. q: s  I
51. end
    : N  V% P7 S: k; i& V/ L
52. z=1:length(x0);& c- l8 l0 t, d: o
    
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察0 O. n1 A+ ^& c. ^6 Q. |4 z- A
    
54. text(2,x0(2),'History data: real line')2 Y! _7 E- X6 ?1 u0 Y& Z1 g$ Y
    
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')( I9 a* M, M- O3 I6 e; {5 ~
    
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2]
   6 j& s3 `/ Q) X. t( o3 }4 w
2. Warning: Input arguments must be scalar.
   - A' a) J8 K% a
3. > In fungry1 at 4' e: e6 @$ o4 N) A/ r' o
   
4. Warning: Input arguments must be scalar.) ^0 O! I. Z3 y. ?' Y% G
   
5. > In fungry1 at 5
   . a9 Y) W6 N1 k& T# R  U. Z
6. Warning: Matrix is singular to working precision.
   8 Y. K, Z! K& T  s# L* [; g  h
7. > In fungry1 at 17
   ' J3 S' t  @% X
8. 
   $ Z+ \( ^! k/ {4 n\" r( Y2 f& @0 I! a
9. HatA =
   $ v, f. ]# j1 l+ _2 Y
10. 4 r9 m/ k$ w3 D
    
11.      0
    ! A- k: v+ V* `% b/ l+ i. `) D
12.      01 r# Z: }; h, r3 V
    
13. 
    , e  }) ?9 ~7 n3 \
14. 2 f- [( K) I6 W0 |  I
    
15. p =5 K( _5 ^& h2 S/ q2 x
    
16. 9 H1 q9 ]\" [$ E) j& p
    
17.      0
    , s9 ^& z* l4 ]! c$ r
18. 
    ! y7 ^' m: a7 n+ d3 r  w\" d
19. 8 D* B6 z8 p' `, F
    
20. ans =
    ! P$ V& l$ P* C\" p9 S
21. 
    % G2 T$ e+ n5 x! a# n
22.      0
    1 I  J( u5 A! Y* z1 ]
23. 
    1 [  ~2 T6 ~# X8 O+ l
24. The model is bad and try again9 h2 x* q( p' u
    
25. ??? Index exceeds matrix dimensions.  ?, ?$ T& ^/ ]2 |+ x: S
    
26. 8 C5 g( z$ S& R, r
    
27. Error in ==> fungry1 at 549 `# Y0 Y; X! C6 B; D
    
28. text(2,x0(2),'History data: real line'): G- M% Z5 N, C( E! t+ R
    
29. 6 I3 ^# M: N8 t
    
30. >>

复制代码

是怎么回事。

Mlearner

求高手解答。

zqyzixin

牛啊,想不到的强帖