Inverse transform sampling and ambiguous Intervals

Question

Let $F_i:\mathbb R\to[0,1]$ be a distribution function$^1$ and $$F_i^{-1}(t):=\inf\left\{x\in\mathbb R:F_i(x)\ge t\right\}\;\;\;\text{for }t\in[0,1].$$

I've got a computer program where only $F_i^{-1}$ (not $F_i$) is directly available. Assume $t_1\in[0,1]$ and $x:=F_1^{-1}(t_1)\in F_2^{-1}([0,1])$. I need to compute a $t_2\in[0,1]$ such that $x=F_2^{-1}(t_2)$. How can I do this?

First of all, we know that $$\left\{t_2\in[0,1]:F_2^{-1}(t_2)=x\right\}=\begin{cases}[F_2(x-),F_2(x)]&\text{, if }F_2\text{ is continuous at }x\text{ or }\forall y<x:F_2(y)<F_2(x-)\\(F_2(x-),F_2(x)]&\text{, otherwise}.\end{cases}$$

Now, I've read (here in section 6.1) the following, but can't really make sense of it: Assuming that $F_2^{-1}([a,b])=\{x\}$ for some $0\le a\le b\le 1$ we can find $t_2$ by sampling $u$ with uniform distribution on $[0,1]$ and set $t_2:=a+(b-a)u$. Why does this work? And if it works, how can we apply it in practice (i.e. how do we find $a,b$)?

Remark: I think this book page (above Example 3.31) is related.

---

$^1$ i.e. $F_i$ is right-continuous and nondecreasing with $F(-\infty):=\lim_{x\to-\infty}F(x)=0$ and $F(\infty):=\lim_{x\to\infty}F(x)=1$.

Answer 1

The bisection method is guaranteed to work even for such discontinuous $F^{-1},$ provided it is suitably implemented.

Here is pseudocode (that actually works in R):

function(f, x, tol=1e-8, ...) {
  u <- 1
  l <- 0
  repeat {
    m <- (u + l) / 2
    if (f(m, ...) - x <= 0) l <- m else u <- m
    if (u - l <= tol) break # (See the end of this post for a better test)
  }
  return(u)
}

The arguments are (1) the name of $F^{-1},$ (2) the value of $x,$ (3) a positive error tolerance (the result will be accurate to this amount), and (4) any other arguments that need to be passed to $F^{-1}.$ I will refer to this function as findroot.

Before proving this works, let's look at how it might be used, again using R.

> findroot(qpois, 2, tol=0, lambda=2)
[1] 0.6766764

qpois is the Poisson percentile function. Thus, with $\lambda=2$ we hope that

$$0.6766764 = e^{-\lambda}(1 + \lambda + \lambda^2/2!)$$

and indeed that's the case. This figure plots part of $F^{-1},$ showing $x$ as a horizontal dashed line and the solution as a vertical red line:

Let's turn to proving this works. Let $\epsilon \ge 0$ be the tolerance. Consider the proposition

$$\mathcal{P}_{x}(l,u):\ F^{-1}(l) - x \le 0 \le F^{-1}(u) - x\quad \text{ and }\quad u^\prime > u \implies F^{-1}(u^\prime) - x \gt 0. $$

If we take the values of $F^{-1}$ at any number greater than $1$ to be $\infty,$ then $\mathcal{P}_{x}(0,1)$ is true. Assuming hypothetically $\mathcal{P}_{x}(l,u)$ at the beginning of the loop, note that $u$ will be decreased to $u^\prime$ only when $F^{-1}(u^\prime) - x \gt 0$ and in any event $F^{-1}$ changes sign between the new $l$ and new $u.$ Thus, $\mathcal{P}_{x}(l,u)$ remains true at the end of the loop. After exiting, $u$ and $l$ are within $\epsilon$ of each other and $\mathcal{P}_{x}(l,u)$ remains true (by induction). Thus, the value $t = u$ returned by findroot enjoys two properties:

$$F^{-1}(t-\epsilon)-x \le 0 \lt F^{-1}(t) - x.$$

That's what it means for $t$ to be within $\epsilon$ of a solution to $x = F^{-1}(t),$ QED.

Notice that after $n$ iterations of the loop, the difference $u-l = 2^{-n}.$ Therefore this procedure terminates after $\lceil -\log_2 \epsilon \rceil$ iterations. That's a reasonably sparing use of calls to $F^{-1}.$

---

In a practical application, the test u - l <= tol is too naive about floating-point roundoff error: if tol is very small (but still positive), this condition might never hold. One way to guarantee termination is to set an upper limit on the number of iterations; $52$ will be fine for double-precision arithmetic. A slightly more flexible solution in R uses zapsmall, as in

    if (zapsmall(c(u - l, 1))[1] <= tol) break

When $u-l$ is indistinguishable from $0$ compared to $1,$ it is set to $0,$ guaranteeing termination of the loop.