[fpc-devel] Numeric error in the calculation of Variance

[Werner Pamler]

Am 23.12.2017 um 13:24 schrieb Anton Shepelev:
> Hello, all
>
> May  I  ask you if there is any chance of fixing the
> error in the Math unit that I reported here:
>
>     https://bugs.freepascal.org/view.php?id=32804
>
> and provided a patch?  Here is a test case:
>
> Program Vartest;
> uses Math;
>
> const Size = 1000000;
> var dataS: array of Single;
>      dataD: array of Double;
>      dataE: array of Extended;
>      i,n: longint;
> begin
>    WriteLn('Each run should return a value near unity.');
>
>    WriteLn('Single:');
>    SetLength( dataS, Size );
>    for n := 1 to 4 do
>    begin
>      for i := 0 to Size - 1 do
>      begin
>        dataS[i] := 10000000 + RandG(0,1);
>      end;
>      WriteLn( Math.Variance( dataS ):5:3 );
>    end;
>
>    WriteLn('Double:');
>    SetLength( dataD, Size );
>    for n := 1 to 4 do
>    begin
>      for i := 0 to Size - 1 do
>      begin
>        dataD[i] := 1000000000000000 + RandG(0,1);
>      end;
>      WriteLn( Math.Variance( dataD ):5:3 );
>    end;
>
>    WriteLn('Extended:');
>    SetLength( dataE, Size );
>    for n := 1 to 4 do
>    begin
>      for i := 0 to Size - 1 do
>      begin
>        dataE[i] := 1000000000000000000 + RandG(0,1);
>      end;
>      WriteLn( Math.Variance( dataE ):5:3 );
>    end;
> end.

Yes, the formula used is  SS-Sqr(S)/N, and this is the bad one. 
Explained somewhere in the Numerical Recipes, but also in Widepedia 
(https://en.wikipedia.org/wiki/Variance#Formulae_for_the_variance): "it 
is not advised for computer calculations as it suffers from catastrophic 
cancellation <https://en.wikipedia.org/wiki/Catastrophic_cancellation> 
if the two components of the equation are similar in magnitude and 
floating point arithmetic is used."
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.freepascal.org/pipermail/fpc-devel/attachments/20171224/83b2c81c/attachment.html>