R语言 ma模型预测 r语言gmm估计

文章目录

• GM(1,1)模型

• 1 GM(1,1)模型概述
• 2 GM(1,1)数据处理方法
• 3 GM(1,1)模型方法的可行性检验
• 4 GM(1,1)预测模型的构建
• 5 GM(1,1) 模型的检验
• 6 GM(1,1) 模型的适用范围
• 7 GM(1,1) 残差模型

• R语言实现
• python

GM(1,1)模型

1 GM(1,1)模型概述

灰色预测经常用来解决数据量较少且不能直接发现规律的数据。对于包含不确定信息的序列，灰色预测方法通过对原始数据进行处理，使之转化为灰色序列，并建立微分方程模型模型

GM(1,1)模型是灰色预测理论的基本模型，也是灰色系统理论中运用最广泛的一种动态预测模型，模型由一个单变量的一阶微分方程构成。GM(1,1)模型适合对“少数据，贫信息”的对象进行中短期预测

2 GM(1,1)数据处理方法

在灰色系统中，能获得的信息非常有限，且不易发现内部规律，因此我们需要一些处理，减弱序列的随机性，使得一般规律可以显现。处理方式主要有以下三种：累加生成、累减生成和加权邻值生成

3 GM(1,1)模型方法的可行性检验

（1）级比检验

由原始序列
$X^{\left( 0 \right)}=\left( x^{\left( 0 \right)}\left( 1 \right) ,x^0\left( 2 \right) ,...,x^0\left( n \right) \right)$

设 $\lambda \left( k \right)$ 为序列$X^{\left( 0 \right)}$

$\lambda \left( k \right) =\frac{x^{\left( 0 \right)}\left( k-1 \right)}{x^{\left( 0 \right)}\left( k \right)}\text{，}k=2,3,\cdots ,n$若满足：

$\lambda \left( k \right) \in \left( e^{\frac{-2}{\text{n}+1}},e^{\frac{2}{n+1}} \right)$
其中 n 为原始序列的数据个数，则$X^{\left( 0 \right)}$

否则，需对数据进行处理，使其落在有效范围内，如平移变换：

$y^{\left( 0 \right)}\left( k \right) =x^{\left( 0 \right)}\left( k \right) +c\text{，}k=1,2,\cdots ,n$

---

（2）光滑度检验

设数据序列 $\left\{ x^{\left( 0 \right)}\left( k \right) \text{，}k=1,2,\cdots ,n \right\}$ 为非负数据序列，则该数据序列为光滑离散序列的充要条件为 $p(k)$ 是 $k$ 的递减函数，即
$p\left( k \right) =\frac{x^{\left( 0 \right)}\left( k \right)}{x^{\left( 1 \right)}\left( k-1 \right)}\text{，}k=2,3,\cdots ,n$
只有原始数据序列满足了这一条件才可以用于建模预测。

4 GM(1,1)预测模型的构建

设原始序列$X^{\left( 0 \right)}=\left( x^{\left( 0 \right)}\left( 1 \right) ,x^0\left( 2 \right) ,...,x^0\left( n \right) \right)$满足上述条件，其中$x^{\left( 0 \right)}\left( k \right) \ge 0$

（1）序列的累加处理
将原始序列进行一次累加生成后处理后，生成 $X^{\left( 0 \right)}$的1-AGO 序列 (累加生成序列)：
$X^{\left( 1 \right)}=\left( x^{\left( 1 \right)}\left( 1 \right) ,x^{\left( 1 \right)}\left( 2 \right) ,...,x^{\left( 1 \right)}\left( n \right) \right)$
其中
$x^{\left( 1 \right)}\left( k \right) =\sum\limits_{i=1}^k{x^{\left( 0 \right)}\left( i \right)}\ \ \ \ k=1,2,......,n$

