bugGNU Scientific Library - Bugs: bug #38548, Rounding issues in gsl-histogram...

 
 

You are not allowed to post comments on this tracker with your current authentification level.

bug #38548: Rounding issues in gsl-histogram with integer numbers

Submitted by:  Christian Leitold <cleitold>
Submitted on:  Tue 19 Mar 2013 03:53:57 PM UTC  
 
Category: Accuracy problemSeverity: 3 - Normal
Operating System: Debian 7.0 amd64Status: None
Assigned to: NoneOpen/Closed: Open
Release: 1.15

Tue 19 Mar 2013 03:53:57 PM UTC, original submission:

A rounding problem occurs in gsl-histogram when you try to make a histogram from pure integer data, where it is appropriate to choose a bin size of 1. I have tracked this down to the file histogram/init.c, where the array range is filled with the bin limits. There, in make_uniform one has

double f1 = ((double) (n-i) / (double) n);
double f2 = ((double) i / (double) n);
range[i] = f1 * xmin + f2 * xmax;

In the case of e. g. xmin = 0.0, xmax = 60.0, n = 60, this (on my
system, gcc 4.6, standard make command from the file archive) leads to

range[31] = 31.000000000000003552713679

where it should be, without rounding error,

range[31] = 31.000000000000000000000000

This then leads to the number 31 being put into bin number 30, instead of 31. A bunch of uniformely distributed integer data would then lead to a peak of double the average height at 30, and no entry at 31. I have attached a sample list of integer data where the problem occurs (varlist.dat). As discussed in the corresponding thread [1] on the mailing list, a proposed solution is

double dx = (xmax-xmin) / (double) n;
double f1 = xmin + i*dx;
double f2 = xmax - (n-i) * dx;
range[i] = (f1 + f2) * 0.5;

As I understand it, in particular after reading [2], there are certain special cases where floating-point arithmetic is just as accurate and error-free as integer arithmetic. This should be fulfilled here, when xmin, xmax, n are real integer values (even though represented as doubles), and thus dx (=1), f1 and f2 are integer values as well.

[1] http://lists.gnu.org/archive/html/bug-gsl/2013-02/msg00006.html
[2] http://en.wikipedia.org/wiki/Double-precision_floating-point_format

Christian Leitold <cleitold>

 

Attached Files
file #27634:  varlist.dat added by cleitold (1kB - application/x-ns-proxy-autoconfig - Sample data, use e. g. xmin = 0, xmax = n = 60)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by cleitold (Submitted the item)
  •  

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follows 1 latest change.

    Date Changed By Updated Field Previous Value => Replaced By
    Tue 19 Mar 2013 03:53:57 PM UTCcleitoldAttached File-=>Added varlist.dat, #27634

    Back to the top


    Powered by Savane 3.1-cleanup