草庐IT

python数学建模--灰色预测

夺笋123 2023-07-23 原文

目录

灰色预测

灰色预测模型是通过少量的、不完全的信息,建立数学模型做出预测的预测方法,基于客观事物的过去和现在的发展规律,借助于科学的方法对未来的发展趋势和状况进行描述和分析,并形成科学的假设和判断。

一阶灰色方程GM(1,1)

对于原始数据序列 x ( 0 ) = [ x ( 0 ) ( 1 ) , x ( 0 ) ( 2 ) , . . . , x ( 0 ) ( n ) ] x^{(0)}=[x^{(0)}(1),x^{(0)}(2),...,x^{(0)}(n)] x(0)=[x(0)(1),x(0)(2),...,x(0)(n)]

建模步骤

  1. 计算级比
    已知序列的级比为
    σ ( k ) = x ( k ) x ( k − 1 ) , k = ( 2 , 3 , . . . , n ) \sigma(k)=\frac{x(k)}{x(k-1)},k=(2,3,...,n) σ(k)=x(k1)x(k),k=(2,3,...,n)
    若所有级比 σ ( k ) \sigma(k) σ(k)都落在区间 ( e − 2 n + 1 , e 2 n + 1 ) (e^{-\frac{2}{n+1}},e^{\frac{2}{n+1}}) (en+12,en+12)内,则序列 x ( 0 ) x^{(0)} x(0)可使用GM(1,1)进行灰色预测,否则应对序列数据 x ( 0 ) x^{(0)} x(0)进行平移变换
    x ( 0 ) ( k ) = x ( 0 ) ( k ) + c , k = 1 , 2 , . . . , n x^{(0)}(k)=x^{(0)}(k)+c,k=1,2,...,n x(0)(k)=x(0)(k)+c,k=1,2,...,n

遗留的问题:原始序列进行累加后为什么不会影响到最终结果?

  1. 累加序列及其均值生成序列
    一次累加序列
    x ( 1 ) = ( x ( 1 ) ( 1 ) , x ( 1 ) ( 2 ) , . . . , x ( 1 ) ( 7 ) ) = ( x ( 0 ) ( 1 ) , x ( 0 ) ( 1 ) + x ( 0 ) ( 2 ) , . . . , x ( 0 ) ( 1 ) + . . . + x ( 0 ) ( 7 ) ) \begin{aligned} x^{(1)}&=(x^{(1)}(1),x^{(1)}(2),...,x^{(1)}(7))\\ &=(x^{(0)}(1),x^{(0)}(1)+x^{(0)}(2),...,x^{(0)}(1)+...+x^{(0)}(7)) \end{aligned} x(1)=(x(1)(1),x(1)(2),...,x(1)(7))=(x(0)(1),x(0)(1)+x(0)(2),...,x(0)(1)+...+x(0)(7))
    一次累加序列的均值生成序列为
    z ( 1 ) = ( z ( 1 ) ( 2 ) , z ( 1 ) ( 3 ) , . . . , z ( 1 ) ( 7 ) ) 其中 z ( 1 ) ( k ) = α x ( 1 ) ( k ) + α x ( 1 ) ( k − 1 ) , k = ( 2 , 3 , . . . , 7 ) z^{(1)}=(z^{(1)}(2),z^{(1)}(3),...,z^{(1)}(7))\\ 其中z^{(1)}(k)=\alpha x^{(1)}(k)+\alpha x^{(1)}(k-1),k=(2,3,...,7) z(1)=(z(1)(2),z(1)(3),...,z(1)(7))其中z(1)(k)=αx(1)(k)+αx(1)(k1),k=(2,3,...,7)

注:均值生成序列中的权重参数取值范围为 0 ≤ α ≤ 1 0\leq \alpha \leq 1 0α1,一般取0.5

  1. 建立模型
    白化方程
    d x ( 1 ) ( t ) d t − a z ( 1 ) ( t ) = b \frac{dx^{(1)}(t)}{dt}-az^{(1)}(t)=b dtdx(1)(t)az(1)(t)=b
    灰微分方程
    x ( 0 ) ( k ) − a z k ( 1 ) = b 即 − a z k ( 1 ) + b = x ( 0 ) ( k ) (1) x^{(0)}(k)-az_k^{(1)}=b\\ 即-az_k^{(1)}+b=x^{(0)}(k)\tag{1} x(0)(k)azk(1)=bazk(1)+b=x(0)(k)(1)

其实在这个地方还是有一些遗留问题的,比如白化方程、灰微分方程各是怎么得到的?建立这两个方程背后的数学原理?

