[exceeding] add full diagnositic plot (qqplot & intensity plots) for gumbel...

[exceeding] add full diagnositic plot (qqplot & intensity plots) for gumbel stationary & add gumbel gof tests from r package

extreme_fit/distribution/gumbel/gnfit_fixed.R

+# Title     : TODO
+# Objective : TODO
+# Created by: erwan
+# Created on: 15/05/2020
+##
+## Exported symobls in package `gnFit`
+##
+
+## Exported package methods
+
+gnfit_fixed <- function (dat, dist, df = NULL, pr = NULL, threshold = NULL)
+{
+    dat <- as.numeric(dat)
+    x <- NULL
+    z <- list()
+    op <- par(mfrow = c(1, 2))
+    if (is.null(pr)) {
+        if (dist == "gev" | dist == "gpd")
+            stop("Enter Parameters!")
+        else if (dist == "t") {
+            if (df > 2) {
+                loc <- mean(dat)
+                sc <- sqrt((df - 2) * var(dat)/df)
+                xdat <- (dat - loc)/sc
+                prob <- pt(xdat, df)
+            }
+            else stop("DF must be > 2")
+        }
+        else if (dist == "gum") {
+            sc <- sqrt(6 * var(dat)/pi^2)
+            loc <- mean(dat) - 0.577 * sc
+            pr <- c(loc, sc, 0)
+            prob <- gevf(pr, dat)
+        }
+        else {
+            loc <- mean(dat)
+            ifelse(dist == "norm", sc <- sd(dat), NA)
+            ifelse(dist == "laplace", sc <- sqrt(var(dat)/2),
+                NA)
+            ifelse(dist == "logis", sc <- sqrt(3 * var(dat)/pi^2),
+                NA)
+            prob <- get(paste0("p", dist))(dat, loc, sc)
+        }
+    }
+    else {
+        if (dist == "gev") {
+            prob <- gevf(pr, dat)
+        }
+        else if (dist == "gum") {
+            pr[3] <- 0
+            prob <- gevf(pr, dat)
+        }
+        else if (dist == "gpd")
+            if (!is.null(threshold)) {
+                u <- threshold
+                dat <- dat[dat > u]
+                prob <- gpdf(pr, u, dat)
+            }
+            else stop("threshold is missing!")
+    }
+    n <- length(dat)
+    k <- seq(1:n)
+    qnor <- qnorm(sort(prob))
+    pnor <- pnorm((qnor - mean(qnor))/sd(qnor))
+    w <- round((sum((pnor - (2 * k - 1)/(2 * n))^2) + 1/(12 *
+        n)) * (1 + 0.5/n), 4)
+    if (w < 0.0275) {
+        pval <- 1 - exp(-13.953 + 775.5 * w - 12542.61 * w^2)
+    }
+    else if (w < 0.051) {
+        pval <- 1 - exp(-5.903 + 179.546 * w - 1515.29 * w^2)
+    }
+    else if (w < 0.092) {
+        pval <- exp(0.886 - 31.62 * w + 10.897 * w^2)
+    }
+    else if (w < 1.1) {
+        pval <- exp(1.111 - 34.242 * w + 12.832 * w^2)
+    }
+    else {
+        warning("p-value is smaller than 7.37e-10")
+    }
+    z$Wpval <- pval
+    A <- (-n - sum((2 * k - 1) * log(pnor) + (2 * n + 1 - 2 *
+        k) * log(1 - pnor))/n) * (1 + 0.75/n + 2.25/n^2)
+    A <- round((1 + 0.75/n + 2.25/n^2) * A, 4)
+    if (A < 0.2) {
+        pval <- 1 - exp(-13.436 + 101.14 * A - 223.73 * A^2)
+    }
+    else if (A < 0.34) {
+        pval <- 1 - exp(-8.318 + 42.796 * A - 59.938 * A^2)
+    }
+    else if (A < 0.6) {
+        pval <- exp(0.9177 - 4.279 * A - 1.38 * A^2)
+    }
+    else if (A < 10) {
+        pval <- exp(1.2937 - 5.709 * A + 0.0186 * A^2)
+    }
+    else {
+        warning("p-value is smaller than 7.37e-10")
+    }
+    z$Apval <- pval
+    z$Cram <- w
+    z$Ander <- A
+    invisible(z)
+}
+
+
+
+
+## Package Data
+
+# none
+
+
+## Package Info
+
+.skeleton_package_title = "Goodness of Fit Test for Continuous Distribution Functions"
+
+.skeleton_package_version = "0.2.0"
+
+.skeleton_package_depends = ""
+
+.skeleton_package_imports = "ismev,rmutil"
+
+
+## Internal
+
+.skeleton_version = 5
+
+
+## EOF

