From ac99ffcc57cae22530f159403d13b6ac22fb7f63 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 30 Oct 2018 18:34:26 +0100
Subject: [PATCH] [GEV] Add gev_fit in R, and corresponding wrapper in Python
---
R/gev_fit/gev_fit.R | 32 +++++++++++++++++++++++++++++++
R/gev_fit/gev_marginal.py | 38 +++++++++++++++++++++++++++++++++++++
R/gev_fit/gev_mle_fit.py | 32 +++++++++++++++++++++++++++++++
test/R/test_gev_fit.py | 40 +++++++++++++++++++++++++++++++++++++++
4 files changed, 142 insertions(+)
create mode 100644 R/gev_fit/gev_fit.R
create mode 100644 R/gev_fit/gev_marginal.py
create mode 100644 R/gev_fit/gev_mle_fit.py
create mode 100644 test/R/test_gev_fit.py
diff --git a/R/gev_fit/gev_fit.R b/R/gev_fit/gev_fit.R
new file mode 100644
index 00000000..124bd803
--- /dev/null
+++ b/R/gev_fit/gev_fit.R
@@ -0,0 +1,32 @@
+library(stats4)
+library(SpatialExtremes)
+
+# Boolean for python call
+call_main = !exists("python_wrapping")
+
+if (call_main) {
+ set.seed(42)
+ N <- 50
+ loc = 0; scale = 1; shape <- 1
+ x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
+ start_loc = 0; start_scale = 1; start_shape = 1
+}
+
+minus_log_likelihood_gev <- function(loc = 0, scale = 1, shape = 0){
+ R = suppressWarnings(dgev(x_gev, loc = loc, scale = scale, shape = shape))
+ -sum(log(R))
+}
+
+mle_gev <- function(loc, scale, shape){
+ mle(minuslogl = minus_log_likelihood_gev, start = list(loc = loc, scale = scale, shape = shape))
+}
+
+main <- function(){
+ res = mle_gev(start_loc, start_scale, start_shape)
+ print(attributes(res))
+ print(attr(res, 'coef'))
+}
+
+if (call_main) {
+ main()
+}
\ No newline at end of file
diff --git a/R/gev_fit/gev_marginal.py b/R/gev_fit/gev_marginal.py
new file mode 100644
index 00000000..7758532f
--- /dev/null
+++ b/R/gev_fit/gev_marginal.py
@@ -0,0 +1,38 @@
+from R.gev_fit.gev_mle_fit import GevMleFit
+
+
+def frechet_unitary_transformation(data, location, scale, shape):
+ """
+ Compute the unitary Frechet transformed data
+ (1 + \zeta \frac{z - \mu}{\sigma}) ^ {\frac{1}{\zeta}}
+ """
+ return (1 + shape * (data - location) / scale) ** (1 / shape)
+
+
+class GevMarginal(object):
+
+ def __init__(self, position, data, estimator=GevMleFit):
+ """
+ Class to fit a GEV a list of data. Compute also the corresponding unitary data
+
+ :param position: Represents the spatio-temporal position of the marginals
+ :param data: array of data corresponding to this position (and potentially its neighborhood)
+ """
+ self.position = position
+ self.data = data
+ self.gev_estimator = estimator(x_gev=data)
+ self.gev_parameters_estimated = [self.location, self.scale, self.shape]
+ self.frechet_unitary_data = frechet_unitary_transformation(data=data, location=self.location,
+ scale=self.scale, shape=self.shape)
+
+ @property
+ def location(self):
+ return self.gev_estimator.location
+
+ @property
+ def scale(self):
+ return self.gev_estimator.scale
+
+ @property
+ def shape(self):
+ return self.gev_estimator.shape
diff --git a/R/gev_fit/gev_mle_fit.py b/R/gev_fit/gev_mle_fit.py
new file mode 100644
index 00000000..b6f05736
--- /dev/null
+++ b/R/gev_fit/gev_mle_fit.py
@@ -0,0 +1,32 @@
+import rpy2.robjects as ro
+import numpy as np
+import rpy2.robjects.numpy2ri as rpyn
+
+
+# Defining some constants
+GEV_SCALE = 'scale'
+GEV_LOCATION = 'loc'
+GEV_SHAPE = 'shape'
+
+
+def mle_gev(x_gev: np.ndarray, start_loc=0, start_scale=1, start_shape=0):
+ assert np.ndim(x_gev) == 1
+ assert start_scale > 0
+ r = ro.r
+ x_gev = rpyn.numpy2ri(x_gev)
+ r.assign('x_gev', x_gev)
+ r.assign('python_wrapping', True)
+ r.source(file="/home/erwan/Documents/projects/spatiotemporalextremes/R/gev_fit/gev_fit.R")
+ res = r.mle_gev(loc=start_loc, scale=start_scale, shape=start_shape)
+ mle_params = dict(r.attr(res, 'coef').items())
+ return mle_params
+
+
+class GevMleFit(object):
+
+ def __init__(self, x_gev: np.ndarray, start_loc=0, start_scale=1, start_shape=0):
+ self.x_gev = x_gev
+ self.mle_params = mle_gev(x_gev, start_loc, start_scale, start_shape)
+ self.shape = self.mle_params[GEV_SHAPE]
+ self.scale = self.mle_params[GEV_SCALE]
+ self.location = self.mle_params[GEV_LOCATION]
\ No newline at end of file
diff --git a/test/R/test_gev_fit.py b/test/R/test_gev_fit.py
new file mode 100644
index 00000000..c2f4b955
--- /dev/null
+++ b/test/R/test_gev_fit.py
@@ -0,0 +1,40 @@
+import unittest
+
+from R.gev_fit.gev_marginal import frechet_unitary_transformation
+from R.gev_fit.gev_mle_fit import GevMleFit
+import rpy2.robjects as ro
+import numpy as np
+
+
+class TestGevFit(unittest.TestCase):
+
+ def test_unitary_mle_gev_fit(self):
+ r = ro.r
+ # Generate some sample from a gev
+ r.library('SpatialExtremes')
+ r("""
+ set.seed(42)
+ N <- 50
+ loc = 0; scale = 1; shape <- 1
+ x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
+ start_loc = 0; start_scale = 1; start_shape = 1
+ """)
+ # Get the MLE estimator
+ estimator = GevMleFit(x_gev=np.array(r['x_gev']),
+ start_loc=np.float(r['start_loc'][0]),
+ start_scale=np.float(r['start_scale'][0]),
+ start_shape=np.float(r['start_shape'][0]))
+ # Compare the MLE estimated parameters to the reference
+ mle_params_estimated = estimator.mle_params
+ mle_params_ref = {'loc': 0.0219, 'scale': 1.0347, 'shape': 0.8290}
+ for key in mle_params_ref.keys():
+ self.assertAlmostEqual(mle_params_ref[key], mle_params_estimated[key], places=3)
+
+ # def test_unitary_frechet_unitary_transformation(self):
+ # data = np.array([0.0, 1.0])
+ # # frechet_unitary_transformation()
+ # self.assertTrue(False)
+
+
+if __name__ == '__main__':
+ unittest.main()
\ No newline at end of file
--
GitLab