From 4d582087c63236490ea06960149f4085e817ea20 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 22 Jan 2019 16:23:29 +0100
Subject: [PATCH] [GEV] modify gevmle_fit, use code from spatial_extreme
---
extreme_estimator/gev/fit_gev.py | 31 ---------------
.../gev/{fit_gev.R => gevmle_fit.R} | 4 ++
extreme_estimator/gev/gevmle_fit.py | 39 +++++++++++++++++++
...test_gev_mle_fit.py => test_gevmle_fit.py} | 14 ++++---
4 files changed, 51 insertions(+), 37 deletions(-)
delete mode 100644 extreme_estimator/gev/fit_gev.py
rename extreme_estimator/gev/{fit_gev.R => gevmle_fit.R} (86%)
create mode 100644 extreme_estimator/gev/gevmle_fit.py
rename test/test_extreme_estimator/test_gev/{test_gev_mle_fit.py => test_gevmle_fit.py} (68%)
diff --git a/extreme_estimator/gev/fit_gev.py b/extreme_estimator/gev/fit_gev.py
deleted file mode 100644
index 0a19e0cf..00000000
--- 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 124bd803..50d5ac61 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 00000000..21bb625d
--- /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 9b948f46..97b2d8ba 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}
--
GitLab