extreme_fit/distribution/gumbel/gumbel_gof.R

@@ -4,9 +4,10 @@
 # Created on: 15/05/2020
 # library(rmutil)
 library(evd)
-library(gnFit)
+library(ismev)
+source("gnfit_fixed.R")
 data = rgumbel(10)
-res = gnfit(data, "gum")
+res = gnfit_fixed(data, "gum")
 print(res)
 # r = rnorm(0)
 # disp(r)

extreme_fit/distribution/gumbel/gumbel_gof.py

+import io
+from contextlib import redirect_stdout
+
 import numpy as np
 
 from extreme_fit.model.result_from_model_fit.abstract_result_from_model_fit import AbstractResultFromModelFit
@@ -5,11 +8,11 @@ from extreme_fit.model.utils import r
 
 
 def cramer_von_mises_and_anderson_darling_tests_pvalues_for_gumbel_distribution(data):
-    res = r.gnfit(data, "gum")
+    res = r.gnfit_fixed(data, "gum")
     res = AbstractResultFromModelFit.get_python_dictionary(res)
     res = {k: np.array(v)[0] for k, v in res.items()}
     return res['Wpval'], res['Apval']
 
 
 if __name__ == '__main__':
-    cramer_von_mises_and_anderson_darling_tests_pvalues_for_gumbel_distribution(np.array([2.0, 3.0]))
+    cramer_von_mises_and_anderson_darling_tests_pvalues_for_gumbel_distribution(np.arange(50))

extreme_fit/model/utils.py

@@ -25,15 +25,15 @@ pandas2ri.activate()
 r.library('SpatialExtremes')
 r.library('data.table')
 r.library('quantreg')
-r.library('gnFit')
 # Desactivate temporarily warnings
 default_filters = warnings.filters.copy()
 warnings.filterwarnings("ignore")
 # Load ismev
 r.library('ismev')
 # Load fevd fixed
-for filename in ['ci_fevd_fixed.R', 'fevd_fixed.R']:
-    fevd_fixed_filepath = op.join(get_root_path(), 'extreme_fit', 'distribution', 'gev', filename)
+for j, filename in enumerate(['ci_fevd_fixed.R', 'fevd_fixed.R', 'gnfit_fixed.R']):
+    folder = 'gev' if j <= 1 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)
 # Reactivate warning

extreme_trend/abstract_gev_trend_test.py

@@ -9,6 +9,8 @@ from scipy.stats.stats import KstestResult
 
 from extreme_data.eurocode_data.utils import EUROCODE_QUANTILE, YEAR_OF_INTEREST_FOR_RETURN_LEVEL
 from extreme_data.meteo_france_data.scm_models_data.crocus.crocus_variables import AbstractSnowLoadVariable
+from extreme_fit.distribution.gumbel.gumbel_gof import \
+    cramer_von_mises_and_anderson_darling_tests_pvalues_for_gumbel_distribution
 from extreme_fit.estimator.margin_estimator.utils import fitted_linear_margin_estimator
 from extreme_fit.distribution.gev.gev_params import GevParams
 from extreme_fit.model.margin_model.utils import \
@@ -72,8 +74,10 @@ class AbstractGevTrendTest(object):
         # significance_level=array([25. , 10. ,  5. ,  2.5,  1. ]))
         index_for_significance_level_5_percent = 2
         quantiles = self.compute_empirical_quantiles(estimator=self.unconstrained_estimator)
