r - 为 R 中的年降水量数据拟合 GEV

Question

我正在使用年降水量最大值，并希望使用 GEV 来拟合我的数据。虽然我能够做到这一点，但我不确定 R 中的 GEV 分析是否检测到我如何分配我的块最大值。理想情况下，我想使用块最大值方法将我的 14 年时间序列分成两个 7 年块（即第 1 年到第 7 年，然后是第 8 年到第 14 年）。我也有 9 个数据集，所以我汇总了 9 个数据集的最大值。例如，我通过对第 1 年到第 14 年的所有年度最大值进行分组来汇总数据，从而产生 14 年的 126 个年度最大值：

library(extRemes)

#Obtaining annual maxima from the 9 datasets

Gmax <- sapply(unstack(subset), function(r){max(values(r))}) 
Gmax1 <- sapply(unstack(subset5), function(r){max(values(r))})
Gmax2 <- sapply(unstack(subset8), function(r){max(values(r))})
Gmax3 <- sapply(unstack(subset11), function(r){max(values(r))})
Gmax4 <- sapply(unstack(subset15), function(r){max(values(r))})
Gmax5 <- sapply(unstack(subset18), function(r){max(values(r))})
Gmax6 <- sapply(unstack(subset21), function(r){max(values(r))})
Gmax7 <- sapply(unstack(subset25), function(r){max(values(r))})
Gmax8 <- sapply(unstack(subset27), function(r){max(values(r))})

#Indexing to select all of the annual maxima for the 14 years for the 9 datasets (126 maxima)

Gcomb <- c(Gmax[1], Gmax1[1], Gmax2[1], Gmax3[1], Gmax4[1], Gmax5[1], Gmax6[1], Gmax7[1], Gmax8[1],
Gmax[2], Gmax1[2], Gmax2[2], Gmax3[2], Gmax4[2], Gmax5[2], Gmax6[2], Gmax7[2], Gmax8[2], Gmax[3],  
Gmax1[3], Gmax2[3], Gmax3[3], Gmax4[3], Gmax5[3], Gmax6[3], Gmax7[3], Gmax8[3], Gmax[4], Gmax1[4], 
Gmax2[4], Gmax3[4], Gmax4[4], Gmax5[4], Gmax6[4], Gmax7[4], Gmax8[4], Gmax[5], Gmax1[5], Gmax2[5], 
Gmax3[5], Gmax4[5], Gmax5[5], Gmax6[5], Gmax7[5], Gmax8[5], Gmax[6], Gmax1[6], Gmax2[6], Gmax3[6], 
Gmax4[6], Gmax5[6], Gmax6[6], Gmax7[6], Gmax8[6], Gmax[7], Gmax1[7], Gmax2[7], Gmax3[7], Gmax4[7], 
Gmax5[7], Gmax6[7], Gmax7[7], Gmax8[7], Gcomb1<-c(Gmax[8], Gmax1[8], Gmax2[8], Gmax3[8], Gmax4[8], 
Gmax5[8], Gmax6[8], Gmax7[8], Gmax8[8], Gmax[9], Gmax1[9], Gmax2[9], Gmax3[9], Gmax4[9], Gmax5[9], 
Gmax6[9], Gmax7[9], Gmax8[9], Gmax[10], Gmax1[10], Gmax2[10], Gmax3[10], Gmax4[10], Gmax5[10], 
Gmax6[10], Gmax7[10], Gmax8[10], Gmax[11], Gmax1[11], Gmax2[11], Gmax3[11], Gmax4[11], Gmax5[11], 
Gmax6[11], Gmax7[11], Gmax8[11], Gmax[12], Gmax1[12], Gmax2[12], Gmax3[12], Gmax4[12], Gmax5[12], 
Gmax6[12], Gmax7[12], Gmax8[12], Gmax[13], Gmax1[13], Gmax2[13], Gmax3[13], Gmax4[13], Gmax5[13], 
Gmax6[13], Gmax7[13], Gmax8[13], Gmax[14], Gmax1[14], Gmax2[14], Gmax3[14], Gmax4[14], Gmax5[14], 
Gmax6[14], Gmax7[14], Gmax8[14]))

> Gcomb

[1] 285.9906 159.0246 216.1329 317.3168 134.0706 240.5591 143.3593 278.6423 300.2213 287.6668 
133.6294 235.7827 348.0946 133.2577
[15] 200.7394 142.7288 255.3043 332.2073 276.2251 153.0331 218.5597 320.7929 123.5906 215.9244 
165.5872 386.4905 259.4115 247.1718
[29] 160.4190 276.9192 313.5304 175.7406 211.9053 132.6234 318.1059 297.1881 254.9781 148.2776 
270.5989 292.8830 155.7216 194.3641
[43] 126.4537 236.9236 383.5452 308.1921 142.6571 244.4956 356.3880 160.6649 243.4610 154.6999 
244.6703 349.7785 234.6858 172.3195
[57] 275.7311 359.1511 115.8130 185.3560 180.2979 210.3161 302.0780 235.7213 157.7454 262.2424 
307.9303 136.0832 224.9627 146.9899
[71] 297.2613 337.8194 242.0865 148.3942 205.8541 342.4027 153.2334 191.5551 137.0983 367.9224     
428.6493 247.1022 165.9219 247.4879
[85] 327.0393 146.6705 202.0553 146.3462 254.7741 346.1007 233.6528 139.6812 317.0052 336.3853 
155.3950 219.7212 144.6838 284.4246
[99] 306.7610 243.9972 140.4338 267.6197 366.8654 159.3714 182.8306 136.9110 284.0912 362.1337 
297.6855 161.9404 395.7047 402.7022
[113] 128.2466 218.4132 135.6528 320.8012 308.9064 264.3488 136.1558 285.4910 324.7307 177.1121 
226.3639 138.8371 280.7834 370.2415

#Fitting a GEV

superfit <- fevd(x=Gcomb, type="GEV")

fevd(x = Gcomb, type = "GEV")

[1] "Estimation Method used: MLE"


Negative Log-Likelihood Value:  727.0276 


Estimated parameters:
location       scale       shape 
206.1101087  72.3405803  -0.1773972 

Standard Error Estimates:
 location     scale     shape 
7.8269021 6.0762930 0.1016432 

Estimated parameter covariance matrix.
         location      scale       shape
location 61.2603966 16.5091185 -0.42074406
scale    16.5091185 36.9213365 -0.41062737
shape    -0.4207441 -0.4106274  0.01033135

AIC = 1460.055 

BIC = 1468.564 

此输出基于我的整个时间序列。我的问题是 R 怎么知道我正在制作 7 年的区块（即第 1 年到第 7 年，以及第 8 年到第 14 年）？

谢谢，任何反馈将不胜感激！

标签： r