matlab、lingo程序代码11-灰色预测

灰色预测

%程序中的变量定义:alpha是包含α、μ值的矩阵;

%ago是预测后累加值矩阵;var是预测值矩阵;

%error是残差矩阵; c是后验差比值

function basicgrey(x,m) %定义函数basicgray(x)

if nargin==1 %m为想预测数据的个数，默认为1

m=1;

end

clc; %清屏,以使计算结果独立显示

if length(x(:,1))==1 %对输入矩阵进行判断,如不是一维列矩阵,进行转置变换x=x';

end

n=length(x); %取输入数据的样本量

x1(:,1)=cumsum(x); %计算累加值,并将值赋与矩阵be

for i=2:n %对原始数列平行移位

Y(i-1,:)=x(i,:);

end

for i=2:n %计算数据矩阵B的第一列数据

z(i,1)=0.5*x1(i-1,:)+0.5*x1(i,:);

end

B=ones(n-1,2); %构造数据矩阵B

B(:,1)=-z(2:n,1);

alpha=inv(B'*B)*B'*Y; %计算参数α、μ矩阵

for i=1:n+m %计算数据估计值的累加数列,如改n+1为n+m可预测后m个值ago(i,:)=(x1(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha(2,:)/alpha(1,:);

end

var(1,:)=ago(1,:);

for i=1:n+m-1 %可预测后m个值

var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原,并计算出下m个预测值end

[P,c,error]=lcheck(x,var); %进行后验差检验

[rela]=relations([x';var(1:n)']); %关联度检验

ago %显示输出预测值的累加数列

alpha %显示输出参数α、μ数列

var %显示输出预测值

error %显示输出误差

P %显示计算小残差概率

c %显示后验差的比值c

rela %显示关联度

judge(P,c,rela) %评价函数显示这个模型是否合格

clear

clc

X=[136 143 165 152 165 181 204 272 319 491 571 605 665 640 628]; x1(1)=X(1);

X1=[];

for i=1:1:14

x1(i+1)=x1(i)+X(i+1);

X1=[X1,x1(i)];

end

X1=[X1,X1(14)+X(15)]

for k=3:1:15

p(k)=X(k)/X1(k-1);

p1(k)=X1(k)/X1(k-1);

end

p,p1

clear k

Z=[];

for k=2:1:15

z(k)=0.5*X1(k)+0.5*X1(k-1);

Z=[Z,z(k)];

end

Z

B=[-Z',ones(14,1)]

Y=[];

clear i

for i=2:1:15

Y=[Y;X(i)];

end

Y

A=inv(B'*B)*B'*Y

clear k

y1=[];

for k=1:1:15

y(k)=(X(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1);

y1=[y1;y(k)];

end

y1

clear k

X2=[];

for k=2:1:15

x2(k)=y1(k)-y1(k-1);

X2=[X2;x2(k)];

end

X2=[y1(1);X2]

m=abs(e)./X'

s=e'*e

n=sum(m)/13

clear k

syms k

y=(X(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1)

Y1=[];

for j=16:1:21

y11=subs(y,k,j)-subs(y,k,j-1);

Y1=[Y1;y11];

end

Y1

function judge(P,c,rela)

%评价指标并显示比较结果

if rela>0.6

'根据经验关联度检验结果为满意（关联度只是参考主要看后验差的结果）'

else

'根据经验关联度检验结果为不满意（关联度只是参考主要看后验差的结果）'

end

if P>0.95&c<0.5

'后验差结果显示这个模型评价为“优”'

else if P>0.8&c<0.5

'后验差结果显示这个模型评价为“合格”'

else if P>0.7&c<0.65

'后验差结果显示这个模型评价为“勉强合格”'

else

'后验差结果显示这个模型评价为“不合格”'

end

end

end

function [P,c,error]=lcheck(x,var)

%进行后验差检验

n=length(x);

for i=1:n

error(i,:)=abs(var(i,:)-x(i,:)); %计算绝对残差

end

c=std(abs(error))/std(x); %调用统计工具箱的标准差函数计算后验差的比值c

s0=0.6745*std(x);

ek=abs(error-mean(error));

for i=1:n

if ek(i,:)

pk=pk+1;

end

end

P=pk/n; %计算小残差概率

%附带的质料里有一部分讲了关联度

function [rela]=relations(x)

%以x(1,:)的参考序列求关联度

[m,n]=size(x);

for i=1:m

for j=n:-1:2

x(i,j)=x(i,j)/x(i,1);

end

end

for i=2:m

x(i,:)=abs(x(i,:)-x(1,:)); %求序列差

end

c=x(2:m,:);

Max=max(max(c)); %求两极差

Min=min(min(c));

p=0.5; %p称为分辨率,0

for j=1:n

r(i,j)=(Min+p*Max)/(c(i,j)+p*Max); %计算关联系数

end

end

for i=1:m-1

rela(i)=sum(r(i,:))/n; %求关联度

end