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/25 01:54:27 UTC
svn commit: r978950 - in /mahout/trunk: conf/ core/ core/doc/
core/src/main/java/org/apache/mahout/classifier/
core/src/main/java/org/apache/mahout/classifier/evaluation/
core/src/main/java/org/apache/mahout/classifier/sgd/
core/src/test/java/org/apach...
Author: tdunning
Date: Sat Jul 24 23:54:26 2010
New Revision: 978950
URL: http://svn.apache.org/viewvc?rev=978950&view=rev
Log:
Reverting mouse slip checkin
Removed:
mahout/trunk/conf/cat.props
mahout/trunk/conf/runlogistic.props
mahout/trunk/conf/trainlogistic.props
mahout/trunk/core/doc/
mahout/trunk/core/src/main/java/org/apache/mahout/classifier/MurmurHash.java
mahout/trunk/core/src/main/java/org/apache/mahout/classifier/evaluation/
mahout/trunk/core/src/main/java/org/apache/mahout/classifier/sgd/
mahout/trunk/core/src/test/java/org/apache/mahout/classifier/MurmurHashTest.java
mahout/trunk/core/src/test/java/org/apache/mahout/classifier/evaluation/
mahout/trunk/core/src/test/java/org/apache/mahout/classifier/sgd/
mahout/trunk/core/src/test/resources/logP.csv
mahout/trunk/core/src/test/resources/r.csv
mahout/trunk/core/src/test/resources/sgd.csv
mahout/trunk/examples/src/main/java/org/apache/mahout/classifier/sgd/
mahout/trunk/examples/src/main/resources/donut-test.csv
mahout/trunk/examples/src/main/resources/donut.csv
mahout/trunk/examples/src/main/resources/test-data.csv
mahout/trunk/examples/src/test/java/org/apache/mahout/classifier/
mahout/trunk/math/src/test/java/org/apache/mahout/math/jet/stat/GammaTest.java
Modified:
mahout/trunk/conf/driver.classes.props
mahout/trunk/core/pom.xml
mahout/trunk/core/src/test/java/org/apache/mahout/cf/taste/impl/common/FastIDSetTest.java (props changed)
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
Modified: mahout/trunk/conf/driver.classes.props
URL: http://svn.apache.org/viewvc/mahout/trunk/conf/driver.classes.props?rev=978950&r1=978949&r2=978950&view=diff
==============================================================================
--- mahout/trunk/conf/driver.classes.props (original)
+++ mahout/trunk/conf/driver.classes.props Sat Jul 24 23:54:26 2010
@@ -23,6 +23,3 @@ org.apache.mahout.math.hadoop.decomposer
org.apache.mahout.math.hadoop.decomposer.EigenVerificationJob = cleansvd : Cleanup and verification of SVD output
org.apache.mahout.math.hadoop.similarity.RowSimilarityJob = rowsimilarity : Compute the pairwise similarities of the rows of a matrix
org.apache.mahout.cf.taste.hadoop.similarity.item.ItemSimilarityJob = itemsimilarity : Compute the item-item-similarities for item-based collaborative filtering
-org.apache.mahout.classifier.sgd.TrainLogistic = trainlogistic : Train a logistic regression using stochastic gradient descentorg.apache.mahout.classifier.sgd.TrainLogistic = trainlogistic : Train a logistic regression using stochastic gradient descent
-org.apache.mahout.classifier.sgd.RunLogistic = runlogistic : Run a logistic regression model against CSV data
-org.apache.mahout.classifier.sgd.PrintResourceOrFile = cat : Print a file or resource as the logistic regression models would see it
\ No newline at end of file
Modified: mahout/trunk/core/pom.xml
URL: http://svn.apache.org/viewvc/mahout/trunk/core/pom.xml?rev=978950&r1=978949&r2=978950&view=diff
==============================================================================
--- mahout/trunk/core/pom.xml (original)
+++ mahout/trunk/core/pom.xml Sat Jul 24 23:54:26 2010
@@ -1,21 +1,21 @@
<?xml version="1.0" encoding="UTF-8"?>
<!--
- ~ 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.
- -->
+ 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.
+-->
<project xmlns="http://maven.apache.org/POM/4.0.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/maven-v4_0_0.xsd">
<modelVersion>4.0.0</modelVersion>
@@ -268,11 +268,7 @@
<type>test-jar</type>
<scope>test</scope>
</dependency>
- <dependency>
- <groupId>com.google.guava</groupId>
- <artifactId>guava</artifactId>
- <version>r03</version>
- </dependency>
+
</dependencies>
<repositories>
Propchange: mahout/trunk/core/src/test/java/org/apache/mahout/cf/taste/impl/common/FastIDSetTest.java
('svn:mergeinfo' removed)
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=978950&r1=978949&r2=978950&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 Sat Jul 24 23:54:26 2010
@@ -50,8 +50,7 @@ public class NegativeBinomial extends Ab
*/
public NegativeBinomial(int n, double p, RandomEngine randomGenerator) {
setRandomGenerator(randomGenerator);
- this.n = n;
- this.p = p;
+ setNandP(n, p);
this.gamma = new Gamma(n, 1, randomGenerator);
this.poisson = new Poisson(0.0, randomGenerator);
}
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=978950&r1=978949&r2=978950&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 Sat Jul 24 23:54:26 2010
@@ -34,14 +34,42 @@ public class Gamma extends org.apache.ma
protected Gamma() {
}
- /** Returns the Gamma function of the argument. This can be dangerous since a
- * very moderate value of x will lead to overflow. Usually it is better to use
- * logGamma instead. One interesting exception is for negative values of x where
- * logGamma loses the sign information.
- * @param x The value for which we want Gamma(x)
- * @return The value of the gamma function evaluated at x.
- * */
- public static double gamma(double x) {
+ /**
+ * Returns the beta function of the arguments.
+ * <pre>
+ * - -
+ * | (a) | (b)
+ * beta( a, b ) = -----------.
+ * -
+ * | (a+b)
+ * </pre>
+ * @param alpha
+ * @param beta
+ * @return The beta function for given values of alpha and beta.
+ */
+ @Deprecated
+ public static double beta(double alpha, double beta) {
+
+ double y = alpha + beta;
+ y = gamma(y);
+ if (y == 0.0) {
+ return 1.0;
+ }
+
+ if (alpha > beta) {
+ y = gamma(alpha) / y;
+ y *= gamma(beta);
+ } else {
+ y = gamma(beta) / y;
+ y *= gamma(alpha);
+ }
+
+ return (y);
+ }
+
+ /** Returns the Gamma function of the argument. */
+ @Deprecated
+ public static double gamma(double x) throws ArithmeticException {
double[] P = {
1.60119522476751861407E-4,
@@ -132,136 +160,365 @@ public class Gamma extends org.apache.ma
}
/**
- * Power series for incomplete beta integral; formerly named <tt>pseries</tt>. Use when b*x is small and x not too
- * close to 1.
+ * Returns the Incomplete Beta Function evaluated from zero to <tt>xx</tt>; formerly named <tt>ibeta</tt>.
+ *
+ * @param aa the alpha parameter of the beta distribution.
+ * @param bb the beta parameter of the beta distribution.
+ * @param xx the integration end point.
*/
@Deprecated
- static double powerSeries(double a, double b, double x) {
+ public static double incompleteBeta(double aa, double bb, double xx) throws ArithmeticException {
- double ai = 1.0 / a;
- double u = (1.0 - b) * x;
- double v = u / (a + 1.0);
- double t1 = v;
- double t = u;
- double n = 2.0;
- double s = 0.0;
- double z = MACHEP * ai;
- while (Math.abs(v) > z) {
- u = (n - b) * x / n;
- t *= u;
- v = t / (a + n);
- s += v;
- n += 1.0;
+ if (aa <= 0.0 || bb <= 0.0) {
+ throw new
+ ArithmeticException("ibeta: Domain error!");
}
- s += t1;
- s += ai;
- u = a * Math.log(x);
- if ((a + b) < MAXGAM && Math.abs(u) < MAXLOG) {
- t = gamma(a + b) / (gamma(a) * gamma(b));
- s = s * t * Math.pow(x, a);
- } else {
- t = logGamma(a + b) - logGamma(a) - logGamma(b) + u + Math.log(s);
- if (t < MINLOG) {
- s = 0.0;
- } else {
- s = Math.exp(t);
+ if ((xx <= 0.0) || (xx >= 1.0)) {
+ if (xx == 0.0) {
+ return 0.0;
+ }
+ if (xx == 1.0) {
+ return 1.0;
}
+ throw new ArithmeticException("ibeta: Domain error!");
}
- return s;
- }
- /**
- * Returns the Incomplete Gamma function; formerly named <tt>igamma</tt>.
- *
- * @param alpha the shape parameter of the gamma distribution.
- * @param x the integration end point.
- * @return The value of the unnormalized incomplete gamma function.
- */
- public static double incompleteGamma(double alpha, double x){
- if (x <= 0 || alpha <= 0) {
- return 0.0;
+ double t;
+ if ((bb * xx) <= 1.0 && xx <= 0.95) {
+ t = powerSeries(aa, bb, xx);
+ return t;
}
- if (x > 1.0 && x > alpha) {
- return 1.0 - incompleteGammaComplement(alpha, x);
- }
+ double w = 1.0 - xx;
- /* Compute x**a * exp(-x) / gamma(a) */
- double ax = alpha * Math.log(x) - x - logGamma(alpha);
- if (ax < -MAXLOG) {
- return (0.0);
+ /* Reverse a and b if x is greater than the mean. */
+ double xc;
+ double x;
+ double b;
+ double a;
+ boolean flag = false;
+ if (xx > (aa / (aa + bb))) {
+ flag = true;
+ a = bb;
+ b = aa;
+ xc = xx;
+ x = w;
+ } else {
+ a = aa;
+ b = bb;
+ xc = w;
+ x = xx;
}
- ax = Math.exp(ax);
-
- /* power series */
- double r = alpha;
- double c = 1.0;
- double ans = 1.0;
-
- do {
- r += 1.0;
- c *= x / r;
- ans += c;
+ if (flag && (b * x) <= 1.0 && x <= 0.95) {
+ t = powerSeries(a, b, x);
+ if (t <= MACHEP) {
+ t = 1.0 - MACHEP;
+ } else {
+ t = 1.0 - t;
+ }
+ return t;
}
- while (c / ans > MACHEP);
-
- return (ans * ax / alpha);
- }
+ /* Choose expansion for better convergence. */
+ double y = x * (a + b - 2.0) - (a - 1.0);
+ if (y < 0.0) {
+ w = incompleteBetaFraction1(a, b, x);
+ } else {
+ w = incompleteBetaFraction2(a, b, x) / xc;
+ }
- /**
- * Returns the Complemented Incomplete Gamma function; formerly named <tt>igamc</tt>.
- *
- * @param alpha the shape parameter of the gamma distribution.
- * @param x the integration start point.
- */
- public static double incompleteGammaComplement(double alpha, double x) {
+ /* Multiply w by the factor
+ a b _ _ _
+ x (1-x) | (a+b) / ( a | (a) | (b) ) . */
- if (x <= 0 || alpha <= 0) {
- return 1.0;
+ y = a * Math.log(x);
+ t = b * Math.log(xc);
+ if ((a + b) < MAXGAM && Math.abs(y) < MAXLOG && Math.abs(t) < MAXLOG) {
+ t = Math.pow(xc, b);
+ t *= Math.pow(x, a);
+ t /= a;
+ t *= w;
+ t *= gamma(a + b) / (gamma(a) * gamma(b));
+ if (flag) {
+ if (t <= MACHEP) {
+ t = 1.0 - MACHEP;
+ } else {
+ t = 1.0 - t;
+ }
+ }
+ return t;
}
-
- if (x < 1.0 || x < alpha) {
- return 1.0 - incompleteGamma(alpha, x);
+ /* Resort to logarithms. */
+ y += t + logGamma(a + b) - logGamma(a) - logGamma(b);
+ y += Math.log(w / a);
+ if (y < MINLOG) {
+ t = 0.0;
+ } else {
+ t = Math.exp(y);
}
- double ax = alpha * Math.log(x) - x - logGamma(alpha);
- if (ax < -MAXLOG) {
- return 0.0;
+ if (flag) {
+ if (t <= MACHEP) {
+ t = 1.0 - MACHEP;
+ } else {
+ t = 1.0 - t;
+ }
}
+ return t;
+ }
- ax = Math.exp(ax);
+ /** Continued fraction expansion #1 for incomplete beta integral; formerly named <tt>incbcf</tt>. */
+ @Deprecated
+ static double incompleteBetaFraction1(double a, double b, double x) throws ArithmeticException {
- /* continued fraction */
- double y = 1.0 - alpha;
- double z = x + y + 1.0;
- double c = 0.0;
- double pkm2 = 1.0;
- double qkm2 = x;
- double pkm1 = x + 1.0;
- double qkm1 = z * x;
- double ans = pkm1 / qkm1;
+ double k1 = a;
+ double k2 = a + b;
+ double k3 = a;
+ double k4 = a + 1.0;
+ double k5 = 1.0;
+ double k6 = b - 1.0;
+ double k7 = k4;
+ double k8 = a + 2.0;
- double t;
+ double pkm2 = 0.0;
+ double qkm2 = 1.0;
+ double pkm1 = 1.0;
+ double qkm1 = 1.0;
+ double ans = 1.0;
+ double r = 1.0;
+ int n = 0;
+ double thresh = 3.0 * MACHEP;
do {
- c += 1.0;
- y += 1.0;
- z += 2.0;
- double yc = y * c;
- double pk = pkm1 * z - pkm2 * yc;
- double qk = qkm1 * z - qkm2 * yc;
+ double xk = -(x * k1 * k2) / (k3 * k4);
+ double pk = pkm1 + pkm2 * xk;
+ double qk = qkm1 + qkm2 * xk;
+ pkm2 = pkm1;
+ pkm1 = pk;
+ qkm2 = qkm1;
+ qkm1 = qk;
+
+ xk = (x * k5 * k6) / (k7 * k8);
+ pk = pkm1 + pkm2 * xk;
+ qk = qkm1 + qkm2 * xk;
+ pkm2 = pkm1;
+ pkm1 = pk;
+ qkm2 = qkm1;
+ qkm1 = qk;
+
if (qk != 0) {
- double r = pk / qk;
+ r = pk / qk;
+ }
+ double t;
+ if (r != 0) {
t = Math.abs((ans - r) / r);
ans = r;
} else {
t = 1.0;
}
- pkm2 = pkm1;
- pkm1 = pk;
+ if (t < thresh) {
+ return ans;
+ }
+
+ k1 += 1.0;
+ k2 += 1.0;
+ k3 += 2.0;
+ k4 += 2.0;
+ k5 += 1.0;
+ k6 -= 1.0;
+ k7 += 2.0;
+ k8 += 2.0;
+
+ if ((Math.abs(qk) + Math.abs(pk)) > big) {
+ pkm2 *= biginv;
+ pkm1 *= biginv;
+ qkm2 *= biginv;
+ qkm1 *= biginv;
+ }
+ if ((Math.abs(qk) < biginv) || (Math.abs(pk) < biginv)) {
+ pkm2 *= big;
+ pkm1 *= big;
+ qkm2 *= big;
+ qkm1 *= big;
+ }
+ } while (++n < 300);
+
+ return ans;
+ }
+
+ /** Continued fraction expansion #2 for incomplete beta integral; formerly named <tt>incbd</tt>. */
+ @Deprecated
+ static double incompleteBetaFraction2(double a, double b, double x) throws ArithmeticException {
+
+ double k1 = a;
+ double k2 = b - 1.0;
+ double k3 = a;
+ double k4 = a + 1.0;
+ double k5 = 1.0;
+ double k6 = a + b;
+ double k7 = a + 1.0;
+ double k8 = a + 2.0;
+
+ double pkm2 = 0.0;
+ double qkm2 = 1.0;
+ double pkm1 = 1.0;
+ double qkm1 = 1.0;
+ double z = x / (1.0 - x);
+ double ans = 1.0;
+ double r = 1.0;
+ int n = 0;
+ double thresh = 3.0 * MACHEP;
+ do {
+ double xk = -(z * k1 * k2) / (k3 * k4);
+ double pk = pkm1 + pkm2 * xk;
+ double qk = qkm1 + qkm2 * xk;
+ pkm2 = pkm1;
+ pkm1 = pk;
+ qkm2 = qkm1;
+ qkm1 = qk;
+
+ xk = (z * k5 * k6) / (k7 * k8);
+ pk = pkm1 + pkm2 * xk;
+ qk = qkm1 + qkm2 * xk;
+ pkm2 = pkm1;
+ pkm1 = pk;
+ qkm2 = qkm1;
+ qkm1 = qk;
+
+ if (qk != 0) {
+ r = pk / qk;
+ }
+ double t;
+ if (r != 0) {
+ t = Math.abs((ans - r) / r);
+ ans = r;
+ } else {
+ t = 1.0;
+ }
+
+ if (t < thresh) {
+ return ans;
+ }
+
+ k1 += 1.0;
+ k2 -= 1.0;
+ k3 += 2.0;
+ k4 += 2.0;
+ k5 += 1.0;
+ k6 += 1.0;
+ k7 += 2.0;
+ k8 += 2.0;
+
+ if ((Math.abs(qk) + Math.abs(pk)) > big) {
+ pkm2 *= biginv;
+ pkm1 *= biginv;
+ qkm2 *= biginv;
+ qkm1 *= biginv;
+ }
+ if ((Math.abs(qk) < biginv) || (Math.abs(pk) < biginv)) {
+ pkm2 *= big;
+ pkm1 *= big;
+ qkm2 *= big;
+ qkm1 *= big;
+ }
+ } while (++n < 300);
+
+ return ans;
+ }
+
+ /**
+ * Returns the Incomplete Gamma function; formerly named <tt>igamma</tt>.
+ *
+ * @param alpha the shape parameter of the gamma distribution.
+ * @param x the integration end point.
+ * @return The value of the unnormalized incomplete gamma function.
+ */
+ public static double incompleteGamma(double alpha, double x){
+ if (x <= 0 || alpha <= 0) {
+ return 0.0;
+ }
+
+ if (x > 1.0 && x > alpha) {
+ return 1.0 - incompleteGammaComplement(alpha, x);
+ }
+
+ /* Compute x**a * exp(-x) / gamma(a) */
+ double ax = alpha * Math.log(x) - x - logGamma(alpha);
+ if (ax < -MAXLOG) {
+ return (0.0);
+ }
+
+ ax = Math.exp(ax);
+
+ /* power series */
+ double r = alpha;
+ double c = 1.0;
+ double ans = 1.0;
+
+ do {
+ r += 1.0;
+ c *= x / r;
+ ans += c;
+ }
+ while (c / ans > MACHEP);
+
+ return (ans * ax / alpha);
+
+ }
+
+ /**
+ * Returns the Complemented Incomplete Gamma function; formerly named <tt>igamc</tt>.
+ *
+ * @param alpha the shape parameter of the gamma distribution.
+ * @param x the integration start point.
+ */
+ public static double incompleteGammaComplement(double alpha, double x) {
+
+ if (x <= 0 || alpha <= 0) {
+ return 1.0;
+ }
+
+ if (x < 1.0 || x < alpha) {
+ return 1.0 - incompleteGamma(alpha, x);
+ }
+
+ double ax = alpha * Math.log(x) - x - logGamma(alpha);
+ if (ax < -MAXLOG) {
+ return 0.0;
+ }
+
+ ax = Math.exp(ax);
+
+ /* continued fraction */
+ double y = 1.0 - alpha;
+ double z = x + y + 1.0;
+ double c = 0.0;
+ double pkm2 = 1.0;
+ double qkm2 = x;
+ double pkm1 = x + 1.0;
+ double qkm1 = z * x;
+ double ans = pkm1 / qkm1;
+
+ double t;
+ do {
+ c += 1.0;
+ y += 1.0;
+ z += 2.0;
+ double yc = y * c;
+ double pk = pkm1 * z - pkm2 * yc;
+ double qk = qkm1 * z - qkm2 * yc;
+ if (qk != 0) {
+ double r = pk / qk;
+ t = Math.abs((ans - r) / r);
+ ans = r;
+ } else {
+ t = 1.0;
+ }
+
+ pkm2 = pkm1;
+ pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
if (Math.abs(pk) > big) {
@@ -373,10 +630,51 @@ public class Gamma extends org.apache.ma
}
/**
+ * Power series for incomplete beta integral; formerly named <tt>pseries</tt>. Use when b*x is small and x not too
+ * close to 1.
+ */
+ @Deprecated
+ static double powerSeries(double a, double b, double x) throws ArithmeticException {
+
+ double ai = 1.0 / a;
+ double u = (1.0 - b) * x;
+ double v = u / (a + 1.0);
+ double t1 = v;
+ double t = u;
+ double n = 2.0;
+ double s = 0.0;
+ double z = MACHEP * ai;
+ while (Math.abs(v) > z) {
+ u = (n - b) * x / n;
+ t *= u;
+ v = t / (a + n);
+ s += v;
+ n += 1.0;
+ }
+ s += t1;
+ s += ai;
+
+ u = a * Math.log(x);
+ if ((a + b) < MAXGAM && Math.abs(u) < MAXLOG) {
+ t = gamma(a + b) / (gamma(a) * gamma(b));
+ s = s * t * Math.pow(x, a);
+ } else {
+ t = logGamma(a + b) - logGamma(a) - logGamma(b) + u + Math.log(s);
+ if (t < MINLOG) {
+ s = 0.0;
+ } else {
+ s = Math.exp(t);
+ }
+ }
+ return s;
+ }
+
+ /**
* Returns the Gamma function computed by Stirling's formula; formerly named <tt>stirf</tt>. The polynomial STIR is
* valid for 33 <= x <= 172.
*/
- private static double stirlingFormula(double x) {
+ @Deprecated
+ static double stirlingFormula(double x) throws ArithmeticException {
double[] STIR = {
7.87311395793093628397E-4,
-2.29549961613378126380E-4,
@@ -401,306 +699,4 @@ public class Gamma extends org.apache.ma
y = SQTPI * y * w;
return y;
}
-
- /**
- * Returns the beta function of the arguments.
- * <pre>
- * - -
- * | (a) | (b)
- * beta( a, b ) = -----------.
- * -
- * | (a+b)
- * </pre>
- * @param alpha
- * @param beta
- * @return The beta function for given values of alpha and beta.
- */
- @Deprecated
- public static double beta(double alpha, double beta) {
-
- double y = alpha + beta;
- y = gamma(y);
- if (y == 0.0) {
- return 1.0;
- }
-
- if (alpha > beta) {
- y = gamma(alpha) / y;
- y *= gamma(beta);
- } else {
- y = gamma(beta) / y;
- y *= gamma(alpha);
- }
-
- return (y);
- }
-
- /**
- * Returns the Incomplete Beta Function evaluated from zero to <tt>xx</tt>; formerly named <tt>ibeta</tt>.
- *
- * @param aa the alpha parameter of the beta distribution.
- * @param bb the beta parameter of the beta distribution.
- * @param xx the integration end point.
- */
- @Deprecated
- public static double incompleteBeta(double aa, double bb, double xx) throws ArithmeticException {
-
- if (aa <= 0.0 || bb <= 0.0) {
- throw new
- ArithmeticException("ibeta: Domain error!");
- }
-
- if ((xx <= 0.0) || (xx >= 1.0)) {
- if (xx == 0.0) {
- return 0.0;
- }
- if (xx == 1.0) {
- return 1.0;
- }
- throw new ArithmeticException("ibeta: Domain error!");
- }
-
- double t;
- if ((bb * xx) <= 1.0 && xx <= 0.95) {
- t = powerSeries(aa, bb, xx);
- return t;
- }
-
- double w = 1.0 - xx;
-
- /* Reverse a and b if x is greater than the mean. */
- double xc;
- double x;
- double b;
- double a;
- boolean flag = false;
- if (xx > (aa / (aa + bb))) {
- flag = true;
- a = bb;
- b = aa;
- xc = xx;
- x = w;
- } else {
- a = aa;
- b = bb;
- xc = w;
- x = xx;
- }
-
- if (flag && (b * x) <= 1.0 && x <= 0.95) {
- t = powerSeries(a, b, x);
- if (t <= MACHEP) {
- t = 1.0 - MACHEP;
- } else {
- t = 1.0 - t;
- }
- return t;
- }
-
- /* Choose expansion for better convergence. */
- double y = x * (a + b - 2.0) - (a - 1.0);
- if (y < 0.0) {
- w = incompleteBetaFraction1(a, b, x);
- } else {
- w = incompleteBetaFraction2(a, b, x) / xc;
- }
-
- /* Multiply w by the factor
- a b _ _ _
- x (1-x) | (a+b) / ( a | (a) | (b) ) . */
-
- y = a * Math.log(x);
- t = b * Math.log(xc);
- if ((a + b) < MAXGAM && Math.abs(y) < MAXLOG && Math.abs(t) < MAXLOG) {
- t = Math.pow(xc, b);
- t *= Math.pow(x, a);
- t /= a;
- t *= w;
- t *= gamma(a + b) / (gamma(a) * gamma(b));
- if (flag) {
- if (t <= MACHEP) {
- t = 1.0 - MACHEP;
- } else {
- t = 1.0 - t;
- }
- }
- return t;
- }
- /* Resort to logarithms. */
- y += t + logGamma(a + b) - logGamma(a) - logGamma(b);
- y += Math.log(w / a);
- if (y < MINLOG) {
- t = 0.0;
- } else {
- t = Math.exp(y);
- }
-
- if (flag) {
- if (t <= MACHEP) {
- t = 1.0 - MACHEP;
- } else {
- t = 1.0 - t;
- }
- }
- return t;
- }
-
- /** Continued fraction expansion #1 for incomplete beta integral; formerly named <tt>incbcf</tt>. */
- @Deprecated
- private static double incompleteBetaFraction1(double a, double b, double x) throws ArithmeticException {
-
- double k1 = a;
- double k2 = a + b;
- double k3 = a;
- double k4 = a + 1.0;
- double k5 = 1.0;
- double k6 = b - 1.0;
- double k7 = k4;
- double k8 = a + 2.0;
-
- double pkm2 = 0.0;
- double qkm2 = 1.0;
- double pkm1 = 1.0;
- double qkm1 = 1.0;
- double ans = 1.0;
- double r = 1.0;
- int n = 0;
- double thresh = 3.0 * MACHEP;
- do {
- double xk = -(x * k1 * k2) / (k3 * k4);
- double pk = pkm1 + pkm2 * xk;
- double qk = qkm1 + qkm2 * xk;
- pkm2 = pkm1;
- pkm1 = pk;
- qkm2 = qkm1;
- qkm1 = qk;
-
- xk = (x * k5 * k6) / (k7 * k8);
- pk = pkm1 + pkm2 * xk;
- qk = qkm1 + qkm2 * xk;
- pkm2 = pkm1;
- pkm1 = pk;
- qkm2 = qkm1;
- qkm1 = qk;
-
- if (qk != 0) {
- r = pk / qk;
- }
- double t;
- if (r != 0) {
- t = Math.abs((ans - r) / r);
- ans = r;
- } else {
- t = 1.0;
- }
-
- if (t < thresh) {
- return ans;
- }
-
- k1 += 1.0;
- k2 += 1.0;
- k3 += 2.0;
- k4 += 2.0;
- k5 += 1.0;
- k6 -= 1.0;
- k7 += 2.0;
- k8 += 2.0;
-
- if ((Math.abs(qk) + Math.abs(pk)) > big) {
- pkm2 *= biginv;
- pkm1 *= biginv;
- qkm2 *= biginv;
- qkm1 *= biginv;
- }
- if ((Math.abs(qk) < biginv) || (Math.abs(pk) < biginv)) {
- pkm2 *= big;
- pkm1 *= big;
- qkm2 *= big;
- qkm1 *= big;
- }
- } while (++n < 300);
-
- return ans;
- }
-
- /** Continued fraction expansion #2 for incomplete beta integral; formerly named <tt>incbd</tt>. */
- @Deprecated
- private static double incompleteBetaFraction2(double a, double b, double x) throws ArithmeticException {
-
- double k1 = a;
- double k2 = b - 1.0;
- double k3 = a;
- double k4 = a + 1.0;
- double k5 = 1.0;
- double k6 = a + b;
- double k7 = a + 1.0;
- double k8 = a + 2.0;
-
- double pkm2 = 0.0;
- double qkm2 = 1.0;
- double pkm1 = 1.0;
- double qkm1 = 1.0;
- double z = x / (1.0 - x);
- double ans = 1.0;
- double r = 1.0;
- int n = 0;
- double thresh = 3.0 * MACHEP;
- do {
- double xk = -(z * k1 * k2) / (k3 * k4);
- double pk = pkm1 + pkm2 * xk;
- double qk = qkm1 + qkm2 * xk;
- pkm2 = pkm1;
- pkm1 = pk;
- qkm2 = qkm1;
- qkm1 = qk;
-
- xk = (z * k5 * k6) / (k7 * k8);
- pk = pkm1 + pkm2 * xk;
- qk = qkm1 + qkm2 * xk;
- pkm2 = pkm1;
- pkm1 = pk;
- qkm2 = qkm1;
- qkm1 = qk;
-
- if (qk != 0) {
- r = pk / qk;
- }
- double t;
- if (r != 0) {
- t = Math.abs((ans - r) / r);
- ans = r;
- } else {
- t = 1.0;
- }
-
- if (t < thresh) {
- return ans;
- }
-
- k1 += 1.0;
- k2 -= 1.0;
- k3 += 2.0;
- k4 += 2.0;
- k5 += 1.0;
- k6 += 1.0;
- k7 += 2.0;
- k8 += 2.0;
-
- if ((Math.abs(qk) + Math.abs(pk)) > big) {
- pkm2 *= biginv;
- pkm1 *= biginv;
- qkm2 *= biginv;
- qkm1 *= biginv;
- }
- if ((Math.abs(qk) < biginv) || (Math.abs(pk) < biginv)) {
- pkm2 *= big;
- pkm1 *= big;
- qkm2 *= big;
- qkm1 *= big;
- }
- } while (++n < 300);
-
- return ans;
- }
}