mahout-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From gsing...@apache.org
Subject svn commit: r891983 [12/47] - in /lucene/mahout/trunk: ./ core/ core/src/main/java/org/apache/mahout/cf/taste/hadoop/item/ core/src/main/java/org/apache/mahout/clustering/ core/src/main/java/org/apache/mahout/clustering/canopy/ core/src/main/java/org/a...
Date Thu, 17 Dec 2009 23:22:41 GMT
Added: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/BreitWigner.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/BreitWigner.java?rev=891983&view=auto
==============================================================================
--- lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/BreitWigner.java (added)
+++ lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/BreitWigner.java Thu Dec 17 23:22:16 2009
@@ -0,0 +1,92 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
+is hereby granted without fee, provided that the above copyright notice appear in all copies and 
+that both that copyright notice and this permission notice appear in supporting documentation. 
+CERN makes no representations about the suitability of this software for any purpose. 
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.random.engine.RandomEngine;
+
+/** @deprecated until unit tests are in place.  Until this time, this class/interface is unsupported. */
+@Deprecated
+public class BreitWigner extends AbstractContinousDistribution {
+
+  private double mean;
+  private double gamma;
+  private double cut;
+
+  // The uniform random number generated shared by all <b>static</b> methods.
+  private static final BreitWigner shared = new BreitWigner(1.0, 0.2, 1.0, makeDefaultGenerator());
+
+  /**
+   * Constructs a BreitWigner distribution.
+   *
+   * @param cut </tt>cut==Double.NEGATIVE_INFINITY</tt> indicates "don't cut".
+   */
+  public BreitWigner(double mean, double gamma, double cut, RandomEngine randomGenerator) {
+    setRandomGenerator(randomGenerator);
+    setState(mean, gamma, cut);
+  }
+
+  /** Returns a random number from the distribution. */
+  @Override
+  public double nextDouble() {
+    return nextDouble(mean, gamma, cut);
+  }
+
+  /**
+   * Returns a random number from the distribution; bypasses the internal state.
+   *
+   * @param cut </tt>cut==Double.NEGATIVE_INFINITY</tt> indicates "don't cut".
+   */
+  public double nextDouble(double mean, double gamma, double cut) {
+
+    if (gamma == 0.0) {
+      return mean;
+    }
+    double displ;
+    double rval;
+    if (cut == Double.NEGATIVE_INFINITY) { // don't cut
+      rval = 2.0 * randomGenerator.raw() - 1.0;
+      displ = 0.5 * gamma * Math.tan(rval * (Math.PI / 2.0));
+      return mean + displ;
+    } else {
+      double val = Math.atan(2.0 * cut / gamma);
+      rval = 2.0 * randomGenerator.raw() - 1.0;
+      displ = 0.5 * gamma * Math.tan(rval * val);
+
+      return mean + displ;
+    }
+  }
+
+  /**
+   * Sets the mean, gamma and cut parameters.
+   *
+   * @param cut </tt>cut==Double.NEGATIVE_INFINITY</tt> indicates "don't cut".
+   */
+  public void setState(double mean, double gamma, double cut) {
+    this.mean = mean;
+    this.gamma = gamma;
+    this.cut = cut;
+  }
+
+  /**
+   * Returns a random number from the distribution.
+   *
+   * @param cut </tt>cut==Double.NEGATIVE_INFINITY</tt> indicates "don't cut".
+   */
+  public static double staticNextDouble(double mean, double gamma, double cut) {
+    synchronized (shared) {
+      return shared.nextDouble(mean, gamma, cut);
+    }
+  }
+
+  /** Returns a String representation of the receiver. */
+  public String toString() {
+    return this.getClass().getName() + '(' + mean + ',' + gamma + ',' + cut + ')';
+  }
+
+}

Propchange: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/BreitWigner.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/BreitWignerMeanSquare.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/BreitWignerMeanSquare.java?rev=891983&view=auto
==============================================================================
--- lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/BreitWignerMeanSquare.java (added)
+++ lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/BreitWignerMeanSquare.java Thu Dec 17 23:22:16 2009
@@ -0,0 +1,84 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
+is hereby granted without fee, provided that the above copyright notice appear in all copies and 
+that both that copyright notice and this permission notice appear in supporting documentation. 
+CERN makes no representations about the suitability of this software for any purpose. 
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.random.engine.RandomEngine;
+
+/** @deprecated until unit tests are in place.  Until this time, this class/interface is unsupported. */
+@Deprecated
+public class BreitWignerMeanSquare extends BreitWigner {
+
+  private Uniform uniform; // helper
+
+  // The uniform random number generated shared by all <b>static</b> methods.
+  private static final BreitWigner shared = new BreitWignerMeanSquare(1.0, 0.2, 1.0, makeDefaultGenerator());
+
+  /**
+   * Constructs a mean-squared BreitWigner distribution.
+   *
+   * @param cut </tt>cut==Double.NEGATIVE_INFINITY</tt> indicates "don't cut".
+   */
+  public BreitWignerMeanSquare(double mean, double gamma, double cut, RandomEngine randomGenerator) {
+    super(mean, gamma, cut, randomGenerator);
+    this.uniform = new Uniform(randomGenerator);
+  }
+
+  /**
+   * Returns a deep copy of the receiver; the copy will produce identical sequences. After this call has returned, the
+   * copy and the receiver have equal but separate state.
+   *
+   * @return a copy of the receiver.
+   */
+  @Override
+  public Object clone() {
+    BreitWignerMeanSquare copy = (BreitWignerMeanSquare) super.clone();
+    if (this.uniform != null) {
+      copy.uniform = new Uniform(copy.randomGenerator);
+    }
+    return copy;
+  }
+
+  /**
+   * Returns a mean-squared random number from the distribution; bypasses the internal state.
+   *
+   * @param cut </tt>cut==Double.NEGATIVE_INFINITY</tt> indicates "don't cut".
+   */
+  @Override
+  public double nextDouble(double mean, double gamma, double cut) {
+    if (gamma == 0.0) {
+      return mean;
+    }
+    if (cut == Double.NEGATIVE_INFINITY) { // don't cut
+      double val = Math.atan(-mean / gamma);
+      double rval = this.uniform.nextDoubleFromTo(val, Math.PI / 2.0);
+      double displ = gamma * Math.tan(rval);
+      return Math.sqrt(mean * mean + mean * displ);
+    } else {
+      double tmp = Math.max(0.0, mean - cut);
+      double lower = Math.atan((tmp * tmp - mean * mean) / (mean * gamma));
+      double upper = Math.atan(((mean + cut) * (mean + cut) - mean * mean) / (mean * gamma));
+      double rval = this.uniform.nextDoubleFromTo(lower, upper);
+
+      double displ = gamma * Math.tan(rval);
+      return Math.sqrt(Math.max(0.0, mean * mean + mean * displ));
+    }
+  }
+
+  /**
+   * Returns a random number from the distribution.
+   *
+   * @param cut </tt>cut==Double.NEGATIVE_INFINITY</tt> indicates "don't cut".
+   */
+  public static double staticNextDouble(double mean, double gamma, double cut) {
+    synchronized (shared) {
+      return shared.nextDouble(mean, gamma, cut);
+    }
+  }
+
+}

