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

Mlearner

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

luoshichao123

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

Mlearner

luoshichao123 发表于 2012-5-16 07:32
! l2 A: m' P: G3 D% Q  p* T4 x' K7 l这个程序自己编啊，原理很简单的

  i5 Y. \5 |3 @# s3 I* O0 V网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x0) B+ O! [( F\" C  _0 H( e
   
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点9 @. L5 W\" ~7 J5 B
   
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);: b$ U. a2 l, u; T% |7 B
   
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);0 I% t* ~1 G8 n: U4 h8 {
   
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);% w) ~3 i# ^6 M2 G# P) Q
   
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);) J9 M6 t9 i5 ]  @, k# `
   
7. for i=1:length(x0)
   7 {5 u6 S, Q6 n
8.     for j=1:i
   5 t( C% U6 F- a  }0 ?) F$ g$ z9 a
9.         x1(i)=x1(i)+x0(j);\" J/ N! z: P, D6 s5 x# _% T
   
10.     end
    4 E8 ~7 g. A9 w& e
11. end
    ' Z4 q0 K\" E1 y/ C  C8 K2 O
12. for i=1:length(x0)-1
    ' B/ b: p: b$ n6 j3 r$ j
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));/ ~1 K  ~  Z2 v# x
    
14.     B(i,2)=1;
    ( x1 p6 d; K: N+ z) N, ~8 }  C
15.     yn(i)=x0(i+1);
    5 t5 ]\" j: c( X4 [7 Y. e, |
16. end
    ' f7 {0 a2 c\" ]: g% A6 j
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计
    \" g  O% K- f& A4 v1 N. I$ z
18. for k=1:length(x0)+T5 _. N+ s5 M5 T3 [) V
    
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);
    - F2 Z\" L' g% S1 Z0 G! T
20. end: \/ L$ {1 M, U$ X$ K! d3 g$ A
    
21. Hatx0(1)=Hatx1(1);
    * F% C, E! E: z, p3 d) }
22. for k=2:length(x0)+T
    # N# W' z# \) |6 b  X6 w  \
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值
    / c9 F: n; a% i, s( a% F$ h
24. end
    ! h2 [\" i  K+ _7 U: J1 i' e& G+ K0 f7 O
25. for i=1:length(x0) %开始模型检验. _! t; h5 \) b0 }% T  o
    
26.     epsilon(i)=x0(i)-Hatx0(i);8 V1 ^: f7 y$ ~& g2 w- c
    
27.     omega(i)=(epsilon(i)/x0(i))*100;
    6 e$ c6 p7 X# z; A
28. end
    . {\" Z  ~  d% _6 V
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据+ P) T+ J0 N( o8 N  H  J
    
30. c=std(epsilon)/std(x0);p=0;
    9 t( {; h7 d; E* F# j; X
31. for i=1:length(x0), }) j& H8 J7 H; f* |( n1 l
    
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)* i( ^6 j# S: e
    
33.         p=p+1;
    1 g0 B% X8 v\" T\" S& s\" R; A) ^8 S. C4 z
34.     end
    : [# Y! c- w6 t
35. end5 j* ^: V; g2 X9 u$ n\" ^/ X! u
    
36. p=p/length(x0). `, W- B( T  K1 B# X\" `3 ^
    
37. if p>0.95 & c<0.35* |, L# i1 L8 s
    
38.     disp('The model is good,and the forecast is:'),, I) D( x# P) E9 w  x9 w
    
39.     disp(Hatx0(length(x0)+T))3 H$ y& S7 n1 {
    
40. elseif p>0.85 & c<0.5- G3 A% E0 I+ L- A. r
    
41.     disp('The model is eligibility,and the forecast is:'),
    ( E, M+ x8 o4 h! Q) A3 p
42.     disp(Hatx0(length(x0)+T))- C: @8 E( u0 \( r; r; W1 S
    
43. elseif p>0.7 & c>0.65
    \" P  r( a6 b! X
44.     disp('The model is not good,and the forecast is:'),1 V  V2 \1 O! V) c
    
45.     disp(Hatx0(length(x0)+T))% [2 J5 D9 m/ D7 d
    
46. else p<=0.7 & c>0.65
    1 s+ n6 C\" J% w0 i( _
47.     disp('The model is bad and try again')% p% D! q; G: X9 N
    
48. end
    8 }1 G( X7 y2 t! |
49. for i=1:length(x0)- j; P$ z3 a\" T  q+ ~
    
50.     Hatx00(i)=Hatx0(i);
    $ j+ ?, M' v2 x8 \! h
51. end* L* R5 v) l2 R& T: k# |\" N4 A
    
52. z=1:length(x0);
    ' s8 r2 V* m6 @/ g
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察
    ) y% n, Y, w: e\" q* j) U
54. text(2,x0(2),'History data: real line')4 L& |# I2 [! B% s' q) K
    
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')! r8 D! n* j5 E# \1 Q  H! n
    
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2]
   7 T* q: |% f- ]6 B# ~5 N1 r
2. Warning: Input arguments must be scalar.
   0 U$ ~' ?& i  }/ V+ G& X
3. > In fungry1 at 4
   8 v6 F( {5 W6 y5 @6 h5 {
4. Warning: Input arguments must be scalar.
   * H5 V  W2 T+ \, d6 `* \0 \
5. > In fungry1 at 51 ^/ B7 n: X2 o% f0 B# L: s+ E
   
6. Warning: Matrix is singular to working precision.) Z2 v. M# l8 ~. j% `9 b
   
7. > In fungry1 at 17
   % w6 \1 ?9 [) H, L+ h
8. 
   # ]$ o/ D( \* q3 H0 I. }
9. HatA =; `; p+ a# N! _7 w( P9 e  k
   
10. & y# y+ J& ~5 g+ S
    
11.      0+ v4 w3 j+ n7 s; A7 z
    
12.      03 \# n' q! b& j$ D9 P
    
13. ; u- \4 P) r- Y) E2 l
    
14. ( G; ~, _6 ~4 L\" o! Q0 e- U
    
15. p =
    * u\" r/ C2 a9 a) V% X/ Y! j( a' f
16. 
    + d1 a: p: x: U) _' n
17.      0
    3 t& ^8 R! j% X3 z! Z
18. 
    - ]3 s  Z4 O& S: K. N  o\" E: i
19. 
    ! G& V, q  r0 h. L4 o3 y
20. ans =
    6 ^/ N. g  _! \/ r8 g# f5 i2 I' P+ U
21. 
      F3 B$ I' r& Z9 c- K8 T
22.      0  a* V2 V' j6 B% m$ G4 j8 J) y$ Z  G! @
    
23. * [; E0 |6 u7 w0 q6 P
    
24. The model is bad and try again
    ' x; K' x/ e7 I' R
25. ??? Index exceeds matrix dimensions.2 y2 C7 I' Y0 t& l' U& ~' D  }2 F
    
26. 
    ; D4 f* G  R! d7 }7 q
27. Error in ==> fungry1 at 54( m% Z9 n) a9 `4 h- h1 b9 L
    
28. text(2,x0(2),'History data: real line')
    5 K5 g6 c! ~6 v3 @
29. 
    0 N\" B2 |+ Z  s* y\" H# `1 e
30. >>

复制代码

是怎么回事。

Mlearner

求高手解答。

zqyzixin

牛啊,想不到的强帖