Generalized least squares

[mapi1]

This is a first draft of the method I mentioned. I am not happy with it, as i think it could be done easier. I would be really happy for some feedback, especially the comments marked with # ->. Failing test are related to the issue with ar() #62 . Hope this makes some sense, and thanks in advance!

[baggepinnen]

I've had a first look through the code and have some minor comments

• the word residue is commonly used in the complex-analysis sense when working with transfer functions, perhaps the word residuals is more suitable to avoid confusion.
• camelCase function naming deviates slightly from the julia style guide which I mostly try to follow
• the functions num/den are deprecated, denpoly/denvec are encouraged instead.
• Having input arguments determining the number of output arguments is non-standard in julia since it causes type-instability (unless the compiler is smart enough to eliminate the problem with constant propagation). I think it's okay to simply always return everything, never return the optional things or return everything in the form of a named tuple or custom struct.
• I can't really figure out what's going on in the residue function, in particular, the denominator of a transfer function is considered an impluse-reponse? The lines

  denominator = den(Gest[1, 1])[1]
  A = impRes2tf(denominator)

  appears to be a convoluted way of inverting the transfer function?

I am not really the greatest fan of the name Generalized Least squares, since you can imagine a million ways of generalizing any method and the name does not contribute any guidance on in what way it is generalized. I am myself guilty of naming another function plr for pseudo-linear regression which is an equally bad name. Ideally, the functions hould be named after which kind of model it estimates, like ar and arma, but I realize that's not always feasible since there is no particular name (to my knowledge) for this particular model structure. Perhaps some variant of armax? Hopefully you have a better suggestion that would make it obvious what the function does :)

[mapi1]

How does the saying go, 'there are only two hard things in CS, cache invalidation and naming things'

I did a revision of the method on the basis of your remarks and the new ar definition. Two points that remain relevant are:

• The name gls: I agree that the name is sub optimal, as this kind of model can be rewritten as ADy = BDu + e it can be understood as an arx model with the constraint, that A and B share a common factor, so maybe constrainedARX or similar would give a better name.
• The residuals function: I replaced residue , there are now a residuals and a predict function that are inspired by the predict for ar. These are more in the style of the package I would say. Maybe it makes sense to even document and export them? I therefore also updated params to handle MISO tfs, but I think I will need to add some more tests to be sure it still plays well with the rest.

[mapi1]

Not sure if it is best practice to fix it here, but I found that adding inputdelay, actually broke params for some cases when stochastic = trueI added some tests for these cases now. I created a compareTFs function to compare normal tfs to those using Particles, maybe there is a better way to do this?

[baggepinnen]

Comparing probability distributions can be really tricky, mostly because there are many different ways you can imagine doing it. in this case, the posterior distributions over the parameters are all Gaussians, so comparing their means like you did is an okay approach. The stochastic part should really be tested separately which I see has not been done properly before (no worries, that's a story for another time :)

[mapi1]

I found a name for that model structure which makes sense to me, ARXAR. It is not widely established but I think as a function name arxar is clearer about what is going on, namely fitting an ARX model with an AR noise model. So if you think these changes are fine, I would consider this PR good to go. I personally would like to add the 'stochastic' keyword here as well, but I think I have to dive deeper into MonteCarloMeasurements first and figure out the math.

[baggepinnen]

I'm not sure why CI is not running on this PR. I wonder if it has anything to do with it starting out as a draft? Closing and opening usually triggers CI but that trick did not work here.

[baggepinnen]

I merged it to get CI on the master branch, can always fix any eventual error there. Thanks for your hard work and excellent contribution with this PR!