StatisticUtils.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.statistics.inference;

  19. import java.util.Arrays;
  20. import java.util.Collection;
  21. import org.apache.commons.numbers.core.DD;
  22. import org.apache.commons.numbers.core.Precision;
  23. import org.apache.commons.numbers.core.Sum;
  24. import org.apache.commons.statistics.descriptive.Mean;

  26. /**
  27.  * Utility computation methods.
  28.  *
  29.  * @since 1.1
  30.  */
  31. final class StatisticUtils {
  32.     /** No instances. */
  33.     private StatisticUtils() {}

  35.     /**
  36.      * Compute {@code x - y}.
  37.      *
  38.      * <p>If {@code y} is zero the original array is returned, else a new array is created
  39.      * with the difference.
  40.      *
  41.      * @param x Array.
  42.      * @param y Value.
  43.      * @return x - y
  44.      * @throws NullPointerException if {@code x} is null and {@code y} is non-zero
  45.      */
  46.     static double[] subtract(double[] x, double y) {
  47.         return y == 0 ? x : Arrays.stream(x).map(v -> v - y).toArray();
  48.     }

  50.     /**
  51.      * Compute the degrees of freedom as {@code n - 1 - m}.
  52.      *
  53.      * <p>This method is common functionality shared between the Chi-square test and
  54.      * G-test. The pre-conditions for those tests are performed by this method.
  55.      *
  56.      * @param n Number of observations.
  57.      * @param m Adjustment (assumed to be positive).
  58.      * @return the degrees of freedom
  59.      * @throws IllegalArgumentException if the degrees of freedom is not strictly positive
  60.      */
  61.     static int computeDegreesOfFreedom(int n, int m) {
  62.         final int df = n - 1 - m;
  63.         if (df <= 0) {
  64.             throw new InferenceException("Invalid degrees of freedom: " + df);
  65.         }
  66.         return df;
  67.     }

  69.     /**
  70.      * Gets the ratio between the sum of the observed and expected values.
  71.      * The ratio can be used to scale the expected values to have the same sum
  72.      * as the observed values:
  73.      *
  74.      * <pre>
  75.      * sum(o) = sum(e * ratio)
  76.      * </pre>
  77.      *
  78.      * <p>This method is common functionality shared between the Chi-square test and
  79.      * G-test. The pre-conditions for those tests are performed by this method.
  80.      *
  81.      * @param expected Expected values.
  82.      * @param observed Observed values.
  83.      * @return the ratio
  84.      * @throws IllegalArgumentException if the sample size is less than 2; the array
  85.      * sizes do not match; {@code expected} has entries that are not strictly
  86.      * positive; {@code observed} has negative entries; or all the observations are zero.
  87.      */
  88.     static double computeRatio(double[] expected, long[] observed) {
  89.         Arguments.checkValuesRequiredSize(expected.length, 2);
  90.         Arguments.checkValuesSizeMatch(expected.length, observed.length);
  91.         Arguments.checkStrictlyPositive(expected);
  92.         Arguments.checkNonNegative(observed);
  93.         DD e = DD.ZERO;
  94.         DD o = DD.ZERO;
  95.         for (int i = 0; i < observed.length; i++) {
  96.             e = e.add(expected[i]);
  97.             o = add(o, observed[i]);
  98.         }
  99.         if (o.doubleValue() == 0) {
  100.             throw new InferenceException(InferenceException.NO_DATA);
  101.         }
  102.         // sum(o) / sum(e)
  103.         final double ratio = o.divide(e).doubleValue();
  104.         // Allow a sum within 1 ulp of 1.0
  105.         return Precision.equals(ratio, 1.0, 0) ? 1.0 : ratio;
  106.     }

  108.     /**
  109.      * Adds the value to the sum.
  110.      *
  111.      * @param sum Sum.
  112.      * @param v Value.
  113.      * @return the new sum
  114.      */
  115.     private static DD add(DD sum, long v) {
  116.         // The condition here is a high probability branch if the sample is
  117.         // frequency counts which are typically in the 32-bit integer range,
  118.         // i.e. all the upper bits are zero.
  119.         return (v >>> Integer.SIZE) == 0 ?
  120.             sum.add(v) :
  121.             sum.add(DD.of(v));
  122.     }

  124.     // Specialised statistic methods not directly supported by o.a.c.statistics.descriptive

  126.     /**
  127.      * Returns the arithmetic mean of the entries in the input arrays,
  128.      * or {@code NaN} if the combined length of the arrays is zero.
  129.      *
  130.      * <p>Supports a combined length above the maximum array size.
  131.      *
  132.      * <p>A two-pass, corrected algorithm is used, starting with the definitional formula
  133.      * computed using the array of stored values and then correcting this by adding the
  134.      * mean deviation of the data values from the arithmetic mean. See, e.g. "Comparison
  135.      * of Several Algorithms for Computing Sample Means and Variances," Robert F. Ling,
  136.      * Journal of the American Statistical Association, Vol. 69, No. 348 (Dec., 1974), pp.
  137.      * 859-866.
  138.      *
  139.      * @param samples Values.
  140.      * @return the mean of the values or NaN if length = 0
  141.      */
  142.     static double mean(Collection<double[]> samples) {
  143.         final double mean = samples.stream()
  144.             .map(Mean::of)
  145.             .reduce(Mean::combine)
  146.             .orElseGet(Mean::create)
  147.             .getAsDouble();
  148.         // Second-pass correction.
  149.         // Note: The correction may not be finite in the event of extreme values.
  150.         // In this case the calling method computation will fail when the mean
  151.         // is used and we do not check for overflow here.
  152.         final long n = samples.stream().mapToInt(x -> x.length).sum();
  153.         return mean + samples.stream()
  154.             .flatMapToDouble(Arrays::stream).map(v -> v - mean).sum() / n;
  155.     }

  157.     /**
  158.      * Returns the mean of the (signed) differences between corresponding elements of the
  159.      * input arrays.
  160.      *
  161.      * <pre>
  162.      * sum(x[i] - y[i]) / x.length
  163.      * </pre>
  164.      *
  165.      * <p>This method avoids intermediate array allocation.
  166.      *
  167.      * @param x First array.
  168.      * @param y Second array.
  169.      * @return mean of paired differences
  170.      * @throws IllegalArgumentException if the arrays do not have the same length.
  171.      */
  172.     static double meanDifference(double[] x, double[] y) {
  173.         final int n = x.length;
  174.         if (n != y.length) {
  175.             throw new InferenceException(InferenceException.VALUES_MISMATCH, n, y.length);
  176.         }
  177.         // STATISTICS-84: Use a single-pass extended precision sum.
  178.         final Sum sum = Sum.create();
  179.         for (int i = 0; i < n; i++) {
  180.             sum.add(x[i] - y[i]);
  181.         }
  182.         return sum.getAsDouble() / n;
  183.     }

  185.     /**
  186.      * Returns the variance of the (signed) differences between corresponding elements of
  187.      * the input arrays, or {@code NaN} if the arrays are empty.
  188.      *
  189.      * <pre>
  190.      * var(x[i] - y[i])
  191.      * </pre>
  192.      *
  193.      * <p>Returns the bias-corrected sample variance (using {@code n - 1} in the denominator).
  194.      * Returns 0 for a single-value (i.e. length = 1) sample.
  195.      *
  196.      * <p>This method avoids intermediate array allocation.
  197.      *
  198.      * <p>Uses a two-pass algorithm. Specifically, these methods use the "corrected
  199.      * two-pass algorithm" from Chan, Golub, Levesque, <i>Algorithms for Computing the
  200.      * Sample Variance</i>, American Statistician, vol. 37, no. 3 (1983) pp.
  201.      * 242-247.
  202.      *
  203.      * @param x First array.
  204.      * @param y Second array.
  205.      * @param mean the mean difference between corresponding entries
  206.      * @return variance of paired differences
  207.      * @throws IllegalArgumentException if the arrays do not have the same length.
  208.      * @see #meanDifference(double[], double[])
  209.      */
  210.     static double varianceDifference(double[] x, double[] y, double mean) {
  211.         final int n = x.length;
  212.         if (n != y.length) {
  213.             throw new InferenceException(InferenceException.VALUES_MISMATCH, n, y.length);
  214.         }
  215.         if (n == 1) {
  216.             return 0;
  217.         }
  218.         // As per o.a.c.statistics.descriptive.Variance
  219.         double s = 0;
  220.         double ss = 0;
  221.         for (int i = 0; i < n; i++) {
  222.             final double dx = (x[i] - y[i]) - mean;
  223.             s += dx;
  224.             ss += dx * dx;
  225.         }
  226.         // sum-of-squared deviations = sum(x^2) - sum(x)^2 / n
  227.         // Divide by (n-1) for sample variance
  228.         return (ss - (s * s / n)) / (n - 1);
  229.     }
  230. }
Created with JaCoCo 0.8.12.202403310830