Propchange: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/BreitWignerMeanSquare.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/ChiSquare.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/ChiSquare.java?rev=891983&view=auto
==============================================================================
--- lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/ChiSquare.java (added)
+++ lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/ChiSquare.java Thu Dec 17 23:22:16 2009
@@ -0,0 +1,174 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
+is hereby granted without fee, provided that the above copyright notice appear in all copies and 
+that both that copyright notice and this permission notice appear in supporting documentation. 
+CERN makes no representations about the suitability of this software for any purpose. 
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.random.engine.RandomEngine;
+import org.apache.mahout.math.jet.stat.Probability;
+
+/** @deprecated until unit tests are in place.  Until this time, this class/interface is unsupported. */
+@Deprecated
+public class ChiSquare extends AbstractContinousDistribution {
+
+  private double freedom;
+
+  // cached vars for method nextDouble(a) (for performance only)
+  private double freedom_in = -1.0;
+  private double b;
+  private double vm;
+  private double vd;
+
+  // The uniform random number generated shared by all <b>static</b> methods.
+  private static final ChiSquare shared = new ChiSquare(1.0, makeDefaultGenerator());
+
+  /**
+   * Constructs a ChiSquare distribution. Example: freedom=1.0.
+   *
+   * @param freedom degrees of freedom.
+   * @throws IllegalArgumentException if <tt>freedom &lt; 1.0</tt>.
+   */
+  public ChiSquare(double freedom, RandomEngine randomGenerator) {
+    setRandomGenerator(randomGenerator);
+    setState(freedom);
+  }
+
+  /** Returns the cumulative distribution function. */
+  public double cdf(double x) {
+    return Probability.chiSquare(freedom, x);
+  }
+
+  /** Returns a random number from the distribution. */
+  @Override
+  public double nextDouble() {
+    return nextDouble(this.freedom);
+  }
+
+  /**
+   * Returns a random number from the distribution; bypasses the internal state.
+   *
+   * @param freedom degrees of freedom. It should hold <tt>freedom &lt; 1.0</tt>.
+   */
+  public double nextDouble(double freedom) {
+/******************************************************************
+ *                                                                *
+ *        Chi Distribution - Ratio of Uniforms  with shift        *
+ *                                                                *
+ ******************************************************************
+ *                                                                *
+ * FUNCTION :   - chru samples a random number from the Chi       *
+ *                distribution with parameter  a > 1.             *
+ * REFERENCE :  - J.F. Monahan (1987): An algorithm for           *
+ *                generating chi random variables, ACM Trans.     *
+ *                Math. Software 13, 168-172.                     *
+ * SUBPROGRAM : - anEngine  ... pointer to a (0,1)-Uniform        *
+ *                engine                                          *
+ *                                                                *
+ * Implemented by R. Kremer, 1990                                 *
+ ******************************************************************/
+
+    double u, v, z, zz, r;
+
+    //if( a < 1 )  return (-1.0); // Check for invalid input value
+
+    if (freedom == 1.0) {
+      while (true) {
+        u = randomGenerator.raw();
+        v = randomGenerator.raw() * 0.857763884960707;
+        z = v / u;
+        if (z < 0) {
+          continue;
+        }
+        zz = z * z;
+        r = 2.5 - zz;
+        if (z < 0.0) {
+          r += zz * z / (3.0 * z);
+        }
+        if (u < r * 0.3894003915) {
+          return (z * z);
+        }
+        if (zz > (1.036961043 / u + 1.4)) {
+          continue;
+        }
+        if (2.0 * Math.log(u) < (-zz * 0.5)) {
+          return (z * z);
+        }
+      }
+    } else {
+      if (freedom != freedom_in) {
+        b = Math.sqrt(freedom - 1.0);
+        vm = -0.6065306597 * (1.0 - 0.25 / (b * b + 1.0));
+        vm = (-b > vm) ? -b : vm;
+        double vp = 0.6065306597 * (0.7071067812 + b) / (0.5 + b);
+        vd = vp - vm;
+        freedom_in = freedom;
+      }
+      while (true) {
+        u = randomGenerator.raw();
+        v = randomGenerator.raw() * vd + vm;
+        z = v / u;
+        if (z < -b) {
+          continue;
+        }
+        zz = z * z;
+        r = 2.5 - zz;
+        if (z < 0.0) {
+          r += zz * z / (3.0 * (z + b));
+        }
+        if (u < r * 0.3894003915) {
+          return ((z + b) * (z + b));
+        }
+        if (zz > (1.036961043 / u + 1.4)) {
+          continue;
+        }
+        if (2.0 * Math.log(u) < (Math.log(1.0 + z / b) * b * b - zz * 0.5 - z * b)) {
+          return ((z + b) * (z + b));
+        }
+      }
+    }
+  }
+
+  /** Returns the probability distribution function. */
+  public double pdf(double x) {
+    if (x <= 0.0) {
+      throw new IllegalArgumentException();
+    }
+    double logGamma = Fun.logGamma(freedom / 2.0);
+    return Math.exp((freedom / 2.0 - 1.0) * Math.log(x / 2.0) - x / 2.0 - logGamma) / 2.0;
+  }
+
+  /**
+   * Sets the distribution parameter.
+   *
+   * @param freedom degrees of freedom.
+   * @throws IllegalArgumentException if <tt>freedom &lt; 1.0</tt>.
+   */
+  public void setState(double freedom) {
+    if (freedom < 1.0) {
+      throw new IllegalArgumentException();
+    }
+    this.freedom = freedom;
+  }
+
+  /**
+   * Returns a random number from the distribution.
+   *
+   * @param freedom degrees of freedom.
+   * @throws IllegalArgumentException if <tt>freedom &lt; 1.0</tt>.
+   */
+  public static double staticNextDouble(double freedom) {
+    synchronized (shared) {
+      return shared.nextDouble(freedom);
+    }
+  }
+
+  /** Returns a String representation of the receiver. */
+  public String toString() {
+    return this.getClass().getName() + '(' + freedom + ')';
+  }
+
+}

