| polylog {copula} | R Documentation |
Polylogarithm \mathrm{Li_s(z)} and Debye Functions
Description
Compute the polylogarithm function \mathrm{Li_s(z)},
initially defined as the power series,
\mathrm{Li}_s(z) = \sum_{k=1}^\infty {z^k \over 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(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)
Arguments
z |
numeric or complex vector |
s |
complex number; current implementation is aimed at
|
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 |
is.logmlog |
logical; if TRUE, the specified argument |
asymp.w.order |
currently only default is implemented. |
n.sum |
for |
x |
numeric vector, may contain |
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) = {z \over 1-z}, \ \
\mathrm{Li}_{-1}(z) = {z \over (1-z)^2}, \ \
\mathrm{Li}_{-2}(z) = {z \,(1+z) \over (1-z)^3},
\mathrm{Li}_{-3}(z) = {z \,(1+4z+z^2) \over (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 \,{\partial \over \partial z} \right)^n \frac{z}{1-z}
= \sum_{k=0}^n k! \,S(n+1,k+1) \left({z \over {1-z}} \right)^{k+1}
\ \ (n=0,1,2,\ldots),
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) = {1 \over (1-z)^{n+1}} \sum_{k=0}^{n-1}
\left\langle {n \atop k} \right\rangle z^{n-k} =
\frac{z \sum_{k=0}^{n-1} \left\langle {n \atop k} \right
\rangle z^k}{(1-z)^{n+1}},
\qquad (n=1,2,3,\ldots) ~,
where \left\langle {n \atop k} \right\rangle are the
Eulerian numbers; see also Eulerian.
Value
numeric/complex vector as z, or x, respectively.
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)