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

Mlearner

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

luoshichao123

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

Mlearner

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

网上下了个，总出错。

1. function GM1=fungry1(x0) %输入原始数据x0
   3 O7 C' U) H. j: O$ Z1 X
2. T=input('T=');%从键盘输入从最后一个历史数据算起的第T时点  R0 p; }% k. u2 X
   
3. x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);
   * r* X; u- Q/ z4 U5 d2 @! E8 _. K0 r' B
4. yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);
   * [$ d  {9 W& X) G' e\" D
5. Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);, L/ Y& W\" s% h9 ~( j
   
6. epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);
   0 f% |( o/ e: [
7. for i=1:length(x0)
   & u  c/ g) b, Q7 B* j0 w
8.     for j=1:i
   % @+ d; m$ r2 B, j
9.         x1(i)=x1(i)+x0(j);$ m4 O1 M/ o' z4 r3 H6 M\" m6 r
   
10.     end5 A! c9 F/ Q5 J' Z3 u; }
    
11. end8 `! k& Q1 B2 K) |, ?+ o. r0 P; `
    
12. for i=1:length(x0)-1: c% ~. A% S# E. s
    
13.     B(i,1)=(-1/2)*(x1(i)+x1(i+1));
    - s2 y0 F2 N2 \8 k
14.     B(i,2)=1;! x. h7 v9 I6 n
    
15.     yn(i)=x0(i+1);2 H3 Y+ W6 u9 z, ]$ ~
    
16. end
    8 y- h% V- N8 h, L% e
17. HatA=(inv(B'*B))*B'*yn % GM(1,1)模型参数估计
    \" V. Y3 `0 V# {+ T0 {
18. for k=1:length(x0)+T
      d# l; `, k# X4 ~& i. h
19.     Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);6 A5 \! U( `\" [: A
    
20. end1 X6 _- T9 E! L( l1 _! e  a6 K/ l
    
21. Hatx0(1)=Hatx1(1);
    , a% n2 {- w* [
22. for k=2:length(x0)+T
    ) A( t& {! |' A: m/ c
23.     Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值9 K0 n9 H/ e( K( F# v3 N
    
24. end3 N1 ^6 w, B; i( q! [: f! M0 d& S
    
25. for i=1:length(x0) %开始模型检验
    & F, J; P5 I$ e+ U6 I7 a/ l# C) V
26.     epsilon(i)=x0(i)-Hatx0(i);
    # n8 @; E& F/ ?
27.     omega(i)=(epsilon(i)/x0(i))*100;
    - X  o\" ^+ K* o8 ]\" g
28. end
    \" W: s! [\" o& c) x7 K8 o
29. % x0;Hatx0;epsilon;omega;  %必要时去掉％得到各种数据
    ) P7 l3 w; Y3 ^: _; V
30. c=std(epsilon)/std(x0);p=0;7 c& V) ]( Q* {' e
    
31. for i=1:length(x0)
    4 V0 z: x\" U! |$ d: }
32.     if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)
    $ x; b; a6 G. J& d! H
33.         p=p+1;! d. b- ^3 o0 }4 \* _& k2 H
    
34.     end
    & \5 c+ H- x2 m6 m+ J; ~6 N
35. end1 Q# D9 }8 d$ Q$ a' q) V
    
36. p=p/length(x0)
    2 G0 N4 U# v* i: t( S, F  U
37. if p>0.95 & c<0.35
    8 b\" _+ z8 B% D7 N, t\" x* W  L# w! M
38.     disp('The model is good,and the forecast is:'),
    1 q\" J, P% L5 ]& a! X5 z\" @: s: A
39.     disp(Hatx0(length(x0)+T))8 Y+ n% F2 n; D* R  d/ j8 v! Y
    
40. elseif p>0.85 & c<0.5
    ( H( ]\" O# }7 c8 U/ Y& I2 x5 z
41.     disp('The model is eligibility,and the forecast is:'),2 f5 |; \; A5 f+ P1 k' Q\" _  ^3 R
    
42.     disp(Hatx0(length(x0)+T))
    7 A5 o\" y# c! m\" }& r7 c- w
43. elseif p>0.7 & c>0.65
      ^4 r' `$ B; [. F$ m
44.     disp('The model is not good,and the forecast is:'),
    2 R% g. R: e4 z, ^( W9 Y4 F
45.     disp(Hatx0(length(x0)+T))
    . v) D/ T; T! E( ^
46. else p<=0.7 & c>0.65
    , S2 M8 i# n) t  u9 M
47.     disp('The model is bad and try again')2 C( y7 Z3 h# ?3 i( Q\" R! `
    
48. end
    & g. ~$ |, X\" O) `$ k# p
49. for i=1:length(x0)0 g) H) q) Z8 [9 V0 r
    
50.     Hatx00(i)=Hatx0(i);: ?; W0 a5 h! |: _# D
    
51. end6 d; `* H+ M4 ?( G. u* t
    
52. z=1:length(x0);8 v5 O: }+ E1 T' N# V( J
    
53. plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察
    7 P* K4 K9 q; P# G, X
54. text(2,x0(2),'History data: real line')8 ^. D  T% Y6 a4 J
    
55. text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')9 S' h3 S3 I1 X( Q5 v+ a
    
56. end

复制代码

试着输入fungry1(6)出现

1. T=[2 3 4 5 3 2], x' }5 X( `& q! d, ]0 x+ K6 i. p6 r
   
2. Warning: Input arguments must be scalar.
   ; {& x* }: p, X0 j. N1 s7 D; ~
3. > In fungry1 at 4
   $ u; i$ d4 }) b1 Q
