Craynerd
Posts: 37
Joined: Mon Nov 25, 2013 9:09 pm

How to calculate moving average

Mon Feb 17, 2014 10:50 pm

Hi guys

I'm fairly new to python but have a script that is calculating the number of pulses every second - anywhere between 0 and 150.

However, a snapshot of one of these counts would be too random and they need smoothing with an average.

With my limited knowledge, I could count second one, count second two, ...three... Store then in separate variables and then every 30 seconds average all thirty variables, reset and start again. Ultimately, my output would be every 30 seconds for the average.

I've been told that there is a way to calculate a running/moving average, i.e 30 values stored, average taken. New input comes in every second, oldest value gets kicked out, new average taken. Consequently it is a smoother running average of the last 30 inputs and if needed the output can be taken every second.

Has anyone got a similar code or example they could send me? Or better still adapt my code!

Code: Select all

while True:
    time.sleep(1)
    count = wind_tally()   # brings through the new 1hz count
   print count - old-count   # prints latest count or could be count1 if averaging 
    old-count = count
In the code above, the print would become the new variable storage.
Count1 = count-old. ..next into count2... etc


Any help appreciated.

Chris

tenochtitlanuk
Posts: 156
Joined: Fri Jul 06, 2012 8:51 pm
Location: Taunton, Somerset, UK
Contact: Website

Re: How to calculate moving average

Mon Feb 17, 2014 11:13 pm

Simple way is to add 1/30 of the NEW value each time to 29/30 of the recorded current average.

For elegance you could take the first say 10 values and average them, and then start the loop described.
Of course such smoothing algorithms will adversely affect your ability to 'see' genuine sudden changes.

tenochtitlanuk
Posts: 156
Joined: Fri Jul 06, 2012 8:51 pm
Location: Taunton, Somerset, UK
Contact: Website

Re: How to calculate moving average

Mon Feb 17, 2014 11:21 pm

Example, throwing a totally random input in range 0<-->100, so the running average should settle near 50...
Note I use Python3, so 'print' uses brackets...

Code: Select all

import random, time

runAve =0

while True:
    time.sleep(1)
    count =random.randrange( 0, 100)
    print ( count, runAve)
    runAve =9/10*runAve +1/10*count

simplesi
Posts: 2327
Joined: Fri Feb 24, 2012 6:19 pm
Location: Euxton, Lancashire, UK
Contact: Website

Re: How to calculate moving average

Tue Feb 18, 2014 12:56 am

I sometimes (to avoid daft edge results) stick the string of values into an array - sort it and then average the middle values

(Bit like how the Olympic judges scoring is done - effectively throws away the high and low values)

I do this for Ultrasonic distance measurements where I might get a stupidly high or low value

Simon
Seeking help with Scratch and I/O stuff for Primary age children
http://cymplecy.wordpress.com/ @cymplecy on twitter

BudBennett
Posts: 89
Joined: Fri May 17, 2013 2:45 pm
Location: Westcliffe, Colorado, USA

Re: How to calculate moving average

Tue Feb 18, 2014 1:30 am

Try this:

Code: Select all

from collections import deque
class rolling_avg :
    ''' returns a simple rolling average of n most recent values '''

    def __init__(self, n=5,debug=False):
        "determine lengh of roll at instantiation"
        self.n = n
        self.xqueue = deque('')
        self.debug = debug
        
    def rolling_avg(self,x):
        # if the queue is empty then fill it with values of x
        if(self.xqueue == deque([])):
            for i in range(self.n):
                self.xqueue.append(x)
        self.xqueue.append(x)
        self.xqueue.popleft()
        avg = 0
        for i in self.xqueue:
            avg += i
        avg = avg/float(self.n)
        if self.debug:
            print("Rolling Avg:")
            for i in self.xqueue:
                print(i)
            print("avg: %f" % avg)
        return avg
Just instantiate it with the number of sample you need to average:

myRollingAvg = rolling_avg(n=5)

Bud

tenochtitlanuk
Posts: 156
Joined: Fri Jul 06, 2012 8:51 pm
Location: Taunton, Somerset, UK
Contact: Website

Re: How to calculate moving average

Tue Feb 18, 2014 1:54 pm

