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

Mlearner

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

zqyzixin

牛啊,想不到的强帖

Mlearner

求高手解答。

Mlearner

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

网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x0% l- m4 K, |) K. w
   
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点
   & |! Y0 I4 q; q
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);
   5 q7 V- G1 ?' f; y8 Z3 U7 E+ i: N
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);5 d. B3 E* o0 Q9 W) b+ S% v/ x) F+ l
   
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);
   1 X. K6 P2 a9 e( T, U2 a' l
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);6 }; ^, ~0 i6 W- g* F  v
   
7. for i=1:length(x0)
   - O- _( H9 {1 G  Z/ o
8.     for j=1:i
   2 R6 d, G. v! d
9.         x1(i)=x1(i)+x0(j);
   7 B& N! I( i& o; v
10.     end
    \" J4 h- k\" i! ~6 d
11. end: `/ c' L3 {0 ^& }4 |8 H
    
12. for i=1:length(x0)-1
    / e1 E) I! L. o* t6 a' g; m
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));
    2 i+ a; B, b7 s+ H; \
14.     B(i,2)=1;2 n& Z  ]) j$ p* \2 g1 l, C) E
    
15.     yn(i)=x0(i+1);
    9 N8 o$ Z2 }3 c: Q- j
16. end2 f( p\" Y1 X6 N0 V+ @& u
    
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计
    ( D( z, t7 ~* I) ^
18. for k=1:length(x0)+T  ]\" |2 e% ]\" B% w! ?/ e: {  u
    
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);6 v- X8 X# |; O) q9 ~# |  Q: c. T
    
20. end
    9 D& m! F\" I% P2 V. \# I* P
21. Hatx0(1)=Hatx1(1);
    % r7 g( `  A  I5 c  _0 s* s
22. for k=2:length(x0)+T0 }9 T, k) ^; f
    
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值
    % j. d4 V6 h! d& |! x! C* `- ?# E3 `
24. end
    ; g$ L8 U8 Q) E% Y1 M
25. for i=1:length(x0) %开始模型检验7 P. F5 V& S( M0 A  Q9 H5 b
    
26.     epsilon(i)=x0(i)-Hatx0(i);. n8 \5 {+ i) Y9 @$ o+ r* @# n! r
    
27.     omega(i)=(epsilon(i)/x0(i))*100;\" x2 @( I2 s; @3 x. ~. h
    
28. end
    \" s+ X! U/ V' ~8 U
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据
    - w7 j7 N: D5 O& R; C  f- n
30. c=std(epsilon)/std(x0);p=0;9 M\" f3 u8 ^5 y, D' E4 L* X
    
31. for i=1:length(x0)
    \" R8 ^7 Z. u$ ?: n; A
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)
    4 K  N- d1 h2 e
33.         p=p+1;9 N5 L3 r# @& ~$ Z  u4 F. P
    
34.     end/ F% I0 \0 ~( o! k
    
35. end
    # F/ k1 G/ T* p0 A, v
36. p=p/length(x0)6 a' }\" W2 N) [, L) k
    
37. if p>0.95 & c<0.35' R+ F6 N5 V( `. S( }
    
38.     disp('The model is good,and the forecast is:'),
    & H# d\" ?/ m* B6 P: k3 y
39.     disp(Hatx0(length(x0)+T))8 p* U/ r. Z3 |/ T. I
    
40. elseif p>0.85 & c<0.52 E1 E7 H  J- ]+ `8 F4 K  M; |
    
41.     disp('The model is eligibility,and the forecast is:'),4 P* `& g2 b1 T
    
42.     disp(Hatx0(length(x0)+T))- E3 N( v% v+ t9 S8 g9 y( \  P
    
43. elseif p>0.7 & c>0.65% @- }' b. b# s; j, p$ x: M
    
44.     disp('The model is not good,and the forecast is:'),0 |  o- C' Z+ O
    
45.     disp(Hatx0(length(x0)+T))2 x# ?% E' d+ y  C
    
46. else p<=0.7 & c>0.65: p4 U: M) C9 B! e6 v3 h2 M( t8 A
    
47.     disp('The model is bad and try again')
    , {4 V4 O- W7 V9 g) j1 p& P  F
48. end
      A1 @\" e+ d: g/ ~\" O, b% g
49. for i=1:length(x0)
    ( J: e' `9 n2 }
50.     Hatx00(i)=Hatx0(i);! r\" f7 S4 C: T% [+ j
    
51. end
    6 G9 B& h$ c6 h/ H; B9 V' R
52. z=1:length(x0);
    7 }' ]5 [( F+ m+ ?/ Y
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察
    , i8 R- x' X. e: u* R8 y: x
54. text(2,x0(2),'History data: real line')* \* S% z! A$ j
    
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')
    8 h\" S% l. C: q) L
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2]; R, D& _9 ]9 w7 {. C( ~$ {: F4 g
   
2. Warning: Input arguments must be scalar.
   $ i6 u  b/ E9 f* n
3. > In fungry1 at 4
   5 }$ E4 e' N9 T  X8 p: h4 S
4. Warning: Input arguments must be scalar.7 s! @$ `! }3 \0 j
   
5. > In fungry1 at 5( K. }  x( p3 j  \7 o/ o7 c
   
6. Warning: Matrix is singular to working precision.
   \" s# i% ^& I5 Q$ p& Q
7. > In fungry1 at 176 @; p& Q% h, U* B: h. F% W  A
   
8. * m' X7 \' S  G8 u- ^2 p' g2 {
   
9. HatA =: ^0 d0 ]6 V/ w2 O3 z! B
   
10. + D! L, e/ m% i
    
11.      0
    $ m# o/ f\" ]; C, `, p7 @$ ?7 |
12.      0) p7 o5 E: ~2 m: L$ ^
    
13. : u0 u% T\" {' z( g, l- h' L4 E6 @
    
14. 
    - T2 o! B8 p* o
15. p =) P6 E* J3 z0 ]8 N4 p8 s\" a4 z) V
    
16.   \5 O5 l  x5 f# y3 N: u& W1 h
    
17.      0
      a) _* g/ T, K. L9 x# K8 Y  e8 \) O
18. 
    . e\" n: i4 }1 F. u
19. 
    ) e& n2 |; |. c\" H7 ]: M7 [
20. ans =
    8 a8 ~6 O2 X\" ?5 y7 A
21. 
    1 W; N; F6 g\" f\" v. ~8 F
22.      00 P0 w1 Y6 A# O$ d
    
23. 
    % u2 t7 @& k0 x+ C
24. The model is bad and try again
    * l3 S% C6 p3 Z9 L+ T8 y8 [* F
25. ??? Index exceeds matrix dimensions.2 E2 L( o0 P; x2 d0 @7 g
    
26. 
    ' `6 s/ H' R) i6 v) [& v
27. Error in ==> fungry1 at 54
    . X, y+ w6 K  J8 v0 v  Y2 m
28. text(2,x0(2),'History data: real line')9 t; {6 f, ]+ O% d+ t. b
    
29. 5 Q  Z( w/ Q( d+ F( \4 u6 F
    
30. >>

复制代码

是怎么回事。

luoshichao123

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