| ... | ... | @@ -45,26 +45,27 @@ Une [version plus aboutie est mise en ligne dans un autre projet](https://gitlab |
|
|
|
|
|
|
|
|
|
|
|
# CODE **ProbaCruesMaxAn_SurNannees.py**
|
|
|
|
Il s'agit de probabilités discrètes, calculées sur des valeurs de p entiers, avec deux paramètres N et T que l'on souhaite faire varier.
|
|
|
|
On peut donc représenter ces probabilités en fonction de p soit pas **plot** avec des markers non reliés, soit par une fonction de type "step". On va présenter ici le premier cas car les courbes tracées avec **plot** sont plus faciles à modifier, et on présentera à la fin une variante avec une graphique en escalier.
|
|
|
|
Il s'agit de probabilités discrètes, calculées sur des valeurs **k** entières, avec deux paramètres **N** et **T** que l'on souhaite faire varier.
|
|
|
|
On peut donc représenter ces probabilités en fonction de k soit par **plot** avec des markers non reliés, soit par une fonction de type "step". On va présenter ici le premier cas car les courbes tracées avec **plot** sont plus faciles à modifier, et on présentera à la fin une variante avec une graphique en escalier.
|
|
|
|
|
|
|
|
## 1. La fonction et ses paramètres
|
|
|
|
On rappelle que la période de retour T est l’inverse de la fréquence au dépassement : T= 1/f.
|
|
|
|
On raisonne avec les crues maximales annuelles, on ne retient de chaque année que la plus forte occurence de crue. Une crue centennale est dépassée en moyenne une fois tous les 100 ans et a donc une chance sur 100 d’être dépassée une année donnée : T=100 ans et f=1/100.
|
|
|
|
On rappelle que la période de retour **T** est l’inverse de **f**, la fréquence au dépassement : **T= 1/f**.
|
|
|
|
On raisonne avec les crues maximales annuelles, on ne retient de chaque année que la plus forte occurrence de crue. Une crue centennale est dépassée en moyenne une fois tous les 100 ans et a donc une chance sur 100 d’être dépassée une année donnée : T=100 ans et f=1/100.
|
|
|
|
|
|
|
|
Si on illustre par un arbre des possibles, pour avoir p crues supérieures au niveau fixé il faut retenir les réalisations pour lesquelles on a exactement p crues supérieures au niveau voulu, et leur nombre est égal au nombre de façons de sélectionner p positions parmi N, soit C(n,p). Pour ces C(n,p) réalisations, on a p fois un tirage de probabilité f et le reste du temps, soit (N-p) tirages, un tirage de probabilité complémentaire (1-f).
|
|
|
|
Si on illustre par un arbre des possibles, le nombre de réalisations pour lesquelles on a exactement **k** crues supérieures au niveau voulu est égal au nombre de façons de sélectionner **k** positions parmi **N**, soit **C(N,k)**. Pour ces **C(N,k)** réalisations, on a **k** fois un tirage de probabilité **f** et le reste du temps, soit **(N-p)** tirages, un tirage de probabilité complémentaire **(1-f)**.
|
|
|
|
|
|
|
|
La probabilité d’avoir sur N années données exactement p crues qui dépassent le niveau de la crue de période de retour T=1/f se calcule donc par la formule : C(n,p).f^p.(1-f)^(N-p).
|
|
|
|
La probabilité d’avoir sur **N** années données exactement **k** crues qui dépassent le niveau de la crue de période de retour T=1/f se calcule donc par la formule : **C(N,k).f^k.(1-f)^(N-k).**
|
|
|
|
|
|
|
|
``` python
|
|
|
|
from scipy.special import comb
|
|
|
|
def proba_crue(n, p, f):
|
|
|
|
""" Probabilité d’avoir exactement p crues de fréquence f en n années : C(n,p)*f^p*(1-f)^(N-p)
|
|
|
|
def proba_crue(n, k, f):
|
|
|
|
""" Probabilité d’avoir exactement k crues de fréquence f en n années : C(n,k)*f^p*(1-f)^(N-k)
|
|
|
|
"""
|
|
|
|
return comb(n, p)*(f**p)*((1-f)**(n-p))
|
|
|
|
return comb(n, k)*(f**k)*((1-f)**(n-k))
|
|
|
|
```
|
|
|
|
On représente la probabilité ainsi calculée en fonction de p, pour N et T fixés.
|
|
|
|
A chaque fois que l'on a une probabilité pour p, on peut aussi calculer par cumul la probabilité d'avoir entre 0 et p crues de période de retour supérieure à N : quand cette probabilité cumulée dépasse un seuil, ici 0.999 on arrête la boucle sur p. les entiers p à partir de 0 jusqu’à ce que ‘probabilité pour p crues ou moins de p’ dépasse 0,999 (courbe orange sur la copie d’écran suivante).
|
|
|
|
On représente la probabilité ainsi calculée en fonction de k, pour N et T fixés.
|
|
|
|
A chaque fois que l'on a une probabilité pour k, on peut aussi calculer par cumul la probabilité d'avoir entre 0 et k crues de période de retour supérieure à N (courbe orange sur la copie d’écran suivante).
|
|
|
|
Il n'est pas utile d'explorer tous les k de 0 à N : quand la probabilité cumulée dépasse un seuil, ici 0.999, on arrête la boucle sur k. .
|
|
|
|
|
|
|
|
``` python
|
|
|
|
def calcul_courbes(n,f):
|
| ... | ... | @@ -72,24 +73,24 @@ def calcul_courbes(n,f): |
|
|
|
cumul_au_plus = 0
|
|
|
|
liste_probas_exactement=[]
|
|
|
|
liste_probas_auplus = []
|
|
|
|
p= 0
|
|
|
|
k= 0
|
|
|
|
while cumul_au_plus < 0.999:
|
|
|
|
proba = proba_crue(n,p,f)
|
|
|
|
proba = proba_crue(n,k,f)
|
|
|
|
liste_probas_exactement.append(proba)
|
|
|
|
cumul_au_plus += proba
|
|
|
|
liste_probas_auplus.append(cumul_au_plus)
|
|
|
|
p+=1
|
|
|
|
k+=1
|
|
|
|
|
|
|
|
return p-1, liste_probas_exactement, liste_probas_auplus
|
|
|
|
return k-1, liste_probas_exactement, liste_probas_auplus
|
|
|
|
```
|
|
|
|
|
|
|
|
On définit *N* et *T* dès le départ, avec une valeur par défaut.
|
|
|
|
On définit **N** et **T** dès le départ, avec une valeur par défaut.
|
|
|
|
|
|
|
|
On calcule les probabilités avex ces valeurs par défaut, et on recueille les résultats dans 3 variables qui seront utilisées pour le tracé.
|
|
|
|
On calcule les probabilités avec ces valeurs par défaut, et on recueille les résultats dans 3 variables qui seront utilisées pour le tracé.
|
|
|
|
``` python
|
|
|
|
p, liste_probas_exactement, liste_probas_auplus = calcul_courbes(n, 1 / T)
|
|
|
|
k, liste_probas_exactement, liste_probas_auplus = calcul_courbes(n, 1 / T)
|
|
|
|
```
|
|
|
|
Deux sliders vont permettre ensuite de modifier N et T ; on les utilise ensuite pour calculer et tracer des résultats via une fonction ; il faut donc que la modification soit répercutée dans tout le code. Cela va nous permettre de manipuler la notion de portée des variables.
|
|
|
|
Deux **sliders** vont permettre ensuite de modifier **N** et **T** ; on va relancer les calculs et tracer des résultats dès qu'ils seront modifiés. Ici, on a choisi de répercuter la modification aux vraibles N et T définies dans le corps du programme. Ce n'est pas forcément la meilleure façon de faire, mais pour un petit code cela ne paose pas de problème et cela va nous permettre de manipuler la notion de portée des variables.
|
|
|
|
|
|
|
|
|
|
|
|
## 2. La figure
|
| ... | ... | @@ -120,20 +121,20 @@ plt.show() |
|
|
|
En haut, on a donc une vignette "ax" dans laquelle on va tracer deux courbes avec "plot", avec les résultats qui viennent d'être calculés. On nomme le premier objet retourné par la fonction "plot", puisqu'on aura besoin d'agir sur lui : on remplacera les données (vecteurs des x et vecteurs des y) après chaque action sur les sliders.
|
|
|
|
|
|
|
|
``` python
|
|
|
|
courbe_exactement, = ax.plot(range(p + 1), liste_probas_exactement, label="exactement p crues",ls='None', marker='o', markersize=10))
|
|
|
|
courbe_cumul, = ax.plot(range(p + 1), liste_probas_auplus, label="p crues ou moins de p",ls='None', marker='o', markersize=10)).
|
|
|
|
courbe_exactement, = ax.plot(range(p + 1), liste_probas_exactement, label="exactement k crues",ls='None', marker='o', markersize=10))
|
|
|
|
courbe_cumul, = ax.plot(range(p + 1), liste_probas_auplus, label="k crues ou moins de k",ls='None', marker='o', markersize=10)).
|
|
|
|
```
|
|
|
|
|
|
|
|
### définition des sliders et liaison avec une fonction "par widget"mise à jour de la courbe"
|
|
|
|
### définition des sliders et liaison avec une fonction par widget de "mise à jour de la courbe"
|
|
|
|
|
|
|
|
Sous la vignette "ax" il y a une vignette "ax_espace" dans laquelle on n'ajoute rien, puis les vignettes des deux sliders.
|
|
|
|
Sous la vignette "ax" il y a une vignette intercalaire "ax_espace" dans laquelle on n'ajoute rien, puis les vignettes des deux sliders.
|
|
|
|
Les arguments utilisés ici pour les sliders sont :
|
|
|
|
|
|
|
|
nom_du_slider = Slider(
|
|
|
|
ax_ou_le_slider_est_installe, "titre", valmin=valeur_mini, valmax = valeur_max, valstep=incrément, valfmt=format précédé de % : '%0.2f', valinit=valeur de départ, color="nom_de_couleur")
|
|
|
|
|
|
|
|
|
|
|
|
Comme on veut des entiers, on définit pour valstep un incrément égal à 1, et comme ce sont des flottants on définit un affichage avec zéro chiffre après la virgule.
|
|
|
|
Comme on veut des entiers, on définit pour valstep un incrément égal à 1, et comme la valeur prise sera quand même un flottant on l'affiche avec zéro chiffre après la virgule.
|
|
|
|
|
|
|
|
On définit enfin la fonction qui sera appelée si on change le curseur du slider :
|
|
|
|
nom_du_slider.on_changed(fonction_appellée)
|
| ... | ... | @@ -155,47 +156,45 @@ slider_T.on_changed(update_graphique) |
|
|
|
|
|
|
|
### définition de la fonction commune associée aux deux widgets
|
|
|
|
|
|
|
|
On a associé la même fonction, update_graphique, aux deux sliders car finalement ils vont tous les deux déclancher la même suite d'actions : recalcul et retracé.
|
|
|
|
On a associé la même fonction, update_graphique, aux deux sliders car finalement ils vont tous les deux déclencher la même suite d'actions : recalcul et retracé.
|
|
|
|
|
|
|
|
Le slider_N (respectivement slider_T) récupère une valeur numérique et la transmet en tant qu'argument *val* à la fonction qui lui a été liée.
|
|
|
|
Les fonctions liées à un widget ne peuvent prendre qu'un argument, imposé par le type du widget. Pour un slider, l'argument ne peut être que la valeur du curseur, val par convention.
|
|
|
|
|
|
|
|
Cependant, on va appeler une fonction qui va retracer les courbes en fonction de N et T ; il faut donc connaître ces deux variables dans cette fonction, et en modifier une (celle du curseur actif) pour tout le code.
|
|
|
|
|
|
|
|
On a choisi de préciser dans la fonction que l'on travaille avec les deux variables N et T au niveau global:
|
|
|
|
global N, T
|
|
|
|
|
|
|
|
Il n'est pas possible d'utiliser l'argument passé val car dans la fonction appelée on ne sait pas quel est le curseur qui a été activé, cela peut être slider_N ou slider_T.
|
|
|
|
Cependant, on appelle cette fonction soit depuis slider_N soit depuis slider_T : il n'est pas possible d'utiliser l'argument passé val car dans la fonction appelée on ne sait pas quel est le curseur qui a été activé, cela peut être slider_N ou slider_T.
|
|
|
|
Heureusement, on peut appeler directement la valeur de chaque curseur en tant qu'attribut ([ source utilisée ](http://www.math.buffalo.edu/~badzioch/MTH337/PT/PT-matplotlib_slider/PT-matplotlib_slider.html)):
|
|
|
|
|
|
|
|
N = slider_N.val
|
|
|
|
n = slider_N.val
|
|
|
|
T = slider_T.val
|
|
|
|
|
|
|
|
Enfin, on a choisi que les sliders modifient la valeur de N et de T pour tout le code. Pour cela, on précise dans la fonction que l'on travaille avec les deux variables N et T au niveau global, avec l'instruction
|
|
|
|
**global** N, T
|
|
|
|
|
|
|
|
``` python
|
|
|
|
def update_graphique(val):
|
|
|
|
global n, T
|
|
|
|
n = int(slider_n.val) # on précise que ce sont des entiers pour l'utilisation dans C(n,p)
|
|
|
|
T = int(slider_T.val)
|
|
|
|
|
|
|
|
p, liste_probas_exactement, liste_probas_auplus = calcul_courbes(n,1/T)
|
|
|
|
fig.suptitle(f"Démo : probabilité d'avoir p crues > T={T}ans en {n} années")
|
|
|
|
courbe_exactement.set_data(range(p + 1), liste_probas_exactement)
|
|
|
|
courbe_cumul.set_data(range(p + 1), liste_probas_auplus)
|
|
|
|
ax.set_xlim(0, p+1)
|
|
|
|
k, liste_probas_exactement, liste_probas_auplus = calcul_courbes(n,1/T)
|
|
|
|
fig.suptitle(f"Démo : probabilité d'avoir k crues > T={T}ans en {n} années")
|
|
|
|
courbe_exactement.set_data(range(k + 1), liste_probas_exactement)
|
|
|
|
courbe_cumul.set_data(range(k + 1), liste_probas_auplus)
|
|
|
|
ax.set_xlim(0, k+1)
|
|
|
|
fig.canvas.draw_idle()
|
|
|
|
```
|
|
|
|
|
|
|
|
Le recalcul ne présente pas de difficulté ; on modifie ensuite les données des courbes avec set_data, ce qui revient à leur passer le nouveau vecteur des x (entiers de 0 à p) et des y (probabilités pour p ou cumulées jusque p) . On modifie aussi le titre et l'étendue de l'axe des x.
|
|
|
|
fig.canvas.draw_idle() force les objets modifiés à se retracer.
|
|
|
|
|
|
|
|
:warning: les objets tracés avec **plot** sont très faciles à mettre à jour avex **set_data**. Ce n'est malheureusement pas le cas de tous les types de graphiques. Dans certains cas, il sera même peut-être plus simple d'effacer la courbe (on doit donc avoir la permission de le faire, par exemple avec global) et de la retracer.
|
|
|
|
:warning: les objets tracés avec **plot** sont très faciles à mettre à jour avec **set_data**. Ce n'est malheureusement pas le cas de tous les types de graphiques. Dans certains cas, il sera même peut-être plus simple d'effacer la courbe (on doit donc avoir la permission de le faire, par exemple avec global) et de la retracer.
|
|
|
|
|
|
|
|
#### remarque sur la portée des variables
|
|
|
|
Ni fig ni ax ni les courbes n'ont été passés en argument, pourtant ils sont connus depuis l'intérieur de la fonction. Comme c'est un petit code, on sait que l'on va manipuler les bons objets, il n'y a pas d'ambiguité. Dans la fonction, on ne pourrait pas écraser l'un ou l'autre objet, par contre si on applique sur eux une méthode (suptitle, set_data, set_xlim...) la modification de leurs attributs correspondant est répercutée.
|
|
|
|
Ni **fig** ni **ax** ni les courbes n'ont été passés en argument, pourtant ils sont connus depuis l'intérieur de la fonction. Comme c'est un petit code, on sait que l'on va manipuler les bons objets, il n'y a pas d'ambiguité. Dans la fonction, on ne pourrait pas écraser l'un ou l'autre objet, par contre si **on applique sur eux une méthode** (suptitle, set_data, set_xlim...) la modification de leurs attributs correspondant est répercutée.
|
|
|
|
|
|
|
|
## Variante avec une courbe "en escalier"
|
|
|
|
|
|
|
|
On a réécrit un code mais en gardant la possibilité de représenter la courbe soit avec **plot** soit avec un escalier en utilisant **hlines** et **fillbetween**.
|
|
|
|
On a réécrit un code avec la possibilité de représenter la courbe soit avec **plot** soit avec un escalier en utilisant **hlines** et **fillbetween**.
|
|
|
|
|
|
|
|
[exemple de graphique en escalier](https://gitlab.irstea.fr/christine.poulard/atelier-matplotlib/-/blob/master/ScE_-_Hydrologie_-_D%C3%A9mo_probas_crue_sur_N_ann%C3%A9es.png)
|
|
|
|
Pour passer d'une représentation à l'autre, il faut changer la valeur du booléen TRACE_PAR_PLOT (on rappelle que les valeurs du booléen sont False et True, avec une majuscule).
|
| ... | ... | @@ -223,8 +222,8 @@ Au passage, on en profite pour choisir quelle fonction update_graphique sera li |
|
|
|
|
|
|
|
``` python
|
|
|
|
if TRACE_PAR_PLOT:
|
|
|
|
courbe_exactement, = ax.plot(range(p + 1), liste_probas_exactement, label="exactement p crues",ls='None', marker='o', markersize=10)
|
|
|
|
courbe_cumul, = ax.plot(range(p + 1), liste_probas_auplus, label="p crues ou moins de p", ls='None', marker='o', markersize=10)
|
|
|
|
courbe_exactement, = ax.plot(range(k + 1), liste_probas_exactement, label="exactement k crues",ls='None', marker='o', markersize=10)
|
|
|
|
courbe_cumul, = ax.plot(range(k + 1), liste_probas_auplus, label="k crues ou moins de k", ls='None', marker='o', markersize=10)
|
|
|
|
update_graphique = update_graphique_plot
|
|
|
|
else:
|
|
|
|
stair_plot(range(p + 1),liste_probas_exactement, couleur="blue")
|
| ... | ... | @@ -235,10 +234,10 @@ else: |
|
|
|
### fontion update, cas graphique_stairs
|
|
|
|
|
|
|
|
Les calculs sont les mêmes, mais au lieu de modifier une courbe on efface le contenu de la vignette avec "clear" et on retrace.
|
|
|
|
Plus intéressant : on veut que les étiquettes de p soient au milieu des classes. Pour cela, on utilise les Locator et Formattor :
|
|
|
|
Plus intéressant : on veut que les étiquettes de k soient au milieu des classes. Pour cela, on utilise les Locator et Formattor :
|
|
|
|
- si p> 15 on place les ticks de 5 en 5 en commençant par 0.5 et les ticks secondaires de 1 en 1 depuis 0.5
|
|
|
|
- si p<= 15 on place les ticks de 1 en 1 en commençant par 0.5 et on enlève les ticks secondaires
|
|
|
|
- évidemment, on ne veut pas que les étiquettes soient "0.5, 1,5..." ; donc on les formatte au format entier.
|
|
|
|
- évidemment, on ne veut pas que les étiquettes soient "0.5, 1,5..." ; donc on les formate au format entier.
|
|
|
|
|
|
|
|
|
|
|
|
``` python
|
| ... | ... | @@ -247,19 +246,19 @@ def update_graphique_stairs(val): |
|
|
|
n = int(slider_n.val)
|
|
|
|
T = int(slider_T.val)
|
|
|
|
|
|
|
|
p, liste_probas_exactement, liste_probas_auplus = calcul_courbes(n, 1 / T)
|
|
|
|
k, liste_probas_exactement, liste_probas_auplus = calcul_courbes(n, 1 / T)
|
|
|
|
fig.suptitle(f"Démo : probabilité d'avoir p crues > T={T}ans en {n} années")
|
|
|
|
|
|
|
|
ax.clear()
|
|
|
|
stair_plot(range(p + 1), liste_probas_exactement, couleur="blue")
|
|
|
|
stair_plot(range(p + 1), liste_probas_auplus, couleur="sienna")
|
|
|
|
stair_plot(range(k + 1), liste_probas_exactement, couleur="blue")
|
|
|
|
stair_plot(range(k + 1), liste_probas_auplus, couleur="sienna")
|
|
|
|
|
|
|
|
ax.set_xlim(0, p + 1)
|
|
|
|
ax.set_xlim(0, k + 1)
|
|
|
|
if p > 15:
|
|
|
|
ax.xaxis.set_major_locator(plt.FixedLocator([0.5 + i for i in range(0, p, 5)]))
|
|
|
|
ax.xaxis.set_minor_locator(plt.FixedLocator([0.5 + i for i in range(0, p, 1)]))
|
|
|
|
ax.xaxis.set_major_locator(plt.FixedLocator([0.5 + i for i in range(0, k, 5)]))
|
|
|
|
ax.xaxis.set_minor_locator(plt.FixedLocator([0.5 + i for i in range(0, k, 1)]))
|
|
|
|
else:
|
|
|
|
ax.xaxis.set_major_locator(plt.FixedLocator([0.5 + i for i in range(0, p, 1)]))
|
|
|
|
ax.xaxis.set_major_locator(plt.FixedLocator([0.5 + i for i in range(0, k, 1)]))
|
|
|
|
ax.xaxis.set_minor_locator(plt.NullLocator())
|
|
|
|
ax.xaxis.set_major_formatter(ticker.FormatStrFormatter("%d"))
|
|
|
|
fig.canvas.draw_idle()
|
| ... | ... | @@ -287,7 +286,7 @@ def T_emp_parametre(n_annees, a , b): |
|
|
|
# instructions qui affichent des résultats dans la console
|
|
|
|
# 1) version simple, résultats numériques non formattés = trop de chiffres après la virgule
|
|
|
|
print("Paramétré par curseur ", plotting_positions_param)
|
|
|
|
# 2) version formatée avec des "f strings": on rappelle a, b et N, a et b étant formattées comme des flottants avec 1 chiffre après la virgule ; on crée une liste de valeurs formattées par {freq:.2f} et on joint les strings avec la chaîne " / "
|
|
|
|
# 2) version formatée avec des "f strings": on rappelle a, b et N, a et b étant formatées comme des flottants avec 1 chiffre après la virgule ; on crée une liste de valeurs formatées par {freq:.2f} et on joint les strings avec la chaîne " / "
|
|
|
|
print(f"Périodes de retour empiriques, a={a:.1f}, b={b:.1f}, N={N} : ", " / ".join([f"{freq:.2f}" for freq in plotting_positions_param]))
|
|
|
|
|
|
|
|
return plotting_positions_param
|
| ... | ... | |
