oseiskar / simdkalman

Python Kalman filters vectorized as Single Instruction, Multiple Data
https://simdkalman.readthedocs.io/
MIT License
176 stars 36 forks source link

Regression with time-varying coefficients #14

Closed cs224 closed 4 years ago

cs224 commented 4 years ago

Hello,

I am trying to implement regression with time-varying coefficients as described in Time Series Analysis by State Space Methods: Second Edition section 3.6.1. For that purpose I use the last example given at lectures: kalman-filters.

I somehow do not manage to replicate the results of pykalman and cannot tell why. I cannot see that I am doing something wrong.

I've summarized my results in the following notebook: https://nbviewer.jupyter.org/gist/cs224/031305297abc17a5546d4f0c2073716c

First I replicate the desired results with pykalman and then I try to do the same with simdkalman. simdkalman only returns my first state parameter and a copy of the first state parameter instead of returning the filter results for both state parameters.

I'd prefer to use simdkalman as it seems to be faster and actively developed.

Thanks a lot, Christian

oseiskar commented 4 years ago

Hi. Just checked the documentation and noticed that this part is not really covered, but:

If you have any time-varying matrices in the model, then you must use the primitives module

This is the most comprehensive example, which demonstrates how to use time-varying noise. In your case, the time-varying parameter is the observation model (matrix H, which is called Z in your notebook) so you'll have to modify the example accordingly.

If you do not care the about the Kalman-smoothed estimates, you can skip the "backward pass" and just take the filtered values that are appended to the m_list in the example

cs224 commented 4 years ago

Ok, I tried to follow your advice here: https://nbviewer.jupyter.org/gist/cs224/49d3ba15219b82d6244b4fcb296f240d Look at the very bottom. Is that how you would have done it yourself or would you do it differently?

While the slope component matches very exactly the pykalman result, the intercept component is different. The difference is small (0.00189), but noticeable. What do you think where this difference comes from? Which result is more correct, the pykalman or the simdkalman result? Or did I make a mistake somwhere along the way that would explain the difference?

oseiskar commented 4 years ago

I turns out my earlier examples were a bit sloppy. Great that you paid attention to this, thanks.

With a single-line fix to your notebook

-    m, P = simdkalman.primitives.predict(m, P, state_transition, process_noise)
+    if i > 0: m, P = simdkalman.primitives.predict(m, P, state_transition, process_noise)

the results match perfectly. (Also note that both the slope and intercept components were different, the slope diff was just below the np.allclose default thresholds.)

The difference was how the initial value is interpreted and at which point you do the first "predict" operation. The fixed version (which also used by both simdkalman's and pykalman's filtering routines in case of fixed coefficients) makes more sense in this case, as there is not much point to apply the predict step to the initial value before the first update.

As a side note, this is not universally true. For example, in various sensor fusion use cases, the predict operations are not all matched with an update operation, i.e., there are a lot more predict steps than update steps. Therefore I would not say that the existing example was strictly "incorrect". It's mostly a matter of convention and what makes sense for the problem at hand.

cs224 commented 4 years ago

Thank you very much for your correction and detailed explanation! I'll close this issue now as I think you've fully resolved all my open questions. Thanks again! Christian