《灰色預測(GM)的MATLAB實現》
一、 灰色模型GM(1,1)
1. 問題
請以下表的資料為依據,預測2005-2014年長江的汙水排放量(單位:億噸)。
1995-2004年的長江汙水排放量
年份 | 1995 | 1996 | 1997 | 1998 | 1999 | 2000 | 2001 | 2002 | 2003 | 2004 |
---|---|---|---|---|---|---|---|---|---|---|
汙水量/億噸 | 174 | 179 | 183 | 189 | 207 | 234 | 220.5 | 256 | 270 | 285 |
2. 分析
此問題為一個複雜的非線性系統,樣本資料量少,但需要預測的時間較長,且汙水排放量的變化規律是一個不確定的系統。如果使用神經網路演算法很難取得理想的效果,故考慮採用GM預測來預測未來的汙水排放量。
3. MATLAB實現原始碼
GM(1,1).m
%建立符號變數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;
執行結果:
預測資料為:
G =
1 至 14 列
174.0000 172.8090 183.9355 195.7785 208.3839 221.8010 236.0820 251.2825 267.4616 284.6825 303.0122 322.5221 343.2881 365.3912
15 至 20 列
388.9175 413.9585 440.6118 468.9812 499.1772 531.3174
相對殘差Q檢驗:
Q =
0.0234
方差比C檢驗:
C =
0.1870
小誤差概率P檢驗:
P =
1
4. MATLAB繪製的曲線圖
二、 灰色Verhulst模型(即Logistic模型)
1. 問題
將一定量的大腸桿菌菌種接種在液體培養基中,在一定條件下進行培養,觀察其生長繁殖規律。細菌懸液的濃度與混濁度成正比,故可用分光亮度計測定細菌懸液的光密度來推知菌液的濃度。每隔5h記錄OD600的值,得到下表。請你預測大腸桿菌的數量。
時間點均勻取樣/5h | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
細菌培養液吸光度/OD600 | 0.025 | 0.023 | 0.029 | 0.044 | 0.084 | 0.164 | 0.332 | 0.521 | 0.97 | 1.6 | 2.45 | 3.11 | 3.57 | 3.76 | 3.96 | 4 | 4.46 | 4.4 | 4.49 | 4.76 | 5.01 |
2. 分析
此問題涉及生物的生長和繁殖規律,其曲線一般呈S型或變異S型,故考慮使用GM Verhulst模型來預測。
3. MATLAB實現原始碼
(以下程式預測到時間點41)
GM_Verhulst.m
%建立符號變數a(發展係數)和b(灰作用量)
syms a b;
c = [a b]';
%原始數列 A
A = [0.025, 0.023, 0.029, 0.044, 0.084, 0.164, 0.332, 0.521, 0.97, 1.6, 2.45, 3.11, 3.57, 3.76, 3.96, 4, 4.46, 4.4, 4.49, 4.76, 5.01];
n = length(A);
%對原始數列 A 做累減得到數列 B
for i = 2:n
H(i) = A(i) - A(i - 1);
end
H(1) = [];
%對原始數列 A 做緊鄰均值生成
for i = 2:n
C(i) = (A(i) + A(i-1))/2;
end
C(1) = [];
%構造資料矩陣
D = [-C; C.^2];
Y = H; Y = Y';
%使用最小二乘法計算引數 a(發展係數)和b(灰作用量)
c = inv(D*D')*D*Y;
c = c';
a = c(1); b = c(2);
%得到預測出的資料
F = []; F(1) = A(1);
for i = 2:(n+n)
F(i) = (a*A(1))/(b*A(1)+(a - b*A(1))*exp(a*(i-1)));
end
disp('預測資料為:');
F
%繪製曲線圖
t1 = 0:n-1;
t2 = 0:2*n-1;
plot(t1, A, 'ro'); hold on;
plot(t2, F);
xlabel('時間點均勻取樣/5h'); ylabel('細菌培養液吸光度/OD600');
legend('實際數量','預測數量');
title('大腸桿菌培養S形增長曲線');
grid on;
執行結果:
預測資料為:
F =
1 至 14 列
0.0250 0.0416 0.0691 0.1143 0.1880 0.3059 0.4900 0.7658 1.1551 1.6603 2.2492 2.8555 3.4051 3.8485
15 至 28 列
4.1738 4.3963 4.5412 4.6326 4.6891 4.7236 4.7445 4.7571 4.7647 4.7692 4.7720 4.7736 4.7746 4.7752
29 至 42 列
4.7755 4.7757 4.7759 4.7759 4.7760 4.7760 4.7760 4.7760 4.7760 4.7760 4.7760 4.7760 4.7760 4.7760