001/* 002 * Copyright (C) 2012 The Guava Authors 003 * 004 * Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except 005 * in compliance with the License. You may obtain a copy of the License at 006 * 007 * http://www.apache.org/licenses/LICENSE-2.0 008 * 009 * Unless required by applicable law or agreed to in writing, software distributed under the License 010 * is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express 011 * or implied. See the License for the specific language governing permissions and limitations under 012 * the License. 013 */ 014 015package com.google.common.math; 016 017import static com.google.common.base.Preconditions.checkState; 018import static com.google.common.primitives.Doubles.isFinite; 019import static java.lang.Double.NaN; 020import static java.lang.Double.isNaN; 021 022import com.google.common.annotations.Beta; 023import com.google.common.annotations.GwtIncompatible; 024 025/** 026 * A mutable object which accumulates paired double values (e.g. points on a plane) and tracks some 027 * basic statistics over all the values added so far. This class is not thread safe. 028 * 029 * @author Pete Gillin 030 * @since 20.0 031 */ 032@Beta 033@GwtIncompatible 034public final class PairedStatsAccumulator { 035 036 // These fields must satisfy the requirements of PairedStats' constructor as well as those of the 037 // stat methods of this class. 038 private final StatsAccumulator xStats = new StatsAccumulator(); 039 private final StatsAccumulator yStats = new StatsAccumulator(); 040 private double sumOfProductsOfDeltas = 0.0; 041 042 /** Adds the given pair of values to the dataset. */ 043 public void add(double x, double y) { 044 // We extend the recursive expression for the one-variable case at Art of Computer Programming 045 // vol. 2, Knuth, 4.2.2, (16) to the two-variable case. We have two value series x_i and y_i. 046 // We define the arithmetic means X_n = 1/n \sum_{i=1}^n x_i, and Y_n = 1/n \sum_{i=1}^n y_i. 047 // We also define the sum of the products of the differences from the means 048 // C_n = \sum_{i=1}^n x_i y_i - n X_n Y_n 049 // for all n >= 1. Then for all n > 1: 050 // C_{n-1} = \sum_{i=1}^{n-1} x_i y_i - (n-1) X_{n-1} Y_{n-1} 051 // C_n - C_{n-1} = x_n y_n - n X_n Y_n + (n-1) X_{n-1} Y_{n-1} 052 // = x_n y_n - X_n [ y_n + (n-1) Y_{n-1} ] + [ n X_n - x_n ] Y_{n-1} 053 // = x_n y_n - X_n y_n - x_n Y_{n-1} + X_n Y_{n-1} 054 // = (x_n - X_n) (y_n - Y_{n-1}) 055 xStats.add(x); 056 if (isFinite(x) && isFinite(y)) { 057 if (xStats.count() > 1) { 058 sumOfProductsOfDeltas += (x - xStats.mean()) * (y - yStats.mean()); 059 } 060 } else { 061 sumOfProductsOfDeltas = NaN; 062 } 063 yStats.add(y); 064 } 065 066 /** 067 * Adds the given statistics to the dataset, as if the individual values used to compute the 068 * statistics had been added directly. 069 */ 070 public void addAll(PairedStats values) { 071 if (values.count() == 0) { 072 return; 073 } 074 075 xStats.addAll(values.xStats()); 076 if (yStats.count() == 0) { 077 sumOfProductsOfDeltas = values.sumOfProductsOfDeltas(); 078 } else { 079 // This is a generalized version of the calculation in add(double, double) above. Note that 080 // non-finite inputs will have sumOfProductsOfDeltas = NaN, so non-finite values will result 081 // in NaN naturally. 082 sumOfProductsOfDeltas += 083 values.sumOfProductsOfDeltas() 084 + (values.xStats().mean() - xStats.mean()) 085 * (values.yStats().mean() - yStats.mean()) 086 * values.count(); 087 } 088 yStats.addAll(values.yStats()); 089 } 090 091 /** Returns an immutable snapshot of the current statistics. */ 092 public PairedStats snapshot() { 093 return new PairedStats(xStats.snapshot(), yStats.snapshot(), sumOfProductsOfDeltas); 094 } 095 096 /** Returns the number of pairs in the dataset. */ 097 public long count() { 098 return xStats.count(); 099 } 100 101 /** Returns an immutable snapshot of the statistics on the {@code x} values alone. */ 102 public Stats xStats() { 103 return xStats.snapshot(); 104 } 105 106 /** Returns an immutable snapshot of the statistics on the {@code y} values alone. */ 107 public Stats yStats() { 108 return yStats.snapshot(); 109 } 110 111 /** 112 * Returns the population covariance of the values. The count must be non-zero. 113 * 114 * <p>This is guaranteed to return zero if the dataset contains a single pair of finite values. It 115 * is not guaranteed to return zero when the dataset consists of the same pair of values multiple 116 * times, due to numerical errors. 117 * 118 * <h3>Non-finite values</h3> 119 * 120 * <p>If the dataset contains any non-finite values ({@link Double#POSITIVE_INFINITY}, {@link 121 * Double#NEGATIVE_INFINITY}, or {@link Double#NaN}) then the result is {@link Double#NaN}. 122 * 123 * @throws IllegalStateException if the dataset is empty 124 */ 125 public double populationCovariance() { 126 checkState(count() != 0); 127 return sumOfProductsOfDeltas / count(); 128 } 129 130 /** 131 * Returns the sample covariance of the values. The count must be greater than one. 132 * 133 * <p>This is not guaranteed to return zero when the dataset consists of the same pair of values 134 * multiple times, due to numerical errors. 135 * 136 * <h3>Non-finite values</h3> 137 * 138 * <p>If the dataset contains any non-finite values ({@link Double#POSITIVE_INFINITY}, {@link 139 * Double#NEGATIVE_INFINITY}, or {@link Double#NaN}) then the result is {@link Double#NaN}. 140 * 141 * @throws IllegalStateException if the dataset is empty or contains a single pair of values 142 */ 143 public final double sampleCovariance() { 144 checkState(count() > 1); 145 return sumOfProductsOfDeltas / (count() - 1); 146 } 147 148 /** 149 * Returns the <a href="http://mathworld.wolfram.com/CorrelationCoefficient.html">Pearson's or 150 * product-moment correlation coefficient</a> of the values. The count must greater than one, and 151 * the {@code x} and {@code y} values must both have non-zero population variance (i.e. {@code 152 * xStats().populationVariance() > 0.0 && yStats().populationVariance() > 0.0}). The result is not 153 * guaranteed to be exactly +/-1 even when the data are perfectly (anti-)correlated, due to 154 * numerical errors. However, it is guaranteed to be in the inclusive range [-1, +1]. 155 * 156 * <h3>Non-finite values</h3> 157 * 158 * <p>If the dataset contains any non-finite values ({@link Double#POSITIVE_INFINITY}, {@link 159 * Double#NEGATIVE_INFINITY}, or {@link Double#NaN}) then the result is {@link Double#NaN}. 160 * 161 * @throws IllegalStateException if the dataset is empty or contains a single pair of values, or 162 * either the {@code x} and {@code y} dataset has zero population variance 163 */ 164 public final double pearsonsCorrelationCoefficient() { 165 checkState(count() > 1); 166 if (isNaN(sumOfProductsOfDeltas)) { 167 return NaN; 168 } 169 double xSumOfSquaresOfDeltas = xStats.sumOfSquaresOfDeltas(); 170 double ySumOfSquaresOfDeltas = yStats.sumOfSquaresOfDeltas(); 171 checkState(xSumOfSquaresOfDeltas > 0.0); 172 checkState(ySumOfSquaresOfDeltas > 0.0); 173 // The product of two positive numbers can be zero if the multiplication underflowed. We 174 // force a positive value by effectively rounding up to MIN_VALUE. 175 double productOfSumsOfSquaresOfDeltas = 176 ensurePositive(xSumOfSquaresOfDeltas * ySumOfSquaresOfDeltas); 177 return ensureInUnitRange(sumOfProductsOfDeltas / Math.sqrt(productOfSumsOfSquaresOfDeltas)); 178 } 179 180 /** 181 * Returns a linear transformation giving the best fit to the data according to <a 182 * href="http://mathworld.wolfram.com/LeastSquaresFitting.html">Ordinary Least Squares linear 183 * regression</a> of {@code y} as a function of {@code x}. The count must be greater than one, and 184 * either the {@code x} or {@code y} data must have a non-zero population variance (i.e. {@code 185 * xStats().populationVariance() > 0.0 || yStats().populationVariance() > 0.0}). The result is 186 * guaranteed to be horizontal if there is variance in the {@code x} data but not the {@code y} 187 * data, and vertical if there is variance in the {@code y} data but not the {@code x} data. 188 * 189 * <p>This fit minimizes the root-mean-square error in {@code y} as a function of {@code x}. This 190 * error is defined as the square root of the mean of the squares of the differences between the 191 * actual {@code y} values of the data and the values predicted by the fit for the {@code x} 192 * values (i.e. it is the square root of the mean of the squares of the vertical distances between 193 * the data points and the best fit line). For this fit, this error is a fraction {@code sqrt(1 - 194 * R*R)} of the population standard deviation of {@code y}, where {@code R} is the Pearson's 195 * correlation coefficient (as given by {@link #pearsonsCorrelationCoefficient()}). 196 * 197 * <p>The corresponding root-mean-square error in {@code x} as a function of {@code y} is a 198 * fraction {@code sqrt(1/(R*R) - 1)} of the population standard deviation of {@code x}. This fit 199 * does not normally minimize that error: to do that, you should swap the roles of {@code x} and 200 * {@code y}. 201 * 202 * <h3>Non-finite values</h3> 203 * 204 * <p>If the dataset contains any non-finite values ({@link Double#POSITIVE_INFINITY}, {@link 205 * Double#NEGATIVE_INFINITY}, or {@link Double#NaN}) then the result is {@link 206 * LinearTransformation#forNaN()}. 207 * 208 * @throws IllegalStateException if the dataset is empty or contains a single pair of values, or 209 * both the {@code x} and {@code y} dataset have zero population variance 210 */ 211 public final LinearTransformation leastSquaresFit() { 212 checkState(count() > 1); 213 if (isNaN(sumOfProductsOfDeltas)) { 214 return LinearTransformation.forNaN(); 215 } 216 double xSumOfSquaresOfDeltas = xStats.sumOfSquaresOfDeltas(); 217 if (xSumOfSquaresOfDeltas > 0.0) { 218 if (yStats.sumOfSquaresOfDeltas() > 0.0) { 219 return LinearTransformation.mapping(xStats.mean(), yStats.mean()) 220 .withSlope(sumOfProductsOfDeltas / xSumOfSquaresOfDeltas); 221 } else { 222 return LinearTransformation.horizontal(yStats.mean()); 223 } 224 } else { 225 checkState(yStats.sumOfSquaresOfDeltas() > 0.0); 226 return LinearTransformation.vertical(xStats.mean()); 227 } 228 } 229 230 private double ensurePositive(double value) { 231 if (value > 0.0) { 232 return value; 233 } else { 234 return Double.MIN_VALUE; 235 } 236 } 237 238 private static double ensureInUnitRange(double value) { 239 if (value >= 1.0) { 240 return 1.0; 241 } 242 if (value <= -1.0) { 243 return -1.0; 244 } 245 return value; 246 } 247}