\documentclass[article,nojss]{jss} %% NOTA BENE: More definitions --> further down %%%%%%%%%%%% % \author{Martin M\"achler \\ ETH Zurich% \\ May 2011 {\tiny (\LaTeX'ed \today)}%---- for now } \title{Numerically Stable Frank Copula Functions via Multiprecision: \proglang{R} Package \pkg{Rmpfr}} % \def\mythanks{a version of this paper, for \pkg{nacopula} 0.4\_4, has been published % in JSS, \url{http://www.jstatsoft.org/v39/i09}.} %% for pretty printing and a nice hypersummary also set: \Plainauthor{Martin M\"achler} %% comma-separated \Plaintitle{Numerically Stable Frank Copula Functions via Multiprecision: R Package 'Rmpfr'} \Shorttitle{Numerically Stable Frank via Multiprecision in R} % %\VignetteIndexEntry{Franknacopula} %\VignetteDepends{nacopula} %\VignetteDepends{Rmpfr} \SweaveOpts{engine=R,eps=FALSE,pdf=TRUE,width=7,height=5,strip.white=true,keep.source=TRUE} %% an abstract and keywords \Abstract{ The package \pkg{nacopula} .... ... Archimedean copulas The package \pkg{Rmpfr} .... } \Keywords{Archimedean copulas, Frank Copula, Multiple Precision, R, MPFR, Rmpfr} %% at least one keyword must be supplied %% publication information %% NOTE: Typically, this can be left commented and will be filled out by the technical editor %% \Volume{13} %% \Issue{9} %% \Month{September} %% \Year{2004} %% \Submitdate{2004-09-29} %% \Acceptdate{2004-09-29} %% The address of (at least) one author should be given %% in the following format: \Address{ Martin M\"achler\\ Seminar f\"ur Statistik, HG G~16\\ ETH Zurich\\ 8092 Zurich, Switzerland\\ E-mail: \email{maechler@stat.math.ethz.ch}\\ URL: \url{http://stat.ethz.ch/people/maechler} } %% It is also possible to add a telephone and fax number %% before the e-mail in the following format: %% Telephone: +43/1/31336-5053 %% Fax: +43/1/31336-734 %% for those who use Sweave please include the following line (with % symbols): %% MM: this is "substituted" by jss.cls: %% need no \usepackage{Sweave.sty} \usepackage[american]{babel}%for American English \usepackage{amsmath}%sophisticated mathematical formulas with amstex (includes \text{}) \usepackage{mathtools}%fix amsmath deficiencies \usepackage{amssymb}%sophisticated mathematical symbols with amstex (includes \mathbb{}) % \usepackage{amsthm}%theorem environments \usepackage{bm}%for bold math symbols: \bm (= bold math) \usepackage{enumitem}%for automatic numbering of new enumerate environments % This is already in jss above -- but withOUT the fontsize=\small part !! \DefineVerbatimEnvironment{Sinput}{Verbatim}{fontsize=\small,fontshape=sl} \DefineVerbatimEnvironment{Soutput}{Verbatim}{fontsize=\small} \DefineVerbatimEnvironment{Scode}{Verbatim}{fontsize=\small,fontshape=sl} %% makes space between Sinput and Soutput smaller: \fvset{listparameters={\setlength{\topsep}{0pt}}}% !! quite an effect! %% ??? FIXME but it also influences all other lists, itemize, ... ???? FIXME %% \setkeys{Gin}{width=\textwidth}% Sweave.sty has {width=0.8\textwidth} \newcommand*{\R}{\proglang{R}}%{\textsf{R}} %% Marius' commands \newcommand*{\eps}{\varepsilon} %NON-STANDARD{bbm}:\newcommand*{\I}{\mathbbm{1}} % \newcommand*{\I}{\mathbf{1}} % \newcommand*{\IN}{\mathbb{N}} % \newcommand*{\IK}{\mathbb{K}} % \newcommand*{\IR}{\mathbb{R}} % \newcommand*{\IC}{\mathbb{C}} % \newcommand*{\IP}{\Prob} % \newcommand*{\IE}{\E} % \newcommand*{\V}{\operatorname*{V}} %- \abs{ab} --> | ab | ``absolut Betrag'' \newcommand{\abs}[1] {\left| #1 \right|} \newcommand*{\Li}[1]{\mathrm{Li}_{#1}}% polylog / polylogarithm % \renewcommand*{\S}{\operatorname*{S}} % \newcommand*{\tS}{\operatorname*{\tilde{S}}} % \newcommand*{\ran}{\operatorname*{ran}} %\newcommand*{\sgn}{\operatorname*{sgn}} \DeclareMathOperator{\sign}{sign} \newcommand*{\psiD}{\psi^\prime} \newcommand*{\psii}{{\psi^{-1}}} \newcommand*{\psiis}[1]{{\psi_{#1}^{-1}}} % \renewcommand*{\L}{\mathcal{L}} % \newcommand*{\Li}{\mathcal{L}^{-1}} % \newcommand*{\LS}{\mathcal{LS}} % \newcommand*{\LSi}{\LS^{-1}} \newcommand{\tr}{\ensuremath{^\mathsf{T}}}% or ^{\intercal} \renewcommand*{\O}{\mathcal{O}} % \newcommand*{\Geo}{\operatorname*{Geo}} % \newcommand*{\Exp}{\operatorname*{Exp}} % \newcommand*{\Sibuya}{\operatorname*{Sibuya}} % \newcommand*{\Log}{\operatorname*{Log}} % \newcommand*{\U}{\operatorname*{U}} % \newcommand*{\B}{\operatorname*{B}} % \newcommand*{\NB}{\operatorname*{NB}} % \newcommand*{\N}{\operatorname*{N}} \DeclareMathOperator{\var}{var} \DeclareMathOperator{\Var}{Var} \DeclareMathOperator{\Cov}{Cov} \DeclareMathOperator{\cov}{cov} \DeclareMathOperator{\Cor}{Corr} \DeclareMathOperator{\cor}{corr} % \newcommand*{\Var}{\operatorname*{Var}} % \newcommand*{\Cov}{\operatorname*{Cov}} % \newcommand*{\Cor}{\operatorname*{Cor}} \hyphenation{Ar-chi-me-dean} %% journal specific aliases \newcommand*{\setcapwidth}[1]{} %% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{document} %% include your article here, just as usual %% Note that you should use the \pkg{}, \proglang{} and \code{} commands. % \section[About Java]{About \proglang{Java}} %% Note: If there is markup in \(sub)section, then it has to be escape as above. %% Note: These are explained in '?RweaveLatex' : <>= op.orig <- options(width = 75, SweaveHooks= list(fig=function() par(mar=c(5.1, 4.1, 1.1, 2.1))), useFancyQuotes = FALSE, ## for JSS, but otherwise MM does not like it: ## prompt="R> ", continue=" ")# 2 (or 3) blanks: use same length as 'prompt' copDDir <- system.file('doc', package='nacopula') copSrcDDir <- { if(Sys.getenv("USER") == "maechler") '~/R/Pkgs/nacopula/inst/doc' else file.path(dirname(tempdir()),"nacopula-saveDir") } if(!file.exists(copSrcDDir)) stopifnot(dir.create(copSrcDDir)) Sys.setenv(LANGUAGE = "en") if(.Platform$OS.type != "windows") Sys.setlocale("LC_MESSAGES","C") @ %\section[Introduction]{Introduction \small~\footnote{\mythanks}} %\section{Introduction} \section{The diagonal density of Frank's copula} The ``diagonal density'' of a copula $C(\cdot)$ is the density of $\max(U_1,U_2,\dots,U_d)$, where \linebreak[4] $\bm U = (U_1,U_2,\dots,U_d)\tr \sim C$. The (cumulative) distribution function, by definition, \begin{equation} F^D(u) := P\left[ \max(U_1,U_2,\dots,U_d) \le u \right] = P\left[U_1\le u, U_2,\le u, \dots, U_d\le u\right] = C(u,u,\dots,u), \label{eq:diagF} \end{equation} evaluates the copula only on the diagonal and is therefore called \emph{``the diagonal of $C$''}. Its density $f^D(u) := \frac{d}{du}F^D(u)$, is therefore called the diagonal density of $C$. For Archimedean copulas, i.e., where \begin{align} C(\bm{u})=C(\bm{u};\psi)=\psi(\psii(u_1)+\dots+\psii(u_d)),\ \bm{u}\in[0,1]^d,\label{ac} \end{align} the diagonal density is \begin{eqnarray} f^D(u) &=& \frac{d}{du}F^D(u) = \frac{d}{du}\psi\left(\sum\nolimits_{j=1}^d \psii(u)\right) = \frac{d}{du}\psi(d\cdot \psii(u)) = \nonumber\\ &=& \psi^\prime\!\left(d\cdot \psii(u)\right) \cdot d \cdot \frac{d}{du} \psii(u)\nonumber\\ &=& d\cdot \psi^\prime\!\left(d\cdot \psii(u)\right) \cdot \left[\psii\right]^\prime(u). \label{eq:fD} \end{eqnarray} For this reason, the \pkg{nacopula} package's \code{dDiag()} function for computing the diagonal density $f^D(u)$ makes use of the following %% ==> lines 261 ff from ~/R/Pkgs/nacopula/R/estimation.R %% ================================ <>= nacopula:::dDiagA <>= writeLines(head(capture.output(print( <> )), -1)) @ where the three functions \begin{align} \mathtt{psiDabs(t,thet)} &= \bigl|\psi_\theta^\prime(t)\bigr|,\label{def:psiDabs}\\ \mathtt{psiInv(u,thet)} &= \psi_\theta^{-1}(u), \ \mathrm{and}\label{def:psiInv}\\ \mathtt{psiInvD1abs(u,thet)} &= \bigl|[\psi_\theta^{-1}]^\prime(u) \bigr| \label{def:psiInvD1abs} \end{align} are all provided by the slots of the corresponding Archimedean copula family. For the following explorations, we need a definition of \code{dDiagA} which is more flexible as it does not work with the copula family object but gets the three functions as arguments, <>= dDiagA <- function(u, th, d, psiInv, psiDabs, psiInvD1abs, log = FALSE) { stopifnot(is.finite(th), d > 0, is.function(psiInv), is.function(psiDabs), is.function(psiInvD1abs)) if(log) { log(d) + psiDabs(d*psiInv(u,th), th, log = TRUE) + psiInvD1abs(u, th, log = TRUE) } else { d* psiDabs(d*psiInv(u,th), th) * psiInvD1abs(u,th) } } @ Now, for the Frank copula (see Hofert and M\"achler (2011)),% FIXME \cite{Hofert:Maechler:201....} \begin{eqnarray} \psi_\theta(t) &=& - \frac{1}{\theta}\log\left(1 - (1 - e^{-\theta})e^{-t}\right), \label{eq:psiFrank} \ \ \ \theta > 0, \ \ \mathrm{\ \ hence,} \\ \psi_\theta^{-1}(u) &=& - \log\left(\frac{e^{-u\theta} - 1}{e^{-\theta} - 1}\right), \label{eq:psiInvFrank} \ \ \mathrm{\ \ and} \\ %% psiDabs(t,thet) (-1)^k {\psi_\theta}^{(k)}(t) = \bigl|{\psi_\theta}^{(k)}(t)\bigr| &=& \frac{1}{\theta} \Li{k-1}((1-e^{-\theta})e^{-t}), \label{eq:psiDabs} \ \ \ \ \ \ \mathrm{and} \\ %% psiInvD1abs.1(u,theta) \bigl|[\psi_\theta^{-1}]^\prime(u) \bigr| &=& \theta/(e^{u\cdot\theta} - 1), \label{eq:psiInvD1} \end{eqnarray} where $\Li{s}(z)$ is the \emph{polylogarithm of order} $s$ at $z$, defined as (analytic continuation of) $\sum_{k=1}^\infty \frac{z^k}{k^s}$. When $s = -n$, $n\in\mathbb{N}$, one has \begin{equation} \label{eq:Li-n} \Li{-n}(z) = \left(z \cdot \frac{\partial}{\partial z} \right)^n \frac{z}{1-z}, \end{equation} and we note that here, only the first derivative is needed, $ - {\psi_\theta}^\prime(t) = \left|{\psi_\theta}^\prime(t)\right|$, and hence only \begin{equation} \mathtt{polylog(z,\: s=0)} = \Li{0}(z) = z / (1 - z). \label{eq:polylog-0} \end{equation} First note that numerically, $e^{-a} - 1$ suffers from cancellation when $0 < a \ll 1 $, and the \R\ (and C) function \code{expm1(-a)} is advisably used instead of \code{exp(-a) - 1}. For this reason, in \pkg{nacopula}, I had replaced the original \code{psiInv.0(u, th)} for $\psi_\theta^{-1}(u)$ by \code{psiInv.1()}, making use of \code{expm1()}. These and the derivative (\ref{eq:psiDabs}) originally were <>= psiInv.0 <- function(u,theta) -log( (exp(-theta*u)-1) / (exp(-theta)-1) ) psiInv.1 <- function(u,theta) -log(expm1(-u*theta) / expm1(-theta)) psiInvD1abs.1 <- function(u, theta, log = FALSE) if(log) log(theta)-log(expm1(u*theta)) else theta/expm1(u*theta) @ and the general $k$-th derivative (\ref{eq:psiInvD1}), simplified, for $k= 1$ (\code{degree = 1}), <>= require("nacopula")# for polylog() psiDabs.1 <- function(t, theta, log=FALSE) { p <- -expm1(-theta) Li. <- polylog(log(p) - t, s = 0, log=log, method="negI-s-Eulerian", is.log.z=TRUE) if(log) Li. - log(theta) else Li. / theta } @ where we however now assume that the \code{polylog()} function for \code{s=0} would basically use the direct formula (\ref{eq:polylog-0}), such that we use <>= psiDabs.2 <- function(t, theta, log=FALSE) { w <- log(-expm1(-theta)) - t Li. <- if(log) w - log(-expm1(w)) else -exp(w)/expm1(w) if(log) Li. - log(theta) else Li. / theta } @ \section{Computing the ``diagonal MLE''} The most important use of the diagonal density is to compute the ``diagonal maximum likelihood estimator'', \code{dmle}, ${\hat\theta}^D$ which for a sample of observations $\bm u_1,\bm u_2,\ldots, \bm u_n$, ($\bm u_i \in [0,1]^d$) is defined as minimizer of the negative log-likelihood, \begin{align} {\hat\theta}^D &= \arg\min_\theta - l(\theta; \bm u_1,\ldots,\bm u_n),\ \ \mathrm{where}\\ l(\theta; \bm u_1,\ldots,\bm u_n) &= \sum_{i=1}^n \log f^D(\tilde u_i)\ \ \mathrm{and}\\ \tilde{u_i} &= \max_{j=1,\dots,d} u_{i,j}. \label{eq:dmle} \end{align} In our exploration of the \code{dmle} estimator, we found cases with numerical problems, at first already in evaluating the logarithm of the diagonal density $\log f^D = \log f^D_\theta(u)$\code{= dDiag(u, theta, *, log=TRUE)} for non-small $\theta$ and ``large'' $u$, i.e., $u \approx 1$: <>= curve(dDiagA(x, th = 38, d = 2, psiInv = psiInv.0, psiDabs=psiDabs.1, psiInvD1abs=psiInvD1abs.1, log = TRUE), ylab = "dDiagA(x, theta= 38, *, log=TRUE)", 0, 1, col = 4, n = 1000) ## and using the slightly better psiInv.1 does not help here: curve(dDiagA(x, th = 38, d = 2, psiInv = psiInv.1, psiDabs=psiDabs.2, psiInvD1abs=psiInvD1abs.1, log = TRUE), add = TRUE, col = 2, n=1000) legend("bottom", c("psiInv.0()","psiInv.1()"),col=c(4,2), lty=1, bty="n") @ However, it's not hard to see that indeed our initial computation of Frank's $\psii$, i.e, (\ref{eq:psiInvFrank}), $\psiis{\theta}(u) = -\log\left(\frac{1-e^{-u\theta}}{1-e^{-\theta}}\right)$, suffers from ``division cancellation'' for ``large'' $\theta$ ($\theta=38$ in ex.) when computed directly with \code{psiInv.0()}, see above, and that the improvement of using \code{expm1(-t)} instead of \code{exp(-t) - 1}, as used in \code{psiInv.1()}, see above, helps only for $t\approx 0$. However, we can rewrite $\psii$ as \begin{equation} \label{eq:psii:log:1} \psiis{\theta}(u) = -\log\left(1 -\frac{e^{-u\theta}-e^{-\theta}}{1-e^{-\theta}}\right), \end{equation} which when using \code{log1p(e)} for $\log(1+e)$ is much better numerically: <>= psiInv.2 <- function(u,theta) -log1p((exp(-u*theta)-exp(-theta)) / expm1(-theta)) curve(dDiagA(x, th = 38, d = 2, psiInv = psiInv.2, psiDabs=psiDabs.2, psiInvD1abs=psiInvD1abs.1, log = TRUE), ylab = "dDiagA(x, theta= 38, *, log=TRUE)", 0, 1, col = 4, n = 1000) ## previously curve(dDiagA(x, th = 38, d = 2, psiInv = psiInv.1, psiDabs=psiDabs.2, psiInvD1abs=psiInvD1abs.1, log = TRUE), add = TRUE, col = "darkgray", lwd=2, lty=3, n=1000) @ Unfortunately, this is not enough to get numerically stable evaluations of the negative log-likelihood $l()$: <>= d <- 5 (theta <- copFrank@tauInv(tau = 0.75)) cop <- onacopulaL("Frank", list(theta, 1:d)) set.seed(1); for(l in 1:4) U <- rnacopula(n = 100, cop) U. <- sort(apply(U, 1, max)) # build the max <>= mlogL <- function(theta) -sum(dDiagA(U., theta, d=d, psiInv = psiInv.2, psiDabs=psiDabs.2, psiInvD1abs=psiInvD1abs.1, log = TRUE)) @ Now, plot this negative log likelihood function in an interval $\theta$, close to what is proposed \pkg{nacopula}'s \code{initOpt()} function, defining a utility function that we'll reuse later: <>= p.mlogL <- function(th, mlogL, col= "red2", lwd = 1, lty = 1, add= FALSE) { stopifnot(is.numeric(th), is.function(mlogL)) nll <- vapply(th, mlogL, 0.) if(add) lines(nll ~ th, col=col, lwd=lwd, lty=lty) else plot(nll ~ th, xlab=expression(theta), ylab = expression(- logLik(theta, .)), type = "l", col=col, lwd=lwd, lty=lty) invisible(nll) # return invisibly } thet <- seq(11, 99, by = 1/4) p.mlogL(thet, mlogL) require("Rmpfr")## compute the same with *high* accuracy ... ## using three different precisions: MPrecBits <- c(160, 128, 96) mkNm <- function(bits) sprintf("%03d.bits", bits) ## As it takes a while seconds, cache the result: fnam <- sprintf("mlogL_mpfr_%s.rda", Sys.info()[["machine"]]) if (!file.exists(fn <- file.path(copDDir,fnam))) { print(system.time( nllMP <- lapply(MPrecBits, function(pBit) { nlM <- thM <- mpfr(thet, precBits = pBit) ## (vapply() does not work for "Rmpfr":) for(i in seq_along(thet)) nlM[i] <- mlogL(thM[i]) nlM }) )) ## 91.226 0.013 91.506 [nb-mm icore 5] names(nllMP) <- mkNm(MPrecBits) save(nllMP, file = file.path(copSrcDDir, fnam)) } else load(fn) colB <- c("blue3","violetred4","tan3") ltyB <- c(5:3) lwdB <- c(2,2,2) for(i in seq_along(nllMP)) { lines(thet, as.numeric(nllMP[[i]]), col=colB[i], lty = ltyB[i], lwd = lwdB[i]) } leg <- c("double prec.", sprintf("mpfr(*, precBits = %d)", MPrecBits)) legend("top", leg, col= c("red3",colB), lty=c(1, ltyB), lwd=c(1,lwdB), bty="n") @ So, clearly, high-precision computations can solve the numerical problems, if the precision is high enough. E.g., for $\theta = 100$, it needs more than 128 bits precision. Let's look at the phenomenon in more details now. The flesh in the \code{mlogL()} computation is (up to the constant $\log(d)$, $d=5$), only the sum of the two terms <>= op <- options(prompt = " ") <>= psiDabs(d*psiInv(u,th), th, log = TRUE) + psiInvD1abs(u, th, log = TRUE) @ currently, with the three functions <<3f, eval=false>>= psiInv = psiInv.2; psiDabs = psiDabs.2; psiInvD1abs = psiInvD1abs.1 <>= options(op) @ where we have already tried to ensure that the \code{psiInv()} function is ok, but now can confirm it, e.g., for $\theta = 50$. Further note, that using high-precision arithmetic, we can also ``partially afford'' to use the simplistic \code{psiInv.0()} function instead the more stable \code{psiInv.2()} one: <>= stopifnot(all.equal(psiInv.2(U., 50 ), psiInv.2(U., mpfr(50, 96))), all.equal(psiInv.0(U., mpfr(50, 200)), pI.U <- psiInv.2(U., mpfr(50, 200)), tol=1e-50) ) @ However, we can observe dramatic differences in \code{psiDabs.2()} ($= |\psi^\prime(.)|$): <>= psD.n <- psiDabs.2(as.numeric(pI.U), 40) psD.M <- as.numeric(psiDabs.2(pI.U, mpfr(40, 200))) matplot(U., cbind(psD.n, psD.M), type="l", log="y") legend("top", c("double prec.", "mpfr(*, precBits = 200)"), col= 1:2, lty=1:2, bty="n") @ where we see a very large difference (note the \emph{$\log$} scale!) for $u \approx 1$, i.e., for very small \linebreak[3] \code{pI.U= psiInv.2(U., theta)}$ = \psi_\theta^{-1}(u)$. <>= u0 <- 2^-(100:1) psD.n <- psiDabs.2(u0, 40) psD.M <- as.numeric(psiDabs.2(u0, mpfr(40, 200))) matplot(u0, cbind(psD.n, psD.M), type="l", log="xy") legend("top", c("double prec.", "mpfr(*, precBits = 200)"), col= 1:2, lty=1:2, bty="n") @ Further investigation shows that the culprit is really the use of \code{log(-expm1(-theta))} inside \code{psiDabs.2()} which underflows for large theta, and hence should be replaced by the generally accurate %% as shown by Maechler (2011) --- FIXME: finally right that small paper!! %% ========================= ~/R/MM/NUMERICS/log1-exp.R <>= log1mexpm <- function(a) { stopifnot(a >= 0) r <- a tst <- a <= log(2) r[ tst] <- log(-expm1(-a[ tst])) r[!tst] <- log1p(-exp(-a[!tst])) r } @ so that we rather compute $|\psi^\prime(.)|$ via <>= psiDabs.3 <- function(t, theta, log=FALSE) { w <- log1mexpm(theta) - t Li. <- if(log) w - log1mexpm(-w) else -exp(w)/expm1(w) if(log) Li. - log(theta) else Li. / theta } @ Does this already solve the ``diagonal likelihood'' problem? Let's see, using a <>= p.mlogL(th = seq(11, 99, by = 1/4), mlogL = (mlogL2 <- function(theta) -sum(dDiagA(U., theta, d=d, psiInv = psiInv.2, psiDabs = psiDabs.3, psiInvD1abs=psiInvD1abs.1, log = TRUE))), lwd = 2) @ Yes, indeed, using a numerically stable version for \code{psiDabs()} did solve the numerical problem of computing, the ``diagonal likelihood'', and hence the \code{dmle()} (``diagonal MLE''). Well, if we really insist, there are more problems, but probably not really practical: <>= thet <- 9:1000 nll <- p.mlogL(thet, mlogL = mlogL2, lwd=2) (th0 <- thet[i0 <- max(which(is.finite(nll)))]) abline(v = th0, col="red2", lty="15", lwd=2) @ where we see that for $\theta > \Sexpr{th0}$, \code{mlogL()} is not finite, e.g., <>= dDiagA(0.999, 715, d = d, psiInv = psiInv.2, psiDabs = psiDabs.3, psiInvD1abs = psiInvD1abs.1, log = TRUE) @ which after closer inspection is from the \code{psiInvD1abs(u, th, log = TRUE)} part of \code{dDiagA()} and that, see (\ref{def:psiInvD1abs}), uses \code{log(theta)-log(expm1(u*theta))}, where clearly already \code{expm1(u*theta)) = expm1(0.999 * 715)} numerically overflows to \code{Inf}: <>= psiInvD1abs.1(0.999, th = 715, log=TRUE) @ However, as $\log(\mathrm{expm1}(y)) = \log(e^y -1) = \log(e^y(1 - e^{-y})) = y + \log(1 - e^{-y})$, the ``numerical stable'' solution is to replace \code{log(expm1(y))} by \code{y + log1mexpm(y)}, such that we will use <>= psiInvD1abs.2 <- function(u, theta, log = FALSE) if(log) log(theta)- {y <- u*theta; y + log1mexpm(y)} else theta/expm1(u*theta) @ Unfortunately, this improves the situation for large $\theta$ only slightly: <>= plot(nll ~ thet, xlab=expression(theta), ylab = expression(- logLik(theta, .)), type = "l", col="red2", lwd=2) abline(v = th0, col="red2", lty="15", lwd=2) nll3 <- p.mlogL(thet, mlogL = function(theta) -sum(dDiagA(U., theta, d=d, psiInv= psiInv.2, psiDabs= psiDabs.3, psiInvD1abs = psiInvD1abs.2, log = TRUE)), col = "blue2", lwd=3, lty=2, add = TRUE) nll3[thet == 800] @ where we can see, that this time, the numerical overflow to $-\infty$ happens in \code{psiDabs(d*psiInv(u,th), th, log = TRUE)}, as \code{psiInv(u,th) = psiInv(0.999, th=800)}$=$\Sexpr{format(psiInv.2(u=0.999,th=800))} underflows to zero --- by necessity: the correct value is smaller than the smallest representable double precision number: <>= pI <- psiInv.2(u=0.999, th= mpfr(800, 200)) cat(sapply(list(pI, .Machine$double.xmin), format, digits = 7), "\n") @ The only solution here is to allow passing the \emph{logarithm} $\log(t)$ instead of $t$ to \code{psiDabs()}, i.e., we start computing \code{log(psiInv(.))} directly (evading the underflow to 0 there). \dots \dots to be continued. \section{Session Information} <>= toLatex(sessionInfo()) <>= options(op.orig) @ \section{Conclusion} %\bibliography{nacopula}% Rmpfr \end{document}