# 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. –  Martin Sleziak Jul 22 '12 at 16:27

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 Nov 25 '12 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. –  Rafal Golcz Dec 20 '13 at 9:02

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 Jul 22 '12 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. –  Ross Millikan Jul 22 '12 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 Nov 25 '12 at 5:32

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.

-

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. –  Joffan Feb 7 at 6:31
wait a minute, I am working on it. –  Wangyan Li Feb 7 at 6:48
@Joffan it's okay now, have a look at it –  Wangyan Li Feb 7 at 6:57
Thanks, much better –  Joffan Feb 7 at 7:15

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

-