Let $F(x;\theta)$ be a cumulative distribution function and $\beta>0$. I need to evaluate


but, for some values of $\theta$, R returns $F(x;\theta)=1$, and then I have some problems calculating $\rho$. The CDF is not exactly 1 at those parameter values. Is there a trick to obtain a stable calculation of $\rho$?

  • 2
    $\begingroup$ For most built-in functions, R supports options to (a) return the complementary CDF and/or (b) return its logarithm. Can you avail yourself of them? if not, then apply L'Hopital's Rule to estimate $\rho.$ $\endgroup$ – whuber Dec 18 '18 at 19:01
  • $\begingroup$ @whuber so, are you suggesting to calculate $F = 1 - (1-F)$? $\endgroup$ – Ratio Dec 18 '18 at 19:04
  • 2
    $\begingroup$ Sort of--but there will be algebraic cancellation. If you let $z=1-F,$ then to second order in $z,$ $$\rho = \frac{1}{\beta}z^{-1} + \frac{\beta-1}{2\beta} + \frac{\beta^2-1}{12\beta}z + O(z^2).$$ Thus, depending on the size of $\beta$ (you need $1/\beta$ and $z/\beta$ to be small), all you really need is a decent estimate of $z.$ $\endgroup$ – whuber Dec 18 '18 at 19:08

Let $W_x(\theta) \equiv \log F(x|\theta)$ and note that $W_x(\theta) \leqslant 0$ since it is a log-probability. (Since you have specified that $F(x|\theta)$ is close to one, this means that $W_x(\theta)$ is close to zero.) Then you can have:

$$\begin{equation} \begin{aligned} \log \rho &= \beta \log F(x;\theta) - \log ( F(x;\theta)-F(x;\theta)^{\beta+1}) \\[6pt] &= (\beta-1) \log F(x;\theta) - \log ( 1 - F(x;\theta)^\beta) \\[6pt] &= (\beta-1) \log F(x;\theta) - \log ( 1 - \exp(\beta \log F(x;\theta))) \\[6pt] &= (\beta-1) W_x(\theta) - \log ( 1 - \exp(\beta W_x(\theta))) \\[6pt] &= (\beta-1) W_x(\theta) - \text{log1mexp}(-\beta W_x(\theta)), \\[6pt] \end{aligned} \end{equation}$$

where $\text{log1mexp}(a) \equiv \log (1-e^{-a})$. Accurate calculation of expression hinges on accurate calculation of the $\text{log1mexp}$ function for argument values $a$ near zero. Accurate calculation of this function is considered in Machler (2012), and is implemented in the VGAM package in R. This is the required code:


#Function to compute log-rho
#Input W and beta
LOG_RHO <- function(W, beta) { (beta-1)*W - VGAM::log1mexp(beta*abs(W)); }

#Example with W close to zero (so F is close to one)
W    <- -10^(-8);
beta <- 100;

LOG_RHO(W, beta);
[1] 13.81551

exp(LOG_RHO(W, beta));
[1] 999999.5

You can see that this allows you to calculate your function even for values of $F(x|\theta)$ that are extremely close to one. This should serve your purposes well. If you would like to know more about computational issues relating to log-probabilties, see related questions here and here.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.