Propchange: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/ChiSquare.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Distributions.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Distributions.java?rev=891983&view=auto
==============================================================================
--- lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Distributions.java (added)
+++ lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Distributions.java Thu Dec 17 23:22:16 2009
@@ -0,0 +1,322 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
+is hereby granted without fee, provided that the above copyright notice appear in all copies and 
+that both that copyright notice and this permission notice appear in supporting documentation. 
+CERN makes no representations about the suitability of this software for any purpose. 
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.random.engine.RandomEngine;
+
+/** @deprecated until unit tests are in place.  Until this time, this class/interface is unsupported. */
+@Deprecated
+public class Distributions {
+
+  private Distributions() {
+  }
+
+  /**
+   * Returns the probability distribution function of the discrete geometric distribution. <p> <tt>p(k) = p *
+   * (1-p)^k</tt> for <tt> k &gt;= 0</tt>. <p>
+   *
+   * @param k the argument to the probability distribution function.
+   * @param p the parameter of the probability distribution function.
+   */
+  public static double geometricPdf(int k, double p) {
+    if (k < 0) {
+      throw new IllegalArgumentException();
+    }
+    return p * Math.pow(1 - p, k);
+  }
+
+  /**
+   * Returns a random number from the Burr II, VII, VIII, X Distributions. <p> <b>Implementation:</b> Inversion method.
+   * This is a port of <tt>burr1.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND /
+   * WIN-RAND</A> library. C-RAND's implementation, in turn, is based upon <p> L. Devroye (1986): Non-Uniform Random
+   * Variate Generation, Springer Verlag, New York. <p>
+   *
+   * @param r  must be &gt; 0.
+   * @param nr the number of the burr distribution (e.g. 2,7,8,10).
+   */
+  public static double nextBurr1(double r, int nr, RandomEngine randomGenerator) {
+/******************************************************************
+ *                                                                *
+ *        Burr II, VII, VIII, X Distributions - Inversion         *
+ *                                                                *
+ ******************************************************************
+ *                                                                *
+ * FUNCTION :   - burr1 samples a random number from one of the   *
+ *                Burr II, VII, VIII, X distributions with        *
+ *                parameter  r > 0 , where the no. of the         *
+ *                distribution is indicated by a pointer          *
+ *                variable.                                       *
+ * REFERENCE :  - L. Devroye (1986): Non-Uniform Random Variate   *
+ *                Generation, Springer Verlag, New York.          *
+ * SUBPROGRAM : - drand(seed) ... (0,1)-uniform generator with    *
+ *                unsigned long integer *seed.                    *
+ *                                                                *
+ ******************************************************************/
+
+    double y = Math.exp(Math.log(randomGenerator.raw()) / r);
+    switch (nr) {
+      // BURR II
+      case 2:
+        return (-Math.log(1 / y - 1));
+
+      // BURR VII
+      case 7:
+        return (Math.log(2 * y / (2 - 2 * y)) / 2);
+
+      // BURR VIII
+      case 8:
+        return (Math.log(Math.tan(y * Math.PI / 2.0)));
+
+      // BURR X
+      case 10:
+        return (Math.sqrt(-Math.log(1 - y)));
+    }
+    return 0;
+  }
+
+  /**
+   * Returns a random number from the Burr III, IV, V, VI, IX, XII distributions. <p> <b>Implementation:</b> Inversion
+   * method. This is a port of <tt>burr2.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND
+   * / WIN-RAND</A> library. C-RAND's implementation, in turn, is based upon <p> L. Devroye (1986): Non-Uniform Random
+   * Variate Generation, Springer Verlag, New York. <p>
+   *
+   * @param r  must be &gt; 0.
+   * @param k  must be &gt; 0.
+   * @param nr the number of the burr distribution (e.g. 3,4,5,6,9,12).
+   */
+  public static double nextBurr2(double r, double k, int nr, RandomEngine randomGenerator) {
+/******************************************************************
+ *                                                                *
+ *      Burr III, IV, V, VI, IX, XII Distribution - Inversion     *
+ *                                                                *
+ ******************************************************************
+ *                                                                *
+ * FUNCTION :   - burr2 samples a random number from one of the   *
+ *                Burr III, IV, V, VI, IX, XII distributions with *
+ *                parameters r > 0 and k > 0, where the no. of    *
+ *                the distribution is indicated by a pointer      *
+ *                variable.                                       *
+ * REFERENCE :  - L. Devroye (1986): Non-Uniform Random Variate   *
+ *                Generation, Springer Verlag, New York.          *
+ * SUBPROGRAM : - drand(seed) ... (0,1)-Uniform generator with    *
+ *                unsigned long integer *seed.                    *
+ *                                                                *
+ ******************************************************************/
+    double u = randomGenerator.raw();
+    double y = Math.exp(-Math.log(u) / r) - 1.0;
+    switch (nr) {
+      case 3:               // BURR III
+        return (Math.exp(-Math.log(y) / k));      // y^(-1/k)
+
+      case 4:               // BURR IV
+        y = Math.exp(k * Math.log(y)) + 1.0;         // y^k + 1
+        y = k / y;
+        return (y);
+
+      case 5:               // BURR V
+        y = Math.atan(-Math.log(y / k));           // arctan[log(y/k)]
+        return (y);
+
+      case 6:               // BURR VI
+        y = -Math.log(y / k) / r;
+        y = Math.log(y + Math.sqrt(y * y + 1.0));
+        return (y);
+
+      case 9:               // BURR IX
+        y = 1.0 + 2.0 * u / (k * (1.0 - u));
+        y = Math.exp(Math.log(y) / r) - 1.0;         // y^(1/r) -1
+        return Math.log(y);
+
+      case 12:               // BURR XII
+        return Math.exp(Math.log(y) / k);        // y^(1/k)
+    }
+    return 0;
+  }
+
+  /**
+   * Returns a cauchy distributed random number from the standard Cauchy distribution C(0,1). <A
+   * HREF="http://www.cern.ch/RD11/rkb/AN16pp/node25.html#SECTION000250000000000000000"> math definition</A> and <A
+   * HREF="http://www.statsoft.com/textbook/glosc.html#Cauchy Distribution"> animated definition</A>. <p> <tt>p(x) = 1/
+   * (mean*pi * (1+(x/mean)^2))</tt>. <p> <b>Implementation:</b> This is a port of <tt>cin.c</tt> from the <A
+   * HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library. <p>
+   */
+  public static double nextCauchy(RandomEngine randomGenerator) {
+    return Math.tan(Math.PI * randomGenerator.raw());
+  }
+
+  /** Returns an erlang distributed random number with the given variance and mean. */
+  public static double nextErlang(double variance, double mean, RandomEngine randomGenerator) {
+    int k = (int) ((mean * mean) / variance + 0.5);
+    k = (k > 0) ? k : 1;
+    double a = k / mean;
+
+    double prod = 1.0;
+    for (int i = 0; i < k; i++) {
+      prod *= randomGenerator.raw();
+    }
+    return -Math.log(prod) / a;
+  }
+
+  /**
+   * Returns a discrete geometric distributed random number; <A HREF="http://www.statsoft.com/textbook/glosf.html#Geometric
+   * Distribution">Definition</A>. <p> <tt>p(k) = p * (1-p)^k</tt> for <tt> k &gt;= 0</tt>. <p> <b>Implementation:</b>
+   * Inversion method. This is a port of <tt>geo.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND
+   * / WIN-RAND</A> library.
+   *
+   * @param p must satisfy <tt>0 &lt; p &lt; 1</tt>. <p>
+   */
+  public static int nextGeometric(double p, RandomEngine randomGenerator) {
+/******************************************************************
+ *                                                                *
+ *              Geometric Distribution - Inversion                *
+ *                                                                *
+ ******************************************************************
+ *                                                                *
+ * On generating random numbers of a discrete distribution by     *
+ * Inversion normally sequential search is necessary, but in the  *
+ * case of the Geometric distribution a direct transformation is  *
+ * possible because of the special parallel to the continuous     *
+ * Exponential distribution Exp(t):                               *
+ *    X - Exp(t): G(x)=1-exp(-tx)                                 *
+ *        Geo(p): pk=G(k+1)-G(k)=exp(-tk)*(1-exp(-t))             *
+ *                p=1-exp(-t)                                     *
+ * A random number of the Geometric distribution Geo(p) is        *
+ * obtained by k=(long int)x, where x is from Exp(t) with         *
+ * parameter t=-log(1-p).                                         *
+ *                                                                *
+ ******************************************************************
+ *                                                                *
+ * FUNCTION:    - geo samples a random number from the Geometric  *
+ *                distribution with parameter 0<p<1.              *
+ * SUBPROGRAMS: - drand(seed) ... (0,1)-Uniform generator with    *
+ *                unsigned long integer *seed.                    *
+ *                                                                *
+ ******************************************************************/
+    double u = randomGenerator.raw();
+    return (int) (Math.log(u) / Math.log(1.0 - p));
+  }
+
+  /**
+   * Returns a lambda distributed random number with parameters l3 and l4. <p> <b>Implementation:</b> Inversion method.
+   * This is a port of <tt>lamin.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND /
+   * WIN-RAND</A> library. C-RAND's implementation, in turn, is based upon <p> J.S. Ramberg, B:W. Schmeiser (1974): An
+   * approximate method for generating asymmetric variables, Communications ACM 17, 78-82. <p>
+   */
+  public static double nextLambda(double l3, double l4, RandomEngine randomGenerator) {
+    double l_sign;
+    if ((l3 < 0) || (l4 < 0)) {
+      l_sign = -1.0;                          // sign(l)
+    } else {
+      l_sign = 1.0;
+    }
+
+    double u = randomGenerator.raw();                           // U(0/1)
+    double x = l_sign * (Math.exp(Math.log(u) * l3) - Math.exp(Math.log(1.0 - u) * l4));
+    return x;
+  }
+
+  /**
+   * Returns a Laplace (Double Exponential) distributed random number from the standard Laplace distribution L(0,1). <p>
+   * <b>Implementation:</b> Inversion method. This is a port of <tt>lapin.c</tt> from the <A
+   * HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library. <p>
+   */
+  public static double nextLaplace(RandomEngine randomGenerator) {
+    double u = randomGenerator.raw();
+    u = u + u - 1.0;
+    if (u > 0) {
+      return -Math.log(1.0 - u);
+    } else {
+      return Math.log(1.0 + u);
+    }
+  }
+
+  /**
+   * Returns a random number from the standard Logistic distribution Log(0,1). <p> <b>Implementation:</b> Inversion
+   * method. This is a port of <tt>login.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND
+   * / WIN-RAND</A> library.
+   */
+  public static double nextLogistic(RandomEngine randomGenerator) {
+    double u = randomGenerator.raw();
+    return (-Math.log(1.0 / u - 1.0));
+  }
+
+  /**
+   * Returns a power-law distributed random number with the given exponent and lower cutoff.
+   *
+   * @param alpha the exponent
+   * @param cut   the lower cutoff
+   */
+  public static double nextPowLaw(double alpha, double cut, RandomEngine randomGenerator) {
+    return cut * Math.pow(randomGenerator.raw(), 1.0 / (alpha + 1.0));
+  }
+
+  /**
+   * Returns a random number from the standard Triangular distribution in (-1,1). <p> <b>Implementation:</b> Inversion
+   * method. This is a port of <tt>tra.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND
+   * / WIN-RAND</A> library. <p>
+   */
+  public static double nextTriangular(RandomEngine randomGenerator) {
+/******************************************************************
+ *                                                                *
+ *     Triangular Distribution - Inversion: x = +-(1-sqrt(u))     *
+ *                                                                *
+ ******************************************************************
+ *                                                                *
+ * FUNCTION :   - tra samples a random number from the            *
+ *                standard Triangular distribution in (-1,1)      *
+ * SUBPROGRAM : - drand(seed) ... (0,1)-Uniform generator with    *
+ *                unsigned long integer *seed.                    *
+ *                                                                *
+ ******************************************************************/
+
+    double u = randomGenerator.raw();
+    if (u <= 0.5) {
+      return (Math.sqrt(2.0 * u) - 1.0);                      /* -1 <= x <= 0 */
+    } else {
+      return (1.0 - Math.sqrt(2.0 * (1.0 - u)));
+    }                 /*  0 <= x <= 1 */
+  }
+
+  /**
+   * Returns a weibull distributed random number. Polar method. See Simulation, Modelling & Analysis by Law & Kelton,
+   * pp259
+   */
+  public static double nextWeibull(double alpha, double beta, RandomEngine randomGenerator) {
+    // Polar method.
+    // See Simulation, Modelling & Analysis by Law & Kelton, pp259
+    return Math.pow(beta * (-Math.log(1.0 - randomGenerator.raw())), 1.0 / alpha);
+  }
+
+  /**
+   * Returns a zipfian distributed random number with the given skew. <p> Algorithm from page 551 of: Devroye, Luc
+   * (1986) `Non-uniform random variate generation', Springer-Verlag: Berlin.   ISBN 3-540-96305-7 (also 0-387-96305-7)
+   *
+   * @param z the skew of the distribution (must be &gt;1.0).
+   */
+  public static int nextZipfInt(double z, RandomEngine randomGenerator) {
+    /* Algorithm from page 551 of:
+    * Devroye, Luc (1986) `Non-uniform random variate generation',
+    * Springer-Verlag: Berlin.   ISBN 3-540-96305-7 (also 0-387-96305-7)
+    */
+    double b = Math.pow(2.0, z - 1.0);
+    double constant = -1.0 / (z - 1.0);
+
+    int result;
+    while (true) {
+      double u = randomGenerator.raw();
+      double v = randomGenerator.raw();
+      result = (int) (Math.floor(Math.pow(u, constant)));
+      double t = Math.pow(1.0 + 1.0 / result, z - 1.0);
+      if (v * result * (t - 1.0) / (b - 1.0) <= t / b) {
+        break;
+      }
+    }
+    return result;
+  }
+}

