# 6章
# 6.2.1項
library(evd)
# ?fgev # ブラウザのヘルプが立ち上がる。ページ下からindexに移る
library(ismev)
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
# ?ismev # ブラウザのヘルプが立ち上がる。ページ下からindexに移る
demo(exchange.rate) # デモの一つを実演
##
##
## demo(exchange.rate)
## ---- ~~~~~~~~~~~~~
##
## > data(euroex)
##
## > length(euroex)
## [1] 975
##
## > plot(euroex, type = "l", xlab = "Day Number", ylab = "Exchange Rate")

##
## > euroex.ldr <- log(euroex[-1]/euroex[-975])
##
## > plot(euroex.ldr, type = "l", xlab = "Day Number", ylab = "Log-Daily Return")

##
## > euroex.ldr <- 100*euroex.ldr
##
## > mrl.plot(euroex.ldr)

##
## > sum(euroex.ldr > 0.9)
## [1] 39
##
## > gpd.fitrange(euroex.ldr, -1, 1.4, nint = 100)

##
## > euroex.gpd <- gpd.fit(euroex.ldr, 0.9, npy = 250)
## $threshold
## [1] 0.9
##
## $nexc
## [1] 39
##
## $conv
## [1] 0
##
## $nllh
## [1] -9.420511
##
## $mle
## [1] 0.3534629 -0.2015665
##
## $rate
## [1] 0.04004107
##
## $se
## [1] 0.07277418 0.13338195
##
##
## > gpd.diag(euroex.gpd)

##
## > gpd.profxi(euroex.gpd, -0.5, 0.3)
## If routine fails, try changing plotting interval

##
## > gpd.prof(euroex.gpd, m = 10, npy = 250, 1.75, 3.5)
## If routine fails, try changing plotting interval

##
## > trend <- as.matrix((-500:473)/250)
##
## > euroex.gpd2 <- gpd.fit(euroex.ldr, 0.9, npy = 250, ydat = trend,
## + sigl = 1, siglink = exp)
## $model
## $model[[1]]
## [1] 1
##
## $model[[2]]
## NULL
##
##
## $link
## [1] "c(exp, identity)"
##
## $threshold
## [1] 0.9
##
## $nexc
## [1] 39
##
## $conv
## [1] 0
##
## $nllh
## [1] -12.46977
##
## $mle
## [1] -0.9325977 0.2067568 -0.3842224
##
## $rate
## [1] 0.04004107
##
## $se
## [1] 0.20806649 0.05980706 0.14626536
##
##
## > gpd.diag(euroex.gpd2)

# 6.2.2項
library(evd)
pgev(0.5) # GEV(0,1,0)の分布関数のx=0.5での値
## [1] 0.5452392
exp(-exp(-0.5)) # 上の分布関数を定義通り計算する
## [1] 0.5452392
curve(dnweibull(x, loc=1, scale=2, shape=3), xlim=c(-4,2))

# Weibull(1,2,3)分布の密度関数のグラフ
curve(pnweibull(x, loc=1, scale=2, shape=3), xlim=c(-4,2))

# Weibull(1,2,3)分布の分布関数のグラフ
qgumbel(0.7, scale=2) # Gumbel(0,2)の分布関数F(x)=0.7となるx
## [1] 2.061861
pgumbel(2.061861, scale=2) # 上の分位数を確認
## [1] 0.7
rfrechet(5) # Frechet(0,1,1)分布に従う乱数を5個
## [1] 2.8737969 0.6251997 0.8635976 2.7185670 0.4182139
# 6.2.3項
library(evd)
set.seed(5) # 乱数の種を指定
x <- rgev(100) # GEV(0,1,0)の乱数を100個発生
plot(sort(x)) # xを昇順に並べてグラフで様子を見る

(x.gev.evd <- fgev(x)) # GEVモデルで推定する
##
## Call: fgev(x = x)
## Deviance: 296.3341
##
## Estimates
## loc scale shape
## -0.075216 0.906971 -0.002022
##
## Standard Errors
## loc scale shape
## 0.10210 0.07392 0.07472
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 20
## Gradient Evaluations: 6
fitted(x.gev.evd)
## loc scale shape
## -0.075216206 0.906970565 -0.002021981
plot(x.gev.evd) # 推定結果のEDA




# 以下は1つの画面に4つのプロットを出す方法
oldpar <- par() # 現在の設定を退避する
par(mfrow=c(2,2)) # 画面を2*2型で横を先に入れる
plot(x.gev.evd) # プロットを描画

par(oldpar) # 元の設定に戻す
x.prof <- profile(x.gev.evd)
## [1] "profiling loc"
## [1] "profiling scale"
## [1] "profiling shape"
oldpar <- par()
par(mfrow=c(1,3)) # 3つのグラフを並べる
plot(x.prof, ci=c(0.95, 0.99)) # プロファイル尤度のグラフを描く

par(oldpar)
confint(x.gev.evd, level=0.95) # ワルドの区間推定、0.95のときはlevelを省略可
## 2.5 % 97.5 %
## loc -0.2753347 0.1249023
## scale 0.7620936 1.0518475
## shape -0.1484777 0.1444337
confint(x.prof) # プロファイル区間推定
## lower upper
## loc -0.2714348 0.1306632
## scale 0.7775528 1.0713340
## shape -0.1388744 0.1556046
fgev(x, prob=0.01) # 100年に1回起こる
##
## Call: fgev(x = x, prob = 0.01)
## Deviance: 296.3341
##
## Estimates
## quantile scale shape
## 4.077570 0.907012 -0.002052
##
## Standard Errors
## quantile scale shape
## 0.68545 0.07391 0.07475
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 30
## Gradient Evaluations: 5
fgev(x, shape=0) # グンベルモデルを使う
##
## Call: fgev(x = x, shape = 0)
## Deviance: 296.3349
##
## Estimates
## loc scale
## -0.07617 0.90631
##
## Standard Errors
## loc scale
## 0.09539 0.07085
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 5
## Gradient Evaluations: 1
# 6.2.4項
library(ismev)
x.gev.ismev <- gev.fit(x)# GEVモデルで推定する
## $conv
## [1] 0
##
## $nllh
## [1] 148.1671
##
## $mle
## [1] -0.075240212 0.907032604 -0.002073445
##
## $se
## [1] 0.10210600 0.07393291 0.07472648
gev.diag(x.gev.ismev)

# 6.3節
library(evd) # evdパッケージをロードする
library(ismev) # ismevパッケージをロードする
set.seed(100)
x <- rexp(10000)
dim(x) <- c(50,200) # 50年分で1年あたり200個のデータとする
x[1:5,1:5] # 行列の1部がどのようなものか眺める
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.9242116 0.02620275 0.9977575 1.25242206 0.12050840
## [2,] 0.7238372 1.53172826 0.3179299 1.18030233 1.29447220
## [3,] 0.1046449 1.12569914 0.2120438 0.68628497 0.07671477
## [4,] 3.0973623 0.96308593 0.8480349 0.09184619 1.53113115
## [5,] 0.6248052 0.03491950 0.4356617 0.23421824 0.60266001
(xmax <- apply(x, 1, max)) # 各行で最大値を取る
## [1] 6.203711 4.808009 5.750058 5.561197 6.659164 6.529056 6.832202
## [8] 6.728025 3.800758 8.948393 6.941624 7.866865 4.967155 5.452945
## [15] 5.321147 6.890748 5.560308 3.830903 5.800984 5.546578 7.013092
## [22] 5.458236 5.047301 6.298012 6.219634 5.322889 8.072023 5.213717
## [29] 4.669476 6.929013 6.846274 10.823676 6.467423 8.123058 6.911492
## [36] 3.920887 6.701608 5.767668 6.757418 8.763617 7.831077 6.125121
## [43] 5.549881 5.231606 6.163502 7.513349 4.704594 4.570383 5.697831
## [50] 5.166852
(x.gev.evd <- fgev(xmax)) # evdによるGEVモデル
##
## Call: fgev(x = xmax)
## Deviance: 167.4616
##
## Estimates
## loc scale shape
## 5.6192 1.1587 -0.0801
##
## Standard Errors
## loc scale shape
## 0.18092 0.12632 0.08523
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 27
## Gradient Evaluations: 6
oldpar <- par()
par(mfrow=c(2,2)) # 4つのグラフを並べる
plot(x.gev.evd) # そのEDA

x.gev.ismev <- gev.fit(xmax) # ismevによるGEVモデル
## $conv
## [1] 0
##
## $nllh
## [1] 83.73079
##
## $mle
## [1] 5.61940312 1.15877607 -0.08010996
##
## $se
## [1] 0.18093013 0.12632839 0.08523348
gev.diag(x.gev.ismev) # そのEDA

