浮点数学结果略有不同（从C到golang）

Question

I am working on developing a library of technical indicators directly in golang. It is, among other things, an exercise in learning golang.

I've been validating the results of my algorithms by building test cases with data generated with TA-Lib (or rather the ruby wrapper around TA-Lib).

This has been working fine until I got to the implementation of Bollinger Bands. My implementation seems to work fine, but differs at the 14-15th decimal place.

I've read Floating point math in different programming languages and suspect that this might be the culprit (I am doing the calculations in a slightly different order).

Edited to add:

The question above deals with a very simple manifestation of the floating point math. It's much harder to confirm that a longer piece of code is actually hitting this problem.

How can I confirm that it is just variations in floating point math because of the order?

/ End Edit

Am I correct in my understanding?

Here is my implementation:

package ta

import (
  "math"
)

func BollingerBands(values []float64, period int) ([]float64, []float64, []float64) {
  deviationsUp := 2.0
  deviationsDown := 2.0

  middleBand := Sma(values, period)
  offset := len(values)-len(middleBand)
  var upperBand []float64
  var lowerBand []float64
  for idx, v := range middleBand {
    backIdx := offset+idx-period+1
    curIdx := offset+idx+1
    if backIdx < 0 {
      backIdx = 0
    }
    stdDev := SliceStdDev(values[backIdx:curIdx])
    upperBand = append(upperBand, v + (stdDev * deviationsUp))
    lowerBand = append(lowerBand, v - (stdDev * deviationsDown))
  }
  return upperBand, middleBand, lowerBand
}

// Sma produces the Simple Moving Average for the
// supplied array of float64 values for a given period
func Sma(values []float64, period int) []float64{
  var result []float64
  for index,_ := range values {
    indexPlusOne := index+1
    if(indexPlusOne>=period) {
      avg := Mean(values[indexPlusOne-period:indexPlusOne])
      result = append(result, avg)
    }
  }
  return result
}

// SliceMean returns the Mean of the slice of float64
func SliceMean(values []float64) float64 {
  var total float64=0
    for _,element := range values {
        total += element
    }
  return total / float64(len(values))
}

// SliceVariance returns the variance of the slice of float64.
func SliceVariance(values []float64) float64 {
    if 0 == len(values) {
        return 0.0
    }
    m := SliceMean(values)
    var sum float64
    for _, v := range values {
        d := v - m
        sum += d * d
    }
    return sum / float64(len(values))
}

// SliceStdDev returns the standard deviation of the slice of float64.
func SliceStdDev(values []float64) float64 {
    return math.Sqrt(SliceVariance(values))
}

Which results in the following values for the upper band <[]float64 | len:6, cap:8>: [94.92564730599291, 94.50588827974477, 92.12752961253167, 101.58367006802706, 114.64331379078675, 120.58088881180322]

Using ruby:

require 'indicator/mixin'
x = [26.0, 54.0, 8.0, 77.0, 61.0, 39.0, 44.0, 91.0, 98.0, 17.0]
y = x.indicator(:bbands_5)
# {:out_real_upper_band=>[94.9256473059929, 94.50588827974477, 92.12752961253167, 101.58367006802709, 114.64331379078678, 120.58088881180323, nil, nil, nil, nil] <SNIP>}

Answer 1

I think the algorithms are different. For example variance:

/* Do the MA calculation using tight loops. */
/* Add-up the initial periods, except for the last value. */
periodTotal1 = 0;
periodTotal2 = 0;
trailingIdx = startIdx-nbInitialElementNeeded;

i=trailingIdx;
if( optInTimePeriod > 1 )
{
   while( i < startIdx ) {
      tempReal = inReal[i++];
      periodTotal1 += tempReal;
      tempReal *= tempReal;
      periodTotal2 += tempReal;
   }
}

/* Proceed with the calculation for the requested range.
 * Note that this algorithm allows the inReal and
 * outReal to be the same buffer.
 */
outIdx = 0;
do
{
   tempReal = inReal[i++];

   /* Square and add all the deviation over
    * the same periods.
    */

   periodTotal1 += tempReal;
   tempReal *= tempReal;
   periodTotal2 += tempReal;

   /* Square and add all the deviation over
    * the same period.
    */

   meanValue1 = periodTotal1 / optInTimePeriod;
   meanValue2 = periodTotal2 / optInTimePeriod;

   tempReal = inReal[trailingIdx++];
   periodTotal1 -= tempReal;
   tempReal *= tempReal;
   periodTotal2 -= tempReal;

   outReal[outIdx++] = meanValue2-meanValue1*meanValue1;
} while( i <= endIdx );

That doesn't look like your variance. If you were to reproduce the algorithms so that they did the exact same operations then the Go version should produce the same result. Go is just doing standard, IEEE 754 floating point arithmetic.

As to the question "does order matter?" It definitely does. Since floating point arithmetic is inexact you will lose information as you do the calculations. Most of the time it doesn't make much of a difference, but sometimes algorithms can be very susceptible to these changes. (so rearranging your formula algebraically may not lead to the same answer in real code)

You often find in libraries like these that algorithms have been designed to account for these issues and so they often don't look like the naive implementation. For example mean is usually a trivial function, but here's how its calculated in the GSL:

double
FUNCTION (gsl_stats, mean) (const BASE data[], const size_t stride, const size_t size)
{
  /* Compute the arithmetic mean of a dataset using the recurrence relation 
     mean_(n) = mean(n-1) + (data[n] - mean(n-1))/(n+1)   */

  long double mean = 0;
  size_t i;

  for (i = 0; i < size; i++)
    {
      mean += (data[i * stride] - mean) / (i + 1);
    }

  return mean;
}

So unless you match the algorithms exactly your answers will be subtly different. (which doesn't necessarily mean you're program is wrong)

One solution often used for this is to do equality comparisons within a very small number (math.Abs(expected-result) < ɛ, where you define ɛ: const ɛ = 0.0000001) rather than using ==.