This is an automated email from the ASF dual-hosted git repository. aherbert pushed a commit to branch master in repository https://gitbox.apache.org/repos/asf/commons-statistics.git
commit 8afda419d31fbf065f2b1794b8d10e9a985c2f56 Author: Alex Herbert <[email protected]> AuthorDate: Tue Dec 26 08:24:50 2023 +0000 Remove use of Math.pow to compute higher-order moments A performance benchmark shows using pow(x, 3) and pow(x, 4) is an order of magnitude slower than compound multiplication. Multiplication will be at most 1 ULP error per term. --- .../commons/statistics/descriptive/Statistics.java | 3 +- .../descriptive/SumOfCubedDeviations.java | 15 ++++- .../descriptive/SumOfFourthDeviations.java | 18 ++++- .../statistics/descriptive/SkewnessTest.java | 3 +- .../jmh/descriptive/MomentPerformance.java | 78 +++++++++++++++++++++- 5 files changed, 109 insertions(+), 8 deletions(-) diff --git a/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Statistics.java b/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Statistics.java index f43218c..657ca7d 100644 --- a/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Statistics.java +++ b/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Statistics.java @@ -90,6 +90,7 @@ final class Statistics { // is smaller than the squared precision of the mean (m1). // Precision is set to 15 decimal digits // (1e-15 ~ 4.5 eps where eps = 2^-52). - return m2 <= Math.pow(1e-15 * m1, 2); + final double meanPrecision = 1e-15 * m1; + return m2 <= meanPrecision * meanPrecision; } } diff --git a/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/SumOfCubedDeviations.java b/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/SumOfCubedDeviations.java index bcc5565..fc1b69f 100644 --- a/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/SumOfCubedDeviations.java +++ b/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/SumOfCubedDeviations.java @@ -134,12 +134,23 @@ class SumOfCubedDeviations extends SumOfSquaredDeviations { // and opposite sign. So the sum-of-cubed deviations is zero. if (values.length > LENGTH_TWO) { for (final double x : values) { - s += Math.pow(x - xbar, 3); + s += pow3(x - xbar); } } return new SumOfCubedDeviations(s, ss); } + /** + * Compute {@code x^3}. + * Uses compound multiplication. + * + * @param x Value. + * @return x^3 + */ + private static double pow3(double x) { + return x * x * x; + } + /** * Updates the state of the statistic to reflect the addition of {@code value}. * @@ -212,7 +223,7 @@ class SumOfCubedDeviations extends SumOfSquaredDeviations { final double n1n2 = n1 + n2; final double dm = 2 * (halfDiffOfMean / n1n2); sumCubedDev += (sumSquaredDev * n2 - other.sumSquaredDev * n1) * dm * 3 + - (n2 - n1) * (n1 * n2) * Math.pow(dm, 3) * n1n2; + (n2 - n1) * (n1 * n2) * pow3(dm) * n1n2; } } } diff --git a/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/SumOfFourthDeviations.java b/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/SumOfFourthDeviations.java index 5e1a639..522aa34 100644 --- a/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/SumOfFourthDeviations.java +++ b/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/SumOfFourthDeviations.java @@ -117,11 +117,23 @@ class SumOfFourthDeviations extends SumOfCubedDeviations { // Note: This handles n=1. double s = 0; for (final double x : values) { - s += Math.pow(x - xbar, 4); + s += pow4(x - xbar); } return new SumOfFourthDeviations(s, sc); } + /** + * Compute {@code x^4}. + * Uses compound multiplication. + * + * @param x Value. + * @return x^4 + */ + private static double pow4(double x) { + final double x2 = x * x; + return x2 * x2; + } + /** * Updates the state of the statistic to reflect the addition of {@code value}. * @@ -195,7 +207,7 @@ class SumOfFourthDeviations extends SumOfCubedDeviations { sumFourthDev += (sumCubedDev - other.sumCubedDev) * halfDiffOfMean * 4 + (sumSquaredDev + other.sumSquaredDev) * (halfDiffOfMean * halfDiffOfMean) * 6 + - Math.pow(halfDiffOfMean, 4) * n1 * 2; + pow4(halfDiffOfMean) * n1 * 2; } else { final double n1n2 = n1 + n2; final double dm = 2 * (halfDiffOfMean / n1n2); @@ -204,7 +216,7 @@ class SumOfFourthDeviations extends SumOfCubedDeviations { sumFourthDev += (sumCubedDev * n2 - other.sumCubedDev * n1) * dm * 4 + (n2 * n2 * sumSquaredDev + n1 * n1 * other.sumSquaredDev) * (dm * dm) * 6 + - (n1 * n2) * (n1n2 * n1n2 - 3 * (n1 * n2)) * Math.pow(dm, 4) * n1n2; + (n1 * n2) * (n1n2 * n1n2 - 3 * (n1 * n2)) * pow4(dm) * n1n2; } } } diff --git a/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/SkewnessTest.java b/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/SkewnessTest.java index 125ba22..fcc550c 100644 --- a/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/SkewnessTest.java +++ b/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/SkewnessTest.java @@ -166,7 +166,8 @@ final class SkewnessTest extends BaseDoubleStatisticTest<Skewness> { // Matlab will return NaN for bias correction when n<=2. // This implementation matches the behaviour of Matlab. builder.accept(Arguments.of(new double[] {1, 2}, 0, nan, tol)); - builder.accept(Arguments.of(new double[] {1, 3, 7, 9, 11}, -0.15798755143759588, -0.23551393640880602, tol)); + builder.accept(Arguments.of(new double[] {1, 3, 7, 9, 11}, -0.15798755143759588, -0.23551393640880602, + DoubleTolerances.ulps(4))); return builder.build(); } } diff --git a/commons-statistics-examples/examples-jmh/src/main/java/org/apache/commons/statistics/examples/jmh/descriptive/MomentPerformance.java b/commons-statistics-examples/examples-jmh/src/main/java/org/apache/commons/statistics/examples/jmh/descriptive/MomentPerformance.java index ad55c5f..7b77433 100644 --- a/commons-statistics-examples/examples-jmh/src/main/java/org/apache/commons/statistics/examples/jmh/descriptive/MomentPerformance.java +++ b/commons-statistics-examples/examples-jmh/src/main/java/org/apache/commons/statistics/examples/jmh/descriptive/MomentPerformance.java @@ -145,9 +145,12 @@ public class MomentPerformance { public static class FunctionSource { /** Name of the source. */ @Param({MEAN, ROLLING_MEAN, SAFE_ROLLING_MEAN, SCALED_ROLLING_MEAN, - INLINE_SAFE_ROLLING_MEAN, INLINE_SAFE_ROLLING_MEAN_EXT + INLINE_SAFE_ROLLING_MEAN, INLINE_SAFE_ROLLING_MEAN_EXT, // Same speed as the ROLLING_MEAN, i.e. the DoubleConsumer is not an overhead //INLINE_ROLLING_MEAN + // Higher moments + "SumOfCubed", "SumOfCubedPow", + "SumOfFourth", "SumOfFourthPow", }) private String name; @@ -180,6 +183,14 @@ public class MomentPerformance { function = MomentPerformance::arrayInlineSafeRollingFirstMoment; } else if (INLINE_SAFE_ROLLING_MEAN_EXT.equals(name)) { function = MomentPerformance::arrayInlineSafeRollingFirstMomentExt; + } else if ("SumOfCubed".equals(name)) { + function = MomentPerformance::arraySumOfCubed; + } else if ("SumOfCubedPow".equals(name)) { + function = MomentPerformance::arraySumOfCubedPow; + } else if ("SumOfFourth".equals(name)) { + function = MomentPerformance::arraySumOfFourth; + } else if ("SumOfFourthPow".equals(name)) { + function = MomentPerformance::arraySumOfFourthPow; } else { throw new IllegalStateException("Unknown function: " + name); } @@ -492,6 +503,71 @@ public class MomentPerformance { return correctMeanKahan(data, m1); } + /** + * Create the sum-of-cubed deviations from the mean. + * + * @param data Data. + * @return the statistic + */ + static double arraySumOfCubed(double[] data) { + final double m = arrayInlineSafeRollingFirstMoment(data); + double s = 0; + for (final double x : data) { + final double dx = x - m; + s += dx * dx * dx; + } + return s; + } + + /** + * Create the sum-of-cubed deviations from the mean using the + * {@link Math#pow(double, double)} function. + * + * @param data Data. + * @return the statistic + */ + static double arraySumOfCubedPow(double[] data) { + final double m = arrayInlineSafeRollingFirstMoment(data); + double s = 0; + for (final double x : data) { + s += Math.pow(x - m, 3); + } + return s; + } + + /** + * Create the sum-of-fourth deviations from the mean. + * + * @param data Data. + * @return the statistic + */ + static double arraySumOfFourth(double[] data) { + final double m = arrayInlineSafeRollingFirstMoment(data); + double s = 0; + for (final double x : data) { + double dx = x - m; + dx *= dx; + s += dx * dx; + } + return s; + } + + /** + * Create the sum-of-fourth deviations from the mean using the + * {@link Math#pow(double, double)} function. + * + * @param data Data. + * @return the statistic + */ + static double arraySumOfFourthPow(double[] data) { + final double m = arrayInlineSafeRollingFirstMoment(data); + double s = 0; + for (final double x : data) { + s += Math.pow(x - m, 4); + } + return s; + } + /** * Create the mean from a stream of {@code double} values. *
