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, 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 (6.2.1), page 258).

     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 (6.1.1), page 255)

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

     'factorial(x)' is x! and identical to 'gamma(x+1)' and
     'lfactorial' is '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') is more generally computing the 'deriv'-th
     derivative of psi(x).

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


     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. 
      'choose(*,k)' uses direct arithmetic (instead of '[l]gamma'
     calls) for small 'k', for speed and accuracy reasons.

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

_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:

     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, ...
     x <- sort(c(seq(-3,4, length=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")

     x <- seq(.1, 4, length = 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)
     }

     ## "Extended" Pascal triangle:
     fN <- function(n) formatC(n, wid=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)