dynamicslab / pysindy

A package for the sparse identification of nonlinear dynamical systems from data
https://pysindy.readthedocs.io/en/latest/
Other
1.41k stars 308 forks source link

Ensembling with weighted mean_squared_error w.r.t the test or training set #339

Closed SM-CAU closed 1 year ago

SM-CAU commented 1 year ago

Hello, I am trying out 'ensembling' of my time-series data of a Magnetoelectric-Sensor.

  1. By default model.print generates the median of the model co-efficients. But I want it to generate from the weighted average based on the mean-squared-error w.r.t the test_set. How can I do that? The function I have used for calculating weighted avg is attached.

weighted_mse

akaptano commented 1 year ago

What's wrong with this function? It looks to me that the line avg = ... is already computing the final coefficients as the weighted average of the ensemble models, using the mean-squared-error of the test set. So is the issue just setting these coefficients in the model? If so, try model.coef = avg or model.optimizer.coef = model.coef_ (I'm forgetting which one works right now)

SM-CAU commented 1 year ago

Thanks a lot. yes, the problem is setting them in the model so that model.print can use it to show the final equation. As this is a post processing, I meanafter I have printed the model with ensembling, how can I get the final equation using this weighted avg. This seems promising to me as this allows the model to see more of the dataset and generate based on that. So, I use your suggetion after the post processing and print it again? Do I insert this into ensemble_aggregator or just update the modelcoef and print the model again?

Problem 2: Tuning model Hyperparameters.: I used regularization (alpha) in the optimizer and trying to find the optimum value for 'threshold and regularization parameter 'alpha'. But for the threshold I come across a confusing trend in the generated plots and can't decide on the thresholds Is there any bug in the implementation? simulated xdot xdot

threshold_scan = np.linspace(0, 0.2, 20) coefs = [] rmse = mean_squared_error(data_x_train, np.zeros(data_x_train.shape), squared=False)

x_train_added_noise = x_train + np.random.normal(0, rmse / 10.0, x_train.shape)

for i, threshold in enumerate(threshold_scan): sparse_regression_optimizer = ps.STLSQ(threshold=threshold) model_sine = ps.SINDy( optimizer=sparse_regression_optimizer,

optimizer=sr3,

differentiation_method= sfd,
feature_names=["x", "xdot", "t"],
discrete_time=False,
feature_library=generalized_library
)
model_sine.fit(data_x_train, 1/sampling_rate, quiet=True)
#model_sine.fit(data_x_train, 1/sampling_rate, ensemble=True, replace=False, n_subset=1000, n_models=3, quiet=True, unbias=False, ensemble_aggregator=np.median)
coefs.append(model_sine.coefficients())

plot_pareto(coefs, sparse_regression_optimizer, model_sine, threshold_scan, data_x_train, t_train)

Jacob-Stevens-Haas commented 1 year ago

Hey Sayeed, here's how you could go about doing it: EnsembleOptimizer has a paramter ensemble_aggregator, which is a function to replace median. You can see in the code that it only can take a list of coefficients. It's easy enough to do that with a partial function. However, the challenge is that the optimizer doesn't have access to model.predict

What you could do is sublcass EnsembleOptimizer and add an attribute parent_model_features. You could then modify your subclass's _reduce() method to create a model bag_model with parent_model_features and set the coefficients manually for each fitted bag. Then run bag_model.predict() and compare to test data to determine weights.

Also FWIW, if you're using test data to choose weights in combining models, you might instead call it validation data instead of test data.

As to your second question - I don't believe there's a bug. The example you gave, with 10e9 error in derivative, is wild and not a part of the pareto frontier that I would expect to display the traditional elbow. Your example is also incomplete. If you can submit a MWE that demonstrates this more accurately, do it in a separate issue.

SM-CAU commented 1 year ago

Hello Stevens, Sorry for my late reply as last few weeks was too exhausting. Thanks a lot for your answer. Ensembling was latter point omitted as my model could predict without this, and we settled down for that. About the 10e9 error, its coming from one issue: In Sindy model.fit, the training data I fed is stacked measurement data from a sensor, their derivatives (up to this, this is required by the system dynamics of a duffing system which is a 2nd order differential equation) and ,ironically and I think inevitably as I see no other way around, time to make the driving force term of a duffing equation (F cos wt) be expressed with time. So, I had to pass time in the training data as a feature. The error is coming from them, I guess, because when I simulate both the measurement data and their derivatives (for train, validation, test all of them) forward in time with the model it fits nice with a R2 score of 0.99 and 1.

Jacob-Stevens-Haas commented 1 year ago

Ok, no worries. Do you have any remaining questions? The general answer is that ensembling based on error in the test/training set is a substantial feature that doesn't exist.

SM-CAU commented 1 year ago

Thank you so much for the reply. I have another question for which I have opened another issue 'https://github.com/dynamicslab/pysindy/issues/355'. I would be grateful if you can help me find the answer for that issue. However, I would be happy to have an expert opinion on one thing: Say I have a timeseries data what we assume to follow the Duffing system. Then I trained one portionof mydataset, validated with another portion, and tested for rest of the part. If I now make a new measurement by changing the excitation frequency (F cos (wt), this is the R.H.S of a duffing equation) w, should my developed model be able to predict the dataset without re-training on thenew dataset? As I understand it, it should not as the driving force is changes , the change should be reflected in the part the driving force is driving to. If so, what would be the gain for the built model. I a lil bit confused to find a good answer for that. However, the built algorithms are uploaded in the issue i have referred here. It would be helpful if you have time to have a look and make a comment onthe uilt algorithm. Thank you very much, again.

Jacob-Stevens-Haas commented 1 year ago

It depends how you train it. If w and t are both control inputs in the SINDy model, and you input the new w as a control variable in the new dataset, then yes, it will work. If the model has to discover w, then if you supply a new dataset with a different driving frequency, the model will struggle to simulate that system.

If that's what you're looking for, I think there's a paper by @Ohjeah on detecting parameter changes in SINDy models.

SM-CAU commented 1 year ago

Thank you so much for your reply. it was very helpful.

On Thu, Jul 13, 2023, 9:07 PM Jacob Stevens-Haas @.***> wrote:

It depends how you train it. If w and t are both control inputs in the SINDy model, and you input the new w as a control variable in the new dataset, then yes, it will work. If the model has to discover w, then if you supply a new dataset with a different driving frequency, the model will struggle to simulate that system.

If that's what you're looking for, I think there's a paper by @Ohjeah https://github.com/Ohjeah on detecting parameter changes in SINDy models.

— Reply to this email directly, view it on GitHub https://github.com/dynamicslab/pysindy/issues/339#issuecomment-1634754717, or unsubscribe https://github.com/notifications/unsubscribe-auth/A6BP2P34BJAEJQCFU3M7CMLXQBBQDANCNFSM6AAAAAAYGLCUV4 . You are receiving this because you modified the open/close state.Message ID: @.***>