# 7章
# 7.2.1項
library(evd)
dgpd(3, loc=1, scale=2, shape=0.5) # GPの密度関数の値
## [1] 0.1481481
1/2*(1+0.5*(3-1)/2)^(-1/0.5-1) # 上の値をg(y)に数値を代入して直接計算する
## [1] 0.1481481
pgpd(5,loc=2, scale=1.5, shape=0.2) # GPの分布関数の値
## [1] 0.8140656
curve(pgpd(x,loc=2, scale=1.5, shape=0.2),xlim=c(1,10)) # 分布関数のグラフ

qgpd(0.8, loc=1, scale=2, shape=0.7) # G(x)=0.8となる点(分位数)
## [1] 6.957627
pgpd(6.957627, loc=1, scale=2, shape=0.7) # 上の分位数を確認
## [1] 0.8
rgpd(5, loc=1, scale=2, shape=0.3) # GPに従う乱数を5個
## [1] 3.005424 1.595850 3.556109 1.370340 10.660415
# 7.2.2項
library(evd)
set.seed(100)
tmp <- rgpd(10000) # GP(0,1,0)に従う乱数を10000個
hist(tmp) # 乱数の分布の確認

mrlplot(tmp) # 標本平均超過プロット

oldpar <- par()
par(mfrow=c(1,2))
tcplot(tmp, tlim=c(0,7)) # GPモデルによる母数推定プロット

par(mfcol=c(1,3))
tcplot(tmp, model="pp", tlim=c(0,7)) # PPモデルによる母数推定プロット

sum(tmp>7) # 7を超える乱数の数
## [1] 11
(tmp.g1 <- fpot(tmp, 1)) # GPモデル u=1
##
## Call: fpot(x = tmp, threshold = 1)
## Deviance: 7441.09
##
## Threshold: 1
## Number Above: 3684
## Proportion Above: 0.3684
##
## Estimates
## scale shape
## 0.97816 0.03201
##
## Standard Errors
## scale shape
## 0.02309 0.01692
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 35
## Gradient Evaluations: 5
(tmp.g7 <- fpot(tmp, 7)) # GPモデル u=7
##
## Call: fpot(x = tmp, threshold = 7)
## Deviance: 27.08549
##
## Threshold: 7
## Number Above: 11
## Proportion Above: 0.0011
##
## Estimates
## scale shape
## 1.8014 -0.3574
##
## Standard Errors
## scale shape
## 0.7273 0.2893
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 24
## Gradient Evaluations: 9
(tmp.p1 <- fpot(tmp, model="pp", 1)) # PPモデル u=1
##
## Call: fpot(x = tmp, threshold = 1, model = "pp")
## Deviance: -45694.36
##
## Threshold: 1
## Number Above: 3684
## Proportion Above: 0.3684
##
## Estimates
## loc scale shape
## 9.76431 1.14605 0.01756
##
## Standard Errors
## loc scale shape
## 0.4060 0.1136 0.0136
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 344
## Gradient Evaluations: 59
# (tmp.p7 <- fpot(tmp, model="pp", 7)) # PPモデル u=7
# (tmp.p7 <- fpot(tmp, model="pp", 7, method="Nelder-Mead")) # PPモデル u=7
tmp.g1.pro <- profile(tmp.g1) # GP u=1のプロファイル
## [1] "profiling scale"
## [1] "profiling shape"
plot(tmp.g1.pro) # 上のプロファイルのプロット
tmp.g7.pro <- profile(tmp.g7) # GP u=7のプロファイル 計算不可
## [1] "profiling scale"
## [1] "profiling shape"
confint(tmp.g1) # GP u=1のワルド区間推定
## 2.5 % 97.5 %
## scale 0.932891325 1.0234209
## shape -0.001156472 0.0651684
confint(tmp.g7) # GP u=7のワルド区間推定
## 2.5 % 97.5 %
## scale 0.3759408 3.2267971
## shape -0.9244318 0.2096428
confint(tmp.p1) # PP u=1のワルド区間推定
## 2.5 % 97.5 %
## loc 8.968529533 10.56009489
## scale 0.923420945 1.36867622
## shape -0.009097813 0.04422255
confint(tmp.g1.pro) # GP u=1のプロファイル区間推定
## lower upper
## scale 0.9337310958 1.02418655
## shape 0.0001977867 0.06643153
par(mfrow=c(2,2))

