优化的灰色模型
优化的是μ值和C值,若有错误希望大家多多指出
MATLAB 算法实现
clear;clc;
ysh =[0.851,0.66733333,0.603,0.779,0.857666667,0.374666667,0.669,0.813333333,0.252,0.536333333,0.190666667]; %输入的数据
n = length(ysh);
AGO = zeros(1,n);
AGO(1,1) = ysh(1,1);
for k = 2:n
AGO(1,k) = AGO(1,k-1) + ysh(1,k);
end %一次累加序列
for i = 1:n
P(n) = ysh(n-1)/AGO(n);
if P(n) <0.5
continue
else
disp('不满足光滑性检验')
end
end %光滑性检验
for i = 1:n
O(n) = AGO(n-1)/AGO(n);
if 1<=O(n)<=1.5
continue
else
disp('不满足指数规律检验要求')
end
end %指数规律检验
for i = 2:n
S(n) = ysh(n)/ysh(n-1);
if exp(-2/(n+1))<=S(n)<=exp(-2/(n+1))
continue
else
disp('不满足级比检验要求')
end
end %级比检验
i = 0;
M = []; A = [];C = [];S = [];F = [];
for m = 0:0.001:1 %确定权重系数起始值、递增值和终止值
i = i + 1;
z = zeros(1,n-1); %计算背景值
for k = 1:n - 1
z(1,k) = m*AGO(1,k+1) + (1-m)*AGO(1,k); %邻值生成序列
end
B = zeros(n-1,2);
Y = zeros(n-1,1); %计算数据矩阵 B 和数据向量 Y
for k = 1:n-1
B(k,1) = -z(1,k);
B(k,2) = 1;
Y(k,1) = ysh(1,k+1);
end
N = zeros(2,1);
N = inv(B'*B)*B'*Y;
a = N(1,1);
u = N(2,1); %计算a,u值
c1 = 0;c2 = 0;
for k = 2:n
c1 = c1 + ysh(1,k)*exp(-a*(k-1));
c2 = c2 + exp(-2*a*(k-1));
end
c = ((ysh(1,1)-u/a)/(1-exp(a)) + c1)/((1-exp(a))^(-2)+c2); %计算改进模型参数 C
ycAGO = zeros(1,n);
for k = 0:n-1
ycAGO(1,k+1) = c*(1-exp(a))^(-1)*exp(-a*k)+u/a; %计算预测值一次累加生成(1AGO)
end
yc = zeros(1,n);
for k = 1:n-1
yc(1,1) = ycAGO(1,1);
yc(1,k+1) = c * exp(-a*k); %yc即预测值
end
s = sum((ysh - yc).^2); %s即离差平方和
M(1,i) = m;
A(1,i) = a;
U(1,i) = u;
C(1,i) = c;
S(1,i) = s;
F(1,i) = m;
F(2,i) = s;
end
Y = sortrows(F',2) %依据残差大小进行排序
yc
plot(M,S) %绘制离方平均和与权重系数的关系图
结果
由结果可知,当权重系数为0.500的时候,离差平方和最小,所以选择权重系数为0.500(下面的离差平方和相等的权重系数也可以)
## 若不满足级比检验的要求,则应该对原始数据进行适当的平移:Y(0k) = Y(0k) + C,直到满足要求为止。`clear;clc;
ysh = [0.5997,0.597,0.559,0.542,0.509,0.437,0.409,0.395,0.391,0.393333333,0.371666667]; %输入的数据
n = length(ysh);
AGO = zeros(1,n);
AGO(1,1) = ysh(1,1);
for k = 2:n
AGO(1,k) = AGO(1,k-1) + ysh(1,k);
end
for j = -10:0.01:10
aysh(n) = ysh(n) + j;
aS(n) = aysh(n)/aysh(n-1);
for i = 2:n
if exp(-2/(n+1))<=aS(n)<=exp(-2/(n+1))
i
j
disp(‘满足’)
else
j;
end
end
end
从上述代码的结果中找到适当的平移数C##
选择权重系数为0.500,进行灰色模型的构建,具体代码如下所示:
clear;clc;
ysh = [0.851,0.66733333,0.603,0.779,0.857666667,0.374666667,0.669,0.813333333,0.252,0.536333333,0.190666667];
n = length(ysh);
AGO = zeros(1,n);
AGO(1,1) = ysh(1,1);
for k = 2:n
AGO(1,k) = AGO(1,k-1) + ysh(1,k);
end
z = zeros(1,n-1);
for k = 1:n - 1
z(1,k) = 0.5*AGO(1,k+1) + 0.5*AGO(1,k); %权值选择为0.5
end
B = zeros(n-1,2);
Y = zeros(n-1,1);
for k = 1:n-1
B(k,1) = -z(1,k);
B(k,2) = 1;
Y(k,1) = ysh(1,k+1);
end
N = zeros(2,1);
N = inv(B'*B)*B'*Y;
a = N(1,1);
u = N(2,1);
c1 = 0;c2 = 0;
for k = 2:n
c1 = c1 + ysh(1,k)*exp(-a*(k-1));
c2 = c2 + exp(-2*a*(k-1));
end
c = ((ysh(1,1)-u/a)/(1-exp(a)) + c1)/((1-exp(a))^(-2)+c2);
ycAGO = zeros(1,n);
for k = 0:n-1
ycAGO(1,k+1) = c*(1-exp(a))^(-1)*exp(-a*k)+u/a;
end
yc = zeros(1,n);
for k = 1:n-1
yc(1,1) = ycAGO(1,1);
yc(1,k+1) = c * exp(-a*k);
end
s = sum((ysh - yc).^2);
yc %预测值
H = yc(1:11);
%计算残差序列
epsilon = ysh - yc;
%法一:相对残差Q检验
%计算相对误差序列
delta = abs(epsilon./ysh);
%计算相对误差Q
disp('相对残差Q检验:')
Q = mean(delta)
%法二:方差比C检验
disp('方差比C检验:')
C = std(epsilon, 1)/std(ysh, 1)
%法三:小误差概率P检验
S1 = std(ysh, 1);
tmp = find(abs(epsilon - mean(epsilon))< 0.6745 * S1);
disp('小误差概率P检验:')
P = length(tmp)/n
%绘制曲线图
t1 = 1:11;
t2 = 1:11;
plot(t1, ysh,'ro'); hold on;
plot(t2, yc, ':b*');
xlabel('xx'); ylabel('yy');
legend('实际值','预测值');
grid on;
结果如下图所示
可以根据三个检验的值来判断模型的精度,具体如下图所示:
参考文献> 1.杨旭. 改进的灰色预测GM(1,1)模型的MATLAB实现[J]. 江苏科技信息, 2014(7):69-70.
2.组合灰色模型在黑土区玉米产量预测中的应用