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

Question

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

• http://bilgin.esme.org/BitsBytes/KalmanFilterforDummies.aspx
• https://scipy-cookbook.readthedocs.io/items/KalmanFiltering.html

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

Answer 1

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

Answer 2

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

Answer 3

To give some flavor, see this list of books:

http://www.cs.unc.edu/~welch/kalman/kalmanBooks.html

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.

Answer 4

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.

Answer 5

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