From 33f31bc8bd11c1edfc138579ea7bcf6f8443bfe4 Mon Sep 17 00:00:00 2001 From: Le Roux Erwan <erwan.le-roux@irstea.fr> Date: Thu, 13 Dec 2018 15:35:44 +0100 Subject: [PATCH] [THESIS REPORT] add quantile plot --- example.R | 35 ----------------------------------- thesis_report/gev_plot.py | 22 +++++++++++++++++++--- 2 files changed, 19 insertions(+), 38 deletions(-) delete mode 100644 example.R diff --git a/example.R b/example.R deleted file mode 100644 index 0083b6ea..00000000 --- a/example.R +++ /dev/null @@ -1,35 +0,0 @@ -library(SpatialExtremes) -# Title : TODO -# Objective : TODO -# Created by: erwan -# Created on: 27/11/18 -# define the coordinate of each location -set.seed(42) -n.site <- 30 -locations <- matrix(runif(2*n.site, 0, 10), ncol = 2) -colnames(locations) <- c("lon", "lat") -print(locations) -##Simulate a max-stable process - with unit Frechet margins -data <- rmaxstab(40, locations, cov.mod = "whitmat", nugget = 0, range = 3, -smooth = 0.5) -##Now define the spatial model for the GEV parameters -param.loc <- -10 + 2 * locations[,2] -param.scale <- 5 + 2 * locations[,1] + locations[,2]^2 -param.shape <- rep(0.2, n.site) -##Transform the unit Frechet margins to GEV -for (i in 1:n.site) -data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i], param.shape[i]) -##Define a model for the GEV margins to be fitted -##shape ~ 1 stands for the GEV shape parameter is constant -##over the region -loc.form <- loc ~ lat -scale.form <- scale ~ lon + I(lat^2) -shape.form <- shape ~ 1 -##Fit a max-stable process using the Schlather -res = fitmaxstab(data, locations, "whitmat", loc.form, scale.form, shape.form) -print(res['fitted.values']) -## Model without any spatial structure for the GEV parameters -## Be careful this could be *REALLY* time consuming -# res = fitmaxstab(data, locations, "whitmat", fit.margin=TRUE) -# print(res['fitted.values']) - diff --git a/thesis_report/gev_plot.py b/thesis_report/gev_plot.py index d2d6b9e5..a5916cc0 100644 --- a/thesis_report/gev_plot.py +++ b/thesis_report/gev_plot.py @@ -28,7 +28,7 @@ def max_stable_plot(): lim = 10 ** (power_n) x = np.linspace(0, lim, 100) loc, scale, shape = 1, 1, 1 - for n in [10**i for i in range(power_n)]: + for n in [10 ** i for i in range(power_n)]: label = 'n={}'.format(n) y = np.array(r.pgev(x, loc, scale, shape)) y **= n @@ -40,6 +40,22 @@ def max_stable_plot(): plt.show() +def quantile_function_plot(): + p = np.linspace(5e-1, 1e-5, 100) + loc, scale = 2, 1 + shapes = [-1, 0, 1] + for shape in shapes: + label = '$\zeta= {} $'.format(shape) + funct = [r.qgev, r.qgpd][0] + y = funct(1 - p, loc, scale, shape) + plt.loglog(1/p, y, basex=10, label=label) + plt.legend() + plt.xlabel('$1/p$') + plt.ylabel('$q(1- p|\mu, \sigma, \zeta)$') + plt.show() + + if __name__ == '__main__': -# gev_plot() - max_stable_plot() + # gev_plot() + # max_stable_plot() + quantile_function_plot() \ No newline at end of file -- GitLab