xsort <- apply(x, 1, sort, decreasing=T) # 各行を降順にソートする
xsort[1:5,1:5] # xで行だったものがxsortで列になってしまっている
## [,1] [,2] [,3] [,4] [,5]
## [1,] 6.203711 4.808009 5.750058 5.561197 6.659164
## [2,] 5.879858 4.623220 4.126996 4.081248 6.202558
## [3,] 5.474357 4.220684 3.707195 3.965980 5.334864
## [4,] 4.712133 4.219024 3.475196 3.646878 4.687164
## [5,] 4.686745 4.172299 3.379892 3.378743 4.568518
dim(xsort)
## [1] 200 50
xsort2 <- t(xsort) # 転置を取る
xsort2[1:5, 1:5]
## [,1] [,2] [,3] [,4] [,5]
## [1,] 6.203711 5.879858 5.474357 4.712133 4.686745
## [2,] 4.808009 4.623220 4.220684 4.219024 4.172299
## [3,] 5.750058 4.126996 3.707195 3.475196 3.379892
## [4,] 5.561197 4.081248 3.965980 3.646878 3.378743
## [5,] 6.659164 6.202558 5.334864 4.687164 4.568518
x.rgev <- rlarg.fit(xsort2) # rGEVモデルを当てはめる
## $conv
## [1] 0
##
## $nllh
## [1] -32968.85
##
## $mle
## [1] 5.46469134 1.07780268 0.01712443
##
## $se
## [1] 0.108336910 0.047455518 0.009929399
# rlarg.diag(x.rgev) # xsortの当てはまりを見る 途中で抜ける
x.5gev <- rlarg.fit(xsort2,5) # 上位5個とする
## $conv
## [1] 0
##
## $nllh
## [1] 95.58674
##
## $mle
## [1] 5.573835818 1.106366565 0.002981914
##
## $se
## [1] 0.13091336 0.08522798 0.05990017
# rlarg.diag(x.5gev)
error <- rlarg.fit(xsort)
## $conv
## [1] 0
##
## $nllh
## [1] -26358.38
##
## $mle
## [1] 5.0911854 2.3533967 0.3955472
##
## $se
## [1] 0.13533996 0.10946564 0.01252596
# rlarg.diag(error)
par(oldpar)
# 6.4節
library(evd)
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
## 1893 94 87 80 72 71 71 70 69 66 65
## 1894 67 66 65 64 63 62 61 60 60 59
## 1895 105 90 87 80 80 77 77 77 77 75
## 1896 115 115 114 106 87 86 78 77 76 75
## 1897 97 89 85 81 81 77 72 72 72 72
## 1898 103 92 86 86 85 82 82 80 78 78
## 1899 90 73 72 70 68 67 66 66 65 64
## 1900 98 97 96 93 90 90 86 85 83 80
## 1901 101 92 92 92 91 87 84 80 80 78
## 1902 83 73 66 66 66 65 64 63 63 62
## 1903 130 89 88 87 85 84 83 81 79 77
## 1904 97 96 89 87 83 83 80 79 78 74
## 1905 100 87 85 83 82 80 79 78 76 75
## 1906 118 102 100 96 94 93 90 89 86 84
## 1907 89 84 84 81 73 72 72 72 70 70
## 1908 77 72 68 67 65 63 62 62 62 61
## 1909 95 73 73 71 71 68 67 67 66 66
## 1910 115 105 95 90 89 87 86 83 82 81
## 1911 95 89 88 86 86 85 82 78 77 77
## 1912 87 87 84 83 81 79 74 74 73 73
## 1913 95 85 83 78 75 72 71 71 70 69
## 1914 114 94 93 82 82 80 79 79 78 76
## 1915 99 94 92 90 87 86 82 78 78 78
## 1916 136 102 101 101 99 94 94 93 93 89
## 1917 107 103 102 93 92 89 89 86 85 83
## 1918 82 81 78 77 77 77 76 76 72 71
## 1919 106 97 95 92 84 83 82 82 81 81
## 1920 111 100 85 84 84 83 82 80 74 71
## 1921 65 65 63 63 61 60 60 59 58 58
## 1922 106 NA NA NA NA NA NA NA NA NA
## 1923 105 97 81 80 78 75 73 73 72 71
## 1924 70 70 70 69 68 68 66 65 65 64
## 1925 96 92 92 87 86 86 79 79 77 76
## 1926 94 92 92 90 90 87 87 86 84 84
## 1927 111 108 102 100 95 90 88 84 80 80
## 1928 110 97 93 91 90 90 89 84 82 80
## 1929 90 89 88 78 78 75 73 73 70 70
## 1930 100 90 86 85 85 84 83 81 79 78
## 1931 103 99 98 96 94 89 86 85 84 79
## 1932 78 78 74 73 73 72 71 70 70 69
## 1933 121 113 106 105 102 89 89 88 86 85
## 1934 116 113 91 91 91 89 88 88 86 81
## 1935 115 107 105 101 93 91 NA NA NA NA
## 1936 147 106 93 90 87 87 87 84 82 81
## 1937 119 107 107 106 105 102 98 95 94 94
## 1938 114 97 85 83 82 81 79 76 74 69
## 1939 89 86 82 81 80 80 78 78 77 77
## 1940 102 101 98 97 96 94 94 91 90 89
## 1941 99 98 96 95 94 94 89 88 87 86
## 1942 91 91 87 83 83 81 78 75 75 73
## 1943 97 88 82 79 78 78 76 76 74 71
## 1944 106 95 94 90 89 84 84 83 82 80
## 1945 105 102 98 88 86 84 84 80 80 80
## 1946 136 104 103 101 100 91 88 84 83 82
## 1947 126 108 101 99 98 98 96 96 94 93
## 1948 132 126 119 107 101 98 92 86 85 83
## 1949 104 102 102 101 93 92 90 88 88 86
## 1950 117 96 91 89 88 86 86 86 86 85
## 1951 151 117 114 109 106 104 100 99 99 99
## 1952 116 104 103 98 91 91 90 90 87 86
## 1953 107 102 98 98 92 89 89 89 84 84
## 1954 112 100 95 94 94 90 86 82 81 81
## 1955 97 96 96 95 94 92 90 90 90 88
## 1956 95 91 90 85 85 84 82 82 82 80
## 1957 119 107 100 98 98 97 93 92 91 90
## 1958 124 114 113 110 108 104 104 102 100 99
## 1959 118 117 108 107 105 102 96 96 94 94
## 1960 145 126 123 116 114 110 108 108 107 106
## 1961 122 108 104 100 100 95 94 93 92 91
## 1962 114 110 108 107 106 104 99 99 98 95
## 1963 118 116 114 112 110 109 104 103 103 102
## 1964 107 104 104 103 102 98 92 90 90 90
## 1965 110 108 106 102 101 101 100 98 98 98
## 1966 194 127 126 104 103 102 102 99 98 97
## 1967 138 118 118 107 100 96 95 93 92 92
## 1968 144 132 123 114 112 110 108 107 107 106
## 1969 138 120 116 114 108 106 104 103 103 103
## 1970 123 122 119 110 105 99 99 97 96 96
## 1971 122 116 116 109 104 101 100 100 99 98
## 1972 120 118 113 111 96 92 91 91 90 89
## 1973 114 111 99 98 97 97 96 92 90 90
## 1974 96 95 95 93 92 90 90 89 88 88
## 1975 125 110 109 103 102 101 101 97 90 89
## 1976 124 122 114 109 108 108 104 104 102 100
## 1977 120 102 100 98 96 96 95 94 91 90
## 1978 132 114 110 107 105 102 100 100 100 99
## 1979 166 140 131 130 122 118 116 115 115 112
## 1980 134 114 111 109 107 106 104 103 102 99
## 1981 138 136 130 128 119 110 107 104 104 104
## 1982 135 128 120 117 114 113 111 108 106 104
## 1983 122 105 100 100 100 99 99 98 98 98
## 1984 119 118 115 111 105 105 103 102 101 101
## 1985 124 103 103 95 95 94 91 91 91 91
## 1986 159 114 103 99 99 97 96 95 95 93
## 1987 137 130 109 107 100 100 100 97 97 97
## 1988 103 97 92 91 87 86 86 86 84 84
## 1989 102 99 97 95 94 91 91 91 89 89
## 1990 128 116 113 111 111 107 104 103 99 95
## 1991 127 107 105 103 101 95 95 95 94 92
## 1992 142 132 124 122 118 113 110 105 103 102
## 1993 123 112 111 110 103 102 100 99 97 94
## 1994 111 105 100 100 98 96 95 91 90 90
## 1995 112 107 104 100 100 95 94 93 91 91
## 1996 134 128 115 115 111 111 111 110 110 107
## 1997 126 123 120 118 112 111 111 110 107 104
## 1998 124 113 112 112 105 104 104 104 102 102
## 1999 121 116 116 112 110 107 105 104 103 103
## 2000 144 124 116 111 111 110 108 108 103 103
## 2001 124 115 112 111 110 108 107 103 103 103
## 2002 147 126 123 123 122 121 115 113 113 111
## 2003 109 108 107 102 101 100 100 100 98 98
## 2004 137 128 124 115 115 108 107 106 104 104
## 2005 132 101 99 97 96 95 95 94 94 93
## 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
# ?venice2
dim(venice2)
## [1] 125 10
plot(1887:2011,venice2[,1],type="l") # 折れ線グラフでデータの様子を眺める

(venice.gev <- fgev(venice2[,1]))
##
## Call: fgev(x = venice2[, 1])
## Deviance: 1111.223
##
## Estimates
## loc scale shape
## 105.2995 19.3543 -0.1463
##
## Standard Errors
## loc scale shape
## 1.87769 1.27804 0.04176
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 27
## Gradient Evaluations: 11
oldpar <- par()
par(mfrow=c(2,2)) # 4つのグラフを並べる
plot(venice.gev)

fgev(venice2[,1], prob=0.01)
##
## Call: fgev(x = venice2[, 1], prob = 0.01)
## Deviance: 1111.223
##
## Estimates
## quantile scale shape
## 170.0808 19.3711 -0.1469
##
## Standard Errors
## quantile scale shape
## 5.47737 1.28114 0.04166
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 44
## Gradient Evaluations: 22
library(ismev)
venice.5gev <- rlarg.fit(venice2,5)
## $conv
## [1] 0
##
## $nllh
## [1] 1850.09
##
## $mle
## [1] 116.5859859 14.9862451 -0.1542552
##
## $se
## [1] 1.10176039 0.44165912 0.01582627
# rlarg.diag(venice.5gev)
venice.4gev <- rlarg.fit(venice2,4)
## $conv
## [1] 0
##
## $nllh
## [1] 1593.704
##
## $mle
## [1] 115.5945120 15.7247156 -0.1594054
##
## $se
## [1] 1.1898084 0.4998447 0.0177804
# rlarg.diag(venice.4gev)
venice.3gev <- rlarg.fit(venice2,3)
## $conv
## [1] 0
##
## $nllh
## [1] 1296.104
##
## $mle
## [1] 113.7315222 16.4475401 -0.1577057
##
## $se
## [1] 1.28758998 0.59582815 0.02149216
# rlarg.diag(venice.3gev)
venice.2gev <- rlarg.fit(venice2,2)
## $conv
## [1] 0
##
## $nllh
## [1] 962.4824
##
## $mle
## [1] 110.4548562 17.6343495 -0.1535678
##
## $se
## [1] 1.45851405 0.78571744 0.02847785
# rlarg.diag(venice.2gev)
fox
## berlin wright
## 1918 6.05 16.3
## 1919 2.67 13.1
## 1920 5.15 16.6
## 1921 2.45 14.2
## 1922 5.92 20.1
## 1923 6.05 13.7
## 1924 4.02 15.5
## 1925 2.52 8.3
## 1926 3.44 9.1
## 1927 3.17 13.3
## 1928 5.92 15.1
## 1929 6.62 20.6
## 1930 3.00 6.6
## 1931 1.14 3.1
## 1932 1.91 9.9
## 1933 2.60 8.9
## 1934 1.91 6.7
## 1935 4.34 11.1
## 1936 4.34 6.3
## 1937 3.26 13.5
## 1938 6.19 18.0
## 1939 4.91 18.2
## 1940 4.72 17.5
## 1941 3.54 16.6
## 1942 2.74 19.8
## 1943 5.08 21.3
## 1944 2.29 10.8
## 1945 3.46 15.8
## 1946 6.90 21.3
## 1947 3.16 11.0
## 1948 4.54 10.3
## 1949 2.00 6.4
## 1950 4.63 10.9
# ?fox
dim(fox)
## [1] 33 2
(fox.gev <- fgev(fox[,1]))
##
## Call: fgev(x = fox[, 1])
## Deviance: 120.806
##
## Estimates
## loc scale shape
## 3.3804 1.4492 -0.2317
##
## Standard Errors
## loc scale shape
## 0.3030 0.2332 0.1949
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 26
## Gradient Evaluations: 9
plot(fox.gev)

data(glass)
glass
## [1] 0.55 0.74 0.77 0.81 0.84 0.93 1.04 1.11 1.13 1.24 1.25 1.27 1.28 1.29 1.30
## [16] 1.36 1.39 1.42 1.48 1.48 1.49 1.49 1.50 1.50 1.51 1.52 1.53 1.54 1.55 1.55
## [31] 1.58 1.59 1.60 1.61 1.61 1.61 1.61 1.62 1.62 1.63 1.64 1.66 1.66 1.66 1.67
## [46] 1.68 1.68 1.69 1.70 1.70 1.73 1.76 1.76 1.77 1.78 1.81 1.82 1.84 1.84 1.89
## [61] 2.00 2.01 2.24
(glass.gev <- fgev(glass))
##
## Call: fgev(x = glass)
## Deviance: 34.2449
##
## Estimates
## loc scale shape
## 1.4080 0.3475 -0.3860
##
## Standard Errors
## loc scale shape
## 0.04683 0.03149 0.04972
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 58
## Gradient Evaluations: 15
par(mfrow=c(2,2))
plot(glass.gev)

par(oldpar)
# 6.5.1項
getwd() # 作業ディレクトリの確認
## [1] "C:/Users/saigo/Documents/R/EVTbook"
Tokyo <- read.csv("TokyoTemp.csv")
head(Tokyo)
## date temp 品質情報 均質番号
## 1 1990/1/1 9.8 8 1
## 2 1990/1/2 9.4 8 1
## 3 1990/1/3 9.1 8 1
## 4 1990/1/4 9.3 8 1
## 5 1990/1/5 10.6 8 1
## 6 1990/1/6 16.3 8 1
summer.n <- grep("/[789]/", Tokyo[,1]) # 7,8,9月の行番号を取り出す
length(summer.n); summer.n[1:20]
## [1] 2668
## [1] 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## [20] 201
Tokyo2 <- Tokyo[summer.n,2] # 7,8,9月の気温を取り出す
head(Tokyo2)
## [1] 22.1 25.9 25.3 25.7 30.7 29.6
2668/29
## [1] 92
dim(Tokyo2) <- c(92,29) # 行列に整形する
colnames(Tokyo2) <- 1990:2018 # 各列に名前を付ける
head(Tokyo2)
## 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004
## [1,] 22.1 27.8 30.5 20.7 27.6 26.3 30.7 30.2 30.0 28.1 33.5 36.7 25.5 23.6 29.2
## [2,] 25.9 26.7 24.8 21.0 32.4 27.6 31.2 29.7 32.9 29.0 33.1 30.5 27.2 26.8 28.0
## [3,] 25.3 26.1 26.1 27.0 35.6 27.2 31.0 34.0 34.4 27.1 32.8 32.5 26.7 26.7 29.5
## [4,] 25.7 28.6 29.4 26.6 32.6 24.6 29.7 34.4 36.1 31.0 31.7 35.3 31.5 29.3 29.5
## [5,] 30.7 31.5 27.7 21.6 35.3 22.9 25.8 37.7 31.6 24.4 28.3 34.9 28.4 29.0 29.4
## [6,] 29.6 24.5 27.2 25.8 33.1 23.3 31.9 35.8 29.2 23.6 28.5 33.0 30.6 24.6 32.3
## 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018
## [1,] 27.9 27.9 26.9 25.5 25.7 29.5 33.9 24.7 26.1 29.6 22.4 30.6 25.3 32.0
## [2,] 29.4 29.9 25.6 26.9 22.6 31.0 28.9 27.0 27.9 30.3 25.0 31.8 32.0 33.5
## [3,] 22.7 29.4 27.1 27.3 26.3 28.8 31.1 27.7 27.4 27.5 23.4 35.4 32.5 32.7
## [4,] 22.9 29.3 22.7 31.6 27.6 31.6 34.3 29.0 27.3 23.4 25.9 33.8 29.9 31.2
## [5,] 29.3 24.1 30.4 31.4 26.3 30.3 32.3 29.4 28.9 24.3 21.9 25.6 31.4 28.7
## [6,] 23.7 28.3 30.1 30.7 24.7 30.6 32.5 28.8 33.7 28.5 21.1 26.9 30.7 25.0
Tokyo3 <- apply(Tokyo2, 2, sort, decreasing=T) # 各列でソートする
head(Tokyo3)
## 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004
## [1,] 35.9 35.6 35.2 32.9 39.1 36.4 38.7 37.7 36.1 34.8 37.8 38.1 35.8 34.3 39.5
## [2,] 35.9 35.4 35.1 32.9 36.6 36.3 36.6 35.8 35.5 34.2 34.9 36.7 35.7 34.3 38.1
## [3,] 34.6 35.2 34.7 32.3 35.9 35.9 34.8 35.0 35.2 34.2 34.5 36.3 35.7 34.1 36.5
## [4,] 34.5 35.0 34.2 32.3 35.6 35.8 34.5 34.5 34.4 34.2 34.5 36.1 35.6 34.1 36.3
## [5,] 34.3 34.6 34.1 32.3 35.5 35.6 34.5 34.4 34.4 34.2 34.5 35.7 35.3 33.9 35.1
## [6,] 34.2 34.5 33.8 32.1 35.3 35.4 33.7 34.3 34.2 34.2 34.1 35.6 35.2 33.7 35.1
## 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018
## [1,] 35.8 36.1 37.5 35.3 34.2 37.2 36.1 35.7 38.3 36.1 37.7 37.7 37.1 39.0
## [2,] 35.6 35.4 37.0 34.9 33.9 36.3 35.2 35.6 37.4 35.7 35.9 36.7 35.0 37.3
## [3,] 35.0 35.0 37.0 34.8 33.8 36.3 35.1 35.4 36.8 35.6 35.8 35.4 34.9 36.0
## [4,] 34.9 34.8 36.4 34.5 33.6 36.1 34.8 35.2 35.8 35.5 35.5 34.3 34.9 35.8
## [5,] 34.8 34.4 35.9 34.5 33.6 35.9 34.6 35.1 35.7 35.3 35.3 34.2 34.8 35.6
## [6,] 34.3 34.3 35.7 34.5 33.2 35.9 34.5 35.0 35.4 34.7 35.2 34.0 33.9 35.6
tTokyo3 <- t(Tokyo3) # 関数の仕様に合わせて転置する
plot(1990:2018, tTokyo3[,1], type="l") # 各年の最高気温の折れ線グラフ

