copula (version 1.1-6)

polylog: Polylogarithm ${\mathrm{Li}}_{\mathrm{s}}(\mathrm{z})$ and Debye Functions

Description

Compute the polylogarithm function ${\mathrm{Li}}_{\mathrm{s}}(\mathrm{z})$, initially defined as the power series, $${\mathrm{Li}}_s(z)=\sum _{k=1}^{\mathrm{\infty }}\frac{z^k}{k^s},$$ for $|z|<1$, and then more generally (by analytic continuation) as $${\mathrm{Li}}_1(z)=-\log (1-z),$$ and $${\mathrm{Li}}_{s+1}(z)={\int }_0^z\frac{{\mathrm{Li}}_s(t)}{t}dt.$$

Currently, mainly the case of negative integer $s$ is well supported, as that is used for some of the Archimedean copula densities.

For $s=2$, ${\mathrm{Li}}_2(\mathrm{z})$ is also called ‘dilogarithm’ or “Spence's function”. The "default" method uses the dilog or complex_dilog function from package gsl, respectively when $s=2$.

Also compute the Debye_n functions, for $n=1$ and $n=2$, in a slightly more general manner than the gsl package functions debye_1 and debye_2 (which cannot deal with non-finite x.)

Usage

polylog(z, s,
        method = c("default", "sum", "negI-s-Stirling",
                   "negI-s-Eulerian", "negI-s-asymp-w"),
        logarithm = FALSE, is.log.z = FALSE, is.logmlog = FALSE,
        asymp.w.order = 0, n.sum)
debye1(x)
debye2(x)

Value

numeric/complex vector as z, or x, respectively.

Arguments

z

numeric or complex vector

s

complex number; current implementation is aimed at $s\in \{0,-1,\dots \}$

method

a string specifying the algorithm to be used.

logarithm

logical specified to return log(Li.(.)) instead of Li.(.)

is.log.z

logical; if TRUE, the specified z argument is really $w=\log (z)$; that is, we compute ${\mathrm{Li}}_s(\exp (w))$, and we typically have $w<0$, or equivalently, $z<1$.

is.logmlog

logical; if TRUE, the specified argument z is $lw=\log (-w)=\log (-\log (z))$ (where as above, $w=\log (z)$).

asymp.w.order

currently only default is implemented.

n.sum

for method="sum" only: the number of terms used.

x

numeric vector, may contain Inf, NA, and negative values.

Details

Almost entirely taken from https://en.wikipedia.org/wiki/Polylogarithm:

For integer values of the polylogarithm order, the following explicit expressions are obtained by repeated application of $z\frac{\partial }{\partial z}$ to ${\mathrm{Li}}_1(z)$: $${\mathrm{Li}}_1(z)=-\log (1-z),{\mathrm{Li}}_0(z)=\frac{z}{1-z},{\mathrm{Li}}_{-1}(z)=\frac{z}{(1-z{)}^2},{\mathrm{Li}}_{-2}(z)=\frac{z(1+z)}{(1-z{)}^3},$$ ${\mathrm{Li}}_{-3}(z)=\frac{z(1+4z+z^2)}{(1-z{)}^4}$, etc.

Accordingly, the polylogarithm reduces to a ratio of polynomials in z, and is therefore a rational function of z, for all nonpositive integer orders. The general case may be expressed as a finite sum:

$${\mathrm{Li}}_{-n}(z)={\left(z\frac{\partial }{\partial z}\right)}^n\frac{z}{1-z}=\sum _{k=0}^{n}k!S(n+1,k+1){\left(\frac{z}{1-z}\right)}^{k+1}(n=0,1,2,\dots ),$$ where $S(n,k)$ are the Stirling numbers of the second kind.

