Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Problem with estimate_delay? #38

Closed
finmod opened this issue Oct 12, 2017 · 13 comments
Closed

Problem with estimate_delay? #38

finmod opened this issue Oct 12, 2017 · 13 comments

Comments

@finmod
Copy link

finmod commented Oct 12, 2017

Just trying to run the nlts example for continuous time. The command:

using DynamicalSystems, PyPlot

ds = Systems.lorenz() # Max lyapunov is around 0.90
dt = 0.05
x = timeseries(ds, 1000.0; dt = dt)[:, 1]
τ1 = estimate_delay(x) #gives 7

Returns error:

UndefVarError: x not defined

Stacktrace:
[1] estimate_delay(::Array{Float64,1}) at C:\Users\Denis.julia\v0.6\DynamicalSystems\src\delay_coords.jl:186
[2] include_string(::String, ::String) at .\loading.jl:515

@Datseris
Copy link
Member

Datseris commented Oct 13, 2017

Yes that was a typo. I was writing s instead of x at some points.

Thank you very much for pointing this out. It is now fixed on master branch
(you can use Pkg.checkout() to get there)

A new release will be coming soon but though.

@Datseris
Copy link
Member

I've just tagged a new release. It will probably be merged in one day.

@finmod
Copy link
Author

finmod commented Oct 13, 2017

Success! Nice work. I have a follow up question for you: are you planning to extend DynamicalSysytems.jl to parameter estimation from noisy data? For instance recover all three parameters of Rossler from the three synthetic noisy timeseries.

I have been participating to JuliaDiffDocs Roadmap on this topic but I am frustrated by the lack of progress in Julia with respect to what has been achieved in Matlab or Python.

@Datseris
Copy link
Member

Keep in mind that Julia is much, much younger than Python or Matlab. And regardless, in my eyes at least, already way outmatches both of them. In addition, due to Julia being open-sourced and easy to contribute, if something bothers you you can simply improve it immediatelly; just contribute code!

Can you elaborate a bit more on what you mean by parameter estimation for synthetic data?

Specifically, at which parameters you are referring to? The parameters of the Reconstruction or parameters of a model a user provides?

For the first one you don't need any interface. Just start reconstruct with different parameters and see what fits best.

There is no interface for the latter and I will also not do it. I have other methods I want to develop that are pure to nonlinear dynamics and chaos, and also sophisticated enough to require an expert to write them. Parameter estimation is already possible from many different packages, including DiffEq, and also the nonlinear optimization packages. You can combine these with the reconstruct from my package and have pretty decent results!

Of course if you want to contribute a pretty interface, I would gladly accept.

@ChrisRackauckas
Copy link
Member

ChrisRackauckas commented Oct 13, 2017

Success! Nice work. I have a follow up question for you: are you planning to extend DynamicalSysytems.jl to parameter estimation from noisy data? For instance recover all three parameters of Rossler from the three synthetic noisy timeseries. I have been participating to JuliaDiffDocs Roadmap on this topic but I am frustrated by the lack of progress in Julia with respect to what has been achieved in Matlab or Python.

I think @finmod is mentioning the JuliaDiffEq work here and although it's only a side project of a side project since most of what I have been working on is the solvers themselves, it's pretty completely false that there has been a lack of progress here. Let's go through the current feature set and then see what we are missing.

We already have Normal log-likelihood (L2) estimation (non-log is just a bad idea that gets worse results because of numerical issues) with choices of regularization and your choice of optimization algorithm. Other than sensitivity analysis (which is just a faster way of computing gradients w.r.t. parameters than autodiff in this case), this means we have everything AMIGO2 has. And we are benchmarking very well too. Then we also have DiffEqBayes.jl for Bayesian estimation which is under construction but already working. What we don't have are arbitrary likelihood fits on SDEs, but I don't know of any other package which has that anyways.

On the equations that AMIGO2 is testing on we are doing just as well (usually better due to higher accuracy and better global optimizers). So let's not talk about generalities since that doesn't make sense. More specifically, what we have does really well in comparison to other packages for non-chaotic models.

What we don't cover is chaotic models. That is a very new subject and I don't know of anything which is truly robust, but the best you can get according to the literature is:

http://www.sciencedirect.com/science/article/pii/S0375960116000839
http://www.sciencedirect.com/science/article/pii/S1877750315000149

Essentially DiffEqParamEstim.jl can already do that if someone implements an optimizer that does PSO or ADA. PSO already exists because of DiffEqParamEstim + Evolutionary.jl, which I know you found because of SciML/DiffEqParamEstim.jl#48 and I told you that you need to Pkg.checkout("Evolutionary") to make it work because the author of that library isn't tagging a new version.

I still don't think that's going to be the holy grail though since you're using a tougher test problem than what I usually see in the papers: a long-time integration of a Lorenz system (this is the quintessential super hard numerical problem anyways). But there is a clear trend with those chaotic estimation papers: they use the same strategy we use in DiffEq but with newer and fancier global optimizers. In that case, to get the newest tools you'd likely want to talk to (/help out) the optimization community build more global optimizers. In that sense, @Datseris is doing a lot of great work but unless you know of a completely different method for doing this, it doesn't sound like the next developments are within the scope of DynamicalSystems.jl

Of course @finmod, feel free to chime in with new methods as well because now that we've exhausted what we see in AMIGO2 (sans multiple shooting) I'm not aware of other commonly used methods.

@Datseris
Copy link
Member

