爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 8555|回复: 10

[源程序] Bayesian时间序列突变点分析,趋势分析,和季节去除分解 (Matlab, Python,和R)

[复制链接]

新浪微博达人勋

发表于 2022-6-10 13:26:01 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 zhaokg 于 2024-2-24 08:02 编辑

借贵宝地分享自己开发的一个贝叶斯间序列分析算法BEAST。它主要用于把时间序列分解为趋势(trend)和季节周期波动(seasonality), 同时对这2项也进行突变点检测。如果时间序列没有季节项,也可以只进行趋势分析和突变点检测 (changepoint detection and trend analysis)。  这个算法和程序已经被一些文章和教材用过或者介绍过,其中包括一篇文献用它来分析降雨数据(White, J.H., Walsh, J.E. and Thoman Jr, R.L., 2021. Using Bayesian statistics to detect trends in Alaskan precipitation. International Journal of Climatology, 41(3), pp.2045-2059.)我觉的BEAST可能对分析某些气象数据或许有些用处。 这个算法在这个网址有更多介绍 https://github.com/zhaokg/Rbeast.

相比传统的突变点分析方法(比如bfast或者是M-K),它的优势是不仅告诉那儿发生突变,并且给出每个时间点发出突变的概率;也可以告诉发生几个突变点的概率(比如有1个突变点的概率,有2个的概率,....)。


这个算法可以在R, Matlab,或者 Python中调用。在R中,已经做成一个package,可以直接用install.packages("Rbeast")安装

  1. install.packages("Rbeast")
  2. library(Rbeast)
  3. data(Nile)                                       #  annual streamflow of the Nile River   
  4. out = beast(Nile, season='none')      #  'none': trend-only data without seasonlaity   
  5. print(out)                  
  6. plot(out)
复制代码




在Matlab中,可以不用手动下载, 直接可以用下面前2行代码在线安装 (如果下面自动download 不work, 可以直接去https://github.com/zhaokg/Rbeast/tree/master/Matlab手动下载:
  1. beastPath = 'c:\beast'                                                 % 本地安装目录
  2. eval(webread('http://b.link/rbeast',weboptions('cert','')))  % 自动安装在beastPath指定目录
  3. load('Nile.mat')                                                          % 尼罗河年径流量
  4. out = beast(Nile, 'season', 'none','start', 1871)               % trend-only data without seasonality
  5. printbeast(out)
  6. plotbeast(out)
  7. %下图为这个例子的结果
  8. %第一个截屏是 突变点个数的概率,比如说有54%的概率存在1个突变点
  9. %第二个截屏是 画出的主要结果
复制代码


NileProb.png



                               
登录/注册后可看大图


Matlab中另外一个例子

  1. load('googletrend.mat')               % 美国自2004年每月在Google搜索'beach'关键词的频度:这个时间序列有一个周期项
  2. o = beast( beach )                       % Apply BEAST to the 'beach time series: beach is a data vector only; the time
  3.                                                  % info can be supplied using the start and deltat
  4.                                                  % keywords, as in the next commented line
  5. %o = beast( beach,'start',[2004,1],'deltat',1/12 )
  6. printbeast(o)          % print the changepoints detected
  7. plotbeast(o)            % plot the results: o.season.Y and o.trend.Y are the seasonal and trend compoents
复制代码


                               
登录/注册后可看大图

在Python中,已经做成一个package放在https://pypi.org/project/Rbeast/, 可以直接用pip install Rbeast安装:


  1. import Rbeast as rb
  2. nile, year = rb.load_example('nile')
  3. o            = rb.beast( nile, start=1871, season='none')
  4. rb.plot(o, title='Annual streamflow of the Nile River')
  5. rb.print(o)
复制代码


密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2022-6-30 20:28:22 | 显示全部楼层
感谢楼主介绍,赞一个
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2022-8-9 18:27:07 | 显示全部楼层
谢谢分享,学习了
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2022-6-11 13:31:16 | 显示全部楼层
感谢楼主介绍,赞一个
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-1-15 16:03:25 | 显示全部楼层
感谢楼主分享!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-6-15 09:43:46 | 显示全部楼层
matlab 的安装链接安装不上去
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2023-6-15 21:34:01 | 显示全部楼层
本帖最后由 zhaokg 于 2024-9-11 12:39 编辑
deemo7 发表于 2023-6-15 09:43
matlab 的安装链接安装不上去

谢谢,我在我这边试了下,可以看装。请务必选一个你有write 权限的安装路径. 如果不行的话, 请到https://github.com/zhaokg/Rbeast/tree/master/Matlab手动下载所有文件。

                               
登录/注册后可看大图

Untitled.png

密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-6-16 07:30:46 | 显示全部楼层
谢谢
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2023-7-19 21:26:17 | 显示全部楼层
感谢楼主,这个怎么保存成图片呢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2023-8-17 07:00:57 | 显示全部楼层
zq20180902 发表于 2023-7-19 21:26
感谢楼主,这个怎么保存成图片呢

刚看到您的贴子, 如果在matlab中的话, 用matlab自带的exportgraphics可以。
  1. load('Nile.mat')                                                           
  2. out  = beast(Nile, 'season')
  3. plotbeast(out)
  4. exportgraphics(gca, 'results.pdf')
  5. exportgraphics(gca, 'results.png')
复制代码


如果需要做自己customized的图片的话, 可以i用输出 out中的结果自己做i图。 不知有没有回答您的问题。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表