# Evaluating problematic function when cdf is close to one?

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

$$\rho=\frac{F(x;\theta)^\beta}{F(x;\theta)-F(x;\theta)^{\beta+1}},$$

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$$?

• 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.$ – whuber Dec 18 '18 at 19:01
• @whuber so, are you suggesting to calculate $F = 1 - (1-F)$? – Ratio Dec 18 '18 at 19:04
• 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.$ – 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{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}

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:

library(VGAM);

#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.