Skip to content
Snippets Groups Projects
Commit 8ffee955 authored by FROGE Ewen's avatar FROGE Ewen
Browse files

Upload New File

parent 85d8348c
No related branches found
No related tags found
No related merge requests found
#%%
# This script generates a Lorenz system time series, splits it into training and testing datasets,
# and applies an analog forecasting method using nearest neighbors. The training data (bib) is used
# to build a database of past states, while the testing data (signal) is predicted based on analogs
# found in the training set. The script then visualizes the innovation (prediction error) and its PDF.
# Modules
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from Compute_Innovation import create_database, prediction_neighbors
# Parameters
p = 3 # Size of analogs
k = 70 # Number of analogs
N = 2**16 # Signal size
overlapping = False # Whether to allow overlapping analogs in database
weighting = True # Whether to weight analogs based on distance
normalizing = True # Whether to remove analogs mean
def lorenz_system(t, state, sigma=10.0, beta=8/3, rho=28.0):
"""
Defines the Lorenz system equations.
"""
x, y, z = state
dx = sigma * (y - x)
dy = x * (rho - z) - y
dz = x * y - beta * z
return [dx, dy, dz]
def generate_lorenz_signal(n_steps, dt=0.01):
"""
Generates a Lorenz system time series (x component).
"""
t_span = (0, n_steps * dt)
initial_state = [1.0, 1.0, 1.0]
t_eval = np.linspace(t_span[0], t_span[1], n_steps)
sol = solve_ivp(lorenz_system, t_span, initial_state, t_eval=t_eval, method='RK45')
return sol.y[0] # Extract x-component
data = generate_lorenz_signal(2 * N)
bib = (data[:N] - np.mean(data[:N])) / np.std(data[:N])
signal = (data[N:] - np.mean(data[N:])) / np.std(data[N:])
db = create_database(bib, p, overlapping=overlapping)
vals, variances = prediction_neighbors(
signal, db, k, weighting=weighting, normalize=normalizing
)
innovation = signal - vals # First p values are NaN due to alignment
# Plot results
fig, axs = plt.subplots(2, 2, figsize=(10, 8))
axs[0, 0].plot(innovation)
axs[0, 0].set_title("Innovation")
axs[0, 1].plot(signal)
axs[0, 1].set_title("Signal")
axs[1, 0].hist(innovation[~np.isnan(innovation)], bins=50, density=True)
axs[1, 0].set_title("PDF of Innovation")
axs[1, 1].hist(signal, bins=50, density=True)
axs[1, 1].set_title("PDF of Signal")
plt.tight_layout()
plt.show()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment