Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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
share|improve this question
    
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

4 Answers 4

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
share|improve this answer
    
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.

share|improve this answer
    
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:

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.

share|improve this answer

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.