Text Only | Text with Attachments

QB64.org Forum

Active Forums => Programs => Topic started by: codeguy on October 11, 2018, 07:23:33 pm

Title: numerically stable average
Post by: codeguy on October 11, 2018, 07:23:33 pm
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
  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
  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&)
  27. a# = sum# / (finish - start + 1)
  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
  49. average = aa_int + aa_temp / StatN
  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.
Text Only | Text with Attachments