Model types and algorithms

[ElOceanografo]

As we add more model types and algorithms, I've started thinking about how best to conceptualize and organize them. Right now, we've been adding a new type and associated methods for each filtering algorithm. The existing ones are summarized in this table:

Type	Supertype	State evolution	Observation	Noise	Input	Filtering algorithm
LinearGaussianSSM	AbstractLinearGaussian	Matrix	Matrix	Gaussian	No	Kalman
LinearGaussianCISSM	AbstractLinearGaussian	Matrix	Matrix	Gaussian	Yes	Kalman
NonlinearGaussianSSM	AbstractGaussianSSM	Function	Function	Gaussian	No	EKF
AdditiveNonLinUKFSSM	AbstractSSM	Function	Function	Gaussian (?)	No	UKF

In addition, the following are on the wish list:

Type	Supertype	State evolution	Observation	Noise	Input	Filtering algorithm
NonlinearGaussianSSM	AbstractSSM	Function	Function	Approx. Gauss.	Yes	Ensemble Kalman
NonlinearSSM	AbstractSSM	Function	Function	Anything	Yes	Particle Filter

And finally, we want to add time-varying evolution and observation matrices.

There's some overlap between these types in terms of the information they encapsulate. There is also some overlap in which filtering algorithms could be used for each type of model. E.g., a nonlinear Gaussian model can be filtered with the EKF, UKF, EnKF, or PF.

--	Linear	Nonlinear
Gaussian	Any	All but KF
Non-Gaussian	PF (maybe EnKF)	PF (maybe EnKF)

I can see users might want to try more than one filter on the same model--for instance, they've got a nonlinear Gaussian model, try the EKF first, decide it doesn't work well enough, and switch to the UKF or EnKF.

I see three basic options:

1. One type per algorithm, with type-specific update, filter, and smooth methods...basically what we have now. Maybe write convert functions to switch easily between model types.
2. Reduce the number of redundant/aliased types, and have algorithm-specific functions, e.g. just one NonlinearGaussianSSM type with filter_kf, filter_ekf, etc.
3. Combine redundant/aliased types, keep one update/filter/smooth function that takes an optional argument for the algorithm to use, e.g. filter(my_lingaussm, y, x0, method="EnKF").

Thoughts?

[JonnyCBB]

Definitely a good point to fix this up. Thanks for bringing it up.

If I think about what the problem set up is: we have a system where we think we know the process with some level of confidence and we have have observations with some observation error. This, at least to me, is the basic model for which a filtering algorithm can be applied. Any suitable algorithm! Therefore it would seem sensible to me that you could create the model without specifying its type, unlike what we do now. With that said, it narrows it down to option 2. or 3.

Here is where I can't really decide on the same basis as above. My personal preference is option 3. I like the idea that you have a model and you want to filter it (or smooth it or update it or whatever) and you can just write something like filter(model, y, x0, method="EKF"). A bit like in _JuliaOpt/Optim.lj_ e.g. optimize(f, [0.0, 0.0], method = :nelder_mead). I think this is more consistent for a user API personally.

However if there are any objections I would love to hear them.

[ElOceanografo]

Option 3 (b): define a type for each filter, and have methods that dispatch on both the model and filter type. So if your model is, say, a NonlinearGaussianSSM, you could do either e.g. filter(model, data, EKF()) or filter(model, data, UKF()). For a particle filter, it could be filter(model, data, PF(nparticles)).

I've been thinking about it and I think I like this best...

[ElOceanografo]

Yes, we would have a linear Gaussian type, a nonlinear Gaussian type, and an "anything goes" type, name TBD, that can only be filtered with the ensemble Kalman or particle filters.

I've just started work on porting @codles implementation in TimeModels, which hopefully will let these types be time-varying or invariant, with or without control input and feedback.

On Sun, Nov 22, 2015 at 3:32 PM, Jonny Brooks-Bartlett notifications@github.com wrote:

I like it. This solution should also avoid the horrible UKFParameters type that I introduced for the UKF. So with the new type proposition it could read something like filter(model, data, UKF(a, b, k)).

To be clear:

There are types for the model — Reply to this email directly or view it on GitHub.

[JonnyCBB]

I like it. This solution should also avoid the horrible UKFParameters type that I introduced for the UKF. So with the new type proposition it could read something like filter(model, data, UKF(a, b, k)) should the user want to specify the UKF parameters (we can use the default parameter values if they are not explicitly passed to the UKF type).

To be clear:

• There are types for the model, which the user should have distinguished already.
  • LinearGaussianSSM
  • NonlinearGaussianSSM
  • NonlinearSSM
• And there are types for the filtering algorithms, which I guess would extend to to the smoothing and maybe parameter estimation algorithms. For the filtering we would have
  • KF or BKF (for the standard linear Kalman filter)
  • EKF
  • UKF
  • etc.

Since this proposition would mean using types for the algorithms as well as the models we should decide on the available types for these or at the very least, decide how to determine whether a new type is needed.

Types for filtering algorithms

The possible types for the filter algorithms are straightforward I guess - implementation of a new filter algorithm requires a new type. We just have to be aware of when a particular algorithm is a special case of a more general algorithm.

Types for models

The possible types for the models may be less straightforward if my understanding isn't too naive. Since there are types for the different algorithms I think it suffices to have model types that differ in the types of input that are required. For example we need different types for LinearGaussianSSM and NonlinearGaussianSSM because the former requires process and observation matrices whereas the latter requires process and observation functions. My understanding (although I'll need to do some reading to confirm this) is that particle filters require the explicit probability distributions for the process (prior) and the observation (likelihood) so would also need a different type NonlinearSSM. I really have no idea about the EnKF so again I have some reading to do.

I hope this makes sense.

[rob-luke]

I like options 2 and 3. Should I hold off on implementing time-varying models until these changes are made? I haven't gotten around to starting it yet.

[JonnyCBB]

Yes I would hold off implementing the time-varying stuff for now. Sounds like @ElOceanografo has started implementing some time varying stuff. But I imagine the changes proposed here are likely to change the file structure as well as the algorithm implementations.

[JonnyCBB]

@ElOceanografo it sounds like we're happy to go with the proposed changes. We should probably get going with these changes sooner rather than later so that the code base foundation is stable for other additions such as additional filtering algorithms and time varying matrices.

Are you going to start with the changes? If not, I think I may have some time over the Xmas break to have a good go at this and I'm happy to do it. I just want to make sure that we're not duplicating efforts.

[ElOceanografo]

Cool. I've started doing some of this refactoring on my local branch, when I get it working (hopefully in the next week or so) I'll push it to the /dev branch.

[ElOceanografo]

Okay, I've pushed my tentative changes to the dev branch. All tests that were passing before are still passing. I've reorganized a fair amount, but most of the actual business code is unchanged. Comment period is open!

[ElOceanografo]

Merged the new implementation into master 684012f9aca5894fd98ba0aa2a061c5542fa3a81