plot(tmp.g1) # GPモデル u=1 のEDA

plot(tmp.g7) # GPモデル u=7 のEDA

plot(tmp.p1) # PPモデル u=1 のEDA

par(oldpar)
# 7.2.4項
library(evd)
library(ismev)
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
set.seed(100)
tmp <- rgpd(10000) # GP(0,1,0)に従う乱数を10000個
hist(tmp) # 乱数の分布の確認

mrl.plot(tmp) # 標本平均超過プロット

gpd.fitrange(tmp, umin=1, umax=8) # GPモデルによる母数推定プロット

pp.fitrange(tmp, umin=1, umax=8) # PPモデルによる母数推定プロット

tmp.g1.is <- gpd.fit(tmp,1) # GPモデル u=1
## $threshold
## [1] 1
##
## $nexc
## [1] 3684
##
## $conv
## [1] 0
##
## $nllh
## [1] 3720.545
##
## $mle
## [1] 0.9782663 0.0319127
##
## $rate
## [1] 0.3684
##
## $se
## [1] 0.02309609 0.01691508
tmp.g7.is <- gpd.fit(tmp,7) # GPモデル u=7
## $threshold
## [1] 7
##
## $nexc
## [1] 11
##
## $conv
## [1] 0
##
## $nllh
## [1] 13.54275
##
## $mle
## [1] 1.8016879 -0.3575407
##
## $rate
## [1] 0.0011
##
## $se
## [1] 0.7274643 0.2893074
tmp.p1.is <- pp.fit(tmp,1) # PPモデル u=1
## $threshold
## [1] 1
##
## $npy
## [1] 365
##
## $nexc
## [1] 3684
##
## $conv
## [1] 0
##
## $nllh
## [1] -10651.89
##
## $mle
## [1] 6.19195346 1.14491210 0.03213587
##
## $se
## [1] 0.16243564 0.07823430 0.01682988
tmp.p7.is <- pp.fit(tmp,7) # PPモデル u=7
## $threshold
## [1] 7
##
## $npy
## [1] 365
##
## $nexc
## [1] 11
##
## $conv
## [1] 0
##
## $nllh
## [1] 34.58077
##
## $mle
## [1] 5.055955 2.496722 -0.357519
##
## $se
## [1] 1.2789577 1.6466989 0.2891916
gpd.diag(tmp.g1.is) # GPモデル u=1 のEDA

gpd.diag(tmp.g7.is) # GPモデル u=7 のEDA

pp.diag(tmp.p1.is) # PPモデル u=1 のEDA

pp.diag(tmp.p7.is) # PPモデル u=7 のEDA

# 7.3節
library(evd)
set.seed(100)
tmpn <- rnorm(2000) # 正規乱数2000個
hist(tmpn) # 乱数の分布の確認

mrlplot(tmpn) # 標本平均超過プロット

max(tmpn) # 乱数の最大値
## [1] 3.304151
mrlplot(tmpn, tlim=c(-1,3)) # 範囲を指定する

oldpar <- par()
par(mfrow=c(2,1))
tcplot(tmpn, tlim=c(-1,3)) # GPモデルによる区間推定

par(mfrow=c(3,1))
# tcplot(tmpn, model="pp", tlim=c(-1,3)) # PPモデルによる区間推定
tcplot(tmpn, model="pp", tlim=c(-1,2)) # PPモデルによる区間推定 範囲を取り直す

sum(tmpn>1.5) # 1.5を超えるデータサイズ
## [1] 126
(tmpn.pot <- fpot(tmp, 1.5)) # GPモデル u=1.5での推定
##
## Call: fpot(x = tmp, threshold = 1.5)
## Deviance: 4535.206
##
## Threshold: 1.5
## Number Above: 2239
## Proportion Above: 0.2239
##
## Estimates
## scale shape
## 0.9624 0.0511
##
## Standard Errors
## scale shape
## 0.02996 0.02289
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 22
## Gradient Evaluations: 6
tmpn.prof <- profile(tmpn.pot) # そのプロファイル
## [1] "profiling scale"
## [1] "profiling shape"
confint(tmpn.pot) # ワルド信頼区間
## 2.5 % 97.5 %
## scale 0.903666456 1.02112159
## shape 0.006242669 0.09596399
par(mfrow=c(2,2))
plot(tmpn.pot) # GPモデルのEDA

