原文链接：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 shape
mom 1.9538076495 0.2423393
mle 2.0345084386 0.2053905
pwmu 2.0517348996 0.2043644
pwmb 2.0624399910 0.2002131
pickands 2.3693985422 -0.0708419
med 2.2194363549 0.1537701
mdpd 2.0732577511 0.1809110
mple 2.0499646631 0.1960452
ad2r 0.0005539296 27.5964097

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


> fit(x, thresh = 1, shape = 0, est = "mle")
Estimator: MLE
Deviance: 322.686
AIC: 324.686
Varying Threshold: FALSEThreshold Call: 1
Number Above: 100
Proportion Above: 1
Estimates
scale
1.847
Standard Error Type: observed
Standard Errors
scale
0.1847
Asymptotic Variance Covariance
scale
scale 0.03410
Optimization Information
Convergence: successful
Function Evaluations: 7
Gradient Evaluations: 1
> fitgpd(x, thresh = 1, scale = 2, est = "mle")
Estimator: MLE
Deviance: 323.3049
AIC: 325.3049
Varying Threshold: FALSE
Threshold Call: 1
Number Above: 100
Proportion Above: 1
Estimates
shape
0.0003398
Standard Error Type: observed
Standard Errors
shape
0.06685
Asymptotic Variance Covariance
shape
shape 0.004469
Optimization Information
Convergence: successful
Function Evaluations: 5
Gradient Evaluations: 1
If 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: MLE
Deviance: -176.1669
AIC: -172.1669
Varying Threshold: TRUE
Threshold Call: 1:2
Number Above: 500
Proportion Above: 1
Estimates
scale shape
0.3261 -0.0556
Standard Error Type: observed
Standard Errors
scale shape
0.02098 0.04632
Asymptotic Variance Covariance
scale shape
scale 0.0004401 -0.0007338
shape -0.0007338 0.0021451
Optimization Information
Convergence: successful
Function Evaluations: 62
Gradient Evaluations: 11

scale1	shape1	scale2	shape2	α
6.784e-02	5.303e-02	2.993e-02	3.718e-02	2.001e-06

Asymptotic Variance Covariance
scale1 shape1 scale2 shape2 alpha
scale1 4.602e-03 -2.285e-03 1.520e-06 -1.145e-06 -3.074e-11
shape1 -2.285e-03 2.812e-03 -1.337e-07 4.294e-07 -1.843e-11
scale2 1.520e-06 -1.337e-07 8.956e-04 -9.319e-04 8.209e-12
shape2 -1.145e-06 4.294e-07 -9.319e-04 1.382e-03 5.203e-12
alpha -3.074e-11 -1.843e-11 8.209e-12 5.203e-12 4.003e-12
Optimization Information
Convergence: successful
Function Evaluations: 150
Gradient Evaluations: 21

双变量情况

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

vgpd(cbind(x, y), c(0, 2), model = "log")
> Mlog
Estimator: MLE
Dependence Model and Strenght:
Model : Logistic
lim_u Pr[ X_1 > u | X_2 > u] = NA
Deviance: 1313.260
AIC: 1323.260
Marginal Threshold: 0 2
Marginal Number Above: 500 500
Marginal Proportion Above: 1 1
Joint Number Above: 500
Joint Proportion Above: 1
Number of events such as {Y1 > u1} U {Y2 > u2}: 500
Estimates
scale1 shape1 scale2 shape2 alpha
0.9814 0.2357 0.5294 -0.2835 0.9993
Standard 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 : Logistic
lim_u Pr[ X_1 > u | X_2 > u] = NA
Deviance: 1312.842
AIC: 1320.842
Marginal Threshold: 0 2
Marginal Number Above: 500 500
Marginal Proportion Above: 1 1
Joint Number Above: 500
Joint Proportion Above: 1
Number of events such as {Y1 > u1} U {Y2 > u2}: 500
Estimates
scale1 shape1 scale2 shape2
0.9972 0.2453 0.5252 -0.2857
Standard Errors
scale1 shape1 scale2 shape2
0.07058 0.05595 0.02903 0.03490
Asymptotic Variance Covariance

