R语言有极值(EVT)依赖结构的马尔可夫链(MC)对洪水极值分析

原文链接:http://tecdat.cn/?p=17375

为了帮助客户使用POT模型,本指南包含有关使用此模型的实用示例。本文快速介绍了极值理论(EVT)、一些基本示例,最后则通过案例对河流的极值进行了具体的统计分析。

EVT的介绍

单变量情况

假设存在归一化常数an> 0和bn使得:

根据极值类型定理(Fisher和Tippett,1928年),G必须是Fr'echet,Gumbel或负Weibull分布。Jenkinson(1955)指出,这三个分布可以合并为一个参数族:广义极值(GEV)分布。GEV具有以下定义的分布函数:

根据这一结果,Pickands(1975)指出,当阈值接近目标变量的端点µend时,阈值阈值的标准化超额的极限分布是广义Pareto分布(GPD)。也就是说,如果X是一个随机变量,则:

基本用法
随机数和分布函数

首先,让我们从基本的东西开始。将R用于随机数生成和分布函数。


> rgpd(5, loc = 1, scale = 2, shape = -0.2)[1] 1.523393 2.946398 2.517602 1.199393 2.541937> rgpd(6, c(1, -5), 2, -0.2)[1] 1.3336965 -4.6504749 3.1366697 -0.9330325 3.5152161 -4.4851408> rgpd(6, 0, c(2, 3), 0)[1] 3.1139689 6.5900384 0.1886106 0.9797699 3.2638614 5.4755026> pgpd(c(9, 15, 20), 1, 2, 0.25)[1] 0.9375000 0.9825149 0.9922927> qgpd(c(0.25, 0.5, 0.75), 1, 2, 0)[1] 1.575364 2.386294 3.772589> dgpd(c(9, 15, 20), 1, 2, 0.25)[1] 0.015625000 0.003179117 0.001141829

使用选项lower.tail = TRUE或lower.tail = FALSE分别计算不超过或超过概率;
指定分位数是否超过概率分别带有选项lower.tail = TRUE或lower.tail = FALSE;
指定是分别使用选项log = FALSE还是log = TRUE计算密度或对数密度。

阈值选择图

此外,可以使用Fisher信息来计算置信区间。


> x <- runif(10000)> par(mfrow = c(1, 2))

结果如图所示。我们可以清楚地看到,将阈值设为0.98是合理的选择。

可以将置信区间添加到该图,因为经验均值可以被认为是正态分布的(中心极限定理)。但是,对于高阈值,正态性不再成立,此外,通过构造,该图始终会收敛到点(xmax; 0)。
这是另一个综合示例。


