3

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

Ben
  • 91,027
  • 3
  • 150
  • 376
Ratio
  • 31
  • 1
  • 2
    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
  • 2
    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

1 Answers1

1

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:

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.

Ben
  • 91,027
  • 3
  • 150
  • 376