（2）计算紧邻均值序列$Z^{\left( 1 \right)}=\left( z^{\left( 1 \right)}\left( 2 \right) ,z^{\left( 1 \right)}\left( 3 \right) ,...,z^{\left( 1 \right)}\left( n \right) \right)$
其中$z^{\left( 1 \right)}\left( k \right) =\frac{1}{2}\left( x^{\left( 1 \right)}\left( k \right) +x^{\left( 1 \right)}\left( k-1 \right) \right)$

（3）建立相关模型
建立一阶微分线性方程，即灰色微分方程，得到 GM(1,1)模型的均值形式：$x^{\left( 0 \right)}\left( k \right) +az^{\left( 1 \right)}\left( k \right) =b$

通过 GM(1,1) 模型相应的白化微分方程：
$\frac{\text{dx}^{\left( 1 \right)}}{dt}+ax^{\left( 1 \right)}=b$

其中 a 表示发展系数，反映了$\hat{x}^{\left( 1 \right)}$ 和$\hat{x}^{\left( 0 \right)}$

（4）计算 a、b

设 $\hat{a}$ 为待估参数向量，令 $\hat{a}=\left[ \begin{array}{c} \text{a}\\ b\\ \end{array} \right]$，利用最小二乘法求解，可得：
$\hat{a}=\left( B^TB \right) ^{-1}B^TY$
其中B、Y 分别为
$\text{B}=\left[ \begin{matrix} -z^{\left( 1 \right)}\left( 2 \right)& 1\\ -z^{\left( 1 \right)}\left( 3 \right)& 1\\ \vdots& \vdots\\ -z^{\left( 1 \right)}\left( n \right)& 1\\ \end{matrix} \right] \text{，Y}=\left[ \begin{array}{c} x^{\left( 0 \right)}\left( 2 \right)\\ x^{\left( 0 \right)}\left( 3 \right)\\ \vdots\\ x^{\left( 0 \right)}\left( \text{n} \right)\\ \end{array} \right]$

GM(1,1)模型中的参数 $-a$ 为发展系数，$b$

$-a$ 反应了 $X^{(1)}$ 及原始序列 $X^{(0)}$

作用量在物理学中是一个比较抽象的概念，它表示一个动力系统内在的演化趋势。

对原始数据序列进行计算求出的GM(1,1)模型中的 $b$

（5）建立 GM(1,1) 的时间响应式，即预测模型
$\left\{ \begin{array}{l} \hat{x}^{\left( 1 \right)}\left( k+1 \right) =\left( x^{\left( 0 \right)}\left( 1 \right) -\frac{b}{a} \right) e^{-ak}+\frac{b}{a},\text{k}=1,2\cdots ,\text{n}\\ \\ \hat{x}^{\left( 0 \right)}\left( k+1 \right) =\hat{x}^{\left( 1 \right)}\left( k+1 \right) -\hat{x}^{\left( 1 \right)}\left( k \right)\\ \end{array} \right.$

5 GM(1,1) 模型的检验

检验灰色预测模型效果主要包括三种方法：残差检验、关联度检验和后验差检验方法，以此来检验模型的准确性。

（1） 残差检验

按预测模型计算 $\hat{x}^{\left( 1 \right)}\left( k \right)$，并将 $\hat{x}^{\left( 1 \right)}\left( i \right)$累减生成 $\hat{x}^{\left( 0 \right)}\left( k \right)$，然后计算原始序列 $x^{\left( 0 \right)}\left( k \right)$ 与 $\hat{x}^{\left( 0 \right)}\left( k \right)$的绝对误差序列即相对误差序列。

$\left\{ \begin{array}{l} \text{残差：}\varepsilon \left( k \right) =x^{\left( 0 \right)}\left( k \right) -\hat{x}^{\left( 0 \right)}\left( k \right)\\ \\ \text{相对误差：}\Delta _k=\frac{\left| \varepsilon \left( k \right) \right|}{x^{\left( 0 \right)}\left( k \right)}\\ \\ \text{平均相对误差：}\Delta =\frac{1}{n-1}\sum\limits_{k=2}^n{\Delta _k}\\ \end{array} \right.$

若对相应序列残差有 $\Delta _k<10\%$ ，则代表用 GM(1,1) 模型的模拟精度较高，若 $\Delta _k<20\%$

GM(1,1) 模型建模精度：
$p=\left( 1-\varDelta \right) \times 100\%$

（2）关联度检验法

绝对误差：
$\varDelta ^{\left( 0 \right)}\left( k \right) =\left| x^{\left( 0 \right)}\left( k \right) -\hat{x}^{\left( 0 \right)}\left( k \right) \right|$

得出最大值、最小值

$\left\{ \begin{array}{l} \varDelta _{\max}^{\left( 0 \right)}=\max \varDelta ^{\left( 0 \right)}\left( k \right)\\ \\ \varDelta _{\min}^{\left( 0 \right)}=\min \varDelta ^{\left( 0 \right)}\left( k \right)\\ \end{array} \right.$
计算关联系数
$\eta \left( k \right) =\frac{\varDelta ^{\left( 0 \right)}_{\max}+\rho \varDelta ^{\left( 0 \right)}_{\min}}{\varDelta ^{\left( 0 \right)}\left( k \right) +\rho \varDelta ^{\left( 0 \right)}_{\max}}\text{，}\rho =0.5\text{，}k=1,2,\cdots ,n$
得出关联度
$r=\frac{1}{n}\sum_{k=1}^n{\eta \left( k \right)}$
根据经验，当 p =0.5 时，关联度 r >0.6 时可认为模型的预测效果较好。

（3）后验差检验
原始序列的均值
$\bar{x}^{\left( 0 \right)}=\frac{1}{n}\sum_{k=1}^n{x^{\left( 0 \right)}\left( k \right)}$

计算的原始序列的标准差：
$S_1=\sqrt{\frac{\sum{\left[ x^{\left( 0 \right)}\left( k \right) -\bar{x}^{\left( 0 \right)} \right]}^2}{n-1}}$

残差序列的均值：
$\varDelta ^{\left( 0 \right)}=\frac{1}{n}\sum_{k=1}^n{\left[ x^{\left( 0 \right)}\left( k \right) -\hat{x}^{\left( 0 \right)}\left( k \right) \right]}$
计算绝对误差序列的标准差：
$S_2=\sqrt{\frac{\sum{\left[ \Delta ^{\left( 0 \right)}\left( k \right) -\Delta ^{\left( 0 \right)} \right]}^2}{n-1}}$
计算后验差比值为：：
$C=\frac{S_2}{S_1}$
计算小误差概率：
$P=p\left\{ \left| \Delta ^{\left( 0 \right)}\left( k \right) -\Delta ^{\left( 0 \right)} \right|<0.6745S_1 \right\}$

GM(1,1) 通常用后验差检验方法来评价预测结果的好坏，主要根据其中的后验差比值（C）和小误差概率（P）这两个数值来检验模型。根据 C 和 P 的大小可以将灰色模型的预测精度分为以下几个等级，如下表

后验差比值 C	小误差概率 P	预测精度等级
<0.35	>0.95	好（一级）
<0.50	>0.80	合格（二级）
<0.65	>0.70	勉强合格（三级）
>=0.65	<=0.65	不合格 （四级）

6 GM(1,1) 模型的适用范围

建立于灰色理论基础上的GM(1,1)模型具有所需样本量小，无需考虑数据分布规律，模型预测精度较高，且易于计算检验，决定了GM(1,1)模型被广泛的应用于经济生活各个领域。但这并非说明GM(1,1)模型应用具有随意性与其他任何数学模型相同GM(1,1)模型也有其适用范围，若超出其适用范围使用GM(1,1)模型往往很难得到满意的的结果。

GM(1,1)模型的参数与GM(1,1)模型适用性之间的关系得到如下结论

（1）当 -a ≤ 0.3时，GM(1,1)模型可用于中长期预测；

（2）当 0.3 < -a ≤ 0.5时，GM(1,1)模型可用于短期预测，中长期预测需谨慎；

（3）当 0.5 < -a ≤ 0.8时，用GM(1,1)模型作短期预测慎用；

（4）当 0.8 < -a ≤ 1.0时，需采用残差修正GM(1,1)模型；

（5）当 1 < -a 时，不宜使用GM(1,1)模型。

也即

模型预测时长判断标准

7 GM(1,1) 残差模型

当原始数据序列$x^{\left( 0 \right)}\left( k \right)$建立的GM(1,1)模型检验不合格时，可以用GM(1,1)残差模型来修正。如果原始序列建立的GM(1,1)模型不够精确，也可以用GM(1,1)残差模型来提高精度。

若用原始序列 $x^{\left( 0 \right)}\left( k \right)$

$\hat{x}^{\left( 1 \right)}\left( k+1 \right) =\left( x^{\left( 0 \right)}\left( 1 \right) -\frac{b}{a} \right) e^{-ak}+\frac{b}{a},\text{k}=1,2\cdots ,\text{n}$

可获得生成序列$x^{\left( 1 \right)}\left( k \right)$的预测值，定义残差序列$e^{\left( 0 \right)}\left( k \right) =x^{\left( 1 \right)}\left( k \right) -\hat{x}^{\left( 1 \right)}\left( k \right)$。则对应的残差序列为：$e^{\left( 0 \right)}\left( k \right) =\left\{ e^{\left( 0 \right)}\left( 1 \right) ,e^{\left( 0 \right)}\left( 2 \right) ,\cdots ,e^{\left( 0 \right)}\left( n \right) \right\}$

计算其生成序列$e^{\left( 1 \right)}\left( k \right)$，并据此建立相应的GM(1,1)模型：

$\hat{e}^{\left( 1 \right)}\left( k+1 \right) =\left( e^{\left( 0 \right)}\left( 1 \right) -\frac{b_e}{a_e} \right) e^{-ak}+\frac{b_e}{a_e},\text{k}=1,2\cdots ,\text{n}$
得修正模型
$\hat{x}^{\left( 1 \right)}\left( k+1 \right) =\left( x^{\left( 0 \right)}\left( 1 \right) -\frac{b}{a} \right) e^{-ak}+\frac{b}{a}+\delta \left( k-i \right) \left( -a_e \right) \left[ e^{\left( 0 \right)}\left( 1 \right) -\frac{b_e}{a_e} \right] e^{-ak}$
其中
$\delta \left( k-i \right) =\left\{ \begin{array}{l} 1\text{，}k\ge i\\ 0\text{，}k<i\\ \end{array} \right.$

为修正系数。

应用此模型时要考虑：

（1）一般不是使用全部残差数据来建立模型，而只是利用了部分残差。

（2）修正模型所代表的是差分微分方程，其修正作用与$\delta \left( k-i \right)$

R语言实现

灰色系统模型GM(1,1)的R语言实现

灰色预测GM11，R语言源程序

python

'''灰色预测函数'''
def GM11(x0): #自定义灰色预测函数
    import numpy as np
    x1 = x0.cumsum() # 生成累加序列
    z1 = (x1[:len(x1)-1] + x1[1:])/2.0 # 生成紧邻均值（MEAN）序列，比直接使用累加序列好，共 n-1 个值
    z1 = z1.reshape((len(z1),1))
    B = np.append(-z1, np.ones_like(z1), axis = 1)    # 生成 B 矩阵
    Y = x0[1:].reshape((len(x0)-1, 1))    # Y 矩阵
    [[a],[u]] = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y)    #计算参数
    f = lambda k: (x0[0]-u/a)*np.exp(-a*(k-1))-(x0[0]-u/a)*np.exp(-a*(k-2))    #还原值
    delta = np.abs(x0 - np.array([f(i) for i in range(1,len(x0)+1)]))    # 计算残差
    C = delta.std()/x0.std()
    P = 1.0*(np.abs(delta - delta.mean()) < 0.6745*x0.std()).sum()/len(x0)
    return f, a, u, x0[0], C, P #返回灰色预测函数、a、b、首项、方差比、小残差概率