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()

Bayesian Ridge

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()

Bilateral LSTM e GRU con dati di estensimetro

Avevo gia' provato ad usare LSTM per dati derivanti da un estensimetro qui ... adesso ci riprovo con la variante Bilateral LSTM e GRU prendendo spunto da questa pagina 

Questi sono i dati raw. I dati sono campionati ogni 20 minuti e sono presenti diverse fasi di un movimento di versante

Prima di iniziare un paio di considerazioni

In questo grafico sono riassunti gli spostamenti registrati dall'estensimetro, la velocita' di movimento ed i dati pluviometrici

Dalla cross correlazione tra il dato meteo e quello di spostamento si un lag di 27 giorni

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

df_b = pd.read_csv('est2.csv', sep=',',header=0, parse_dates = ['Date'], date_format = '%Y-%m-%d %H:%M:%S')

pluvio = pd.read_csv('arpa.csv', sep=';',header=0, parse_dates = ['Date'], date_format = '%d/%m/%Y %H:%M:%S')

df_b_res = df_b.resample('24h', on='Date').mean()
pluvio_res = pluvio.resample('24h', on='Date').mean()

print(df_b_res)
print(pluvio)

df_b_res['Diff'] = df_b_res['Est'].diff()
df_b_res['Inv'] = 1/df_b_res['Diff']

time = pd.date_range('2024-02-01', periods=214, freq='D')
print(time)
print(pluvio_res.iloc[:,0])
print(df_b_res.iloc[:,0])

df = pd.DataFrame({'series1': df_b_res.iloc[:,0], 'series2': pluvio_res.iloc[:,0]}, index=time)
cross_corr = df['series1'].corr(df['series2'])
def cross_correlation(series1, series2, max_lag):
    corr = []
    for lag in range(-max_lag, max_lag + 1):
        corr.append(df['series1'].shift(lag).corr(df['series2']))
    return corr
max_lag = 40  # You can adjust the maximum lag
lags = range(-max_lag, max_lag + 1)
correlation_values = cross_correlation(df['series1'], df['series2'], max_lag)

# Find the lag that maximizes the correlation
best_lag = lags[np.argmax(correlation_values)]
print(f"The best lag is: {best_lag} days")

#print(df_b_res)

plt.figure(figsize=(11,7))
plt.plot(df_b_res['Est'],"-b",label="Est. (0.1 mm)")
plt.plot(df_b_res['Diff']*35,"-r",label="35*Vel. (0.1 mm/day)")
plt.plot(pluvio_res['Prec'],"-g",label="Rain (mm)")


plt.title('Velocita')
plt.xlabel('Date')
plt.ylabel('Estensimetro 0.1mm')
plt.legend(loc="lower left")

plt.show()

plt.figure(figsize=(11,7))
plt.plot(df_b_res['Inv'],"+b",label="Est.")
plt.title('1/Velocita')
plt.xlabel('Date')
plt.ylabel('1/V day/0.1mm')
plt.legend(loc="lower left")

plt.show()

Un paio di considerazioni

1) LSTM e' molto sensibile all'overfitting del modello.Nel caso specifico valori di 5-6 epochs sono gia' ottimali per la convergenza del modello

2) se il modello non converge la predizione magari ha una forma e' corretta ma e' shiftata (tipo il grafico sottostante

3) si deve giocare sul valore della finestra sia come units che come time_steps (oltre che sulle epochs) per ottenere un modello ottimale

import numpy as np

import pandas as pd

import tensorflow as tf

from tensorflow.keras.models import Sequential

from tensorflow.keras.layers import LSTM, Dense, Bidirectional

from sklearn.preprocessing import MinMaxScaler

import matplotlib.pyplot as plt

data = pd.read_csv('est2.csv')

data = data['Est'].values.reshape(-1, 1)

scaler = MinMaxScaler(feature_range=(0, 1))

data_scaled = scaler.fit_transform(data)

def create_dataset(data, time_step=1):

X, y = [], []

for i in range(len(data) - time_step):

X.append(data[i:i + time_step, 0])

y.append(data[i + time_step, 0])

return np.array(X), np.array(y)

time_step = 200

X, y = create_dataset(data_scaled, time_step)

X = X.reshape(X.shape[0], X.shape[1], 1)

train_size = int(len(X) * 0.8)

X_train, X_test = X[:train_size], X[train_size:]

y_train, y_test = y[:train_size], y[train_size:]

model = Sequential()

model.add(Bidirectional(LSTM(units=200, return_sequences=True), input_shape=(time_step, 1)))

model.add(Bidirectional(LSTM(units=200)))

model.add(Dense(units=1))

model.compile(optimizer='adam', loss='mean_squared_error')

model.fit(X_train, y_train, epochs=5, batch_size=32)

predictions = model.predict(X_test)

predictions_rescaled = scaler.inverse_transform(predictions)

y_test_rescaled = scaler.inverse_transform(y_test.reshape(-1, 1))

plt.figure(figsize=(10,6))

plt.plot(y_test_rescaled, label='Actual Data')

plt.plot(predictions_rescaled, label='Predicted Data')

plt.title('Bidirectional LSTM Time Series Forecasting')

plt.xlabel('Time')

plt.ylabel('Value')

plt.legend()

plt.show()

forecast_steps = 1000 # Number of time steps to predict

input_sequence = X_test[-1]

forecast = []

for _ in range(forecast_steps):

predicted_value = model.predict(input_sequence.reshape(1, time_step, 1))

forecast.append(predicted_value[0, 0])

input_sequence = np.append(input_sequence[1:], predicted_value, axis=0)

forecast_rescaled = scaler.inverse_transform(np.array(forecast).reshape(-1, 1))

predictions_rescaled = scaler.inverse_transform(predictions)

y_test_rescaled = scaler.inverse_transform(y_test.reshape(-1, 1))

# Plot the results

plt.figure(figsize=(10,6))

plt.plot(forecast_rescaled, label='Forecasted Data')

plt.title('Bidirectional LSTM Time Series Forecasting')

plt.xlabel('Time')

plt.ylabel('Value')

plt.legend()

plt.show()

risultati sul set di validazione

questi sono i dati di forecast

Ho provato a chiedere a ChatGPT che mi ha risposto con questo codice

import numpy as np

import pandas as pd

import tensorflow as tf

from tensorflow.keras.models import Model

from tensorflow.keras.layers import LSTM, Dense, Input, Bidirectional

from sklearn.preprocessing import MinMaxScaler

from sklearn.model_selection import train_test_split

import matplotlib.pyplot as plt

def load_data(file_path):

df = pd.read_csv(file_path)

data = df['Est'].values

scaler = MinMaxScaler(feature_range=(0, 1))

data_scaled = scaler.fit_transform(data.reshape(-1, 1))

return data_scaled, scaler

def create_dataset(data, time_step=200):

X, y = [], []

for i in range(len(data) - time_step):

X.append(data[i:i + time_step, 0])

y.append(data[i + time_step, 0])

X = np.array(X)

y = np.array(y)

X = np.reshape(X, (X.shape[0], X.shape[1], 1))

return X, y

def create_bilateral_lstm_model(input_shape):

inputs = Input(shape=input_shape)

forward_lstm = LSTM(200, return_sequences=False)(inputs)

backward_lstm = LSTM(200, return_sequences=False, go_backwards=True)(inputs)

merged = tf.keras.layers.concatenate([forward_lstm, backward_lstm])

output = Dense(1)(merged)

model = Model(inputs=inputs, outputs=output)

return model

def train_and_forecast(file_path, time_step=200, forecast_days=200):

data_scaled, scaler = load_data(file_path)

X, y = create_dataset(data_scaled, time_step)

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=False)

model = create_bilateral_lstm_model(input_shape=(X_train.shape[1], 1))

model.compile(optimizer='adam', loss='mean_squared_error')

model.fit(X_train, y_train, epochs=5, batch_size=32, validation_data=(X_test, y_test))

y_pred = model.predict(X_test)

y_pred_rescaled = scaler.inverse_transform(y_pred)

y_test_rescaled = scaler.inverse_transform(y_test.reshape(-1, 1))

plt.figure(figsize=(12, 6))

plt.plot(y_test_rescaled, label="True Data")

plt.plot(y_pred_rescaled, label="Predicted Data")

plt.title("True vs Predicted Data (Test Set)")

plt.legend()

plt.show()

last_sequence = data_scaled[-time_step:].reshape(1, time_step, 1)

future_predictions = []

for _ in range(forecast_days):

forecast = model.predict(last_sequence)

future_predictions.append(forecast[0, 0])

last_sequence = np.append(last_sequence[:, 1:, :], forecast.reshape(1, 1, 1), axis=1)

future_predictions_rescaled = scaler.inverse_transform(np.array(future_predictions).reshape(-1, 1))

plt.figure(figsize=(12, 6))

plt.plot(range(len(data_scaled)), scaler.inverse_transform(data_scaled), label="Historical Data")

plt.plot(range(len(data_scaled), len(data_scaled) + forecast_days), future_predictions_rescaled, label="Forecasted Data", linestyle="--")

