View source: R/psigamma-deriv.R
dpsifn | R Documentation |
Log Gamma derivatives, Psi Gamma functions. dpsifn()
is an R
interface to the R API function R_dpsifn()
.
dpsifn(x, m, deriv1 = 0L, k2 = FALSE)
x |
numeric vector. |
m |
number of derivatives to return, an integer >= 0. |
deriv1 |
“start” derivative .... |
k2 |
a |
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’.
A numeric l_x \times m
matrix
(where
l_x=
length(x)
) of scaled \psi^{(k)}(x)
values. The matrix has attributes
underflow |
of |
ierr |
length- |
Martin Maechler (R interface); R Core et al., see
digamma
.
See those in psigamma
digamma
, trigamma
, psigamma
.
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) )
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.