library(ismev) # ismevパッケージをロードする
Tokyo_model_r1 <- rlarg.fit(tTokyo3,1) # r=1のモデルを当てる
## $conv
## [1] 0
##
## $nllh
## [1] 53.94325
##
## $mle
## [1] 36.142854 1.652356 -0.399661
##
## $se
## [1] 0.3437395 0.2590920 0.1505488
# rlarg.diag(Tokyo_model_r1) # モデルの適合性を図示する
Tokyo_model_r2 <- rlarg.fit(tTokyo3,2)
## $conv
## [1] 0
##
## $nllh
## [1] 72.30083
##
## $mle
## [1] 36.3428173 1.2503843 -0.2383249
##
## $se
## [1] 0.21884633 0.11282737 0.09497959
# rlarg.diag(Tokyo_model_r2)
library(evd)
fgev(tTokyo3[,1], prob=0.01)
##
## Call: fgev(x = tTokyo3[, 1], prob = 0.01)
## Deviance: 107.8865
##
## Estimates
## quantile scale shape
## 39.6198 1.6523 -0.3996
##
## Standard Errors
## quantile scale shape
## 0.5016 0.2591 0.1506
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 35
## Gradient Evaluations: 15
# 6.5.2項
rain <- read.csv("TateyamaAme.csv")
dim(rain)
## [1] 838 5
rain2 <- na.omit(rain) # 初めの方はデータがないようなので削除する
head(rain2)
## X rainfall 現象なし情報 品質情報 均質番号
## 221 May-68 4 0 4 1
## 222 Jun-68 8 0 8 1
## 223 Jul-68 4 0 8 1
## 224 Aug-68 16 0 8 1
## 225 Sep-68 10 0 8 1
## 226 Oct-68 6 0 8 1
raints <- ts(rain2[,2],start=c(1968,5),frequency=12) # 時系列の形式に直す
raints
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1968 4.0 8.0 4.0 16.0 10.0 6.0 3.5 8.5
## 1969 1.5 13.0 7.0 6.0 4.5 11.0 5.0 18.5 6.0 18.0 6.0 3.5
## 1970 5.5 4.0 4.0 3.0 4.5 6.5 6.0 7.0 17.5 5.0 7.0 17.5
## 1971 2.5 2.5 7.0 3.5 4.0 4.0 5.0 14.5 12.5 7.0 0.5 4.0
## 1972 4.0 9.5 1.5 7.0 5.0 6.5 8.0 7.5 20.0 3.5 6.5 9.0
## 1973 5.5 3.0 2.0 8.5 6.0 3.5 8.0 5.5 11.5 15.0 6.0 0.5
## 1974 6.5 2.0 6.5 4.0 5.5 8.5 10.5 6.5 9.0 4.5 6.0 4.5
## 1975 9.5 5.0 13.5 2.5 4.0 8.5 5.5 8.0 8.5 5.0 15.0 5.0
## 1976 1.5 6.5 8.0 5.5 5.0 12.5 5.5 7.0 7.5 14.5 10.5 6.5
## 1977 2.5 1.0 6.0 6.0 4.0 8.5 15.5 8.5 7.5 3.0 4.0 5.5
## 1978 1.0 2.0 4.0 3.5 7.0 5.0 3.5 2.0 8.0 4.0 3.0 4.5
## 1979 9.5 10.5 5.0 4.0 10.0 11.0 5.0 13.5 5.0 18.5 10.5 3.0
## 1980 5.0 4.0 5.0 4.0 6.5 9.5 9.0 10.5 7.0 12.5 4.0 7.5
## 1981 4.0 3.5 4.5 8.0 8.5 8.0 24.0 5.0 3.5 7.0 7.0 2.0
## 1982 7.0 4.0 2.0 12.0 6.0 6.5 7.5 6.0 14.0 3.5 6.0 3.5
## 1983 4.5 2.0 2.5 4.5 3.0 4.5 15.5 7.0 3.5 10.0 4.5 6.0
## 1984 1.5 4.0 3.0 2.5 8.5 7.0 7.5 9.0 3.5 7.5 2.5 4.5
## 1985 1.5 6.0 3.0 7.0 4.0 6.0 5.0 9.0 5.0 3.5 10.0 1.5
## 1986 7.5 1.0 5.0 3.0 3.5 5.0 14.0 11.0 12.5 3.5 5.5 7.5
## 1987 4.0 2.0 3.0 2.5 5.0 6.0 11.5 5.0 10.5 6.0 11.0 3.0
## 1988 2.5 4.0 3.5 6.0 3.5 7.5 5.5 8.5 14.5 6.0 1.0 0.5
## 1989 5.0 4.5 6.0 4.0 7.5 8.0 10.0 23.5 9.5 7.0 12.5 2.0
## 1990 7.0 3.5 7.0 11.5 2.5 5.0 6.5 7.0 10.0 14.0 4.5 7.0
## 1991 3.5 3.0 2.5 4.0 2.5 11.0 6.0 9.5 7.0 7.0 7.0 4.0
## 1992 1.5 1.5 4.0 6.5 7.0 7.0 3.5 6.0 6.5 16.0 6.5 4.5
## 1993 2.5 9.0 3.0 3.5 2.5 14.0 13.5 12.5 5.0 13.0 15.0 10.5
## 1994 1.0 6.0 8.5 7.5 4.5 5.5 1.5 14.0 17.0 6.0 12.5 3.5
## 1995 9.0 2.0 5.5 10.5 10.0 14.5 13.0 2.0 3.5 12.0 15.0 1.5
## 1996 7.0 3.0 6.5 3.5 5.0 7.0 7.5 6.5 13.0 8.0 3.0 5.0
## 1997 2.5 6.5 5.5 7.5 11.5 4.5 8.0 6.0 5.5 7.0 4.5 6.5
## 1998 2.0 6.0 3.5 10.5 8.5 6.0 7.0 3.5 9.0 7.5 11.5 2.0
## 1999 2.5 2.0 4.0 4.0 2.5 2.5 15.5 13.0 8.0 16.5 4.0 0.5
## 2000 8.0 2.5 16.5 4.0 5.5 9.5 8.0 15.0 6.5 4.5 6.5 3.5
## 2001 11.0 4.5 3.0 3.5 7.5 3.0 1.5 8.0 9.5 8.5 9.5 5.0
## 2002 15.0 1.5 3.5 1.5 3.5 7.5 7.0 20.0 11.5 14.5 2.5 24.0
## 2003 5.0 4.0 3.0 4.0 3.5 3.0 10.5 7.5 4.0 4.0 18.5 2.0
## 2004 2.5 3.5 5.0 3.5 3.0 7.0 8.5 8.0 11.5 21.0 10.5 4.0
## 2005 2.0 7.5 4.5 4.5 5.0 6.5 7.0 11.5 6.5 9.0 7.0 0.5
## 2006 15.5 11.0 4.0 9.5 5.0 5.0 9.5 14.5 20.0 7.5 2.5 9.5
## 2007 4.5 8.0 16.0 10.0 9.0 5.0 20.0 2.5 10.0 7.0 2.0 6.0
## 2008 1.5 2.0 5.5 8.0 5.0 9.5 5.0 15.0 8.5 20.0 5.5 4.5
## 2009 7.5 6.0 5.0 3.0 23.0 6.0 10.5 11.5 5.5 5.0 8.0 7.0
## 2010 14.5 8.0 3.5 9.5 3.5 4.0 11.5 4.5 7.0 11.0 13.5 7.5
## 2011 0.5 4.0 4.0 11.5 7.0 4.0 13.5 7.5 11.0 11.5 7.5 12.0
## 2012 2.0 3.0 3.0 8.5 6.0 8.5 8.0 2.0 8.0 14.5 12.0 6.5
## 2013 5.0 5.0 5.0 7.5 3.0 4.5 6.5 2.5 13.0 13.5 14.0 3.5
## 2014 9.0 4.0 6.0 4.0 7.5 5.5 3.0 6.5 11.5 15.5 4.0 12.0
## 2015 4.0 3.5 7.0 4.5 4.0 8.5 6.5 14.0 24.5 10.5 5.5 3.0
## 2016 2.0 6.5 6.0 5.0 4.5 16.0 7.5 11.0 12.5 2.5 8.0 5.0
## 2017 3.0 1.5 2.0 7.5 8.0 4.0 12.0 16.5 14.0 7.0 16.0 2.5
## 2018 5.0 1.5 7.0 3.5 5.5 9.5 14.0 7.0 12.0 8.5 9.0 3.0
## 2019 1.5 3.5 3.0 5.5 6.0 9.0 10.0 7.5 16.0 18.0
plot(raints)