plt.title(f"Forecasted Data for Next {forecast_days} Steps")

plt.legend()

plt.show()

return future_predictions_rescaled

file_path = 'est2.csv'

forecast = train_and_forecast(file_path, time_step=200, forecast_days=200)

print("Forecasted future values:")

print(forecast)

ingrandimento dei soli dati di forevast

Visto che il modello converge adesso vediamo come si comporta quando si ha un evento inatteso (ovvero l'innescarsi del movimento di versante)..ho tagliato i dati sul primo movimento che corrisponda circa con la coordinata di ascissa 2000

l'addestramento del modello non e' male (anche se non perfetto)

il problema e' che quando entra in modalita' forecast sbaglia completamente indicando un movimento esattamente opposto a quello della frana

Con GRU le cose non sono andate molto bene

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

from sklearn.preprocessing import MinMaxScaler, StandardScaler

import warnings

from scipy import stats

import tensorflow as tf

from tensorflow import keras

from tensorflow.keras import Sequential, layers, callbacks

from tensorflow.keras.layers import Dense, LSTM, Dropout, GRU, Bidirectional

tf.random.set_seed(1234)

print("Num GPUs Available: ", len(tf.config.list_physical_devices('GPU')))

df = pd.read_csv("est2.csv", parse_dates = ["Date"])

# Define a function to draw time_series plot

def timeseries (x_axis, y_axis, x_label):

plt.figure(figsize = (10, 6))

plt.plot(x_axis, y_axis, color ="black")

plt.xlabel(x_label, {'fontsize': 12})

plt.ylabel('Est', {'fontsize': 12})

plt.show()

print(df)

timeseries(df.index, df['Est'], 'Time (day)')

df = df.drop('Date', axis = 1)

train_size = int(len(df)*0.8)

train_data = df.iloc[:train_size]

test_data = df.iloc[train_size:]

scaler = MinMaxScaler().fit(train_data)

train_scaled = scaler.transform(train_data)

test_scaled = scaler.transform(test_data)

def create_dataset (X, look_back = 1):

Xs, ys = [], []

for i in range(len(X)-look_back):

v = X[i:i+look_back]

Xs.append(v)

ys.append(X[i+look_back])

return np.array(Xs), np.array(ys)

LOOK_BACK = 300

X_train, y_train = create_dataset(train_scaled,LOOK_BACK)

X_test, y_test = create_dataset(test_scaled,LOOK_BACK)

# Print data shape

print("X_train.shape: ", X_train.shape)

print("y_train.shape: ", y_train.shape)

print("X_test.shape: ", X_test.shape)

print("y_test.shape: ", y_test.shape)

# Create GRU model

def create_gru(units):

model = Sequential()

# Input layer

model.add(GRU (units = units, return_sequences = True,

input_shape = [X_train.shape[1], X_train.shape[2]]))

model.add(Dropout(0.2))

# Hidden layer

model.add(GRU(units = units))

model.add(Dropout(0.2))

model.add(Dense(units = 1))

#Compile model

model.compile(optimizer='adam',loss='mse')

return model

model_gru = create_gru(300)

def fit_model(model):

early_stop = keras.callbacks.EarlyStopping(monitor = 'val_loss',patience = 10)

history = model.fit(X_train, y_train, epochs = 5,

validation_split = 0.2,

batch_size = 32, shuffle = False,

callbacks = [early_stop])

return history

history_gru = fit_model(model_gru)

y_test = scaler.inverse_transform(y_test)

y_train = scaler.inverse_transform(y_train)

def plot_loss (history, model_name):

plt.figure(figsize = (10, 6))

plt.plot(history.history['loss'])

plt.plot(history.history['val_loss'])

plt.title('Model Train vs Validation Loss for ' + model_name)

plt.ylabel('Loss')

plt.xlabel('epoch')

plt.legend(['Train loss', 'Validation loss'], loc='upper right')

plt.show()

plot_loss (history_gru, 'GRU')

# Make prediction

def prediction(model):

prediction = model.predict(X_test)

prediction = scaler.inverse_transform(prediction)

return prediction

prediction_gru = prediction(model_gru)

# Plot test data vs prediction

def plot_future(prediction, model_name, y_test):

plt.figure(figsize=(10, 6))

range_future = len(prediction)

plt.plot(np.arange(range_future), np.array(y_test), label='Test data')

plt.plot(np.arange(range_future), np.array(prediction),label='Prediction')

plt.title('Test data vs prediction for ' + model_name)

plt.legend(loc='upper left')

plt.xlabel('Time (day)')

plt.ylabel('Est')

plt.show()

plot_future(prediction_gru, 'GRU', y_test)

# bilstm 5 epoch

#model_bilstm.save("bilstm.keras")

model_gru.save("gru.keras")

Filtraggio dati estensimetro con Kalman

 Ho ripreso in mano il discorso filtro Kalman applicato pero' a dati reali di un estensimetro (i dati sono espressi in decimi di millimetro)

 

 

 

 

 

 

 

 

 

import numpy as np
import matplotlib.pyplot as plt


dati = np.genfromtxt('est.csv', delimiter=';')
dist = dati[:,2]
z = dist.T
print(z)
dimensione = dist.shape
#print(dimensione[0])
n_iter = dimensione[0]
Q = 1e-5 # process variance

# allocate space for arrays
xhat=np.zeros(dimensione)      # a posteri estimate of x
P=np.zeros(dimensione)         # a posteri error estimate
xhatminus=np.zeros(dimensione) # a priori estimate of x
Pminus=np.zeros(dimensione)    # a priori error estimate
K=np.zeros(dimensione)         # gain or blending factor
R = 0.1**2 # estimate of measurement variance, change to see effect

x = 23.6
xhat[0] = 23.6
P[0] = 1.0


for k in range(1,n_iter):
    # time update
    xhatminus[k] = xhat[k-1]
    Pminus[k] = P[k-1]+Q

    # measurement update
    K[k] = Pminus[k]/( Pminus[k]+R )
    xhat[k] = xhatminus[k]+K[k]*(z[k]-xhatminus[k])
    P[k] = (1-K[k])*Pminus[k]

sp = np.diff(xhat, prepend=xhat[0])

plt.figure()
plt.plot(z,'k+',label='noisy measurements')
plt.plot(xhat,'b-',label='a posteri estimate')
plt.legend()
plt.title('Kalman filter ', fontweight='bold')
plt.xlabel('Time')
plt.ylabel('Distance')


plt.figure()
valid_iter = range(1,n_iter) # Pminus not valid at step 0
plt.plot(valid_iter,Pminus[valid_iter],label='a priori error estimate')
plt.title('Estimated $\it{\mathbf{a \ priori}}$ error vs. iteration step', fontweight='bold')
plt.xlabel('Iteration')
plt.ylabel('$(Voltage)^2$')
plt.setp(plt.gca(),'ylim',[0,.01])

plt.figure()
plt.plot(sp,'b-')
plt.title('Speed', fontweight='bold')
plt.xlabel('Time')
plt.ylabel('Speed (0.1mm/day)')

plt.show()

 

 

 

Bayesian linear regression

Volevo provare un approccio differente dal solito per la regressione lineare

Il problema e' relativo alla terza fase del creep in frane la cui legge tra il tempo e l'inverso della velocita' e' approssimata a lineare

Questi dati sono ripresi dall' articolo Tommaso Carlà,Emanuele Intrieri,Federico Di Traglia,Teresa Nolesini,Giovanni Gigli,Nicola Casagli "Guidelines on the use of inverse velocity method as a tool for setting alarm thresholds and forecasting landslides and structure collapses" Landslides DOI 10.1007/s10346-016-0731-5 per la frana di Montebeni. La freccia verde indica la data dell'evento reale e la freccia rossa la data stimata

 La terza fase e' stata modellata tramite regressione a minimi quadrati applicando alcuni filtri

Ho provato a trattare gli stessi dati tramite regressione lineare bayesiana con Metropolis-Hastings Sampler

Alla fine i risultati sono molto simili. C'e' da notare che 

1) essendo il modello Montecarlo non deterministico i risultati non sono mai identici tra una run e la successiva. Ciao' comporta anche uno spostamento dell'intercetta della retta con l'asse del tempo spostando il TOF (time of failure)

