gm11.m

灰色预测模型称为CM模型

%
% IIi为预测点
% x为预测序列（一行）
% X0X为预测值（对有IIi而言）
% Q0Q为预测精度
% pp为小误差概率
% C为后验误差比
% pp>0.95  或 C<0.35  为 1 级(好)
% pp>0.80  或 C<0.50  为 2 级(合格)
% pp>0.70  或 C<0.65  为 3 级(勉强合格)
% pp<=0.70 或 C>=0.65 为 4 级(不合格)
%
% AUA=[a;u;ua;xua];
% XEG=[xg;eg];
% PPC=[pp;C];
function [X0X,Q0Q,PPC,pddstr,XIi,AUA,XEG]=GM11(x,IIi);
%%%%%%%%%%%%%%%%%%%%%%%

%例子
%x=[3.711 3.723 3.716 3.721 3.728]
%x=[2.67 3.13 3.25 3.36 3.56 3.72]
% x=[2.97 3.23 3.29 3.46 3.59 3.71]
%x=[43.45 47.05 52.75  57.14 62.64 68.52]
% x=[3.38 4.27 4.55 4.69 5.59]
% x=[4.24 4.33 5.20 6.42 7.32 8.53 8.82 10.72]
% 
% x=[3.38 4.27 4.55 4.69 5.59]
% IIi=[  5 6 7 8 9 10 11]
% IIi=[]

%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                %%
%%        一、数据处理            %%
%%                                %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%1建立X1生成数列
format long
a12=size(x);
n=a12(2);
x1=ones(1,n);
x1(1)=x(1);
for i=2:n
    x1(i)=x1(i-1)+x(i);
end
x1;
%2构造数据距阵B和数据向量Yn:
B=ones(n-1,2);
for i=1:n-1
    B(i,1)=(-0.5)*(x1(i)+x1(i+1));
end
Yn=ones(n-1,1);
for i=2:n
    Yn(i-1,1)=x(i);
