numerically stable average

[codeguy]

In this post are two methods for calculating the average of a dataset. The standard method proves to be problematic for some cases, where the element values are large or the range is large. These two algorithms agree UNTIL n becomes large, then the truncation errors start creeping in for the standard method, using a running sum of all the elements in the range. The numerically stable method avoids this error quite well. In this particular instance, the numerically stable version lags the standard method by 50% or so, but the fact that overflow is completely avoided may be of benefit in some cases when dealing with large datasets or numbers that are large or perhaps even both.

Code: QB64: [Select]

1. ntest& = 0
2. DO
3.     ntest& = ntest& + 1
4.     REDIM a(0 TO ntest&) AS DOUBLE
5.     FOR q& = LBOUND(a) TO UBOUND(a)
6.         a(q&) = ntest& / 2 + .125
7.     NEXT
8.     s0! = TIMER(.001)
9.     AverageArrayNS a(), LBOUND(a), ntest&, average#
10.     f0! = TIMER(.001)
11.     PRINT "numerically stable "; ntest&; average#; f0! - s0!
12.     s1! = TIMER(.001)
13.     AverageStandardX a(), LBOUND(a), ntest&, averagx#
14.     f1! = TIMER(.001)
15.     PRINT "standard not stable"; ntest&; averagx#; f1! - s1!
16.     ntest& = ntest& * 2
17. LOOP
18.  
19. '****************
20. '* this method is subject to truncation/overflow errors -- fails at large N.
21. '****************
22. SUB AverageStandardX (a() AS DOUBLE, start AS LONG, finish AS LONG, a#)
23. sum# = 0
24. FOR c& = start TO finish
25.     sum# = sum# + a(c&)
26. NEXT
27. a# = sum# / (finish - start + 1)
28. END SUB
29.  
30. '****************
31. '* a numerically stable way to calculate the average of elements in an array without the dreaded overflow.
32. '****************
33. SUB AverageArrayNS (a() AS DOUBLE, start AS LONG, finish AS LONG, average AS DOUBLE)
34. DIM aa_int AS DOUBLE: aa_int = 0
35. DIM aa_temp AS DOUBLE: aa_temp = 0
36. DIM StatN AS LONG: StatN = finish - start + 1
37.  
38. FOR c& = start TO finish
39.     aa_temp = a(c&) + aa_temp
40.     DO
41.         IF aa_temp < StatN THEN
42.             EXIT DO
43.         ELSE
44.             aa_int = aa_int + 1
45.             aa_temp = aa_temp - StatN
46.         END IF
47.     LOOP
48. NEXT
49. average = aa_int + aa_temp / StatN
50. END SUB
51.  

Just another perhaps boring but useful algo to add to your code arsenal if overflow is to be avoided at all costs and results MUST be accurate.