Propchange: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Distributions.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Empirical.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Empirical.java?rev=891983&view=auto
==============================================================================
--- lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Empirical.java (added)
+++ lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Empirical.java Thu Dec 17 23:22:16 2009
@@ -0,0 +1,180 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
+is hereby granted without fee, provided that the above copyright notice appear in all copies and 
+that both that copyright notice and this permission notice appear in supporting documentation. 
+CERN makes no representations about the suitability of this software for any purpose. 
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.random.engine.RandomEngine;
+
+/** @deprecated until unit tests are in place.  Until this time, this class/interface is unsupported. */
+@Deprecated
+public class Empirical extends AbstractContinousDistribution {
+
+  private double[] cdf; // cumulative distribution function
+  private int interpolationType;
+
+  private static final int LINEAR_INTERPOLATION = 0;
+  private static final int NO_INTERPOLATION = 1;
+
+  /**
+   * Constructs an Empirical distribution. The probability distribution function (pdf) is an array of positive real
+   * numbers. It need not be provided in the form of relative probabilities, absolute probabilities are also accepted.
+   * The <tt>pdf</tt> must satisfy both of the following conditions <ul> <li><tt>0.0 &lt;= pdf[i] :
+   * 0&lt;=i&lt;=pdf.length-1</tt> <li><tt>0.0 &lt; Sum(pdf[i]) : 0&lt;=i&lt;=pdf.length-1</tt> </ul>
+   *
+   * @param pdf               the probability distribution function.
+   * @param interpolationType can be either <tt>Empirical.NO_INTERPOLATION</tt> or <tt>Empirical.LINEAR_INTERPOLATION</tt>.
+   * @param randomGenerator   a uniform random number generator.
+   * @throws IllegalArgumentException if at least one of the three conditions above is violated.
+   */
+  public Empirical(double[] pdf, int interpolationType, RandomEngine randomGenerator) {
+    setRandomGenerator(randomGenerator);
+    setState(pdf, interpolationType);
+  }
+
+  /** Returns the cumulative distribution function. */
+  public double cdf(int k) {
+    if (k < 0) {
+      return 0.0;
+    }
+    if (k >= cdf.length - 1) {
+      return 1.0;
+    }
+    return cdf[k];
+  }
+
+  /**
+   * Returns a deep copy of the receiver; the copy will produce identical sequences. After this call has returned, the
+   * copy and the receiver have equal but separate state.
+   *
+   * @return a copy of the receiver.
+   */
+  @Override
+  public Object clone() {
+    Empirical copy = (Empirical) super.clone();
+    if (this.cdf != null) {
+      copy.cdf = this.cdf.clone();
+    }
+    return copy;
+  }
+
+  /** Returns a random number from the distribution. */
+  @Override
+  public double nextDouble() {
+    double rand = randomGenerator.raw();
+    if (this.cdf == null) {
+      return rand;
+    } // Non-existing pdf
+
+    // binary search in cumulative distribution function:
+    int nBins = cdf.length - 1;
+    int nbelow = 0;     // largest k such that I[k] is known to be <= rand
+    int nabove = nBins; // largest k such that I[k] is known to be >  rand
+
+    while (nabove > nbelow + 1) {
+      int middle = (nabove + nbelow + 1) >> 1; // div 2
+      if (rand >= cdf[middle]) {
+        nbelow = middle;
+      } else {
+        nabove = middle;
+      }
+    }
+    // after this binary search, nabove is always nbelow+1 and they straddle rand:
+
+    if (this.interpolationType == NO_INTERPOLATION) {
+      return ((double) nbelow) / nBins;
+    } else if (this.interpolationType == LINEAR_INTERPOLATION) {
+      double binMeasure = cdf[nabove] - cdf[nbelow];
+      // binMeasure is always aProbFunc[nbelow],
+      // but we don't have aProbFunc any more so we subtract.
+
+      if (binMeasure == 0.0) {
+        // rand lies right in a bin of measure 0.  Simply return the center
+        // of the range of that bin.  (Any value between k/N and (k+1)/N is
+        // equally good, in this rare case.)
+        return (nbelow + 0.5) / nBins;
+      }
+
+      double binFraction = (rand - cdf[nbelow]) / binMeasure;
+      return (nbelow + binFraction) / nBins;
+    } else {
+      throw new InternalError();
+    } // illegal interpolation type
+  }
+
+  /** Returns the probability distribution function. */
+  public double pdf(double x) {
+    throw new UnsupportedOperationException("not implemented");
+    //if (x < 0 || x > cdf.length-2) return 0.0;
+    //int k = (int) x;
+    //return cdf[k-1] - cdf[k];
+  }
+
+  /** Returns the probability distribution function. */
+  public double pdf(int k) {
+    if (k < 0 || k >= cdf.length - 1) {
+      return 0.0;
+    }
+    return cdf[k - 1] - cdf[k];
+  }
+
+  /**
+   * Sets the distribution parameters. The <tt>pdf</tt> must satisfy both of the following conditions <ul> <li><tt>0.0
+   * &lt;= pdf[i] : 0 &lt; =i &lt;= pdf.length-1</tt> <li><tt>0.0 &lt; Sum(pdf[i]) : 0 &lt;=i &lt;= pdf.length-1</tt>
+   * </ul>
+   *
+   * @param pdf               probability distribution function.
+   * @param interpolationType can be either <tt>Empirical.NO_INTERPOLATION</tt> or <tt>Empirical.LINEAR_INTERPOLATION</tt>.
+   * @throws IllegalArgumentException if at least one of the three conditions above is violated.
+   */
+  public void setState(double[] pdf, int interpolationType) {
+    if (interpolationType != LINEAR_INTERPOLATION &&
+        interpolationType != NO_INTERPOLATION) {
+      throw new IllegalArgumentException("Illegal Interpolation Type");
+    }
+    this.interpolationType = interpolationType;
+
+    if (pdf == null || pdf.length == 0) {
+      this.cdf = null;
+      //throw new IllegalArgumentException("Non-existing pdf");
+      return;
+    }
+
+    // compute cumulative distribution function (cdf) from probability distribution function (pdf)
+    int nBins = pdf.length;
+    this.cdf = new double[nBins + 1];
+
+    cdf[0] = 0;
+    for (int ptn = 0; ptn < nBins; ++ptn) {
+      double prob = pdf[ptn];
+      if (prob < 0.0) {
+        throw new IllegalArgumentException("Negative probability");
+      }
+      cdf[ptn + 1] = cdf[ptn] + prob;
+    }
+    if (cdf[nBins] <= 0.0) {
+      throw new IllegalArgumentException("At leat one probability must be > 0.0");
+    }
+    for (int ptn = 0; ptn < nBins + 1; ++ptn) {
+      cdf[ptn] /= cdf[nBins];
+    }
+    // cdf is now cached...
+  }
+
+  /** Returns a String representation of the receiver. */
+  public String toString() {
+    String interpolation = null;
+    if (interpolationType == NO_INTERPOLATION) {
+      interpolation = "NO_INTERPOLATION";
+    }
+    if (interpolationType == LINEAR_INTERPOLATION) {
+      interpolation = "LINEAR_INTERPOLATION";
+    }
+    return this.getClass().getName() + '(' + ((cdf != null) ? cdf.length : 0) + ',' + interpolation + ')';
+  }
+
+}

Propchange: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Empirical.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/EmpiricalWalker.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/EmpiricalWalker.java?rev=891983&view=auto
==============================================================================
--- lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/EmpiricalWalker.java (added)
+++ lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/EmpiricalWalker.java Thu Dec 17 23:22:16 2009
@@ -0,0 +1,332 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
+is hereby granted without fee, provided that the above copyright notice appear in all copies and 
+that both that copyright notice and this permission notice appear in supporting documentation. 
+CERN makes no representations about the suitability of this software for any purpose. 
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.random.engine.RandomEngine;
+
+/** @deprecated until unit tests are in place.  Until this time, this class/interface is unsupported. */
+@Deprecated
+public class EmpiricalWalker extends AbstractDiscreteDistribution {
+
+  private int K;
+  private int[] A;
+  private double[] F;
+
+  private double[] cdf; // cumulative distribution function
+
+  /*
+  * James Theiler, jt@lanl.gov, the author of the GSL routine this port is based on, describes his approach as follows:
+  *
+  * Based on: Alastair J Walker, An efficient method for generating
+  * discrete random variables with general distributions, ACM Trans
+  * Math Soft 3, 253-256 (1977).  See also: D. E. Knuth, The Art of
+  * Computer Programming, Volume 2 (Seminumerical algorithms), 3rd
+  * edition, Addison-Wesley (1997), p120.
+
+  * Walker's algorithm does some preprocessing, and provides two
+  * arrays: floating point F[k] and integer A[k].  A value k is chosen
+  * from 0..K-1 with equal likelihood, and then a uniform random number
+  * u is compared to F[k].  If it is less than F[k], then k is
+  * returned.  Otherwise, A[k] is returned.
+
+  * Walker's original paper describes an O(K^2) algorithm for setting
+  * up the F and A arrays.  I found this disturbing since I wanted to
+  * use very large values of K.  I'm sure I'm not the first to realize
+  * this, but in fact the preprocessing can be done in O(K) steps.
+
+  * A figure of merit for the preprocessing is the average value for
+  * the F[k]'s (that is, SUM_k F[k]/K); this corresponds to the
+  * probability that k is returned, instead of A[k], thereby saving a
+  * redirection.  Walker's O(K^2) preprocessing will generally improve
+  * that figure of merit, compared to my cheaper O(K) method; from some
+  * experiments with a perl script, I get values of around 0.6 for my
+  * method and just under 0.75 for Walker's.  Knuth has pointed out
+  * that finding _the_ optimum lookup tables, which maximize the
+  * average F[k], is a combinatorially difficult problem.  But any
+  * valid preprocessing will still provide O(1) time for the call to
+  * gsl_ran_discrete().  I find that if I artificially set F[k]=1 --
+  * ie, better than optimum! -- I get a speedup of maybe 20%, so that's
+  * the maximum I could expect from the most expensive preprocessing.
+  * Folding in the difference of 0.6 vs 0.75, I'd estimate that the
+  * speedup would be less than 10%.
+
+  * I've not implemented it here, but one compromise is to sort the
+  * probabilities once, and then work from the two ends inward.  This
+  * requires O(K log K), still lots cheaper than O(K^2), and from my
+  * experiments with the perl script, the figure of merit is within
+  * about 0.01 for K up to 1000, and no sign of diverging (in fact,
+  * they seemed to be converging, but it's hard to say with just a
+  * handful of runs).
+
+  * The O(K) algorithm goes through all the p_k's and decides if they
+  * are "smalls" or "bigs" according to whether they are less than or
+  * greater than the mean value 1/K.  The indices to the smalls and the
+  * bigs are put in separate stacks, and then we work through the
+  * stacks together.  For each small, we pair it up with the next big
+  * in the stack (Walker always wanted to pair up the smallest small
+  * with the biggest big).  The small "borrows" from the big just
+  * enough to bring the small up to mean.  This reduces the size of the
+  * big, so the (smaller) big is compared again to the mean, and if it
+  * is smaller, it gets "popped" from the big stack and "pushed" to the
+  * small stack.  Otherwise, it stays put.  Since every time we pop a
+  * small, we are able to deal with it right then and there, and we
+  * never have to pop more than K smalls, then the algorithm is O(K).
+
+  * This implementation sets up two separate stacks, and allocates K
+  * elements between them.  Since neither stack ever grows, we do an
+  * extra O(K) pass through the data to determine how many smalls and
+  * bigs there are to begin with and allocate appropriately.  In all
+  * there are 2*K*sizeof(double) transient bytes of memory that are
+  * used than returned, and K*(sizeof(int)+sizeof(double)) bytes used
+  * in the lookup table.
+
+  * Walker spoke of using two random numbers (an integer 0..K-1, and a
+  * floating point u in [0,1]), but Knuth points out that one can just
+  * use the integer and fractional parts of K*u where u is in [0,1].
+  * In fact, Knuth further notes that taking F'[k]=(k+F[k])/K, one can
+  * directly compare u to F'[k] without having to explicitly set
+  * u=K*u-int(K*u).
+
+  * Usage:
+
+  * Starting with an array of probabilities P, initialize and do
+  * preprocessing with a call to:
+
+  *    gsl_rng *r;
+  *    gsl_ran_discrete_t *f;
+  *    f = gsl_ran_discrete_preproc(K,P);
+
+  * Then, whenever a random index 0..K-1 is desired, use
+
+  *    k = gsl_ran_discrete(r,f);
+
+  * Note that several different randevent struct's can be
+  * simultaneously active.
+
+  * Aside: A very clever alternative approach is described in
+  * Abramowitz and Stegun, p 950, citing: Marsaglia, Random variables
+  * and computers, Proc Third Prague Conference in Probability Theory,
+  * 1962.  A more accesible reference is: G. Marsaglia, Generating
+  * discrete random numbers in a computer, Comm ACM 6, 37-38 (1963).
+  * If anybody is interested, I (jt) have also coded up this version as
+  * part of another software package.  However, I've done some
+  * comparisons, and the Walker method is both faster and more stingy
+  * with memory.  So, in the end I decided not to include it with the
+  * GSL package.
+
+  * Written 26 Jan 1999, James Theiler, jt@lanl.gov
+  * Adapted to GSL, 30 Jan 1999, jt
+
+  */
+
+  /**
+   * Constructs an Empirical distribution. The probability distribution function (pdf) is an array of positive real
+   * numbers. It need not be provided in the form of relative probabilities, absolute probabilities are also accepted.
+   * The <tt>pdf</tt> must satisfy both of the following conditions <ul> <li><tt>0.0 &lt;= pdf[i] :
+   * 0&lt;=i&lt;=pdf.length-1</tt> <li><tt>0.0 &lt; Sum(pdf[i]) : 0&lt;=i&lt;=pdf.length-1</tt> </ul>
+   *
+   * @param pdf               the probability distribution function.
+   * @param interpolationType can be either <tt>Empirical.NO_INTERPOLATION</tt> or <tt>Empirical.LINEAR_INTERPOLATION</tt>.
+   * @param randomGenerator   a uniform random number generator.
+   * @throws IllegalArgumentException if at least one of the three conditions above is violated.
+   */
+  public EmpiricalWalker(double[] pdf, int interpolationType, RandomEngine randomGenerator) {
+    setRandomGenerator(randomGenerator);
+    setState(pdf, interpolationType);
+    setState2(pdf);
+  }
+
+  /** Returns the cumulative distribution function. */
+  public double cdf(int k) {
+    if (k < 0) {
+      return 0.0;
+    }
+    if (k >= cdf.length - 1) {
+      return 1.0;
+    }
+    return cdf[k];
+  }
+
+  /**
+   * Returns a deep copy of the receiver; the copy will produce identical sequences. After this call has returned, the
+   * copy and the receiver have equal but separate state.
+   *
+   * @return a copy of the receiver.
+   */
+  @Override
+  public Object clone() {
+    EmpiricalWalker copy = (EmpiricalWalker) super.clone();
+    if (this.cdf != null) {
+      copy.cdf = this.cdf.clone();
+    }
+    if (this.A != null) {
+      copy.A = this.A.clone();
+    }
+    if (this.F != null) {
+      copy.F = this.F.clone();
+    }
+    return copy;
+  }
+
+  /** Returns a random integer <tt>k</tt> with probability <tt>pdf(k)</tt>. */
+  @Override
+  public int nextInt() {
+    double u = this.randomGenerator.raw();
+    u *= this.K;
+    int c = (int) u;
+    u -= c;
+    double f = this.F[c];
+    if (f == 1.0) {
+      return c;
+    }
+    if (u < f) {
+      return c;
+    } else {
+      return this.A[c];
+    }
+  }
+
+  /** Returns the probability distribution function. */
+  public double pdf(int k) {
+    if (k < 0 || k >= cdf.length - 1) {
+      return 0.0;
+    }
+    return cdf[k - 1] - cdf[k];
+  }
+
+  /**
+   * Sets the distribution parameters. The <tt>pdf</tt> must satisfy all of the following conditions <ul> <li><tt>pdf !=
+   * null && pdf.length &gt; 0</tt> <li><tt>0.0 &lt;= pdf[i] : 0 &lt; =i &lt;= pdf.length-1</tt> <li><tt>0.0 &lt;
+   * Sum(pdf[i]) : 0 &lt;=i &lt;= pdf.length-1</tt> </ul>
+   *
+   * @param pdf probability distribution function.
+   * @throws IllegalArgumentException if at least one of the three conditions above is violated.
+   */
+  public void setState(double[] pdf, int interpolationType) {
+    if (pdf == null || pdf.length == 0) {
+      throw new IllegalArgumentException("Non-existing pdf");
+    }
+
+    // compute cumulative distribution function (cdf) from probability distribution function (pdf)
+    int nBins = pdf.length;
+    this.cdf = new double[nBins + 1];
+
+    cdf[0] = 0;
+    for (int i = 0; i < nBins; i++) {
+      if (pdf[i] < 0.0) {
+        throw new IllegalArgumentException("Negative probability");
+      }
+      cdf[i + 1] = cdf[i] + pdf[i];
+    }
+    if (cdf[nBins] <= 0.0) {
+      throw new IllegalArgumentException("At leat one probability must be > 0.0");
+    }
+    // now normalize to 1 (relative probabilities).
+    for (int i = 0; i < nBins + 1; i++) {
+      cdf[i] /= cdf[nBins];
+    }
+    // cdf is now cached...
+  }
+
+  /**
+   * Sets the distribution parameters. The <tt>pdf</tt> must satisfy both of the following conditions <ul> <li><tt>0.0
+   * &lt;= pdf[i] : 0 &lt; =i &lt;= pdf.length-1</tt> <li><tt>0.0 &lt; Sum(pdf[i]) : 0 &lt;=i &lt;= pdf.length-1</tt>
+   * </ul>
+   *
+   * @param pdf probability distribution function.
+   * @throws IllegalArgumentException if at least one of the three conditions above is violated.
+   */
+  public void setState2(double[] pdf) {
+    int size = pdf.length;
+    int k;
+    double pTotal = 0;
+
+    //if (size < 1) {
+    //  throw new IllegalArgumentException("must have size greater than zero");
+    //}
+    /* Make sure elements of ProbArray[] are positive.
+    * Won't enforce that sum is unity; instead will just normalize
+    */
+    for (k = 0; k < size; ++k) {
+      //if (pdf[k] < 0) {
+      //throw new IllegalArgumentException("Probabilities must be >= 0: "+pdf[k]);
+      //}
+      pTotal += pdf[k];
+    }
+
+    /* Begin setting up the internal state */
+    this.K = size;
+    this.F = new double[size];
+    this.A = new int[size];
+
+    // normalize to relative probability
+    double[] E = new double[size];
+    for (k = 0; k < size; ++k) {
+      E[k] = pdf[k] / pTotal;
+    }
+
+    /* Now create the Bigs and the Smalls */
+    double mean = 1.0 / size;
+    int nSmalls = 0;
+    int nBigs = 0;
+    for (k = 0; k < size; ++k) {
+      if (E[k] < mean) {
+        ++nSmalls;
+      } else {
+        ++nBigs;
+      }
+    }
+    Stack Bigs = new Stack(nBigs);
+    Stack Smalls = new Stack(nSmalls);
+    for (k = 0; k < size; ++k) {
+      if (E[k] < mean) {
+        Smalls.push(k);
+      } else {
+        Bigs.push(k);
+      }
+    }
+    /* Now work through the smalls */
+    int b;
+    while (Smalls.size() > 0) {
+      int s = Smalls.pop();
+      if (Bigs.size() == 0) {
+        /* Then we are on our last value */
+        this.A[s] = s;
+        this.F[s] = 1.0;
+        break;
+      }
+      b = Bigs.pop();
+      this.A[s] = b;
+      this.F[s] = size * E[s];
+
+      double d = mean - E[s];
+      E[s] += d;              /* now E[s] == mean */
+      E[b] -= d;
+      if (E[b] < mean) {
+        Smalls.push(b); /* no longer big, join ranks of the small */
+      } else if (E[b] > mean) {
+        Bigs.push(b); /* still big, put it back where you found it */
+      } else {
+        /* E[b]==mean implies it is finished too */
+        this.A[b] = b;
+        this.F[b] = 1.0;
+      }
+    }
+    while (Bigs.size() > 0) {
+      b = Bigs.pop();
+      this.A[b] = b;
+      this.F[b] = 1.0;
+    }
+  }
+
+  /** Returns a String representation of the receiver. */
+  public String toString() {
+    return this.getClass().getName() + '(' + ((cdf != null) ? cdf.length : 0) + ')';
+  }
+}

Propchange: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/EmpiricalWalker.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Exponential.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Exponential.java?rev=891983&view=auto
==============================================================================
--- lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Exponential.java (added)
+++ lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Exponential.java Thu Dec 17 23:22:16 2009
@@ -0,0 +1,72 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
+is hereby granted without fee, provided that the above copyright notice appear in all copies and 
+that both that copyright notice and this permission notice appear in supporting documentation. 
+CERN makes no representations about the suitability of this software for any purpose. 
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.random.engine.RandomEngine;
+
+/** @deprecated until unit tests are in place.  Until this time, this class/interface is unsupported. */
+@Deprecated
+public class Exponential extends AbstractContinousDistribution {
+
+  private double lambda;
+
+  // The uniform random number generated shared by all <b>static</b> methods.
+  private static final Exponential shared = new Exponential(1.0, makeDefaultGenerator());
+
+  /** Constructs a Negative Exponential distribution. */
+  public Exponential(double lambda, RandomEngine randomGenerator) {
+    setRandomGenerator(randomGenerator);
+    setState(lambda);
+  }
+
+  /** Returns the cumulative distribution function. */
+  public double cdf(double x) {
+    if (x <= 0.0) {
+      return 0.0;
+    }
+    return 1.0 - Math.exp(-x * lambda);
+  }
+
+  /** Returns a random number from the distribution. */
+  @Override
+  public double nextDouble() {
+    return nextDouble(lambda);
+  }
+
+  /** Returns a random number from the distribution; bypasses the internal state. */
+  public double nextDouble(double lambda) {
+    return -Math.log(randomGenerator.raw()) / lambda;
+  }
+
+  /** Returns the probability distribution function. */
+  public double pdf(double x) {
+    if (x < 0.0) {
+      return 0.0;
+    }
+    return lambda * Math.exp(-x * lambda);
+  }
+
+  /** Sets the mean. */
+  public void setState(double lambda) {
+    this.lambda = lambda;
+  }
+
+  /** Returns a random number from the distribution with the given lambda. */
+  public static double staticNextDouble(double lambda) {
+    synchronized (shared) {
+      return shared.nextDouble(lambda);
+    }
+  }
+
+  /** Returns a String representation of the receiver. */
+  public String toString() {
+    return this.getClass().getName() + '(' + lambda + ')';
+  }
+
+}

Propchange: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Exponential.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/ExponentialPower.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/ExponentialPower.java?rev=891983&view=auto
==============================================================================
--- lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/ExponentialPower.java (added)
+++ lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/ExponentialPower.java Thu Dec 17 23:22:16 2009
@@ -0,0 +1,113 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
+is hereby granted without fee, provided that the above copyright notice appear in all copies and 
+that both that copyright notice and this permission notice appear in supporting documentation. 
+CERN makes no representations about the suitability of this software for any purpose. 
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.random.engine.RandomEngine;
+
+/** @deprecated until unit tests are in place.  Until this time, this class/interface is unsupported. */
+@Deprecated
+public class ExponentialPower extends AbstractContinousDistribution {
+
+  private double tau;
+
+  // cached vars for method nextDouble(tau) (for performance only)
+  private double s, sm1, tau_set = -1.0;
+
+  // The uniform random number generated shared by all <b>static</b> methods.
+  private static final ExponentialPower shared = new ExponentialPower(1.0, makeDefaultGenerator());
+
+  /**
+   * Constructs an Exponential Power distribution. Example: tau=1.0.
+   *
+   * @throws IllegalArgumentException if <tt>tau &lt; 1.0</tt>.
+   */
+  public ExponentialPower(double tau, RandomEngine randomGenerator) {
+    setRandomGenerator(randomGenerator);
+    setState(tau);
+  }
+
+  /** Returns a random number from the distribution. */
+  @Override
+  public double nextDouble() {
+    return nextDouble(this.tau);
+  }
+
+  /**
+   * Returns a random number from the distribution; bypasses the internal state.
+   *
+   * @throws IllegalArgumentException if <tt>tau &lt; 1.0</tt>.
+   */
+  public double nextDouble(double tau) {
+
+    if (tau != tau_set) { // SET-UP
+      s = 1.0 / tau;
+      sm1 = 1.0 - s;
+
+      tau_set = tau;
+    }
+
+    // GENERATOR
+    double x;
+    double v;
+    double u;
+    do {
+      u = randomGenerator.raw();                             // U(0/1)
+      u = (2.0 * u) - 1.0;                                     // U(-1.0/1.0)
+      double u1 = Math.abs(u);
+      v = randomGenerator.raw();                             // U(0/1)
+
+      if (u1 <= sm1) { // Uniform hat-function for x <= (1-1/tau)
+        x = u1;
+      } else { // Exponential hat-function for x > (1-1/tau)
+        double y = tau * (1.0 - u1);
+        x = sm1 - s * Math.log(y);
+        v *= y;
+      }
+    }
+
+    // Acceptance/Rejection
+    while (Math.log(v) > -Math.exp(Math.log(x) * tau));
+
+    // Random sign
+    if (u < 0.0) {
+      return x;
+    } else {
+      return -x;
+    }
+  }
+
+  /**
+   * Sets the distribution parameter.
+   *
+   * @throws IllegalArgumentException if <tt>tau &lt; 1.0</tt>.
+   */
+  public void setState(double tau) {
+    if (tau < 1.0) {
+      throw new IllegalArgumentException();
+    }
+    this.tau = tau;
+  }
+
+  /**
+   * Returns a random number from the distribution.
+   *
+   * @throws IllegalArgumentException if <tt>tau &lt; 1.0</tt>.
+   */
+  public static double staticNextDouble(double tau) {
+    synchronized (shared) {
+      return shared.nextDouble(tau);
+    }
+  }
+
+  /** Returns a String representation of the receiver. */
+  public String toString() {
+    return this.getClass().getName() + '(' + tau + ')';
+  }
+
+}

