Online Smoothing via update() and smooth_current()?

[mikediamore]

Hello,

This library has been quite helpful! I was wondering how to properly call the update and smooth_current functions to update the smoother in an online fashion.

For example, if I have the following to smooth 10000 normal random variables

data = np.random.normal(size=10000)

 ##KF specificiation taken from examples
kf = simdkalman.KalmanFilter(
    state_transition = [[1,1],[0,1]],        # matrix A
    process_noise = np.diag([0.1, 0.01]),    # Q
    observation_model = np.array([[1,0]]),   # H
    observation_noise = 1.0)

kf_results = kf.smooth(data)

And I wanted to get the smoothed value for the 10,001st term would I do the following?

#Snippet2

m = kf_results.states.mean[-1]
P = kf_results.states.mean[-1]
y = np.random.normal() #10,001st term
post_mean, post_cov, loglike = kf.update(m = m, P = P, y=y)
smooth_mean, smooth_cov, smooth_gain = kf.smooth_current(m = m, P = P, ms=post_mean Ps=post_cov)
new_smoothed_observation = smooth_mean
m = post_mean
P = posterior_cov

I'm finding that the above doesn't do a very good job of recreating the smoothing on the whole series especially if I continue this process over longer time periods. For example:

Here are the results from the last 500 points of the kf.smooth(data) call

And here are the results from using the following for loop, adapting #Snippet 2 above

m = kf_results.states.mean[-1]
P = kf_results.states.cov[-1]
smoothed_current_obs = []
for index in range(9500,10000):
    y = data[index] #get next obs to smooth
    post_mean, post_cov, loglike = kf.update(m = m, P = P, y=y.reshape(-1,1))
    smooth_mean, smooth_cov, smooth_gain = kf.smooth_current(m = m, P = P, ms=post_mean, Ps=post_cov)
    smoothed_current_obs.append(smooth_mean[0][0])
    m = post_mean
    P = post_cov    

Is it possible to recreate the output of kf.smooth() using kf.update() and kf.smooth_current()?

Here's the full code example to make it easier to run

import numpy as np
import simdkalman
import pandas as pd

data = np.random.normal(size=10000)

 ##KF specificiation taken from examples
kf = simdkalman.KalmanFilter(
    state_transition = [[1,1],[0,1]],        # matrix A
    process_noise = np.diag([0.1, 0.01]),    # Q
    observation_model = np.array([[1,0]]),   # H
    observation_noise = 1.0)

kf_results = kf.smooth(data)
smoothed_obs = kf_results.observations.mean

#Make Plots
df = pd.DataFrame([data,smoothed_obs]).T
df.columns = ['data','smoothed_obs']
df.iloc[-500:].plot()

m = kf_results.states.mean[-1]
P = kf_results.states.cov[-1]
smoothed_current_obs = []
for index in range(9500,10000):
    y = data[index] #get next obs to smooth
    post_mean, post_cov, loglike = kf.update(m = m, P = P, y=y.reshape(-1,1))
    smooth_mean, smooth_cov, smooth_gain = kf.smooth_current(m = m, P = P, ms=post_mean, Ps=post_cov)
    smoothed_current_obs.append(smooth_mean[0][0])
    m = post_mean
    P = post_cov  

df2 = pd.DataFrame([data[9500:10000],smoothed_current_obs]).T
df2.columns = ['data','smoothed_current_data']
df2.plot()

[oseiskar]

Hi,

Great to hear that you have found the library useful! Unfortunately, what you are attempting is not possible. This is not due to this library but how Kalman filters and smoothers work. Let me explain...

If you have a sequence of T values, [y₁, y₂, ..., y_T], the output of the Kalman filter at element t, 1 ≤ t ≤ T depends on the past values only, that is, y₁, y₂, ..., y_{t - 1}, y_t. The state of the Kalman filter can be efficiently updated online using only the previous state (m_t-1, P_t-1) and the new value y_t (using kf.update in this library).

Conversely, the value of the Kalman smoother (the output of kf.smooth) at sequence element t depends on all the past and the future values: y₁, y₂, ..., y_t, ..., y_T. It cannot be computed online becase every value in the sequence affects every output value of the Kalman smoother.

Here's an example (code follows the image) where the Kalman filter (dashed line) and smoother (solid line) for the same model have been applied to the full data with 150 elements (black) and the first 100 elements (blue).

A few things to notice:

• The output of the Kalman filter computed on the first 100 elements coincides with the first 100 elements of the Kalman filter computed on the full data (blue dashed line is exactly on top of the black dashed line)
• This is not the case for the Kalman smoother. The smoothed version of the first 100 elements differs from the first 100 elements of the smoothed full data
• Kalman filter and smoother outputs are always equal at the end of the available data (dashed and solid lines meet at the end)

import numpy as np
import simdkalman
import pandas as pd
import matplotlib.pyplot as plt

np.random.seed(0)
data = np.random.normal(size=150)

kf = simdkalman.KalmanFilter(
    state_transition = [[1,1],[0,1]],        # matrix A
    process_noise = np.diag([0.05, 0.002]),  # Q
    observation_model = np.array([[1,0]]),   # H
    observation_noise = 20.0)

kf_results = kf.smooth(data)

t = range(data.size)

def kf_smooth(y):
    return kf.smooth(y).observations.mean

def kf_filter(y):
    return kf.compute(y, 0, filtered=True).filtered.observations.mean

plt.plot(t, data, 'kx', alpha=0.2, label='data')
plt.plot(t[:100], data[:-50], 'bx', alpha=0.4, label='first 100')
plt.plot(t, kf_smooth(data), 'k-', label='smoothed full')
plt.plot(t, kf_filter(data), 'k--', label='filtered full')
plt.plot(t[:100], kf_smooth(data[:100]), 'b-', label='smoothed first 100')
plt.plot(t[:100], kf_filter(data[:100]), 'b--', label='filtered first 100')
plt.legend()
plt.show()

I hope this clarifies the problem.

BR, Otto