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

Mlearner

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

zqyzixin

牛啊,想不到的强帖

Mlearner

求高手解答。

Mlearner

luoshichao123 发表于 2012-5-16 07:32
+ x9 e7 O  Y$ @8 b3 |4 t/ y5 e这个程序自己编啊，原理很简单的

; ~$ \% P8 F. \& w2 }( f网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x0  y5 {$ d. p1 l8 S' W# S% y
   
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点
   7 e% n\" ~- C1 {& ^! b
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);
   5 Q4 x5 t% O: s$ d# z) K6 }
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);! n/ T/ t3 g2 s4 Y* K* [4 B8 w) Z
   
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);
   / s* e7 I5 Q0 z# Q/ v, a
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);/ Z& b$ k% p  b
   
7. for i=1:length(x0)
   # @2 j& r3 o2 h: j# d
8.     for j=1:i6 \: W\" l3 J1 D2 z4 h
   
9.         x1(i)=x1(i)+x0(j);
   2 k' @: D- e+ b8 q
10.     end8 L% e\" W; f* Q) J# K/ U( j& O
    
11. end
    . i, K+ }6 o) u6 Q& n
12. for i=1:length(x0)-17 Z& f  ]6 @3 H  b\" {( l
    
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));9 {' i2 u' @. }! G\" y7 V/ S, W
    
14.     B(i,2)=1;, U0 b1 w\" S6 G' k
    
15.     yn(i)=x0(i+1);# B\" @, ]' q8 G4 F5 j
    
16. end6 H/ I( D  u5 ~' n. u' I0 r5 p
    
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计) [: i\" `\" M, Q' b! R* P6 l
    
18. for k=1:length(x0)+T- Q+ ^5 I. p5 t/ Z% l\" B6 j
    
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);
    + L\" d8 ?. p9 i$ z& p
20. end) F8 R. t\" F  V+ Q) n0 C' Z* b
    
21. Hatx0(1)=Hatx1(1);
    7 L0 n) h0 n  h$ x7 I' `' @
22. for k=2:length(x0)+T2 |5 E. q, l% y- v2 T4 c3 V
    
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值1 ^# m\" r: m3 o3 U* \9 }
    
24. end
    7 U/ f+ N! C2 m% E( w\" ?$ u
25. for i=1:length(x0) %开始模型检验/ ?3 f, o1 @) Z0 G+ o- c2 \
    
26.     epsilon(i)=x0(i)-Hatx0(i);
    ! |, ^7 C+ n) p5 x8 K
27.     omega(i)=(epsilon(i)/x0(i))*100;9 v/ I9 V& Y6 Q3 a$ a, V3 z: r
    
28. end
    / a6 |+ Y9 t# g+ w
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据
    . T; f0 W' A# I) \- _2 u
30. c=std(epsilon)/std(x0);p=0;
    + S, O6 r$ s& t6 H\" Z* j
31. for i=1:length(x0)
    ) V5 d1 ~% }\" s, B) i3 v
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)
    - j' x. u2 S6 v) }  o; f3 z$ F
33.         p=p+1;* D$ N& F' ^; k( L' a4 a! j$ T
    
34.     end  M, A3 s$ S$ P; ~
    
35. end
    \" B1 L7 [$ |1 V- m$ C) a# K
36. p=p/length(x0)
    : R$ }& T( X  e/ V
37. if p>0.95 & c<0.35\" W/ u( N9 y% y1 ]: z/ Z
    
38.     disp('The model is good,and the forecast is:'),
    8 D0 ~, ~4 D( P( C9 z( P
39.     disp(Hatx0(length(x0)+T))
    5 q, s, \- _0 e* y. ]
40. elseif p>0.85 & c<0.5' W4 D# z- _* E5 g
    
41.     disp('The model is eligibility,and the forecast is:'),5 b8 \6 a$ o7 A
    
42.     disp(Hatx0(length(x0)+T))
    & }, E9 |+ n+ |6 M3 r; I* `
43. elseif p>0.7 & c>0.65. f\" t# `/ G9 K# I% p- ^
    
44.     disp('The model is not good,and the forecast is:'),, n' n0 z6 \8 _5 t0 g\" k: e- K
    
45.     disp(Hatx0(length(x0)+T))
    7 N/ `% ?# t% q; N\" D1 X4 Z
46. else p<=0.7 & c>0.65% s# z# u% R3 O
    
47.     disp('The model is bad and try again')
    ) V4 D- c: N8 i8 V5 H/ E
48. end9 b7 e- E$ _. j8 w
    
49. for i=1:length(x0)6 _& V- H# o  M; T. o- M( h
    
50.     Hatx00(i)=Hatx0(i);6 \, W3 u6 L% J3 W9 s
    
51. end- U& C8 B# @0 c1 m  ?% N8 D
    
52. z=1:length(x0);! v, z' l- X8 [6 w0 T  n- L! k
    
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察
    5 Y+ V0 x  f6 C  p# u
54. text(2,x0(2),'History data: real line')\" Y7 v2 K' x. c0 b5 b$ j; ]- L
    
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')9 B# w* T3 h9 s# n2 Z
    
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2]3 ?- m: ^  A\" {( ?# g- N2 ]/ p
   
2. Warning: Input arguments must be scalar./ P' n% O& N$ ^
   
3. > In fungry1 at 4
   1 e2 E! {9 J+ }$ o\" V9 D
4. Warning: Input arguments must be scalar.
   % G7 t& X, |2 w0 Q3 K
5. > In fungry1 at 5
   3 c$ R2 x5 y% ^8 i  A3 }; A
6. Warning: Matrix is singular to working precision.  x/ f* d4 j# D' n+ G7 \
   
7. > In fungry1 at 17
   ( _( W\" G* B/ \: o
8. + P, j; v1 u$ U\" L5 w, c
   
9. HatA =8 c. S8 ~. a+ K6 X
   
10. . D0 o; Q\" P3 q2 z# O$ m( x
    
11.      0
    & t( z0 j. A: R# }
12.      0- o1 s. \( V% x2 |4 I
    
13. . N# n1 {6 `  J
    
14. + g  K' D( y5 e; N! x7 C. ~
    
15. p =
    - f' t* D  a. c- M
16. 
    / [\" C( t, P2 U  h' X
17.      0
    ) Z. D6 ]\" T  F' u5 {) m7 p
18. 
    5 R- u* t7 \) ^# s
19. 
    # G4 u4 _4 `# R. Y9 Q
20. ans =  h2 @3 N) p\" R6 Y5 o4 ?* k. ~; N
    
21. 
    6 C- J# a, {- E7 E
22.      0; U\" C2 \4 m. z% t! b: v9 `! @
    
23. : h3 o& @8 ^/ K$ H* ?% J) P
    
24. The model is bad and try again  {  i% e1 E. U% W\" l; L0 u! n
    
25. ??? Index exceeds matrix dimensions.
    7 c. w) s7 Q\" c* `, i
26. 8 _5 t, l8 @# F1 p
    
27. Error in ==> fungry1 at 54; j* u: u: y6 y4 z$ ?2 H; ~
    
28. text(2,x0(2),'History data: real line')
    & R7 C- t. v8 @6 e
29. - `% x0 c% F1 {5 g3 m
    
30. >>

复制代码

是怎么回事。

luoshichao123

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