将(1)式化为矩阵形式为

[ − z ( 1 ) ( 2 ) 1 − z ( 1 ) ( 3 ) 1 ⋮ 1 − z ( 1 ) ( k ) 1 ] [ a b ] = [ x ( 0 ) ( 2 ) x ( 0 ) ( 3 ) ⋮ x ( 0 ) ( k ) ] \begin{bmatrix} -z^{(1)}(2)&1\\ -z^{(1)}(3)&1\\ \vdots&1\\ -z^{(1)}(k)&1\\ \end{bmatrix} \begin{bmatrix} a\\b \end{bmatrix}=\begin{bmatrix} x^{(0)}(2)\\ x^{(0)}(3)\\ \vdots\\ x^{(0)}(k) \end{bmatrix} z(1)(2)z(1)(3)z(1)(k)1111 [ab]= x(0)(2)x(0)(3)x(0)(k)
X u = Y Xu=Y Xu=Y,根据最小二乘法得到参数矩阵 u u u的估计值为
u = ( X T X ) − 1 X T Y u=(X^TX)^{-1}X^TY u=(XTX)1XTY

代入数据求得参数矩阵 u u u,由此得到白化方程与灰微分方程的参数 a , b a,b a,b
那么白化方程的时间响应式:
x ( 1 ) ( k + 1 ) = ( x ( 0 ) ( 1 ) − b a ) e − a k + b a ,其中 k = 0 , 1 , . . . x^{(1)}(k+1)=(x^{(0)}(1)-\frac{b}{a})e^{-ak}+\frac{b}{a},其中k=0,1,... x(1)(k+1)=(x(0)(1)ab)eak+ab,其中k=0,1,...
将k代入上式,得到一次累减序列 x m i n u s ( 1 ) = [ x m i n u s ( 1 ) ( 1 ) , x m i n u s ( 1 ) ( 2 ) , . . . , x m i n u s ( 1 ) ( n ) ] x^{(1)}_{minus}=[x^{(1)}_{minus}(1),x_{minus}^{(1)}(2),...,x_{minus}^{(1)}(n)] xminus(1)=[xminus(1)(1),xminus(1)(2),...,xminus(1)(n)]

那么预测序列为
x p r e d ( 0 ) = [ x p r e d ( 0 ) ( 1 ) , x p r e d ( 0 ) ( 2 ) , . . . , x p r e d ( 0 ) ( k ) ] , k = 1 , 2 , . . . 其中 x p r e d ( 0 ) ( k ) = x m i n u s ( 0 ) ( k ) − x m i n u s ( 0 ) ( k − 1 ) x p r e d ( 0 ) ( 1 ) = x m i n u s ( 1 ) ( 1 ) x^{(0)}_{pred}=[x^{(0)}_{pred}(1),x^{(0)}_{pred}(2),...,x^{(0)}_{pred}(k)],k=1,2,...\\其中x^{(0)}_{pred}(k)=x^{(0)}_{minus}(k)-x^{(0)}_{minus}(k-1)\\ x^{(0)}_{pred}(1)=x^{(1)}_{minus}(1) xpred(0)=[xpred(0)(1),xpred(0)(2),...,xpred(0)(k)],k=1,2,...其中xpred(0)(k)=xminus(0)(k)xminus(0)(k1)xpred(0)(1)=xminus(1)(1)
5. 误差检验
(1)相对误差检验
λ ( k ) = ∣ x ( 0 ) ( k ) − x p r e d ( 0 ) ( k ) ∣ x ( 0 ) ( k ) , k = 1 , 2 , . . . , n \lambda(k)=\frac{|x^{(0)}(k)-x^{(0)}_{pred}(k)|}{x^{(0)}(k)},k=1,2,...,n λ(k)=x(0)(k)x(0)(k)xpred(0)(k),k=1,2,...,n

λ ( k ) < 0.2 时 \lambda (k)<0.2时 λ(k)<0.2,可认为模型达到了一般要求
λ ( k ) < 0.1 时 \lambda (k)<0.1时 λ(k)<0.1,可认为模型达到了较高要求
(2)级比偏差值检验
ϵ ( k ) = ∣ 1 − ( 1 − 0.5 a 1 + 0.5 a ) σ ( k ) ∣ , k = 2 , 3 , . . . , n \epsilon(k)=|1-(\frac{1-0.5a}{1+0.5a})\sigma (k)|,k=2,3,...,n ϵ(k)=∣1(1+0.5a10.5a)σ(k),k=2,3,...,n

