关于数据挖掘: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 的阈值是正当的。

点击题目查阅往期内容

极值实践 EVT、POT 超阈值、GARCH 模型剖析股票指数 VaR、条件 CVaR：多元化投资组合预测危险测度剖析

左右滑动查看更多

01

02

03

04

拟合 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: FALSE


Threshold 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.939282

conf.inf.scale conf.sup.scale
1.633362 3.126838

conf.inf.scale conf.sup.scale
2.138851 3.074945

conf.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.2

shape
0.9974115
> prob
npy retper prob
1 1.707897 226.1982 0.9974115

因而，2000 年 6 月产生的最大事件的重现期大概为 240 年。

conf.inf conf.sup
25.56533 90.76633

点击文末 “浏览原文”

获取全文残缺材料。

本文选自《 R 语言有极值（EVT）依赖构造的马尔可夫链 (MC) 对洪水极值剖析》。

点击题目查阅往期内容

R 语言极值实践：希尔 HILL 统计量尾部指数参数估计可视化
极值实践 EVT、POT 超阈值、GARCH 模型剖析股票指数 VaR、条件 CVaR：多元化投资组合预测危险测度剖析
R 语言 POT 超阈值模型和极值实践 EVT 剖析
R 语言极值推断：狭义帕累托散布 GPD 应用极大似然预计、轮廓似然预计、Delta 法
R 语言极值实践 EVT：基于 GPD 模型的火灾损失散布剖析
R 语言有极值（EVT）依赖构造的马尔可夫链 (MC) 对洪水极值剖析
R 语言 POT 超阈值模型和极值实践 EVT 剖析
R 语言混合正态分布极大似然预计和 EM 算法
R 语言多项式线性模型：最大似然预计二次曲线
R 语言 Wald 测验 vs 似然比测验
R 语言 GARCH-DCC 模型和 DCC（MVT）建模预计
R 语言非参数办法：应用核回归平滑预计和 K -NN(K 近邻算法)分类预测心脏病数据
matlab 实现 MCMC 的马尔可夫转换 ARMA – GARCH 模型预计
R 语言基于 Bootstrap 的线性回归预测置信区间预计办法
R 语言随机搜寻变量抉择 SSVS 预计贝叶斯向量自回归（BVAR）模型
Matlab 马尔可夫链蒙特卡罗法（MCMC）预计随机稳定率（SV，Stochastic Volatility）模型
Matlab 马尔可夫区制转换动静回归模型预计 GDP 增长率 R 语言极值推断：狭义帕累托散布 GPD 应用极大似然预计、轮廓似然预计、Delta 法