R中的极值分析入门–第3部分:峰值阈值法

嘿,

welcome to part 3 of our short introduction to 极值分析 using the extRemes package in R.

讨论了块最大值方法 最后一次,我们现在来看一下阈值超出方法。

根据Coles(2001),如果没有间隙的完整(时间)序列可用,则阈值方法比块极大方法更有效,因为所有超过特定阈值的值都可以用作模型拟合的基础。在某些情况下,将分布拟合到块最大值数据是一种浪费的方法,因为每个块仅使用一个值进行建模,而阈值过量方法则可能提供有关极限的更多信息。

但是,类似于在块最大值方法中选择块大小,对于部分持续时间模型的阈值的选择也要在偏差(低阈值)和方差(高阈值)之间进行权衡。

Coles(2001)描述了两种不同的阈值选择方法。首先,有一种基于平均剩余寿命图的探索方法。在实际模型拟合之前应用此技术。其次,另一种方法是评估参数估计的稳定性。因此,该模型拟合的敏感性分析是在一系列不同阈值范围内执行的。

但是,选择适当的阈值可能是使用部分持续时间序列执行极值分析的最关键部分。 Scarrott和MacDonald在其2012年的文章中对阈值估计的方法进行了很好的概述 极值阈值估计和不确定性量化的回顾 (REVSTAT 10(1):33–59)。

找到合适的阈值后,超过此阈值的极值的最终子集将用于拟合广义的Pareto分布。

根据Pickands-Balkema-de Haan定理,超过阈值的值的分布可以近似为 广义pareto分布.

The following code shows a short practical example of fitting a 广义pareto分布 to a time series of precipitation 数据 using the extRemes package in R. The sample 数据 set features precipitation 数据 of Oberwang (Salzburg, Austria) from 1981 through 2014 and is is provided by the hydrographic services of Austria via via H. The function read_ehyd() for importing the 数据 set can be found at 使用R从eHYD读取数据.

# load packages
library(extRemes)
library(xts)

# get 数据 from H
ehyd_url <- "http://ehyd.gv.at/eHYD/MessstellenExtraData/nlv?id=105700&file=2"
precipitation_xts <- read_ehyd(ehyd_url)

# mean residual life plot:
mrlplot(precipitation_xts, main="Mean Residual Life Plot")
# The mean residual life plot depicts the Thresholds (u) vs Mean Excess flow.
# The idea is to find the lowest threshold where the plot is nearly linear;
# taking into account the 95% confidence bounds.

# fitting the GPD model over a range of thresholds
threshrange.plot(precipitation_xts, r = c(30, 45), nint = 16)
# ismev implementation is faster:
# ismev::gpd.fitrange(precipitation_xts, umin=30, umax=45, nint = 16)
# set threshold
th <- 40

# maximum likelihood estimation
pot_mle <- fevd(as.vector(precipitation_xts), method = "MLE", type="GP", threshold=th)
# diagnostic plots
plot(pot_mle)
rl_mle <- return.level(pot_mle, conf = 0.05, return.period= c(2,5,10,20,50,100))

# L-moments estimation
pot_lmom <- fevd(as.vector(precipitation_xts), method = "Lmoments", type="GP", threshold=th)
# diagnostic plots
plot(pot_lmom)
rl_lmom <- return.level(pot_lmom, conf = 0.05, return.period= c(2,5,10,20,50,100))

# return level plots
par(mfcol=c(1,2))
# return level plot w/ MLE
plot(pot_mle, type="rl",
     main="Return Level Plot for Oberwang w/ MLE",
     ylim=c(0,200), pch=16)
loc <- as.numeric(return.level(pot_mle, conf = 0.05,return.period=100))
segments(100, 0, 100, loc, col= 'midnightblue',lty=6)
segments(0.01,loc,100, loc, col='midnightblue', lty=6)

# return level plot w/ LMOM
plot(pot_lmom, type="rl",
     main="Return Level Plot for Oberwang w/ L-Moments",
     ylim=c(0,200))
loc <- as.numeric(return.level(pot_lmom, conf = 0.05,return.period=100))
segments(100, 0, 100, loc, col= 'midnightblue',lty=6)
segments(0.01,loc,100, loc, col='midnightblue', lty=6)

# comparison of return levels
results <- t(data.frame(mle=as.numeric(rl_mle),
                        lmom=as.numeric(rl_lmom)))
colnames(results) <- c(2,5,10,20,50,100)
round(results,1)

return_levels_pot
This is an example that illustrates nicely why the approach based on L-moments may be preferred over maximum likelihood estimation, as the right image proves clearly that the influence of the outlier is much smaller when using L-moments estimation. In addition to these classical estimation methods, extRemes offers Generalized Maximum Likelihood Estimation (GMLE, 马丁斯和施丁格,2000年)和贝叶斯估计方法(吉兰(Gilleland)和卡兹(Katz),2016年)。

参考:Coles,Stuart(2001)。极值统计建模简介。统计中的Springer系列。伦敦:施普林格。

关于作者

马蒂亚斯在维也纳自然资源与生命科学大学学习了环境信息管理,并获得了环境统计博士学位。他的论文的重点是罕见(极端)事件的统计建模,作为对关键基础设施进行漏洞评估的基础。他目前在奥地利国家气象和地球物理服务局(ZAMG)和BOKU大学山区风险工程研究所工作。他目前专注于(统计)不良天气事件和自然灾害以及减少灾害风险的评估。他的主要兴趣是环境现象的统计建模以及用于数据科学,地理信息和遥感的开源工具。

9条留言

