[Hans Blom, 950927c]
(Bill Powers (950927.0645 MDT))
You get another WOW for that.
Thank you. Did you get my statement that it was _knowledge acqui-
sition_ and the subsequent _use_ of this knowledge that provides the
improvement in control?
But I guess I need some kindergarten lessons in how your system
works ... I realize that I'm asking you to teach Kalman Filter
theory (if that's what it is) from scratch, but I would really
appreciate being given a helping hand here.
What can I assume that you know already? Are you familiar with
statistics, in particular (normal or Gaussian) probabilities? I'll
assume little and give just the mechanics.
The basis of everything is that any variable x can be regarded as
having two components:
^ ~
x = x + x
where the first is the (known) expected value of x, also written as
E (x), and the second is the (unknown) uncertainty that is inherent
in our expectation. This second term will always be unknown, but we
might know (or estimate) its second moment or variance
E (x^2)
which is actually the square of the RMS error. In my program this
variance is called pxx. For example, if we know that a variable may
be anywhere between 4 and 18, we can take its expected value to be 11
and its standard deviation to be 7^2 = 49.
An example. Suppose that the model is
x = a * u + b
where a and b are unknown. Therefore x is unknown as well. We decom-
pose as follows:
^ ~ ^ ~ ^ ~
x + x = (a + a) * u + b + b
We do not decompose u, because it is fully known. Now take the "hat"-
and the non-"hat" parts separately:
^ ^ ^
x = a * u + b
~ ~ ~
x = a * u + b
The first gives the best prediction ("expected value") of x based on
the known value of u and the expected values of a and b; squaring the
second gives its expected variance:
~ ~ ~ ~ ~ ~ ~
E ( x^2) = E ((a * u + b)^2) = E (a^2 * u^2 + 2 * a * u * b + b^2)
This is simply plugging the expression for x in and squaring. Now
apply the E-operator to terms separately:
~ ~ ~ ~
= E (a^2 * u^2) + E (2 * a * u * b) + E (b^2)
and bring _known_ values (u and 2 are known) before the E-operator:
~ ~ ~ ~
= u^2 * E (a^2) + 2 * u * E (a * b) + E (b^2)
Now substitute names:
= u^2 * paa + 2 * u * pab + pbb
~ ~
using the notation E (x * y) = pxy, as the program does.
In the program, we assume that x, k and d are such random variables
that have to be decomposed. We have to follow the way in which the
expected values of x, k and d propagate, as well as their (co)vari-
ances pxx, pkk, pdd, pxk, pxd and pkd. Decomposition and the above
procedure gives the prediction part:
pxx := pkk * sqr (u) + 2.0 * pdk * u + pdd;
pxk := pkk * u + pdk;
pxd := pdk * u + pdd;
pkk := pkk + pmm;
pdd := pdd + pnn;
which you may now recognize. The terms pmm and pnn are taken to be
zero average normally distributed variables, so actually we would
have to write
~ ~
m = m, n = n
because the expected values of m and n are (assumed to be) zero. We
also know (or assume) that the random variables m and n are not
correlated with anything else (they are truly random), so terms like
pmx and pkn will always be zero and are therefore not found in the
program.
The next set of equations describes what happens when we have _two_
sources of information for _one_ variable. You may remember this from
Physics 101: given two measurements of the length of a stick with two
yardsticks of unequal but known precision, compute your best estimate
of the length of that stick and the precision of this estimate.
The second source of information in our case is, of course, the
observation. For x, for instance, the prediction yielded an expected
value and a pxx. But from the observation y we know that x = y +/- 0,
so the observation processing step gives x := y and pxx := 0. Things
get more complicated when we have two _inexact_ information sources.
In the general case, we have
^ ~ ^ ~
q = q1 + q1 = q2 + q2, for example 3 +/- 5 and 5 +/- 2
which we want to reduce to an equivalent
^ ~
q = q3 + q3, for example 4.5 +/- 1.7
This is what the second set of equations does. Try to derive the
formulas yourself. Note that the uncertainty decreases through the
observation.
k := k + (y - x) * pxk / pxx;
d := d + (y - x) * pxd / pxx;
pkk := pkk - sqr (pxk) / pxx;
pdd := pdd - sqr (pxd) / pxx;
pdk := pdk - pxd * pxk / pxx;
Maybe these are easier to understand when you rewrite them:
k * pxx := k * pxx + (y - x) * pxk
pkk * pxx := pkk * pxx - pxk * pxk
Finally, we do not know the values of pmm and pnn, but we can
estimate them from a running average. Since we modeled k as
k [i] := k [i-1] + m [i],
we have
pmm = E (m [i]^2) = E ((k [i] - k [i-1])^2) which can be
approximated to
pmm := pmm + (sqr (k - kold) - pmm) / 20.0; kold := k;
pnn := pnn + (sqr (d - dold) - pnn) / 20.0; dold := d;
Practice shows that a rough estimate works fine.
I realize that this is probably far too short to be really compre-
hensible, but the alternative would be, indeed, a course in Kalman
Filter theory. Which I won't do...
Greetings,
Hans