Bayesian linear regression Vajont

 Nello stesso articolo sono presentati anche i dati della frana del Vajont avvenuta nella notte edl 9 ottobre 1963 (digitalizzando il grafico il punto risulta essere alle 22:30 del 9/10/1963, l'orario esatto e' le 22:39)

Le previsioni risultano essere

Metodo minimi quadrati : 07/10/1963 14:30 circa

Bayesian regression :  07/10/1963 15:10 circa

 In questo caso i due metodi sono sostanzialmente sovrapponibili e sbagliano entrambi la data dell'evento

 

Grafico originale dell'articolo

 grafici originali dell'invaso del Vajont. Si vede che le condizioni al contorno non sono stabili e questo genera i flessi nella serie di misure

 

 

Ascissa dell'intercetta con l'asse x: x = 2438310.10911
Errore stimato per l'ascissa dell'intercetta: σ_x = 83814.68063
Valore di RMSE: 0.00745
Valore di R^2: 0.99006
Coefficiente angolare (m): -0.00372
Intercetta (b): 9082.54821
Coefficiente angolare (m): -0.00372 ± 0.00009
Intercetta (b): 9082.54821 ± 220.76079

 

 

Mean Squared Error (MSE): 2.22082811850103e-06
R^2 Score: 0.999641381986245
Coefficiente: [-0.003733]
Intercetta: 9102.213990030079
Intercetta sull'asse x: 2438310.1325 ± 2608880.4621
Coefficiente angolare (m): -0.0037 ± 0.0040
Intercetta sull'asse y (b): 9102.2140 ± 0.0087 

se si prendono in considerazione solo gli ultimi 6 punti si ha che la stima dell'evento e'

Metodo minimi quadrati : 10/10/1963 05:30 circa

Bayesian regression : 10/10/1963 12:00 circa

 

 

 

Mean Squared Error (MSE): 1.1478977160850574e-06
R^2 Score: 0.8087633959041971
Coefficiente: [-0.00254709]
Intercetta: 6210.613078964549
Intercetta sull'asse x: 2438312.9889 ± 2810932.1096
Coefficiente angolare (m): -0.0025 ± 0.0029
Intercetta sull'asse y (b): 6210.6131 ± 0.0008

Ascissa dell'intercetta con l'asse x: x = 2438312.73126
Errore stimato per l'ascissa dell'intercetta: σ_x = 203901.17932
Valore di RMSE: 0.00055
Valore di R^2: 0.98621
Coefficiente angolare (m): -0.00271
Intercetta (b): 6599.21777
Coefficiente angolare (m): -0.00271 ± 0.00016
Intercetta (b): 6599.21777 ± 390.21805

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import BayesianRidge
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import train_test_split

# Crea un DataFrame dai dati
data = {
    "Tempo": [

2438249.6968,
2438254.7664,
2438259.6571,
2438264.6074,
2438270.0944,
2438274.3290,
2438280.2336,
2438285.4225,
2438290.1342,
2438295.3827,
2438300.0348,
2438303.3151,
2438304.1501,
2438305.2833,
2438306.4165,
2438307.4304,
2438308.2654,
2438309.2793,
2438310.4125

    ],
    "1/V": [

0.2248,
0.2060,
0.1903,
0.1770,
0.1649,
0.1310,
0.1038,
0.0802,
0.0591,
0.0452,
0.0300,
0.0246,
0.0228,
0.0210,
0.0161,
0.0143,
0.0119,
0.0095,
0.0065


    ]
}
df = pd.DataFrame(data)

# Variabili indipendente (Tempo) e dipendente (1/V)
X = df[["Tempo"]]
y = df["1/V"]

# Dividi il dataset in training e test
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Inizializza il modello Bayesian Ridge
model = BayesianRidge()

# Addestra il modello
model.fit(X_train, y_train)

# Effettua previsioni
y_pred = model.predict(X_test)

# Calcola metriche
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)

# Stampa i risultati
print("Mean Squared Error (MSE):", mse)
print("R^2 Score:", r2)
print("Coefficiente:", model.coef_)
print("Intercetta:", model.intercept_)

# Incertezze (varianza -> deviazione standard)
sigma_m = np.sqrt(1 / model.lambda_)  # Deviazione standard di m
sigma_b = np.sqrt(1 / model.alpha_)   # Deviazione standard di b

# Calcolo punto di intercetta sull'asse x
m = model.coef_[0]         # Coefficiente angolare (slope)
b = model.intercept_       # Intercetta sull'asse y
x_intercept = -b / m

# Propagazione errore per l'intercetta sull'asse x
sigma_x = np.sqrt((sigma_b / m) ** 2 + (b * sigma_m / m**2) ** 2)

# Risultati
print(f"Intercetta sull'asse x: {x_intercept:.4f} ± {sigma_x:.4f}")
print(f"Coefficiente angolare (m): {m:.4f} ± {sigma_m:.4f}")
print(f"Intercetta sull'asse y (b): {b:.4f} ± {sigma_b:.4f}")



# Visualizza il fit del modello
plt.figure(figsize=(10, 6))
plt.scatter(df["Tempo"], df["1/V"], color="blue", label="Dati originali")
plt.plot(df["Tempo"], model.predict(df[["Tempo"]]), color="red", label="Fit del modello")
intercept_x = -model.intercept_ / model.coef_[0]
plt.axvline(x=intercept_x, color="green", linestyle="--", label=f"Intersezione x={intercept_x:.2f}")
plt.xlabel("Tempo")
plt.ylabel("1/V")
plt.title("Regressione Lineare Bayesiana: Tempo vs 1/V")
plt.legend()
plt.grid()
plt.show()