数学建模常用Matlab/Lingo/c代码总结系列——灰色预测

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]
e=X'-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


%程序中的变量定义:alpha是包含α、μ值的矩阵;%ago是预测后累加值矩阵;var是预测值矩阵;%error是残差矩阵; c是后验差比值function basicgrey(x,m) %定义函数basicgray(x)if nargin==1 %m为想预测数据的个数，默认为1 m=1;endclc; %清屏,以使计算结果独立显示if length(x(:,1))==1 %对输入矩阵进行判断,如不是一维列矩阵,进行转置变换 x=x';endn=length(x); %取输入数据的样本量x1(:,1)=cumsum(x);
 %计算累加值,并将值赋与矩阵be for i=2:n %对原始数列平行移位 Y(i-1,:)=x(i,:);endfor i=2:n %计算数据矩阵B的第一列数据 z(i,1)=0.5*x1(i-1,:)+0.5*x1(i,:);endB=ones(n-1,2); %构造数据矩阵BB(:,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,:);endvar(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 %显示后验差的比值crela %显示关联度judge(P,c,rela)
 %评价函数 显示这个模型是否合格

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));
pk=0;
for i=1:n
    if ek(i,:)<s0
        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<p<1,一般取p=0.5   
for i=1:m-1
    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