Loi exponentielle d’espérance 1.5 (de fonction de répartition F(z) = 1 − exp(−),
pour z > 0),
Loi de Pareto (de fonction de répartition F(z) = 1 − z−0.5 pour z > 1).
Question 1Pour chacune des lois ci-dessus: Trouver par le calcul les suitesrenormalisantes μnet σnqui assure la convergence du maximum. En déduire l’indice deforme des extrêmes ξ. Donner la loi généralisée des extrêmes qui approche la loi du maximumrenormalisé.
Question 2Simuler N = 100 échantillons du maximum de n = 12, 365, 1000 et
10000 valeurs pour les trois lois ci-dessus, puis tracer les histogrammes des maximumscorrespondants et la densité de la loi des extrêmes théoriques. On utilisera les fonctionsunif.sce, expo.sce, pareto.sce. Discuter de la qualité d’ajustement. (Pour créer et afficher lesgraphiques dans la fenêtre graphique i, utiliser la commande xset(``window'',i).)
Le vecteur unif(:,2) contient l’échantillon de N maxima de n = 365 simulations de la loi
uniforme.
1.2 Estimation dans le modèle GEV
On reprend les simulations du paragraphe précédent.
Question 3Écrire la vraisemblance et réaliser l’inférence par maximum de vraisemblance.Donner un intervalle de confiance pour μ,σ et ξ. On utilisera la fonction gev_fitdéfiniedans le fichier gev.sce.
La fonction gev_fit (y, init) calcule les estimateurs du maximum de vraisemblance de μ,σ
et ξ de l’échantillon y ainsi que leur intervalle de confiance. Le paramètre optionnel init
permet de spécifier les valeurs initiales de μ,σ et ξ utilisées pour l’optimisation de la
log-vraisemblance.
Question 4Relancer 10 fois les simulations, et comparer les estimations avec les valeursthéoriques.
2 Données Réelles
2.1 Estimation statistique de la crue centennale à Bar sur Seine (modèle POT)
A Bar sur Seine, sur 26 ans, 74 valeurs de débits de la Seine ont dépassé 250 (dixièmes de mm
équivalent lame d’eau sur le bassin versant). Voici les valeurs de dépassements et le nombre de
dépassement sur chaque année :
On ne considère dans un premier temps que les 15 premières années de données.
Question 5Réaliser l’inférence du modèle POT par maximum de vraisemblance et donnerun intervalle de confiance pour λ,σ et ξ. Pour l’estimation et l’intervalle de confiance deσ,ξ, utiliser fonction pot_fitdéfinie dans le fichier pot.sce.
[sigma, xi, cov] = pot_fit(x_15);
Question 6Peut on accepter l’hypothèse d’un dépassement de forme exponentielle ξ = 0?
Question 7Donner un intervalle de confiance pour le débit centennal à Bar sur Seine. Onpourra utiliser la fonction pot_quantile(y, f, p).
pot_quantile(x_15, f(1:15), 0.01)
Question 8Même travail avec toutes les données disponibles. Que constatez vous?
2.2 Estimation statistique sur un modèle de maximum par paquets
Les maxima annuels sur les 31 années de 1950 à 1980 de la Seine à Bar sur Seine sont
Question 9Écrire la vraisemblance et réaliser l’inférence par maximum de vraisemblance.Donner un intervalle de confiance pour pour μ,σ et ξ. De nouveau, on utilisera la fonctiongev_fit.
Question 10Donner un intervalle de confiance pour le débit maximal annuel centennal àBar sur Seine.
On utilisera la fonction gev_quantile qui prend en argument y, p, mu, sigma ,xi où y est
l’échantillon des maxima annuels, p=0.01 dans notre cas et mu, sigma ,xi sont les valeurs des
estimateurs trouvés à la question précédente.
Question 11On suppose μ et σ sont connus, donnés par l’évaluation à la question9.
Le coût unitaire de construction (surélévation) est de 100. On suppose que le coût desdommages est égal à la puissance 1.2 du débit non contenu par la digue et que la durée duprojet est de 50 ans. En utilisant un calcul par simulation de Monte-Carlo, quelle hauteurrecommandez vous si ξ = 0, si ξ = 0.1, si ξ = 0.3. Est-il plus pertinent de répondre à cettequestion avec le modèle GEV ou POT?
Pour cette question, on utilisera la fonction simul_gev(mu, sigma, xi, z) qui calcule par une
méthode de Monte Carlo (avec 10000 simulations) l’espérance de (D − z)+1.2 lorsque la variable
aléatoire D (représentant le débit) suit respectivement un modèle GEV de paramètres
mu, sigma, xi La fonction fournit également la précision de la méthode de Monte Carlo à
95%.
On pourra aussi utiliser la boucle:
xbasc(); z=[100:10:600]; for i=1:length(z) y(i)=50*simul_gev(mu, sigma, xi, z(i))+ 100*z(i); end plot(z,y);