2) sia la regressione a minimi quadrati che quella bayesiana sono molto sensibili al dataset (anche perche' alla fine ci sono solo una vendita di misure). Basta omettere un solo dato all'inizio che la correlazione e' sensibilmente differente. Questo comporta la difficolta' di definire in modo univoco il passaggio dalla fase 2 alla fase 3 del creep  

I dati sono

Tempo,1/V
14, 0.4586046511627907
28, 0.44000000000000006
42, 0.30604651162790697
56, 0.27813953488372095
70, 0.28930232558139535
84, 0.3227906976744187
98, 0.21860465116279076
112, 0.17023255813953395
126, 0.17953488372092963
140, 0.1367441860465109
154, 0.13302325581395344
168, 0.11441860465116252
182, 0.08837209302325501
196, 0.0958139534883719
210, 0.14046511627906993
224, 0.13860465116279058
238, 0.09581395348837188
252, 0.04744186046511625
266, 0.051162790697674376
280, 0.02883720930232564
294, 0.030697674418604465
308, 0.010232558139534276

Per il codice sottostante ho modificato un esempio preso da chatgpt

import sys

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


# Load CSV data
def load_data(csv_path):
    data = pd.read_csv(csv_path)
    return data


# Define the likelihood function (Gaussian likelihood)
def likelihood(y, X, beta):
    predictions = np.dot(X, beta)
    return -0.5 * np.sum((y - predictions) ** 2)


# Define the prior function (Normal prior)
def prior(beta, sigma=10):
    return -0.5 * np.sum(beta ** 2) / (sigma ** 2)


# Define the proposal distribution (random walk)
def proposal(beta, step_size=0.1):
    return beta + np.random.normal(0, step_size, size=beta.shape)


# Metropolis-Hastings sampler
def metropolis_hastings(X, y, initial_beta, n_samples, step_size=0.1):
    beta = initial_beta
    samples = []

    for _ in range(n_samples):
        # Propose new beta
        beta_new = proposal(beta, step_size)

        # Compute likelihood and prior for current and proposed beta
        current_likelihood = likelihood(y, X, beta) + prior(beta)
        new_likelihood = likelihood(y, X, beta_new) + prior(beta_new)

        # Metropolis-Hastings acceptance step
        acceptance_ratio = min(1, np.exp(new_likelihood - current_likelihood))

        # Accept or reject the proposed beta
        if np.random.rand() < acceptance_ratio:
            beta = beta_new

        samples.append(beta)

    return np.array(samples)


# Prediction function based on posterior samples
def predict(X, posterior_samples):
    predictions = np.dot(X, posterior_samples.T)
    return predictions.mean(axis=1), predictions.std(axis=1)


def predictLeastSquare(X, beta):
    # Add intercept term to the new data
    X_with_intercept = np.hstack([np.ones((X.shape[0], 1)), X])
    return X_with_intercept @ beta

# Forecast function for new input data
def forecast(X_new, posterior_samples):
    # X_new is the new data for which we want to forecast
    return predict(X_new, posterior_samples)


def linear_least_squares(X, y):
    # Add intercept term (column of ones)
    X_with_intercept = np.hstack([np.ones((X.shape[0], 1)), X])

    # Compute the coefficients using the normal equation: β = (X^T X)^(-1) X^T y
    beta = np.linalg.inv(X_with_intercept.T @ X_with_intercept) @ X_with_intercept.T @ y
    return beta

# Main function
def run_bayesian_linear_regression(csv_path, n_samples=1000000, step_size=0.1):
    # Load data from CSV
    data = load_data(csv_path)

    # Assuming the last column is the target 'y' and the rest are features 'X'
    X = data.iloc[:, :-1].values
    y = data.iloc[:, -1].values

    beta = linear_least_squares(X, y)
    y_pred = predictLeastSquare(X, beta)

    # Add a column of ones for the intercept term (beta_0)
    X = np.hstack([np.ones((X.shape[0], 1)), X])

    # Initialize the parameters (beta_0, beta_1, ...)
    initial_beta = np.zeros(X.shape[1])

    # Perform Metropolis-Hastings sampling
    posterior_samples = metropolis_hastings(X, y, initial_beta, n_samples, step_size)

    # Get predictions
    mean_predictions, std_predictions = predict(X, posterior_samples)


    # Plot results
    plt.figure(figsize=(10, 6))
    plt.plot(y, label="True values")
    plt.plot(mean_predictions, label="Bayesian Regr")
    plt.xlabel("Time")
    plt.ylabel("1/V")
    plt.fill_between(range(len(y)), mean_predictions -  std_predictions,
                     mean_predictions + std_predictions, alpha=0.3, label="Err")
    plt.plot(y_pred, label="Least Square", color='red')

    plt.legend()
    plt.show()

    return posterior_samples


# Usage example
csv_path = 'dati.csv'  # Replace with the path to your CSV file
posterior_samples = run_bayesian_linear_regression(csv_path)


new_data = pd.DataFrame({
    'feature1': [30, 35]  # Replace with your features
})

X_new = np.hstack([np.ones((new_data.shape[0], 1)), new_data.values])

# Forecast values using the posterior samples
mean_forecast, std_forecast = forecast(X_new, posterior_samples)

# Output the forecasts
print("Forecasted values (mean):", mean_forecast)
print("Forecasted values (95% CI):", mean_forecast - 1.96 * std_forecast, mean_forecast + 1.96 * std_forecast)