Special                 package:base                 R Documentation

_S_p_e_c_i_a_l _F_u_n_c_t_i_o_n_s _o_f _M_a_t_h_e_m_a_t_i_c_s

_D_e_s_c_r_i_p_t_i_o_n:

     Special mathematical functions related to the beta and gamma
     functions.

_U_s_a_g_e:

     beta(a, b)
     lbeta(a, b)

     gamma(x)
     lgamma(x)
     psigamma(x, deriv = 0)
     digamma(x)
     trigamma(x)

     choose(n, k)
     lchoose(n, k)
     factorial(x)
     lfactorial(x)

_A_r_g_u_m_e_n_t_s:

    a, b: non-negative numeric vectors.

    x, n: numeric vectors.

k, deriv: integer vectors.

_D_e_t_a_i_l_s:

     The functions 'beta' and 'lbeta' return the beta function and the
     natural logarithm of the beta function,

              B(a,b) = (Gamma(a)Gamma(b))/(Gamma(a+b)).

     The formal definition is

                 integral_0^1 t^(a-1) (1-t)^(b-1) dt

     (Abramowitz and Stegun section 6.2.1, page 258).  Note that it is
     only defined in R for non-negative 'a' and 'b', and is infinite if
     either is zero.

     The functions 'gamma' and 'lgamma' return the gamma function
     Gamma(x) and the natural logarithm of _the absolute value of_ the
     gamma function.  The gamma function is defined by (Abramowitz and
     Stegun section 6.1.1, page 255)

                  integral_0^Inf t^(a-1) exp(-t) dt

     for all real 'x' except zero and negative integers (when 'NaN' is
     returned). 'factorial(x)' (x! for non-negative integer 'x') is
     defined to be 'gamma(x+1)' and 'lfactorial' to be 'lgamma(x+1)'.

     The functions 'digamma' and 'trigamma' return the first and second
     derivatives of the logarithm of the gamma function. 'psigamma(x,
     deriv)' ('deriv >= 0') computes the 'deriv'-th derivative of
     psi(x).

  'digamma(x)' = psi(x) = d/dx {ln Gamma(x)} = Gamma'(x) / Gamma(x)

     This is often called the 'polygamma' function, e.g. in Abramowitz
     and Stegun (section 6.4.1, page 260); and its higher derivatives
     ('deriv = 2:4') have occasionally been called 'tetragamma',
     'pentagamma', and 'hexagamma'.

     The functions 'choose' and 'lchoose' return binomial coefficients
     and their logarithms.  Note that 'choose(n,k)' is defined for all
     real numbers n and integer k.  For k >= 1 as n(n-1)...(n-k+1) /
     k!, as 1 for k = 0 and as 0 for negative k. Non-integer values of
     'k' are rounded to an integer, with a warning. 
      'choose(*,k)' uses direct arithmetic (instead of '[l]gamma'
     calls) for small 'k', for speed and accuracy reasons.  Note the
     function 'combn' (package 'utils') for enumeration of all possible
     combinations.

     The 'gamma', 'lgamma', 'digamma' and 'trigamma' functions are
     generic: methods can be defined for them individually or via the
     'Math' group generic.

_S_o_u_r_c_e:

     'gamma', 'lgamma', 'beta' and 'lbeta' are based on C translations
     of Fortran subroutines by W. Fullerton of Los Alamos Scientific
     Laboratory (now available as part of SLATEC).

     'digamma', 'trigamma' and 'psigamma' are based on

     Amos, D. E. (1983). A portable Fortran subroutine for derivatives
     of the psi function, Algorithm 610, _ACM Transactions on
     Mathematical Software_ *9(4)*, 494-502.

_R_e_f_e_r_e_n_c_e_s:

     Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) _The New S
     Language_. Wadsworth & Brooks/Cole. (For 'gamma' and 'lgamma'.)

     Abramowitz, M. and Stegun, I. A. (1972) _Handbook of Mathematical
     Functions._ New York: Dover. Chapter 6: Gamma and Related
     Functions.

_S_e_e _A_l_s_o:

     'Arithmetic' for simple, 'sqrt' for miscellaneous mathematical
     functions and 'Bessel' for the real Bessel functions.  Note that
     'gammaCody(x)' is considerably faster than 'gamma(x)' but slightly
     less accurate and (potentially) less reliable.

     For the incomplete gamma function see 'pgamma'.

_E_x_a_m_p_l_e_s:

     require(graphics)

     choose(5, 2)
     for (n in 0:10) print(choose(n, k = 0:n))

     factorial(100)
     lfactorial(10000)

     ## gamma has 1st order poles at 0, -1, -2, ...
     ## this will generate loss of precision warnings, so turn off
     op <- options("warn")
     options(warn = -1)
     x <- sort(c(seq(-3,4, length.out=201), outer(0:-3, (-1:1)*1e-6, "+")))
     plot(x, gamma(x), ylim=c(-20,20), col="red", type="l", lwd=2,
          main=expression(Gamma(x)))
     abline(h=0, v=-3:0, lty=3, col="midnightblue")
     options(op)

     x <- seq(.1, 4, length.out = 201); dx <- diff(x)[1]
     par(mfrow = c(2, 3))
     for (ch in c("", "l","di","tri","tetra","penta")) {
       is.deriv <- nchar(ch) >= 2
       nm <- paste(ch, "gamma", sep = "")
       if (is.deriv) {
         dy <- diff(y) / dx # finite difference
         der <- which(ch == c("di","tri","tetra","penta")) - 1
         nm2 <- paste("psigamma(*, deriv = ", der,")",sep='')
         nm  <- if(der >= 2) nm2 else paste(nm, nm2, sep = " ==\n")
         y <- psigamma(x, deriv=der)
       } else {
         y <- get(nm)(x)
       }
       plot(x, y, type = "l", main = nm, col = "red")
       abline(h = 0, col = "lightgray")
       if (is.deriv) lines(x[-1], dy, col = "blue", lty = 2)
     }
     par(mfrow = c(1, 1))

     ## "Extended" Pascal triangle:
     fN <- function(n) formatC(n, width=2)
     for (n in -4:10) cat(fN(n),":", fN(choose(n, k= -2:max(3,n+2))), "\n")

     ## R code version of choose()  [simplistic; warning for k < 0]:
     mychoose <- function(r,k)
         ifelse(k <= 0, (k==0),
                sapply(k, function(k) prod(r:(r-k+1))) / factorial(k))
     k <- -1:6
     cbind(k=k, choose(1/2, k), mychoose(1/2, k))

     ## Binomial theorem for n=1/2 ;
     ## sqrt(1+x) = (1+x)^(1/2) = sum_{k=0}^Inf  choose(1/2, k) * x^k :
     k <- 0:10 # 10 is sufficient for ~ 9 digit precision:
     sqrt(1.25)
     sum(choose(1/2, k)* .25^k)

