matlab实践（七）：GM（1,1）灰色预测模型

最新推荐文章于 2024-02-05 08:21:52 发布

从零开始的奋豆

最新推荐文章于 2024-02-05 08:21:52 发布

阅读量2.1k

点赞数 6

分类专栏： MATLAB 文章标签： matlab 算法开发语言

本文链接：https://blog.csdn.net/m0_68926749/article/details/133362954

版权

MATLAB 专栏收录该内容

12 篇文章 9 订阅

订阅专栏

1.过程

第一步：设原始数据序列： $X^{(0)}=\{x^{(0)}(1),x^{(0)}(2),\cdots ,x^{(0)}(N)\}$
用AGO生成一阶累加生成模块得到： $X^{(1)}=\{x^{(1)}(1),x^{(1)}(2),\cdots ,x^{(1)}(N)\}$

第二步：构造累加矩阵 $B$ 与常数项向量 $Y_{N}$ ，即

$B=\begin{pmatrix} 0.5*(x^{(1)}(1)+x^{(1)}(2)) & 1\\ 0.5*(x^{(1)}(2)+x^{(1)}(3))& 1\\ \cdots &\cdots \\ 0.5*(x^{(1)}(N-1)+x^{(1)}(N)) &1 \end{pmatrix}$

$Y_N=\left[{x_1}^{(0)}(2),{x_1}^{(0)}(3),\cdots,{x_1}^{(0)}(N)\right]^T$

$\hat{a}=\begin{bmatrix} a\\ u \end{bmatrix}$

上述方程组中， $Y_{N}$ 和 $B$ 为已知量，a为待定参数。由于变量只有a和u二个，而方程个数却有N-1个，而N-1>2，故方程组无解。但可用最小二乘法得到最小二乘解。

第三步：用最小二乘法解灰参数

$\hat{a}=(B^{T}B)^{-1}B^{T}Y_{N}$

第四步：将灰参数代入时间参数

$\hat{x}^{(1)}(t+1)=(x^{(0)}(1)-\frac{u}{a})e^{-at}+\frac{u}{a}$

第五步：对 $\hat{x}^{(1)}$ 求导还原得到 $\hat{x}^{(0)}$

$X^{(0)}(t+1)=\hat{x}^{(1)}(t+1)-\hat{x}^{(1)}(t)$

第六步：计算拟合误差

$\xi (t)=x^{(0)}(t)-\hat{x}^{(0)}(t)$

第七步：灰色GM（1，N）模型的检验分为三个方面：关联度检验、残差检验、后验差检验。后验差检验是残差分析统计特性的检验，模型诊断及应用模型进行预报。后验差比值C：残差方差S1与数据方差S2之比，即有：
$C=\frac{S_{1}}{S_{2}}$

计算残差方差S1
$S_{1}^{2}=\sum_{t=1}^{m}(x^{(0)}(t)-\hat{x}^{(0)}(t))^2$

及数据方差S2
$S_{2}^{2}=\frac{1}{m-1}\sum_{t=1}^{m}(q^{(0)}(t)-\bar{q}^{(0)}(t))^2$

小误差概率

$P=\{|q^{(0)}(t)-\bar{q}^{(0)}(t)|<0.6745*S_{1}\}$

若对于给定的 $C_{0}$ ＞0，当 $C$ ＜ $C_{0}$ 时，称模型为均方差比合格模型；如对给定的 $P_{0}$ ＞0，当 $P$ ＞ $P_{0}$ 时，称模型为小残差概率合格模型。根据后验比C和小误差概率P对模型进行判断，预测模型的精度等级见表二所示：

2.代码实现

