Bayesian linear regression (2)

Per migliorare il post precedente ho inserito le date in formato giuliano (dati digitalizzati con https://apps.automeris.io/wpd4/)

L'evento e' avvenuto il 28 dicembre 2002 ore 10 corrispondente al tempo giuliano 2452636,91522

Base 1-2

Il metodo ai minimi quadrati ha stimato l'evento per il 24 dicembre 2002 alle 01 (circa)

Il metodo di regressione bayesiana ha stimato l'evento per  25 dicembre 2002 ore 6 (circa)

Base  15-13

Il metodo ai minimi quadrati ha stimato l'evento per il 21 dicembre 2002 alle 05 (circa)

Il metodo di regressione bayesiana ha stimato l'evento per  23 dicembre 2002 ore 18 (circa)

Al di la' di quale sia il metodo migliore e' da mettere in evidenzia il valore di σ_x

 

Base 1-2

Metodo minimi quadrati

Ascissa dell'intercetta con l'asse x: x = 2452632.56444
Errore stimato per l'ascissa dell'intercetta: σ_x = 136727.00244
Valore di RMSE: 0.07491
Valore di R^2: 0.97426
Coefficiente angolare (m): -0.01137
Intercetta (b): 27882.92583
Coefficiente angolare (m): -0.01137 ± 0.00045
Intercetta (b): 27882.92583 ± 1099.10446

 

Regressione lineare bayesiana

Mean Squared Error (MSE): 0.009615241747958183
R^2 Score: 0.9475633566354236
Coefficiente: [-0.01111325]
Intercetta: 27256.72908552682
Intercetta sull'asse x: 2452633.7423 ± 2474970.5684
Coefficiente angolare (m): -0.0111 ± 0.0112
Intercetta sull'asse y (b): 27256.7291 ± 0.0710

Base 15-13

Metodo minimi quadrati

Ascissa dell'intercetta con l'asse x: x = 2452629.70635
Errore stimato per l'ascissa dell'intercetta: σ_x = 290006.65020
Valore di RMSE: 0.03778
Valore di R^2: 0.88275
Coefficiente angolare (m): -0.00254
Intercetta (b): 6217.58752
Coefficiente angolare (m): -0.00254 ± 0.00021
Intercetta (b): 6217.58752 ± 519.84916

 

Regressione lineare bayesiana

Mean Squared Error (MSE): 0.002548090926111798
R^2 Score: 0.8769276021004735
Coefficiente: [-0.00238676]
Intercetta: 5853.843312339611
Intercetta sull'asse x: 2452632.2719 ± 2861671.0348
Coefficiente angolare (m): -0.0024 ± 0.0028
Intercetta sull'asse y (b): 5853.8433 ± 0.0348

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": [
        2452489.13,2452498.78,2452508.23,2452521.87,2452530.69,2452537.62,2452548.53,2452549.79,2452555.88,
    2452562.81,2452568.69,2452577.50,2452585.48,2452590.73,2452599.75,2452605.63,2452610.67,2452622.64,2452628.72
    ],
    "1/V": [
        1.71,1.43,1.45,1.27,1.14,1.2,1.01,0.88,0.89,0.7,0.55,0.55,0.62,0.55,0.4,0.26,0.21,0.13,0.12
    ]
}
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()