bugGNU Octave - Bugs: bug #53199, hist returns different results...

 
 

bug #53199: hist returns different results depending on 2nd 1:n/n input argument

Submitter:  -X- <jsh>
Submitted:  Wed 21 Feb 2018 06:17:00 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 4.2.1 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 01 Mar 2018 09:15:59 AM UTC, comment #9: 

@Rik: You said that the error is 8 eps, which indeed would be larger than expected and larger than necessary, but my point was that for floating-point numbers the discretization step grows with the magnitude of the number. In this case the number with the deviation is 16 (actually 15.9999...), and the deviation is 0.5 eps(16), that is, 1 eps(15.9999...), which is the best possible accuracy after the ideal result.

I did not know that Octave now can do scalar operations on ranges (sorry, Debian Stable here with Octave 4.0.3), but see it in the recent bug reports. Seemingly 4.2.1 on https://octave-online.net/ also cannot do it yet. Yes, of course it is sensible to do these calculations on ranges. But first I would not call the problem of the reporter a bug, and I would not be sure that with the range solution it is equivalent to Matlab in every possible situation without trying a very large number of possibilities. For instance, consider


[y,x]=hist(20*[rand(10,1);0;1],6)


You get a bin size of 3.333..., and you would ideally have a bin division at exactly 10. However, you generate a range of half-integers (that is, with an increment of 1 and a starting value of 0.5, both of which can be exactly represented), but then multiply with 3.333..., which cannot be exactly represented. Perhaps here the result differs depending on whether you do this multiplication on the range or on the corresponding vector (definitely there will be equivalent examples where they differ), and I could not say that Matlab will always give the result that corresponds to the operation on the range, because, as far as I know, Matlab does not have a range type. So the patch is fine, but I would not claim that it fixes a Matlab incompatibility in the general case.

Michael Leitner <mleitner>
Mon 26 Feb 2018 04:59:37 PM UTC, comment #8: 

I don't see how the analysis is wrong.  With every numeric calculation there is a chance for discretization error and rounding-off to the nearest floating point number.  In general, in a series of calculations, these are not statistically guaranteed to cancel out.  Instead, the error tends to grow by eps or a few eps for each calculation.

In this case, Octave has a specific Range data type which is designed for maximum accuracy in a string of serial additions such as


range(i) = base + increment*(i-1)


Replacing the Matrix type with a Range type is enough to get Matlab-equivalent results.


-      x = [0.5:n]'/n;
-      x = (max_val - min_val) * x  + min_val * ones (size (x));
+      ## Use range type to preserve accuracy
+      x = (0.5:n) * (1/n);
+      x = (max_val - min_val) * x + min_val;
+      x = x.';  # Convert to matrix;


This function is fairly old and hasn't received much attention over the years so I overhauled it in this changeset (http://hg.savannah.gnu.org/hgweb/octave/rev/1fa1869650cc).  The changeset fixes this bug, corrects handling of non-finite values in the auto-determination of bins, addresses problems with the validation of inputs, and adds BIST tests for all of these issues.

For the original reporter, this version of hist will be available in the 3.4.0 release.  However, since it is just an m-file you can retrieve it from Octave's Mercurial repository and just copy the file in to your distribution if you want to see the changes.


Rik <rik5>
Group administrator
Thu 22 Feb 2018 10:24:22 AM UTC, comment #7: 

I do not agree with this analysis. The relevant part of the code is lines 108--109 in hist.m. For n=55, min_val=0 and max_val=55 you would get in perfect arithmetic the range (0.5:55). With floating-point arithmetic you get deviations


((.5:55)/55)*55-(.5:55)


This is due to the division and subsequent multiplication by 55, which as intermediate result gives numbers that cannot be exactly expressed (the final result of course is expressible). I have deviations at entry 4, 8, 15 and 16. In all but the third case these are half the value of the corresponding epsilon eps(4), eps(8), eps(16), that is, the resulting number is e.g. 15.99999.... (which has an eps that is half the one of 16). So these are not large deviations, but have to be expected (in fact I am surprised that the result is always so accurate, in most cases even exact).

Of course you could change the order of expressions and fuse the division in line 108 and the multiplication in line 109 and only afterwards multiply with the range. Perhaps this is the actual code Matlab uses, so this would then give (on the same architecture!) results that are perfectly compatible with Matlab. It would not hurt to make this change. But actually, this problem is due to the user misunderstanding the arguments of hist() -- with a sample of integers from 0 to 55 you would of course want 56 bins, not 55. Using 55 bins lets all values fall exactly on bin boundaries, so you will be susceptible to rounding errors, no matter what you do.

