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

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