Propchange: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/ExponentialPower.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Fun.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Fun.java?rev=891983&view=auto
==============================================================================
--- lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Fun.java (added)
+++ lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Fun.java Thu Dec 17 23:22:16 2009
@@ -0,0 +1,430 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
+is hereby granted without fee, provided that the above copyright notice appear in all copies and 
+that both that copyright notice and this permission notice appear in supporting documentation. 
+CERN makes no representations about the suitability of this software for any purpose. 
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.math.Arithmetic;
+
+/**
+ * Contains various mathematical helper methods.
+ *
+ * <b>Implementation:</b> High performance implementation. <dt>This is a port of <tt>gen_fun.cpp</tt> from the <A
+ * HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
+ */
+class Fun {
+
+  private Fun() {
+  }
+
+  private static double _fkt_value(double lambda, double z1, double z2, double x_value) {
+
+    double y_value = Math.cos(z1 * x_value) / (Math.pow((x_value * x_value + z2 * z2), (lambda + 0.5)));
+    return y_value;
+  }
+
+  public static double bessel2_fkt(double lambda, double beta) {
+
+    double[] b0 = {-1.5787132, -0.6130827, 0.1735823, 1.4793411,
+        2.6667307, 4.9086836, 8.1355339,
+    };
+    double[] b05 = {-1.9694802, -0.7642538, 0.0826017, 1.4276355,
+        2.6303682, 4.8857787, 8.1207968,
+    };
+    double[] b1 = {-2.9807345, -1.1969943, -0.1843161, 1.2739241,
+        2.5218256, 4.8172216, 8.0765633,
+    };
+    double[] b2 = {-5.9889676, -2.7145389, -1.1781269, 0.6782201,
+        2.0954009, 4.5452152, 7.9003173,
+    };
+    double[] b3 = {-9.6803440, -4.8211925, -2.6533185, -0.2583337,
+        1.4091915, 4.0993448, 7.6088310,
+    };
+    double[] b5 = {-18.1567152, -10.0939408, -6.5819139, -2.9371545,
+        -0.6289005, 2.7270412, 6.6936799,
+    };
+    double[] b8 = {-32.4910195, -19.6065943, -14.0347298, -8.3839439,
+        -4.9679730, -0.3567823, 4.5589697,
+    };
+
+    if (lambda == 0.0) {
+      if (beta == 0.1) {
+        return (b0[0]);
+      }
+      if (beta == 0.5) {
+        return (b0[1]);
+      }
+      if (beta == 1.0) {
+        return (b0[2]);
+      }
+      if (beta == 2.0) {
+        return (b0[3]);
+      }
+      if (beta == 3.0) {
+        return (b0[4]);
+      }
+      if (beta == 5.0) {
+        return (b0[5]);
+      }
+      if (beta == 8.0) {
+        return (b0[6]);
+      }
+    }
+
+    if (lambda == 0.5) {
+      if (beta == 0.1) {
+        return (b05[0]);
+      }
+      if (beta == 0.5) {
+        return (b05[1]);
+      }
+      if (beta == 1.0) {
+        return (b05[2]);
+      }
+      if (beta == 2.0) {
+        return (b05[3]);
+      }
+      if (beta == 3.0) {
+        return (b05[4]);
+      }
+      if (beta == 5.0) {
+        return (b05[5]);
+      }
+      if (beta == 8.0) {
+        return (b05[6]);
+      }
+    }
+
+    if (lambda == 1.0) {
+      if (beta == 0.1) {
+        return (b1[0]);
+      }
+      if (beta == 0.5) {
+        return (b1[1]);
+      }
+      if (beta == 1.0) {
+        return (b1[2]);
+      }
+      if (beta == 2.0) {
+        return (b1[3]);
+      }
+      if (beta == 3.0) {
+        return (b1[4]);
+      }
+      if (beta == 5.0) {
+        return (b1[5]);
+      }
+      if (beta == 8.0) {
+        return (b1[6]);
+      }
+    }
+
+    if (lambda == 2.0) {
+      if (beta == 0.1) {
+        return (b2[0]);
+      }
+      if (beta == 0.5) {
+        return (b2[1]);
+      }
+      if (beta == 1.0) {
+        return (b2[2]);
+      }
+      if (beta == 2.0) {
+        return (b2[3]);
+      }
+      if (beta == 3.0) {
+        return (b2[4]);
+      }
+      if (beta == 5.0) {
+        return (b2[5]);
+      }
+      if (beta == 8.0) {
+        return (b2[6]);
+      }
+    }
+
+    if (lambda == 3.0) {
+      if (beta == 0.1) {
+        return (b3[0]);
+      }
+      if (beta == 0.5) {
+        return (b3[1]);
+      }
+      if (beta == 1.0) {
+        return (b3[2]);
+      }
+      if (beta == 2.0) {
+        return (b3[3]);
+      }
+      if (beta == 3.0) {
+        return (b3[4]);
+      }
+      if (beta == 5.0) {
+        return (b3[5]);
+      }
+      if (beta == 8.0) {
+        return (b3[6]);
+      }
+    }
+
+    if (lambda == 5.0) {
+      if (beta == 0.1) {
+        return (b5[0]);
+      }
+      if (beta == 0.5) {
+        return (b5[1]);
+      }
+      if (beta == 1.0) {
+        return (b5[2]);
+      }
+      if (beta == 2.0) {
+        return (b5[3]);
+      }
+      if (beta == 3.0) {
+        return (b5[4]);
+      }
+      if (beta == 5.0) {
+        return (b5[5]);
+      }
+      if (beta == 8.0) {
+        return (b5[6]);
+      }
+    }
+
+    if (lambda == 8.0) {
+      if (beta == 0.1) {
+        return (b8[0]);
+      }
+      if (beta == 0.5) {
+        return (b8[1]);
+      }
+      if (beta == 1.0) {
+        return (b8[2]);
+      }
+      if (beta == 2.0) {
+        return (b8[3]);
+      }
+      if (beta == 3.0) {
+        return (b8[4]);
+      }
+      if (beta == 5.0) {
+        return (b8[5]);
+      }
+      if (beta == 8.0) {
+        return (b8[6]);
+      }
+    }
+
+
+    int i;
+    double erg;
+    double sum;
+    if ((beta - 5.0 * lambda - 8.0) >= 0.0) {
+      double my = 4.0 * lambda * lambda;
+      double c = -0.9189385 + 0.5 * Math.log(beta) + beta;
+      sum = 1.0;
+      i = 1;
+      double prod = 0.0;
+      double diff = 8.0;
+      double value = 1.0;
+      while (true) { //while (!NULL) {
+        if ((factorial(i) * Math.pow((8.0 * beta), i)) > 1.0e250) {
+          break;
+        }
+        if (i > 10) {
+          break;
+        }
+        if (i == 1) {
+          prod = my - 1.0;
+        } else {
+          value += diff;
+          prod *= (my - value);
+          diff *= 2.0;
+        }
+        sum += prod / (factorial(i) * Math.pow((8.0 * beta), i));
+        i++;
+      }
+      erg = c - Math.log(sum);
+      return (erg);
+    }
+
+    double pi = Math.PI;
+    if ((lambda > 0.0) && ((beta - 0.04 * lambda) <= 0.0)) {
+      if (lambda < 11.5) {
+        erg = -Math.log(gamma(lambda)) - lambda * Math.log(2.0) + lambda * Math.log(beta);
+        return (erg);
+      } else {
+        erg = -(lambda + 1.0) * Math.log(2.0) - (lambda - 0.5) * Math.log(lambda) + lambda + lambda * Math.log(beta) -
+            0.5 * Math.log(0.5 * pi);
+        return (erg);
+      }
+    }
+
+
+    // otherwise numerical integration of the function defined above
+
+    double x = 0.0;
+
+    double new_value;
+    double x1;
+    double step;
+    if (beta < 1.57) {
+      double fx = (fkt2_value(lambda, beta, x)) * 0.01;
+      double y = 0.0;
+      while (true) { //while (!NULL) {
+        y += 0.1;
+        if ((fkt2_value(lambda, beta, y)) < fx) {
+          break;
+        }
+      }
+      step = y * 0.001;
+      x1 = step;
+      sum = (0.5 * (10.0 * step + fkt2_value(lambda, beta, x1))) * step;
+      double first_value = sum;
+      double epsilon = 0.01;
+      while (true) { //while (!NULL) {
+        x = x1;
+        x1 += step;
+        new_value = (0.5 * (fkt2_value(lambda, beta, x) + fkt2_value(lambda, beta, x1))) * step;
+        sum += new_value;
+        if ((new_value / first_value) < epsilon) {
+          break;
+        }
+      }
+      erg = -Math.log(2.0 * sum);
+      return (erg);
+    } else {
+      double z1 = beta / 1.57;
+      sum = 0.0;
+      double period = pi / z1;
+      step = 0.1 * period;
+      double border = 100.0 / ((lambda + 0.1) * (lambda + 0.1));
+      int nr_per = (int) Math.ceil((border / period)) + 20;
+      x1 = step;
+      int j;
+      double z2 = 1.57;
+      for (i = 1; i <= nr_per; i++) {
+        for (j = 1; j <= 10; j++) {
+          new_value = (0.5 * (_fkt_value(lambda, z1, z2, x) + _fkt_value(lambda, z1, z2, x1))) * step;
+          sum += new_value;
+          x = x1;
+          x1 += step;
+        }
+      }
+      for (j = 1; j <= 5; j++) {
+        new_value = (0.5 * (_fkt_value(lambda, z1, z2, x) + _fkt_value(lambda, z1, z2, x1))) * step;
+        sum += new_value;
+        x = x1;
+        x1 += step;
+      }
+      double first_sum = sum;
+      for (j = 1; j <= 10; j++) {
+        new_value = (0.5 * (_fkt_value(lambda, z1, z2, x) + _fkt_value(lambda, z1, z2, x1))) * step;
+        sum += new_value;
+        x = x1;
+        x1 += step;
+      }
+      double second_sum = sum;
+      sum = 0.5 * (first_sum + second_sum);
+      erg = gamma(lambda + 0.5) * Math.pow((2.0 * z2), lambda) / (Math.sqrt(pi) * Math.pow(z1, lambda)) * sum;
+      erg = -Math.log(2.0 * erg);
+      return (erg);
+    }
+  }
+
+  /** Modified Bessel Functions of First Kind - Order 0. */
+  public static double bessi0(double x) {
+    double ax, ans;
+    double y;
+
+    if ((ax = Math.abs(x)) < 3.75) {
+      y = x / 3.75;
+      y *= y;
+      ans = 1.0 + y * (3.5156229 + y * (3.0899424 + y * (1.2067492
+          + y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
+    } else {
+      y = 3.75 / ax;
+      ans = (Math.exp(ax) / Math.sqrt(ax)) * (0.39894228 + y * (0.1328592e-1
+          + y * (0.225319e-2 + y * (-0.157565e-2 + y * (0.916281e-2
+          + y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1
+          + y * 0.392377e-2))))))));
+    }
+    return ans;
+  }
+
+  /** Modified Bessel Functions of First Kind - Order 1. */
+  public static double bessi1(double x) {
+    double ax, ans;
+    double y;
+
+    if ((ax = Math.abs(x)) < 3.75) {
+      y = x / 3.75;
+      y *= y;
+      ans = ax * (0.5 + y * (0.87890594 + y * (0.51498869 + y * (0.15084934
+          + y * (0.2658733e-1 + y * (0.301532e-2 + y * 0.32411e-3))))));
+    } else {
+      y = 3.75 / ax;
+      ans = 0.2282967e-1 + y * (-0.2895312e-1 + y * (0.1787654e-1
+          - y * 0.420059e-2));
+      ans = 0.39894228 + y * (-0.3988024e-1 + y * (-0.362018e-2
+          + y * (0.163801e-2 + y * (-0.1031555e-1 + y * ans))));
+      ans *= (Math.exp(ax) / Math.sqrt(ax));
+    }
+    return x < 0.0 ? -ans : ans;
+  }
+
+  /** Returns <tt>n!</tt>. */
+  public static long factorial(int n) {
+    return Arithmetic.longFactorial(n);
+  }
+
+  private static double fkt2_value(double lambda, double beta, double x_value) {
+
+    double y_value = cosh(lambda * x_value) * Math.exp(-beta * cosh(x_value));
+    return y_value;
+  }
+
+  private static double cosh(double x) {
+    return (Math.exp(x) + Math.exp(-x)) / 2.0;
+  }
+
+
+  /** Returns the gamma function <tt>gamma(x)</tt>. */
+  public static double gamma(double x) {
+    x = logGamma(x);
+    //if (x > Math.log(Double.MAX_VALUE)) return Double.MAX_VALUE;
+    return Math.exp(x);
+  }
+
+  /** Returns a quick approximation of <tt>log(gamma(x))</tt>. */
+  public static double logGamma(double x) {
+
+    if (x <= 0.0 /* || x > 1.3e19 */) {
+      return -999;
+    }
+
+    double z;
+    for (z = 1.0; x < 11.0; x++) {
+      z *= x;
+    }
+
+    double r = 1.0 / (x * x);
+    double c6 = -1.9175269175269175e-03;
+    double c5 = 8.4175084175084175e-04;
+    double c4 = -5.9523809523809524e-04;
+    double c3 = 7.9365079365079365e-04;
+    double c2 = -2.7777777777777777e-03;
+    double c1 = 8.3333333333333333e-02;
+    double g = c1 + r * (c2 + r * (c3 + r * (c4 + r * (c5 + r + c6))));
+    double c0 = 9.1893853320467274e-01;
+    g = (x - 0.5) * Math.log(x) - x + c0 + g / x;
+    if (z == 1.0) {
+      return g;
+    }
+    return g - Math.log(z);
+  }
+}