(tmpn.pot.p <- fpot(tmpn, model="pp", 1.5)) # PPモデル u=1.5での推定
##
## Call: fpot(x = tmpn, threshold = 1.5, model = "pp")
## Deviance: -901.635
##
## Threshold: 1.5
## Number Above: 126
## Proportion Above: 0.063
##
## Estimates
## loc scale shape
## 3.1690 0.1624 -0.2804
##
## Standard Errors
## loc scale shape
## 0.1273 0.0447 0.0752
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 228
## Gradient Evaluations: 75
confint(tmpn.pot.p) # ワルド区間推定
## 2.5 % 97.5 %
## loc 2.91956604 3.4184339
## scale 0.07483215 0.2500401
## shape -0.42780199 -0.1330205
plot(tmpn.pot.p) # PPモデルのEDA

set.seed(100)
tmpe <- rexp(10000) # 指数乱数10000個
par(oldpar)
mrlplot(tmpe) # 標本平均超過プロット

mrlplot(tmpe, tlim=c(3,8)) # 区間を取り直す

sum(tmpe>4) # 4を超えるデータの数
## [1] 204
(tmpe.pot <- fpot(tmpe,4)) # GPモデル u=4で推定
##
## Call: fpot(x = tmpe, threshold = 4)
## Deviance: 459.8189
##
## Threshold: 4
## Number Above: 204
## Proportion Above: 0.0204
##
## Estimates
## scale shape
## 1.18710 -0.04451
##
## Standard Errors
## scale shape
## 0.11380 0.06556
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 18
## Gradient Evaluations: 6
confint(tmpe.pot) # そのワルド信頼区間
## 2.5 % 97.5 %
## scale 0.9640492 1.4101521
## shape -0.1730035 0.0839897
par(mfrow=c(2,2))
plot(tmpe.pot) # GPモデルのEDA

# (tmpe.pot.p <- fpot(tmpe,model="pp",4)) # PPモデル u=4で推定
(tmpe.pot.p <- fpot(tmpe,model="pp",3)) # PPモデルでu=3に取り直す
##
## Call: fpot(x = tmpe, threshold = 3, model = "pp")
## Deviance: -3960.987
##
## Threshold: 3
## Number Above: 487
## Proportion Above: 0.0487
##
## Estimates
## loc scale shape
## 9.961181 1.118511 -0.001839
##
## Standard Errors
## loc scale shape
## 0.75358 0.27757 0.04706
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 361
## Gradient Evaluations: 98
confint(tmpe.pot.p) # そのワルド信頼区間
## 2.5 % 97.5 %
## loc 8.48418650 11.43817618
## scale 0.57448771 1.66253460
## shape -0.09407089 0.09039203
plot(tmpe.pot.p) # PPモデルのEDA

