refactor a bit. add fevd files (work in progress) to use the GMLE method in a Bayesian mode.

extreme_estimator/margin_fits/gev/fevd.R

+# Title     : TODO
+# Objective : TODO
+# Created by: erwan
+# Created on: 04/10/2019
+library(extRemes)
+library(stats4)
+library(SpatialExtremes)
+
+set.seed(42)
+N <- 1000
+loc = 0; scale = 1; shape <- 0.1
+x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
+start_loc = 0; start_scale = 1; start_shape = 1
+# res = fevd(x_gev, method='GMLE')
+
+fevdPriorMy <- function (theta, q, p, log = FALSE){
+    x = theta["shape"] + 0.5
+
+    print(theta)
+    print(theta["location"])
+    print(dunif(theta["location"]))
+    print(theta[0])
+    dfun <- function(th) dbeta(th[1], shape1 = th[2], shape2 = th[3],
+        log = log)
+    th <- cbind(theta, q, p)
+    res <- apply(th, 1, dfun)
+    return(prod(res))
+}
+
+fevdPriorMyMy <- function (theta, q, p, log = FALSE){
+    print(theta)
+    print(q)
+    print(p)
+    x = theta[length(theta)]
+    # + 0.5 enables to shift the Beta law in the interval [-0.5, 0.5]
+    res = dbeta(x + 0.5, q, p, log = TRUE)
+    return(res)
+}
+
+
+print(pbeta(1.0, 1, 1))
+print(pbeta(0.5, 1, 1))
+print(fevdPriorMy(2.0, 0.0, 0.0))
+
+res = fevd(x_gev, method='Bayesian', priorFun="fevdPriorMyMy", priorParams=list(q=c(6), p=c(9)), iter=5000)
+print(res)
+# res = fevd(x_gev, method='Bayesian')
+# print(res)
+
+priorFun="shapePriorBeta"
+shapePriorBeta
+#
+#
+# print(shapePriorBeta(0.0, 6, 9))
+# priorParams=list(q=c(1, 1, 6), p=c(1, 1, 9))
+# p.i <- do.call(priorFun, c(list(1.0), priorParams))
+# print(p.i)
+
+# priorFun <- "shapePriorBeta"
+# priorParams <- list(q = 6, p = 9)
+# priorFun <- "fevdPriorDefault"
+# priorParams <- list(q = 6, p = 9)
+# e = do.call(priorFun, c(list(0.0), priorParams))
+# print(e)
+#
+# print(res$method)
+# print(res$priorFun)
+# print(res$priorParams)
+# m = res$results
+# print(m[2,1])
+# print(class(res$chain.info))
+# print(res$chain.info[[1]])
+# # summary(res)
+# print(attributes(res))
+# print('here')
+# print(attr(res, 'chain.info'))
+# print(attr(res, "method"))
+# print(attr(res, "x"))
+# print(attr(res, "priorParams"))
+
+# print(res.method)
+
+
+# p.i <- do.call(shapePriorBeta, c(list(theta =  c(-0.12572432087762, -0.0567634605386987, 0.133782230298093)), priorParams=list(q = 6, p = 9)))
+# print(p.i)
+# a = fevd(x_gev, method='Bayesian', priorFun="shapePriorBeta", priorParams=list(q = 6, p = 9))
+
+# priorParams=list(v=c(0.1, 10, 0.1)),
+#     initial=list(location=0, scale=0.1, shape=-0.5)),
+
+# print(a)
+#
+# # S3 method for fevd.bayesian
+# summary(a, FUN = "mean", burn.in = 499)
+
+# print(a.results)
+
+# Bayesian method is using a normal distribution functions for the shape parameter
+# GMLE distribution is using a Beta distribution for the shape parameter

extreme_estimator/margin_fits/gev/fevd.py

+import numpy as np
+
+from extreme_estimator.margin_fits.gev.gev_fit import GevFit
+from extreme_estimator.margin_fits.gev.gev_params import GevParams
+
+
+class IsmevGevFit(GevFit):
+
+    def __init__(self, x_gev: np.ndarray, y=None, mul=None):
+        super().__init__(x_gev)
+        self.y = y
+        self.mul = mul
+        self.res = fevd_gev_fit(x_gev, y, mul)
+
+    # @property
+    # def gev_params(self) -> GevParams:
+    #     assert self.y is None
+    #     gev_params_dict = dict(zip(GevParams.PARAM_NAMES, self.res['mle']))
+    #     return GevParams.from_dict(gev_params_dict)
+
+
+def fevd_gev_fit(x, y, mul):
+    pass
\ No newline at end of file

extreme_estimator/margin_fits/margin_fits_utils.py

@@ -34,9 +34,6 @@ def ismev_gev_fit(x_gev, y, mul):
     For non-stationary fitting it is recommended that the covariates within the generalized linear models are
     (at least approximately) centered and scaled (i.e.the columns of ydat should be approximately centered and scaled).
     """
-    # print('Mean x={}, variance x={}'.format(np.mean(x_gev), np.var(x_gev)))
-    # print('Mean y={}, variance y={}'.format(np.mean(y), np.var(y)))
-
     xdat = ro.FloatVector(x_gev)
     gev_fit = r('gev.fit')
     y = y if y is not None else get_null()

test/test_extreme_estimator/test_margin_fits/test_gev/test_gevmle_fit.py

@@ -29,6 +29,8 @@ class TestGevMleFit(unittest.TestCase):
         ismev_ref = {'loc': 0.0219, 'scale': 1.0347, 'shape': 0.8295}
         self.fit_estimator(estimator, ismev_ref)
 
+    # def test_s
+
     def fit_estimator(self, estimator, ref):
         # Compare the MLE estimated parameters to the reference
         mle_params_estimated = estimator.gev_params