灰色预测GM(1,1)-Matlab实现

ops/2024/11/26 6:13:20/

灰色预测是通过少量的，不完全的信息，建立数学模型并给出预测的一种预测方法

是处理小样本预测问题的有效工具

适用条件

 (1)数据是以年份度量的非负数据（如果是月份或者季度数据一般要用时间序列模型）

 (2)数据能通过准指数规律的检验（确定原始数据是否可以使用灰色预测模型）

 (3)数据的期数较短（>3&&<=10）且其他数据之间的关联性不强（如果期数较长，一般用传统的时间序列模型比较合适）

GM(1,1)模型

只含有一个变量的一阶微分方程模型

模型评价

检验模型对原始数据的拟合程度

 代码

gm11脚本（函数文件）

function[result,x0_hat,relative_residuals,eta]=gm11(x0,predict_num)
n=length(x0);
x1=cumsum(x0);
z1=(x1(1:end-1)+x1(2:end))/2;
y=x0(2:end);x=z1;
k=((n-1)*sum(x.*y)-sum(x)*sum(y))/((n-1)*sum(x.*x)-sum(x)*sum(x));
b=(sum(x.*x)*sum(y)-sum(x)*sum(x.*y))/((n-1)*sum(x.*x)-sum(x)*sum(x));
a=-k;
x0_hat=zeros(n,1);x0_hat(1)=x0(1);
for m=1:n-1x0_hat(m+1)=(1-exp(a))*(x0(1)-b/a)*exp(-a*m);
end
result=zeros(predict_num,1);
for i=1:predict_numresult(i)=(1-exp(a))*(x0(1)-b/a)*exp(-a*(n+i-1));
end
absolute_residuals=x0(2:end)-x0_hat(2:end);
relative_residuals=abs(absolute_residuals)./x0(2:end);
class_ratio=x0(2:end)./x0(1:end-1);
eta=abs(1-(1-0.5*a)/(1+0.5*a)*(1./class_ratio));
end

主脚本

