关于数据挖掘:R语言有极值EVT依赖结构的马尔可夫链MC对洪水极值分析附代码数据

27次阅读

共计 9372 个字符,预计需要花费 24 分钟才能阅读完成。

浏览全文: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 法

正文完
 0