-        test_res = anderson(quantiles, dist='gumbel_r')  # type: AndersonResult
-        return test_res.statistic < test_res.critical_values[index_for_significance_level_5_percent]
+        # test_res = anderson(quantiles, dist='gumbel_r')  # type: AndersonResult
+        # return test_res.statistic < test_res.critical_values[index_for_significance_level_5_percent]
+        test = cramer_von_mises_and_anderson_darling_tests_pvalues_for_gumbel_distribution(quantiles)
+        return test[1] > self.SIGNIFICANCE_LEVEL
 
     @property
     def degree_freedom_chi2(self) -> int:

projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py

@@ -89,8 +89,8 @@ def major_result():
     #                                  ModelSubsetForUncertainty.non_stationary_gumbel_and_gev]
     model_subsets_for_uncertainty = None
     # study_classes = [CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][::-1]
-    # altitudes = paper_altitudes
-    altitudes = [900, 1800, 2700][:1]
+    altitudes = paper_altitudes
+    # altitudes = [900, 1800, 2700][:1]
     for study_class in study_classes:
         intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
                             uncertainty_methods, study_class,

projects/exceeding_snow_loads/section_results/plot_selection_curves.py

@@ -20,8 +20,8 @@ def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNo
     ax = create_adjusted_axes(1, 1)
 
     selected_counter = merge_counter([v.selected_trend_test_class_counter for v in altitude_to_visualizer.values()])
-    selected_and_significative_counter = merge_counter([v.selected_and_significative_trend_test_class_counter for v in altitude_to_visualizer.values()])
-    # selected_and_significative_counter = merge_counter([v.selected_and_anderson_goodness_of_fit_trend_test_class_counter for v in altitude_to_visualizer.values()])
+    # selected_and_significative_counter = merge_counter([v.selected_and_significative_trend_test_class_counter for v in altitude_to_visualizer.values()])
+    selected_and_significative_counter = merge_counter([v.selected_and_anderson_goodness_of_fit_trend_test_class_counter for v in altitude_to_visualizer.values()])
     # selected_and_significative_counter = merge_counter([v.selected_and_kstest_goodness_of_fit_trend_test_class_counter for v in altitude_to_visualizer.values()])
     total_of_selected_models = sum(selected_counter.values())
     l = sorted(enumerate(visualizer.non_stationary_trend_test), key=lambda e: selected_counter[e[1]])
@@ -32,8 +32,8 @@ def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNo
     labels = permute(['${}$'.format(t.label) for t in visualizer.non_stationary_trend_test], permutation)
 
     print(l)
-    print(select_list)
-    print(selected_and_signifcative_list)
+    print(sum(select_list), select_list)
+    print(sum(selected_and_signifcative_list), selected_and_signifcative_list)
     # [(5, <    class 'data.trend_analysis.univariate_test.extreme_trend.trend_test_two_parameters.gev_trend_test_two_parameters.GevLocationAgainstGumbel'> ), (6, < class 'data.trend_analysis.univariate_test.extreme_trend.trend_test_two_parameters.gev_trend_test_two_parameters.GevScaleAgainstGumbel' > ), (2, < class 'data.trend_analysis.univariate_test.extreme_trend.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelScaleTrendTest' > ), (1, < class 'data.trend_analysis.univariate_test.extreme_trend.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelLocationTrendTest' > ), (7, < class 'data.trend_analysis.univariate_test.extreme_trend.trend_test_three_parameters.gev_trend_test_three_parameters.GevLocationAndScaleTrendTestAgainstGumbel' > ), (3, < class 'data.trend_analysis.univariate_test.extreme_trend.trend_test_two_parameters.gumbel_test_two_parameters.GumbelLocationAndScaleTrendTest' > ), (4, < class 'data.trend_analysis.univariate_test.extreme_trend.trend_test_one_parameter.gumbel_trend_test_one_parameter.GevStationaryVersusGumbel' > ), (0, < class 'data.trend_analysis.univariate_test.extreme_trend.trend_test_one_parameter.gumbel_trend_test_one_parameter.GumbelVersusGumbel' > )]
     # [32.64462809917355, 24.380165289256198, 12.396694214876034, 9.50413223140496, 9.090909090909092, 5.785123966942149, 3.71900826446281, 2.479338842975207]
     # [0, 13.223140495867769, 7.851239669421488, 8.264462809917354, 4.958677685950414, 2.479338842975207, 0.8264462809917356, 2.0661157024793386]