Michael Leitner <mleitner>
Thu 22 Feb 2018 06:34:35 AM UTC, comment #6: 

This is due to round-off error.  If I subtract the bin centers generated by Matlab from those generated by Octave the difference is usually 0.  Occasionally, however, there are differences on the order of a few eps.  In the worst case, there is a difference of 8 eps or 1.8e-15.  It's pretty small, but that is the cause of the problem.

I think there must be a better way within hist to calculate the bin centers without getting so much round-off error.

Rik <rik5>
Group administrator
Thu 22 Feb 2018 12:58:25 AM UTC, comment #5: 

See attachment

(file #43382)

-X- <jsh>
Thu 22 Feb 2018 12:48:44 AM UTC, comment #4: 

Could you save the variables to a .mat file and upload them to this report?  I will look and see if I can detect anything.


%% Run code as before, and then
save 53199.mat a cnt1 ctr1 cnt2 ctr2



Rik <rik5>
Group administrator
Thu 22 Feb 2018 12:37:14 AM UTC, comment #3: 

When I look closer the bin centers seems to be the same, but Matlab will somehow use the bins differently (maybe simply due to rounding), so the final count is different.

-X- <jsh>
Thu 22 Feb 2018 12:30:25 AM UTC, comment #2: 

Tested in 2011/2017a and a 15.5 bin center does seem to be generated:


>> a = [  1,  2,  3,  4, 0;
       5,  4,  6,  7, 8;
       9, 12, 11, 10, 0;
      13, 16, 15, 14, 0;
      17, 20, 19, 18, 0;
      21, 22, 23,  2, 0;
      24, 27, 26, 25, 0;
      28, 31, 30, 29, 0;
      32, 35, 34, 33, 0;
      36, 39, 38, 37, 0;
      40, 43, 42, 41, 0;
      44, 47, 46, 45, 0;
      48, 51, 50, 49, 0;
      52, 55, 54, 53, 0];
a = a(:);
b = min (a)
n = max(a)
[cnt1, ctr1] = hist(a, 1:n);
ctr1
[cnt2, ctr2] = hist(a, n);
ctr2
isequal(cnt1,cnt2)

b =

     0


n =

    55


ctr1 =

  Columns 1 through 16

     1     2     3     4     5     6     7     8     9    10    11    12    13    14    15    16

  Columns 17 through 32

    17    18    19    20    21    22    23    24    25    26    27    28    29    30    31    32

  Columns 33 through 48

    33    34    35    36    37    38    39    40    41    42    43    44    45    46    47    48

  Columns 49 through 55

    49    50    51    52    53    54    55


ctr2 =

  Columns 1 through 9

    0.5000    1.5000    2.5000    3.5000    4.5000    5.5000    6.5000    7.5000    8.5000

  Columns 10 through 18

    9.5000   10.5000   11.5000   12.5000   13.5000   14.5000   15.5000   16.5000   17.5000

  Columns 19 through 27

   18.5000   19.5000   20.5000   21.5000   22.5000   23.5000   24.5000   25.5000   26.5000

  Columns 28 through 36

   27.5000   28.5000   29.5000   30.5000   31.5000   32.5000   33.5000   34.5000   35.5000

  Columns 37 through 45

   36.5000   37.5000   38.5000   39.5000   40.5000   41.5000   42.5000   43.5000   44.5000

  Columns 46 through 54

   45.5000   46.5000   47.5000   48.5000   49.5000   50.5000   51.5000   52.5000   53.5000

  Column 55

   54.5000


ans =

  logical

   1


-X- <jsh>
Wed 21 Feb 2018 05:04:49 PM UTC, comment #1: 

There is a difference between the two calling syntaxes.  When you call hist with a vector you directly set the bin centers.  When you call hist with just the number of bins to use it determines the bin centers based on the range of the data.  The two are not the same.  Use the second output argument to see the calculated bin centers.  Try the following code:


a = [  1,  2,  3,  4, 0;
       5,  4,  6,  7, 8;
       9, 12, 11, 10, 0;
      13, 16, 15, 14, 0;
      17, 20, 19, 18, 0;
      21, 22, 23,  2, 0;
      24, 27, 26, 25, 0;
      28, 31, 30, 29, 0;
      32, 35, 34, 33, 0;
      36, 39, 38, 37, 0;
      40, 43, 42, 41, 0;
      44, 47, 46, 45, 0;
      48, 51, 50, 49, 0;
      52, 55, 54, 53, 0];
a = a(:);
b = min (a)
n = max(a)
[cnt1, ctr1] = hist(a, 1:n);
ctr1
[cnt2, ctr2] = hist(a, n);
ctr2


In Octave this produces


octave:19> diary on
octave:20> a = [  1,  2,  3,  4, 0;
>        5,  4,  6,  7, 8;
>        9, 12, 11, 10, 0;
>       13, 16, 15, 14, 0;
>       17, 20, 19, 18, 0;
>       21, 22, 23,  2, 0;
>       24, 27, 26, 25, 0;
>       28, 31, 30, 29, 0;
>       32, 35, 34, 33, 0;
>       36, 39, 38, 37, 0;
>       40, 43, 42, 41, 0;
>       44, 47, 46, 45, 0;
>       48, 51, 50, 49, 0;
>       52, 55, 54, 53, 0];
octave:21> a = a(:);
octave:22> b = min (a)
b = 0
octave:23> n = max(a)
n =  55
octave:24> [cnt1, ctr1] = hist(a, 1:n);
octave:25> ctr1
ctr1 =

 Columns 1 through 19:

    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18   19

 Columns 20 through 38:

   20   21   22   23   24   25   26   27   28   29   30   31   32   33   34   35   36   37   38

 Columns 39 through 55:

   39   40   41   42   43   44   45   46   47   48   49   50   51   52   53   54   55

octave:26> [cnt2, ctr2] = hist(a, n);
octave:27> ctr2
ctr2 =

 Columns 1 through 8:

    0.50000    1.50000    2.50000    3.50000    4.50000    5.50000    6.50000    7.50000

 Columns 9 through 16:

    8.50000    9.50000   10.50000   11.50000   12.50000   13.50000   14.50000   15.50000

 Columns 17 through 24:

   16.50000   17.50000   18.50000   19.50000   20.50000   21.50000   22.50000   23.50000

 Columns 25 through 32:

   24.50000   25.50000   26.50000   27.50000   28.50000   29.50000   30.50000   31.50000

 Columns 33 through 40:

   32.50000   33.50000   34.50000   35.50000   36.50000   37.50000   38.50000   39.50000

 Columns 41 through 48:

   40.50000   41.50000   42.50000   43.50000   44.50000   45.50000   46.50000   47.50000

 Columns 49 through 55:

   48.50000   49.50000   50.50000   51.50000   52.50000   53.50000   54.50000

octave:28> diary off


Could you try the same code in Matlab and also report the version of Matlab you are using?  Octave used to be Matlab compatible, but I notice they are now trying to deprecate the hist() function and move to the histogram function instead.

Rik <rik5>
Group administrator
Wed 21 Feb 2018 06:17:00 AM UTC, original submission:  

In this case below hist returns different results depending if the second argument is 'n' or '1:n' where n is the number of count boxes. Specifically, if only the 'n' argument is supplied then the a=15 count is missing.


a = [  1,  2,  3,  4, 0;
       5,  4,  6,  7, 8;
       9, 12, 11, 10, 0;
      13, 16, 15, 14, 0;
      17, 20, 19, 18, 0;
      21, 22, 23,  2, 0;
      24, 27, 26, 25, 0;
      28, 31, 30, 29, 0;
      32, 35, 34, 33, 0;
      36, 39, 38, 37, 0;
      40, 43, 42, 41, 0;
      44, 47, 46, 45, 0;
      48, 51, 50, 49, 0;
      52, 55, 54, 53, 0];

n = max(a(:));
cnt1 = hist( a(:), 1:n )
cnt2 = hist( a(:),   n )
cnt1 - cnt2
assert( isequal(cnt1,cnt2) )


-X- <jsh>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #43382:  53199.mat added by jsh (602B - application/octet-stream - Added output)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mleitner (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by jsh (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-02-26 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2018-02-22 rik5 StatusWorks For Me Confirmed
    2018-02-22 jsh Attached File- Added 53199.mat, #43382
    2018-02-21 rik5 StatusNone Works For Me

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code