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=B\Y;
% 在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('原始數列','模型數列');

運行結果，程序輸出的圖形以下所示：

版權聲明：本文爲博主原創文章，未經博主容許不得轉載。