灰色预测的matlab代码怎么写？

clc

clear all

% 本程序主要用来计算根据灰色理论建立的模型的预测值。

% 应用的数学模型是 GM(1,1)。

% 橘咐芦原始数据的处理方法是一次累加法。

y=[1662.87 2163.4 1965.35 2472.48 2900.66 简族3034.93 2755.5 3207 3462]%已知数据

n=length(y)

yy=ones(n,1)

yy(1)=y(1)

for i=2:n

yy(i)=yy(i-1)+y(i)

end

B=ones(n-1,2)

for i=1:(n-1)

B(i,1)=-(yy(i)+yy(i+1))/2

B(i,2)=1

end

BT=B'

for j=1:n-1

YN(j)=y(j+1)

end

YN=YN'

A=inv(BT*B)*BT*YN

a=A(1)

u=A(2)

t=u/a

t_test=4  %需要预测个数

i=1:t_test+n

yys(i+1)=(y(1)-t).*exp(-a.*i)+t

yys(1)=y(1)

for j=n+t_test:-1:2

ys(j)=yys(j)-yys(j-1)

end

x=1:n

xs=2:n+t_test

yn=ys(2:n+t_test)

plot(x,y,'^r',xs,yn,'*-b')

det=0

for i=2:n

det=det+abs(yn(i)-y(i))

end

det=det/(n-1)

disp(['百分绝对误差为：',num2str(det),'%'])

disp(['预测值为： ',num2str(ys(n+1:n+t_test))])

输出结果：

百分绝对误差为：圆带228.3113%

预测值为： 3710.152      3978.2142      4265.6442      4573.8413

这是我曾经写过的一个灰色预测的程序：第一个文件为函数，需要在调用时输入原始数据x0和预测周期T, 第二个文件用于计算灰色关联度，使用时直接修改相关参数和原始数据。

--------------------------------------------------------------------------

第一个文件（用于灰色建模）：grymdl.m

--------------------------------------------------------------------------

function GM=grymdl(x0,T)

% 输入原始数据x0

% T为从最后一个历史数据算起的第T时点

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

yn=zeros(length(x0)-1,1)Hatx0=zeros(1,length(x0)+T)

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)模型参数估计

for k=1:length(x0)+T

    Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1)

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

x0

HatA

Hatx0

epsilon

omega

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

end

p=p/length(x0)

if p>=0.95

    M1=1

elseif p>=0.8

    M1=2

elseif p>=0.7

    M1=3

else

    M1=4

end

if c<=0.35

    M2=1

elseif c<=0.5

    M2=2

elseif c<=0.65

    M2=3

else

    M2=4

end

M=max(M1,M2)

if M==1

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

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

elseif M==2

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

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

elseif M==3

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

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

else

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

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

end

for i=1:length(x0)

  袜凯  Hatx00(i)=Hatx0(i)

end

z=1:length(x0)

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

text(2,x0(2),'History data: real line')

text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')

GM=Hatx0(length(x0)+T)

--------------------------------------------------------------------------------

第二个文件（用于计算灰色关联度）：grydgr.m

--------------------------------------------------------------------------------

x=[26.4,40.83,5,18.5,27,4.9

    32.56,50.36,6.5,23.8,34.9,5.97

    41.16,63.66,7.9,30.6,45,6.83

    51.82,80.15,9.6,39.43,60.76,8.21

    70.58,109.16,11.3,52.1,81.8,9.78

    111.47,136.12,15,73.92,118.12,12.3

    127.16,166.7,18,118.17,184.51,16.75

    165.99,196.91,20.66,156.65,240.83,21.49

    246.76,257.56,24.15,218.3,333.4,27.1] %原始数据

delta=zeros(size(x,1),size(x,2)-1) %初始化绝对差

yita=zeros(size(delta,1),size(delta,2)) %初始化关联系数

for i=1:size(x,2)

    x(:,i)=x(:,i)./x(1,i)  %无量纲化处理

end

for i=1:(size(x,2)-1)

    delta(:,i)=abs(x(:,1)-x(:,i+1))  %求解delta

end

delta_min=min(min(delta,[],1)) %求解最小二级差

delta_max=max(max(delta,[],1)) %求解最大二级差

rou=0.5 %设定分辨系数为0.5

for i=1:size(delta,1)

    for j=1:size(delta,2)

        yita(i,j)=(delta_min+rou*delta_max)/(delta(i,j)+rou*delta_max) %计算关联系数

    end

end

r=sum(yita)./size(yita,1) %计算灰色关联

2014.4.11~4.20 预测价格为 5.9535

>>%灰色预测模型键搏程序

clear

syms a b

c=[a b]'

A=[5.76 5.92 5.91 5.92 5.93 5.93 5.93 5.94 5.94 5.95 5.95 5.94 ]% 原始序列

B=cumsum(A)%累加

n=length(A)

for i=1:(n-1)

C(i)=(B(i)+B(i+1))/2

end

%计算待定参数

D=A

D(1)=[]

D=D'

E=[-Cones(1,n-1)]

c=inv(E*E')*E*D

c=c'

a=c(1)

b=c(2)

%预测 往后预测1个数据（搏亮胡即4.11-4.20价格）

F=[]F(1)=A(1)

for i=2:(n+1)

F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a

end

G=[]G(1)=A(1)

for i=2:(n+1)

G(i)=F(i)-F(i-1)

end

t1=2002:2013

t2=2002:2014

G

plot(t1,A,'基拦o',t2,G)

G =

5.76005.91555.91895.92245.92585.92935.93275.93625.93965.94315.94655.95005.9535

---