时间序列分析----灰度模型

最新推荐文章于 2023-11-16 15:56:55 发布

暗夜猎手-大魔王

最新推荐文章于 2023-11-16 15:56:55 发布

阅读量5.2k

点赞数 7

分类专栏：算法与数据结构

本文链接：https://blog.csdn.net/u014106644/article/details/90405842

版权

算法与数据结构专栏收录该内容

108 篇文章 11 订阅

订阅专栏

灰度预测法是一种对含有不确定因素的系统的预测方法，灰色系统是位于白色系统和黑色系统之间的一种系统。

白色系统指的是一个系统内部的特征是完全已知的，使用者不仅知道系统的输入-输出关系，还知道实现输入-输出的具体方式，譬如函数表达式，微分方程的变化公式，或者物理学的基本定律。比方说牛顿第二定律F=ma, 使用者只需要知道物体的质量和加速度，就可以通过牛顿第二定律求出所使用的力F的具体值。或者说当物体的质量固定之后，已知不同的加速度，就可以求出不同的力的值，力和加速度之间有明确的关系表达式。这种系统就称为白色系统。

黑色系统和白色系统就有鲜明的对比，使用者完全不清楚黑色系统的内部特征，只是知道一些输入和相应的输出。使用者需要通过这一系列的输入和输出的对应关系来判断这个黑色系统的内部特征，对于没有办法打开的黑箱，使用者需要通过输入和输出来建立模型来了解这个黑箱。比方说：高校在进行招生的时候，其实对考生是几乎不了解的，这个时候就需要通过高考或者自主招生考试来判断学生的知识掌握程度和运用知识解决问题的能力。通过测试的结果来判断学生的综合能力的大小，进而判断是不是符合高校的专业需求。

而灰度系统恰好位于白色系统和黑色系统之间，灰度系统内部一部分是已知的，但是另一部分却是未知的。比方说：经济学中的模型，通过数学或者统计学的方法，可以判断某些因素之间确实有关系，但是却不能够完全的判断整个系统的情况，对于很多未知的情况，这些模型就有一定的局限性。于是，灰度系统就需要通过判断系统各个因素之间的发展规律，进行相应的关联分析。通过原始的数列来生成规律性更强的数列，通过生成的数列来建立相应的微分方程模型，从而预测数列的未来发展趋势。灰度模型使用等时间距观测到的数据值来构造灰色预测模型，从而达到能够预测未来某一时刻的数值的目的。

灰度模型首先将原始数据进行某种方式的生成处理，从杂乱无章的现象中发现内在规律，常见的生成处理方式有累加，累减以及均值处理等。使用累加生成方式，将原始数据各时刻数据依个累加得到新的生成数据，由于原始数据非负，已知生成数据必为非减递增数列，通过累加可以看出灰量积累过程发展态势，使原始数据中蕴含的规律加以显化，从而更好对于数据中规律进行归纳和利用。

灰色生成数列
灰色系统理论认为，尽管客观表象复杂，但总是有整体功能的，因此必然蕴含某种内在规律。关键在于如何选择适当的方式去挖掘和利用它。灰色系统时通过对原始数据的整理来寻求其变化规律的，这是一种就数据寻求数据的现实规律的途径，也就是灰色序列的生产。一切灰色序列都能通过某种生成弱化其随机性，显现其规律性。数据生成的常用方式有累加生成、累减生成和加权累加生成。
累加生成（AGO）

原始序列： $x^{(0)}=(x^0(1),x^0(2),...,x^0(n))$

累加生成序列：

$x^1(k)=\displaystyle\sum_{i=1}^{k}x^0(i),k=1,2,..n$

$x^{(1)}=(x^1(1),x^1(2),...,x^1(n))$

得到的新数列为数列 $x^{(0)}$ 的1次累加生成数列。类似的有

$x^{(r)}(k)=\displaystyle\sum_{i=1}^{k}x^{(r-1)}(i),k=1,2,...n$

累加生成能使任意非负数列、摆动的与非摆动的，转化为非减的、递增的。
累减生成(IAGO)

原始数列： $x^{(1)}=(x^1(1),x^1(2),...,x^1(n))$

累减生成：

$x^{(0)}(k)=x^{(1)}(k)-x^{(1)}(k-1),k=2,3...n$

加权邻值生成

原始序列为： $x^{(0)}=(x^0(1),x^0(2),...,x^0(n))$

加权临值生成：

$z^{(0)}(k)=\alpha x^{(0)}(k)+(1-\alpha )x^{(0)}(k-1), k=2,3...n$

灰度模型预测计算流程图
在这里插入图片描述
对于原始序列进行级比检验

原始时间序列： $x^{(0)}=(x^0(1),x^0(2),...,x^0(n))$

