[GEV] modify gevmle_fit, use code from spatial_extreme

extreme_estimator/gev/fit_gev.R → extreme_estimator/gev/gevmle_fit.R

@@ -22,9 +22,13 @@ mle_gev <- function(loc, scale, shape){
 }
 
 main <- function(){
+    # My custom mle fit
     res = mle_gev(start_loc, start_scale, start_shape)
     print(attributes(res))
     print(attr(res, 'coef'))
+    # mle fit from spatial extremes
+    param = gevmle(x_gev, method = "Nelder")
+    print(param)
 }
 
 if (call_main) {

extreme_estimator/gev/fit_gev.py → extreme_estimator/gev/gevmle_fit.py

@@ -7,25 +7,33 @@ import os.path as op
 from extreme_estimator.gev_params import GevParams
 from extreme_estimator.extreme_models.utils import get_associated_r_file
 
+from extreme_estimator.extreme_models.utils import r
 
-def gev_mle_fit(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=get_associated_r_file(python_filepath=op.abspath(__file__)))
-    res = r.mle_gev(loc=start_loc, scale=start_scale, shape=start_shape)
-    mle_params = dict(r.attr(res, 'coef').items())
-    return mle_params
+
+# def gev_mle_fit(x_gev: np.ndarray, start_loc=0.0, start_scale=1.0, start_shape=1.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=get_associated_r_file(python_filepath=op.abspath(__file__)))
+#     print(start_loc, start_scale, start_shape)
+#     res = r.mle_gev(loc=start_loc, scale=start_scale, shape=start_shape)
+#     mle_params = dict(r.attr(res, 'coef').items())
+#     print('mle params', mle_params)
+#     return mle_params
+
+def spatial_extreme_gevmle_fit(x_gev):
+    res = r.gevmle(x_gev, method="Nelder")
+    return dict(zip(res.names, res))
 
 
 class GevMleFit(object):
 
-    def __init__(self, x_gev: np.ndarray, start_loc=0, start_scale=1, start_shape=0):
+    def __init__(self, x_gev: np.ndarray):
         self.x_gev = x_gev
-        self.mle_params = gev_mle_fit(x_gev, start_loc, start_scale, start_shape)
+        self.mle_params = spatial_extreme_gevmle_fit(x_gev)
         self.shape = self.mle_params[GevParams.GEV_SHAPE]
         self.scale = self.mle_params[GevParams.GEV_SCALE]
         self.loc = self.mle_params[GevParams.GEV_LOC]

test/test_extreme_estimator/test_gev/test_gev_mle_fit.py → test/test_extreme_estimator/test_gev/test_gevmle_fit.py

@@ -3,12 +3,12 @@ import unittest
 import numpy as np
 
 from extreme_estimator.extreme_models.utils import r, set_seed_r
-from extreme_estimator.gev.fit_gev import GevMleFit
+from extreme_estimator.gev.gevmle_fit import GevMleFit
 
 
 class TestGevMleFit(unittest.TestCase):
 
-    def test_unitary_gev_mle_fit(self):
+    def setUp(self) -> None:
         set_seed_r()
         r("""
         N <- 50
@@ -16,11 +16,13 @@ class TestGevMleFit(unittest.TestCase):
         x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
         start_loc = 0; start_scale = 1; start_shape = 1
         """)
+
+    def test_gevmle_fit(self):
         # 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]))
+        estimator = GevMleFit(x_gev=np.array(r['x_gev']))
+        self.fit_estimator(estimator)
+
+    def fit_estimator(self, estimator):
         # 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}