Propchange: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Fun.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Gamma.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Gamma.java?rev=891983&view=auto
==============================================================================
--- lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Gamma.java (added)
+++ lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Gamma.java Thu Dec 17 23:22:16 2009
@@ -0,0 +1,280 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
+is hereby granted without fee, provided that the above copyright notice appear in all copies and 
+that both that copyright notice and this permission notice appear in supporting documentation. 
+CERN makes no representations about the suitability of this software for any purpose. 
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.random.engine.RandomEngine;
+import org.apache.mahout.math.jet.stat.Probability;
+
+/** @deprecated until unit tests are in place.  Until this time, this class/interface is unsupported. */
+@Deprecated
+public class Gamma extends AbstractContinousDistribution {
+
+  private double alpha;
+  private double lambda;
+
+  // The uniform random number generated shared by all <b>static</b> methods.
+  private static final Gamma shared = new Gamma(1.0, 1.0, makeDefaultGenerator());
+
+  /**
+   * Constructs a Gamma distribution. Example: alpha=1.0, lambda=1.0.
+   *
+   * @throws IllegalArgumentException if <tt>alpha &lt;= 0.0 || lambda &lt;= 0.0</tt>.
+   */
+  public Gamma(double alpha, double lambda, RandomEngine randomGenerator) {
+    setRandomGenerator(randomGenerator);
+    setState(alpha, lambda);
+  }
+
+  /** Returns the cumulative distribution function. */
+  public double cdf(double x) {
+    return Probability.gamma(alpha, lambda, x);
+  }
+
+  /** Returns a random number from the distribution. */
+  @Override
+  public double nextDouble() {
+    return nextDouble(alpha, lambda);
+  }
+
+  /** Returns a random number from the distribution; bypasses the internal state. */
+  public double nextDouble(double alpha, double lambda) {
+/******************************************************************
+ *                                                                *
+ *    Gamma Distribution - Acceptance Rejection combined with     *
+ *                         Acceptance Complement                  *
+ *                                                                *
+ ******************************************************************
+ *                                                                *
+ * FUNCTION:    - gds samples a random number from the standard   *
+ *                gamma distribution with parameter  a > 0.       *
+ *                Acceptance Rejection  gs  for  a < 1 ,          *
+ *                Acceptance Complement gd  for  a >= 1 .         *
+ * REFERENCES:  - J.H. Ahrens, U. Dieter (1974): Computer methods *
+ *                for sampling from gamma, beta, Poisson and      *
+ *                binomial distributions, Computing 12, 223-246.  *
+ *              - J.H. Ahrens, U. Dieter (1982): Generating gamma *
+ *                variates by a modified rejection technique,     *
+ *                Communications of the ACM 25, 47-54.            *
+ * SUBPROGRAMS: - drand(seed) ... (0,1)-Uniform generator with    *
+ *                unsigned long integer *seed                     *
+ *              - NORMAL(seed) ... Normal generator N(0,1).       *
+ *                                                                *
+ ******************************************************************/
+    double a = alpha;
+
+    // Check for invalid input values
+
+    if (a <= 0.0) {
+      throw new IllegalArgumentException();
+    }
+    if (lambda <= 0.0) {
+      throw new IllegalArgumentException();
+    }
+
+    double gds;
+    double b = 0.0;
+    if (a < 1.0) { // CASE A: Acceptance rejection algorithm gs
+      b = 1.0 + 0.36788794412 * a;              // Step 1
+      while (true) {
+        double p = b * randomGenerator.raw();
+        if (p <= 1.0) {                       // Step 2. Case gds <= 1
+          gds = Math.exp(Math.log(p) / a);
+          if (Math.log(randomGenerator.raw()) <= -gds) {
+            return (gds / lambda);
+          }
+        } else {                                // Step 3. Case gds > 1
+          gds = -Math.log((b - p) / a);
+          if (Math.log(randomGenerator.raw()) <= ((a - 1.0) * Math.log(gds))) {
+            return (gds / lambda);
+          }
+        }
+      }
+    } else {        // CASE B: Acceptance complement algorithm gd (gaussian distribution, box muller transformation)
+      double ss = 0.0;
+      double s = 0.0;
+      double d = 0.0;
+      double aa = -1.0;
+      if (a != aa) {                        // Step 1. Preparations
+        aa = a;
+        ss = a - 0.5;
+        s = Math.sqrt(ss);
+        d = 5.656854249 - 12.0 * s;
+      }
+      // Step 2. Normal deviate
+      double v12;
+      double v1;
+      do {
+        v1 = 2.0 * randomGenerator.raw() - 1.0;
+        double v2 = 2.0 * randomGenerator.raw() - 1.0;
+        v12 = v1 * v1 + v2 * v2;
+      } while (v12 > 1.0);
+      double t = v1 * Math.sqrt(-2.0 * Math.log(v12) / v12);
+      double x = s + 0.5 * t;
+      gds = x * x;
+      if (t >= 0.0) {
+        return (gds / lambda);
+      }         // Immediate acceptance
+
+      double u = randomGenerator.raw();
+      if (d * u <= t * t * t) {
+        return (gds / lambda);
+      } // Squeeze acceptance
+
+      double q0 = 0.0;
+      double si = 0.0;
+      double c = 0.0;
+      double aaa = -1.0;
+      if (a != aaa) {                           // Step 4. Set-up for hat case
+        aaa = a;
+        double r = 1.0 / a;
+        double q9 = 0.0001710320;
+        double q8 = -0.0004701849;
+        double q7 = 0.0006053049;
+        double q6 = 0.0003340332;
+        double q5 = -0.0003349403;
+        double q4 = 0.0015746717;
+        double q3 = 0.0079849875;
+        double q2 = 0.0208333723;
+        double q1 = 0.0416666664;
+        q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) *
+            r + q3) * r + q2) * r + q1) * r;
+        if (a > 3.686) {
+          if (a > 13.022) {
+            b = 1.77;
+            si = 0.75;
+            c = 0.1515 / s;
+          } else {
+            b = 1.654 + 0.0076 * ss;
+            si = 1.68 / s + 0.275;
+            c = 0.062 / s + 0.024;
+          }
+        } else {
+          b = 0.463 + s - 0.178 * ss;
+          si = 1.235;
+          c = 0.195 / s - 0.079 + 0.016 * s;
+        }
+      }
+      double v;
+      double q;
+      double a9 = 0.104089866;
+      double a8 = -0.112750886;
+      double a7 = 0.110368310;
+      double a6 = -0.124385581;
+      double a5 = 0.142873973;
+      double a4 = -0.166677482;
+      double a3 = 0.199999867;
+      double a2 = -0.249999949;
+      double a1 = 0.333333333;
+      if (x > 0.0) {                        // Step 5. Calculation of q
+        v = t / (s + s);                  // Step 6.
+        if (Math.abs(v) > 0.25) {
+          q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
+        } else {
+          q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
+              v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
+        }                  // Step 7. Quotient acceptance
+        if (Math.log(1.0 - u) <= q) {
+          return (gds / lambda);
+        }
+      }
+
+      double e7 = 0.000247453;
+      double e6 = 0.001353826;
+      double e5 = 0.008345522;
+      double e4 = 0.041664508;
+      double e3 = 0.166666848;
+      double e2 = 0.499999994;
+      double e1 = 1.000000000;
+      while (true) {                          // Step 8. Double exponential deviate t
+        double sign_u;
+        double e;
+        do {
+          e = -Math.log(randomGenerator.raw());
+          u = randomGenerator.raw();
+          u = u + u - 1.0;
+          sign_u = (u > 0) ? 1.0 : -1.0;
+          t = b + (e * si) * sign_u;
+        } while (t <= -0.71874483771719); // Step 9. Rejection of t
+        v = t / (s + s);                  // Step 10. New q(t)
+        if (Math.abs(v) > 0.25) {
+          q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
+        } else {
+          q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
+              v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
+        }
+        if (q <= 0.0) {
+          continue;
+        }           // Step 11.
+        double w;
+        if (q > 0.5) {
+          w = Math.exp(q) - 1.0;
+        } else {
+          w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) *
+              q + e1) * q;
+        }                            // Step 12. Hat acceptance
+        if (c * u * sign_u <= w * Math.exp(e - 0.5 * t * t)) {
+          x = s + 0.5 * t;
+          return (x * x / lambda);
+        }
+      }
+    }
+  }
+
+  /** Returns the probability distribution function. */
+  public double pdf(double x) {
+    if (x < 0) {
+      throw new IllegalArgumentException();
+    }
+    if (x == 0) {
+      if (alpha == 1.0) {
+        return 1.0 / lambda;
+      } else {
+        return 0.0;
+      }
+    }
+    if (alpha == 1.0) {
+      return Math.exp(-x / lambda) / lambda;
+    }
+
+    return Math.exp((alpha - 1.0) * Math.log(x / lambda) - x / lambda - Fun.logGamma(alpha)) / lambda;
+  }
+
+  /**
+   * Sets the mean and variance.
+   *
+   * @throws IllegalArgumentException if <tt>alpha &lt;= 0.0 || lambda &lt;= 0.0</tt>.
+   */
+  public void setState(double alpha, double lambda) {
+    if (alpha <= 0.0) {
+      throw new IllegalArgumentException();
+    }
+    if (lambda <= 0.0) {
+      throw new IllegalArgumentException();
+    }
+    this.alpha = alpha;
+    this.lambda = lambda;
+  }
+
+  /**
+   * Returns a random number from the distribution.
+   *
+   * @throws IllegalArgumentException if <tt>alpha &lt;= 0.0 || lambda &lt;= 0.0</tt>.
+   */
+  public static double staticNextDouble(double alpha, double lambda) {
+    synchronized (shared) {
+      return shared.nextDouble(alpha, lambda);
+    }
+  }
+
+  /** Returns a String representation of the receiver. */
+  public String toString() {
+    return this.getClass().getName() + '(' + alpha + ',' + lambda + ')';
+  }
+
+}

Propchange: lucene/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Gamma.java
------------------------------------------------------------------------------
    svn:eol-style = native



Mime
View raw message