GM(1,1)模型

clc;
clear;


em=[31.5319,45.468,81.6114,112.0086,106.944,139.2239,165.291,103.4906,130.5011,113.8424,125.8621,94.3148,54.6247,80.1303,64.7934,38.012,79.4888,50.9654,22.6694,35.5949,37,39,47,45.5969,46.0514,46.3908,12.1464,11.382];


% %%%%%%%%%%%%%%%%%%%%%%%%%采用GM(1,1)模型预测图%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 x0_length=length(em);
%定义1-AGO数组
x_1=zeros(size(em));
%给数组Y赋值
Y=emi(2:x0_length);
%定义数组z
z_1=zeros(1,x0_length-1);
%预测数组
x_prediction=zeros(1,x0_length);
for i=1:x0_length
    %计算1-AGO数组
    x_1(i)=sum(emi(1:i));
    if i>1
        %给数组z赋值
        z_1(i-1)=-(x_1(i-1)+x_1(i))/2;
    end
end
%B矩阵赋值
B2=ones(1,x0_length-1);
B=cat(1,z_1,B2);
%进行最小二乘估计，计算参数a,b
a=inv(B*B')*B*Y';
%给第一个赋值
x_prediction(1)=em(1);
%计算估计值
for i=1:x0_length-1
    x_prediction(i+1)=(1-exp(a(1,1)))*(emi(1)-a(2,1)/a(1,1))*exp(-a(1,1)*i);
end
%绘制图像
step=1:1:x0_length;
[H]=figure(1);
set(H,'Position',[300 300 640 480])
% subplot(1,3,1)
plot(step,emi)
hold on
plot(step,x_prediction,':b');
xlabel('时间(s)')
ylabel('dBm')
% ylabel('测量值(V)')
% title('\fontsize {8}GM(1,1)')
[L]=legend('测量值','估计值');
set(L, 'Box', 'off')