gm11_1949_2013_2018.zip_1_GM（1）_forecasting_产量预测_粮食预测

clc;clear; %建立符号变量a(发展系数)和b(灰作用量) syms a b; c = [a b]'; %原始数列 A A = [11318.40,13212.90,14368.90,16392.50,16684.10,16952.80,18394.60,19275.60,19504.50,19766.30,16969.20,14385.70,13650.90,15441.40,16574.10,18088.70,19452.50,21400.90,21782.30,20906.00,21097.30,23995.50,25014.00,24048.00,26493.50,27527.00,28451.50,28630.50,28272.50,30476.50,33211.50,32055.50,32502.00,35450.00,38727.50,40730.50,37910.80,39151.20,40297.70,39408.10,40754.90,44624.30,43529.30,44265.80,45648.80,44510.10,46661.80,50453.50,49417.10,51229.53,50838.58,46217.52,45263.67,45705.75,43069.53,46946.95,48402.19,49804.23,50413.86,53434.29,53940.86,55911.31,58849.33,61222.62,63048.20]; 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+5) G(i) = F(i) - F(i-1); %得到预测出来的数据 end disp('预测数据为：'); G lengg=length(G) %模型检验 H = G(1:65); %计算残差序列 epsilon = A - H; Eg=abs(epsilon) %法一：相对残差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 = 1949:2013; t2 = 1949:2018; LLt1 = length(t1) plot(t1, A,'r-'); hold on; plot(t2, G, 'g-'); xlabel('年份'); ylabel('粮食产量/万吨'); legend('实际量','预测量'); title('粮食产量变化曲线'); grid on;