How can I sample from a distribution with incomputable CDF?

Question

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?

Answer 1

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.

Answer 2

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

Answer 3

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)

Answer 4

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.