DPQ (version 0.5-3)

dpsifn: Psi Gamma Functions Workhorse from R's API

Description

Log Gamma derivatives, Psi Gamma functions. dpsifn() is an R interface to the R API function R_dpsifn().

Usage

dpsifn(x, m, deriv1 = 0L, k2 = FALSE)

Value

A numeric \(l_x \times m\)

matrix (where

\(l_x=\)length(x)) of scaled \(\psi^{(k)}(x)\)

values. The matrix has attributes

underflow

of \(l_x\) integer counts of the number of under- and over-flows, in computing the corresponding i-th matrix column for x[i].

ierr

length-\(l_x\) integer vector of error codes, where 0 is normal/successful.

Arguments

x

numeric vector.

m

number of derivatives to return, an integer >= 0.

deriv1

“start” derivative ....

k2

a logical specifying if kode = 2 should be applied.

Author

Martin Maechler (R interface); R Core et al., see digamma.

Details

dpsifn() is the underlying “workhorse” of R's own digamma, trigamma and (generalized) psigamma functions.

It is useful, e.g., when several derivatives of \(\log\Gamma=\)lgamma are desired. It computes and returns length-m sequence \((-1)^{k+1} / \Gamma(k+1) \cdot \psi^{(k)}(x)\) for \(k = n, n+1,\ldots, n+m-1\), where \(n=\)deriv1, and \(\psi^{(k)}(x)\) is the k-th derivative of \(\psi(x)\), i.e., psigamma(x,k). For more details, see the comments in src/nmath/polygamma.c.

References

See those in psigamma

See Also

digamma, trigamma, psigamma.

Examples

x <- seq(-3.5, 6, by=1/4)
dpx <- dpsifn(x, m = if(getRversion() >= 4.2) 7 else 5)
dpx # in R <= 4.2.1, see that sometimes the 'nz' (under-over-flow count) was uninitialized !!
j <- -1L+seq_len(nrow(dpx)); (fj <- (-1)^(j+1)*gamma(j+1))
## mdpsi <- cbind(di =   digamma(x),      -dpx[1,],
## 	       tri=  trigamma(x),       dpx[2,],
## 	       tetra=psigamma(x,2),  -2*dpx[3,],
## 	       penta=psigamma(x,3),   6*dpx[4,],
## 	       hexa =psigamma(x,4), -24*dpx[5,],
## 	       hepta=psigamma(x,5), 120*dpx[6,],
## 	       octa =psigamma(x,6),-720*dpx[7,])
## cbind(x, ie=attr(dpx,"errorCode"), round(mdpsi, 4))
str(psig <- outer(x, j, psigamma))
dpsi <- t(fj * (`attributes<-`(dpx, list(dim=dim(dpx)))))
if(getRversion() >= 4.2) {
      print( all.equal(psig, dpsi, tol=0) )# -> see 1.185e-16
  stopifnot( all.equal(psig, dpsi, tol=1e-15) )
} else { # R <= 4.1.x; dpsifn(x, ..) *not* ok for x < 0
  i <- x >= 0
      print( all.equal(psig[i,], dpsi[i,], tol=0) )# -> see 1.95e-16
  stopifnot( all.equal(psig[i,], dpsi[i,], tol=1e-15) )
}