> x <- rnorm(10000)plot(x, u.range = c(1, quantile(x, probs = 0.995)), col =

L-矩图

L-矩是概率分布和数据样本的摘要统计量。它们类似于普通矩{它们提供位置,离散度,偏度,峰度以及概率分布或数据样本形状的其他方面的度量值{但是是从有序数据值的线性组合中计算出来的(因此有前缀L)。

这是一个简单的例子。

> x <- c(1 - abs(rnorm(200, 0, 0.2)), rgpd(100, 1, 2, 0.25))

我们发现该图形在真实数据上的性能通常很差。

色散指数图

在处理时间序列时,色散指数图特别有用。EVT指出,超出阈值的超出部分可以通过GPD近似。但是,EVT必须通过泊松过程来表示这些超额部分的发生。

对于下一个示例,我们使用POT包中包含的数据集。此外,由于洪水数据是一个时间序列,因此具有很强的自相关性,因此我们必须“提取”极端事件,同时保持事件之间的独立性。


5, clust.max = TRUE)> diplot(events, u.range = c(2, 20))

色散指数图如图所示。从该图可以看出,大约5的阈值是合理的。

拟合GPD

单变量情况

可以根据多个估算器拟合GPD。
MLE是一种特殊情况,因为它是唯一允许变化阈值的情况。
我们在此给出一些教学示例。




scale shapemom 1.9538076495 0.2423393mle 2.0345084386 0.2053905pwmu 2.0517348996 0.2043644pwmb 2.0624399910 0.2002131pickands 2.3693985422 -0.0708419med 2.2194363549 0.1537701mdpd 2.0732577511 0.1809110mple 2.0499646631 0.1960452ad2r 0.0005539296 27.5964097

MLE,MPLE和MGF估计允许比例或形状参数。例如,如果我们要拟合指数分布:





> fit(x, thresh = 1, shape = 0, est = "mle")Estimator: MLEDeviance: 322.686AIC: 324.686Varying Threshold: FALSE



Threshold Call: 1Number Above: 100Proportion Above: 1Estimatesscale1.847Standard Error Type: observedStandard Errorsscale0.1847Asymptotic Variance Covariancescalescale 0.03410Optimization InformationConvergence: successfulFunction Evaluations: 7Gradient Evaluations: 1> fitgpd(x, thresh = 1, scale = 2, est = "mle")Estimator: MLEDeviance: 323.3049AIC: 325.3049Varying Threshold: FALSEThreshold Call: 1Number Above: 100Proportion Above: 1Estimatesshape0.0003398Standard Error Type: observedStandard Errorsshape0.06685Asymptotic Variance Covarianceshapeshape 0.004469Optimization InformationConvergence: successfulFunction Evaluations: 5Gradient Evaluations: 1If now, we want to fit a GPD with a varying threshold, just do:> x <- rgpd(500, 1:2, 0.3, 0.01)> fitgpd(x, 1:2, est = "mle")Estimator: MLEDeviance: -176.1669AIC: -172.1669Varying Threshold: TRUEThreshold Call: 1:2Number Above: 500Proportion Above: 1Estimatesscale shape0.3261 -0.0556Standard Error Type: observedStandard Errorsscale shape0.02098 0.04632Asymptotic Variance Covariancescale shapescale 0.0004401 -0.0007338shape -0.0007338 0.0021451Optimization InformationConvergence: successfulFunction Evaluations: 62Gradient Evaluations: 11
scale1 shape1 scale2 shape2 α
6.784e-02 5.303e-02 2.993e-02 3.718e-02 2.001e-06

Asymptotic Variance Covariancescale1 shape1 scale2 shape2 alphascale1 4.602e-03 -2.285e-03 1.520e-06 -1.145e-06 -3.074e-11shape1 -2.285e-03 2.812e-03 -1.337e-07 4.294e-07 -1.843e-11scale2 1.520e-06 -1.337e-07 8.956e-04 -9.319e-04 8.209e-12shape2 -1.145e-06 4.294e-07 -9.319e-04 1.382e-03 5.203e-12alpha -3.074e-11 -1.843e-11 8.209e-12 5.203e-12 4.003e-12Optimization InformationConvergence: successfulFunction Evaluations: 150Gradient Evaluations: 21

双变量情况

拟合双变量POT。所有这些模型均使用最大似然估计量进行拟合。


vgpd(cbind(x, y), c(0, 2), model = "log")> MlogEstimator: MLEDependence Model and Strenght:Model : Logisticlim_u Pr[ X_1 > u | X_2 > u] = NADeviance: 1313.260AIC: 1323.260Marginal Threshold: 0 2Marginal Number Above: 500 500Marginal Proportion Above: 1 1Joint Number Above: 500Joint Proportion Above: 1Number of events such as {Y1 > u1} U {Y2 > u2}: 500Estimatesscale1 shape1 scale2 shape2 alpha0.9814 0.2357 0.5294 -0.2835 0.9993Standard Errors

在摘要中,我们可以看到lim_u Pr [X_1> u | X_2> u] = 0.02。这是Coles等人的χ统计量。(1999)。对于参数模型,我们有:

对于自变量,χ= 0,而对于完全依存关系,χ=1。在我们的应用中,值0.02表示变量是独立的{这是显而易见的。从这个角度来看,可以固定一些参数。


vgpd(cbind(x, y), c(0, 2), model = "log"Dependence Model and Strenght:Model : Logisticlim_u Pr[ X_1 > u | X_2 > u] = NADeviance: 1312.842AIC: 1320.842Marginal Threshold: 0 2Marginal Number Above: 500 500Marginal Proportion Above: 1 1Joint Number Above: 500Joint Proportion Above: 1Number of events such as {Y1 > u1} U {Y2 > u2}: 500Estimatesscale1 shape1 scale2 shape20.9972 0.2453 0.5252 -0.2857Standard Errorsscale1 shape1 scale2 shape20.07058 0.05595 0.02903 0.03490Asymptotic Variance Covariance


scale1 shape1 scale2 shape2scale1 4.982e-03 -2.512e-03 -3.004e-13 3.544e-13shape1 -2.512e-03 3.130e-03 1.961e-13 -2.239e-13scale2 -3.004e-13 1.961e-13 8.427e-04 -8.542e-04shape2 3.544e-13 -2.239e-13 -8.542e-04 1.218e-03Optimization InformationConvergence: successfulFunction Evaluations: 35Gradient Evaluations: 9

注意,由于所有双变量极值分布都渐近相关,因此Coles等人的χ统计量。(1999)始终等于1。
检测依赖强度的另一种方法是绘制Pickands的依赖函数:

> pickdep(Mlog)

水平线对应于独立性,而其他水平线对应于完全相关性。请注意,通过构造,混合模型和非对称混合模型无法对完美的依赖度变量建模。

使用马尔可夫链对依赖关系结构进行建模

超越的马尔可夫链进行超过阈值的峰分析的经典方法是使GPD拟合最大值。但是,由于仅考虑群集最大值,因此存在数据浪费。主要思想是使用马尔可夫链对依赖关系结构进行建模,而联合分布显然是多元极值分布。这个想法是史密斯等人首先提出的。(1997)。在本节的其余部分,我们将只关注一阶马尔可夫链。因此,所有超出的可能性为:

对于我们的应用程序,我们模拟具有极值依赖结构的一阶马尔可夫链。


mcgpd(mc, 2, "log")Estimator: MLEDependence Model and Strenght:Model : Logisticlim_u Pr[ X_1 > u | X_2 > u] = NADeviance: 1334.847AIC: 1340.847Threshold Call:Number Above: 998Proportion Above: 1Estimatesscale shape alpha0.8723 0.1400 0.5227Standard Errorsscale shape alpha0.08291 0.04730 0.02345Asymptotic Variance Covariancescale shape alphascale 0.0068737 -0.0016808 -0.0009066shape -0.0016808 0.0022369 -0.0004412alpha -0.0009066 -0.0004412 0.0005501Optimization InformationConvergence: successfulFunction Evaluations: 70Gradient Evaluations: 15

置信区间

拟合统计模型后,通常会给出置信区间。当前,只有mle,pwmu,pwmb,矩估计量可以计算置信区间。




conf.inf.scale conf.sup.scale2.110881 2.939282

conf.inf.scale conf.sup.scale1.633362 3.126838

conf.inf.scale conf.sup.scale2.138851 3.074945

conf.inf.scale conf.sup.scale2.149123 3.089090

对于形状参数置信区间:



conf.inf conf.sup2.141919 NA
conf.inf conf.sup0.05757576 0.26363636

分位数的置信区间也可用。



conf.inf conf.sup2.141919 NA conf.inf conf.sup0.05757576 0.26363636


conf.inf conf.sup8.64712 12.22558 conf.inf conf.sup8.944444 12.833333




conf.inf conf.sup8.64712 12.22558 conf.inf conf.sup8.944444 12.833333

轮廓置信区间函数既返回置信区间,又绘制轮廓对数似然函数。

模型检查

要检查拟合的模型,用户必须调用函数图。




> plot(fitted, npy = 1)

图显示了执行获得的图形窗口。

聚类技术

在处理时间序列时,超过阈值的峰值可能会出现问题。确实,由于时间序列通常是高度自相关的,因此选择高于阈值可能会导致相关事件。
该函数试图在满足独立性标准的同时识别超过阈值的峰。为此,该函数至少需要两个参数:阈值u和独立性的时间条件。群集标识如下:
1.第一次超出会启动第一个集群;
2.在阈值以下的第一个观察结果将“终止集群”,除非tim.cond不成立;
3.下一个超过tim.cond的集群将启动新集群;
4.根据需要进行迭代。
一项初步研究表明,如果两个洪水事件不在8天之内,则可以认为两个洪水事件是独立的,请注意,定义tim.cond的单位必须与所分析的数据相同。
返回一个包含已识别集群的列表。

clu(dat, u = 2, tim.cond = 8/365)

其他

返回周期

将返回周期转换为非超出概率,反之亦然。它既需要返回期,也需要非超越概率。此外,必须指定每年平均事件数。


npy retper prob1 1.8 50 0.9888889npy retper prob1 2.2 1.136364 0.6

无偏样本L矩

函数samlmu计算无偏样本L矩。

l_1 l_2 t_3 t_4 t_5
0.455381591 0.170423740 0.043928262 -0.005645249 -0.009310069
3.7.3

河流阈值分析

在本节中,我们提供了对河流阈值的全面和详细的分析。

时间序列的移动平均窗口

从初始时间序列ts计算“平均”时间序列。这是通过在初始时间序列上使用长度为d的移动平均窗口来实现的。




> plot(dat, type = "l", col = "blue")> lines(tsd, col = "green")

执行过程如图所示。图描绘了瞬时洪水数量-蓝线。

由于这是一个时间序列,因此我们必须选择一个阈值以上的独立事件。首先,我们固定一个相对较低的阈值以“提取”更多事件。因此,其中一些不是极端事件而是常规事件。这对于为GPD的渐近逼近选择一个合理的阈值是必要的。


> par(mfrow = c(2, 2))> plot(even[, "obs"])> abline(v = 6, col = "green"

根据图,阈值6m3 = s应该是合理的。平均剩余寿命图-左上方面板-表示大约10m3 = s的阈值应足够。但是,所选阈值必须足够低,以使其上方具有足够的事件以减小方差,而所选阈值则不能过低,因为它会增加偏差3。
因此,我们现在可以\重新提取阈值6m3 = s以上的事件,以获得对象event1。注意,可以使用极值指数的估计。



> events1 <-365, clust.max = TRUE)> nptime obsMin. :1970 Min. : 0.0221st Qu.:1981 1st Qu.: 0.236Median :1991 Median : 0.542Mean :1989 Mean : 1.0243rd Qu.:1997 3rd Qu.: 1.230Max. :2004 Max. :44.200NA's : 1.000

结果给出了估计器的名称(阈值,阈值以上的观测值的数量和比例,参数估计,标准误差估计和类型,渐近方差-协方差矩阵和收敛性诊断。
图显示了拟合模型的图形诊断。可以看出,拟合模型“ mle”似乎是合适的。假设我们想知道与100年返回期相关的返回水平。


> rp2pnpy retper prob1 1.707897 100 0.9941448>scale36.44331

考虑到不确定性,图描绘了与100年返回期相关的分位数的分布置信区间。


conf.inf conf.sup25.56533 90.76633

有时有必要知道指定事件的估计返回周期。让我们对较大事件进行处理。


> maxEvent[1] 44.2

shape0.9974115> probnpy retper prob1 1.707897 226.1982 0.9974115

因此,2000年6月发生的最大事件的重现期大约为240年。



conf.inf conf.sup25.56533 90.76633
(0)

相关推荐