matlab之灰度预测

概述

这次来学习一下比较简单的灰度预测吧，学习数学建模的道路还真是挺艰难的

1.

clear
syms a b;
c=[a b]';
A=[89677,99215,109655,120333,135823,159878,182321,209407,246619,300670];
B=cumsum(A);  % 原始数据累加
n=length(A);
for i=1:(n-1)
    C(i)=(B(i)+B(i+1))/2;  % 生成累加矩阵
end
% 计算待定参数的值
D=A;D(1)=[];
D=D';
E=[-C;ones(1,n-1)];
c=inv(E*E')*E*D;
c=c';
a=c(1);b=c(2);
% 预测后续数据
F=[];F(1)=A(1);
for i=2:(n+10)
    F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a;
end
G=[];G(1)=A(1);
for i=2:(n+10)
    G(i)=F(i)-F(i-1); %得到预测出来的数据
end 
t1=1999:2008;
t2=1999:2018;
G
plot(t1,A,'ko', 'LineWidth',2)
hold on
plot(t2,G,'k', 'LineWidth',2)
xlabel('年份', 'fontsize',12)
ylabel('利润/(元/年)','fontsize',12)
set(gca,  'LineWidth',2);

2.

%建立符号变量a(发展系数)和b(灰作用量)
syms a b;
c = [a b]';

%原始数列 A
A = [174, 179, 183, 189, 207, 234, 220.5, 256, 270, 285];
n = length(A);

%对原始数列 A 做累加得到数列 B
B = cumsum(A);

%对数列 B 做紧邻均值生成
for i = 2:n
    C(i) = (B(i) + B(i - 1))/2; 
end
C(1) = [];

%构造数据矩阵 
B = [-C;ones(1,n-1)];
Y = A; Y(1) = []; Y = Y';

%使用最小二乘法计算参数 a(发展系数)和b(灰作用量)
c = inv(B*B')*B*Y;
c = c';
a = c(1); b = c(2);

%预测后续数据
F = []; F(1) = A(1);
for i = 2:(n+10)
    F(i) = (A(1)-b/a)/exp(a*(i-1))+ b/a;
end

%对数列 F 累减还原,得到预测出的数据
G = []; G(1) = A(1);
for i = 2:(n+10)
    G(i) = F(i) - F(i-1); %得到预测出来的数据
end

disp('预测数据为：');
G

%模型检验

H = G(1:10);
%计算残差序列
epsilon = A - H;

%法一：相对残差Q检验
%计算相对误差序列
delta = abs(epsilon./A);
%计算相对误差Q
disp('相对残差Q检验：')
Q = mean(delta)

%法二：方差比C检验
disp('方差比C检验：')
C = std(epsilon, 1)/std(A, 1)

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

%绘制曲线图
t1 = 1995:2004;
t2 = 1995:2014;

plot(t1, A,'ro'); hold on;
plot(t2, G, 'g-');
xlabel('年份'); ylabel('污水量/亿吨');
legend('实际污水排放量','预测污水排放量');
title('长江污水排放量增长曲线');
grid on;

3.

x0 = [71.1 72.4 72.4 72.1 71.4 72 71.6]'; %这里是列向量，相当于原始数据中因变量
n = length(x0);
lamda = x0(1:n-1)./x0(2:n) %计算级比
range = minmax(lamda') %计算级比的范围
x1 = cumsum(x0)
B = [-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)]; %这是构造的数据矩阵B
Y = x0(2:n); %数据向量Y
u = B\Y  %拟合参数u(1)=a,u(2)=b
syms x(t)
x = dsolve(diff(x)+u(1)*x==u(2),x(0)==x0(1)); %建立模型求解
xt = vpa(x,6) %以小数格式显示微分方程的解
prediction1 = subs(x,t,[0:n-1]); %求已知数据的预测值
prediction1 = double(prediction1); %符号数转换成数值类型，以便做差分运算
prediction = [x0(1),diff(prediction1)] %差分运算，还原数据
epsilon = x0'-prediction %计算残差
delta = abs(epsilon./x0') %计算相对残差
rho = 1-(1-0.5*u(1))/(1+0.5*u(1))*lamda'%计算级比偏差值，u(1)=a

4.

5.

参考链接：