scale1 shape1 scale2 shape2
scale1 4.982e-03 -2.512e-03 -3.004e-13 3.544e-13
shape1 -2.512e-03 3.130e-03 1.961e-13 -2.239e-13
scale2 -3.004e-13 1.961e-13 8.427e-04 -8.542e-04
shape2 3.544e-13 -2.239e-13 -8.542e-04 1.218e-03
Optimization Information
Convergence: successful
Function Evaluations: 35
Gradient Evaluations: 9

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

> pickdep(Mlog)

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

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

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

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

mcgpd(mc, 2, "log")
Estimator: MLE
Dependence Model and Strenght:
Model : Logistic
lim_u Pr[ X_1 > u | X_2 > u] = NA
Deviance: 1334.847
AIC: 1340.847
Threshold Call:
Number Above: 998
Proportion Above: 1
Estimates
scale shape alpha
0.8723 0.1400 0.5227
Standard Errors
scale shape alpha
0.08291 0.04730 0.02345
Asymptotic Variance Covariance
scale shape alpha
scale 0.0068737 -0.0016808 -0.0009066
shape -0.0016808 0.0022369 -0.0004412
alpha -0.0009066 -0.0004412 0.0005501
Optimization Information
Convergence: successful
Function Evaluations: 70
Gradient Evaluations: 15

置信区间

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


conf.inf.scale conf.sup.scale
2.110881 2.939282conf.inf.scale conf.sup.scale
1.633362 3.126838conf.inf.scale conf.sup.scale
2.138851 3.074945conf.inf.scale conf.sup.scale
2.149123 3.089090

对于形状参数置信区间：

conf.inf conf.sup
2.141919 NA

conf.inf conf.sup
0.05757576 0.26363636

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

conf.inf conf.sup
2.141919 NA

conf.inf conf.sup
0.05757576 0.26363636

 conf.inf conf.sup
8.64712 12.22558

conf.inf conf.sup
8.944444 12.833333


conf.inf conf.sup
8.64712 12.22558

conf.inf conf.sup
8.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 prob
1 1.8 50 0.9888889
npy retper prob
1 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)
> np
time obs
Min. :1970 Min. : 0.022
1st Qu.:1981 1st Qu.: 0.236
Median :1991 Median : 0.542
Mean :1989 Mean : 1.024
3rd Qu.:1997 3rd Qu.: 1.230
Max. :2004 Max. :44.200
NA's : 1.000

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

> rp2p
npy retper prob
1 1.707897 100 0.9941448
>
scale
36.44331

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

conf.inf conf.sup
25.56533 90.76633

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

> maxEvent
[1] 44.2shape
0.9974115
> prob
npy retper prob
1 1.707897 226.1982 0.9974115

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

conf.inf conf.sup
25.56533 90.76633

最受欢迎的见解

1.R语言基于ARMA-GARCH-VaR模型拟合和预测实证研究

2.R语言时变参数VAR随机模型

3.R语言时变参数VAR随机模型

4.R语言基于ARMA-GARCH过程的VAR拟合和预测

5.GARCH（1,1），MA以及历史模拟法的VaR比较

6.R语言时变参数VAR随机模型

7.R语言实现向量自动回归VAR模型

8.R语言随机搜索变量选择SSVS估计贝叶斯向量自回归（BVAR）模型

9.R语言VAR模型的不同类型的脉冲响应分析

拓端tecdat|R语言有极值（EVT）依赖结构的马尔可夫链(MC)对洪水极值分析相关推荐

