急求GM(1,N)的matlab程序！！！！

% GM(1,1)模型计胡饥吵算及检验、作图。文肢伍件名fungry1.m

function GM1=fungry1(x0) %输入原始数据x0

T=input('T=')

x1=zeros(1,length(x0))

B=zeros(length(x0)-1,2)

yn=zeros(length(x0)-1,1)

Hatx0=zeros(1,length(x0)-1,2)

Hatx00=zeros(1,length(x0))

Hatx1=zeros(1,length(x0)+T)

epsilon=zeros(length(x0),1)

omega=zeros(length(x0),1)

for i=1:length(x0)

for j=1:i

x1(i)=x1(i)+x0(j)

end

end

for i=1:length(x0)-1

B(i,1)=(-1/2)*(x1(i)+x1(i+1))

B(i,2)=1

yn(i)=x0(i+1)

end

HatA=(inv(B'*B))*B'*yn% GM(1,1)模型参数估计

a=HatA(1)

b=HatA(2)

for k=1:length(x0)+T

Hatx1(k)=(x0(1)-b/a)*exp(-a*(k-1))+b/a

end

Hatx0(1)=Hatx1(1)

for k=2:length(x0)+T

Hatx0(k)=Hatx1(k)-Hatx1(k-1)% 累减还原得到历史数据的模拟值

end

for i=1:length(x0) % 开始模型检验

epsilon(i)=x0(i)-Hatx0(i)

omega(i)=(epsilon(i)/x0(i))*100

end

c=std(epsilon)/std(x0)p=0

for i=1:length(x0)

if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)

p=p+1

end

p=p/length(x0)

if p>0.95 &c<0.35

disp('裤侍The model is good,and the forecast is:')

disp(Hatx0(length(x0)+T))

elseif p>0.85 &c<0.5

disp('The model is eligibility,and the forecast is:')

disp(Hatx0(length(x0)+T))

elseif p>0.70 &c<0.65

disp('The model is not good,and the forecast is:')

disp(Hatx0(length(x0)+T))

else p<=0.70 &c>0.65

disp('The model is bad and try again')

end

for i=1:length(x0)

Hatx00(i)=Hatx0(i)

end

z=1:length(x0)

plot(z,x0,'-',z,Hatx00,':') %将原始数据和模拟值画在一个图上帮助观察

end

主要是用regress函数来进行孝晌：给你举个例子来说明吧。

x=[0 1 2 3 4 ]'y=[1.0 1.3 1.5,2.0 2.3]'

x=[ones(5,1),x]%给出两个数组元素

[b,bint,r,rint,stats]=regress(y,x,0.05) ％对x和y进行一元知做线性回归，并得到相关系数，其中，stats中第一个数搭慎衡即为相关系数，大于0.9就认为拟合很好。

结果：stats =

0.9829 171.94740.00100.0063

即为0.9829.

---