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.

     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)