拓端tecdat|R语言逻辑回归（Logistic回归）模型分类预测病人冠心病风险
最近我们被客户要求撰写关于冠心病风险的研究报告,包括一些图形和统计输出. 相关视频:R语言逻辑回归(Logistic回归)模型分类预测病人冠心病风险逻辑回归Logistic模型原理和R语言分类预测冠 ...
拓端tecdat|R语言用LOESS(局部加权回归)季节趋势分解（STL）进行时间序列异常检测
最近我们被客户要求撰写关于LOESS(局部加权回归)的研究报告,包括一些图形和统计输出. 这篇文章描述了一种对涉及季节性和趋势成分的时间序列的中点进行建模的方法.我们将对一种叫做STL的算法进行研究, ...
拓端tecdat|R语言向量误差修正模型 (VECMs)分析长期利率和通胀率影响关系
最近我们被客户要求撰写关于向量误差修正模型的研究报告,包括一些图形和统计输出. 向量自回归模型估计的先决条件之一是被分析的时间序列是平稳的.但是,经济理论认为,经济变量之间在水平上存在着均衡关系,可以 ...
拓端tecdat|R语言线性回归和时间序列分析北京房价影响因素可视化案例
最近我们被客户要求撰写关于北京房价影响因素的研究报告,包括一些图形和统计输出. 目的房价有关的数据可能反映了中国近年来的变化: 人们得到更多的资源(薪水),期望有更好的房子人口众多独生子女政策: ...
拓端tecdat荣获掘金社区入驻新人奖
2021年7月,由掘金发起了"入驻成长礼"颁奖活动.本次活动邀请到知名开发者.服务机构代表等业界人士. 据了解,掘金社区"新入驻创作者礼"主要对已经积累了一定历 ...
拓端tecdat荣获2022年度51CTO博主之星
相信技术,传递价值,这是51CTO每一个技术创作者的动力与信念,2022 年度,拓端tecdat 作为新锐的数据分析咨询公司,在51CTO平台上,不断的输出优质的技术文章,分享前沿创新技术,输出最佳生 ...
R语言主成分回归（PCR）、多元线性回归特征降维分析光谱数据和汽车油耗、性能数据...
原文链接:http://tecdat.cn/?p=24152 什么是PCR?(PCR = PCA + MLR)(点击文末"阅读原文"获取完整代码数据). • PCR是处理许多 x ...
R语言稀疏主成分分析、因子分析、KMO检验和Bartlett球度检验分析上市公司财务指标数据...
全文链接:http://tecdat.cn/?p=31080 R中的主成分分析(PCA)和因子分析是统计分析技术,也称为多元分析技术(点击文末"阅读原文"获取完整代码数据). 当可 ...
R语言用户自定义函数的语法结构、编写自定义统计值计算函数（使用ifelse结构计算均值和标准差等）、编写自定义日期格式化（format）函数（switch函数使用不同分枝格式化日期数据）、应用自定函数
R语言用户自定义函数的语法结构.编写自定义统计值计算函数(使用ifelse结构计算均值和标准差等).编写自定义日期格式化(format)函数(switch函数使用不同分枝格式化日期数据).应用自定函数 ...
R语言data.table进行滚动数据连接，滚动连接通常用于分析涉及时间的数据（例如商业销售活动和对应的广告投放的安排之之间的关系）实战：实战和动画说明滚动数据连接的形式及方法
R语言data.table进行滚动数据连接,滚动连接通常用于分析涉及时间的数据(例如商业销售活动和对应的广告投放的安排之之间的关系)实战:实战和动画说明滚动数据连接的形式及方法目录

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

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

原文出处：拓端数据部落公众号

EVT的介绍

单变量情况

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

阈值选择图

L-矩图

色散指数图

拟合GPD

单变量情况

双变量情况

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

置信区间

模型检查

聚类技术

其他

返回周期

无偏样本L矩

河流阈值分析

时间序列的移动平均窗口

拓端tecdat|R语言有极值（EVT）依赖结构的马尔可夫链(MC)对洪水极值分析相关推荐

最新文章

热门文章

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

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

原文出处：拓端数据部落公众号

EVT的介绍

单变量情况

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

阈值选择图

L-矩图

色散指数图

拟合GPD

单变量情况

双变量情况

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

置信区间

模型检查

聚类技术

其他

返回周期

无偏样本L矩

河流阈值分析

时间序列的移动平均窗口

拓端tecdat|R语言有极值（EVT）依赖结构的马尔可夫链(MC)对洪水极值分析相关推荐

最新文章

热门文章

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