clc;clear;
format short
X=input('请输入原始数据:','s');
X=str2num(X);
[m1 m2]=size(X);
k0=input('请输入所要预测的阶数:');
%GM(1,1)模型
for i=1:m1
    n=i;
    x0=X(i,:);
    disp('1.原始数据:');
    Y='';
    for z=1:m2
        Y=strcat(Y,'(',num2str(x0(z)),')');    
    end
    disp(Y);
    
    % 1. 利用一次累加（1-AGO）生成新数列,m2是X的列数，triu返回上三角矩阵，x0是X的每行，
    E=triu(ones(m2));
    x1=x0*E;
    disp('2.一次累加（1-AGO）生成的数据:');
    Y='';
    for z=1:m2
        Y=strcat(Y,'(',num2str(x1(z)),')');    
    end
    disp(Y);
    
    % 2. 计算出发展系数a，灰作用量u
    b1=x1;b1(1)=[];
    b2=x1;b2(m2)=[];
    b=-0.5*(b1+b2);
    B=[b;ones(1,m2-1)];
    B=B'; 
    y0=x0;y0(1)=[];
    y0=y0'; 
    A=((inv(B'*B))*B')*y0; 
    a=A(1);u=A(2);
    disp(a);
    disp(u);
    % 3. 确立模型且求出模拟值
    u_a=u/a;
    format long
    syms t
    ppp=vpa((x0(1)-u_a)*exp(-t*a)+u_a);
    
    for k=0:m2+k0-1
        x2(k+1)=(x0(1)-u_a)*exp(-k*a)+u_a; 
    end 
    
    x3=x2;
    x3(m2+k0)=[];
    x4=[0 x3];   
    x5=x2-x4;
%     disp('3.一次累加（1-AGO）生成的数据的模拟值:')
%     Y='';
%     for z=1:m2+k0
%         Y=strcat(Y,'(',num2str(x2(z)),')');    
%     end
%     disp(Y);
%     disp('4.原始数据的模拟值:')
%     Y='';
%     for z=1:m2+k0
%         Y=strcat(Y,'(',num2str(x5(z)),')');    
%     end
%     disp(Y);
    
    % 4. 模型检验（算出的值到等级参照表中检查其精度等级）
    
  
    x6=x5(1:m2);     
    Q=x0-x6;
    s1=std(Q);
    s2=std(x0);
    C=s1/s2;
    w1=1:m2;
    w1=[ones(m2,1) w1'];
    w2=Q';
    [bb,bint,r1,rint,stats]=regress(w2,w1);
    
    rcoplot(r1,rint)
    C1=strcat('后验差比(均方差比值): C=',num2str(C));disp(C1);
    if  C<=0.35
        disp('  由于C<=0.35，则此模型精度等级为1级（好）。');
    else if  C<=0.5
           disp('  由于0.35<C<=0.5，则此模型精度等级为2级（合格）。');
        else if C<=0.65
                disp('  由于0.5<C<=0.65，则此模型精度等级为3级（勉强）。');
            else 
                 disp('  由于C>0.65，则此模型精度等级为4级（不合格）。');
            end
        end
    end
    
    %计算相对误差序列R
    R=Q./x0;
%     disp('相对误差序列:')
%     Y='';
%     for z=1:m2
%         Y=strcat(Y,'(',num2str(R(z)),')');    
%     end
%     disp(Y);
%     disp('  （1）如果相对误差为0.01则模型精度等级为1级（好）');
%     disp('  （2）如果相对误差为0.05则模型精度等级为2级（合格）');
%     disp('  （3）如果相对误差为0.10则模型精度等级为3级（勉强）');
%     disp('  （4）如果相对误差为0.20则模型精度等级为4级（不合格）');
   
    %计算小误差概率P
    Qmean=mean(Q);
    D=abs(Q-Qmean);
    p0=0.6745*s2;
    t=0; 
    for j=1:m2 
        if  D(j)<p0 
            t=t+1;
        end 
    end 
    P=t/m2;%计算小误差概率P
    P1=strcat('小误差概率: P=',num2str(P));disp(P1)
    if  P>=0.95
        disp('  由于P>=0.95，则此模型精度等级为1级（好）。');
    else if  P>=0.8
           disp('  由于0.80<=P<0.95，则此模型精度等级为2级（合格）。');
        else if P>=0.7
                disp('  由于0.70<=P<0.8，则此模型精度等级为3级（勉强）。');
            else 
                 disp('  由于P<0.70，则此模型精度等级为4级（不合格）。');
            end
        end
    end
    
    % 5.模型适用范围
    a1=strcat('发展系数: a=',num2str(a));disp(a1);
    u1=strcat('  灰作用量: u=',num2str(u));disp(u1);
    if  -a<0.3
        disp('  由于-a<0.3，则此模型适合用于中长期预测。');
    else if  -a<0.5
           disp('  由于0.3<-a<0.5，则此模型适合用于短期预测，中长期预测慎用。');
        else if -a<0.8
                disp('  由于0.5<-a<0.8，则此模型作短期预测应十分谨慎。');
            else if -a<1
                    disp('  由于0.8<-a<1.0，则应采用残差修正 GM(1,1)模型。');
                end
            end
        end
    end
end
x5=x5'
Q=Q'
R=R'
[mm1 mm2]=size(x5);
xx5=x5(m2+1:mm1)
tt=[x0',x6',Q,R]

3.结果

GM（1,1）模型为

从零开始的奋豆

关注

6
点赞
踩
16

收藏

觉得还不错? 一键收藏
打赏
1
评论
matlab实践（七）：GM（1,1）灰色预测模型

设原始数据序列：构造累加矩阵与常数项向量，即上述方程组中，和为已知量，a为待定参数。由于变量只有a和u二个，而方程个数却有N-1个，而N-1>2，故方程组无解。但可用最小二乘法得到最小二乘解。用最小二乘法解灰参数将灰参数代入时间参数对求导还原得到计算拟合误差。后验差检验是残差分析统计特性的检验，模型诊断及应用模型进行预报。后验差比值C：残差方差S1与数据方差S2之比，即有：计算残差方差S1及数据方差S2小误差概率若对于给定的＞0，当。
复制链接

扫一扫