To a rather low-level Python programmer, 'deque' and 'instantiation' are frightening words! I've a bit to go yet, tho' the method looks interesting! ;)
I found this kind of smoothing to work very nicely as follows. It turns your screen into effectively a scrolling x/y display. Good practice in 'bit-slicing', and 'overloading' the '*' operator....
I use here a generated sin+noise function, but it works particularly nicely with the LDR+capacitor+GPIO to make an impressive display. ( 2 components and three or four lines of code to make a Pi have a simple AtoD!)
With short time intervals, show you have 4 (?) fingers by moving them in front of the LDR, or with much longer intervals graph the brightness in your room through a day or more. ( I used to do this printing to a dot-matrix printer- you could see a permanent record of when lights in my classroom went off and on- such as the early-morning cleaners. You can easily do this too in Python, of course.)
Note that smoothing something like a sin-wave+noise you'll get smaller amplitudes and a phase-shift.

Code: Select all

import math, random, time

x      = 0    # this will be an incrementing angle in degrees...
runAve =40    # as good as any to start mid-range! Better would be use first 'noisy' reading...

while x <=1.5 *360:
    time.sleep( 0.2)
    #   generate a 40-offset sine wave, amp 30, with random added noise of +-8 units.
    #   remember Python uses radians rather than degrees...
    count  =40 +30 *math.sin( x *math.pi /180) +16 *random.random() -8
    
    #   create a running average of 10% new & 90% old
    runAve = 9 /10 *runAve +1 /10 *count
    
    #   Construct a string of 80 spaces using the 'overloaded' multiply operator
    #   and modify it so 'C' ('count') and 'A' (average) are at the appropriate place
    #       using 'string-slicing'.
    op ="|" +" " *79
    op =op[ 0:int( count  -1)] +"C" +op[ int( count  +1):80]
    op =op[ 0:int( runAve -1)] +"A" +op[ int( runAve +1):80]
    print( op, "|", " ", x, " degrees.")
    
    x +=15
Output. 'A' -average, smoothed. 'C' =current noisy reading.

Code: Select all

"""
|                                      AC                                      |   0  degrees.
|                                        A               C                     |   30  degrees.
|                                         A             C                      |   60  degrees.
|                                           A                 C                |   90  degrees.
|                                              A                       C       |   120  degrees.
|                                                A          C                  |   150  degrees.
|                                C             A                               |   180  degrees.
|               C                           A                                  |   210  degrees.
|                 C                      A                                     |   240  degrees.
|  C                                 A                                         |   270  degrees.
|         C                        A                                           |   300  degrees.
|                              C  A                                            |   330  degrees.
|                                  A  C                                        |   360  degrees.
|                                   A              C                           |   390  degrees.
|                                      A                         C             |   420  degrees.
|                                         A                   C                |   450  degrees.
|                                          A              C                    |   480  degrees.
|                                           A       C                          |   510  degrees.
|                                  C       A                                   |   540  degrees.
|                          C              A                                    |   570  degrees.
|    C                                A                                        |   600  degrees.
|  C                               A                                           |   630  degrees.
|     C                         A                                              |   660  degrees.
|                              CA                                              |   690  degrees.
|                                A      C                                      |   720  degrees.
|                                 A          C                                 |   750  degrees.
|                                    A                         C               |   780  degrees.
|                                      A                      C                |   810  degrees.
|                                        A                C                    |   840  degrees.
|                                          A             C                     |   870  degrees.
|                                      C   A                                   |   900  degrees.
|                   C                   A                                      |   930  degrees.
|               C                     A                                        |   960  degrees.
|       C                          A                                           |   990  degrees.
|            C                   A                                             |   1020  degrees.
|                      C        A                                              |   1050  degrees.
|                               A                                              |   1080  degrees.
>>> 
"""

PiGraham
Posts: 4208
Joined: Fri Jun 07, 2013 12:37 pm
Location: Waterlooville

Re: How to calculate moving average

Tue Feb 18, 2014 2:32 pm

BudBennett wrote:Try this:

Code: Select all