sum(tmpe>3) # 3を超えるデータの数
## [1] 487
gpd.fit(tmpe,4) # ismevパッケージのGPモデル
## $threshold
## [1] 4
##
## $nexc
## [1] 204
##
## $conv
## [1] 0
##
## $nllh
## [1] 229.9095
##
## $mle
## [1] 1.18691478 -0.04447095
##
## $rate
## [1] 0.0204
##
## $se
## [1] 0.11377471 0.06555223
pp.fit(tmpe,3) # ismevパッケージのPPモデル
## $threshold
## [1] 3
##
## $npy
## [1] 365
##
## $nexc
## [1] 487
##
## $conv
## [1] 0
##
## $nllh
## [1] -368.3076
##
## $mle
## [1] 6.247467693 1.125539987 -0.001714457
##
## $se
## [1] 0.16963382 0.11342593 0.04803009
par(oldpar)
# 7.4節------------------------------------
library(evd)
library(ismev)
# venice2
str(venice2)
## 'data.frame': 125 obs. of 10 variables:
## $ 1 : num 94 90 97 107 87 86 94 67 105 115 ...
## $ 2 : num 93 84 75 104 72 77 87 66 90 115 ...
## $ 3 : num 90 84 74 85 70 74 80 65 87 114 ...
## $ 4 : num 86 78 72 81 67 70 72 64 80 106 ...
## $ 5 : num 85 75 72 79 66 70 71 63 80 87 ...
## $ 6 : num 82 75 68 72 66 69 71 62 77 86 ...
## $ 7 : num 81 72 68 72 65 69 70 61 77 78 ...
## $ 8 : num 80 72 68 70 64 68 69 60 77 77 ...
## $ 9 : num 79 69 67 70 63 68 66 60 77 76 ...
## $ 10: num 76 69 67 67 62 66 65 59 75 75 ...
head(venice2); tail(venice2) # はじめと終わりのデータを見る
## 1 2 3 4 5 6 7 8 9 10
## 1887 94 93 90 86 85 82 81 80 79 76
## 1888 90 84 84 78 75 75 72 72 69 69
## 1889 97 75 74 72 72 68 68 68 67 67
## 1890 107 104 85 81 79 72 72 70 70 67
## 1891 87 72 70 67 66 66 65 64 63 62
## 1892 86 77 74 70 70 69 69 68 68 66
## 1 2 3 4 5 6 7 8 9 10
## 2006 112 105 103 100 98 98 96 93 93 93
## 2007 109 102 101 100 100 99 96 95 94 93
## 2008 156 118 116 116 115 110 109 109 108 107
## 2009 145 144 133 131 123 120 119 118 117 116
## 2010 136 124 122 122 120 119 117 114 114 114
## 2011 112 96 92 89 89 88 87 87 86 86
ve <- c(as.matrix(venice2)) # 行列にした後にベクトルにする
plot(ve)

vets <- ts(ve, start=c(1887,1), end=c(2011,12), frequency=12) # 時系列にする
plot(vets) # 時系列のプロット

# mrlplot(ve) # 標本平均超過プロット
tcplot(ve, tlim=c(90,140)) # 母数の推定プロット


vern <- na.omit(ve) # NAを取り除く
sum(vern>110); sum(vern>120) # u=110, 120を超えるデータサイズ
## [1] 206
## [1] 82
(ve.gp <- fpot(vern, 110)) # evdでGPモデルによる推定
##
## Call: fpot(x = vern, threshold = 110)
## Deviance: 1417.329
##
## Threshold: 110
## Number Above: 206
## Proportion Above: 0.1665
##
## Estimates
## scale shape
## 11.62882 -0.01337
##
## Standard Errors
## scale shape
## 1.06812 0.05986
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 15
## Gradient Evaluations: 4
oldpar <- par()
par(mfrow=c(2,2))
plot(ve.gp) # EDA

# (ve.pp <- fpot(vern, model="pp", 110)) # evdでPPモデルによる推定
ve.gp2 <- gpd.fit(vern, 110) # ismevでGPモデルによる推定
## $threshold
## [1] 110
##
## $nexc
## [1] 206
##
## $conv
## [1] 0
##
## $nllh
## [1] 708.6644
##
## $mle
## [1] 11.63026229 -0.01360121
##
## $rate
## [1] 0.1665319
##
## $se
## [1] 1.06789211 0.05977572
gpd.diag(ve.gp2) # EDA

ve.pp2 <- pp.fit(vern, 110) # ismevでPPモデルによる推定
## $threshold
## [1] 110
##
## $npy
## [1] 365
##
## $nexc
## [1] 206
##
## $conv
## [1] 0
##
## $nllh
## [1] 445.5966
##
## $mle
## [1] 242.6276336 120.1823238 0.6509444
##
## $se
## [1] 2.000446e-06 2.000446e-06 2.000444e-06
pp.diag(ve.pp2) # EDA

plot(fpot(vern, 110, npp=12)) # 1年12月として推定のEDA