wind <- read.csv("TateyamaKaze.csv")
dim(wind)
## [1] 838 6
wind2 <- na.omit(wind) # 初めの方はデータがないようなので削除する
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) # 時系列の形式に直す
windts
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1968 9.8 11.3 13.0 12.8 7.5 11.8 13.5 12.3
## 1969 13.3 10.5 14.7 15.3 11.8 13.3 13.2 13.5 11.2 9.5 10.5 14.8
## 1970 16.8 14.5 13.5 14.7 11.3 10.3 12.8 10.2 9.5 8.7 14.5 13.0
## 1971 12.0 11.0 11.8 12.3 12.2 12.5 11.7 11.7 17.7 12.3 12.5 12.0
## 1972 13.3 13.2 14.7 12.8 13.0 12.7 12.8 14.7 15.0 11.5 13.5 14.2
## 1973 14.0 14.3 12.3 12.8 10.7 10.0 9.3 10.2 9.8 11.8 13.0 13.0
## 1974 13.2 16.5 16.7 12.7 10.8 10.8 10.8 10.7 8.8 9.7 13.2 11.3
## 1975 10.7 12.1 10.9 12.4 11.8 10.5 12.2 10.8 8.5 10.9 10.1 10.1
## 1976 11.4 10.4 10.8 11.5 10.4 10.5 9.7 9.4 8.7 13.0 12.5 11.4
## 1977 9.6 11.4 10.9 11.5 10.6 7.6 10.0 8.8 11.4 7.4 9.9 10.0
## 1978 11.1 11.0 12.7 12.3 9.5 10.9 6.7 9.3 8.4 8.9 11.6 10.3
## 1979 11.2 10.5 12.2 12.0 11.6 10.4 10.5 10.2 8.5 18.1 9.8 9.0
## 1980 11.5 11.3 11.0 12.4 11.5 11.2 10.0 10.4 8.1 13.2 9.9 11.2
## 1981 11.0 10.5 13.0 13.2 10.6 10.5 10.9 15.2 9.2 13.5 10.3 11.5
## 1982 11.6 11.0 12.7 12.1 11.5 13.1 10.6 16.3 16.1 12.6 13.9 13.2
## 1983 10.3 11.6 12.7 13.1 11.6 12.7 11.7 11.2 10.3 9.9 13.9 11.3
## 1984 9.9 10.3 12.0 12.1 11.5 11.2 7.6 9.3 10.7 10.3 10.5 9.1
## 1985 11.5 10.7 10.8 10.5 8.7 10.6 18.4 15.5 9.7 10.1 10.3 11.0
## 1986 11.1 11.4 12.6 11.3 12.5 13.4 9.3 8.9 10.0 10.0 7.6 12.7
## 1987 11.8 12.5 14.3 11.3 12.7 10.0 11.1 8.8 10.6 11.2 10.0 10.7
## 1988 11.2 12.2 12.0 12.1 10.1 12.6 11.7 8.8 9.3 9.9 13.3 11.0
## 1989 10.5 9.2 11.0 12.0 11.7 8.7 10.8 12.4 13.8 8.5 10.9 10.7
## 1990 11.2 12.8 12.4 11.8 10.5 10.6 10.0 11.8 15.7 12.0 12.8 12.6
## 1991 9.9 12.5 12.0 13.9 10.5 11.0 12.3 10.5 11.8 9.0 12.2 11.6
## 1992 11.5 11.4 11.0 12.2 10.3 9.0 11.7 11.0 8.6 10.8 10.4 11.3
## 1993 11.8 12.1 12.9 14.6 13.8 10.2 12.2 15.2 10.6 10.4 10.8 11.6
## 1994 13.1 15.6 11.0 11.5 9.7 9.6 9.7 7.5 9.7 10.0 8.2 11.4
## 1995 10.8 8.9 10.3 13.2 9.9 9.6 10.7 9.8 17.9 10.3 12.3 13.9
## 1996 11.2 10.7 11.2 10.9 11.1 12.0 10.5 12.6 20.5 11.0 10.2 12.0
## 1997 12.1 11.7 11.7 10.4 11.3 13.9 12.8 11.0 9.7 10.6 10.9 10.8
## 1998 11.3 12.3 12.5 12.5 11.5 12.2 9.0 10.3 17.8 12.4 9.9 10.4
## 1999 9.5 11.4 10.0 12.0 12.9 13.2 8.9 8.0 9.8 11.8 12.4 10.0
## 2000 10.2 11.7 13.0 11.8 11.2 12.7 16.0 8.5 9.1 10.5 12.7 12.8
## 2001 12.5 8.9 11.8 9.7 9.1 12.9 10.0 11.5 16.8 9.7 9.8 13.1
## 2002 12.1 10.5 11.8 13.3 8.9 11.2 12.2 9.6 7.7 18.6 11.0 11.1
## 2003 13.4 9.3 12.0 11.7 11.0 11.5 11.0 10.6 10.9 10.1 11.0 11.9
## 2004 11.9 12.1 12.2 12.9 10.1 10.7 9.8 12.4 12.5 15.6 12.4 17.7
## 2005 11.0 12.3 11.5 12.4 9.6 11.4 12.6 18.0 10.6 9.8 9.3 12.3
## 2006 10.9 12.0 14.1 11.0 13.3 8.7 12.6 9.8 7.3 11.7 10.5 10.7
## 2007 14.2 11.8 11.6 12.5 10.0 10.5 10.1 10.1 18.2 14.9 9.8 11.3
## 2008 9.7 12.2 10.2 13.5 12.9 8.9 8.8 9.6 11.3 9.3 11.5 11.6
## 2009 10.7 12.8 14.1 15.3 12.8 11.8 14.4 14.0 9.7 15.1 12.7 15.9
## 2010 13.9 12.9 14.3 13.3 10.4 13.2 10.6 10.8 11.5 10.4 13.6 15.3
## 2011 12.5 11.6 13.4 13.0 12.8 13.4 13.0 10.9 17.3 11.6 14.4 10.4
## 2012 11.2 13.4 12.3 14.1 11.2 17.6 10.8 9.8 14.4 12.4 13.5 14.1
## 2013 15.4 12.3 12.9 16.7 12.0 12.5 12.8 12.5 15.7 20.1 14.4 14.6
## 2014 14.6 12.3 14.9 10.4 11.9 12.2 12.2 11.0 11.0 21.2 11.7 14.1
## 2015 13.0 11.7 14.5 12.7 13.4 10.7 13.8 10.3 11.0 13.0 13.7 16.9
## 2016 13.9 12.9 10.0 13.8 13.4 12.3 12.7 17.7 12.6 13.7 8.8 13.4
## 2017 13.1 16.5 11.0 13.1 14.0 13.3 10.4 11.6 14.0 19.9 12.2 13.6
## 2018 15.4 13.7 15.9 13.4 13.7 12.0 14.9 12.9 13.1 19.8 7.8 10.8
## 2019 12.3 13.0 14.1 13.7 11.7 12.7 10.3 11.7 28.4 20.7
plot(windts)

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
tail(wind2) # データの末尾
## 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
dim(wind2) # データの次元
## [1] 618 6
618/12
## [1] 51.5
windmat <- matrix(c(rep(NA,4), wind2[,2], rep(NA,2)), 52, 12, byrow=T) # データを行列形式にする
rownames(windmat) <- 1968:2019 # 行の名前を年の名前にする
(windmax <- apply(windmat, 1, max)) # 各年の最大値を取る NAができてしまう
## 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983
## NA 15.3 16.8 17.7 15.0 14.3 16.7 12.4 13.0 11.5 12.7 18.1 13.2 15.2 16.3 13.9
## 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999
## 12.1 18.4 13.4 14.3 13.3 13.8 15.7 13.9 12.2 15.2 15.6 17.9 20.5 13.9 17.8 13.2
## 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015
## 16.0 16.8 18.6 13.4 17.7 18.0 14.1 18.2 13.5 15.9 15.3 17.3 17.6 20.1 21.2 16.9
## 2016 2017 2018 2019
## 17.7 19.9 19.8 NA
(windmax <- apply(windmat, 1, max, na.rm=T)) # 各年の最大値を取り直す NAなし
## 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983
## 13.5 15.3 16.8 17.7 15.0 14.3 16.7 12.4 13.0 11.5 12.7 18.1 13.2 15.2 16.3 13.9
## 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999
## 12.1 18.4 13.4 14.3 13.3 13.8 15.7 13.9 12.2 15.2 15.6 17.9 20.5 13.9 17.8 13.2
## 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015
## 16.0 16.8 18.6 13.4 17.7 18.0 14.1 18.2 13.5 15.9 15.3 17.3 17.6 20.1 21.2 16.9
## 2016 2017 2018 2019
## 17.7 19.9 19.8 28.4
plot(1968:2019,windmax, type="l")

