From a47f98657ae951205949c8c968b772f4d013206c Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Fri, 4 Oct 2019 17:04:16 +0200
Subject: [PATCH] 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 | 99 +++++++++++++++++++
extreme_estimator/margin_fits/gev/fevd.py | 23 +++++
.../margin_fits/margin_fits_utils.py | 3 -
.../test_gev/test_gevmle_fit.py | 2 +
4 files changed, 124 insertions(+), 3 deletions(-)
create mode 100644 extreme_estimator/margin_fits/gev/fevd.R
create mode 100644 extreme_estimator/margin_fits/gev/fevd.py
diff --git a/extreme_estimator/margin_fits/gev/fevd.R b/extreme_estimator/margin_fits/gev/fevd.R
new file mode 100644
index 00000000..d1018f17
--- /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 00000000..8089836e
--- /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 31de2d17..bc41ff95 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 82616977..15d003e4 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
--
GitLab