灰色预测原理及实例(附代码)

2023-10-31

灰色预测

引言

古人说:“凡事预则立,不预则废。”办任何事情之前,必须先调查研究,摸清情况,深思熟虑,有科学的预见,周密的计划,这样才能达到预期的成功。

所谓预测,就是人们根据可获得的历史和现实数据,资料,运用一定的科学方法和手段,对人类社会、政治、军事、科学技术等发展趋势作出科学推测,以指导未来行动的方向,减少处理未来事件的盲目性。
灰色预测基于人们对系统演化不确定性特征的认识,运用序列算子对原始数据进行生成、处理,挖掘系统演化规律,建立灰色系统模型,对系统的未来状态作出科学的定量预测。同时,对于一个具体问题,究竟应该选择什么样的预测模型,应以充分的定性分析结论作为依据。模型的选择不是一成不变的,一个模型要经过多种检验才能判定其是否合理,是否有效。只有通过检验的模型才能用作预测模型。

灰色预测简介

通过对原始数据的生成处理来寻求系统变动的规律。生成数据序列较强的规律性,可以用它来建立相应的微分方程模型,从而预测事物未来的发展趋势和未来状态。

灰色预测的类型
  • 时间序列预测
  • 灾变预测
  • 波形预测
  • 系统预测

本篇文章例子为时间序列预测,其他类型使用方面较少并且比较深入,大家可以去借鉴别人文章。

最简单的模型:GM(1,1)
GM(1,1)模型

G:Grey(灰色);
M:模型;
(1,1):只含有一个变量的一阶微分方程模型;

实例
年份 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004
污水/亿吨 174 179 183 189 207 234 220.5 256 270 285

题目要求:根据给定表中数据预测未来10年的长江污水量

原理及求解

数据处理方法:

1.累加生成

原始序列为X(x=1 2 3…)
生成序列为Y(y=1 2 3…)
则:
Y(1)=X(1)
Y(2)=Y(1)+X(2)
Y(3)=Y(2)+X(3)

所谓的累加生成,就是将同一序列中的数据逐次相加以生成新的数据的一种手段,累加前的数列成为原始数列。累加后的数列称为原始数列。累加后的数列成为生成数列。累加生成是使灰色系统变白的一种方法,它在灰色系统理论中占有极其重要的地位。通过累加生成可以看出灰量累积过程的发展态势,使杂乱无章的原始数据中蕴含的积分特性或规律加以显化。

2.累减生成

原始序列为X(x=1 2 3…)
生成序列为Y(y=1 2 3…)
则:
Y(1)=X(1)
Y(2)=X(2)-X(1)
Y(3)=X(3)-X(2)

累减生成是在获取增量信息时常用的生成,多数情况下累减生成对累加生成起还原作用,即累减生成是累加生成的逆运算。

3.均值生成

原始序列为X(x=1 2 3…)
生成序列为Y(y=1 2 3…)
Y(1)=X(1)
Y(2)=aX(2)+(1-a)X(1)
Y(3)=aX(3)+(1-a)X(2)

在收集数据的时候,由于一些不易克服的困难导致数据序列出现空缺或无法使用的异常数据,需要在数据预处理中解决。均值生成是常用的构造新数据、填补老数据空穴、生成新数列的方法。

求解步骤框图

流程图:

Created with Raphaël 2.3.0 开始 输入数据 级比检验 累加生成 邻均值生成 构造数据矩阵及数据向量 计算发展系数及灰作用量 求解时间响应函数并进行预测 精确度检验 结束 思考 yes no yes no
求解步骤
Step1:

首先建立时间序列如下:

  • X=(x(1),x(2),…x(n))
  • 求级比L(k):
  • L(k)=x(k-1)/x(k)
  • 当所有的L(k)落在【exp(-2/(n+1)),exp(2/(n+1))】区间之内,认为可以做比较满意的GM(1,1)建模。(并不是说由有的数据没落入区间之内就不能建模,只是落在区间之内建模效果比较好)。