from collections import deque
class rolling_avg :
    ''' returns a simple rolling average of n most recent values '''

    def __init__(self, n=5,debug=False):
        "determine lengh of roll at instantiation"
        self.n = n
        self.xqueue = deque('')
        self.debug = debug
        
    def rolling_avg(self,x):
        # if the queue is empty then fill it with values of x
        if(self.xqueue == deque([])):
            for i in range(self.n):
                self.xqueue.append(x)
        self.xqueue.append(x)
        self.xqueue.popleft()
        avg = 0
        for i in self.xqueue:
            avg += i
        avg = avg/float(self.n)
        if self.debug:
            print("Rolling Avg:")
            for i in self.xqueue:
                print(i)
            print("avg: %f" % avg)
        return avg
Just instantiate it with the number of sample you need to average:

myRollingAvg = rolling_avg(n=5)

Bud
It is more efficient to keep a sum value and add the new head value, subtract the old tail value and divide by the current length of the queue to get the average. No need to iterate over the queue or pre-fill it.
This is particularly beneficial for long queues since the calculation doesn't change with queue length. A thousand sample average takes the same time to update for a new sample as a 5 sample queue.

User avatar
jojopi
Posts: 3354
Joined: Tue Oct 11, 2011 8:38 pm

Re: How to calculate moving average

Tue Feb 18, 2014 4:00 pm

tenochtitlanuk wrote:Simple way is to add 1/30 of the NEW value each time to 29/30 of the recorded current average.
Not only is this the simplest method, but I would argue that it is one of the best. It is much more mathematically justified than it may at first seem. It is how the Linux kernel calculates the load average values in /proc/loadavg. It smooths anomalous samples but responds relatively quickly to genuine trends. It is an exponentially-weighted moving average: http://en.wikipedia.org/wiki/Moving_ave ... ng_average

There is a slight controversy over the selection of the weighting coefficient. With alpha = 1-exp(1-1/N), or the approximation 1/N as above, the last N samples account for 63.2% of the current average, which fits with the concept of the "time constant" in other fields. However, because all previous samples affect the result slightly, the average age of the contributions is N.

To make the average age agree with (N+1)/2, as in a simple N-point average, it is common to set alpha = 2/(N+1). Then the last N samples account for 86.5% of the current value.

Craynerd
Posts: 37
Joined: Mon Nov 25, 2013 9:09 pm

Re: How to calculate moving average

Tue Feb 18, 2014 6:32 pm

Well thank you for your replies and I have asked the questions and certainly got lots of suggestions and sample codes.

I admit, a lot of this has gone over my head.

As ultimately I'm integrating this code into a script that is already doing an awful lot (a weather station), I do need to try and select the most efficient way of doing it. Tenochtitlanuk and jojopi seem to be suggesting a method but I don't understand if any of these examples are demonstrating this logic. Could you please post an example of this? I'd really appreciate it.

Chris

jamesh
Raspberry Pi Engineer & Forum Moderator
Raspberry Pi Engineer & Forum Moderator
Posts: 27461
Joined: Sat Jul 30, 2011 7:41 pm

Re: How to calculate moving average

Tue Feb 18, 2014 6:46 pm

I think what they are saying is this (not python)

n = number samples so far
new = new sample

if (n = 0)
average = new
else
average = (average*(n-1))/n + new/n

n = n+1
if n>30 then n = 30
Principal Software Engineer at Raspberry Pi (Trading) Ltd.
Contrary to popular belief, humorous signatures are allowed.
I've been saying "Mucho" to my Spanish friend a lot more lately. It means a lot to him.

Craynerd
Posts: 37
Joined: Mon Nov 25, 2013 9:09 pm

Re: How to calculate moving average

Wed Feb 19, 2014 4:08 pm

Call me thick by all means but I still don't see how you are holding the numbed and excluding the last one!! I'd still appreciate help if anyone can lend some.
Chris

jamesh
Raspberry Pi Engineer & Forum Moderator
Raspberry Pi Engineer & Forum Moderator
Posts: 27461
Joined: Sat Jul 30, 2011 7:41 pm

Re: How to calculate moving average

Wed Feb 19, 2014 4:55 pm

It's not strictly an exact average of the last n entries, but an approximation. So you dont store any numbers, just keep the approximate average as you go along. I think.
Principal Software Engineer at Raspberry Pi (Trading) Ltd.
Contrary to popular belief, humorous signatures are allowed.
I've been saying "Mucho" to my Spanish friend a lot more lately. It means a lot to him.

