# 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)