ϵ ( k ) < 0.2 时 \epsilon(k)<0.2时 ϵ(k)<0.2,可认为模型达到了一般要求
ϵ ( k ) < 0.1 时 \epsilon(k)<0.1时 ϵ(k)<0.1,可认为模型达到了较高要求

应用及其求解步骤

已知 x ( 0 ) = ( 25723 , 30379 , 34473 , 38485 , 40514 , 42400 , 48337 ) , 试建立 G M ( 1 , 1 ) 模型 已知x^{(0)}=(25723,30379,34473,38485,40514,42400,48337),试建立GM(1,1)模型 已知x(0)=(25723,30379,34473,38485,40514,42400,48337),试建立GM(1,1)模型
建模开始

求级比

import numpy as np
x0=np.array([25723,30379,34473,38485,40514,42400,48337])
sigma=x0[:-1]/x0[1:]
>>> array([0.84673623, 0.88124039, 0.89575159, 0.94991855, 0.95551887,0.87717484])
bound1=[sigma.min(),sigma.max()]
>>> [0.8467362322657098, 0.9555188679245283]
bound2=[np.exp(-2/(x0.size+1)),np.exp(2/(x0.size+1))]
>>> [0.7788007830714049, 1.2840254166877414]

由于bound1 ∈ \in bound2,故级比符合要求

一次累加序列

x ( 1 ) = ( x ( 1 ) ( 1 ) , x ( 1 ) ( 2 ) , . . . , x ( 1 ) ( 7 ) ) = ( x ( 0 ) ( 1 ) , x ( 0 ) ( 1 ) + x ( 0 ) ( 2 ) , . . . , x ( 0 ) ( 1 ) + . . . + x ( 0 ) ( 7 ) ) \begin{aligned} x^{(1)}&=(x^{(1)}(1),x^{(1)}(2),...,x^{(1)}(7))\\ &=(x^{(0)}(1),x^{(0)}(1)+x^{(0)}(2),...,x^{(0)}(1)+...+x^{(0)}(7))\\ \end{aligned} x(1)=(x(1)(1),x(1)(2),...,x(1)(7))=(x(0)(1),x(0)(1)+x(0)(2),...,x(0)(1)+...+x(0)(7))

x1=np.cumsum(x0)
>>> array([ 25723,  56102,  90575, 129060, 169574, 211974, 260311],dtype=int32)

一次累加序列的均值生成序列为
z ( 1 ) = ( z ( 1 ) ( 2 ) , z ( 1 ) ( 3 ) , . . . , z ( 1 ) ( 7 ) ) 其中 z ( 1 ) ( k ) = 0.5 x ( 1 ) ( k ) + 0.5 x ( 1 ) ( k − 1 ) , k = ( 2 , 3 , . . . , 7 ) z^{(1)}=(z^{(1)}(2),z^{(1)}(3),...,z^{(1)}(7))\\ 其中z^{(1)}(k)=0.5x^{(1)}(k)+0.5x^{(1)}(k-1),k=(2,3,...,7) z(1)=(z(1)(2),z(1)(3),...,z(1)(7))其中z(1)(k)=0.5x(1)(k)+0.5x(1)(k1),k=(2,3,...,7)

z1=(x1[:-1]+x1[1:])/2.0
>>> array([ 40912.5,  73338.5, 109817.5, 149317. , 190774. , 236142.5])

求参数矩阵 u u u

X = [ − z ( 1 ) ( 2 ) 1 − z ( 1 ) ( 3 ) 1 ⋮ 1 − z ( 1 ) ( k ) 1 ] Y = [ x ( 0 ) ( 2 ) x ( 0 ) ( 3 ) ⋮ x ( 0 ) ( k ) ] X=\begin{bmatrix} -z^{(1)}(2)&1\\ -z^{(1)}(3)&1\\ \vdots&1\\ -z^{(1)}(k)&1\\ \end{bmatrix}Y=\begin{bmatrix} x^{(0)}(2)\\ x^{(0)}(3)\\ \vdots\\ x^{(0)}(k) \end{bmatrix} X= z(1)(2)z(1)(3)z(1)(k)1111 Y= x(0)(2)x(0)(3)x(0)(k)

X=np.vstack([-z1,np.ones(x0.size-1)]).T
Y=x0[1:].T
u=np.matmul(np.matmul((np.linalg.inv(np.matmul(X.T,X))),X.T),Y)
>>> array([-8.42648092e-02,  2.78584508e+04])

时间响应式

