Pettitt text
import numpy as np
import pandas as pd
import [Link] as plt
def pettitt_test(series):
"""
Test de Pettitt non paramétrique pour détecter une rupture dans une série temporelle.
Retourne la statistique K, la p-value approximative, la position tau et la série U(t).
"""
series = [Link](series)
n = len(series)
U = [Link](n)
for t in range(n):
somme = 0
for i in range(t + 1):
for j in range(t + 1, n):
if series[i] > series[j]:
somme += 1
elif series[i] < series[j]:
somme -= 1
U[t] = somme
K = [Link]([Link](U))
tau = [Link]([Link](U))
# Approximation de la p-value (validée pour n >= 10)
p_value = 2 * [Link]((-6 * (K**2)) / (n**3 + n**2))
return K, p_value, tau, U
# 🔹 Chargement des données
df = pd.read_csv(r"C:\Users\MON PC\Desktop\traitement\
Q95_annuel_hanning_smoothed.csv") # ← fichier corrigé
# 🔹 Vérification et conversion
years = df["date"].astype(int).values # Les dates sont des entiers : 1993, 1994, ...
precip = df["precip"].astype(float).values # Valeurs numériques
# 🔹 Application du test de Pettitt
K, p_value, tau, U = pettitt_test(precip)
n = len(precip)
# 🔹 Paramètres
alpha = 0.05 # Seuil de signification classique (5%)
rupture = p_value < alpha
date_tau = years[tau]
# 🔹 Affichage des résultats
print(f"Statistique K : {K:.2f}")
print(f"P-value : {p_value:.4f}")
print(f"Seuil α : {alpha}")
print(f"Position τ : {tau} (année : {date_tau})")
print(f"Nombre de données n : {n}")
if rupture:
print("\n✅ Cas 1 : Rupture significative détectée.")
print(f"➡ Rupture estimée en {date_tau}.")
else:
print("\n❌ Cas 2 : Pas de rupture significative au seuil α = 0.05.")
# 🔹 Graphique 1 : U(t)
[Link](figsize=(11, 4.5))
[Link](years, U, marker="o", ms=4, lw=1.2, color="black", alpha=0.9)
[Link](0, linestyle=":", linewidth=1.2, color="gray")
if rupture:
[Link](date_tau, linestyle="--", linewidth=1.6, color="black")
[Link](f"τ = {date_tau}\nK={K:.0f}, p={p_value:.3g}",
xy=(date_tau, U[tau]),
xytext=(10, 15), textcoords="offset points",
bbox=dict(boxstyle="round", fc="white", alpha=0.8))
[Link]("Évolution de la statistique U(t) du test de Pettitt")
[Link]("Années")
[Link]("U(t)")
[Link](alpha=0.5)
[Link](years, rotation=45)
plt.tight_layout()
[Link]()
# 🔹 Graphique 2 : Précipitations
[Link](figsize=(11, 4.5))
[Link](years, precip, marker="o", ms=4, lw=1.2, color="darkorange", alpha=0.9)
if rupture:
[Link](date_tau, linestyle="--", linewidth=1.8, color="darkorange")
mean1 = [Link](precip[:tau])
mean2 = [Link](precip[tau:]) # Note : inclut la valeur à tau dans le second segment ou
non selon convention
[Link](mean1, xmin=years[0], xmax=date_tau,
colors="red", linestyles="--", linewidth=1.5, label=f"Moyenne avant = {mean1:.2f}")
[Link](mean2, xmin=date_tau, xmax=years[-1],
colors="green", linestyles="--", linewidth=1.5, label=f"Moyenne après =
{mean2:.2f}")
[Link](f"Rupture: {date_tau}\nK={K:.0f}, p={p_value:.3g}",
xy=(date_tau, [Link](precip)),
xytext=(10, -25), textcoords="offset points",
bbox=dict(boxstyle="round", fc="white", alpha=0.8))
[Link]("Précipitations minimales annuelles")
[Link]("Années")
[Link]("Précipitations (unité non précisée)")
[Link]()
[Link](alpha=0.5)
[Link](years, rotation=45)
plt.tight_layout()
[Link]()
Kendall test
import pandas as pd
import pymannkendall as mk
import numpy as np
import [Link] as plt
# Charger les données
df = pd.read_csv(r"C:\Users\MON PC\Desktop\traitement\minmax\
min_abo_annuel_hanning_smoothed.csv")
df = df.sort_values("date")
# Test de Mann-Kendall
result = mk.original_test(df['precip'])
print(result)
alpha = 0.05
print("Statistique U_MK :", result.z)
print("p-value :", result.p)
print("Seuil α :", alpha)
if result.p > alpha:
print("\nCas 1 : On accepte H₀ → Absence de tendance significative.")
else:
if result.z > 0:
print("\nCas 2 : On rejette H₀ → Tendance à la hausse (U_MK positif).")
elif result.z < 0:
print("\nCas 2 : On rejette H₀ → Tendance à la baisse (U_MK négatif).")
else:
print("\nCas 2 : On rejette H₀ → Tendance neutre (U_MK ≈ 0).")
# =================== VISUALISATION ===================
# Préparer l’axe temps
if [Link](df["date"].dtype, [Link]) or [Link](df["date"].dtype, [Link]):
years = df["date"].astype(int).to_numpy()
else:
dparse = pd.to_datetime(df["date"], errors="coerce")
years = ([Link](pd.to_numeric(df["date"],
errors="coerce"))).astype(int).to_numpy()
y = df["precip"].to_numpy()
x = [Link](float)
# -------- Figure 1 : Série + Sen’s slope --------
sen_slope = getattr(result, "slope", [Link])
if [Link](sen_slope):
intercept = [Link](y - sen_slope * x)
y_hat = sen_slope * x + intercept
else:
coeffs = [Link](x, y, 1)
sen_slope, intercept = coeffs[0], coeffs[1]
y_hat = sen_slope * x + intercept
sig_text = "significative" if result.p <= alpha else "non significative"
# -------- Fonction UF/UB --------
def mk_progressive_curves(values):
z = [Link](values, dtype=float)
n = len(z)
# UF
S = [Link](n)
for i in range(1, n):
S[i] = S[i-1] + [Link]([Link](z[i] - z[:i]))
i = [Link](n)
EF = i*(i-1)/4.0
VF = i*(i-1)*(2*i+5)/72.0
UF = [Link](n); mask = VF > 0
UF[mask] = (S[mask] - EF[mask]) / [Link](VF[mask])
# UB
Sr = [Link](n)
zr = z[::-1]
for i in range(1, n):
Sr[i] = Sr[i-1] + [Link]([Link](zr[i] - zr[:i]))
ER = i*(i-1)/4.0
VR = i*(i-1)*(2*i+5)/72.0
UR = [Link](n); maskr = VR > 0
UR[maskr] = (Sr[maskr] - ER) / [Link](VR)
UB = -UR[::-1]
return UF, UB
UF, UB = mk_progressive_curves(y)
# Détection de la date de tendance
trend_date = None
for i in range(1, len(years)):
if [Link](UF[i] - UB[i]) != [Link](UF[i-1] - UB[i-1]):
trend_date = years[i]
break
# Si pas d’intersection trouvée → prendre max(|UF-UB|)
if trend_date is None and result.p <= alpha:
idx = [Link]([Link](UF - UB))
trend_date = years[idx]
# --------- Premier graphique ---------
[Link](figsize=(11,4.6))
[Link](years, y, marker="o", lw=1.2, ms=4, color="tab:blue", alpha=0.9, label="Série")
[Link](years, y_hat, linestyle="--", lw=1.6, color="tab:orange",
label=f"Tendance Sen (slope={sen_slope:.3g})")
if trend_date and result.p <= alpha:
# Ligne verticale à la date de tendance
[Link](trend_date, linestyle="--", color="black")
# Moyennes avant et après
mean_before = [Link](y[years < trend_date])
mean_after = [Link](y[years >= trend_date])
[Link](mean_before, xmin=[Link](), xmax=trend_date,
colors="red", linestyles="--", linewidth=1.5, label=f"Moyenne avant
({mean_before:.1f})")
[Link](mean_after, xmin=trend_date, xmax=[Link](),
colors="green", linestyles="--", linewidth=1.5, label=f"Moyenne après
({mean_after:.1f})")
# Annotation
[Link](f"Début tendance : {trend_date}",
xy=(trend_date, [Link](y)),
xytext=(10, -30), textcoords="offset points",
bbox=dict(boxstyle="round", fc="white", alpha=0.8))
print(f"📌 La tendance significative apparaît autour de {trend_date}.")
print(f" Moyenne avant : {mean_before:.2f}")
print(f" Moyenne après : {mean_after:.2f}")
[Link](f"Mann–Kendall : tendance {[Link]} ({sig_text}), p={result.p:.3g}")
[Link]("Année"); [Link]("Précipitations")
[Link](alpha=0.25); [Link]()
plt.tight_layout(); [Link]()
# --------- Second graphique : UF/UB ---------
[Link](figsize=(11,4.6))
[Link](years, UF, lw=1.3, color="tab:green", label="UF (progressif)")
[Link](years, UB, lw=1.3, color="tab:red", label="UB (régressif)")
[Link](0, color="gray", lw=1, linestyle=":")
[Link](1.96, color="gray", lw=1, linestyle="--")
[Link](-1.96, color="gray", lw=1, linestyle="--")
if trend_date and result.p <= alpha:
[Link](trend_date, linestyle="--", color="black")
[Link](f"Début tendance : {trend_date}",
xy=(trend_date, 0), xytext=(10,10), textcoords="offset points",
bbox=dict(boxstyle="round", fc="white", alpha=0.8))
[Link]("Mann–Kendall : courbes UF / UB (α ≈ 5%)")
[Link]("Année"); [Link]("Statistique normalisée")
[Link](alpha=0.25); [Link]()
plt.tight_layout(); [Link]()
SPI
import pandas as pd
from [Link] import gamma, norm
import [Link] as plt
# Charger le fichier CSV contenant les précipitations annuelles
chemin_fichier = (r"C:\Users\MON PC\Desktop\traitement\moyenne_lanta_annuel.csv")
data = pd.read_csv(chemin_fichier)
data = [Link]()
data = data.sort_values('date')
# Ajouter une colonne 'Année' si elle n'existe pas
if 'Année' not in [Link]:
data['date'] = pd.to_datetime(data['date'])
data['Année'] = data['date'].[Link]
# Regrouper les précipitations par année
annual_precip = [Link]('Année')['precip'].sum()
# Calculer le SPI annuel
alpha, loc, beta = [Link](annual_precip, floc=0)
cdf = [Link](annual_precip, a=alpha, loc=loc, scale=beta)
spi_annual = [Link](cdf)
# Créer un DataFrame pour le SPI annuel
spi_annual_df = [Link]({
'Année': annual_precip.index,
'SPI_Annuel': spi_annual
})
# Affichage : histogramme SPI annuel avec ligne à SPI=0
[Link](figsize=(10, 5))
[Link](spi_annual_df['Année'], spi_annual_df['SPI_Annuel'], color='#3498db')
[Link](0, color='black', linestyle='--', linewidth=1)
[Link]('Année')
[Link]('SFI Annuel')
[Link]('SFI LANTA')
[Link](spi_annual_df['Année'], rotation=45)
plt.tight_layout()
[Link]()