[contrasting] add quadratic & cross terms Gumbel altitudinal models & fix a...

[contrasting] add quadratic & cross terms Gumbel altitudinal models & fix a recurrent bug that happenning when checking the AIC value

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

extreme_fit/distribution/gev/summary_fevd_fixed.R

+# 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)
+}
+

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

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)

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)
 
 

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)