diff --git a/extreme_estimator/gev/fit_gev.py b/extreme_estimator/gev/fit_gev.py
deleted file mode 100644
index 0a19e0cfe5f50209629bc03a5697777cef103ffb..0000000000000000000000000000000000000000
--- a/extreme_estimator/gev/fit_gev.py
+++ /dev/null
@@ -1,31 +0,0 @@
-import rpy2.robjects as ro
-import numpy as np
-import rpy2.robjects.numpy2ri as rpyn
-import os.path as op
-
-# Defining some constants
-from extreme_estimator.gev_params import GevParams
-from extreme_estimator.extreme_models.utils import get_associated_r_file
-
-
-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
-
-
-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 = gev_mle_fit(x_gev, start_loc, start_scale, start_shape)
-        self.shape = self.mle_params[GevParams.GEV_SHAPE]
-        self.scale = self.mle_params[GevParams.GEV_SCALE]
-        self.loc = self.mle_params[GevParams.GEV_LOC]
diff --git a/extreme_estimator/gev/fit_gev.R b/extreme_estimator/gev/gevmle_fit.R
similarity index 86%
rename from extreme_estimator/gev/fit_gev.R
rename to extreme_estimator/gev/gevmle_fit.R
index 124bd803133e436eb1c889759221a2f9d249b87c..50d5ac61a189cfa26a9f48137242e6c384d191c7 100644
--- a/extreme_estimator/gev/fit_gev.R
+++ b/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) {
diff --git a/extreme_estimator/gev/gevmle_fit.py b/extreme_estimator/gev/gevmle_fit.py
new file mode 100644
index 0000000000000000000000000000000000000000..21bb625d2623c5335b0c06a83128b1b9ae83af6c
--- /dev/null
+++ b/extreme_estimator/gev/gevmle_fit.py
@@ -0,0 +1,39 @@
+import rpy2.robjects as ro
+import numpy as np
+import rpy2.robjects.numpy2ri as rpyn
+import os.path as op
+
+# Defining some constants
+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.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):
+        self.x_gev = x_gev
+        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]
diff --git a/test/test_extreme_estimator/test_gev/test_gev_mle_fit.py b/test/test_extreme_estimator/test_gev/test_gevmle_fit.py
similarity index 68%
rename from test/test_extreme_estimator/test_gev/test_gev_mle_fit.py
rename to test/test_extreme_estimator/test_gev/test_gevmle_fit.py
index 9b948f469902f474fbe941d432ec522455a3c7c8..97b2d8ba0999b712deb62ea7475c0069bb9aea16 100644
--- a/test/test_extreme_estimator/test_gev/test_gev_mle_fit.py
+++ b/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}