Step2:

将原始数据时间序列进行累加生成

  • 设累加生成之后的序列为X’
Step3:

对累加生成之后的序列进行邻均值生成
邻均值生成是对等时距数列,用相邻数据的平均值构造生成新的数据。
设新生成的邻均值序列为Z,则:

  • Z(k)=(X’(k)+X’(k+1))/2 ,k=1,2,…n
Step4:

构造数据矩阵B及数据向量Y
在这里插入图片描述

Step5:

计算发展系数a及灰作用量b

  • 计算方法如下:在这里插入图片描述
Step6:

建立模型求解时间响应函数并进行预测

  • GM(1,1)模型X(k)+aZ(k)=b的白化方程:
    在这里插入图片描述
    的时间响应函数为:
    在这里插入图片描述

Step7:精确度检验

  • 相对残差检验
  • 方差比检验
  • 小误差概率检验
  • 相对残差检验法:
  • 设实际数据序列为X(k)
  • 设模拟数据序列为X’(k)
  • 则残差&为:
    &(k)=X(k)-X‘(k)
  • 相对误差*k=|&(k)|/X(k)
  • 方差比检验法:
  • 设残差序列为&(k),其序列标准差为A。
  • 原始序列为X(k),其序列标准差为B。
  • 则方差比C=A/B。
  • 小误差概率检验法:
  • 设小误差概率为p,残差序列平均值为*&:
  • P={|&(k)-*&|<0.6745B}
  • 则p=P/原始数据序列长度
附小误差概率p及方差比检验标准(可在题目无要求精度时相对比较):
预测精度等级 p(小误差概率) C(方差比)
一级 >0.95 <0.35
二级 >0.8 <0.5
三级 >0.7 <0.65
四级 >0.6 >=0.65(不合格)
一般用作模型预测的话,精度等级要达到二级以上。
The End : MATLAB求解代码!!!
clc;clear;
%建立符号变量a(发展系数)b(灰作用量)
syms a b;
c = [a b]';

%原始数列
A = [174, 179, 183, 189, 207, 234, 220.5, 256, 270, 285];
%级比检验
n = length(A);
min=exp(-2/(n+1));
max=exp(2/(n+1));
for i=2:n
    ans(i)=A(i-1)/A(i);
end
ans(1)=[];
for i=1:(n-1)
    if ans(i)<max&ans(i)>min
    else
        fprintf('第%d个级比不在标准区间内',i)
        disp(' ');
    end
end

%对原始数列 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:n);
%计算残差序列
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;
好了,这就是本篇博客的全部内容了,首先非常感谢您能看到本篇文章最后,那么如果看到它的你还有什么地方不懂,可以单独来问我,I'm goldsun,期待下次相遇。
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