Equivalent formulae applicable to negative integer orders are (Wood 1992, § 6) ... $${\mathrm{Li}}_{-n}(z)=\frac{1}{(1-z{)}^{n+1}}\sum _{k=0}^{n-1}⟨\genfrac{}{}{0ex}{}{n}{k}⟩z^{n-k}=\frac{z\sum _{k=0}^{n-1}⟨\genfrac{}{}{0ex}{}{n}{k}⟩z^k}{(1-z{)}^{n+1}},(n=1,2,3,\dots ),$$ where $⟨\genfrac{}{}{0ex}{}{n}{k}⟩$ are the Eulerian numbers; see also Eulerian.

References

Wikipedia (2011) Polylogarithm, https://en.wikipedia.org/wiki/Polylogarithm.

Wood, D. C. (June 1992). The Computation of Polylogarithms. Technical Report 15-92. Canterbury, UK: University of Kent Computing Laboratory. https://www.cs.kent.ac.uk/pubs/1992/110/.

Apostol, T. M. (2010), "Polylogarithm", in the NIST Handbook of Mathematical Functions, https://dlmf.nist.gov/25.12

Lewin, L. (1981). Polylogarithms and Associated Functions. New York: North-Holland. ISBN 0-444-00550-1.

For Debye functions: Levin (1981) above, and
Wikipedia (2014) Debye function, https://en.wikipedia.org/wiki/Debye_function.

See Also

The polylogarithm is used in MLE for some Archimedean copulas; see emle;

The Debye functions are used for tau or rho computations of the Frank copula.

Examples

## The dilogarithm,  polylog(z, s = 2) = Li_2(.) -- mathmatically defined on C \ [1, Inf)
## so x -> 1 is a limit case:
polylog(z = 1, s = 2)
## in the limit, should be equal to
pi^2 / 6

## Default method uses  GSL's dilog():
rLi2 <- curve(polylog(x, 2), -5, 1, n= 1+ 6*64, col=2, lwd=2)
abline(c(0,1), h=0,v=0:1, lty=3, col="gray40")
## "sum" method gives the same for |z| < 1 and large number of terms:
ii <- which(abs(rLi2$x) < 1)
stopifnot(all.equal(rLi2$y[ii],
            polylog(rLi2$x[ii], 2, "sum", n.sum = 1e5),
          tolerance = 1e-15))


z1 <- c(0.95, 0.99, 0.995, 0.999, 0.9999)
L   <- polylog(         z1,  s=-3,method="negI-s-Euler") # close to Inf
LL  <- polylog(     log(z1), s=-3,method="negI-s-Euler",is.log.z=TRUE)
LLL <- polylog(log(-log(z1)),s=-3,method="negI-s-Euler",is.logmlog=TRUE)
all.equal(L, LL)
all.equal(L, LLL)

p.Li <- function(s.set, from = -2.6, to = 1/4, ylim = c(-1, 0.5),
                 colors = c("orange","brown", palette()), n = 201, ...)
{
    s.set <- sort(s.set, decreasing = TRUE)
    s <- s.set[1] # <_ for auto-ylab
    curve(polylog(x, s, method="negI-s-Stirling"), from, to,
          col=colors[1], ylim=ylim, n=n, ...)
    abline(h=0,v=0, col="gray")
    for(is in seq_along(s.set)[-1])
        curve(polylog(x, s=s.set[is], method="negI-s-Stirling"),
              add=TRUE, col = colors[is], n=n)
    s <- rev(s.set)
    legend("bottomright",paste("s =",s), col=colors[2-s], lty=1, bty="n")
}

## yellow is unbearable (on white):
palette(local({p <- palette(); p[p=="yellow"] <- "goldenrod"; p}))

## Wikipedia page plot (+/-):
p.Li(1:-3, ylim= c(-.8, 0.6), colors = c(2:4,6:7))

## and a bit more:
p.Li(1:-5)

## For the range we need it:
ccol <- c(NA,NA, rep(palette(),10))
p.Li(-1:-20, from=0, to=.99, colors=ccol, ylim = c(0, 10))
## log-y scale:
p.Li(-1:-20, from=0, to=.99, colors=ccol, ylim = c(.01, 1e7),
     log = "y", yaxt = "n")
if(require(sfsmisc)) eaxis(2) else axis(2)