14
$\begingroup$

Suppose the cumulative distribution function $F$ is given but not invertible to use the inverse transform sampling technique (to compute $X=F^{-1}(Y)$). Do we have other alternative methods? I would appreciate to know the name of all possible methods...

$\endgroup$
2
  • $\begingroup$ Could you clarify in your question what you mean by "not invertible"? Is it "the inverse does not exist everywhere" or is it "the inverse is not easily computed" $\endgroup$
    – Xi'an
    Aug 10, 2020 at 6:49
  • 1
    $\begingroup$ First of all thanks for your answer. Actually i am interested to learn about existing approaches to deal with the CDF, which is either not invertible, or hard to invert. Thankfully, i got different answers about both situations, as well as the name of some approaches, about which i can search more. I got more than what i wanted to know. $\endgroup$
    – Rob
    Aug 10, 2020 at 8:42

3 Answers 3

5
$\begingroup$

In low dimensions a good alternative is to use rejection sampling from the pdf $f_X$ (in high dimensions this becomes very inefficient).

Say $f_X$ is your pdf for some random variable $X$, which you want to sample from in the interval $I=[x_\mathrm{min}, x_\mathrm{max}]$. Then you can draw samples $x_i$ uniformly from $I$ and accept/reject them with the probability $f_X(x_i)$, i.e. you draw another uniformly distributed random number $u_i\in[0, \max(f_X(x)|x \in I)]$ and if $u_i \lt f_X(x_i)$ you accept that sample point otherwise you reject it.

$\endgroup$
22
$\begingroup$

The inverse cdf method operates even when the cdf is not invertible, using the generalised inverse$$F^-(u)=\sup\{x;\ F(x)\le u\}\tag{1}$$ which is always defined for $u\in(0,1)$.

When there is no solution in $x$ to the equation $$F(x)=u$$ it means that $F$ has a jump between a value less than $u$, $u-\epsilon$ and a value more than $u$, $u+\eta$. Hence the distribution has a point mass with mass $\eta+\epsilon$ at a point $x_0$, with $$F(x_0^{-})=u-\epsilon\quad\text{and}\quad F(x_0)=u+\eta$$ In that case (1) leads to$$F^-(u)=x_0$$

Similarly, if the equation $$F(x)=u$$ has an infinite number of solutions, say $x\in [a,b)$, it means that the cdf is constant over this interval and hence that all values in (a,b) have zero probability to occur. In that case, $$F^-(u)=b$$ [which is of course a convention since the event $U=F(a)$ has probability zero to occur.

$\endgroup$
6
  • 11
    $\begingroup$ +1. This solution is illustrated in my post at stats.stackexchange.com/a/18439/919. I found it, along with a closely related post of yours, with a site search on inverse CDF jump. Another nice explanation is given by vqv at stats.stackexchange.com/a/5607/919. $\endgroup$
    – whuber
    Aug 8, 2020 at 16:46
  • 1
    $\begingroup$ Nice answer (+1). $\endgroup$
    – BruceET
    Aug 8, 2020 at 16:52
  • $\begingroup$ @whuber I'm trying to understand the proof... by $F(x_0^-)$ do we mean $F^-(x_0)$? $\endgroup$
    – Jake
    Apr 4, 2022 at 22:40
  • 1
    $\begingroup$ No, by $x_0^-$ I mean the limit inf of the $x$'s less than $x_0$. $\endgroup$
    – Xi'an
    Apr 5, 2022 at 5:08
  • $\begingroup$ @Xi'an How would you write this notation mathematically? I haven't seen it before. is it something like $F(x_0^-)=F(\text{inf}\{x:x<x_0 \})$. Where does the limit part come in? $\endgroup$
    – Jake
    Apr 5, 2022 at 18:40
1
$\begingroup$

Methods can be quite different depending on the distribution you want to simulate. There many good books on simulation methods and lots of info about simulation on the Internet. Methods often exploit relationships among distributions (such as Poisson, exponential, gamma, chi-squared, F, beta---see discussion in comments). Ultimately, almost all current computer simulation ultimately uses standard uniform output from a pseudo-random generator.

Sometimes, methods depend on the current state of technology. For example, the normal CDF is one of the many CDFs that cannot be expressed in closed form.

  • The first method of simulating normal variates in the 1950's, when computation beyond basic arithmetic was expensive, was to use $Z = \sum_{i=1}^{12} U_i - 6 \stackrel{aprx}{\sim}\mathsf{Norm}(0,1),$ by the CLT, where $U_i \stackrel{iid}{\sim}\mathsf{Unif}(0,1).$
  • Subsequently, the Box-Muller method was used to obtain two independent standard normal random variables from to independent standard uniform ones.
  • Currently, it is common to use a very accurate rational approximation to the standard normal quantile function to get one standard normal deviate from one standard uniform deviate. I believe the function runif in R uses Michael Wichura's rational approximation, which is accurate up double-precision computer representation.
$\endgroup$
4
  • 3
    $\begingroup$ Your discussion doesn't look relevant to the question, because every distribution you mention has an invertible CDF. I suspect you might have interpreted "not invertible" as meaning "it might be expensive to compute the inverse." $\endgroup$
    – whuber
    Aug 8, 2020 at 16:48
  • $\begingroup$ Well, except for some gammas. But yes, sometimes the methods using tricks are faster than quantile. $\endgroup$
    – BruceET
    Aug 8, 2020 at 16:50
  • 3
    $\begingroup$ I don't follow--all Gamma CDFs are invertible and numerically inverting them is not expensive. $\endgroup$
    – whuber
    Aug 8, 2020 at 16:52
  • $\begingroup$ Sorry, I was thinking about the use of $-\sum_{i=1}^n \ln U_i \sim \mathsf{Gamma}(n, 1),$ with $U_i$ iid standard uniform as among the 'tricks'. Should have given the example instead of using imprecise language. $\endgroup$
    – BruceET
    Aug 8, 2020 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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