Discrepancy between R^2 score and simulation results

[Thirie]

Hi everyone!

I'm trying to use SINDy on EEG data and currently, I'm trying to make use of the discrete_time=True option, in combination with coef ensembling. I'am using a custom library with sines and cosines with frequencies of 1-13Hz. I'm struggling to obtain correct results. When i use the model.score function, i get a R^2 score of 0.99, which is unrealistically high in my opinion. When i try to plot the model.simulate results, I unfortunately see that it does indeed not seem to capture the underlying dynamics, but only fits a low frequency component (see pictures). True is the pre-processed data, mean is the mean of all the coefs found in ensambling and stable mean is the mean of the coefs which yield a stable system.

This is my power spectrum:

This is my simulation combined with the true data (the x-axis is 1 second):

I therefore have some questions:

• Why is it that i only manage to fit a low frequency line to the high frequency data? I already used a low pass filter of 13Hz on the EEG data and have a sampling frequency of 500Hz
• Is the method i use correct for determining score and simulation results? If so, why do i get a high R^2 with poor simulations?
• I found it necessary to fit a second model without ensambling to swap out the coefficients and simulate different conditions. Is this correct, or is there an oversight in my code?

I attached the code i used below:

ode_lib2 = ps.CustomLibrary(
    library_functions=library_functions,
    function_names=library_function_names,
    interaction_only=False,
    #library_ensemble=True,
    include_bias=False  # Include bias term
)


optimizer = ps.SR3(
    threshold=0.0005, thresholder="l1", max_iter=1000, normalize_columns=True, tol=1e-1
)
original_model = ps.SINDy(feature_library=ode_lib2, optimizer=optimizer, differentiation_method=ps.SmoothedFiniteDifference(smoother_kws=smoother_kws), discrete_time=True)
original_model.fit(normalized_low_freq, t=time_list, multiple_trajectories=True, ensemble=True, n_models=20)
ensemble_coefs = copy.deepcopy(original_model.coef_list)
print("Normal form model: ")
original_model.print()
print("R2 score", original_model.score(normalized_low_freq, t=len(time_vector), multiple_trajectories=True))
print("MSE score",original_model.score(normalized_low_freq, t=len(time_vector), multiple_trajectories=True, metric=mean_squared_error))
#print("Normal fit coef ", optimizer.coef_)

#Fitting model second time, otherwise I cannot swap the coefficients
original_model = ps.SINDy(feature_library=ode_lib2, optimizer=optimizer, differentiation_method=ps.SmoothedFiniteDifference(smoother_kws=smoother_kws), discrete_time=True)
original_model.fit(normalized_low_freq, t=time_list, multiple_trajectories=True)

ensemble_ar = copy.deepcopy(np.asarray(ensemble_coefs))
#%% Model stable coefs
begin_index=1000 
end_index=1500
number_of_steps=end_index-begin_index

t_test=time_list[0][begin_index:end_index]
x_test_sim=normalized_low_freq[0][begin_index:end_index,:]

 stable_ensemble_coefs, x_test_sims = integration_metric_discrete(ensemble_ar, optimizer, t_test, normalized_low_freq, original_model, number_of_steps)

optimizer.coef_ = mean_coef
#print("Mean coef ", optimizer.coef_)
sim_mean = original_model.simulate(normalized_low_freq[0][begin_index, :], number_of_steps, integrator="odeint")
x_test_sim_mean = replace_negatives_with_nan(sim_mean)
    
##%% Simulate median coef
optimizer.coef_ = median_coef
#print("Median coef ", optimizer.coef_)
sim_median = original_model.simulate(normalized_low_freq[0][begin_index, :], number_of_steps, integrator="odeint")
x_test_sim_median = replace_negatives_with_nan(sim_median)
    
    
##%% Simulate stable models mean coef
optimizer.coef_ = stable_mean_coef
#print("Median coef ", optimizer.coef_)
sim_stable_mean = original_model.simulate(normalized_low_freq[0][begin_index, :], number_of_steps, integrator="odeint")
x_test_sim_stable_mean = replace_negatives_with_nan(sim_stable_mean)
    
##%% Simulate stable models median coef
optimizer.coef_ = stable_median_coef
#print("Median coef ", optimizer.coef_)
sim_stable_median = original_model.simulate(normalized_low_freq[0][begin_index, :], number_of_steps, integrator="odeint")
x_test_sim_stable_median = replace_negatives_with_nan(sim_stable_median)

# Compute the 2.5 and 97.5 percentile trajectories
bottom_line_raw = np.percentile(x_test_sims, 2.5, axis=0)
bottom_line = replace_negatives_with_nan(bottom_line_raw)
top_line_raw = np.percentile(x_test_sims, 97.5, axis=0)
top_line = replace_negatives_with_nan(top_line_raw)

t_test=time_list[0][begin_index:end_index]
x_test_sim=normalized_low_freq[0][begin_index:end_index,:]

plt.figure(figsize=(14, 15))
for i in range(3):
    plt.subplot(3, 1, i + 1)
    plt.plot(t_test, x_test_sim[:, i], 'k', label='True')
    plt.plot(t_test, x_test_sim_mean[:, i], 'b', label='Mean')
    plt.plot(t_test, x_test_sim_median[:, i], 'r', alpha=0.5, label='Median')
    plt.plot(t_test, x_test_sim_stable_median[:, i], 'y', alpha=0.5, label='Stable Median')
    plt.plot(t_test, x_test_sim_stable_mean[:, i], 'm', label='Stable Mean')
    plt.plot(t_test, bottom_line[:, i], 'g', alpha=0.6, label='95th percentile')
    plt.plot(t_test, top_line[:, i], 'g', alpha=0.6)
    plt.ylim(0.3, 0.8)
    ax = plt.gca()
    ax.fill_between(t_test, bottom_line[:, i], top_line[:, i], color='g', alpha=0.3)
    plt.grid(True)
    if i != 2:
        ax.set_xticklabels([])
    else:
        plt.xlabel('t', fontsize=20)
    plt.ylabel(feature_names[i], fontsize=20)
    plt.xticks(fontsize=16)
    plt.yticks(fontsize=16)
    ax.yaxis.set_label_coords(-0.05, 0.75 - 0.1 * (i + 1))
    if i == 0:
        ax.legend(bbox_to_anchor=(1.01, 1.05), fontsize=18)
plt.tight_layout()
plt.show()

plt.figure(figsize=(14, 5))
for x in x_test_sims:
    plt.plot(x)
    plt.ylim(0, 1)

[Jacob-Stevens-Haas]

Scoring is based on "if I have the correct point, how closely can I predict the next point?" Simulation is much more complex, as errors from previous steps are accrued nonlinearly. It's possible to have a model with a high R2 that diverges in simulation - sometimes even to infinity. Work on discovering stable models is ongoing.

[Thirie]

Thank you for your answer Jacob!
I struggle to understand how this is related to the optimal model discovery though. If the score is based on predicting the next point in time, given the current point, why wouldn't SINDy return x(n+1) = 1*x(n) given the sample frequency is high? The next point would almost always be very close to the last point when the sample frequency would be 500 Hz for example. And wouldn't this be the same when the signal is continuous?
I hope to hear from you!