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

Mlearner

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

luoshichao123

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

Mlearner

luoshichao123 发表于 2012-5-16 07:32
7 E6 o9 _" d- e( u3 N2 X; g这个程序自己编啊，原理很简单的

网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x0% F/ O, a# ~/ ]$ K2 m- M1 M( P, Q: `
   
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点' M6 T3 {3 t5 ^\" {
   
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);
   ) O8 W5 e( g9 `. u/ R# e
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);
   , u% I, k7 I; v/ j
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);
   - Y' [\" l  J; \& K; k
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);
   2 s/ ~) |- V% |\" E0 D$ o$ q\" x- w/ b
7. for i=1:length(x0)
   ' k, b7 Y- k. r\" n\" J: v/ B* B
8.     for j=1:i
   . d8 y1 r( q3 a8 o1 @. `
9.         x1(i)=x1(i)+x0(j);
   5 }( Q1 X. T9 v6 j2 I+ U( ]. r' z
10.     end
    ! u2 @. ^! u% U7 I8 f: a! `' \5 p' Z6 i7 P
11. end
    0 @* U  B' P( T0 }' T( F
12. for i=1:length(x0)-19 d! R/ M( [. D4 e9 r
    
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));# d( O3 W4 Q9 R4 o: Y
    
14.     B(i,2)=1;
    , q& A# \1 j% B: \  _6 V. R2 |
15.     yn(i)=x0(i+1);
    0 J. D- w, Y: m4 _& X* O; u' |
16. end5 ~7 Y+ V' E( ?5 g- h8 r! x
    
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计
    ) s% y9 O\" X) k/ T# V# t' A7 j
18. for k=1:length(x0)+T
    4 i. P2 H' w) y! ~
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);$ ?# G9 K& @  l
    
20. end1 ^& P4 Y! n+ f: Y. c
    
21. Hatx0(1)=Hatx1(1);( @9 ^2 ]: f2 v2 w0 e/ z
    
22. for k=2:length(x0)+T6 d/ c: n7 P* u; U; ]
    
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值- ?; I& K) z& r/ R- T* t7 V
    
24. end
    * M' N, \, _\" c9 R' O& u- N. _
25. for i=1:length(x0) %开始模型检验
    % a! V( K$ e$ S# N) A2 I1 _
26.     epsilon(i)=x0(i)-Hatx0(i);
    ; n& N) t8 z( h0 n5 R4 x, T/ j
27.     omega(i)=(epsilon(i)/x0(i))*100;3 C6 v( D1 ^, Y% k1 x2 Z6 S; X; K
    
28. end
    , E) ~  P$ m5 U0 i
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据
    ) S8 @9 n: C/ j( I* z
30. c=std(epsilon)/std(x0);p=0;* y# t1 F3 [2 f1 R, s( t( `\" K
    
31. for i=1:length(x0)) ?2 ^0 e9 S& L, \6 U4 \$ x
    
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)4 S) @3 a' F0 O% n
    
33.         p=p+1;
    \" \3 F; G- K$ u* a; _
34.     end
    ; F, V! l( t7 g
35. end
    : ~3 X6 e5 M; h% K' W, H
36. p=p/length(x0), e! @/ p/ Y3 g
    
37. if p>0.95 & c<0.35
    6 F- z\" E- h  I
38.     disp('The model is good,and the forecast is:'),& {$ i! H' F( i4 \6 r; N+ D! M: l
    
39.     disp(Hatx0(length(x0)+T))1 C% n! j$ X3 @) J5 I
    
40. elseif p>0.85 & c<0.5
    5 F- S) Q/ q& ^, H* u( P& s8 s; F
41.     disp('The model is eligibility,and the forecast is:'),
      t: B2 o$ R2 }: t3 d
42.     disp(Hatx0(length(x0)+T))) e% M0 e$ _3 y! L; W
    
43. elseif p>0.7 & c>0.65
    + R7 F6 C7 W# n
44.     disp('The model is not good,and the forecast is:'),
    / e+ X: O\" \+ E
45.     disp(Hatx0(length(x0)+T))& d% j, G% o) o
    
46. else p<=0.7 & c>0.65
    4 R! Y4 M' l8 v
47.     disp('The model is bad and try again')
    . a# @9 v' P! A* W( v+ p% q* s( t
48. end
    1 l' v2 j# ^2 O! G
49. for i=1:length(x0)
    ; t% m' @; F  O& u+ X
50.     Hatx00(i)=Hatx0(i);
    \" [3 D/ d1 t; _  H* [6 B! `- C
51. end& u; [) K0 P7 X! ]. k
    
52. z=1:length(x0);
    5 C  r( ^8 D6 _9 h, U% ?6 a
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察
    : r! L- W- y5 d& D+ G
54. text(2,x0(2),'History data: real line')
    / z% i4 S5 n: T9 u
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')( M: T! J2 C0 y% D: \; g0 h
    
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2]
   $ k: ]. f) E4 o4 u2 T
2. Warning: Input arguments must be scalar.
   6 I7 `1 J3 l5 F2 a1 [
3. > In fungry1 at 45 D8 v5 T% O6 B3 R, q* l
   
4. Warning: Input arguments must be scalar.
   $ {2 r1 A4 G3 t, d& r# V
5. > In fungry1 at 5
   - z% T8 R5 a9 o: Q& N
6. Warning: Matrix is singular to working precision.
   0 \8 J) H  Z3 @2 Z& I+ M; l/ @
7. > In fungry1 at 17
   5 [$ W3 Z\" p: G( l5 |+ h; Y
8. 
   ) c: x- K7 ~$ ]# G0 _
9. HatA =
   # R$ \1 N4 W2 u) n: w
10. 
    : G/ V( d4 }\" o1 T* `0 C# w, N# m
11.      0& r) P1 X! S! F8 z
    
12.      0
    0 k$ u6 M$ F8 C2 L- g
13. 6 B; g* t6 j0 q( K# F8 s/ Y
    
14. 5 v# Y, {& w) `/ b2 f
    
15. p =\" s  q! J, p. v\" m4 H
    
16. 
    % v7 V5 J1 v' U7 e9 k3 f' p
17.      0
    ! \% y\" H. ~# f& d0 H: g- n
18. 7 w* c7 ]1 I) @! Y; I; X
    
19. * g/ Z6 a0 y) L) [
    
20. ans =
    * J- M- p( g1 d2 a; R8 u2 z# T2 m
21. 6 g5 n- ?\" h0 {
    
22.      0\" b9 a2 y$ S4 Y1 k9 G) j
    
23. 4 Z$ ^3 C3 Z+ c/ i1 v\" U2 B
    
24. The model is bad and try again
    : k/ ?- B9 l7 v3 r
25. ??? Index exceeds matrix dimensions.! h; l; ~! E& n. ]- W
    
26. 
    9 a. G0 N# V\" Z4 E5 s
27. Error in ==> fungry1 at 54
    ; \7 Q& C5 i  [9 ^! O; E% N
28. text(2,x0(2),'History data: real line')6 H! `$ F2 _1 O- g/ b
    
29. $ C5 i; u# V& {0 A3 y) q5 o( w
    
30. >>

复制代码

是怎么回事。

Mlearner

求高手解答。

zqyzixin

牛啊,想不到的强帖