# How can I sample from a distribution with incomputable CDF?

Semi-computer science simulation related problem here.

I have a distribution where

P(x) = $$\frac{(e^b-1) e^{b (n-x)}}{e^{b n+b}-1}$$

for some constants b and n, and x is an integer such that $$0\leq x \leq n$$.

Now, I need to sample from this distribution. It has an invertible CDF, so it's possible to do this directly in theory. The problem is that the numbers involved are LARGE. So large in fact, that they both overflow conventionally formatted variables, and take at least minutes (at some point I gave up...) to compute using arbitrary precision formats. Basically, the inverse CDF still involves a term of $$e^{b(n+1)}$$, for $$350 < n < 3500$$. Despite this, the output numbers will still be in the range $$0-n$$, so it seems like there should be a way of doing this.

What I'm looking for is a way of approximately sampling from this distribution that is computable. Are there alternative methods of sampling? What are they?

• Have you considered normalizing or scaling your data in some way so to reduce the domain? Commented Apr 9, 2013 at 0:23

The CDF is readily invertible. A formula for the inversion leads to what has to be one of the simplest and most expedient possible solutions.

Begin by observing that the probability of outcome $k$, $0 \le k \le n$, is proportional to $e^{-b k}$. Thus, if we generate a uniform value $q$ between $0$ and $q_{\max}=\sum_{k=0}^{n} e^{-b k}$ = $(1 - e^{-b(n+1)})/(1 - e^{-b})$, we only need find the largest $k$ for which

$$q \ge \sum_{i=0}^{k} e^{-bi} = \frac{1 - e^{-(k+1)b}}{1 -e^{-b}}.$$

Simple algebra gives the solution

$$k = - \text{ceiling}\left(\frac{\log(1 - q (1-e^{-b}))}{b}\right).$$

Here is an R implementation constructed like all the other random-number generators: its first argument specifies how many iid values to generate and the rest of the arguments name the parameters ($b$ as b and $n$ as n.max):

rgeom.truncated <- function(n=1, b, n.max) {
a <- 1 - exp(-b)
q.max <- (1 - exp(-b*(n.max+1))) / a
q <- runif(n, 0, q.max)
return(-ceiling(log(1 - q*a) / b))
}


As an example of its use, let's generate a million random variates according to this distribution:

b <- 0.001
n.max <- 3500
n.sim <- 10^6
set.seed(17)
system.time(sim <- rgeom.truncated(n.sim, b,n.max))


($0.10$ seconds were needed.)

h <- hist(sim+1, probability=TRUE, breaks=50, xlab="Outcome+1")
pmf <- exp(-b * (0: n.max)); pmf <- pmf / sum(pmf)
lines(0:n.max, pmf, col="Red", lwd=2)