User avatar
jojopi
Posts: 3354
Joined: Tue Oct 11, 2011 8:38 pm

Re: How to calculate moving average

Wed Feb 19, 2014 5:31 pm

An exponentially weighted average is an average of all the previous data points, but weighted so that the most recent values contribute the most, and the contributions of older and older data values decay exponentially.

Consider the simplest case where alpha = 1/2. Then the current average is 1/2 of the most recent reading, plus 1/4 of the previous reading, plus 1/8 of the reading before that, plus 1/16 of the reading before that, all the way back to the start. (Notice this adds to 1, so it is a true average.)

One of the advantages is that we do not need to retain the individual values. If we multiply the current average by 1/2, we have automatically shifted all the contributions back one place, and then we can add 1/2 of the next reading.

This works just as well with other alpha values. (For alpha=1/30, the contribution of the kth term is (1/30)*(29/30)**k, so multiply by 29/30 moves them down and we add 1/30 of the next reading.)

I do not think it is necessary or beneficial to change the rules at the beginning, except that we may want to use the first reading as a guess at the long-term average if there is no better guess available. If the initial averages are displaced as a result, then just drop the first few outputs.

In Python:

Code: Select all

N = 30
alpha = 1.0/N  # or: 2.0/(N+1)
average = wind_tally() # or: guess_at_long_term_average()
while True:
    count = wind_tally()
    average = alpha*count + (1-alpha)*average
    print("current = ", count, ", average = ", average)
    time.sleep(1)
Compared to a simple average of the last N values, the exponentially weighted average is more affected by the recent readings, but less affected by readings suddenly dropping out of the window. So it is arguably both more smooth and faster to respond to genuine trends.

Craynerd
Posts: 37
Joined: Mon Nov 25, 2013 9:09 pm

Re: How to calculate moving average

Wed Feb 19, 2014 7:47 pm

Ok, it makes sense what you are advising I do but maybe I should add more context. This is for a windSpeed measurement. Every 30 seconds or user set time, all data is exported to a website for processing, graphs and such. Do I really want an exponentially weighted average for this - if the wind has truly dropped do I not want to see this immediately and have this data processed.

Then again maybe not, as you said, to show the overall trend perhaps it should be the exponentially weighted average!!!

Chris

PiGraham
Posts: 4208
Joined: Fri Jun 07, 2013 12:37 pm
Location: Waterlooville

Re: How to calculate moving average

Wed Feb 19, 2014 9:46 pm

Craynerd wrote:Ok, it makes sense what you are advising I do but maybe I should add more context. This is for a windSpeed measurement. Every 30 seconds or user set time, all data is exported to a website for processing, graphs and such. Do I really want an exponentially weighted average for this - if the wind has truly dropped do I not want to see this immediately and have this data processed.

Then again maybe not, as you said, to show the overall trend perhaps it should be the exponentially weighted average!!!

Chris
If you want to see the full change immediately, don't filter it with a long time-constant filter such as these infinite impulse response (IIR) or 'exponentially weighted average' filters, or a moving average filter. The reason to filter is to reduce noise and reveal trends at some time scale.

You can tune the filter response time to your needs. An IIR with 50% feedback settles to a new value very quickly. A short-buffer moving average changes quickly.

IIR filters can be much more complex than the one described here, but I think you just need to display a current value that doesn't jump about too much and you can simply experiment until it seems right.

Ravenous
Posts: 1956
Joined: Fri Feb 24, 2012 1:01 pm
Location: UK

Re: How to calculate moving average

Thu Feb 20, 2014 11:38 am

Craynerd wrote:This is for a windSpeed measurement.
Hope butting in doesn't confuse anyone, but I think I read somewhere that the best way to report wind speed is a simple average (mean) over a two minute interval. (This was when I was into gliding, which was years ago!)

Also a single measurement every 30 secs is not enough to measure gusts - of course if the average is all you want that's fine.

Finally, what device is actually measuring the speed - a mechanical thing with a cup & cone or a windmill, or something else like direct pressure measurement or an ultrasound gizmo? Just out of interest, as a really fast responding system will respond differently to gusts.

Return to “Python”