Author: psteitz
Date: Mon May 25 13:20:07 2009
New Revision: 778416
URL: http://svn.apache.org/viewvc?rev=778416&view=rev
Log:
Added trigamma, javadoc fixes for digamma. JIRA: MATH-267. Patched by Ted
Dunning.
Modified:
commons/proper/math/trunk/src/java/org/apache/commons/math/special/Gamma.java
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=778416&r1=778415&r2=778416&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
Mon May 25 13:20:07 2009
@@ -271,25 +271,28 @@
// 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;
+ 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/>
+ * <p>Computes the digamma function of x.</p>
+ *
+ * <p>This is an independently written implementation of the algorithm
described in
* 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
+ * of run-time. 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
+ * <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.
+ * less than 10^5 and 10^-8 relative for results larger than that.</p>
*
- * @param x argument
- * @return value of the digamma function
+ * @param x the argument
+ * @return digamma(x) to within 10-8 relative or absolute error
whichever is smaller
+ * @see <a href="http://en.wikipedia.org/wiki/Digamma_function";> Digamma
at wikipedia </a>
+ * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf";>
Bernardo's original article </a>
* @since 2.0
*/
public static double digamma(double x) {
@@ -303,11 +306,38 @@
// use method 4 (accurate to O(1/x^8)
double inv = 1 / (x * x);
// 1 1 1 1
- // log(x) - --- - ------ - ------- - -------
+ // 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;
}
+
+ /**
+ * <p>Computes the trigamma function of x. This function is derived by
taking the derivative of
+ * the implementation of digamma.</p>
+ *
+ * @param x the argument
+ * @return trigamma(x) to within 10-8 relative or absolute error
whichever is smaller
+ * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function";> Trigamma
at wikipedia </a>
+ * @see Gamma#digamma(double)
+ * @since 2.0
+ */
+ public static double trigamma(double x) {
+ if (x > 0 && x <= S_LIMIT) {
+ return 1 / (x * x);
+ }
+
+ if (x >= C_LIMIT) {
+ double inv = 1 / (x * x);
+ // 1 1 1 1 1
+ // - + ---- + ---- - ----- + -----
+ // x 2 3 5 7
+ // 2 x 6 x 30 x 42 x
+ return 1 / x + inv / 2 + inv / x * (1.0 / 6 - inv * (1.0 / 30 +
inv / 42));
+ }
+
+ return trigamma(x + 1) + 1 / (x * x);
+ }
}
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=778416&r1=778415&r2=778416&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
Mon May 25 13:20:07 2009
@@ -119,6 +119,31 @@
}
}
+ public void testTrigamma() {
+ double eps = 1e-8;
+ // computed using webMathematica. For example, to compute
trigamma($i) = Polygamma(1, $i), use
+ //
+ //
http://functions.wolfram.com/webMathematica/Evaluated.jsp?name=PolyGamma2&plottype=0&vars={%221%22,%22$i%22}&digits=20
+ double[] data = {
+ 1e-4, 1.0000000164469368793e8,
+ 1e-3, 1.0000016425331958690e6,
+ 1e-2, 10001.621213528313220,
+ 1e-1, 101.43329915079275882,
+ 1, 1.6449340668482264365,
+ 2, 0.64493406684822643647,
+ 3, 0.39493406684822643647,
+ 4, 0.28382295573711532536,
+ 5, 0.22132295573711532536,
+ 10, 0.10516633568168574612,
+ 20, 0.051270822935203119832,
+ 50, 0.020201333226697125806,
+ 100, 0.010050166663333571395
+ };
+ for (int i = data.length - 2; i >= 0; i -= 2) {
+ assertEquals(String.format("trigamma %.0f", data[i]), data[i + 1],
Gamma.trigamma(data[i]), eps);
+ }
+ }
+
private void checkRelativeError(String msg, double expected, double
actual, double tolerance) {
assertEquals(msg, expected, actual, Math.abs(tolerance * actual));
}