4. Warning: Input arguments must be scalar.8 N3 @, v3 d6 c\" I! L
   
5. > In fungry1 at 5# ~* ^& R4 g$ ]
   
6. Warning: Matrix is singular to working precision.; y& N. x# y: d- b7 v! i5 v! v
   
7. > In fungry1 at 176 |+ m# @  j, `$ J$ R$ a1 W. p5 S9 {
   
8. 
   ' \# t4 O, I! E) O& C
9. HatA =- ]$ W; ~/ l3 \6 N! u
   
10. 
    * _$ l# z7 Z7 }  _! D! U* x
11.      0# p% ~8 A4 Q4 g& N2 P
    
12.      0: E3 [3 O) G0 s/ C) P4 X/ T4 V: t+ P
    
13. 4 p1 `  a\" x5 V6 V# H/ H# K' X
    
14. 
    ' E! c- r! ?! V: N
15. p =
    & C7 l6 |9 {1 E
16. 0 l$ @+ c, c' U' y5 s
    
17.      08 |0 {6 D/ P+ K$ k# [
    
18. 
    8 ~; F5 S2 e0 [, |/ [1 a
19. 
    ; [5 n/ X1 ]4 u9 c1 B
20. ans =3 c, n. I) v$ I  @' m! `5 K
    
21. ( n7 ^% t' V: M\" r
    
22.      0$ c# W& }& l  d# c( l
    
23. 
    * p9 |( ^  l4 w5 o  D
24. The model is bad and try again
    4 n9 i, \& U6 C1 c8 l2 t1 }
25. ??? Index exceeds matrix dimensions.  Z$ r; \$ Q3 q1 N! j
    
26. 
    3 q' s0 R& K! m; y* @, v. _) x9 X9 Y
27. Error in ==> fungry1 at 54& N1 v& H5 l3 V0 X8 f/ t) G2 ?6 ]
    
28. text(2,x0(2),'History data: real line'), W9 s' u! S7 s+ S& u. Y
    
29. $ X( z' q  ?% L# q3 @( e
    
30. >>

复制代码

是怎么回事。

Mlearner

求高手解答。

zqyzixin

牛啊,想不到的强帖