You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@mahout.apache.org by td...@apache.org on 2010/07/22 00:37:52 UTC
svn commit: r966447 - in /mahout/trunk/math/src:
main/java/org/apache/mahout/math/jet/math/
main/java/org/apache/mahout/math/jet/random/
main/java/org/apache/mahout/math/jet/stat/
test/java/org/apache/mahout/math/jet/random/
Author: tdunning
Date: Wed Jul 21 22:37:51 2010
New Revision: 966447
URL: http://svn.apache.org/viewvc?rev=966447&view=rev
Log:
MAHOUT-443 - added tests for Gamma distribution and changed API to be more standard.
Along the way, I also checked Polynomial, parts of stat.Gamma and stat.Probability.
Added:
mahout/trunk/math/src/test/java/org/apache/mahout/math/jet/random/GammaTest.java
Modified:
mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/math/Polynomial.java
mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/AbstractDistribution.java
mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Gamma.java
mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/NegativeBinomial.java
mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/stat/Gamma.java
mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/stat/Probability.java
Modified: mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/math/Polynomial.java
URL: http://svn.apache.org/viewvc/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/math/Polynomial.java?rev=966447&r1=966446&r2=966447&view=diff
==============================================================================
--- mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/math/Polynomial.java (original)
+++ mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/math/Polynomial.java Wed Jul 21 22:37:51 2010
@@ -1,4 +1,21 @@
/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*
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
@@ -12,8 +29,6 @@ package org.apache.mahout.math.jet.math;
* Polynomial functions.
*/
-/** @deprecated until unit tests are in place. Until this time, this class/interface is unsupported. */
-@Deprecated
public class Polynomial extends Constants {
/** Makes this class non instantiable, but still let's others inherit from it. */
Modified: mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/AbstractDistribution.java
URL: http://svn.apache.org/viewvc/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/AbstractDistribution.java?rev=966447&r1=966446&r2=966447&view=diff
==============================================================================
--- mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/AbstractDistribution.java (original)
+++ mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/AbstractDistribution.java Wed Jul 21 22:37:51 2010
@@ -1,21 +1,20 @@
-/**
- * Licensed to the Apache Software Foundation (ASF) under one
- * or more contributor license agreements. See the NOTICE file
- * distributed with this work for additional information
- * regarding copyright ownership. The ASF licenses this file
- * to you under the Apache License, Version 2.0 (the
- * "License"); you may not use this file except in compliance
- * with the License. You may obtain a copy of the License at
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
*
- * http://www.apache.org/licenses/LICENSE-2.0
+ * http://www.apache.org/licenses/LICENSE-2.0
*
- * Unless required by applicable law or agreed to in writing,
- * software distributed under the License is distributed on an
- * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
- * KIND, either express or implied. See the License for the
- * specific language governing permissions and limitations
- * under the License.
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
*/
+
/*
Copyright � 1999 CERN - European Organization for Nuclear Research.
Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
@@ -31,8 +30,7 @@ import org.apache.mahout.math.function.U
import org.apache.mahout.math.function.IntFunction;
import org.apache.mahout.math.jet.random.engine.RandomEngine;
-public abstract class AbstractDistribution extends PersistentObject
- implements UnaryFunction, IntFunction {
+public abstract class AbstractDistribution extends PersistentObject implements UnaryFunction, IntFunction {
protected RandomEngine randomGenerator;
Modified: mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Gamma.java
URL: http://svn.apache.org/viewvc/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Gamma.java?rev=966447&r1=966446&r2=966447&view=diff
==============================================================================
--- mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Gamma.java (original)
+++ mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/Gamma.java Wed Jul 21 22:37:51 2010
@@ -1,4 +1,21 @@
/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*
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
@@ -11,68 +28,71 @@ package org.apache.mahout.math.jet.rando
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 {
-
+ // shape
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());
+ // rate
+ private double beta;
/**
- * Constructs a Gamma distribution. Example: alpha=1.0, lambda=1.0.
+ * Constructs a Gamma distribution with a given shape (alpha) and rate (beta).
*
- * @throws IllegalArgumentException if <tt>alpha <= 0.0 || lambda <= 0.0</tt>.
+ * @param alpha The shape parameter.
+ * @param beta The rate parameter.
+ * @param randomGenerator The random number generator that generates bits for us.
+ * @throws IllegalArgumentException if <tt>alpha <= 0.0 || alpha <= 0.0</tt>.
*/
- public Gamma(double alpha, double lambda, RandomEngine randomGenerator) {
+ public Gamma(double alpha, double beta, RandomEngine randomGenerator) {
+ this.alpha = alpha;
+ this.beta = beta;
setRandomGenerator(randomGenerator);
- setState(alpha, lambda);
}
- /** Returns the cumulative distribution function. */
+ /**
+ * Returns the cumulative distribution function.
+ * @param x The end-point where the cumulation should end.
+ */
public double cdf(double x) {
- return Probability.gamma(alpha, lambda, x);
+ return Probability.gamma(alpha, beta, x);
}
/** Returns a random number from the distribution. */
@Override
public double nextDouble() {
- return nextDouble(alpha, lambda);
+ return nextDouble(alpha, beta);
}
- /** 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). *
- * *
- ******************************************************************/
-
- // Check for invalid input values
-
+ /** Returns a random number from the distribution; bypasses the internal state.
+ * *
+ * 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). *
+ * *
+ * @param beta Scale parameter.
+ * @param alpha Shape parameter.
+ * @return A gamma distributed sample.
+ */
+ public double nextDouble(double alpha, double beta) {
if (alpha <= 0.0) {
throw new IllegalArgumentException();
}
- if (lambda <= 0.0) {
+ if (beta <= 0.0) {
throw new IllegalArgumentException();
}
@@ -85,12 +105,12 @@ public class Gamma extends AbstractConti
if (p <= 1.0) { // Step 2. Case gds <= 1
gds = Math.exp(Math.log(p) / alpha);
if (Math.log(randomGenerator.raw()) <= -gds) {
- return (gds / lambda);
+ return (gds / beta);
}
} else { // Step 3. Case gds > 1
gds = -Math.log((b - p) / alpha);
if (Math.log(randomGenerator.raw()) <= ((alpha - 1.0) * Math.log(gds))) {
- return (gds / lambda);
+ return (gds / beta);
}
}
}
@@ -117,12 +137,12 @@ public class Gamma extends AbstractConti
double x = s + 0.5 * t;
gds = x * x;
if (t >= 0.0) {
- return (gds / lambda);
+ return (gds / beta);
} // Immediate acceptance
double u = randomGenerator.raw();
if (d * u <= t * t * t) {
- return (gds / lambda);
+ return (gds / beta);
} // Squeeze acceptance
double q0 = 0.0;
@@ -179,7 +199,7 @@ public class Gamma extends AbstractConti
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);
+ return (gds / beta);
}
}
@@ -219,61 +239,37 @@ public class Gamma extends AbstractConti
} // 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);
+ return (x * x / beta);
}
}
}
}
- /** Returns the probability distribution function. */
+ /** Returns the probability distribution function.
+ * @param x Where to compute the density function.
+ *
+ * @return The value of the gamma density at x.
+ */
public double pdf(double x) {
if (x < 0) {
throw new IllegalArgumentException();
}
if (x == 0) {
if (alpha == 1.0) {
- return 1.0 / lambda;
+ return beta;
+ } else if (alpha < 1) {
+ return Double.POSITIVE_INFINITY;
} else {
- return 0.0;
+ return 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 <= 0.0 || lambda <= 0.0</tt>.
- */
- public void setState(double alpha, double lambda) {
- if (alpha <= 0.0) {
- throw new IllegalArgumentException();
- }
- if (lambda <= 0.0) {
- throw new IllegalArgumentException();
+ return beta * Math.exp(-x * beta);
}
- this.alpha = alpha;
- this.lambda = lambda;
+ return beta * Math.exp((alpha - 1.0) * Math.log(x * beta) - x * beta - Fun.logGamma(alpha));
}
- /**
- * Returns a random number from the distribution.
- *
- * @throws IllegalArgumentException if <tt>alpha <= 0.0 || lambda <= 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 + ')';
+ return this.getClass().getName() + '(' + beta + ',' + alpha + ')';
}
-
}
Modified: mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/NegativeBinomial.java
URL: http://svn.apache.org/viewvc/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/NegativeBinomial.java?rev=966447&r1=966446&r2=966447&view=diff
==============================================================================
--- mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/NegativeBinomial.java (original)
+++ mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/random/NegativeBinomial.java Wed Jul 21 22:37:51 2010
@@ -1,4 +1,21 @@
/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*
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
@@ -34,7 +51,7 @@ public class NegativeBinomial extends Ab
public NegativeBinomial(int n, double p, RandomEngine randomGenerator) {
setRandomGenerator(randomGenerator);
setNandP(n, p);
- this.gamma = new Gamma(n, 1.0, randomGenerator);
+ this.gamma = new Gamma(n, 1, randomGenerator);
this.poisson = new Poisson(0.0, randomGenerator);
}
@@ -97,7 +114,7 @@ public class NegativeBinomial extends Ab
double x = p / (1.0 - p);
//double p1 = p;
- double y = x * this.gamma.nextDouble(n, 1.0);
+ double y = x * this.gamma.nextDouble(n, 1);
return this.poisson.nextInt(y);
}
Modified: mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/stat/Gamma.java
URL: http://svn.apache.org/viewvc/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/stat/Gamma.java?rev=966447&r1=966446&r2=966447&view=diff
==============================================================================
--- mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/stat/Gamma.java (original)
+++ mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/stat/Gamma.java Wed Jul 21 22:37:51 2010
@@ -1,4 +1,21 @@
/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*
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
@@ -26,22 +43,25 @@ public class Gamma extends org.apache.ma
* -
* | (a+b)
* </pre>
+ * @param alpha
+ * @param beta
+ * @return The beta function for given values of alpha and beta.
*/
@Deprecated
- public static double beta(double a, double b) throws ArithmeticException {
+ public static double beta(double alpha, double beta) {
- double y = a + b;
+ double y = alpha + beta;
y = gamma(y);
if (y == 0.0) {
return 1.0;
}
- if (a > b) {
- y = gamma(a) / y;
- y *= gamma(b);
+ if (alpha > beta) {
+ y = gamma(alpha) / y;
+ y *= gamma(beta);
} else {
- y = gamma(b) / y;
- y *= gamma(a);
+ y = gamma(beta) / y;
+ y *= gamma(alpha);
}
return (y);
@@ -411,24 +431,21 @@ public class Gamma extends org.apache.ma
/**
* Returns the Incomplete Gamma function; formerly named <tt>igamma</tt>.
*
- * @param a the parameter of the gamma distribution.
+ * @param alpha the shape parameter of the gamma distribution.
* @param x the integration end point.
+ * @return The value of the unnormalized incomplete gamma function.
*/
- @Deprecated
- public static double incompleteGamma(double a, double x)
- throws ArithmeticException {
-
-
- if (x <= 0 || a <= 0) {
+ public static double incompleteGamma(double alpha, double x){
+ if (x <= 0 || alpha <= 0) {
return 0.0;
}
- if (x > 1.0 && x > a) {
- return 1.0 - incompleteGammaComplement(a, x);
+ if (x > 1.0 && x > alpha) {
+ return 1.0 - incompleteGammaComplement(alpha, x);
}
/* Compute x**a * exp(-x) / gamma(a) */
- double ax = a * Math.log(x) - x - logGamma(a);
+ double ax = alpha * Math.log(x) - x - logGamma(alpha);
if (ax < -MAXLOG) {
return (0.0);
}
@@ -436,7 +453,7 @@ public class Gamma extends org.apache.ma
ax = Math.exp(ax);
/* power series */
- double r = a;
+ double r = alpha;
double c = 1.0;
double ans = 1.0;
@@ -447,28 +464,27 @@ public class Gamma extends org.apache.ma
}
while (c / ans > MACHEP);
- return (ans * ax / a);
+ return (ans * ax / alpha);
}
/**
* Returns the Complemented Incomplete Gamma function; formerly named <tt>igamc</tt>.
*
- * @param a the parameter of the gamma distribution.
+ * @param alpha the shape parameter of the gamma distribution.
* @param x the integration start point.
*/
- @Deprecated
- public static double incompleteGammaComplement(double a, double x) throws ArithmeticException {
+ public static double incompleteGammaComplement(double alpha, double x) {
- if (x <= 0 || a <= 0) {
+ if (x <= 0 || alpha <= 0) {
return 1.0;
}
- if (x < 1.0 || x < a) {
- return 1.0 - incompleteGamma(a, x);
+ if (x < 1.0 || x < alpha) {
+ return 1.0 - incompleteGamma(alpha, x);
}
- double ax = a * Math.log(x) - x - logGamma(a);
+ double ax = alpha * Math.log(x) - x - logGamma(alpha);
if (ax < -MAXLOG) {
return 0.0;
}
@@ -476,7 +492,7 @@ public class Gamma extends org.apache.ma
ax = Math.exp(ax);
/* continued fraction */
- double y = 1.0 - a;
+ double y = 1.0 - alpha;
double z = x + y + 1.0;
double c = 0.0;
double pkm2 = 1.0;
@@ -517,7 +533,7 @@ public class Gamma extends org.apache.ma
}
/** Returns the natural logarithm of the gamma function; formerly named <tt>lgamma</tt>. */
- public static double logGamma(double x) throws ArithmeticException {
+ public static double logGamma(double x) {
double p, q, z;
double[] A = {
Modified: mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/stat/Probability.java
URL: http://svn.apache.org/viewvc/mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/stat/Probability.java?rev=966447&r1=966446&r2=966447&view=diff
==============================================================================
--- mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/stat/Probability.java (original)
+++ mahout/trunk/math/src/main/java/org/apache/mahout/math/jet/stat/Probability.java Wed Jul 21 22:37:51 2010
@@ -1,4 +1,21 @@
/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*
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
@@ -417,54 +434,53 @@ public class Probability extends Constan
/**
* Returns the integral from zero to <tt>x</tt> of the gamma probability density function.
* <pre>
- * x
- * b -
- * a | | b-1 -at
- * y = ----- | t e dt
- * - | |
- * | (b) -
- * 0
+ *
+ * alpha - x
+ * beta | alpha-1 -beta t
+ * y = --------- | t e dt
+ * - |
+ * | (alpha) - 0
* </pre>
* The incomplete gamma integral is used, according to the relation
*
- * <tt>y = Gamma.incompleteGamma( b, a*x )</tt>.
+ * <tt>y = Gamma.incompleteGamma( alpha, beta*x )</tt>.
+ *
+ * See http://en.wikipedia.org/wiki/Gamma_distribution#Probability_density_function
*
- * @param a the paramater a (alpha) of the gamma distribution.
- * @param b the paramater b (beta, lambda) of the gamma distribution.
+ * @param alpha the shape parameter of the gamma distribution.
+ * @param beta the rate parameter of the gamma distribution.
* @param x integration end point.
*/
- @Deprecated
- public static double gamma(double a, double b, double x) {
+ public static double gamma(double alpha, double beta, double x) {
if (x < 0.0) {
return 0.0;
}
- return Gamma.incompleteGamma(b, a * x);
+ return Gamma.incompleteGamma(alpha, beta * x);
}
/**
* Returns the integral from <tt>x</tt> to infinity of the gamma probability density function:
* <pre>
- * inf.
- * b -
- * a | | b-1 -at
- * y = ----- | t e dt
- * - | |
- * | (b) -
- * x
+ * alpha - infinity
+ * beta | alpha-1 -beta t
+ * y = --------- | t e dt
+ * - |
+ * | (alpha) - x
* </pre>
* The incomplete gamma integral is used, according to the relation <p> y = Gamma.incompleteGammaComplement( b, a*x
* ).
*
- * @param a the paramater a (alpha) of the gamma distribution.
- * @param b the paramater b (beta, lambda) of the gamma distribution.
+ * TODO this method is inconsistent with gamma(alpha, beta, x)
+ *
+ * @param alpha the shape parameter of the gamma distribution.
+ * @param beta the rate parameter of the gamma distribution.
* @param x integration end point.
*/
- @Deprecated
- public static double gammaComplemented(double a, double b, double x) {
+ public static double gammaComplemented(double alpha, double beta, double x) {
if (x < 0.0) {
return 0.0;
}
- return Gamma.incompleteGammaComplement(b, a * x);
+ return Gamma.incompleteGammaComplement(alpha, beta * x);
}
/**
Added: mahout/trunk/math/src/test/java/org/apache/mahout/math/jet/random/GammaTest.java
URL: http://svn.apache.org/viewvc/mahout/trunk/math/src/test/java/org/apache/mahout/math/jet/random/GammaTest.java?rev=966447&view=auto
==============================================================================
--- mahout/trunk/math/src/test/java/org/apache/mahout/math/jet/random/GammaTest.java (added)
+++ mahout/trunk/math/src/test/java/org/apache/mahout/math/jet/random/GammaTest.java Wed Jul 21 22:37:51 2010
@@ -0,0 +1,127 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.mahout.math.jet.random;
+
+import org.apache.mahout.math.jet.random.engine.MersenneTwister;
+import org.apache.mahout.math.jet.random.engine.RandomEngine;
+import org.junit.Assert;
+import org.junit.Test;
+
+import java.util.Arrays;
+
+import static org.junit.Assert.assertEquals;
+
+public class GammaTest {
+ @Test
+ public void testNextDouble() {
+ double[] z = new double[100000];
+ RandomEngine gen = new MersenneTwister(1);
+ for (double alpha : new double[]{1, 2, 10, 0.1, 0.01, 100}) {
+ Gamma g = new Gamma(alpha, 1, gen);
+ for (int i = 0; i < z.length; i++) {
+ z[i] = g.nextDouble();
+ }
+ Arrays.sort(z);
+
+ // verify that empirical CDF matches theoretical one pretty closely
+ for (double q : seq(0.01, 1, 0.01)) {
+ double p = z[(int) (q * z.length)];
+ assertEquals(q, g.cdf(p), 0.01);
+ }
+ }
+ }
+
+ @Test
+ public void testCdf() {
+ RandomEngine gen = new MersenneTwister(1);
+
+ // verify scaling for special case of alpha = 1
+ for (double beta : new double[]{1, 0.1, 2, 100}) {
+ Gamma g1 = new Gamma(1, beta, gen);
+ Gamma g2 = new Gamma(1, 1, gen);
+ for (double x : seq(0, 0.99, 0.1)) {
+ assertEquals(String.format("Rate invariance: x = %.4f, alpha = 1, beta = %.1f", x, beta),
+ 1 - Math.exp(-x * beta), g1.cdf(x), 1e-9);
+ assertEquals(String.format("Rate invariance: x = %.4f, alpha = 1, beta = %.1f", x, beta),
+ g2.cdf(beta * x), g1.cdf(x), 1e-9);
+ }
+ }
+
+ // now test scaling for a selection of values of alpha
+ for (double alpha : new double[]{0.01, 0.1, 1, 2, 10, 100, 1000}) {
+ Gamma g = new Gamma(alpha, 1, gen);
+ for (double beta : new double[]{0.1, 1, 2, 100}) {
+ Gamma g1 = new Gamma(alpha, beta, gen);
+ for (double x : seq(0, 0.9999, 0.001)) {
+ Assert.assertEquals(String.format("Rate invariance: x = %.4f, alpha = %.2f, beta = %.1f", x, alpha, beta),
+ g.cdf(x * beta), g1.cdf(x), 0);
+ }
+ }
+ }
+
+ // now check against known values computed using R for various values of alpha
+ checkGammaCdf(0.01, 1, 0.0000000, 0.9450896, 0.9516444, 0.9554919, 0.9582258, 0.9603474, 0.9620810, 0.9635462, 0.9648148, 0.9659329, 0.9669321);
+ checkGammaCdf(0.1, 1, 0.0000000, 0.7095387, 0.7591012, 0.7891072, 0.8107067, 0.8275518, 0.8413180, 0.8529198, 0.8629131, 0.8716623, 0.8794196);
+ checkGammaCdf(1, 1, 0.0000000, 0.1812692, 0.3296800, 0.4511884, 0.5506710, 0.6321206, 0.6988058, 0.7534030, 0.7981035, 0.8347011, 0.8646647);
+ checkGammaCdf(10, 1, 0.000000e+00, 4.649808e-05, 8.132243e-03, 8.392402e-02, 2.833757e-01, 5.420703e-01, 7.576078e-01, 8.906006e-01, 9.567017e-01, 9.846189e-01, 9.950046e-01);
+ checkGammaCdf(100, 1, 0.000000e+00, 3.488879e-37, 1.206254e-15, 1.481528e-06, 1.710831e-02, 5.132988e-01, 9.721363e-01, 9.998389e-01, 9.999999e-01, 1.000000e+00, 1.000000e+00);
+
+// > pgamma(seq(0,0.02,by=0.002),0.01,1)
+// [1] 0.0000000 0.9450896 0.9516444 0.9554919 0.9582258 0.9603474 0.9620810 0.9635462 0.9648148 0.9659329 0.9669321
+// > pgamma(seq(0,0.2,by=0.02),0.1,1)
+// [1] 0.0000000 0.7095387 0.7591012 0.7891072 0.8107067 0.8275518 0.8413180 0.8529198 0.8629131 0.8716623 0.8794196
+// > pgamma(seq(0,2,by=0.2),1,1)
+// [1] 0.0000000 0.1812692 0.3296800 0.4511884 0.5506710 0.6321206 0.6988058 0.7534030 0.7981035 0.8347011 0.8646647
+// > pgamma(seq(0,20,by=2),10,1)
+// [1] 0.000000e+00 4.649808e-05 8.132243e-03 8.392402e-02 2.833757e-01 5.420703e-01 7.576078e-01 8.906006e-01 9.567017e-01 9.846189e-01 9.950046e-01
+// > pgamma(seq(0,200,by=20),100,1)
+// [1] 0.000000e+00 3.488879e-37 1.206254e-15 1.481528e-06 1.710831e-02 5.132988e-01 9.721363e-01 9.998389e-01 9.999999e-01 1.000000e+00 1.000000e+00
+ }
+
+ private void checkGammaCdf(double alpha, double beta, double... values) {
+ Gamma g = new Gamma(alpha, beta, new MersenneTwister());
+ int i = 0;
+ for (double x : seq(0, 2 * alpha, 2 * alpha / 10)) {
+ assertEquals(String.format("alpha=%.2f, i=%d, x=%.2f", alpha, i, x), values[i], g.cdf(x), 1e-7);
+ i++;
+ }
+ }
+
+ private double[] seq(double from, double to, double by) {
+ double[] r = new double[(int) Math.ceil(0.999999 * (to - from) / by)];
+ int i = 0;
+ for (double x = from; x < to - (to - from) * 1e-6; x += by) {
+ r[i++] = x;
+ }
+ return r;
+ }
+
+ @Test
+ public void testPdf() {
+ RandomEngine gen = new MersenneTwister(1);
+ for (double alpha : new double[]{0.01, 0.1, 1, 2, 10, 100}) {
+ for (double beta : new double[]{0.1, 1, 2, 100}) {
+ Gamma g1 = new Gamma(alpha, beta, gen);
+ for (double x : seq(0, 0.99, 0.1)) {
+ double p = Math.pow(beta, alpha) * Math.pow(x, alpha - 1) * Math.exp(-beta * x - org.apache.mahout.math.jet.stat.Gamma.logGamma(alpha));
+ assertEquals(String.format("alpha=%.2f, beta=%.2f, x=%.2f\n", alpha, beta, x), p, g1.pdf(x), 1e-9);
+ }
+ }
+ }
+ }
+}