end
B;
Yn;
%3计算B'B,inv(B'B)和B'Yn；得到参数a和u
%BTB=B'*B;%
%BTb=inv(BTB);
%BTYn=B'*Yn;
%au=BTb*BTYn;
au=inv(B'*B)*B'*Yn;


a=au(1,1);
u=au(2,1);

% fprintf('参数a：%f\n',a);
% fprintf('参数u：%f\n',u);
ua=u/a;
xua=x(1)-ua;
%4得出预测模型
%预测第I的预测值
%有问题

%  XI1=xua*exp((-a)*I)+ua;%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                %%
%%        二、模型检验            %%
%%                                %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%5:残差检验
%51计算X1
X1=ones(1,n);
for i=1:n
    X1(i)=xua*exp((-a)*(i-1))+ua;
end
X1;
%52累减生成X0序列
X0=ones(1,n);
X0(1)=X1(1);
for i=2:n
    X0(i)=X1(i)-X1(i-1);
end
X0;
%53计算绝对误差及相对误差序列
Dd=abs(X0-x);%计算绝对误差序列
DDd=(Dd./x).*100;%计算相对误差序列 x%
% fprintf('相对误差序列：');
% fprintf('%f%%  ',DDd);
% fprintf('\n');
%6进行关联度检验
%61计算序列绝对误差Xx
Xx=ones(1,n);
for i=1:n
    Xx(i)=x(i)-X0(i);
end
% fprintf('序列绝对误差：');
% fprintf('%f ',Xx);
% fprintf('\n');

minXx=min(Dd);
% fprintf('最小差：');
% fprintf('%f ',minXx);
% fprintf('\n');

maxXx=max(Dd);
% fprintf('最大差：');
% fprintf('%f ',maxXx);
% fprintf('\n');
%62计算关联系数P=0.5 (则r=0.6)
P=0.5;
nI=ones(1,n);
for i=1:n
    nI(i)=(minXx+P*maxXx)/(Dd(i)+P*maxXx);
end
%63计算关联度,(P=0.5时，则r=0.6)
r=(1/n)*sum(nI);
% fprintf('请查找P=0.5是的检验准则r=0.6是否大于%f\n',r);
%7后验差检验
%71
xm=mean(x);
%72求的均方差
s1=(sum((x-xm).^2)/(n-1))^(1/2);
%73计算残差的均值
Ddm=mean(Dd);
%74计算残差的均方差
s2=(sum((Dd-Ddm).^2)/(n-1))^(1/2);
%75计算后验误差比C:
C=s2/s1;
% fprintf('验误差比C：%f\n',C);
%76计算小误差概率
%%%%%%%%%%%%%%

pr=abs(Xx-mean(Xx))<0.6745*s1; %m满足条件的样本
pa=size(find(pr==1));
ps=pa(1,2); %m满足条件的样本个数
pb=size(pr);
pS=pb(1,2); %m总样本个数

pp=ps/pS; %小误差概率
% fprintf('小误差概率pp：%f\n',pp);

%%%%%%%%%%%%%%%%
%检验  预测精度ycd

% fprintf('   pp>0.95  或 C<0.35  为 1 级(好)\n   pp>0.80  或 C<0.50  为 2 级(合格)\n   pp>0.70  或 C<0.65  为 3 级(勉强合格)\n   pp<=0.70 或 C>=0.65 为 4 级(不合格)\n\n',C);
if pp>0.95  
    pd=1;
%     fprintf('   因 pp>0.95 且 ');
end
if pp<=0.95&pp>0.80  
    pd=2;
%     fprintf('   因 pp<=0.95&pp>0.80 且 '); 
end
if pp<=0.80&pp>0.70  
    pd=3;
%     fprintf('   因 pp<=0.80&pp>0.70 且 ');
end
if pp<=0.70  
    pd=4; 
%     fprintf('   因 pp<=0.70 且 '); 
end

pd1=0;
if C<0.35  
    pd1=1; 
%     fprintf('C<0.35'); 
end
if C>=0.35&C<0.50 
    pd1=2;
%     fprintf('C>=0.35&C<0.50'); 
end
if C>=0.50&C<0.65 
    pd1=3;
%     fprintf('C>=0.50&C<0.65');
end
if C>=0.65 
    pd1=4; 
%     fprintf('C>=0.65');
end

pdd=max(pd,pd1);

if pdd==1 
    pddstr='1 级（好）';       
%     fprintf(' 故根据经验，预测精度为 1 级（好）\n\n'); 
end
if pdd==2
    pddstr='2 级（合格）' ;   
%     fprintf(' 故根据经验，预测精度为 2 级（合格）\n\n'); 
end
if pdd==3
    pddstr='3 级（勉强合格）'; 
%     fprintf(' 故根据经验，预测精度为 3 级（勉强合格）\n\n'); 
end
if pdd==4
    pddstr='4 级（不合格）';  
%     fprintf(' 故根据经验，预测精度为 4 级（不合格）\n\n');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                %%
%%          三、预测              %%
%%                                %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%8模型经验合格后可用于预测，预测公式为：X0(i+1)=Xx(i+1)-Xx(i)
a1=size(IIi);
n1=a1(2);
X0X=ones(1,n1);
for i=1:n1
    X0X(i)=xua*(exp((-a)*(IIi(i)))-exp((-a)*(IIi(i)-1)));
%    fprintf('预测结果：第%d个值(原数据的第%d个序列号)的预测结果为：%f\n',IIi(i),IIi(i)+1,X0X(i));
end
px=[1:n];
% plot(px,x,'-k',IIi+1,X0X,'-*r');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                %%
%%       三、预测值精度评估       %%
%%                                %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算出的模型值为
xg=ones(1,n);
xg(1)=x(1);
for i=2:n
    xg(i)=xua*(exp((-a)*(i-1))-exp((-a)*(i-2)));
end
% fprintf('计算出的模型值为：');
% fprintf('%f ',xg);
% fprintf('\n');

%残差分别为
eg=x-xg;
% fprintf('残差分别为：');
% fprintf('%f ',eg);
% fprintf('\n');

Q=inv(B'*B);
Q11=Q(1,1);
Q12=Q(1,2);
Q21=Q(1,2);
Q22=Q(2,2);

q0=((eg*eg')/(n-1))^0.5;

q0q=ones(1,n);
for k=0:n-1
    q0q(k+1)=((a*k*x(1)-x(1)-k*u)^2*Q11+Q22+2*(a*k*x(1)-x(1)-k*u)*Q12)^0.5*exp(-a*k)*q0;
end
% fprintf('\n');
% for i=1:n
%     fprintf('原数据的第%d个序列号的预测值为：%f±%f\n',i,x(i),q0q(i));
% end
% fprintf('\n');
%========================
AUA=[a;u;ua;xua];
XEG=[xg;eg];
PPC=[pp;C];
%========================
%所以预测值为
XIi=0;
if isempty(IIi)
    XIi='-';
    X0X='-';
    Q0Q='-';
    return;        
end
fprintf('');
Q0Q=ones(1,n1);
for k=1:n1   
    Q0Q(k)=((a*IIi(k)*x(1)-x(1)-IIi(k)*u)^2*Q11+Q22+2*(a*IIi(k)*x(1)-x(1)-IIi(k)*u)*Q12)^0.5*exp(-a*IIi(k))*q0;
%     fprintf('预测结果：第%d个值(原数据的第%d个序列号)的预测值为：%f±%f\n',IIi(k),IIi(k)+1,X0X(k),Q0Q(k));
end