mean, variance and standard deviation with math lite
Tracy Allen
Posts: 6,664
Catching up with the thread on memory requirements for averaging and ran across this comment from Heater:
Heater, is that a puzzle? There is a recursive algorithm in volume 2 of Knuth 'Art of Computer Programming" page 232, equations 15 and 16. it is part of a discussion of ugly errors that can occur in floating point, how the textbook method of calculating the standard deviation can come up needing the square root of a negative number. Here is the recursive algorithm as Propeller code. This version came out of Knuth's text (which can be quite intimidating!) as implemented by Michael Macdonald, a bright intern worked for me on a project that required multiple instances to track several variables and epochs.
Aside: How do you continuously calculate a "running" standard deviation of a signal using hardly any memory, no square roots and only one multiply per sample point and of course integer maths? That is a neat little trick I saw an 8085 based control system using back 1980 or so (it had previously been implemented as just a bunch of TTL chips).
Heater, is that a puzzle? There is a recursive algorithm in volume 2 of Knuth 'Art of Computer Programming" page 232, equations 15 and 16. it is part of a discussion of ugly errors that can occur in floating point, how the textbook method of calculating the standard deviation can come up needing the square root of a negative number. Here is the recursive algorithm as Propeller code. This version came out of Knuth's text (which can be quite intimidating!) as implemented by Michael Macdonald, a bright intern worked for me on a project that required multiple instances to track several variables and epochs.
PUB addNum(number)
'computes the current average value and the variance
'Method to calculate mean and variance, as shown in TAOCP Vol 2 Pg 232. implemented by Michael McDonald
'Only as accurate as fixed point math can do, which turns out to be fairly accurate actually, though values are still off by a little.
' "Proof" offered as an excel document.
' statistical variables are VARs in this object, separate set of VARS for each instance of the object.
n++ 'Increment the number of samples by one every time a value is added.
if(n == 1)
oldM := number
newM := number
oldS := 0
else
newM := oldM + (number - oldM)/n
newS := oldS + (number - oldM)*(number - newM)
oldM := newM
oldS := newS
PUB getAverage
' Returns the current average value of the data set.
return newM
PUB getVariance
'Returns the variance of the data in the set.
return newS/(n-1)
Comments
That is not how I saw it being done, with almost no maths, back in the late seventies.
Assume the incoming samples have a normal distribution.
Consider the positive half of that bell curve. One can mark a point on the horizontal axis, call it Sp, such that 13 out of every 14 sample values is below Sp and one in 14 is above Sp.
So we do he following:
For every new sample that is above Sp we increment Sp by one.
For every 13 new samples that are below Sp we decrement Sp by one.
Now imagine doing the same for the lower side of the bell curve. Producing a value Sm
Clearly Sp and Sm will drift up and down as new samples arrive. The mean value is (Sm + Sp) / 2. And (Sp - Sm) is some measure of deviation. Call it D.
Now, D is not the standard deviation but multiplying it by some magic number gets you the right value.
My this means we can arrive at the mean and standard deviation using nothing but comparators, counters and adders. This was implemented with good old TTL chips, there were no micro-processors at the time, for monitoring the SD of a high speed production process.
The down side of all this is that if the standard deviation of your samples suddenly increases your estimate will rise rapidly, but if your standard deviation suddenly decreases the estimate takes far longer to fall the desired value.
By the time I arrived on the scene the TLL had been replaced by an 8085 and the algorithm above implemented in software. I suspect that guys doing it knew nothing of Knuth, for sure I did not.
Ah, now I remember. They had done it by using a bunch of counter/timer chips. The incoming data was a measurement of density, the sensor produced a frequency proportional to density. So the A to D conversion was a counter chip. Basically they were hardware guys and had reimplemented the old TTL scheme in counter/timers! Besides that poor old 8085 could hardly keep up with the machine. They had to add a second one, with shared RAM. Then a third for the UI. Then a fourth for something else....
I'll have to puzzle about that one too. How can it test the 1:14 split without storing the history of samples? It is like non-parametric statistics, the median and the quartiles. It is easy if all the samples are stored and sorted in memory, but not so clear how to do it if you can't store and sort. You mention TTL logic counter/timer chips and a frequency proportional to density. Maybe the samples themselves were time intervals between events?
1) Every time a new sample is above some threshold increment the value of the threshold.
2) Count the samples that are below that threshold. When that count reaches 13 decrement the threshold value. AND reset the count to zero.
Do this for both a high and low threshold.
After a while it must be that thirteen out of fourteen must be between those two thresholds. From which a single multiply gets you the SD.
The sensors frequency output is of no consequence. Does the same job as an analog output and an A to D. Except that one has no A to D, one just needs to measure the frequency with a counter.
This was all running at about 40000 samples per minute, about 666 samples per second. 1.5ms per sample. Pretty tight for that old 8085 at about 40000 instructions per second!
1st sample Sk (k=1) is say 1000. I'll set Sp:=1000 and Sm:=1000. Both counters Cp~ and Cm~. Core snippet, sequential samples Sk?
if Sk>Sp Sp++ if Sk<Sp if ++Cp = 13 Cp~ Sp-- if Sk<Sm Sm-- if Sk>Sm if ++Cm = 13 Cm~ Sm++Is that on track? I'm not where I can try it. Just curious. The Knuth method does work, and it does not involve keeping the sum of the squares, but it does require a multiplication. I recall it was possible to come up with odd counterexamples where it disagreed with the standard definition, but I will take Knuth's word that it mostly converges to the standard definition.