1 files changed, 176 insertions, 0 deletions
diff --git a/src/de/lmu/ifi/dbs/elki/math/statistics/ProbabilityWeightedMoments.java b/src/de/lmu/ifi/dbs/elki/math/statistics/ProbabilityWeightedMoments.java
new file mode 100644
index 00000000..c783f5f1
--- /dev/null
+++ b/src/de/lmu/ifi/dbs/elki/math/statistics/ProbabilityWeightedMoments.java
@@ -0,0 +1,176 @@
+package de.lmu.ifi.dbs.elki.math.statistics;
+
+/*
+ This file is part of ELKI:
+ Environment for Developing KDD-Applications Supported by Index-Structures
+
+ Copyright (C) 2013
+ Ludwig-Maximilians-Universität München
+ Lehr- und Forschungseinheit für Datenbanksysteme
+ ELKI Development Team
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Affero General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU Affero General Public License for more details.
+
+ You should have received a copy of the GNU Affero General Public License
+ along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+import de.lmu.ifi.dbs.elki.utilities.datastructures.arraylike.NumberArrayAdapter;
+import de.lmu.ifi.dbs.elki.utilities.documentation.Reference;
+
+/**
+ * Estimate the L-Moments of a sample.
+ * 
+ * Reference:
+ * <p>
+ * J. R. M. Hosking, J. R. Wallis, and E. F. Wood<br />
+ * Estimation of the generalized extreme-value distribution by the method of
+ * probability-weighted moments.<br />
+ * Technometrics 27.3
+ * </p>
+ * 
+ * Also based on:
+ * <p>
+ * J. R. M. Hosking<br />
+ * Fortran routines for use with the method of L-moments Version 3.03<br />
+ * IBM Research.
+ * </p>
+ * 
+ * @author Erich Schubert
+ */
+@Reference(authors = "J.R.M. Hosking, J. R. Wallis, and E. F. Wood", title = "Estimation of the generalized extreme-value distribution by the method of probability-weighted moments.", booktitle = "Technometrics 27.3", url = "http://dx.doi.org/10.1080/00401706.1985.10488049")
+public class ProbabilityWeightedMoments {
+  /**
+   * Compute the alpha_r factors using the method of probability-weighted
+   * moments.
+   * 
+   * @param sorted <b>Presorted</b> data array.
+   * @param nmom Number of moments to compute
+   * @return Alpha moments (0-indexed)
+   */
+  public static <A> double[] alphaPWM(A data, NumberArrayAdapter<?, A> adapter, final int nmom) {
+    final int n = adapter.size(data);
+    final double[] xmom = new double[nmom];
+    double weight = 1. / n;
+    for (int i = 0; i < n; i++) {
+      final double val = adapter.getDouble(data, i);
+      xmom[0] += weight * val;
+      for (int j = 1; j < nmom; j++) {
+        weight *= (n - i - j + 1) / (n - j + 1);
+        xmom[j] += weight * val;
+      }
+    }
+    return xmom;
+  }
+
+  /**
+   * Compute the beta_r factors using the method of probability-weighted
+   * moments.
+   * 
+   * @param sorted <b>Presorted</b> data array.
+   * @param nmom Number of moments to compute
+   * @return Beta moments (0-indexed)
+   */
+  public static <A> double[] betaPWM(A data, NumberArrayAdapter<?, A> adapter, final int nmom) {
+    final int n = adapter.size(data);
+    final double[] xmom = new double[nmom];
+    double weight = 1. / n;
+    for (int i = 0; i < n; i++) {
+      final double val = adapter.getDouble(data, i);
+      xmom[0] += weight * val;
+      for (int j = 1; j < nmom; j++) {
+        weight *= (i - j + 1) / (n - j + 1);
+        xmom[j] += weight * val;
+      }
+    }
+    return xmom;
+  }
+
+  /**
+   * Compute the alpha_r and beta_r factors in parallel using the method of
+   * probability-weighted moments. Usually cheaper than computing them
+   * separately.
+   * 
+   * @param sorted <b>Presorted</b> data array.
+   * @param nmom Number of moments to compute
+   * @return Alpha and Beta moments (0-indexed, interleaved)
+   */
+  public static <A> double[] alphaBetaPWM(A data, NumberArrayAdapter<?, A> adapter, final int nmom) {
+    final int n = adapter.size(data);
+    final double[] xmom = new double[nmom << 1];
+    double aweight = 1. / n, bweight = aweight;
+    for (int i = 0; i < n; i++) {
+      final double val = adapter.getDouble(data, i);
+      xmom[0] += aweight * val;
+      xmom[1] += bweight * val;
+      for (int j = 1, k = 2; j < nmom; j++, k += 2) {
+        aweight *= (n - i - j + 1) / (n - j + 1);
+        bweight *= (i - j + 1) / (n - j + 1);
+        xmom[k + 1] += aweight * val;
+        xmom[k + 1] += bweight * val;
+      }
+    }
+    return xmom;
+  }
+
+  /**
+   * Compute the sample L-Moments using probability weighted moments.
+   * 
+   * @param sorted <b>Presorted</b> data array.
+   * @param nmom Number of moments to compute
+   * @return Array containing Lambda1, Lambda2, Tau3 ... TauN
+   */
+  public static <A> double[] samLMR(A sorted, NumberArrayAdapter<?, A> adapter, final int nmom) {
+    final int n = adapter.size(sorted);
+    if (nmom >= n) {
+      throw new ArithmeticException("Can't compute higher order moments for just" + n + " observations.");
+    }
+    final double[] sum = new double[nmom];
+    // Estimate probability weighted moments (unbiased)
+    for (int i = 0; i < n; i++) {
+      double term = adapter.getDouble(sorted, i);
+      // Robustness: skip bad values
+      if (Double.isInfinite(term) || Double.isNaN(term)) {
+        continue;
+      }
+      sum[0] += term;
+      for (int j = 1, z = i; j < nmom; j++, z--) {
+        term *= z;
+        sum[j] += term;
+      }
+    }
+    // Normalize by "n choose (j + 1)"
+    sum[0] /= n;
+    double z = n;
+    for (int j = 1; j < nmom; j++) {
+      z *= n - j;
+      sum[j] /= z;
+    }
+    for (int k = nmom - 1; k >= 1; --k) {
+      double p = ((k & 1) == 0) ? +1 : -1;
+      double temp = p * sum[0];
+      for (int i = 0; i < k; i++) {
+        double ai = i + 1.;
+        p *= -(k + ai) * (k - i) / (ai * ai);
+        temp += p * sum[i + 1];
+      }
+      sum[k] = temp;
+    }
+    if (nmom > 2 && !(sum[1] > 0)) {
+      throw new ArithmeticException("Can't compute higher order moments for constant data. Sum: " + sum[1]);
+    }
+    // Map lambda3...lambdaN to tau3...tauN
+    for (int i = 2; i < nmom; i++) {
+      sum[i] /= sum[1];
+    }
+    return sum;
+  }
+}