于是,白化方程的时间响应式为
x ( 1 ) ( k + 1 ) = ( x ( 0 ) ( 1 ) − b a ) e − a k + b a = ( 25723 − 27858.45 − 0.0843 ) e 0.0843 k − 27858.45 − 0.0843 = 356328.99 e 0.0843 k − 330605.99 \begin{aligned} x^{(1)}(k+1)&=(x^{(0)}(1)-\frac{b}{a})e^{-ak}+\frac{b}{a}\\ &=(25723-\frac{27858.45}{-0.0843})e^{0.0843k}-\frac{27858.45}{-0.0843}\\ &=356328.99e^{0.0843k}-330605.99 \end{aligned} x(1)(k+1)=(x(0)(1)ab)eak+ab=(257230.084327858.45)e0.0843k0.084327858.45=356328.99e0.0843k330605.99

求预测序列

#生成匿名函数
x_m=lambda k:(x0[0]-(u[1]/u[0]))*np.e**(-u[0]*k)+u[1]/u[0]
# 一次累减序列
x_m1=x_m(np.arange(x0.size+1))
# 预测序列
x_pred=np.diff(x_m1)
>>> array([31327.3567122 , 34081.56224489, 37077.90911705, 40337.68565577,
       43884.05179286, 47742.20361062, 51939.55235392])

模型检验

相对误差
λ ( k ) = ∣ x ( 0 ) ( k ) − x p r e d ( 0 ) ( k ) ∣ x ( 0 ) ( k ) , k = 1 , 2 , . . . , n \lambda(k)=\frac{|x^{(0)}(k)-x^{(0)}_{pred}(k)|}{x^{(0)}(k)},k=1,2,...,n λ(k)=x(0)(k)x(0)(k)xpred(0)(k),k=1,2,...,n

lambda_k=abs(x0[1:]-x_pred[:-1])/x0[1:]
>>> array([0.03121751, 0.01135491, 0.03656206, 0.00435194, 0.03500122,
       0.0123052 ])

相对误差值 λ ( k ) < 0.1 \lambda(k)<0.1 λ(k)<0.1,可认为模型达到了较高要求

级比偏差值
ϵ ( k ) = ∣ 1 − ( 1 − 0.5 a 1 + 0.5 a ) σ ( k ) ∣ , k = 2 , 3 , . . . , n \epsilon(k)=|1-(\frac{1-0.5a}{1+0.5a})\sigma (k)|,k=2,3,...,n ϵ(k)=∣1(1+0.5a10.5a)σ(k),k=2,3,...,n

epsilon_k=abs(1-((1-0.5*(-0.0834))/(1+0.5*(-0.0834)))*sigma)
>>> array([0.07957306, 0.04206604, 0.02629194, 0.03258912, 0.03867683,
       0.04648542])

级比偏差值 ϵ ( k ) < 0.1 \epsilon(k)<0.1 ϵ(k)<0.1,可认为模型达到了较高要求

实际值与预测值比较及可视化

import matplotlib
import matplotlib.pyplot as plt
x_pred1=x_pred.tolist()
x_pred1.insert(0,x0[0])
matplotlib.rcParams['font.family']='SimHei'
plt.plot(range(x0.size),x0,'r-',label='原始序列')
plt.plot(range(x0.size),x_pred1[:-1],'g--',label='预测序列')
plt.legend()
plt.show()

运行结果

二阶灰色方程GM(2,1)

To be continue…

