[GEV] Add gev_fit in R, and corresponding wrapper in Python

R/gev_fit/gev_fit.R

+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

R/gev_fit/gev_marginal.py

+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

R/gev_fit/gev_mle_fit.py

+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

test/R/test_gev_fit.py

+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