灰色预测原理及实例(附代码) 的相关文章

  • MATLAB:生成给定三种颜色的颜色图

    我正在尝试在 MATLAB 中生成给定三种颜色 最高值 零值和最低值 的颜色图 我的思维过程是从最高端到中间循环 并将每个步骤存储到一个 3xN 第一列是 R 第二列是 G 第三列是 B 矩阵 所以我正在使用 fade from high
  • Matlab 中 interp2 的类似 OpenCV Api

    有没有类似的功能 其工作原理与 interp2 x y frame z xd yd linear 0 在 OpenCV 中 功能cv remap 几乎可以满足您的要求 请参阅文档here http docs opencv org modul
  • 与超类和子类构造函数接口

    我在 matlab 文档和之前有关使用 matlab 继承和类构造函数创建接口的问题中找不到帮助 为了使其整洁 放在一个包内 我可以将其压缩如下 而不是拖拽代码 一套 MyPkg有一个超类Super和一些子类Sub1 Sub2 我的大多数属
  • 使用 libsvm 交叉验证后重新训练

    我知道交叉验证用于选择好的参数 找到它们后 我需要在不使用 v 选项的情况下重新训练整个数据 但我面临的问题是 在使用 v 选项训练后 我得到了交叉验证精度 例如 85 没有模型 我看不到 C 和 gamma 的值 在这种情况下我该如何重新
  • 将 Android 应用程序与服务器上的 Matlab 应用程序连接

    我正在 Android 上开发一个应用程序 它将获取图像输入 并将该输入传递到安装 MATLAB 应用程序的服务器 MATLAB 应用程序将计算结果并将其返回到该 Android 应用程序 我想知道我可以使用哪个服务器 如何将 MATLAB
  • MATLAB 中的抗锯齿线和标记

    您好 我在 MATLAB 中有一张图像 我希望这条线是平滑的 看看从 0 4 到 0 8 的线 这太可怕了 当在图中使用 LineSmoothing on 运算符时 我得到了这个 我在线条上做得很好 但它也使标记变得平滑 而且它们太可怕了
  • 如何建立数据流挖掘的滑动窗口模型?

    我们遇到的情况是 流 来自传感器的数据或服务器上的点击流数据 采用滑动窗口算法 我们必须将最后 例如 500 个数据样本存储在内存中 然后 这些样本用于创建直方图 聚合并捕获有关输入数据流中异常的信息 请告诉我如何制作这样的滑动窗 如果您询
  • 在 MATLAB 中将数据拟合到 B 样条

    我正在尝试估计矩阵形式的时间序列数据中的缺失值 列代表时间点 即现在 我想将矩阵的每一行拟合到 B 样条曲线 并用它来估计缺失值 我可以使用 MATLAB 将数据拟合到普通样条曲线 但我完全陷入尝试找出如何拟合数据以创建 B 样条曲线的困境
  • 将单元格转换为双精度

    gt gt C 1 2 CF 2 C 1 2 CF 2 gt gt whos C Name Size Bytes Class Attributes C 2x2 478 cell 我怎样才能转换C into double以便 gt gt C
  • 如何读取 10 位原始图像?其中包含 RGB-IR 数据

    我想知道如何从我的 10 位原始 它有 rgb ir 图像数据 数据中提取 RGB 图像 如何使用 Python 或 MATLAB 进行阅读 拍摄时的相机分辨率为 1280x720 室内照片图片下载 https drive google c
  • 在 MATLAB 中用两个值替换向量值

    我必须创建一个以向量作为输入的函数v和三个标量a b and c 该函数替换了的每个元素v等于a有一个二元素数组 b c 例如 给定v 1 2 3 4 and a 2 b 5 c 5 输出将是 out 1 5 5 3 4 我的第一次尝试是尝
  • 是否有一个函数可以检查矩阵是否对角占优(行占优)

    矩阵是对角占优 http en wikipedia org wiki Diagonally dominant matrix 按行 如果对角线处的值在绝对意义上大于该行中所有其他绝对值的总和 对于列也是如此 只是相反 matlab中有没有函数
  • 傅里叶变换定理 matlab

    我目前正在尝试理解二维傅里叶位移定理 根据我到目前为止所了解到的情况 图像空间中的平移会导致相位差异 但不会导致频率空间中的幅度差异 我试图用一个小例子来演示这一点 但它只适用于行的移位 而不适用于列的移位 这是一个小演示 我只在这里显示幅
  • 在 MATLAB 中创建共享库

    一位研究人员在 MATLAB 中创建了一个小型仿真 我们希望其他人也能使用它 我的计划是进行模拟 清理一些东西并将其变成一组函数 然后我打算将其编译成C库并使用SWIG https en wikipedia org wiki SWIG创建一
  • 频域和空间域的汉明滤波器

    我想通过在 MATLAB 中应用汉明滤波器来消除一维信号中的吉布斯伪影 我所拥有的是k1这是频域中的信号 我可以通过应用 DFT 来获取时域信号k1 s1 ifft ifftshift k1 该信号具有吉布斯伪影 现在 我想通过 A 乘以汉
  • 为什么 MATLAB 在打印大量 (.png) 图形时速度会变慢?

    我正在将大量数字打印为 png 文件 每个图都是数据矩阵中的一列图 我获取 png 文件并将它们串在一起形成动画 我的问题是 前几百张图像打印得很快 但创建每个新图形的时间却迅速增加 从前几百个 png 文件的约 0 2 秒到第 800 个
  • 使用 R2010b 中的符号工具箱来求解和/或 linsolve

    我前几天问了一个问题here https stackoverflow com questions 20317038 matlab linear congruence solver that supports a non prime modu
  • 如何获取MATLAB句柄对象的ID?

    当我尝试使用时出现问题MATLAB 句柄对象 http www mathworks com help techdoc ref handle html作为关键值MATLAB 容器 Map http www mathworks com help
  • 归一化互相关的基础知识

    我正在尝试使用范数校正2 归一化互相关 http en wikipedia org wiki Cross correlation Normalized cross correlation 来自 MATLAB 用于计算发育中胚胎中移动形状的速
  • 在Matlab中对字符进行分组并形成矩阵

    我有 26 个字符 A 到 Z 我将 4 个字符组合在一起 并用空格分隔以下 4 个字符 如下所示 abcd efgh ijkl mnop qrst uvwx yz 我的Matlab编码如下 str abcdefghijklmnopqrst