有关python数学建模--灰色预测的更多相关文章

  1. python - 如何使用 Ruby 或 Python 创建一系列高音调和低音调的蜂鸣声? - 2

    关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题?更新问题,以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点?我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等),但我确实想创建一个输出文件。

  2. ruby-on-rails - 建模收藏夹 - 2

    我希望将Favorite模型添加到我的User和Link模型。业务逻辑用户可以有多个链接(即可以添加多个链接)用户可以收藏多个链接(他们自己的或其他用户的)一个链接可以被多个用户收藏,但只有一个所有者我对如何为这种关联建模以及在模型就位后如何创建用户收藏夹感到困惑?classUser 最佳答案 下面的数据模型怎么样:classUser:destroyhas_many:favorite_links,:through=>:favorites,:source=>:linkendclassLink:destroyhas_many:favor

  3. Python 相当于 Perl/Ruby ||= - 2

    这个问题在这里已经有了答案:关闭10年前。PossibleDuplicate:Pythonconditionalassignmentoperator对于这样一个简单的问题表示歉意,但是谷歌搜索||=并不是很有帮助;)Python中是否有与Ruby和Perl中的||=语句等效的语句?例如:foo="hey"foo||="what"#assignfooifit'sundefined#fooisstill"hey"bar||="yeah"#baris"yeah"另外,类似这样的东西的通用术语是什么?条件分配是我的第一个猜测,但Wikipediapage跟我想的不太一样。

  4. java - 什么相当于 ruby​​ 的 rack 或 python 的 Java wsgi? - 2

    什么是ruby​​的rack或python的Java的wsgi?还有一个路由库。 最佳答案 来自Python标准PEP333:Bycontrast,althoughJavahasjustasmanywebapplicationframeworksavailable,Java's"servlet"APImakesitpossibleforapplicationswrittenwithanyJavawebapplicationframeworktoruninanywebserverthatsupportstheservletAPI.ht

  5. 华为OD机试用Python实现 -【明明的随机数】 2023Q1A - 2

    华为OD机试题本篇题目:明明的随机数题目输入描述输出描述:示例1输入输出说明代码编写思路最近更新的博客华为od2023|什么是华为od,od薪资待遇,od机试题清单华为OD机试真题大全,用Python解华为机试题|机试宝典【华为OD机试】全流程解析+经验分享,题型分享,防作弊指南华为o

  6. python - 如何读取 MIDI 文件、更改其乐器并将其写回? - 2

    我想解析一个已经存在的.mid文件,改变它的乐器,例如从“acousticgrandpiano”到“violin”,然后将它保存回去或作为另一个.mid文件。根据我在文档中看到的内容,该乐器通过program_change或patch_change指令进行了更改,但我找不到任何在已经存在的MIDI文件中执行此操作的库.他们似乎都只支持从头开始创建的MIDI文件。 最佳答案 MIDIpackage会为您完成此操作,但具体方法取决于midi文件的原始内容。一个MIDI文件由一个或多个音轨组成,每个音轨是十六个channel中任何一个上的

  7. 「Python|Selenium|场景案例」如何定位iframe中的元素? - 2

    本文主要介绍在使用Selenium进行自动化测试或者任务时,对于使用了iframe的页面,如何定位iframe中的元素文章目录场景描述解决方案具体代码场景描述当我们在使用Selenium进行自动化测试的时候,可能会遇到一些界面或者窗体是使用HTML的iframe标签进行承载的。对于iframe中的标签,如果直接查找是无法找到的,会抛出没有找到元素的异常。比如近在咫尺的例子就是,CSDN的登录窗体就是使用的iframe,大家可以尝试通过F12开发者模式查看到的tag_name,class_name,id或者xpath来定位中的页面元素,会抛出NoSuchElementException异常。解决

  8. python ffmpeg 使用 pyav 转换 一组图像 到 视频 - 2

    2022/8/4更新支持加入水印水印必须包含透明图像,并且水印图像大小要等于原图像的大小pythonconvert_image_to_video.py-f30-mwatermark.pngim_dirout.mkv2022/6/21更新让命令行参数更加易用新的命令行使用方法pythonconvert_image_to_video.py-f30im_dirout.mkvFFMPEG命令行转换一组JPG图像到视频时,是将这组图像视为MJPG流。我需要转换一组PNG图像到视频,FFMPEG就不认了。pyav内置了ffmpeg库,不需要系统带有ffmpeg工具因此我使用ffmpeg的python包装p

  9. Python 刷Leetcode题库,顺带学英语单词(31) - 2

    ValidPalindromeGivenastring,determineifitisapalindrome,consideringonlyalphanumericcharactersandignoringcases. [#125]Example:"Aman,aplan,acanal:Panama"isapalindrome."raceacar"isnotapalindrome.Haveyouconsiderthatthestringmightbeempty?Thisisagoodquestiontoaskduringaninterview.Forthepurposeofthisproblem

  10. python - 是否可以使用 Ruby 或 Python 禁用 anchor /引用来发出有效的 YAML? - 2

    是否可以在PyYAML或Ruby的Psych引擎中禁用创建anchor和引用(并有效地显式列出冗余数据)?也许我在网上搜索时遗漏了一些东西,但在Psych中似乎没有太多可用的选项,而且我也无法确定PyYAML是否允许这样做.基本原理是我必须序列化一些数据并将其以可读的形式传递给一个不是真正的技术同事进行手动验证。有些数据是多余的,但我需要以最明确的方式列出它们以提高可读性(anchor和引用是提高效率的好概念,但不是人类可读性)。Ruby和Python是我选择的工具,但如果有其他一些相当简单的方法来“展开”YAML文档,它可能就可以了。 最佳答案

随机推荐