What happened?
Statistics.LinearMoments forms its probability-weighted-moment numerators in int, so they overflow for large samples and the method silently returns a wrong L-kurtosis.
In Numerics/Data/Statistics/Statistics.cs (lines 509-520 at 2a0357a), N is a double but i is an int, so these two products are integer arithmetic:
B2 += (i - 2) * (i - 1) / ((N - 2) * (N - 1)) * sortedData[i - 1];
B3 += (i - 3) * (i - 2) * (i - 1) / ((N - 3) * (N - 2) * (N - 1)) * sortedData[i - 1];
The triple product first exceeds int.MaxValue at i = 1293: 1290 * 1291 * 1292 = 2,151,683,880 against a ceiling of 2,147,483,647. The default unchecked context wraps it to a negative number rather than throwing, so B3 accumulates a wrong weight and τ4 comes back corrupt. There is no warning and no exception.
The pair product wraps too, at i = 46,343 (46,341 * 46,342 = 2,147,534,622; the last safe index is 46,342). λ1 and λ2 use no product and are unaffected, and τ3 is unaffected until 46,343, so for realistic sample sizes the error is confined to τ4.
What did you expect to happen?
τ4 to be correct for any sample length, or an exception rather than a silently wrong shape statistic. A Kappa-4 or GEV fit taken off these L-moments will consume the corrupt value as a real statistic.
1293 daily observations is under four years of record, so this is reachable in ordinary hydrologic use rather than only at extreme sizes.
Steps to reproduce
Call Statistics.LinearMoments on the evenly spaced sample x[i] = 1 + 0.5 * i, whose L-skewness and L-kurtosis are both exactly 0 at every length. Driving the library at 2a0357a:
| n |
τ3 |
τ4 |
| 1292 |
-3.1086244689504383E-15 |
-1.7763568394002505E-15 |
| 1293 |
-1.7763568394002505E-15 |
-0.18525251648817065 |
| 1300 |
-1.3766765505351941E-14 |
-1.446418581370934 |
τ4 is correct to rounding at 1292 and wrong at 1293, and the error grows with n.
Suggested fix
Form the numerators in double, matching what N already is:
double di = i;
B2 += (di - 2) * (di - 1) / ((N - 2) * (N - 1)) * sortedData[i - 1];
B3 += (di - 3) * (di - 2) * (di - 1) / ((N - 3) * (N - 2) * (N - 1)) * sortedData[i - 1];
Below the overflow these products are exact integers well under 2^53, so results are bit-identical to the current behavior for every sample that does not overflow. long would move the failure out to roughly 2.1 million points rather than removing it.
Affected namespace, class, or method
Numerics.Data.Statistics.Statistics.LinearMoments
Numerics version, package version, branch, or commit
2a0357a (v2.1.4).
Operating system and .NET details
Found while validating a C++ port of Numerics against the real library, and confirmed by driving the C# directly at the pin above. The same expression in C++ is a UBSan finding, since C++ leaves signed overflow undefined where C# defines the wrap: runtime error: signed integer overflow: 1665390 * 1292 cannot be represented in type 'int'. That is what led us to it.
What happened?
Statistics.LinearMomentsforms its probability-weighted-moment numerators inint, so they overflow for large samples and the method silently returns a wrong L-kurtosis.In
Numerics/Data/Statistics/Statistics.cs(lines 509-520 at2a0357a),Nis adoublebutiis anint, so these two products are integer arithmetic:The triple product first exceeds
int.MaxValueati = 1293:1290 * 1291 * 1292 = 2,151,683,880against a ceiling of2,147,483,647. The default unchecked context wraps it to a negative number rather than throwing, soB3accumulates a wrong weight and τ4 comes back corrupt. There is no warning and no exception.The pair product wraps too, at
i = 46,343(46,341 * 46,342 = 2,147,534,622; the last safe index is 46,342). λ1 and λ2 use no product and are unaffected, and τ3 is unaffected until 46,343, so for realistic sample sizes the error is confined to τ4.What did you expect to happen?
τ4 to be correct for any sample length, or an exception rather than a silently wrong shape statistic. A Kappa-4 or GEV fit taken off these L-moments will consume the corrupt value as a real statistic.
1293 daily observations is under four years of record, so this is reachable in ordinary hydrologic use rather than only at extreme sizes.
Steps to reproduce
Call
Statistics.LinearMomentson the evenly spaced samplex[i] = 1 + 0.5 * i, whose L-skewness and L-kurtosis are both exactly 0 at every length. Driving the library at2a0357a:-3.1086244689504383E-15-1.7763568394002505E-15-1.7763568394002505E-15-0.18525251648817065-1.3766765505351941E-14-1.446418581370934τ4 is correct to rounding at 1292 and wrong at 1293, and the error grows with n.
Suggested fix
Form the numerators in
double, matching whatNalready is:Below the overflow these products are exact integers well under 2^53, so results are bit-identical to the current behavior for every sample that does not overflow.
longwould move the failure out to roughly 2.1 million points rather than removing it.Affected namespace, class, or method
Numerics.Data.Statistics.Statistics.LinearMomentsNumerics version, package version, branch, or commit
2a0357a(v2.1.4).Operating system and .NET details
Found while validating a C++ port of Numerics against the real library, and confirmed by driving the C# directly at the pin above. The same expression in C++ is a UBSan finding, since C++ leaves signed overflow undefined where C# defines the wrap:
runtime error: signed integer overflow: 1665390 * 1292 cannot be represented in type 'int'. That is what led us to it.