length(windmax)
## [1] 52
windmax18 <- windmax[1:51] # 2018年までのデータ
library(evd)
(wind.gev18 <- fgev(windmax18)) # evdによる推定
##
## Call: fgev(x = windmax18)
## Deviance: 233.425
##
## Estimates
## loc scale shape
## 14.8066 2.2299 -0.1784
##
## Standard Errors
## loc scale shape
## 0.3679 0.2768 0.1400
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 18
## Gradient Evaluations: 8
oldpar <- par()
par(mfrow=c(2,2))
plot(wind.gev18)

str(wind.gev18)
## List of 17
## $ estimate : Named num [1:3] 14.807 2.23 -0.178
## ..- attr(*, "names")= chr [1:3] "loc" "scale" "shape"
## $ std.err : Named num [1:3] 0.368 0.277 0.14
## ..- attr(*, "names")= chr [1:3] "loc" "scale" "shape"
## $ fixed : NULL
## $ param : Named num [1:3] 14.807 2.23 -0.178
## ..- attr(*, "names")= chr [1:3] "loc" "scale" "shape"
## $ deviance : num 233
## $ corr : NULL
## $ var.cov : num [1:3, 1:3] 0.1354 0.0293 -0.025 0.0293 0.0766 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "loc" "scale" "shape"
## .. ..$ : NULL
## $ convergence: chr "successful"
## $ counts : Named int [1:2] 18 8
## ..- attr(*, "names")= chr [1:2] "function" "gradient"
## $ message : NULL
## $ data : num [1:51] 13.5 15.3 16.8 17.7 15 14.3 16.7 12.4 13 11.5 ...
## $ tdata : num [1:51] 13.5 15.3 16.8 17.7 15 14.3 16.7 12.4 13 11.5 ...
## $ nsloc : NULL
## $ n : int 51
## $ prob : NULL
## $ loc : Named num 14.8
## ..- attr(*, "names")= chr "loc"
## $ call : language fgev(x = windmax18)
## - attr(*, "class")= chr [1:3] "gev" "uvevd" "evd"
(temp <- wind.gev18$estimate)
## loc scale shape
## 14.8065593 2.2298747 -0.1784235
1-pgev(28.4, loc=temp[1], scale=temp[2], shape=temp[3]) # 上側確率
## shape
## 0
wind3 <- wind2[wind2[,6]==1,] # 6列目の均質番号が1の行だけ取り出す
head(wind3); tail(wind3); dim(wind3)
## 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
## X kaze 品質情報 X.1 品質情報.1 均質番号
## 696 Dec-07 11.3 8 西南西 8 1
## 697 Jan-08 9.7 8 西 8 1
## 698 Feb-08 12.2 8 西南西 8 1
## 699 Mar-08 10.2 8 北 8 1
## 700 Apr-08 13.5 8 東 8 1
## 701 May-08 12.9 8 南南東 8 1
## [1] 481 6
481/12
## [1] 40.08333
windmat2 <- matrix(c(rep(NA,4), wind3[,2], rep(NA,7)), 41, 12, byrow=T) # データを行列形式にする
rownames(windmat2) <- 1968:2008 # 行の名前を年の名前にする
(windmax2 <- apply(windmat2, 1, max, na.rm=T)) # 各年の最大値を取り直す NAなし
## 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983
## 13.5 15.3 16.8 17.7 15.0 14.3 16.7 12.4 13.0 11.5 12.7 18.1 13.2 15.2 16.3 13.9
## 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999
## 12.1 18.4 13.4 14.3 13.3 13.8 15.7 13.9 12.2 15.2 15.6 17.9 20.5 13.9 17.8 13.2
## 2000 2001 2002 2003 2004 2005 2006 2007 2008
## 16.0 16.8 18.6 13.4 17.7 18.0 14.1 18.2 13.5
(wind.gev08 <- fgev(windmax2))
##
## Call: fgev(x = windmax2)
## Deviance: 178.0132
##
## Estimates
## loc scale shape
## 14.250 1.895 -0.100
##
## Standard Errors
## loc scale shape
## 0.3554 0.2701 0.1691
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 22
## Gradient Evaluations: 8
plot(wind.gev08)

(wind.gev18gum <- fgev(windmax18, shape=0))
##
## Call: fgev(x = windmax18, shape = 0)
## Deviance: 234.7833
##
## Estimates
## loc scale
## 14.601 2.092
##
## Standard Errors
## loc scale
## 0.3096 0.2294
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 12
## Gradient Evaluations: 6
plot(wind.gev18gum)