随机推荐

  • uni-app开发微信小程序,button通过数组的length判断disabled无效(数组length === 0写法无效)

    错误写法
  • caffe特征提取/C++数据格式转换

    Caffe生成的数据分为2种格式 Lmdb 和 Leveldb 它们都是键 值对 Key Value Pair 嵌入式数据库管理系统编程库 虽然lmdb的内存消耗是leveldb的1 1倍 但是lmdb的速度比leveldb快10 至15
  • 国产操作系统进入被彻底抛弃的时代

    当倪光南正在不断呼喊支持国产操作系统的时候 国产操作系统却迎来了噩梦 国产操作系统接连倒闭 国产操作系统进入一个被国家彻底抛弃的时代 红旗linux梦断国产操作系统 今年2月中科红旗linux因为缺钱倒闭解散了 一直以来做得最好的国产操作系
  • 图形学数学基础之基本蒙特卡罗尔积分(Monte Carlo Integration)

    作者 i dovelemon 日期 2017 07 29 来源 CSDN 主题 Monte Carlo Integration 引言 好久没有写博客了 最近一直在忙于工作 同时GLB库中关于PBR的渲染算法 一直卡住 无法实现下去 不过在这
  • dd大牛的《背包九讲》

    P01 01背包问题 题目 有N件物品和一个容量为V的背包 第i件物品的费用是c i 价值是w i 求解将哪些物品装入背包可使这些物品的费用总和不超过背包容量 且价值总和最大 基本思路 这是最基础的背包问题 特点是 每种物品仅有一件 可以选
  • 【HCNP路由交换学习指南】学习笔记丨第07章 BGP

    07 BGP BGP 的基本概念 BGP 对等体关系类型 IBGP 水平分割原则 路由黑洞问题及 BGP 同步规则 路由通告 Router ID 报文类型及格式 Open 报文 Update 报文 Keepalive 报文 Notifica
  • PaddleOCR使用笔记之模型训练

    目录 简介 模型训练 步骤一 文本检测模型 detection 1 准备训练数据集 2 下载预训练模型 模型介绍 下载预训练模型 3 开始训练 断点训练 4 模型评估 5 模型测试 6 训练模型转inference模型 步骤二 文本识别模型
  • RabbitMQ保证消息的一致性解决方案

    RabbitMQ保证消息的一致性 一 采用confirm消息确认机制及return返回机制 确保消息发送成功 二 将队列以及消息设置持久化 保证rabbitmq突然宕机消息仍然存在 三 手动确认接收消息方式 消息处理失败拒收重回队列 1 y
  • 后端响应是否成功、信息、操作码响应前端及异常处理

    异常处理流程 1 自定义异常类型 2 自定义错误代码及错误信息 3 对于可预知的异常由程序员在代码中主动抛出 由SpringMVC统一捕获 可预知异常是程序员在代码中手动抛出本系统定义的特定异常类型 由于是程序员抛出的异常 通常异常信息比较
  • java反射(从认识到应用)-黑马笔记

    此文章是观看黑马雷哥关于Java反射的视频做的笔记 如有错请多指教 一 认识反射 反射在JavaAPI中的详解 说白了 反射就是 加载类 并允许以编程方式解剖类中的各种成分 成员变量 方法 构造器等 如图 二 反射还学什么 加载类 获取类的
  • 做题

    在一个文件中定义一个全局变量 n 主函数 main 在另一个文件中定义函数 fn1 在 main 中对 n 赋值 再调用 fn1 在 fn1 中也对 n 赋值 显示 n 最后的值 include using namespace std in
  • style对象和less/scss互相转换,驼峰转中横线,支持嵌套转换

    不知道有没有小伙伴在维护或重构前端项目时修改样式的时候遇到js style和less scss需要互相转换的问题 本人找网上没有比较完善的转换工具 要么嵌套不支持 要么兼容不好 于是自己写了一个 请大家参考 emotion style In
  • CSDN英雄大会之 SOA技术观感

    5号假装英雄去北京参加了CSDN技术英雄大会 见到了很多一直想见的同行高手还有编辑记者 这里就不一一列举了 只是从SOA中间件开发角度列一下相关的内容 1 IBM如下划分SOA与构件 SOA4类关键构件 基础构件WAS MQ 流程构件WPS
  • 设置一个FreemarkerExceptionHandler捕获freemarker页面上的异常

    在Freemarker页面中如果使用 userName 并且userName为空 那么Freemarker页面就会崩掉 需要设置默认值 userName 来避免对象为空的错误 同理 user userName 也应该写成这样 user us
  • WaveOut播放声音死锁问题原因

    1 首先我们复习下造成死锁的几个充要条件 1 互斥 互斥资源 只能被一个进程使用 2 不剥夺 非抢占式调度 不能强行抢用其他进程资源 3 请求和保持 占有着资源不释放 同时申请其他资源 4 环路等待 没什么可说的 在WaveOutReset
  • java.lang.IncompatibleClassChangeError: Found interface org.elasticsearch.common.bytes.BytesReferenc

    项目场景 再学谷粒商城时 练习elasticsearch时出现一下错误 问题描述 原因分析 提示 出现java lang IncompatibleClassChangeError Found interface org elasticsea
  • redis 由浅入深之 redis.conf配置文件

    是否以后台进程运行 默认为no 如果需要以后台进程运行则改为yes daemonize no 如果以后台进程运行的话 就需要指定pid 你可以在此自定义redis pid文件的位置 pidfile var run redis pid 接受连
  • 2021-04-09

    C 实现XN 2图灵机模拟 算法分析 代码实现 算法分析 1 首先需要完成数的进制转换 将二进制数转化成十进制数 如何转化呢 这里采用除二取余的方法 这种方法就是不断得到余数 放进数组arr 最后除以颠倒一下就行了 2 转换成拓展二进制编码
  • Leetcode刷题-522最长特殊序列II

    题目描述 给定字符串列表 strs 返回它们中最长的特殊序列 如果最长特殊序列不存在 返回 1 最长特殊序列定义 该序列为某字符串独有的最长子序列 即不能是其他字符串的子序列 s 的子序列可以通过删去字符串 s 中的某些字符实现 来源 力扣
  • 灰色预测原理及实例(附代码)

    灰色预测 引言 古人说 凡事预则立 不预则废 办任何事情之前 必须先调查研究 摸清情况 深思熟虑 有科学的预见 周密的计划 这样才能达到预期的成功 所谓预测 就是人们根据可获得的历史和现实数据 资料 运用一定的科学方法和手段 对人类社会 政