Writing a trigamma function for Splus

1998-11-14

Finding the formulas

Suitable formulas can be found in Abramowitz & Stegun, Handbook of Mathematical Functions, published by the National Bureau of Standards and reprinted by Dover, $40.95 at the Bookstore.

Formula 6.4.12 is an asymptotic approximation to trigamma(z) for large z

trigamma(z) ~ 1/z + 1/(2 * z^2) + 1/(6 * z^3) - 1/(30 * z^5) + 1/(42 * z^7) - 1/(30 * z^9) + ...

How large does z need to be? If z = 30 the first term will be 0.0333 and the sixth term will be 1.69e-15, so the first six terms will be enough to give at least 12-digit accuracy as long as z > 30.

To use this formula for z < 30, apply the recursion of formula 6.4.6 until the argument exceeds 30.

trigamma(z+1) = trigamma(z) - 1/(z^2)

Programming the formulas in FORTRAN or C

Here are these calculations coded as a subroutine in FORTRAN 77 and saved in a file named trigama.f. C could be used instead of FORTRAN and the the code would look much the same, but the file would be called trigama.c. Use an editor (pico or vi) to write the code. The name of the file does not have to be the same as the name of the subroutine.

The subroutine trigama(z,trigg) has two double-precision real arguments: the first, z, is the input value, the second, trigg, is the output value, the result of the calculation.

s4p3@stats[62] % cat trigama.f
C  Apply (6.4.6) of Abramowitz & Stegun recursively while z < 30,
C  then use (6.4.12) to compute the trigamma function.
 
      subroutine trigama(z,trigg)
      implicit none
      real*8 z,zp,zpr,trig,trigg
      trig=0.d0
      if(z.gt.0.d0) then
        zp=z
        zpr=1.d0/zp
        do while (zp.lt.30.d0)
          trig=trig+(zpr*zpr)
          zp=zp+1.d0
          zpr=1.d0/zp
        end do
        trig=trig+(zpr)+((zpr**2)/2.d0)+((zpr**3)/6.d0)
     +   -((zpr**5)/30.d0)
     +   +((zpr**7)/42.d0)
     +   -((zpr**9)/30.d0)
      end if
      trigg=trig
      return
      end

Compiling the source code

The following command line in UNIX invokes the FORTRAN 77 compiler to produce object code and save it in a file named (by default) trigama.o.

s4p3@stats[61] % f77 -c trigama.f

If you wrote the source code in C, you would use the command line

s4p3@stats[61] % cc -c trigama.c

Calling the FORTRAN or C subroutine from Splus

The final step is to write an Splus function that calls the FORTRAN or C subroutine and returns the value of trigamma. Since the FORTRAN or C routine expects a scalar input argument but we want the Splus function to work for any input (scalar, vector, matrix or list), my Splus function trigamma() defines a function trig() internally and uses sapply() to apply trig() to whatever object zz is given to trigamma().

> trigamma
function(zz)
{
        trig <- function(z)
        {
                if(z > 0) {
                        # Note 1
                        trigg <- as.double(z)
                        # Note 2
                        dyn.open("/home/stats/s4p3/myslib/trigama.o")
                        # Note 3
                        tr <- .Fortran("trigama",
                                as.double(z),
                                trigg = as.double(trigg))
                        # Note 4
                        tr$trigg
                }
                else NA
        }
        sapply(zz, trig)
}

Note 1

Define a numeric object trigg before calling the FORTRAN subroutine; it will be used to receive the result of the calculation.

Note 2

Dynamically open the file that holds the object code for the FORTRAN subroutine. To be safe, give the complete path to the file.

Note 3

To call the FORTRAN subroutine, specify three things: (1) the name of the subroutine (in this case it is the same as the name of the file but that isn't necessary; besides, the file could be a library with several FORTRAN subroutines in it); (2) which object (z) is to be used as the first argument of the FORTRAN subroutine; and (3) which argument (trigg) is to be used as the second argument of the FORTRAN subroutine. It was convenient to give the objects in the Splus the same names as the corresponding variables in the FORTRAN code, but the names don't have to be the same.

Note 4

The result of the calculation is now available within trig() as the Splus object tr$trigg; this is the value returned by trig() and used in sapply().

Testing the trigamma() function in Splus

Values in the following table can be compared with those in Table 6.1 of Abramowitz & Stegun. They agree to 7 decimal places.

> matrix(c(seq(1,2,.1),trigamma(seq(1,2,.1))),ncol=2)
      [,1]      [,2]
 [1,]  1.0 1.6449341
 [2,]  1.1 1.4332992
 [3,]  1.2 1.2673772
 [4,]  1.3 1.1342534
 [5,]  1.4 1.0253566
 [6,]  1.5 0.9348022
 [7,]  1.6 0.8584319
 [8,]  1.7 0.7932328
 [9,]  1.8 0.7369741
[10,]  1.9 0.6879721
[11,]  2.0 0.6449341

---

Back to Statistics 743