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