fpot(vern, 110, npp=12, mper=100) # 100年再現レベル
##
## Call: fpot(x = vern, threshold = 110, npp = 12, mper = 100)
## Deviance: 1417.329
##
## Threshold: 110
## Number Above: 206
## Proportion Above: 0.1665
##
## Estimates
## rlevel shape
## 169.5579 -0.0127
##
## Standard Errors
## rlevel shape
## 7.06050 0.06014
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 36
## Gradient Evaluations: 18
par(oldpar)
# 7.5節---------------------------------------------------------
library(evd)
wind <- read.csv("TateyamaKaze.csv") # データ読み込み
# wind
wind2 <- na.omit(wind) # NAの削除
head(wind2)
## X kaze 品質情報 X.1 品質情報.1 均質番号
## 221 May-68 9.8 4 南西 8 1
## 222 Jun-68 11.3 8 南西 8 1
## 223 Jul-68 13.0 8 南西 8 1
## 224 Aug-68 12.8 8 南南西 8 1
## 225 Sep-68 7.5 8 東北東 8 1
## 226 Oct-68 11.8 8 南西 8 1
windts <- ts(wind2[,2], start=c(1968,5), frequency=12) # 時系列に変換
plot(windts) # 時系列のグラフ

max(wind2[,2]) # 最大の風速はいくらか
## [1] 28.4
tail(wind2); length(wind2[,2]) # データの末尾とデータの長さを見る
## X kaze 品質情報 X.1 品質情報.1 均質番号
## 833 May-19 11.7 8 南南東 8 2
## 834 Jun-19 12.7 8 南西 8 2
## 835 Jul-19 10.3 8 南南西 8 2
## 836 Aug-19 11.7 8 南西 8 2
## 837 Sep-19 28.4 8 南 8 2
## 838 Oct-19 20.7 8 南南西 8 2
## [1] 618
wind3 <- wind2[1:608,2] # 風速の列をとり、最後の10か月分を除く
mrlplot(wind3) # 標本平均超過プロット

sum(wind3 >14); sum(wind3 >12) # u=14, 12以上のデータサイズ
## [1] 71
## [1] 248
# tcplot(wind3, tlim=c(10,20)) # 推定プロット
sum(wind3 >13) # 13を超えるデータサイズ
## [1] 133
(tate.gp <- fpot(wind3, 13, npp=12)) # 1年12月としてGPモデルで推定
##
## Call: fpot(x = wind3, threshold = 13, npp = 12)
## Deviance: 423.6642
##
## Threshold: 13
## Number Above: 133
## Proportion Above: 0.2188
##
## Estimates
## scale shape
## 1.83442 -0.01401
##
## Standard Errors
## scale shape
## 0.2515 0.1062
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 16
## Gradient Evaluations: 5
oldpar <- par()
par(mfrow=c(2,2))
plot(tate.gp) # EDA

(tate.pp <- fpot(wind3, 13, model="pp")) # PPモデルで推定
##
## Call: fpot(x = wind3, threshold = 13, model = "pp")
## Deviance: -611.1686
##
## Threshold: 13
## Number Above: 133
## Proportion Above: 0.2188
##
## Estimates
## loc scale shape
## 21.65726 1.70567 -0.01511
##
## Standard Errors
## loc scale shape
## 1.4851 0.7119 0.1053
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 168
## Gradient Evaluations: 43
plot(tate.pp) # EDA

(tate.gp2 <- fpot(wind3, 13, npp=12, mper=100)) # 100年再現レベル
##
## Call: fpot(x = wind3, threshold = 13, npp = 12, mper = 100)
## Deviance: 423.6642
##
## Threshold: 13
## Number Above: 133
## Proportion Above: 0.2188
##
## Estimates
## rlevel shape
## 22.82969 -0.01401
##
## Standard Errors
## rlevel shape
## 2.0014 0.1062
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 17
## Gradient Evaluations: 5
plot(tate.gp2) # EDA

(temp <- tate.gp$estimate) # 繰り返しの面倒を避ける
## scale shape
## 1.83441766 -0.01400603
# qgpd(28.4, scale=temp[1], shape=temp[2]) # 分位数
pgpd(28.4, scale=temp[1], shape=temp[2]) # 分布関数
## shape
## 1
par(oldpar)