($1$ was added to each value in order to create a better histogram: R's hist procedure has an idiosyncrasy (=bug) in which the first bar is too high when the left endpoint is set at zero.) The red curve is the reference distribution that this simulation attempts to reproduce. Let's evaluate the goodness of fit with a chi-square test:

observed <- table(sim)
expected <- n.sim * pmf
chi.square <- (observed-expected)^2 / expected
pchisq(sum(chi.square), n.max, lower.tail=FALSE)


The p-value is $0.84$: a beautiful fit.

• Great solution. I never knew one could sample this way (that is, relying on samples from $Uni(0, k), k>1$ instead of $Uni(0,1)$), but it is obvious in retrospect. Commented Apr 10, 2013 at 16:46

You're dealing with a truncated geometric distribution with $p = 1-e^{-b}$. There are a variety of ways of approaching this.

I'd advise different options in different situations; some options would involve simulating from a geometric and regenerating when its outside the range, taking the integer part of an appropriate truncated exponential (as here), or using any of several fast techniques tailored to discrete distributions over a finite range. Given that $n$ is large, taking the floor of a truncated exponential is likely to be relatively fast, but whether it's the best choice also depends on $b$.

Here's a related question on math.SE

Before I attempt specific suggestions, what's a typical range of values for $b$?

• Thanks for your answer! b ~ ln(1 + epsilon), where epsilon is an additional parameter > 0. Commented Apr 9, 2013 at 0:28
• So you have converted my question about the typical range of b to one about the typical range of ε. Before I attempt specific suggestions, what's a typical range of values for ε? Commented Apr 9, 2013 at 4:03
• The reason I ask is which particular approaches are more efficient depends on the characteristics of the situation. It sounds like you're already happy with the other answer, so maybe it's not worth worrying about additional potential efficiency. Commented Apr 9, 2013 at 4:12
• @JohnDoucette : If b is almost zero then your distribution is almost uniform over $\{0,\ldots,n\]$ hence you can use the uniform as a proposal in an accept reject algorithm since the upper bound should not be terrible. Commented Apr 9, 2013 at 13:32
• @Xi'an You would need $nb$ fairly small rather than $b\approx 0$ before it would be appropriate to use a uniform distribution, because the acceptance rate is $(1- e^{-(n+1)b})/((n+1)(1-e^{-b}))$ $\approx (1-\exp(-n b))/(n b)$, which will be inefficiently low when $nb \gg 1$.
– whuber
Commented Apr 10, 2013 at 17:10

First, note that $P(x)\propto e^{-bx}$ which, if $x$ were continuous, would be related to an exponential distribution. Then, what you can do is to simulate from a truncated exponential distribution and take the floor() (integer part) of the observations.

The cdf of a truncated exponential is

$$F(x;n,b)= \dfrac{1-e^{-bx}}{1-e^{-bn}}.$$

Then, if we make $F(x;n,b)=u$, we obtain that $x=-\dfrac{1}{b}\log[1-u(1-e^{-bn})]$. If $bn$ is large, then $e^{-bn}\approx 0$ which suggest to approximate $x\approx -\dfrac{1}{b}\log[1-u]$.

rweirdp <- function(ns,n,b){
u <- runif(ns)
samp <- - log(1-u*(1-exp(-n*b)))/b
return(floor(samp))
}

rweirdp(1000,10,1)

• I think this is basically what I was looking for. bn is always very large, proportionate sampling makes sense. Wasn't aware of the mapping, though it's clear in retrospect. Thanks! Commented Apr 9, 2013 at 0:45
• I am glad to see that helped. I think I did not explained properly but this approach produces samples from the exact target distribution. Cheers. Commented Apr 9, 2013 at 8:46
• @Xi'an Aren't the weights the same if one uses the value $e^{-bn}$ and take the integer part? Commented Apr 9, 2013 at 12:19
• @Xi'an I think that term appears in the numerator of $P(x)$, up to a factorisation ... Commented Apr 9, 2013 at 12:48
• @Xi'an Actually, this works provided rweirdp is modified to change n to n+1. (As given here, it will never return a value equal to n: that is the effect of the approximation.) A slightly more rigorous analysis is given in my answer. Although I obtain a different-appearing formula, it is equivalent to the (simpler!) one given here, once the n-->n+1 modification is made.
– whuber
Commented Apr 10, 2013 at 17:06

A way to sample from the target distribution $p(k)\propto \exp\{-bk\}$ is to

1. run a Metropolis-Hastings experiment to determine the (interesting) support of the distribution, i.e. in which subset of $\{0,1,\ldots,n\}$ it concentrates;

metro=function(N,b,n){
x=sample(0:n,N,rep=TRUE)
for (t in 2:N){
x[t]=prop=x[t-1]+sample(c(-1,1),1)

if ((prop<0)||(prop>n)||(log(runif(1))>b*(x[t]-prop)))
x[t]=x[t-1]
}
return(x)
}

2. Use the support thus determined, $\{k_0,\ldots,k_1\}$ say, to compute the exact probabilities as $p(k)\propto \exp\{-bk+bk_0\}$ to avoid overflows.

Update: When thinking more about it, since $p(\cdot)$ is decreasing in k, the effective support of the distribution will always start at $k_0=0$. If $b$ is quite large, this support will end very quickly, in which case $n$ does not matter much as large values of $k$ will never be visited. If $b$ is very small, the pdf is almost flat, which means that one can use a uniform distribution on $\{0,1,\ldots,n\}$ as an accept-reject proposal. And use logs in the acceptance step to avoid overflows.