temp <- wind.gev18gum$estimate
1-pgev(28.4, loc=temp[1], scale=temp[2], shape=0)
## loc
## 0.001364422
1/(1-pgev(28.4, loc=temp[1], scale=temp[2], shape=0))
## loc
## 732.9109
par(oldpar)
# 6.5.3項
olm <- read.csv("results.csv")
dim(olm)
## [1] 2406 8
head(olm)# 初めの6行を見る
## Gender Event Location Year Medal Name Nationality
## 1 M 10000M Men Rio 2016 G Mohamed FARAH USA
## 2 M 10000M Men Rio 2016 S Paul Kipngetich TANUI KEN
## 3 M 10000M Men Rio 2016 B Tamirat TOLA ETH
## 4 M 10000M Men Beijing 2008 G Kenenisa BEKELE ETH
## 5 M 10000M Men Beijing 2008 S Sileshi SIHINE ETH
## 6 M 10000M Men Beijing 2008 B Micah KOGO KEN
## Result
## 1 25:05.17
## 2 27:05.64
## 3 27:06.26
## 4 27:01.17
## 5 27:02.77
## 6 27:04.11
levels(olm[,2])# 2列目が競技なのでその要素の種類を並べる
## [1] "" "10000M Men"
## [3] "10000M Women" "100M Hurdles Women"
## [5] "100M Men" "100M Women"
## [7] "110M Hurdles Men" "1500M Men"
## [9] "1500M Women" "200M Men"
## [11] "200M Women" "20Km Race Walk Men"
## [13] "20Km Race Walk Women" "3000M Steeplechase Men"
## [15] "3000M Steeplechase Women" "400M Hurdles Men"
## [17] "400M Hurdles Women" "400M Men"
## [19] "400M Women" "4X100M Relay Men"
## [21] "4X100M Relay Women" "4X400M Relay Men"
## [23] "4X400M Relay Women" "5000M Men"
## [25] "5000M Women" "50Km Race Walk Men"
## [27] "800M Men" "800M Women"
## [29] "Decathlon Men" "Discus Throw Men"
## [31] "Discus Throw Women" "Hammer Throw Men"
## [33] "Hammer Throw Women" "Heptathlon Women"
## [35] "High Jump Men" "High Jump Women"
## [37] "Javelin Throw Men" "Javelin Throw Women"
## [39] "Long Jump Men" "Long Jump Women"
## [41] "Marathon Men" "Marathon Women"
## [43] "Pole Vault Men" "Pole Vault Women"
## [45] "Shot Put Men" "Shot Put Women"
## [47] "Triple Jump Men" "Triple Jump Women"
men100 <- olm[olm[,2]=="100M Men",]# 男子100m走を取り出す。
dim(men100)
## [1] 82 8
str(men100)# データの構造を見る
## 'data.frame': 82 obs. of 8 variables:
## $ Gender : Factor w/ 5 levels "-0.9","+0.1",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Event : Factor w/ 48 levels "","10000M Men",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ Location : Factor w/ 24 levels "","Amsterdam",..: 18 18 18 7 7 7 23 23 23 6 ...
## $ Year : int 2016 2016 2016 2008 2008 2008 2000 2000 2000 1992 ...
## $ Medal : Factor w/ 4 levels "","B","G","S": 3 4 2 3 4 2 3 4 2 3 ...
## $ Name : Factor w/ 1684 levels "Fritz Erik ELMSテδ\x83?TER",..: 1529 857 82 1529 1301 1589 1070 142 1193 981 ...
## $ Nationality: Factor w/ 106 levels "","20.67","52.27",..: 59 103 24 59 97 103 103 97 15 46 ...
## $ Result : Factor w/ 1940 levels "","0:42.2","0:42.4",..: 1909 1914 1916 1906 1914 1916 1913 1918 218 1917 ...
# min(men100[,8])# 最速のデータを知る
(c8 <- as.numeric(as.character(men100[,8]))) # 文字列型を経て数値型に変換する
## [1] 9.81 9.89 9.91 9.69 9.89 9.91 9.87 9.99 10.04 9.96 10.02 10.04
## [13] 9.99 10.19 10.22 10.06 10.08 10.14 9.90 10.00 10.00 10.20 10.20 10.30
## [25] 10.40 10.40 10.40 10.30 10.40 10.50 10.80 10.90 10.90 10.80 10.80 10.90
## [37] 10.80 NA NA 11.00 11.10 11.20 9.63 9.75 9.79 9.85 9.86 9.87
## [49] 9.84 9.89 9.90 10.25 10.25 10.39 10.14 10.24 10.33 10.00 10.20 10.20
## [61] 10.50 10.50 10.60 10.30 10.40 10.60 10.30 10.30 10.40 10.60 10.70 10.80
## [73] 10.80 10.90 10.90 11.00 11.20 11.20 12.00 12.20 12.60 12.60
min(c8, na.rm=T)# 最速のデータを知る
## [1] 9.63
men100[,8] <- -c8# 極値統計は最大値を求める方法なので最小値は逆転させる
men1002 <- men100[order(men100$Year),]# 年代順に並べる
men100G <- men1002[men1002[,5]=="G",]# 金メダルだけ取り出す
men100S <- men1002[men1002[,5]=="S",]# 銀メダルだけ取り出す
men100B <- men1002[men1002[,5]=="B",]# 銅メダルだけ取り出す
dim(men100G); dim(men100S); dim(men100B)# データサイズの確認をする
## [1] 27 8
## [1] 27 8
## [1] 28 8
men100B # 銅だけ多いのでどこか調べる
## Gender Event Location Year Medal Name
## 150 M 100M Men Athens 1896 B Alajos SZOKOLYI
## 151 M 100M Men Athens 1896 B Francis LANE
## 111 M 100M Men Paris 1900 B Stanley ROWLEY
## 147 M 100M Men St Louis 1904 B William HOGENSON
## 108 M 100M Men London 1908 B Robert KERR
## 144 M 100M Men Stockholm 1912 B Donald LIPPINCOTT
## 105 M 100M Men Antwerp 1920 B Harry EDWARD
## 141 M 100M Men Paris 1924 B Arthur PORRITT
## 102 M 100M Men Amsterdam 1928 B Georg LAMMERS
## 138 M 100M Men Los Angeles 1932 B Arthur JONATH
## 99 M 100M Men Berlin 1936 B Martinus OSENDARP
## 135 M 100M Men London 1948 B Lloyd LABEACH
## 96 M 100M Men Helsinki 1952 B Emmanuel McDONALD BAILEY
## 132 M 100M Men Melbourne / Stockholm 1956 B Hector HOGAN
## 93 M 100M Men Rome 1960 B Peter RADFORD
## 129 M 100M Men Tokyo 1964 B Harry JEROME
## 90 M 100M Men Mexico 1968 B Charlie GREENE
## 126 M 100M Men Munich 1972 B Lennox MILLER
## 87 M 100M Men Montreal 1976 B Valery BORZOV
## 123 M 100M Men Moscow 1980 B Peter PETROV
## 84 M 100M Men Los Angeles 1984 B Ben JOHNSON
## 81 M 100M Men Barcelona 1992 B Dennis MITCHELL
## 120 M 100M Men Atlanta 1996 B Ato BOLDON
## 78 M 100M Men Sydney 2000 B Obadele THOMPSON
## 117 M 100M Men Athens 2004 B Maurice GREENE
## 75 M 100M Men Beijing 2008 B Walter DIX
## 114 M 100M Men London 2012 B Justin GATLIN
## 72 M 100M Men Rio 2016 B Andre DE GRASSE
## Nationality Result
## 150 HUN -12.60
## 151 USA -12.60
## 111 AUS -11.20
## 147 USA -11.20
## 108 CAN NA
## 144 USA -10.90
## 105 GBR -10.90
## 141 NZL -10.80
## 102 GER -10.90
## 138 GER -10.40
## 99 NED -10.50
## 135 PAN -10.60
## 96 GBR -10.40
## 132 AUS -10.60
## 93 GBR -10.30
## 129 CAN -10.20
## 90 USA -10.00
## 126 JAM -10.33
## 87 URS -10.14
## 123 BUL -10.39
## 84 CAN -10.22
## 81 USA -10.04
## 120 TTO -9.90
## 78 BAR -10.04
## 117 USA -9.87
## 75 USA -9.91
## 114 USA -9.79
## 72 CAN -9.91
# 1896年アテネで2人銅と分かる
(men100mat <- cbind(men100G[,c(4,6,8)], men100S[,c(6,8)], men100B[2:28,c(6,8)]))
## Year Name Result Name Result
## 148 1896 Thomas BURKE -12.00 Fritz HOFMANN -12.20
## 109 1900 Frank JARVIS -11.00 Walter TEWKSBURY -11.10
## 145 1904 Archie HAHN -11.00 Nate CARTMELL -11.20
## 106 1908 Reggie WALKER -10.80 John RECTOR NA
## 142 1912 Ralph CRAIG -10.80 Alvah MEYER -10.90
## 103 1920 Charles PADDOCK -10.80 Morris KIRKSEY -10.80
## 139 1924 Harold ABRAHAMS -10.60 Jackson SCHOLZ -10.70
## 100 1928 Percy WILLIAMS -10.80 Jack LONDON -10.90
## 136 1932 Eddie TOLAN -10.30 Ralph METCALFE -10.30
## 97 1936 Jesse OWENS -10.30 Ralph METCALFE -10.40
## 133 1948 Harrison DILLARD -10.30 Barney EWELL -10.40
## 94 1952 Lindy REMIGINO -10.40 Herbert MCKENLEY -10.40
## 130 1956 Bobby MORROW -10.50 Thane BAKER -10.50
## 91 1960 Armin HARY -10.20 Dave SIME -10.20
## 127 1964 Bob HAYES -10.00 Enrique FIGUEROLA -10.20
## 88 1968 Jim HINES -9.90 Lennox MILLER -10.00
## 124 1972 Valery BORZOV -10.14 Robert TAYLOR -10.24
## 85 1976 Hasely CRAWFORD -10.06 Donald QUARRIE -10.08
## 121 1980 Allan WELLS -10.25 Silvio LEONARD SARRIA -10.25
## 82 1984 Carl LEWIS -9.99 Sam GRADDY -10.19
## 79 1992 Linford CHRISTIE -9.96 Frank FREDERICKS -10.02
## 118 1996 Donovan BAILEY -9.84 Frank FREDERICKS -9.89
## 76 2000 Maurice GREENE -9.87 Ato BOLDON -9.99
## 115 2004 Justin GATLIN -9.85 Francis OBIKWELU -9.86
## 73 2008 Usain BOLT -9.69 Richard THOMPSON -9.89
## 112 2012 Usain BOLT -9.63 Yohan BLAKE -9.75
## 70 2016 Usain BOLT -9.81 Justin GATLIN -9.89
## Name Result
## 148 Francis LANE -12.60
## 109 Stanley ROWLEY -11.20
## 145 William HOGENSON -11.20
## 106 Robert KERR NA
## 142 Donald LIPPINCOTT -10.90
## 103 Harry EDWARD -10.90
## 139 Arthur PORRITT -10.80
## 100 Georg LAMMERS -10.90
## 136 Arthur JONATH -10.40
## 97 Martinus OSENDARP -10.50
## 133 Lloyd LABEACH -10.60
## 94 Emmanuel McDONALD BAILEY -10.40
## 130 Hector HOGAN -10.60
## 91 Peter RADFORD -10.30
## 127 Harry JEROME -10.20
## 88 Charlie GREENE -10.00
## 124 Lennox MILLER -10.33
## 85 Valery BORZOV -10.14
## 121 Peter PETROV -10.39
## 82 Ben JOHNSON -10.22
## 79 Dennis MITCHELL -10.04
## 118 Ato BOLDON -9.90
## 76 Obadele THOMPSON -10.04
## 115 Maurice GREENE -9.87
## 73 Walter DIX -9.91
## 112 Justin GATLIN -9.79
## 70 Andre DE GRASSE -9.91
# 各行が各年に当たり、金銀銅について氏名と記録を並べた行列をつくる 銅を一つ外す
men100mat2 <- men100mat[,c(1,3,5,7)] # 年代と金銀銅記録のみの行列をつくる
plot(men100mat[,1], men100mat[,3])# 年代と金メダル記録の散布図をつくる

# 明らかにトレンドがある
library(evd)
# fgev(men100mat2[,2])# 定常モデル
(men100gev <- fgev(men100mat2[,2], method="Nelder-Mead"))# ネルダー・ミード法
##
## Call: fgev(x = men100mat2[, 2], method = "Nelder-Mead")
## Deviance: 32.03137
##
## Estimates
## loc scale shape
## -10.3912 0.5665 -0.7259
##
## Standard Errors
## loc scale shape
## 0.1173 0.1022 0.1425
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 134
oldpar <- par()
par(mfrow=c(2,2))
plot(men100gev)