%%
%输入原始数据并做出时间序列图
clear;clc
year=[2019:1:2023]';%加'表示转置，将横向量转为列向量
x0=[5503,5222,5429,5119,5175]';%画出原始数据的时间序列图
figure(1);%因为我们的图形不止一个，因此要设置编号
plot(year,x0,'o-');grid on;%plot是用于绘制二维线图的函数；%grid on打开图形的网格线，网格线有助于在图表上更加清晰地看到数据点和趋势
set(gca,'xtick',year(1:1:end))%set(gca, 'xtick', ...) 用于设置当前图形的 x 轴刻度，gca 表示获取当前坐标轴；%year(1:1:end) 表示从 year 向量中提取所有的年份数据。1:1:end 是一种索引方式，表示从第一个元素到最后一个元素，间隔为 1。所以 year(1:1:end) 实际上就是 year 向量的所有元素。
xlabel('Years');ylabel('Number of Pet Cats in China (in 10,000s)');%%
%对一次累加后的数据进行准指数规律的检验
disp('准指数规律检验')
n=length(x0);%计算原始数据的长度
x1=cumsum(x0);%cumsum是累加函数
rho=x0(2:end)./x1(1:end-1);%计算光滑度rho
%画出光滑度的图形，并画上0.5的直线，表示临界值
figure(2)
plot(year(2:end),rho,'o-',[year(2),year(end)],[0.5,0.5],'-');grid on;
text(year(end-1)+0.2,0.55,'临界线')%在坐标(year(end-1)+0.2,0.55)上添加文本
set(gca,'xtick',year(2:1:end))
xlabel('年份');ylabel('原始数据的光滑程度');disp(strcat('指标1：光滑比小于0.5的数据占比为',num2str(100*sum(rho<0.5)/(n-1)),'%'))
disp(strcat('指标2：除去前两个时期外，光滑比小于0.5的数据占比为',num2str(100*sum(rho(3:end)<0.5)/(n-3)),'%'))
disp('参考指标1一般要大于60%，指标2要大于90%，你认为本例数据可以通过检验吗？')
Judge=input('你认为可以通过准指数规律的检验吗？可以通过请输入1，进行进一步的计算：');
%%
%将数据分为训练组和实验组（根据原数据量大小n来取，n<7则取最后2年作为试验组，n>7则取最后3年作为实验组）
if n>7test_num=3;
elsetest_num=2;
end
train_x0=x0(1:end-test_num);%训练数据
disp('训练的数据是：')
disp(mat2str(train_x0'))%mat2str可以将矩阵或向量转换为字符串显示
test_x0=x0(end-test_num+1:end);%试验数据
disp('试验数据是：')
disp(mat2str(test_x0'))
%使用GM(1,1)模型对训练数据进行训练
disp(' ')
result1=gm11(train_x0,test_num);
disp(' ')
%绘制对试验数据进行预测的图形
test_year=year(end-test_num+1:end);
figure(3)
plot(test_year,test_x0,'o-',test_year,result1,'*-');grid on;
set(gca,'xtick',year(end-test_num+1):1:year(end))
legend('试验组的真实数据','GM(1,1)预测结果')
xlabel('Years');ylabel('Number of Pet Cats in China (in 10,000s)')predict_num=input('请输入你要往后面预测的期数：');
[result,x0_hat,relative_residuals,eta]=gm11(x0,predict_num);
%%
%输出使用模型预测出来的结果
disp('对原始数据拟合效果：')
for i=1:ndisp(strcat(num2str(year(i)),': ',num2str(x0_hat(i))))
end
disp(strcat('往后预测',num2str(predict_num),'期的得到的结果：'))
for i=1:predict_numdisp(strcat(num2str(year(end)+i),': ',num2str(result(i))))
end
%%
%绘制相对残差和级比偏差的图形
figure(4)
subplot(2,1,1)
plot(year(2:end),relative_residuals,'*-');grid on;%相对残差relative_residuals
legend('相对残差');xlabel('年份');
set(gca,'xtick',year(2:1:end))
subplot(2,1,2)
plot(year(2:end),eta,'o-');grid on;%级比偏差eta
legend('级比偏差');xlabel('年份');
set(gca,'xtick',year(2:1:end))
disp(' ')
disp('下面将输出对原数据拟合的评价结果')
%%
%残差检验
average_relative_residuals=mean(relative_residuals);
disp(strcat('平均相对残差为',num2str(average_relative_residuals)))
if average_relative_residuals<0.1disp('残差检验的结果表明：该模型对原数据的拟合程度非常不错')
elseif average_relative_residuals<0.2disp('残差检验的结果表明：该模型对原数据的拟合程度达到一般要求')
elsedisp('残差检验的结果表明：该模型对原数据的拟合程度不太好，建议使用其他模型进行预测')
end
%%
%级比偏差检验
average_eta=mean(eta);
disp(strcat('平均级比偏差为',num2str(average_eta)))
if average_eta<0.1disp('级比偏差检验的结果表明：该模型对原数据的拟合程度非常不错')
elseif average_eta<0.2disp('级比偏差检验的结果表明：该模型对原数据的拟合程度达到一般要求')
elsedisp('级比偏差检验的结果表明：该模型对原数据的拟合程度不太好，建议使用其他模型进行预测')
end
%%
%绘制最终的预测效果图
figure(5)
plot(year,x0,'-o',year,x0_hat,'-*m',year(end)+1:year(end)+predict_num,result,'-*b');grid on;%m洋红色 b蓝色
hold on;
plot([year(end),year(end)+1],[x0(end),result(1)],'-*b')
legend('Raw data','Fitting data','Predicted data')
set(gca,'xtick',[year(1):1:year(end)+predict_num])
xlabel('Years');ylabel('Number of Pet Dogs in China (in 10,000s)');


http://www.ppmy.cn/ops/136762.html

相关文章

【C语言】传值调用与传址调用：深度解析与实现

【C语言】传值调用与传址调用：深度解析与实现

博客主页&#xff1a; [小ᶻ☡꙳ᵃⁱᵍᶜ꙳] 本文专栏: C语言 文章目录 &#x1f4af;前言&#x1f4af;什么是传值调用和传址调用&#xff1f;1. 传值调用&#xff08;Call by Value&#xff09;2. 传址调用&#xff08;Call by Reference&#xff09; &#x1f4af;传值调…
阅读更多...
Java基础-Java中的常用类（上）

Java基础-Java中的常用类（上）

(创作不易&#xff0c;感谢有你&#xff0c;你的支持&#xff0c;就是我前行的最大动力&#xff0c;如果看完对你有帮助&#xff0c;请留下您的足迹&#xff09; 目录 String类 创建字符串 字符串长度 连接字符串 创建格式化字符串 String 方法 System类 常用方法 方…
阅读更多...
【6.10】位运算-实现四则运算

【6.10】位运算-实现四则运算

前言 从我们开始上学起就了解到&#xff0c;若要进行加减乘除就要使用对应的运算符号…… 甚至在如今的计算机中也是如此。我们通常只知道如何使用这些运算符号&#xff0c;却很少去关注其底层是如何实现的。假如某天在面试中遇到一道题&#xff0c;要求不使用 “&#xff0c;-…
阅读更多...
Scala中身份证的使用

Scala中身份证的使用

package hfd import scala.util.Random //字符串 //知识点 //1.toInt把字符串转成整数 //2.toUpperCase变大写 //3.toLomerCase变小写 //4.substring(起点&#xff0c;终点&#xff0c;不包括)字符串截取 //5.chartAt(下标)得到对应位置的字符&#xff08;不是字符串&#xff…
阅读更多...
数据库的联合查询

数据库的联合查询

数据库的联合查询 简介为什么要使⽤联合查询多表联合查询时MYSQL内部是如何进⾏计算的构造练习案例数据案例&#xff1a;⼀个完整的联合查询的过程 内连接语法⽰例 外连接语法 ⽰例⾃连接应⽤场景示例表连接练习 ⼦查询语法单⾏⼦查询多⾏⼦查询多列⼦查询在from⼦句中使⽤⼦查…
阅读更多...
极简开源Windows桌面定时提醒休息python程序

极简开源Windows桌面定时提醒休息python程序

当我们长期在电脑面前坐太久后&#xff0c;会产生一系列健康风险&#xff0c;包括干眼症&#xff0c;颈椎&#xff0c;腰椎&#xff0c;肌肉僵硬等等。解决方案是在一定的时间间隔内我们需要have a break, 远眺可以缓解干眼症等眼部症状&#xff0c;站起来走动两步&#xff0c;…
阅读更多...
一键部署zabbix-agent2的脚本

一键部署zabbix-agent2的脚本

1、首先下载agent2的安装包 我是X86的centos 7系统&#xff0c;zabbix-agent2-5.0.42-1.el7.x86_64.rpm下载地址&#xff0c;另外很多国产系统统信、中科方德也适用这个版本。 这个网站里面有其他版本的&#xff0c;自行选择下载 https://repo.zabbix.com/zabbix/5.0/rhel/7/…
阅读更多...
jmeter5.6.3安装教程

jmeter5.6.3安装教程

一、官网下载 需要提前配置好jdk的环境变量 jmeter官网&#xff1a;https://jmeter.apache.org/download_jmeter.cgi 选择点击二进制的zip文件 下载成功后&#xff0c;默认解压下一步&#xff0c;更改安装路径就行(我安装在D盘) 实用jmeter的bin目录作为系统变量 然后把这…
阅读更多...
最新文章

Copyright @ 2022~2023