From aeac50d4eccceda568ef37c2c5005c4d9c298bbd Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 30 Jun 2020 17:03:05 +0200
Subject: [PATCH] [contrasting] add quadratic & cross terms Gumbel altitudinal
models & fix a recurrent bug that happenning when checking the AIC value
---
.../gev/main_fevd_mle_two_covariates.R | 9 +-
.../distribution/gev/summary_fevd_fixed.R | 83 +++++++++++++++++++
.../abstract_result_from_extremes.py | 2 +-
extreme_fit/model/utils.py | 4 +-
.../altitudes_fit/main_altitudes_studies.py | 11 +--
.../test_gev_spatio_temporal_extremes_mle.py | 3 +-
6 files changed, 99 insertions(+), 13 deletions(-)
create mode 100644 extreme_fit/distribution/gev/summary_fevd_fixed.R
diff --git a/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R b/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R
index fb450f26..49fcd28b 100644
--- a/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R
+++ b/extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R
@@ -8,6 +8,7 @@ library(stats4)
library(SpatialExtremes)
source('fevd_fixed.R')
source('ci_fevd_fixed.R')
+source('summary_fevd_fixed.R')
# Sample from a GEV
set.seed(42)
N <- 50
@@ -27,9 +28,11 @@ coord = data.frame(coord, stringsAsFactors = TRUE)
# res = fevd_fixed(x_gev, data=coord, method='MLE', verbose=TRUE, use.phi=FALSE)
# res = fevd_fixed(x_gev, data=coord, location.fun= ~T, scale.fun= ~T, method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
# res = fevd_fixed(x_gev, data=coord, location.fun= ~sin(X) + cos(T), method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
-res = fevd_fixed(x_gev, data=coord, location.fun= ~poly(X * T, 1, raw = TRUE), method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
-# res = fevd_fixed(x_gev, data=coord, location.fun= ~poly(X, 1, raw = TRUE) + poly(T, 2, raw = TRUE) , method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
-print(res)
+# res = fevd_fixed(x_gev, data=coord, location.fun= ~poly(X * T, 1, raw = TRUE), method='MLE', type="Gumbel", verbose=FALSE, use.phi=FALSE)
+res = fevd_fixed(x_gev, data=coord, location.fun= ~poly(X, 1, raw = TRUE) + poly(T, 2, raw = TRUE) , method='MLE', type="Gumbel", verbose=FALSE, use.phi=FALSE)
+# print(res)
+print(summary.fevd.mle_fixed(res))
+# print(summary(res)$AIC)
# Some display for the results
# m = res$results
diff --git a/extreme_fit/distribution/gev/summary_fevd_fixed.R b/extreme_fit/distribution/gev/summary_fevd_fixed.R
new file mode 100644
index 00000000..0f41a492
--- /dev/null
+++ b/extreme_fit/distribution/gev/summary_fevd_fixed.R
@@ -0,0 +1,83 @@
+# Title : TODO
+# Objective : TODO
+# Created by: erwan
+# Created on: 30/06/2020
+summary.fevd.mle_fixed <- function (object, ...)
+{
+ x <- object
+ a <- list(...)
+ out <- list()
+ cov.theta <- se.theta <- NULL
+ if (!is.null(a$silent))
+ silent <- a$silent
+ else silent <- FALSE
+ if (!silent) {
+ cat("\n")
+ print(x$call)
+ cat("\n")
+ print(paste("Estimation Method used: ", x$method, sep = ""))
+ cat("\n")
+ }
+ if (!silent)
+ cat("\n", "Negative Log-Likelihood Value: ", x$results$value,
+ "\n\n")
+ theta.hat <- x$results$par
+ theta.names <- names(theta.hat)
+ if (is.element("log.scale", theta.names)) {
+ theta.hat[theta.names == "log.scale"] <- exp(theta.hat[theta.names ==
+ "log.scale"])
+ theta.names[theta.names == "log.scale"] <- "scale"
+ names(theta.hat) <- theta.names
+ phiU <- FALSE
+ }
+ else phiU <- x$par.models$log.scale
+ out$par <- theta.hat
+ np <- length(theta.hat)
+ designs <- setup.design(x)
+ if (!is.null(x$data.pointer))
+ xdat <- get(x$data.pointer)
+ else xdat <- x$x
+ # cov.theta <- parcov.fevd(x)
+ # out$cov.theta <- cov.theta
+ # if (!is.null(cov.theta)) {
+ # se.theta <- sqrt(diag(cov.theta))
+ # names(se.theta) <- theta.names
+ # out$se.theta <- se.theta
+ # }
+ if (!silent) {
+ cat("\n", "Estimated parameters:\n")
+ print(theta.hat)
+ }
+ if (!is.null(se.theta)) {
+ if (!silent) {
+ cat("\n", "Standard Error Estimates:\n")
+ print(se.theta)
+ }
+ theta.hat <- rbind(theta.hat, se.theta)
+ if (is.matrix(theta.hat) && dim(theta.hat)[1] == 2)
+ rownames(theta.hat) <- c("Estimate", "Std. Error")
+ if (is.matrix(theta.hat))
+ colnames(theta.hat) <- theta.names
+ else names(theta.hat) <- theta.names
+ if (!silent && !is.null(cov.theta)) {
+ cat("\n", "Estimated parameter covariance matrix.\n")
+ print(cov.theta)
+ }
+ }
+ nllh <- x$results$value
+ out$nllh <- nllh
+ if (is.element(x$type, c("GEV", "Gumbel", "Weibull", "Frechet")))
+ n <- x$n
+ else {
+ y <- c(datagrabber(x, cov.data = FALSE))
+ n <- sum(y > x$threshold)
+ }
+ out$AIC <- 2 * nllh + 2 * np
+ out$BIC <- 2 * nllh + np * log(n)
+ if (!silent) {
+ cat("\n", "AIC =", out$AIC, "\n")
+ cat("\n", "BIC =", out$BIC, "\n")
+ }
+ invisible(out)
+}
+
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
index 8df785c8..ec86cb4a 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
@@ -24,7 +24,7 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
# Warning print will not work in this part
f = io.StringIO()
with redirect_stdout(f):
- summary = r('summary')(self.result_from_fit)
+ summary = r('summary.fevd.mle_fixed')(self.result_from_fit)
return self.get_python_dictionary(summary)
@property
diff --git a/extreme_fit/model/utils.py b/extreme_fit/model/utils.py
index 3d183036..89a11766 100644
--- a/extreme_fit/model/utils.py
+++ b/extreme_fit/model/utils.py
@@ -31,8 +31,8 @@ warnings.filterwarnings("ignore")
# Load ismev
r.library('ismev')
# Load fevd fixed
-for j, filename in enumerate(['ci_fevd_fixed.R', 'fevd_fixed.R', 'gnfit_fixed.R']):
- folder = 'gev' if j <= 1 else "gumbel"
+for j, filename in enumerate(['ci_fevd_fixed.R', 'fevd_fixed.R', 'summary_fevd_fixed.R', 'gnfit_fixed.R']):
+ folder = 'gev' if j <= 2 else "gumbel"
fevd_fixed_filepath = op.join(get_root_path(), 'extreme_fit', 'distribution', folder, filename)
assert op.exists(fevd_fixed_filepath)
r.source(fevd_fixed_filepath)
diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index 576b2ca2..baa704c6 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -61,21 +61,22 @@ def plot_moments(studies, massif_names=None):
def main():
- altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][4:7]
+ # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][4:7]
# altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
- # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
+ altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
SafranPrecipitation7Days][:]
- study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day, SafranSnowfall3Days, SafranPrecipitation3Days][:1]
+ study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day, SafranSnowfall3Days, SafranPrecipitation3Days][:]
massif_names = None
# massif_names = ['Aravis']
# massif_names = ['Chartreuse', 'Belledonne']
for study_class in study_classes:
+ print('change study class')
studies = AltitudesStudies(study_class, altitudes, season=Season.winter_extended)
- # plot_time_series(studies, massif_names)
- # plot_moments(studies, massif_names)
+ plot_time_series(studies, massif_names)
+ plot_moments(studies, massif_names)
plot_altitudinal_fit(studies, massif_names)
diff --git a/test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py b/test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py
index 6e147115..eae7fd6a 100644
--- a/test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py
+++ b/test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py
@@ -57,8 +57,7 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
self.common_test(model_class)
def test_altitudinal_gumbel_models(self):
- for model_class in ALTITUDINAL_GUMBEL_MODELS[:3]:
- # print(model_class)
+ for model_class in ALTITUDINAL_GUMBEL_MODELS[:]:
self.common_test(model_class)
--
GitLab