Skip to content
GitLab
Select Git revision
  • master default
Find file BlameHistoryPermalink
main_fevd_bayesian.R 2.75 KiB
1 
# Title     : TODO
2 
# Objective : TODO
3 
# Created by: erwan
4 
# Created on: 04/10/2019
5 
library(extRemes)
6 
# library(data.table)
7 
# library(stats4)
8 
library(SpatialExtremes)
9 
source('fevd_fixed.R')
10 
# Sample from a GEV
11 
set.seed(42)
12 
N <- 100
13 
loc = 0; scale = 1; shape <- 0.1
14 
x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
15
16 
# fevdPriorMy <- function (theta, q, p, log = FALSE){
17 
#     x = theta["shape"] + 0.5
18 
#
19 
#     print(theta)
20 
#     print(theta["location"])
21 
#     print(dunif(theta["location"]))
22 
#     print(theta[0])
23 
#     dfun <- function(th) dbeta(th[1], shape1 = th[2], shape2 = th[3],
24 
#         log = log)
25 
#     th <- cbind(theta, q, p)
26 
#     res <- apply(th, 1, dfun)
27 
#     return(prod(res))
28 
# }
29
30
31
32 
# print(pbeta(1.0, 1, 1))
33 
# print(pbeta(0.5, 1, 1))
34 
# print(fevdPriorCustom(2.0, 0.0, 0.0))
35
36
37
38 
# res = fevd(x_gev, method='Bayesian', priorFun="fevdPriorMyMy", priorParams=list(q=c(6), p=c(9)), iter=5000, verbose=TRUE, use.phi=FALSE)
39 
# res = fevd(x_gev, method='GMLE', iter=5000, verbose=TRUE, use.phi=FALSE)
40
41 
# Without covariate
42 
# res = fevd_fixed(x_gev, method='Bayesian', priorFun="fevdPriorCustom", priorParams=list(q=c(6), p=c(9)), iter=5000, verbose=TRUE, use.phi=FALSE)
43
44 
# Add covariate
45 
coord <- matrix(ncol=1, nrow = N)
46 
coord[,1]=seq(1,N,1)
47 
colnames(coord) = c("T")
48 
coord = data.frame(coord, stringsAsFactors = TRUE)
49 
# res = fevd_fixed(x_gev, data=coord, location.fun= ~T, scale.fun= ~T, method='Bayesian', priorFun="fevdPriorCustom", priorParams=list(q=c(6), p=c(9)), iter=5000, verbose=TRUE, use.phi=FALSE)
50 
res = fevd_fixed(x_gev, data=coord, method='Bayesian', priorFun="fevdPriorCustom", priorParams=list(q=c(6), p=c(9)), iter=5000, verbose=TRUE, use.phi=FALSE)
51
52 
# Some display for the results
53 
# print(res)
54 
# print('here')
55 
# print(res$constant.loc)
56 
# print('here2')
57 
# print(res$method)
58 
# print(res$priorFun)
59 
# print(res$priorParams)
60 
# m = res$results
61 
# print(class(res$chain.info))
62 
# print(dim(m))
63 
# print(m[1,])
64 
# print(m[1,1])
65 
# print(res$chain.info[1,])
66
67
68
69 
# Confidence interval
70 
res_ci = ci(res, alpha = 0.05, type = c("return.level"),
7172737475767778798081828384858687888990919293949596979899100101
return.period = 50, FUN = "mean", tscale = FALSE) print(res_ci) # print(attributes(res_ci)) # print(dimnames(res_ci)) # print(res_ci[1]) # covariates = c(1.0, 100.0, 1000.0) # v = make.qcov(res, vals = list(mu1 = covariates, sigma1 = covariates)) # v = make.qcov(res, vals = list(mu1 = c(0.0), sigma1 = c(0.0))) # res_find = findpars(res, FUN = "mean" , tscale = FALSE, qcov = v, # burn.in=4998) # print(res_find) # # res_ci = ci(res, alpha = 0.05, type = c("return.level"), # burn.in=4998, # return.period = 50, FUN = "mean", tscale = FALSE, qcov=v) # print(res_ci) # # r = qgev(0.98, 0.1427229 , 1.120554, -0.1008511) # print(r) # print(attributes(res_ci)) # print(dimnames(res_ci)) # print(res_ci[1]) # print(attr(res_ci), '')