您可以在这篇文章中发表评论。


  • 恭喜!这真的很有帮助。
    我正在研究R和Ι的极端值,已经使用以下方法估计了GEV和Pareto分布的参数(比例,形状)“MLE” and “L moments”方法。此外,我必须使用贝叶斯方法计算GEV分布的形状和比例值。但是我可以’使用贝叶斯方法估计这些参数以进行帕累托分布。更具体地说,我在写
    fevd(data, threshold, type="GP", method="Bayesian")
    我正在从R接收以下消息:
    Error in colMeans(x$chain.info[, 2:(pdim[2] - 1)], na.rm = TRUE):
    'x' must be an array of at least two dimensions
    .

    你能帮我解决我的问题吗?
    提前致谢!

    佐治亚州L 在4年前 回复


    • 嘿,
      重现此错误消息有点困难。我想这不是与extRemes包有关的问题,而是与您的输入数据集有关的问题。
      您的输入数据数据是什么样的?

      使用贝叶斯参数估计的最小工作示例为:

      library(extRemes)
      library(ismev)
      
      # Southwest-England rainfall 数据 from Coles
      data(rain)
      
      # simple threshold (usually one should not use fixed quantiles)
      u <- quantile(rain, 0.9)
      
      # fit GP using Bayesian parameter estimation
      x <- fevd(rain, threshold = u , type='GP', method='Bayesian')
      

      这需要几秒钟,但效果很好。

      问候,
      马蒂亚斯

      马蒂亚斯 在4年前 回复


  • 嗨,谢谢你的这篇文章。非常清晰易懂。

    I’m working on extreme event in R using fevd {extRemes}. I have a big sample (more than 4 million rows). I’我试图拟合我的数据,但我发现了一个巨大的AIC,但它并没有’与情节相符。较大的AIC更适合我的数据,对于最小的AIC,返回期为空。所以我的问题是:我应该使用绘图还是AIC选择模型?

    谢谢你的帮助

    我的R代码:

    B1.fit=fevd(B1[,2], 数据=B1, threshold=B1.th, type="GEV",method ="MLE", units="KVA",time.units="seconds",period = "Seconds") 
    B1.fit1=fevd(B1[,2], 数据=B1, threshold=B1.th,type="GP",method ="MLE", units="KVA",time.units="seconds",period = "Seconds") 
    B1.fit2=fevd(B1[,2], 数据=B1, threshold=B1.th,type="Gumbel", method ="MLE",units="KVA",time.units="seconds",period = "Seconds")
    B1.fit3=fevd(B1[,2], 数据=B1, threshold=B1.th,type="Exponential",method ="MLE", units="KVA",time.units="seconds",period = "Seconds")
    
    fit.AIC=summary(B1.fit, silent=TRUE)$AIC
    fit1.AIC=summary(B1.fit1, silent=TRUE)$AIC
    fit2.AIC=summary(B1.fit2, silent=TRUE)$AIC
    fit3.AIC=summary(B1.fit3, silent=TRUE)$AIC
    
    fit.AIC
    # [1] 39976258
    fit1.AIC
    # [1] 466351.5
    fit2.AIC
    # [1] 13934878
    fit3.AIC
    # [1] 466330.8 
    
    plot(B1.fit)
    plot(B1.fit1)
    plot(B1.fit2)
    plot(B1.fit3)
    

    费里尔·B 在4年前 回复


    • 你好

      在不知道您的数据集的情况下,我只能提供以下建议:

      • 当采用阈值过量方法时,阈值选择可以说是最关键的部分。确保参考诊断图(参数稳定性图,平均过量图)以选择合适的阈值。
      • 您似乎只使用一个数据集。这应该使您可以仔细查询所产生的拟合功能的各种诊断图。我个人认为结果的目视检查是一个非常重要的方面。拟合优度度量可能会导致令人困惑的结论,尤其是在阈值超出方法中。一个模型可以很好地拟合分布的较低尾部的大量数据,但是在较高尾部的性能较差(即高回报期),其GOF指标可能仍然比在较高尾部具有良好性能的模型(后者高)实际上是我们感兴趣的部分)。
      • You might want to try a different fitting method to double-check your results. extRemes provides several useful options, which can be controlled by the argument method.

      我希望这有帮助。
      问候,
      马蒂亚斯

      马蒂亚斯 在4年前 回复


  • 嗨,Matthias,您在这里的帖子真的很棒。很好!我现在也正在研究一个类似的主题,基本上是关于R中的极端事件,我应该使用Block Maxima方法和Peaks Over Thershold方法,以及不同的分布来计算给定数据的风险价值(VaR)和预期不足(ES) 。您能帮我提供密码吗?

    最好的祝福,
    玛蒂娜(Martina)

    玛蒂娜(Martina) 3年前 回复


  • 请问我能得到关于如何在上面绘制的阈值图像上绘制峰的代码吗

    伊曼纽尔·科兰滕(Emmanuel Koranteng) 2年前 回复


  • 我如何从NetCDF(TRMM数据集)文件中读取降水值,以分析每日极端降水值。

    帕万 2年前 回复


  • 您好,非常感谢您对我们的博客很有帮助。

    如何获得mrplot的矩阵?带有阈值,超出次数,平均超出量等列的字段将以无形方式返回。

    先感谢您,

    帕斯帕蒂玛利亚

    玛丽亚·帕斯帕蒂(Maria Paspati) 2年前 回复


  • 我正在研究R中的极端降雨,但是我无法使用R来估计伽玛pareto和伽玛广义pareto分布的参数。能否请我帮助我使用R studio中的代码?

    吉恩 1年前 回复


发表回复

*