[THESIS REPORT] add quantile plot

example.R

-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'])
-

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