diff --git a/extreme_fit/distribution/gumbel/gnfit_fixed.R b/extreme_fit/distribution/gumbel/gnfit_fixed.R
new file mode 100644
index 0000000000000000000000000000000000000000..6822300aa8717bbb9ac3c54a01f605b84faf0bed
--- /dev/null
+++ b/extreme_fit/distribution/gumbel/gnfit_fixed.R
@@ -0,0 +1,131 @@
+# 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
diff --git a/extreme_fit/distribution/gumbel/gumbel_gof.R b/extreme_fit/distribution/gumbel/gumbel_gof.R
index 847a8246b533c52f63f54a5f429477b63fb2b257..69c6924acfaf6cfe17a2e93518e65df6791b6996 100644
--- a/extreme_fit/distribution/gumbel/gumbel_gof.R
+++ b/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)
diff --git a/extreme_fit/distribution/gumbel/gumbel_gof.py b/extreme_fit/distribution/gumbel/gumbel_gof.py
index 39aded04c91523a62ce8825eda9a08be0e2894f4..2e4ba56cd15f772774d441fde9df1d9f1c0adb63 100644
--- a/extreme_fit/distribution/gumbel/gumbel_gof.py
+++ b/extreme_fit/distribution/gumbel/gumbel_gof.py
@@ -1,3 +1,6 @@
+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))
diff --git a/extreme_fit/model/utils.py b/extreme_fit/model/utils.py
index 968b877d8825f95960fc30a254c48b777b3adf71..0047c85fa43ed68d99b34213c1aebf06ce429045 100644
--- a/extreme_fit/model/utils.py
+++ b/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
diff --git a/extreme_trend/abstract_gev_trend_test.py b/extreme_trend/abstract_gev_trend_test.py
index 2f79e7b32f2d52cc8105bf894550a2a3023e6bf0..feea4fd4d577ef15539bf32b302208f942a9b1d3 100644
--- a/extreme_trend/abstract_gev_trend_test.py
+++ b/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:
diff --git a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
index c22c2832ff1417cf7862a5e421e45662c23272f3..61a2fe2773d054ecc8cfe9bcd95f3eec28675ce0 100644
--- a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
+++ b/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,
diff --git a/projects/exceeding_snow_loads/section_results/plot_selection_curves.py b/projects/exceeding_snow_loads/section_results/plot_selection_curves.py
index da02e8ecc517ba0d7993fcaa4daa3c2a6e63cdb3..fde43a425ea11db9bead800bb4748ef87e0be13a 100644
--- a/projects/exceeding_snow_loads/section_results/plot_selection_curves.py
+++ b/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]