Riprendendo il discorso abbandonato qui , ho ripreso i dati ed applicato il bayesian ridge
Questa implementazione e' interessante perche' fornisce anche il calcolo dell'incertezza del valore dell'intercetta della retta di regressione con l'asse x
Bayesian Ridge
Mean Squared Error (MSE): 0.008041457036234597
R^2 Score: 0.664885642861212
Intercetta sull'asse x: 304.9129 ± 497.0633
Coefficiente angolare (m): -0.0011 ± 0.0018
Intercetta sull'asse y (b): 0.3368 ± 0.0341
Per confronto gli stessi dati sono stati elaborati con la regressione lineare ai minimi quadrati
Regressione a minimi quadrati
Ascissa dell'intercetta con l'asse x: x = 292.05703
Errore stimato per l'ascissa dell'intercetta: σ_x = 30.98612
Valore di RMSE: 0.04669
Valore di R^2: 0.86116
Coefficiente angolare (m): -0.00131 ± 0.00012
Intercetta (b): 0.38237 ± 0.02161
----------------------------------------------------------
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": [14, 28, 42, 56, 70, 84, 98, 112, 126, 140, 154, 168, 182, 196, 210, 224, 238, 252, 266, 280, 294, 308],
"1/V": [
0.4586046511627907, 0.44000000000000006, 0.30604651162790697, 0.27813953488372095,
0.28930232558139535, 0.3227906976744187, 0.21860465116279076, 0.17023255813953395,
0.17953488372092963, 0.1367441860465109, 0.13302325581395344, 0.11441860465116252,
0.08837209302325501, 0.0958139534883719, 0.14046511627906993, 0.13860465116279058,
0.09581395348837188, 0.04744186046511625, 0.051162790697674376, 0.02883720930232564,
0.030697674418604465, 0.010232558139534276
]
}
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()
-----------------------------------------------------------
questo il codice per la regressione lineare
-----------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
# Dataset fornito
data = {
"Tempo": [
14, 28, 42, 56, 70, 84, 98, 112, 126, 140, 154, 168, 182, 196, 210,
224, 238, 252, 266, 280, 294, 308
],
"1/V": [
0.4586046511627907, 0.44000000000000006, 0.30604651162790697,
0.27813953488372095, 0.28930232558139535, 0.3227906976744187,
0.21860465116279076, 0.17023255813953395, 0.17953488372092963,
0.1367441860465109, 0.13302325581395344, 0.11441860465116252,
0.08837209302325501, 0.0958139534883719, 0.14046511627906993,
0.13860465116279058, 0.09581395348837188, 0.04744186046511625,
0.051162790697674376, 0.02883720930232564, 0.030697674418604465,
0.010232558139534276
]
}
# Estrazione dei dati
x = np.array(data["Tempo"])
y = np.array(data["1/V"])
# Calcolo dei coefficienti della regressione lineare
x_mean = np.mean(x)
y_mean = np.mean(y)
n = len(x)
numerator = np.sum((x - x_mean) * (y - y_mean))
denominator = np.sum((x - x_mean)**2)
m = numerator / denominator
b = y_mean - m * x_mean
# Calcolo dei residui e stima degli errori
residuals = y - (m * x + b)
sigma_squared = np.sum(residuals**2) / (n - 2)
sigma_m = np.sqrt(sigma_squared / denominator)
sigma_b = np.sqrt(sigma_squared * np.sum(x**2) / (n * denominator))
# Calcolo dell'ascissa dell'intercetta con l'asse x e del suo errore
x_intercept = -b / m
sigma_x_intercept = np.sqrt((sigma_b / m)**2 + (b * sigma_m / m**2)**2)
# Output dei risultati
print(f"Ascissa dell'intercetta con l'asse x: x = {x_intercept:.5f}")
print(f"Errore stimato per l'ascissa dell'intercetta: σ_x = {sigma_x_intercept:.5f}")
y_pred = m * x + b
rmse = np.sqrt(np.mean((y - y_pred)**2))
# Output del risultato
print(f"Valore di RMSE: {rmse:.5f}")
# Calcolo del valore di R^2
ss_total = np.sum((y - y_mean)**2)
ss_residual = np.sum((y - y_pred)**2)
r_squared = 1 - (ss_residual / ss_total)
# Output del risultato
print(f"Valore di R^2: {r_squared:.5f}")
print(f"Coefficiente angolare (m): {m:.5f}")
print(f"Intercetta (b): {b:.5f}")
residuals = y - y_pred
sigma_squared = np.sum(residuals**2) / (n - 2)
# Errore standard per m
sigma_m = np.sqrt(sigma_squared / denominator)
# Errore standard per b
sigma_b = np.sqrt(sigma_squared * np.sum(x**2) / (n * denominator))
# Stampa dei valori con le incertezze
print(f"Coefficiente angolare (m): {m:.5f} ± {sigma_m:.5f}")
print(f"Intercetta (b): {b:.5f} ± {sigma_b:.5f}")
plt.figure(figsize=(10, 6))
plt.scatter(x, y, color='blue', label='Dati originali')
plt.plot(x, y_pred, color='red', label='Retta di regressione')
plt.xlabel('Tempo')
plt.ylabel('1/V')
plt.title('Regressione Lineare ai Minimi Quadrati')
plt.legend()
plt.grid()
plt.show()
