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.checkArgument;
018import static com.google.common.base.Preconditions.checkNotNull;
019import static com.google.common.base.Preconditions.checkState;
020import static java.lang.Double.NaN;
021import static java.lang.Double.doubleToLongBits;
022import static java.lang.Double.isNaN;
023
024import com.google.common.annotations.Beta;
025import com.google.common.annotations.GwtIncompatible;
026import com.google.common.base.MoreObjects;
027import com.google.common.base.Objects;
028import java.io.Serializable;
029import java.nio.ByteBuffer;
030import java.nio.ByteOrder;
031import org.checkerframework.checker.nullness.qual.Nullable;
032
033/**
034 * An immutable value object capturing some basic statistics about a collection of paired double
035 * values (e.g. points on a plane). Build instances with {@link PairedStatsAccumulator#snapshot}.
036 *
037 * @author Pete Gillin
038 * @since 20.0
039 */
040@Beta
041@GwtIncompatible
042public final class PairedStats implements Serializable {
043
044  private final Stats xStats;
045  private final Stats yStats;
046  private final double sumOfProductsOfDeltas;
047
048  /**
049   * Internal constructor. Users should use {@link PairedStatsAccumulator#snapshot}.
050   *
051   * <p>To ensure that the created instance obeys its contract, the parameters should satisfy the
052   * following constraints. This is the callers responsibility and is not enforced here.
053   *
054   * <ul>
055   *   <li>Both {@code xStats} and {@code yStats} must have the same {@code count}.
056   *   <li>If that {@code count} is 1, {@code sumOfProductsOfDeltas} must be exactly 0.0.
057   *   <li>If that {@code count} is more than 1, {@code sumOfProductsOfDeltas} must be finite.
058   * </ul>
059   */
060  PairedStats(Stats xStats, Stats yStats, double sumOfProductsOfDeltas) {
061    this.xStats = xStats;
062    this.yStats = yStats;
063    this.sumOfProductsOfDeltas = sumOfProductsOfDeltas;
064  }
065
066  /** Returns the number of pairs in the dataset. */
067  public long count() {
068    return xStats.count();
069  }
070
071  /** Returns the statistics on the {@code x} values alone. */
072  public Stats xStats() {
073    return xStats;
074  }
075
076  /** Returns the statistics on the {@code y} values alone. */
077  public Stats yStats() {
078    return yStats;
079  }
080
081  /**
082   * Returns the population covariance of the values. The count must be non-zero.
083   *
084   * <p>This is guaranteed to return zero if the dataset contains a single pair of finite values. It
085   * is not guaranteed to return zero when the dataset consists of the same pair of values multiple
086   * times, due to numerical errors.
087   *
088   * <h3>Non-finite values</h3>
089   *
090   * <p>If the dataset contains any non-finite values ({@link Double#POSITIVE_INFINITY}, {@link
091   * Double#NEGATIVE_INFINITY}, or {@link Double#NaN}) then the result is {@link Double#NaN}.
092   *
093   * @throws IllegalStateException if the dataset is empty
094   */
095  public double populationCovariance() {
096    checkState(count() != 0);
097    return sumOfProductsOfDeltas / count();
098  }
099
100  /**
101   * Returns the sample covariance of the values. The count must be greater than one.
102   *
103   * <p>This is not guaranteed to return zero when the dataset consists of the same pair of values
104   * multiple times, due to numerical errors.
105   *
106   * <h3>Non-finite values</h3>
107   *
108   * <p>If the dataset contains any non-finite values ({@link Double#POSITIVE_INFINITY}, {@link
109   * Double#NEGATIVE_INFINITY}, or {@link Double#NaN}) then the result is {@link Double#NaN}.
110   *
111   * @throws IllegalStateException if the dataset is empty or contains a single pair of values
112   */
113  public double sampleCovariance() {
114    checkState(count() > 1);
115    return sumOfProductsOfDeltas / (count() - 1);
116  }
117
118  /**
119   * Returns the <a href="http://mathworld.wolfram.com/CorrelationCoefficient.html">Pearson's or
120   * product-moment correlation coefficient</a> of the values. The count must greater than one, and
121   * the {@code x} and {@code y} values must both have non-zero population variance (i.e. {@code
122   * xStats().populationVariance() > 0.0 && yStats().populationVariance() > 0.0}). The result is not
123   * guaranteed to be exactly +/-1 even when the data are perfectly (anti-)correlated, due to
124   * numerical errors. However, it is guaranteed to be in the inclusive range [-1, +1].
125   *
126   * <h3>Non-finite values</h3>
127   *
128   * <p>If the dataset contains any non-finite values ({@link Double#POSITIVE_INFINITY}, {@link
129   * Double#NEGATIVE_INFINITY}, or {@link Double#NaN}) then the result is {@link Double#NaN}.
130   *
131   * @throws IllegalStateException if the dataset is empty or contains a single pair of values, or
132   *     either the {@code x} and {@code y} dataset has zero population variance
133   */
134  public double pearsonsCorrelationCoefficient() {
135    checkState(count() > 1);
136    if (isNaN(sumOfProductsOfDeltas)) {
137      return NaN;
138    }
139    double xSumOfSquaresOfDeltas = xStats().sumOfSquaresOfDeltas();
140    double ySumOfSquaresOfDeltas = yStats().sumOfSquaresOfDeltas();
141    checkState(xSumOfSquaresOfDeltas > 0.0);
142    checkState(ySumOfSquaresOfDeltas > 0.0);
143    // The product of two positive numbers can be zero if the multiplication underflowed. We
144    // force a positive value by effectively rounding up to MIN_VALUE.
145    double productOfSumsOfSquaresOfDeltas =
146        ensurePositive(xSumOfSquaresOfDeltas * ySumOfSquaresOfDeltas);
147    return ensureInUnitRange(sumOfProductsOfDeltas / Math.sqrt(productOfSumsOfSquaresOfDeltas));
148  }
149
150  /**
151   * Returns a linear transformation giving the best fit to the data according to <a
152   * href="http://mathworld.wolfram.com/LeastSquaresFitting.html">Ordinary Least Squares linear
153   * regression</a> of {@code y} as a function of {@code x}. The count must be greater than one, and
154   * either the {@code x} or {@code y} data must have a non-zero population variance (i.e. {@code
155   * xStats().populationVariance() > 0.0 || yStats().populationVariance() > 0.0}). The result is
156   * guaranteed to be horizontal if there is variance in the {@code x} data but not the {@code y}
157   * data, and vertical if there is variance in the {@code y} data but not the {@code x} data.
158   *
159   * <p>This fit minimizes the root-mean-square error in {@code y} as a function of {@code x}. This
160   * error is defined as the square root of the mean of the squares of the differences between the
161   * actual {@code y} values of the data and the values predicted by the fit for the {@code x}
162   * values (i.e. it is the square root of the mean of the squares of the vertical distances between
163   * the data points and the best fit line). For this fit, this error is a fraction {@code sqrt(1 -
164   * R*R)} of the population standard deviation of {@code y}, where {@code R} is the Pearson's
165   * correlation coefficient (as given by {@link #pearsonsCorrelationCoefficient()}).
166   *
167   * <p>The corresponding root-mean-square error in {@code x} as a function of {@code y} is a
168   * fraction {@code sqrt(1/(R*R) - 1)} of the population standard deviation of {@code x}. This fit
169   * does not normally minimize that error: to do that, you should swap the roles of {@code x} and
170   * {@code y}.
171   *
172   * <h3>Non-finite values</h3>
173   *
174   * <p>If the dataset contains any non-finite values ({@link Double#POSITIVE_INFINITY}, {@link
175   * Double#NEGATIVE_INFINITY}, or {@link Double#NaN}) then the result is {@link
176   * LinearTransformation#forNaN()}.
177   *
178   * @throws IllegalStateException if the dataset is empty or contains a single pair of values, or
179   *     both the {@code x} and {@code y} dataset must have zero population variance
180   */
181  public LinearTransformation leastSquaresFit() {
182    checkState(count() > 1);
183    if (isNaN(sumOfProductsOfDeltas)) {
184      return LinearTransformation.forNaN();
185    }
186    double xSumOfSquaresOfDeltas = xStats.sumOfSquaresOfDeltas();
187    if (xSumOfSquaresOfDeltas > 0.0) {
188      if (yStats.sumOfSquaresOfDeltas() > 0.0) {
189        return LinearTransformation.mapping(xStats.mean(), yStats.mean())
190            .withSlope(sumOfProductsOfDeltas / xSumOfSquaresOfDeltas);
191      } else {
192        return LinearTransformation.horizontal(yStats.mean());
193      }
194    } else {
195      checkState(yStats.sumOfSquaresOfDeltas() > 0.0);
196      return LinearTransformation.vertical(xStats.mean());
197    }
198  }
199
200  /**
201   * {@inheritDoc}
202   *
203   * <p><b>Note:</b> This tests exact equality of the calculated statistics, including the floating
204   * point values. Two instances are guaranteed to be considered equal if one is copied from the
205   * other using {@code second = new PairedStatsAccumulator().addAll(first).snapshot()}, if both
206   * were obtained by calling {@code snapshot()} on the same {@link PairedStatsAccumulator} without
207   * adding any values in between the two calls, or if one is obtained from the other after
208   * round-tripping through java serialization. However, floating point rounding errors mean that it
209   * may be false for some instances where the statistics are mathematically equal, including
210   * instances constructed from the same values in a different order... or (in the general case)
211   * even in the same order. (It is guaranteed to return true for instances constructed from the
212   * same values in the same order if {@code strictfp} is in effect, or if the system architecture
213   * guarantees {@code strictfp}-like semantics.)
214   */
215  @Override
216  public boolean equals(@Nullable Object obj) {
217    if (obj == null) {
218      return false;
219    }
220    if (getClass() != obj.getClass()) {
221      return false;
222    }
223    PairedStats other = (PairedStats) obj;
224    return (xStats.equals(other.xStats))
225        && (yStats.equals(other.yStats))
226        && (doubleToLongBits(sumOfProductsOfDeltas)
227            == doubleToLongBits(other.sumOfProductsOfDeltas));
228  }
229
230  /**
231   * {@inheritDoc}
232   *
233   * <p><b>Note:</b> This hash code is consistent with exact equality of the calculated statistics,
234   * including the floating point values. See the note on {@link #equals} for details.
235   */
236  @Override
237  public int hashCode() {
238    return Objects.hashCode(xStats, yStats, sumOfProductsOfDeltas);
239  }
240
241  @Override
242  public String toString() {
243    if (count() > 0) {
244      return MoreObjects.toStringHelper(this)
245          .add("xStats", xStats)
246          .add("yStats", yStats)
247          .add("populationCovariance", populationCovariance())
248          .toString();
249    } else {
250      return MoreObjects.toStringHelper(this)
251          .add("xStats", xStats)
252          .add("yStats", yStats)
253          .toString();
254    }
255  }
256
257  double sumOfProductsOfDeltas() {
258    return sumOfProductsOfDeltas;
259  }
260
261  private static double ensurePositive(double value) {
262    if (value > 0.0) {
263      return value;
264    } else {
265      return Double.MIN_VALUE;
266    }
267  }
268
269  private static double ensureInUnitRange(double value) {
270    if (value >= 1.0) {
271      return 1.0;
272    }
273    if (value <= -1.0) {
274      return -1.0;
275    }
276    return value;
277  }
278
279  // Serialization helpers
280
281  /** The size of byte array representation in bytes. */
282  private static final int BYTES = Stats.BYTES * 2 + Double.SIZE / Byte.SIZE;
283
284  /**
285   * Gets a byte array representation of this instance.
286   *
287   * <p><b>Note:</b> No guarantees are made regarding stability of the representation between
288   * versions.
289   */
290  public byte[] toByteArray() {
291    ByteBuffer buffer = ByteBuffer.allocate(BYTES).order(ByteOrder.LITTLE_ENDIAN);
292    xStats.writeTo(buffer);
293    yStats.writeTo(buffer);
294    buffer.putDouble(sumOfProductsOfDeltas);
295    return buffer.array();
296  }
297
298  /**
299   * Creates a {@link PairedStats} instance from the given byte representation which was obtained by
300   * {@link #toByteArray}.
301   *
302   * <p><b>Note:</b> No guarantees are made regarding stability of the representation between
303   * versions.
304   */
305  public static PairedStats fromByteArray(byte[] byteArray) {
306    checkNotNull(byteArray);
307    checkArgument(
308        byteArray.length == BYTES,
309        "Expected PairedStats.BYTES = %s, got %s",
310        BYTES,
311        byteArray.length);
312    ByteBuffer buffer = ByteBuffer.wrap(byteArray).order(ByteOrder.LITTLE_ENDIAN);
313    Stats xStats = Stats.readFrom(buffer);
314    Stats yStats = Stats.readFrom(buffer);
315    double sumOfProductsOfDeltas = buffer.getDouble();
316    return new PairedStats(xStats, yStats, sumOfProductsOfDeltas);
317  }
318
319  private static final long serialVersionUID = 0;
320}