# Why use a Kalman filter instead of keeping a running average?

I've been trying to understand Kalman filters. Here are some examples that have helped me so far:

These use the algorithm to estimate some constant voltage. How could using a Kalman filter for this be better than just keeping a running average? Are these examples just oversimplified use cases of the filter?

(If so, what's an example where a running average doesn't suffice?)

EDIT:

For example, consider the following Java program and output. The Kalman output doesn't match the average, but they're very close. Why pick one over the other?

int N = 10; // Number of measurements

// measurements with mean = .5, sigma = .1;
double z[] =
{
0, // place holder to start indexes at 1
0.3708435, 0.4985331, 0.4652121, 0.6829262, 0.5011293,
0.3867151, 0.6391352, 0.5533676, 0.4013915, 0.5864200
};

double Q = .000001, // Process variance
R = .1*.1;// Estimation variance

double[] xhat = new double[N+1],// estimated true value (posteri)
xhat_prime = new double[N+1],   // estimated true value (priori)
p = new double[N+1],    // estimated error (posteri)
p_prime = new double[N+1],// estimated error (priori)
k = new double[N+1];    // kalman gain

double cur_ave = 0;

// Initial guesses
xhat[0] = 0;
p[0] = 1;

for(int i = 1; i <= N; i++) {
// time update
xhat_prime[i] = xhat[i-1];
p_prime[i] = p[i-1] + Q;

// measurement update
k[i] = p_prime[i]/(p_prime[i] + R);
xhat[i] = xhat_prime[i] + k[i]*(z[i] - xhat_prime[i]);
p[i] = (1-k[i])*p_prime[i];

// calculate running average
cur_ave = (cur_ave*(i-1) + z[i])/((double)i);

System.out.printf("%d\t%04f\t%04f\t%04f\n", i, z[i], xhat[i], cur_ave);
}


output:

 Iter      Input      Kalman     Average
1   0.370844    0.367172    0.370844
2   0.498533    0.432529    0.434688
3   0.465212    0.443389    0.444863
4   0.682926    0.503145    0.504379
5   0.501129    0.502742    0.503729
6   0.386715    0.483419    0.484227
7   0.639135    0.505661    0.506356
8   0.553368    0.511628    0.512233
9   0.401392    0.499365    0.499917
10  0.586420    0.508087    0.508567

• I've removed the filters tag, since it is intended for filters in set-theoretical and order-theoretical sense; see the tag description. Commented Jul 22, 2012 at 16:27
• Overview of Kalman filter in pictures bzarg.com/p/how-a-kalman-filter-works-in-pictures Commented Aug 18, 2015 at 4:08
• what if the system has dynamics?
– JMJ
Commented Apr 28, 2019 at 2:26

## 6 Answers

YES it is oversimplified example, more misleading than educating.

If so, what's an example where a running average doesn't suffice?

Any case when signal is changing.

Imagine moving vehicle. Calculating average means we assume signal value from any moment in time to be equally important. Obviously it is wrong. Intuition says, the last measurement is more reliable than the one from an hour before.

A very nice example to experiment with is of the form $\frac{1}{sT + 1}$. It has one state, so the equations won't get complicated.

In discrete time it could look like this:

x[n] = Ax[n-1] + Bu[n] + w[n]
y[n] = Cx[n] + v[n]
A = 0.99    B=1     C=1


There's the code that uses it (I'm sorry it's Matlab, I didn't use Python recently):

%% Initialize space
N = 100;               % nr of iterations
x = zeros(N,1);
u = zeros(N,1);
yv = zeros(N,1);

xprio = zeros(N,1); % a-priori     xk|k-1
xpost = zeros(N,1); % a-posteriori xk|k
Pprio = zeros(N,1);
Ppost = zeros(N,1);
K = zeros(N,1);

%------------------------ Variables to play with:
modelError = -0.04;     % relative model error
Q = 0.01;               % std. deviation of disturbance
R = 0.1;                % std. deviation of measurement noise
x(1) = 0.5;             % initial state of plant
xpost(1) = 1;           % initial estimate (state of Kalman filter)
Ppost(1) = 0.001;       % initial error estimate (state of Kalman filter)
%------------------------

% Plant
Areal = 0.99;
B = 1;
C = 1;
% Model of plant
Amodel = Areal*(1+modelError); % model never describes reality perfectly

