StatisticUtils.java
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.statistics.inference;
import java.util.Arrays;
import java.util.Collection;
import org.apache.commons.numbers.core.DD;
import org.apache.commons.numbers.core.Precision;
import org.apache.commons.numbers.core.Sum;
import org.apache.commons.statistics.descriptive.Mean;
/**
* Utility computation methods.
*
* @since 1.1
*/
final class StatisticUtils {
/** No instances. */
private StatisticUtils() {}
/**
* Compute {@code x - y}.
*
* <p>If {@code y} is zero the original array is returned, else a new array is created
* with the difference.
*
* @param x Array.
* @param y Value.
* @return x - y
* @throws NullPointerException if {@code x} is null and {@code y} is non-zero
*/
static double[] subtract(double[] x, double y) {
return y == 0 ? x : Arrays.stream(x).map(v -> v - y).toArray();
}
/**
* Compute the degrees of freedom as {@code n - 1 - m}.
*
* <p>This method is common functionality shared between the Chi-square test and
* G-test. The pre-conditions for those tests are performed by this method.
*
* @param n Number of observations.
* @param m Adjustment (assumed to be positive).
* @return the degrees of freedom
* @throws IllegalArgumentException if the degrees of freedom is not strictly positive
*/
static int computeDegreesOfFreedom(int n, int m) {
final int df = n - 1 - m;
if (df <= 0) {
throw new InferenceException("Invalid degrees of freedom: " + df);
}
return df;
}
/**
* Gets the ratio between the sum of the observed and expected values.
* The ratio can be used to scale the expected values to have the same sum
* as the observed values:
*
* <pre>
* sum(o) = sum(e * ratio)
* </pre>
*
* <p>This method is common functionality shared between the Chi-square test and
* G-test. The pre-conditions for those tests are performed by this method.
*
* @param expected Expected values.
* @param observed Observed values.
* @return the ratio
* @throws IllegalArgumentException if the sample size is less than 2; the array
* sizes do not match; {@code expected} has entries that are not strictly
* positive; {@code observed} has negative entries; or all the observations are zero.
*/
static double computeRatio(double[] expected, long[] observed) {
Arguments.checkValuesRequiredSize(expected.length, 2);
Arguments.checkValuesSizeMatch(expected.length, observed.length);
Arguments.checkStrictlyPositive(expected);
Arguments.checkNonNegative(observed);
DD e = DD.ZERO;
DD o = DD.ZERO;
for (int i = 0; i < observed.length; i++) {
e = e.add(expected[i]);
o = add(o, observed[i]);
}
if (o.doubleValue() == 0) {
throw new InferenceException(InferenceException.NO_DATA);
}
// sum(o) / sum(e)
final double ratio = o.divide(e).doubleValue();
// Allow a sum within 1 ulp of 1.0
return Precision.equals(ratio, 1.0, 0) ? 1.0 : ratio;
}
/**
* Adds the value to the sum.
*
* @param sum Sum.
* @param v Value.
* @return the new sum
*/
private static DD add(DD sum, long v) {
// The condition here is a high probability branch if the sample is
// frequency counts which are typically in the 32-bit integer range,
// i.e. all the upper bits are zero.
return (v >>> Integer.SIZE) == 0 ?
sum.add(v) :
sum.add(DD.of(v));
}
// Specialised statistic methods not directly supported by o.a.c.statistics.descriptive
/**
* Returns the arithmetic mean of the entries in the input arrays,
* or {@code NaN} if the combined length of the arrays is zero.
*
* <p>Supports a combined length above the maximum array size.
*
* <p>A two-pass, corrected algorithm is used, starting with the definitional formula
* computed using the array of stored values and then correcting this by adding the
* mean deviation of the data values from the arithmetic mean. See, e.g. "Comparison
* of Several Algorithms for Computing Sample Means and Variances," Robert F. Ling,
* Journal of the American Statistical Association, Vol. 69, No. 348 (Dec., 1974), pp.
* 859-866.
*
* @param samples Values.
* @return the mean of the values or NaN if length = 0
*/
static double mean(Collection<double[]> samples) {
final double mean = samples.stream()
.map(Mean::of)
.reduce(Mean::combine)
.orElseGet(Mean::create)
.getAsDouble();
// Second-pass correction.
// Note: The correction may not be finite in the event of extreme values.
// In this case the calling method computation will fail when the mean
// is used and we do not check for overflow here.
final long n = samples.stream().mapToInt(x -> x.length).sum();
return mean + samples.stream()
.flatMapToDouble(Arrays::stream).map(v -> v - mean).sum() / n;
}
/**
* Returns the mean of the (signed) differences between corresponding elements of the
* input arrays.
*
* <pre>
* sum(x[i] - y[i]) / x.length
* </pre>
*
* <p>This method avoids intermediate array allocation.
*
* @param x First array.
* @param y Second array.
* @return mean of paired differences
* @throws IllegalArgumentException if the arrays do not have the same length.
*/
static double meanDifference(double[] x, double[] y) {
final int n = x.length;
if (n != y.length) {
throw new InferenceException(InferenceException.VALUES_MISMATCH, n, y.length);
}
// STATISTICS-84: Use a single-pass extended precision sum.
final Sum sum = Sum.create();
for (int i = 0; i < n; i++) {
sum.add(x[i] - y[i]);
}
return sum.getAsDouble() / n;
}
/**
* Returns the variance of the (signed) differences between corresponding elements of
* the input arrays, or {@code NaN} if the arrays are empty.
*
* <pre>
* var(x[i] - y[i])
* </pre>
*
* <p>Returns the bias-corrected sample variance (using {@code n - 1} in the denominator).
* Returns 0 for a single-value (i.e. length = 1) sample.
*
* <p>This method avoids intermediate array allocation.
*
* <p>Uses a two-pass algorithm. Specifically, these methods use the "corrected
* two-pass algorithm" from Chan, Golub, Levesque, <i>Algorithms for Computing the
* Sample Variance</i>, American Statistician, vol. 37, no. 3 (1983) pp.
* 242-247.
*
* @param x First array.
* @param y Second array.
* @param mean the mean difference between corresponding entries
* @return variance of paired differences
* @throws IllegalArgumentException if the arrays do not have the same length.
* @see #meanDifference(double[], double[])
*/
static double varianceDifference(double[] x, double[] y, double mean) {
final int n = x.length;
if (n != y.length) {
throw new InferenceException(InferenceException.VALUES_MISMATCH, n, y.length);
}
if (n == 1) {
return 0;
}
// As per o.a.c.statistics.descriptive.Variance
double s = 0;
double ss = 0;
for (int i = 0; i < n; i++) {
final double dx = (x[i] - y[i]) - mean;
s += dx;
ss += dx * dx;
}
// sum-of-squared deviations = sum(x^2) - sum(x)^2 / n
// Divide by (n-1) for sample variance
return (ss - (s * s / n)) / (n - 1);
}
}