fgev(men100mat2[,2], method="Nelder-Mead", prob=0)
##
## Call: fgev(x = men100mat2[, 2], prob = 0, method = "Nelder-Mead")
## Deviance: 32.03143
##
## Estimates
## quantile scale shape
## -9.6109 0.5660 -0.7257
##
## Standard Errors
## quantile scale shape
## 0.03848 0.10181 0.14128
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 126
(men100gev.ns <- fgev(men100mat2[,2], nsloc=men100mat2[,1]))# 非定常モデル
##
## Call: fgev(x = men100mat2[, 2], nsloc = men100mat2[, 1])
## Deviance: 75.13933
##
## Estimates
## loc loctrend scale shape
## -1.056e+01 -1.443e-04 4.106e-01 -9.828e-09
##
## Standard Errors
## loc loctrend scale shape
## 9.846e-02 3.378e-07 3.261e-02 2.236e-02
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 101
## Gradient Evaluations: 3
plot(men100gev.ns)# 当てはまりが悪い

(men100gev2.ns <- fgev(men100mat2[,2],method="Nelder-Mead", nsloc=men100mat2[,1]))
##
## Call: fgev(x = men100mat2[, 2], nsloc = men100mat2[, 1], method = "Nelder-Mead")
## Deviance: 30.52477
##
## Estimates
## loc loctrend scale shape
## -1.109e+01 3.695e-04 5.523e-01 -7.520e-01
##
## Standard Errors
## loc loctrend scale shape
## 1.134e-01 2.001e-06 1.002e-01 1.484e-01
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 297
plot(men100gev2.ns)

(men100gev3.ns <- fgev(men100mat2[,2], nsloc=men100mat2[,1]/100)) # 別の非定常
##
## Call: fgev(x = men100mat2[, 2], nsloc = men100mat2[, 1]/100)
## Deviance: 12.11418
##
## Estimates
## loc loctrend scale shape
## -18.7313 0.4295 0.4234 -0.8884
##
## Standard Errors
## loc loctrend scale shape
## 0.091992 0.000002 0.159343 0.431660
##
## Optimization Information
## Convergence: successful
## Function Evaluations: 83
## Gradient Evaluations: 13
plot(men100gev3.ns)# 当てはまりがよい

par(oldpar)
library(ismev)
men100gev.is <- gev.fit(men100mat2[,2])# 定常モデル
## $conv
## [1] 0
##
## $nllh
## [1] 16.01569
##
## $mle
## [1] -10.3912566 0.5666167 -0.7259960
##
## $se
## [1] 0.1173063 0.1022269 0.1424777
gev.diag(men100gev.is)# 当てはまりがよい

# ismevで非定常モデルを作る
time1 <- (men100mat2[,1]-1896)/120# 1896年から2016年を0-1とする
time=cbind(time1, time1^2)
men100t1 <- gev.fit(men100mat2[,2], ydat=time, mul=1) # 位置母数について1次式
## $model
## $model[[1]]
## [1] 1
##
## $model[[2]]
## NULL
##
## $model[[3]]
## NULL
##
##
## $link
## [1] "c(identity, identity, identity)"
##
## $conv
## [1] 0
##
## $nllh
## [1] -6.133671
##
## $mle
## [1] -11.0279238 1.3125189 0.2471653 -0.6985406
##
## $se
## [1] 0.07373383 0.10262283 0.04206908 0.13535455
gev.diag(men100t1)

men100t2 <- gev.fit(men100mat2[,2], ydat=time, mul=c(1,2))
## $model
## $model[[1]]
## [1] 1 2
##
## $model[[2]]
## NULL
##
## $model[[3]]
## NULL
##
##
## $link
## [1] "c(identity, identity, identity)"
##
## $conv
## [1] 0
##
## $nllh
## [1] -11.62056
##
## $mle
## [1] -11.2241918 2.6055195 -1.2584368 0.2172678 -0.8388815
##
## $se
## [1] 0.07888038 0.28575650 0.23638724 0.03613637 0.16742839
gev.diag(men100t2)

men100gev.is$nllh*2+6 # 定常モデルのAIC
## [1] 38.03137
men100t1$nllh*2+8 # 1次モデルのAIC
## [1] -4.267342
men100t2$nllh*2+10 # 2次モデルのAIC
## [1] -13.24112
men100t2$mle # 2次モデルの推定値
## [1] -11.2241918 2.6055195 -1.2584368 0.2172678 -0.8388815
t1 <- (2020-1896)/120
temp <- men100t2$mle # 何度も書く手間を省く
m100qgev <- function(year, p){
t1 <- (year-1896)/120
qgev(p, loc=temp[1]+temp[2]*t1+temp[3]*t1^2, scale=temp[4],
shape=temp[5])}
m100qgev(2020,c(0.5,0.9,0.99))
## [1] -9.806999 -9.655769 -9.622018
plot(2000+4*(1:10), m100qgev(2000+4*(1:10), 0.9))

men100r <- rlarg.fit(men100mat2[,2:4])
## $conv
## [1] 0
##
## $nllh
## [1] -24.25423
##
## $mle
## [1] -10.0469241 0.3063958 -0.7199866
##
## $se
## [1] 0.05233495 0.02515602 0.06870977
# rlarg.diag(men100r)
men100rt1 <- rlarg.fit(men100mat2[,2:4], ydat=time, mul=1)
## $model
## $model[[1]]
## [1] 1
##
## $model[[2]]
## NULL
##
## $model[[3]]
## NULL
##
##
## $link
## [1] "c(identity, identity, identity)"
##
## $conv
## [1] 0
##
## $nllh
## [1] -72.20731
##
## $mle
## [1] -10.9044558 1.3038146 0.1656709 -0.7247666
##
## $se
## [1] 0.04236019 0.06315224 0.01378331 0.06909201
# rlarg.diag(men100rt1)
men100rt2 <- rlarg.fit(men100mat2[,2:4], ydat=time, mul=2)
## $model
## $model[[1]]
## [1] 2
##
## $model[[2]]
## NULL
##
## $model[[3]]
## NULL
##
##
## $link
## [1] "c(identity, identity, identity)"
##
## $conv
## [1] 0
##
## $nllh
## [1] -57.13127
##
## $mle
## [1] -10.6077131 1.1631728 0.2047979 -0.6932878
##
## $se
## [1] 0.06812956 0.16342240 0.01596855 0.06216688
# rlarg.diag(men100rt2)
men100mat3 <- abs(100/men100mat2) # 100/行列の各成分の行列をつくる
men100mat3[,1] <- men100mat2[,1] # 1列目の年号だけ書き換え
men100mat3
## Year Result Result.1 Result.2
## 148 1896 8.333333 8.196721 7.936508
## 109 1900 9.090909 9.009009 8.928571
## 145 1904 9.090909 8.928571 8.928571
## 106 1908 9.259259 NA NA
## 142 1912 9.259259 9.174312 9.174312
## 103 1920 9.259259 9.259259 9.174312
## 139 1924 9.433962 9.345794 9.259259
## 100 1928 9.259259 9.174312 9.174312
## 136 1932 9.708738 9.708738 9.615385
## 97 1936 9.708738 9.615385 9.523810
## 133 1948 9.708738 9.615385 9.433962
## 94 1952 9.615385 9.615385 9.615385
## 130 1956 9.523810 9.523810 9.433962
## 91 1960 9.803922 9.803922 9.708738
## 127 1964 10.000000 9.803922 9.803922
## 88 1968 10.101010 10.000000 10.000000
## 124 1972 9.861933 9.765625 9.680542
## 85 1976 9.940358 9.920635 9.861933
## 121 1980 9.756098 9.756098 9.624639
## 82 1984 10.010010 9.813543 9.784736
## 79 1992 10.040161 9.980040 9.960159
## 118 1996 10.162602 10.111223 10.101010
## 76 2000 10.131712 10.010010 9.960159
## 115 2004 10.152284 10.141988 10.131712
## 73 2008 10.319917 10.111223 10.090817
## 112 2012 10.384216 10.256410 10.214505
## 70 2016 10.193680 10.111223 10.090817
sp.gev <- gev.fit(men100mat3[,2])
## $conv
## [1] 0
##
## $nllh
## [1] 14.43074
##
## $mle
## [1] 9.6306356 0.5101265 -0.6486014
##
## $se
## [1] 0.10609851 0.08827584 0.13796928
gev.diag(sp.gev)

sp.gev.t1 <- gev.fit(men100mat3[,2], ydat=time, mul=1)
## $model
## $model[[1]]
## [1] 1
##
## $model[[2]]
## NULL
##
## $model[[3]]
## NULL
##
##
## $link
## [1] "c(identity, identity, identity)"
##
## $conv
## [1] 0
##
## $nllh
## [1] -10.05206
##
## $mle
## [1] 9.0198272 1.2780681 0.2022104 -0.6010034
##
## $se
## [1] 0.06994963 0.10728719 0.03234943 0.12565232
gev.diag(sp.gev.t1)

sp.gev.t2 <- gev.fit(men100mat3[,2], ydat=time, mul=2)
## $model
## $model[[1]]
## [1] 2
##
## $model[[2]]
## NULL
##
## $model[[3]]
## NULL
##
##
## $link
## [1] "c(identity, identity, identity)"
##
## $conv
## [1] 0
##
## $nllh
## [1] -2.053732
##
## $mle
## [1] 9.2648284 1.1095349 0.2717253 -0.5902236
##
## $se
## [1] 0.08695174 0.18501438 0.04244043 0.11048359
gev.diag(sp.gev.t2)