% Generate noise
w = Q*randn(N,1);
v = R*randn(N,1);

%% Iterate
for k = 2:N
% simulate plant
x(k) = Areal*x(k-1) + B*u(k-1) + w(k);
% measurement
yv(k) = C*x(k) + v(k);

% prediction: predict current state from previous state and control
xprio(k) = Amodel*xpost(k-1)+B*u(k-1);
Pprio(k) = Amodel*Ppost(k-1)*Amodel' + Q;

% correction: use measurements with proper weight (K)
K(k) = Pprio(k)*C * inv(C*Pprio(k)*C' + R);
xpost(k) = xprio(k) + K(k)*(yv(k) - C*xprio(k));
Ppost(k) = (1 - K(k)*C)*Pprio(k);
end

%% Plot results
figure;
subplot(2,1,1);
plot(x,'k');
hold
plot(yv,'kx');
plot(xpost,'r');
legend('x real','x measure','x estimated');

% Important to see how K changes with time
subplot(2,1,2);
plot(K,'b')
legend('K');


There are some tips:

• Always set Q and R greater than zero.
Case $Q = 0$ is VERY BAD example. You say to the filter: "there is no disturbance acting on the plant", so after a while the filter will belief only to its predictions based on model rather than looking at measurements. Mathematically speaking $K_k \to 0$. As we know models don't describe reality perfectly.
• Experiment with some model inaccuracy - modelError
• Change initial guess of the state (xpost(1)) and see how fast it converges for different Q, R, and initial Ppost(1)
• Check how the filter gain K changes over time depending on Q and R
• You have answered part of my question. The other part is "How could using a Kalman filter for this be better than just keeping a running average?" I'm confused as to why a running average is better than a kalman filter in this particular situation when both are supposed to be optimal.
– Robz
Commented Nov 25, 2012 at 5:34
• I must admit I don't know the answer, but I guess that as long as we are talking about statistics (and Kalman Filter produces "a statistically optimal estimate"), we cannot say anything based on example with 10 samples. Commented Dec 20, 2013 at 9:02

In fact, they are the same thing in certain sense, I will show your something behind Kalman filter and you will be surprised.

Consider the following simplest problem of estimation. We are given a series of measurement $z_1, z_2, \cdots, z_k$, of an unknown constant $x$. We assume the additive model \begin{eqnarray} z_i= x + v_i, \; i=1,2, \cdots, k ~~~~~~~~~~~ (1) \end{eqnarray} where $v_i$ are measurement noises. If nothing else is known, then everyone will agree that a reasonable estimate of $x$ given the $k$ measurements can be given by \begin{eqnarray} \hat{x}_k= \frac{1}{k} \sum_{i=1}^{k} z_i ~~~~~~~~~~~ ~~~~~~~~~~~ (2) \end{eqnarray} this is average.

Now we can re-write above eq.(2) by simple algebraic manipulation to get \begin{eqnarray} \hat{x}_k= \hat{x}_{k-1} + \frac{1}{k} (z_k-\hat{x}_{k-1}) ~~~~~~~~~~~ (3) \end{eqnarray} Eq.(3) which is simply Eq.(2) expressed in recursive form has an interesting interpretation. It says that the best estimate of $x$ after $k$ measurement is the best estimate of $x$ after $k-1$ measurements plus a correction term. The correction term is the difference between what you expect to measure based on $k-1$ measurement, i.e., and what you actually measure $z_k$.

If we label the correction $\frac{1}{k}$ as $P_k$, then again simply algebraic manipulation can write the recursive form of $P_k$ as \begin{eqnarray} P_k=P_{k-1}-P_{k-1}(P_{k-1}+1)^{-1}P_{k-1} ~~~~~~~~~~~ (4) \end{eqnarray}

Believe it or not, Eqs.(3-4) can be recognized as the Kalman filtering equations for this simple case.

Any discussion is welcomed.

Reference:

Explaining Filtering (Estimation) in One Hour, Ten Minutes, One Minute, and One Sentence by Yu-Chi Ho

• Can you pull the text from your slides into a rather more page-friendly format? thanks. Commented Feb 7, 2015 at 6:31
• wait a minute, I am working on it. Commented Feb 7, 2015 at 6:48
• @Joffan it's okay now, have a look at it Commented Feb 7, 2015 at 6:57
• Arriving at eq. 4 would be slightly clearer if an intermediate step showed some of the manipulation, like how 1/k becomes ^-1 with some terms reordered. Commented May 16, 2018 at 21:52
• Replace P_k with 1/k and P_{k-1} with 1/(k-1) in (4). Then you will find the rationale. Commented May 17, 2018 at 5:29

To give some flavor, see this list of books:

I have Grewal+Andrews with MatLab, also Grewal+Weill+Andrews about GPS.

That is the fundamental example, GPS. Here is a simplified example, I interviewed for a job where they were writing software for keeping track of all trucks going in and out of a huge delivery yard, for Walmart or the like. They had two types of information: based on putting an RFID device in each truck, they had pretty good information about the direction each truck was going with measurements possible many times per second, but eventually growing in error, as does any essentially ODE approximation. On a much longer time scale, they could take the GPS position of a truck, which gives a very good unbiased location but has a large variance, you get position within 100 meters or something. How to combine these? That's the main use of the Kalman filter, when you have two sources of information giving roughly opposite types of error. My idea, which i would have told them if they had paid me, was to place a device on each semi where the cab meets the trailer, giving the current turning radius. This could have been integrated to give very good short-time information about the direction the truck was heading.

Well, that is what they do with almost anything moving nowadays. The one I thought was cute was farms in India, keeping track of where tractors were. The moving body does not need to be moving rapidly to bring about the same questions. But, of course, the first major use was the NASA Apollo project...My father met Kalman at some point. Dad worked mostly on navigation, initially missiles for the Army, later submarines for the Navy.

A running average is one kind of Kalman filter. Following the notation in your first link $\hat{X}_k=K_kZ_k+(1-K_k)\hat{X}_{k-1}$, a running average sets $K_k=\frac 1k$. If your underlying model is that the parameter of interest doesn't change with time, it is what you get. Other forms are needed if $X$ changes with time.

• I don't see it. K isn't a function of k in any equations for a Kalman filter that I've seen. It's usually something like "p/(p+R)". How is "1/k" what you get? And if the Kalman filter is optimal, why wouldn't it converge to same as the average?
– Robz
Commented Jul 22, 2012 at 19:17
• @Robz: Setting $K_k=1/k$ will give the running average, as it weights each new point by $\frac 1k$ and decreases the weight of all the previous by a factor $\frac {k-1}k$. $K$ is allowed to be a function of time, which is why the subscript. Many models do not use this flexibility, but some do. If the parameter is truly constant, an average of all the data is optimal, but if it is changing, you want to weight the newer data more. Commented Jul 22, 2012 at 19:32
• "K is allowed to be a function of time"--I still don't see it. Looking at any equations anywhere about kalman filters, K is never an explicit function of time. As I understand them, nothing in the kalman filter equations are a function of time, except the transition matrices which can depend on delta time between iterations. The crux of my question is on this contradiction: (1) kalman filters are optimal estimators for linear systems (2) the system I describe here is linear (3) taking an average is optimal in this system (4) the average and the kalman filter do not produce the same results.
– Robz
Commented Nov 25, 2012 at 5:32

As mentioned by a previous poster, you can use the following Kalman filter to implement a running average:

$\hat{X}_k=K_kZ_k+(1-K_k)\hat{X}_{k-1}$,

where $k$ runs from 1 to $N-1$. The discrepancy you observe stems from the fact that you don't use the measurement of $Z_0$ in your calculation. The Kalman filter gives you the same value for the average if you compute the average of $Z$ for $k=1..N-1$, that is, leaving the first measurement out. Alternatively you can do one more iteration by upping $k$ by one, but using $Z_0$ (as $Z_{N}$ does not exist).

Hope this helps.

Lucas

what's an example where a running average doesn't suffice?

Such example can be seen with this series of random 1D positions, created with a constant mean velocity and some instantaneous variance, represented by the plain curve:

These values are all pseudo-random numbers. Imagine the position is reported by a sensor which accuracy is a Gaussian distribution. The measurements are indicated by crosses on the figure above.

Looking at values between epochs 8 and 18, or after epoch 30, we see successive measurements repeatedly smaller or larger than actual values. A running average will be inaccurate for these ranges.

On the other hand, a Kalman filter is able to correct for these situations to some extent. However there is no ready-to-use Kalman filter, a filter must be tuned.

It must be provided with a model for prediction (here it will be a constant velocity model) and it must be tuned for the variance of the input. When tuned, the result is clearly in favor of the Kalman filter:

The filter has no access to the actual values of the continuous curve, still it follows them quite accurately.

How could using a Kalman filter for this be better than just keeping a running average?

There is a significant difference between a running average and a Kalman filter:

• The running average is a convolution with a fixed window and fixed weights.

• A Kalman filter is a dynamic implementation of Bayes theorem.

The Kalman filter above tracks two variables: The position, which is observed by the measurements, and the velocity, which is hidden, as there are no related measurements. Velocity is correlated to position, that's how the filter can track it.

A Kalman filter assumes all variables have Gaussian distributions. It might be far from the truth in certain cases, but it is useful in most cases, like in our case where all numbers are drawn from a normal distribution.

The power of a Kalman filter lies in its ability to combine multiple variables, determine a covariance between them and use this covariance to reduce the variance of the result. There is nothing comparable in the unvariate average.

This bivariate filter is provided a model in the form of a state transition matrix corresponding to a constant velocity motion: $$F = \begin{pmatrix} 1 & dt \\ 0 & 1 \end{pmatrix}$$

used to multiply the bivariate state (position and velocity at epoch k) and predict the next state (position and velocity at epoch k+1). A covariance matrix of the system is also tracked and updated by comparing the prediction and the measurement at each epoch.

The ratio between variance of the filter and variance of the measurement is used to take a weighted average of the filter prediction and the measurement provided. E.g. the filter is able to depart from the measurements after epoch 30 because the confidence in the prediction has been gained before this epoch.

The result from the filter and the running average greatly diverges due to this learning curve not present in the running average.

However the filter confidence in prediction must be moderated to prevent it from constantly diverging from the measurements. This is the role of the covariance matrix Q, This filter has been tuned explicitly for the measurements:

$$Q = \begin{pmatrix} 0.2 & 0.4 \\ 0.4 & 0.8 \end{pmatrix}$$

This matrix is added to the predicted covariance matrix. For the update phase, a coefficient of R = 50 is used when determining the filter gain K.

This tuning holds when more epochs are taken into account, e.g. with 500 samples, the mean error is the same:

As an example, here is result with different values for Q:

The colored band shows the position variance of the system covariance matrix. The thinner the band, the larger the confidence in predictions. Increasing Q has the effect of increasing the system covariance. As the system trusts less and less the prediction it trusts more and more the measurements and the weighted mean starts moving closer to the black curve, q=0.8 is the optimum:

Below is the Python code used for this simulation. It uses FilterPy from Roger Labbe to create the Q noise matrix, but the other equations for a Kalman filter are included in the program code.

import sys; sys.path.append('lib')
import numpy as np
import scipy.ndimage as sn
import matplotlib.pyplot as plt
from filterpy.common import Q_discrete_white_noise

#%% Helpers

# Create path with mean velocity and velocity variance
def rnd_walk(n, x0, *, dt=1, v_mean=1, v_var=1, seed=None):
rng = np.random.default_rng(seed)
std = np.sqrt(v_var)
x = x0
xs = [x]
for _ in range(n-1):
v = rng.normal(loc=v_mean, scale=std)
x += v * dt
xs.append(x)
return np.array(xs)

# Create a series approximating another series, by some variance
def sense(xs, var, *, seed=None):
rng = np.random.default_rng(seed)
std = np.sqrt(var)
ys = rng.normal(loc=xs, scale=std)
return ys

# Create a figure
def figure(mosaic, **kw):
mosaic = np.asarray(mosaic)
kw2 = dict(layout='constrained')
kwargs = kw2| kw
fig, axs = plt.subplot_mosaic(mosaic, **kwargs)
return fig, axs

#%% Kalman Filter

class kalman():
def __init__(self, x, P, /, F, H, R, Q):
# State
self.x = x # position and velocity
self.P = P # system covariance
self.F = F # state transition function(innovation)
self.H = H # measurement function
self.R = R # measurement covariance
self.Q = Q # process (prediction) noise

# Prediction function
def predict(self):
self.x = self.F @ self.x
self.P = self.F @ self.P @ self.F.T + self.Q
return self.x, self.P

# Update function
def update(self, z):
if z is None: return self.x, self.P

# Residual and gain
S = self.H @ self.P @ self.H.T + self.R
K = self.P @ self.H.T @ np.linalg.inv(S)
y = z - self.H @ self.x

# Transition
self.x = self.x + K @ y
self.P = self.P - K @ self.H @ self.P
return self.x, self.P

# Run predict-update cycles for measurements in zs
def run(self, zs):
# History arrays
ph = [self.x.copy()]
fh = [self.x.copy()]
Ph = [self.P.copy()]

# Run for the current epoch
for z in zs[:-1]:
# Predict and update state
p, self.P = self.predict()
self.x, self.P = self.update(z)

ph.append(p.copy())
fh.append(self.x.copy())
Ph.append(self.P.copy())

ph = np.array(ph)
fh = np.array(fh)
Ph = np.array(Ph)

return ph, fh, Ph

#%% Simulation

n = 500 # number of steps
dt = 1 # time between steps
seed = 0 # seed for random generator

# Actual track
x0 = 10 # initial position
v_mean = 1 # mean velocity
v_var = 10 # velocity variance

# Variance of measurements
m_var = 50

# Plot parameters
k = ['k', 'r']
titles = ['Kalman', 'Running average']
colors = ['magenta', 'dodgerblue']

#%% Kalman filter parameters

# Values assumed to initialize the Kalman filter
a_x0 = 10 # initial position
a_v0 = 0 # initial velocity
a_x0_var = 50 # variance of position
a_v0_var = 50 # variance of velocity

# Kalman variable parameters
q_var = 0.8 # for prediction noise Q
z_var = 50 # For measurement noise R

#%% Running average parameters

L = 2 # size of the window

#%% Initialization

# Create actual track and measurements
xs = rnd_walk(n, x0, dt=dt, v_mean=v_mean*dt, v_var=v_var, seed=seed)
zs = sense(xs, var=m_var, seed=seed)
zs_e = zs-xs
zs_v = 2*zs_e.std()

epoch = np.arange(n)

# Create figure
mosaic = [*[[kk] for kk in k], ['e']]
h_ratios = [1, 1, 0.6]
kw = dict(figsize=(8, 3*sum(h_ratios)), height_ratios=h_ratios)
fig, axs = figure(mosaic, **kw)

# Create main subplots, plot measurements
for i, (axid, title) in enumerate(zip(k, titles)):
ax = axs[axid]
ax.set_title(title, loc='left')
ax.set_xlabel('epoch')
ax.set_ylabel('position')
ax.grid()
if i==0: ax_ref = ax
else: ax.sharey(ax_ref)
ax.plot(xs, label='actual', c='k')
ax.scatter(epoch, zs, label='measured', c='k', marker='+', s=50, alpha=0.3)

# Create error subplot
ax = axs['e']
ax.set_title('Error vs. actual')
ax.set_xlabel('epoch')
ax.set_ylabel('error')
ax.axhline(0)
ax.grid()

#%% Filter measurements

# Kalman filter parameters
F = np.array([[1, dt], [0, 1]]) # state transition matrix
Q = Q_discrete_white_noise(dim=2, dt=dt, var=q_var)
H = np.array([[1., 0.]]) # to convert from state to measurement space
R = np.array([[z_var]])

# Kalman initial state
x = np.array([a_x0, a_v0]) # initial state (position and velocity)
P = np.diag([a_x0_var, a_v0_var]) # uncertainty on initial state

# Filter measurements
kf = kalman(x, P, F=F, H=H, R=R, Q=Q)
ph, fh, Ph = kf.run(zs)
ks = fh[:,0]

# Running average
rs = sn.uniform_filter1d(zs, L, mode='nearest', origin=(-L // 2))

# Compute error and its variance
ks_e = ks-xs
rs_e = rs-xs
ks_v = 2*ks_e.std()
rs_v = 2*rs_e.std()

# Plot results
for i, (s, e, v) in enumerate(zip([ks, rs],
[ks_e, rs_e],
[ks_v, rs_v])):

ax = axs[k[i]]
ax.plot(s, c=colors[i], lw=3)

# Plot errors
ax = axs['e']
ax.plot(epoch, e, label=f'{titles[i]} ($$\\sigma^2$$={v:.1f})', c=colors[i])
ax.legend()