Arx and some more

[mapi1]

Hey hey
I would like to contribute some more to this repository, and thought this discussion tool would be a good place to talk about some ideas and existing code that may fits here.

Again the arx method ... I found that for certain problems I need a 'direct input / static gain' which is currently not directly possible. I could think of two different ways for adding this, that are not breaking existing code:

1. adding another keyword like direct with the same structure as nb consisting off Bools indicating if direct input is desired or not and defaults to fill(false, size(nb)). (I tested this in my fork already)
2. making direct input the default but at the same time setting inputdelay to 1 as a fallback.
   The first option offers the most flexibility and makes it explicit but also more complex, while the second option is minimal invasive but maybe there is also a better way doing this?

Existing code that may fit here:

• A detrend(signal/iddata; order::Int = 1) function, that removes the mean and linear/ polynomial trend of order order. This function is very prominent in introductory textbooks on System Identification using matlab.
• Fast least trimmed squares (Rousseeuw et al. - Computing LTS Regression for Large Data Sets, 1999), as a robust version of \ for up to 50% outlier. Not sure if this is still state of the art or if this belongs here in general, as I saw a lot to solve Least squares problems is done in TotalLeastSquares.jl.
• In the 80s the method was know as generalized least squares but I think the name changed or vanished. Basically a method to fit the model Ay = Bu + (1/C)e by iteratively fitting an arx model, estimating an ar model on the residuals and then filtering the IO data (Söderström - Convergence Properties of the Generalised Least Squares Identitication Method, 1974). Not sure if this duplicate to some functionality already provided, but it was the method of choice for a physiological model I experimented with so there it is. If there is a more modern approach to this I would also be happy to hear about this :D

Ideas:

When I started out with System Identification I often struggled with small tasks that where not really covered in the literature but handled by all the convenience functions provided by matlab, like the before mentioned detrend. So with providing more of those in a well documented form, the hurdle for using julia instead of matlab when diving into Control Theory / System Identification could be reduced.

So far with what I could contribute, feedback about how useful this is is welcome!

[baggepinnen]

Hey, I would be very happy to have you as a continued contributor to this package! Let me address the topics you bring up

• It makes sense to support a direct term in arx, with the current definition of inputdelay I guess it would correspond to inputdelay=-1 (which is a bit unfortunate, but could perhaps be documented?). Having inputdelay default to 1 and make inputdelay=0 correspond to the model having a direct term would make more sense, but would technically be breaking. I guess most systems identified in practice does not have a direct term (unless modeled like a direct term due to very slow sampling), so it makes sense to default to identifying a system without direct term. In what way do you see the extra argument direct offer more flexibility than choosing inputdelay?
• A detrend function could perhaps be useful, but most introductory texts do not explain how to use it properly. It's one thing to detrend the inputs, but there must be proper "retrending" of the output in order to make predictions useful. If such a feature could be well supported, perhaps in a similar fashion to how MultivariateStats.jl handle whitening transforms etc., I would be for adding detrend functionality.
• Fast least trimmed squares: This would fit much better in TotalLeastSquares.jl, which already contains the robust irls, iteratively reweighted least squares.
• Generalized least-squares: I'm not familiar with this method. It sounds a bit like it is similar to irls but for dynamical systems, but I'll have to read the paper to get familiar with it. It also shares some similarities with pseudo-linear regression (plr in this toolbox) that makes two fits in order to fit a noise model as well. If you find this method useful, I see no reason not to have an implementation of it in this toolbox :) In general though, I tend to fit statespace models if I want to have a general model with poles, zeros and noise models, they tend to be more general.
• Other convenience functions: feel free to contribute as many of those as you feel are required. I'm up for discussing interfaces etc., many matlab interfaces are poorly designed and some are replaced by newer ones, making the user confused about what function to call since both are available.

[mapi1]

Thats great!

1. Yes, now it would correspondent to inputdelay = -1, which could be documented. Setting inputdelay = 0 correspondent to having a direct term, is not necessarily breaking but probably needs a dedicated version of params2ploy. This may be used to bring inputdelay and a direct term to other functions, like plr? Having a direct keyword allows me at the moment to have a direct term, then some inputdelay and afterwards some delayed terms. In this way I can bring some a priori physiological knowledge about fast and slow pathways into the estimation, but this is probably already a very special case and the ultimate abstraction would probably be to allow Boolean vectors for nb manually selecting a input structure. :D
2. True, I never really spent a thought about retrending so far. If I understand the approach of MultivariateStats.jl right, that would mean the introduction of a Trend type storing the detrended information for later use together with functions like predict or simulate. I will do some experimentation to see how a concept like this could play out. Maybe there is also a smart way to store it together with the iddata, with trend-types inheriting from the respective Data-type containing now detrended data and trend information and sth. like retrend would return the original data. But that is just some wild speculation on my side.
3. I prepare a PR to TotalLeastSquares.jl then. :)
4. Yeah the idea is similar to irls. When you say to fit statespace models, are you talking about n4sid? Thats the method I experimented with so far and achieved quite good results for some parts of the model, but other parts I could not recover what I put in. Is there a way to impose some structure onto the statespace model?
5. I think I will start with the above and then see what comes up :)

[baggepinnen]

1. I think inputdelay = 0 corresponding to a direct term would be most logical. This would apply to pretty much all functions that call get_ARXregressor.

the ultimate abstraction would probably be to allow Boolean vectors for nb manually selecting a input structure

I think the ultimate abstration would be to allow the user to specify a constraint function, but at that point it's hard to provide generic functionality and the user might as well use the lower-level functions and call any solver of choice.

2. A Trend structure would correspond to the MvS.jl approach. Not sure how to handle it though, should the user be responsible for applying it, the model wrapped in another model with the trend information or should it be passed along to simulate etc.? Depending on the detrenting, the model polynomials may be modified after estimation to account for the detrending, in which case no further action is needed, the model is already in the correct format, it was just estimated using detrending. I think this is an attractive property.
3. Yeah, like n4sid or PEM in general. The PEM implementation in this package is not great, but I have a prototype that makes use of LowLevelParticleFilters.jl to provide the state estimation that would potentially be much more general. With PEM, you can provide a regularization function and custom cost functions that can encode any structure you want, but the interface is not great.