计算级比： $\lambda(k)=\frac{x^{(0)}(k-1)}{x^{(0)}(k)},k=2,3..n$

如果所有级别都满足区间 $(e^{\frac{-2}{n+1}},e^{\frac{2}{n+1}})$ ，则可以使用灰度模型预测，否则进行平移变换

$y^{(0)}(k)=x^{(0)}(k)+c$ 使其满足上述级比要求

灰度模型GM(1,1)

对原始序列进行一次累加得累加序列

$x^{(1)}=(x^1(1),x^1(2),...,x^1(n))$ $x^1(k)=\displaystyle\sum_{i=1}^{k}x^0(i),k=1,2,..n$

灰度系统是对离散序列建立的微分方程，其中GM(1,1)建立一阶微分方程模型，其微分方程如下：

$\frac{dx^{(1)}t}{dt}+ax^{(1)}(t)=b$ 其中a为发展系数 b为控制系数

如果 $\delta$ t很小,取 $\delta$ t=1则根据微分的意义有：

$\frac{dx^{(1)}t}{dt}=\frac{x^{(1)}(t+1)-x^{(1)}(t)}{\delta t}=x^{(1)}(t+1)-x^{(1)}(t)=x^{(0)}(t)$

如果 $\delta t=1$ 很小，则在很短的时间内，累加序列 $x^{(1)}(t)$ 以及 $x^{(1)}(t+\delta t)$ 不会出现突变量，可用 $x^{(1)}(t)$ 以及 $x^{(1)}(t+\delta t)$ 平均值作为 $x^{(1)}(t)$ 背景值，取α=0.5对累加时间序列取均值：

$z^{(1)}(t)=x^{(1)}(t)=\frac{1}{2}(x^{(1)}(t)+x^{(1)}(t+1))$

将上式代入微分方程得离散方程为：

$x^{(0)}(t)+az^{(1)}(t)=b, t=1,2,...n$

由最小二乘法进行数据拟合

$a,b)^T=(B^TB)^{-1}B^TY$

其中 $B=\begin{bmatrix} -z^{(1)}(2)&1\\ -z^{(1)}(3)&1\\ ...&1\\ -z^{(1)}(n)&1\end{bmatrix}$ $Y=\begin{bmatrix} x^{(0)}(2),&x^{(0)}(3),...x^{(0)}(n)\end{bmatrix}^T$

得到微分方程解为：

$x^1_p(k)=(x^0(1)-\frac{b}{a})e^{-a(k-1)}+\frac{b}{a}, k=1,2..n$

得到最终解：

$x^0_p(k)=x^1_p(k)-x^1_p(k-1)=(x^0(1)-\frac{b}{a})e^{-a(k-1)}(1-e^{a}),k=2,3..n$

由上式将初始值代入，依次计算预测值，得到预测值序列

$x^0_p(1),x^0_p(2),...,x^0_p(n)$

对预测值序列进行相对残差检验

$\varepsilon(k)=\frac{x^0_p(k-1)}{x^0_p(k)},k=2,3,...n$

若εk <0.2时，模型预测效果一般，若εk <0.1，则说明模型预测效果较好。

对级比偏差值检验

$\rho(k)=(1-\frac{1-0.5a}{1+0.5a}) \lambda(k)$

若 $\rho(k)$ <0.2，模型预测效果一般，若 $\rho(k)$ <0.1，则说明模型预测效果较好。

灰度Verhulst模型

根据原始序列求解累加序列，中值序列，建立非线性常微分方程

微分方程 $\frac{dx^{(1)}t}{dt}+ax^{(1)}(t)=b(x^{(1)}(t))^2$

离散方程 $x^{(0)}(k)+az^{(1)}(k)=b(z^{(1)}(k))^2, k=2,3...n$

得到 $a,b)^T=(B^TB)^{-1}B^TY$

其中 $B=\begin{bmatrix} -z^{(1)}(2)&z^{(1)}(2)^2\\ -z^{(1)}(3)&z^{(1)}(3)^2\\ ...&...\\ -z^{(1)}(n)&z^{(1)}(n)^2\end{bmatrix}$ $Y=\begin{bmatrix} x^{(0)}(2),&x^{(0)}(3),...x^{(0)}(n)\end{bmatrix}^T$

最终解的：

$x^{(1)}_p(k)=\frac{ax^{(0)}(1)}{bx^{(0)}(1)+(a-bx^{(0)}(1))e^{a(k-1)}},k=1,2..n$

$x^{(0)}_p(k)=x^{(1)}_p(k)-x^{(1)}_p(k-1),k=2,3,...n$

灰度模型DGM(2,1)

使用原数列，建立微分方程

$\frac{d^2x^{(0)}t}{dt^2}+a\frac{dx^{(0)}(t)}{dt}=b$

得到 $a,b)^T=(B^TB)^{-1}B^TY$

其中 $B=\begin{bmatrix} -z^{(0)}(2)&1\\ -z^{(1)}(3)&1\\ ...&...\\ -z^{(0)}(n)&1\end{bmatrix}$ $Y=\begin{bmatrix} x^{(0)}(2)-x^{(0)}(1),...x^{(0)}(n)-x^{(0)}(n-1)\end{bmatrix}^T$

解的
$x^{(0)}_p(k)=(\frac{b}{a^2}-\frac{x^{(0)}(1)}{a})(1-e^a)e^{-a(k-1)}+\frac{b}{a}, k=2,3..n$

灰度模型GM(1,1)编程实现

使用python来使用GM(1,1)模型来预测数值，一下是一段时间内cpu利用率的数据，要求来进行拟合，预测某一段时间内的cpu利用率
在这里插入图片描述
利用灰度模型建模，依次使用前N个数据点来预测第N+1个数据点，分别测试在不同N下数据拟合的优劣
代码如下：

# -*- coding: utf-8 -*-
"""
Created on Fri Mar  1 11:03:25 2019

@author: chen
"""
import pandas as pd
import numpy as np

#matplotlib inline
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['font.sans-serif'] = ['FangSong'] # 指定默认字体
matplotlib.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题

class GrayForecast():
#初始化
    def __init__(self, data, datacolumn=None):
        pass
#级比校验
    def level_check(self):
        pass
#GM(1,1)建模
    def GM_11_build_model(self, forecast=5):
        pass
#预测
    def forecast(self, time=5, forecast_data_len=5):
        pass
#打印日志
    def log(self):
        pass
#重置
    def reset(self):
        pass
#作图
    def plot(self):
        pass
    
    def __init__(self, data, datacolumn=None):
    
        print('init start')
        if isinstance(data, pd.core.frame.DataFrame):
            self.data=data
            try:
                self.data.columns = ['数据']
            except:
                if not datacolumn:
                    raise Exception('您传入的dataframe不止一列')
                else:
                    self.data = pd.DataFrame(data[datacolumn])
                    self.data.columns=['数据']
                    #print('11111')
        elif isinstance(data, pd.core.series.Series):
            self.data = pd.DataFrame(data, columns=['数据'])
            print('222222')
        else:
            self.data = pd.DataFrame(data, columns=['数据'])
            print('33333')
    
        self.forecast_list = self.data.copy()
        #print(self.forecast_list['数据'])
        if datacolumn:
            self.datacolumn = datacolumn
        else:
            self.datacolumn = None
        #save arg:
        #        data                DataFrame    数据
        #        forecast_list       DataFrame    预测序列
        #        datacolumn          string       数据的含义    
        
    def level_check(self):
        # 数据级比校验
        n = len(self.data)
        lambda_k = np.zeros(n-1)
        for i in range(n-1):
            lambda_k[i] = self.data.ix[i]["数据"]/self.data.ix[i+1]["数据"]
            if lambda_k[i] < np.exp(-2/(n+1)) or lambda_k[i] > np.exp(2/(n+2)):
                flag = False
        else:
            flag = True
    
        self.lambda_k = lambda_k
    
        if not flag:
            print("级比校验失败，请对X(0)做平移变换")
            return False
        else:
            print("级比校验成功，请继续")
            return True
    
    #save arg:
    #        lambda_k            1-d list
            
    def GM_11_build_model(self, forecast=5):
        #参数time为向后预测的次数，forecast_data_len为每次预测所用末尾数据的条目数
        #print(self.data)
        if forecast > len(self.data):
            raise Exception('您的数据行不够')
        X_0 = np.array(self.forecast_list['数据'].tail(forecast))
        #print('X0')
        #print(X_0)
    #       1-AGO
        X_1 = np.zeros(X_0.shape)
        for i in range(X_0.shape[0]):
            X_1[i] = np.sum(X_0[0:i+1])
    #       紧邻均值生成序列
        Z_1 = np.zeros(X_1.shape[0]-1)
        for i in range(1, X_1.shape[0]):
            Z_1[i-1] = -0.5*(X_1[i]+X_1[i-1])
    
        B = np.append(np.array(np.mat(Z_1).T), np.ones(Z_1.shape).reshape((Z_1.shape[0], 1)), axis=1)
        Yn = X_0[1:].reshape((X_0[1:].shape[0], 1))
    
        B = np.mat(B)
        Yn = np.mat(Yn)
        a_ = (B.T*B)**-1 * B.T * Yn
    
        a, b = np.array(a_.T)[0]
    
        X_ = np.zeros(X_0.shape[0])
        def f(k):
            return (X_0[0]-b/a)*(1-np.exp(a))*np.exp(-a*(k))
        
        self.forecast_list.loc[len(self.forecast_list)] = f(X_.shape[0])
        print(f(X_.shape[0]))
        
    def forecast(self, time=5, forecast_data_len=5):
        for i in range(time):
            self.GM_11_build_model(forecast=forecast_data_len)
        print(self.forecast_list['数据'])
        
    def forecast2(self, num):
        print(self.data[0:num])
        for index in range(len(self.data)-(num-1)):
            self.GM_11_build_model2(index, num)
        #print(self.forecast_list['数据'][0:len(self.forecast_list)])
        
    def GM_11_build_model2(self, index, num):
        #参数time为向后预测的次数，forecast_data_len为每次预测所用末尾数据的条目数
        #print(self.data)
        X_0 = np.array(self.data['数据'][index:index+num])
        #print(index)
        #print('X0')
        #print(X_0)
    #       1-AGO
        X_1 = np.zeros(X_0.shape)
        for i in range(X_0.shape[0]):
            X_1[i] = np.sum(X_0[0:i+1])
    #       紧邻均值生成序列
        Z_1 = np.zeros(X_1.shape[0]-1)
        for i in range(1, X_1.shape[0]):
            Z_1[i-1] = -0.5*(X_1[i]+X_1[i-1])
    
        B = np.append(np.array(np.mat(Z_1).T), np.ones(Z_1.shape).reshape((Z_1.shape[0], 1)), axis=1)
        Yn = X_0[1:].reshape((X_0[1:].shape[0], 1))
    
        B = np.mat(B)
        Yn = np.mat(Yn)
        a_ = (B.T*B)**-1 * B.T * Yn
    
        a, b = np.array(a_.T)[0]
    
        X_ = np.zeros(X_0.shape[0])
        def f(k):
            return (X_0[0]-b/a)*(1-np.exp(a))*np.exp(-a*(k))
        r = f(X_.shape[0])
        self.forecast_list.loc[index+5] = r
        print(r)
            
    def log(self):
        res = self.forecast_list.copy()
        if self.datacolumn:
            res.columns = [self.datacolumn]
        return res
    
    def reset(self):
        self.forecast_list = self.data.copy()
        
    def plot(self):
        self.forecast_list.plot()
        self.data.plot()
        if self.datacolumn:
            plt.ylabel(self.datacolumn)
            plt.legend([self.datacolumn])
            
def calculateE():
    f = open("预测数据.csv")
    df = pd.read_csv(f)
    df.tail()
    print(df)
    print(df['原始数据'])
    print(len(df['3个数据点']))
    print(score(df['3个数据点']-df['原始数据']))
    print(score(df['4个数据点']-df['原始数据']))
    print(score(df['5个数据点']-df['原始数据']))
    print(score(df['6个数据点']-df['原始数据']))
    print(score(df['7个数据点']-df['原始数据']))
    print(score(df['8个数据点']-df['原始数据']))
    
def score(df):
    X = np.array(df)
    XT = np.mat(X).T
    a = np.mat(X)*XT / len(df)
    return a
        
def test():
    f = open("cpu利用率数据.csv")
    df = pd.read_csv(f)
    df.tail()
    print(df)
    gf = GrayForecast(df, '原始数据')
    print('test======')
    gf.forecast2(8)
    gf.log()
    gf.plot()

最终结果
在这里插入图片描述
可知四个数据点时，预测效果较好。

灰度模型GM(1,1) 通过对于原始序列累加，使任意非负数列、摆动的与非摆动的，转化为非减的、递增的，然后利用一阶线性微分方程来建立模型，其实质是猜测累加序列在某一段时间内变化趋势为类指数变化趋势，然后来进行模拟求解。

参考链接：

时间序列模型之灰度模型 https://zr9558.com/2015/09/14/greysystemtheory/

时间序列模型之灰度模型 https://blog.csdn.net/cainiaozr/article/details/48439371

灰色预测模型GM(1,1) 与例题分析 https://blog.csdn.net/qq547276542/article/details/77865341

【数学建模】灰色预测及Python实现 https://www.jianshu.com/p/a35ba96d852b

暗夜猎手-大魔王

关注

7
点赞
踩
51

收藏

觉得还不错? 一键收藏
0
评论
时间序列分析----灰度模型

灰度预测法是一种对含有不确定因素的系统的预测方法，灰色系统是位于白色系统和黑色系统之间的一种系统。白色系统指的是一个系统内部的特征是完全已知的，使用者不仅知道系统的输入-输出关系，还知道实现输入-输出的具体方式，譬如函数表达式，微分方程的变化公式，或者物理学的基本定律。比方说牛顿第二定律F=ma, 使用者只需要知道物体的质量和加速度，就可以通过牛顿第二定律求出所使用的力F的具体值。或者说当物体的...
复制链接

扫一扫

专栏目录