Good that Chris replied here because I got no clue for parameter estimation (apparently it is not the same thing as what I was thinking)!

@finmod
Copy link
Author

finmod commented Oct 15, 2017

Your nlts.jl brought me back to my starting point in 2014 with Ulrich Parlitz and Jan Schumann-Bischoff. Back then parameter estimation of their usual benchmark models (the Hindmarsh-Rose model, the Colpitts oscillator, and the Rössler system) was through a cumbersome mix of Python and C, ADOL-C and all that jazz. This has been my main reason to turn to Julia and discover DifferentialEquations.

The novelty that DynamicalSystems.jl is bringing to this old problem is to turn the problem on its head by infering the dynamical characteristics and parameters from delay reconstruction. As you indicate ,nlts.jl implements the Parlitz chapter of Chaos, Detection and Predictability and promises estimating reconstruction parameters "more coming soon!". This is what my original question was about. The code in delay-coords.jl is already set up for some parameter estimation. I just want to make sure that you are dealing with "Estimability and dependency analysis of model parameters based on delay coordinates, J. Schumann-Bischoff, S. Luther, and U. Parlitz, Phys. Rev. E 94, 032221 – Published 28 September 2016. This has more practical relevance as it hints to multivariate time series and models with several structural parameters (for Lorenz, 3 parameters recovered to compare with the 3 original parameters used to produce the timeseries).

As for the interaction with JuliaDiffEq, I just see two complementary methods: statistical estimation that was followed by Schumann-Bischoff (2011-2015) and the delay reconstruction method followed here. I just wish that the testing models be identical with some serious paramter estimation challenges.

@Datseris
Copy link
Member

@finmod before I say anything, keep in mind that "promise" is a very strong word. I don't promise anything; this is an open-sourced project and I am developing it purely out of hobby. I don't need any of the methods there for my work. So "coming soon" simply means "coming whenever I have too much free time and I feel like doing that, which could be in many years from now". Anybody can contribute by the way, and you seem to be very knowledgable in the field.

Another thing to keep in mind is that there is no novelty per se in this package. All methods implemented are well known and published.

Now, about the "coming soon": The estimate_delay is already implemented with 2 methods; only the first is exported though (the second is to find the 0 of the correlation, which I simply find a bad approach). The other half of the story is a method that estimates the best reconstruction dimension. I know of a method again from discussion with Parlitz that simply projects and finds the false nearest neighbors. But it is not trivial to implement so it won't be coming "soon". Maybe the paper you are citing has a nice way to do it. (I removed the "coming soon" expression from the docs because I realise now that it is a misdirection)

I am making this quick reply because I do not have the time to check the papers you are citing today. When I check them out I will write another message.

P.S.: Have you worked with Parlitz? This is the 2nd year I am tutoring his course so I know him personally well!

@ChrisRackauckas
Copy link
Member

I see, these methods are the methods specific to chaotic systems which is completely different than the optimization approaches usually taken. I read through these delay methods and yeah, for things like the Lorenz attractor this is probably the thing to do. I should probably add a mention in the DiffEq docs.

I don't think that the models that are tested really need to be the same since, as I noted before, I would not expect the "standard" approaches to do well on chaotic systems because the sensitivity due to chaos is exactly what would make the parameter estimation sensitive. It should just be well documented that when one thinks the solution lives on a chaotic attractor that they should look to these other methods.

As always, good stuff @Datseris

@finmod
Copy link
Author

finmod commented Oct 16, 2017

I invite you to spare a few minutes for this recent video lecture of Ulrich Parlitz (https://www.youtube.com/watch?v=CiZQocRu6-U). It shows that the estimation of chaotic dyanmical systems is an eminently practical and a life-saving activity! If you have only 5 minutes, you can start at 22.00 minute with the super fast speed reading of the 2016 article at the 25.00 point.
And yes, you will note that Ulrich Parlitz is a mature man like me and has been around for a while. His message is interesting as well as refreshing.

@ChrisRackauckas
Copy link
Member

That talk is great and explains it well. I think it fits in well with the DynamicalSystems.jl sphere of influence, and I'll make a note in my docs saying that data on chaotic attractors should probably be estimated using this method.

It shows that the estimation of chaotic dyanmical systems is an eminently practical and a life-saving activity!

It is practical but I think this is a good point to have separation of concerns. I'm working on solvers and addon functionality which includes parameter estimation in the general sense, but am letting different domains develop their own domain-specific tools. @Datseris is the dynamical systems domain-expert who can handle this domain, and of course I am lending a hand (but you force-pushed over my commits so it shows I have no contributions haha).

As for swapping over, you can get the ODEProblem from any ContinuousSystem by 1 function call, so that's an easy way to handle it. A reverse function (which assumes t=0 always works) would be nice as well.

@Datseris
Copy link
Member

Yeah the assumption of t=0 comes from the fact that I realized its easier to handle all systems as non-autonomous with one extra dimension.

I don't know if having explicit time can offer optimizations for the solvers of DiffEq, but I kind of hoped not when I made the choice.

Sorry about the force pull I may be expert in dynamical systems but I am a complete noob in git :D

One last comment: @finmod maybe you would like to join or gitter chat channel, since this discussion is way off topic in this current issue page! The gitter channel is here : https://gitter.im/JuliaDynamics/Lobby

@ChrisRackauckas
Copy link
Member

I don't know if having explicit time can offer optimizations for the solvers of DiffEq, but I kind of hoped not when I made the choice.

It doesn't.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants