GM(1,1)程序设计

GM(1,1)是灰色模型中较为常见的模型，下面是程序，x0是数据，可更改。

之前编辑忘了说了，一般就是给定一组数据，自己根据这些数据拟合一个灰色模型，底下的代码可以得到该模型对应的公式。

代码：

% GM(1,1)
% 程序有详尽注释
clc;
clear all;
x0=[92.810 97.660 98.800 99.281 99.537 99.537 99.817 100.000];

n=length(x0);
% 做级比判断，看看是否适合用GM(1,1)建模
lamda=x0(1:n-1)./x0(2:n);
range=minmax(lamda);
% 判定是否适合用一阶灰色模型建模
if range(1,1)<exp(-(2/(n+2)))| range(1,2)>exp(2/(n+2))
    error('级比没有落入灰色模型的范围内，不适合用GM(1,1)建模');
else
    % 空行输出
    disp('            ');
    disp('可用GM(1,1)建模');
end

%做AGO累加处理
x1=cumsum(x0);
for i=2:n
    %计算紧邻均值，也就是白化背景值
    z(i)=0.5*(x1(i)+x1(i-1));
end
B=[-z(2:n)',ones(n-1,1)];
Y=x0(2:n)';
% 矩阵做除法，计算发展系数和灰色作用量
% 注意：这里是右除不是左除
u=BY;
% 在MATLAB中，用大写字母D表示导数,dsolve函数用来求解符号常微分方程
x=dsolve('Dx+a*x=b','x(0)=x0');
% subs函数的作用是替换元素，这里是把常量a，b，x0换成具体u(1),u(2),x(1)数值
x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)});
disp('函数:');
vpa(x,6)
forecast1=subs(x,'t',[0:n-1]);
% digits和vpa函数用来控制计算的有效数字的位数
digits(6)
% y值是AGO形式的(还是累加的)
y=vpa(x);
% 把AGO输出值进行累减
% diff用于对符号表达时进行求导
% 但是如果输入的是向量，则表示计算原向量相邻元素的差
exchange=diff(forecast1);
% 输出灰色模型预测的值
disp('输出预测模型预测的值：');
forecast=[x0(1),exchange]
% 计算残差
 epsilon=x0-forecast;
% 计算相对误差
delta=abs(epsilon./x0);

% 检验模型的误差
% 检验方法一：相对误差Q检验法
disp('相对误差Q检验值：');
Q=mean(delta)

% 检验方法二：方差比C检验法
% 计算标准差函数为std(x,a)
% 如果后面一个参数a取0表示的是除以n-1,如果是1就是最后除以n
disp('方差比C检验值：');
C=std(epsilon,1)/std(x0,1)

% 检验方法三：小误差概率P检验法
S1=std(x0,1);
S1_new=S1*0.6745;
temp_P=find(abs(epsilon-mean(epsilon))<S1_new);
disp('小误差概率P检验值：');
P=length(temp_P)/n

%绘制原始数列和灰色模型预测得出的数列差异折线图
plot(1:n,x0,'ro','markersize',11);
hold on
plot(1:n,forecast,'k-','LineWidth',2.5);
grid on;
axis tight;
xlabel('x');
ylabel('y');
title('保有量比例与时间序列的时间关系');
legend('原始数列','模型数列');

运行结果，程序输出的图形如下所示：