diff --git a/extreme_estimator/margin_fits/gev/fevd.R b/extreme_estimator/margin_fits/gev/fevd.R
new file mode 100644
index 0000000000000000000000000000000000000000..d1018f1725acb7ec7b3bd85bcd0e8569daba5859
--- /dev/null
+++ b/extreme_estimator/margin_fits/gev/fevd.R
@@ -0,0 +1,99 @@
+# 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
diff --git a/extreme_estimator/margin_fits/gev/fevd.py b/extreme_estimator/margin_fits/gev/fevd.py
new file mode 100644
index 0000000000000000000000000000000000000000..8089836e856da6effbc63d77e2389671b8b03001
--- /dev/null
+++ b/extreme_estimator/margin_fits/gev/fevd.py
@@ -0,0 +1,23 @@
+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
diff --git a/extreme_estimator/margin_fits/margin_fits_utils.py b/extreme_estimator/margin_fits/margin_fits_utils.py
index 31de2d174bcea68bf6dac5dd5e3ad48821fe933e..bc41ff956447f852acf85242a95958a7b49f0a05 100644
--- a/extreme_estimator/margin_fits/margin_fits_utils.py
+++ b/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()
diff --git a/test/test_extreme_estimator/test_margin_fits/test_gev/test_gevmle_fit.py b/test/test_extreme_estimator/test_margin_fits/test_gev/test_gevmle_fit.py
index 82616977e72575c9c43f5e859cdbc0de8207ea57..15d003e45c3e2a3588f866f852422e37c501d4d8 100644
--- a/test/test_extreme_estimator/test_margin_fits/test_gev/test_gevmle_fit.py
+++ b/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