Author: psteitz
Date: Sun May 24 06:03:19 2009
New Revision: 778091
URL: http://svn.apache.org/viewvc?rev=778091&view=rev
Log:
Added digamma function.
JIRA: MATH-267
Contributed by Ted Dunning
Modified:
commons/proper/math/trunk/src/java/org/apache/commons/math/special/Gamma.java
commons/proper/math/trunk/src/site/xdoc/changes.xml
commons/proper/math/trunk/src/test/org/apache/commons/math/special/GammaTest.java
Modified:
commons/proper/math/trunk/src/java/org/apache/commons/math/special/Gamma.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/special/Gamma.java?rev=778091&r1=778090&r2=778091&view=diff
==============================================================================
---
commons/proper/math/trunk/src/java/org/apache/commons/math/special/Gamma.java
(original)
+++
commons/proper/math/trunk/src/java/org/apache/commons/math/special/Gamma.java
Sun May 24 06:03:19 2009
@@ -33,6 +33,9 @@
/** Serializable version identifier */
private static final long serialVersionUID = -6587513359895466954L;
+ /** <a
href="http://en.wikipedia.org/wiki/Euler-Mascheroni_constant";>Euler-Mascheroni
constant</a> */
+ public static final double GAMMA = 0.577215664901532860606512090082;
+
/** Maximum allowed numerical error. */
private static final double DEFAULT_EPSILON = 10e-15;
@@ -59,7 +62,7 @@
/** Avoid repeated computation of log of 2 PI in logGamma */
private static final double HALF_LOG_2_PI = 0.5 * Math.log(2.0 * Math.PI);
-
+
/**
* Default constructor. Prohibit instantiation.
*/
@@ -261,4 +264,45 @@
return ret;
}
+
+
+ // limits for switching algorithm in digamma
+ /** C limit */
+ private static final double C_LIMIT = 49;
+ /** S limit */
+ private static final double S_LIMIT = 1e-5;
+
+ /**
+ * <p>Computes the <a
href="http://en.wikipedia.org/wiki/Digamma_function";>digamma function</a>
+ * using the algorithm defined in <br/>
+ * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied
Statistics, 1976.</p>
+ *
+ * <p>Some of the constants have been changed to increase accuracy at the
moderate expense
+ * of run-time performance. The result should be accurate to within 10^-8
absolute tolerance for
+ * x >= 10^-5 and within 10^-8 relative tolerance for x > 0.</p>
+ *
+ * <p> Performance for large negative values of x will be quite expensive
(proportional to
+ * |x|). Accuracy for negative values of x should be about 10^-8 absolute
for results
+ * less than 10^5 and 10^-8 relative for results larger than that.
+ * @param x argument
+ * @return value of the digamma function
+ */
+ public static double digamma(double x) {
+ if (x > 0 && x <= S_LIMIT) {
+ // use method 5 from Bernardo AS103
+ // accurate to O(x)
+ return -GAMMA - 1 / x;
+ }
+
+ if (x >= C_LIMIT) {
+ // use method 4 (accurate to O(1/x^8)
+ double inv = 1 / (x * x);
+ // 1 1 1 1
+ // log(x) - --- - ------ - ------- - -------
+ // 2 x 12 x^2 120 x^4 252 x^6
+ return Math.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 /
120 - inv / 252));
+ }
+
+ return digamma(x + 1) - 1 / x;
+ }
}
Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=778091&r1=778090&r2=778091&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Sun May 24 06:03:19 2009
@@ -39,6 +39,9 @@
</properties>
<body>
<release version="2.0" date="TBD" description="TBD">
+ <action dev="psteitz" type="add" issue="MATH-267" due=to="Ted Dunning">
+ Added digamma function.
+ </action>
<action dev="psteitz" type="add" issue="MATH-136" due=to="John Gant">
Added Spearman's rank correlation (SpearmansCorrelation).
</action>
Modified:
commons/proper/math/trunk/src/test/org/apache/commons/math/special/GammaTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/special/GammaTest.java?rev=778091&r1=778090&r2=778091&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/org/apache/commons/math/special/GammaTest.java
(original)
+++
commons/proper/math/trunk/src/test/org/apache/commons/math/special/GammaTest.java
Sun May 24 06:03:19 2009
@@ -25,10 +25,7 @@
* @version $Revision$ $Date$
*/
public class GammaTest extends TestCase {
- /**
- * Constructor for BetaTest.
- * @param name
- */
+
public GammaTest(String name) {
super(name);
}
@@ -92,4 +89,37 @@
public void testLogGammaPositive() {
testLogGamma(0.6931471805599457, 3.0);
}
+
+ public void testDigammaLargeArgs() {
+ double eps = 1e-8;
+ assertEquals(4.6001618527380874002, Gamma.digamma(100), eps);
+ assertEquals(3.9019896734278921970, Gamma.digamma(50), eps);
+ assertEquals(2.9705239922421490509, Gamma.digamma(20), eps);
+ assertEquals(2.9958363947076465821, Gamma.digamma(20.5), eps);
+ assertEquals(2.2622143570941481605, Gamma.digamma(10.1), eps);
+ assertEquals(2.1168588189004379233, Gamma.digamma(8.8), eps);
+ assertEquals(1.8727843350984671394, Gamma.digamma(7), eps);
+ assertEquals(0.42278433509846713939, Gamma.digamma(2), eps);
+ assertEquals(-100.56088545786867450, Gamma.digamma(0.01), eps);
+ assertEquals(-4.0390398965921882955, Gamma.digamma(-0.8), eps);
+ assertEquals(4.2003210041401844726, Gamma.digamma(-6.3), eps);
+ }
+
+ public void testDigammaSmallArgs() {
+ // values for negative powers of 10 from 1 to 30 as computed by
webMathematica with 20 digits
+ // see functions.wolfram.com
+ double[] expected = {-10.423754940411076795, -100.56088545786867450,
-1000.5755719318103005,
+ -10000.577051183514335, -100000.57719921568107,
-1.0000005772140199687e6, -1.0000000577215500408e7,
+ -1.0000000057721564845e8, -1.0000000005772156633e9,
-1.0000000000577215665e10, -1.0000000000057721566e11,
+ -1.0000000000005772157e12, -1.0000000000000577216e13,
-1.0000000000000057722e14, -1.0000000000000005772e15, -1e+16,
+ -1e+17, -1e+18, -1e+19, -1e+20, -1e+21, -1e+22, -1e+23,
-1e+24, -1e+25, -1e+26,
+ -1e+27, -1e+28, -1e+29, -1e+30};
+ for (double n = 1; n < 30; n++) {
+ checkRelativeError(String.format("Test %.0f: ", n), expected[(int)
(n - 1)], Gamma.digamma(Math.pow(10.0, -n)), 1e-8);
+ }
+ }
+
+ private void checkRelativeError(String msg, double expected, double
actual, double tolerance) {
+ assertEquals(msg, expected, actual, Math.abs(tolerance * actual));
+ }
}
