Sampling from discrete distribution - running time of accept reject

Question

I have a question about the running time of the accept/reject method to sample from discrete distributions.

I'm sampling from two distributions built like this:

• draw $n_{outcomes}$ number from an exponential distribution (Pareto with $\alpha = 2$ in the second case)

• compute the sum $Z$

• divide each number by $Z$

That is $p_i = \frac{R_i}{\sum_i R_i}$, $R_i$ random number drawn from Exponential or Pareto distribution.

For the exponential case:

def create_discrete_dist_exp(n_items):
    # creates a collection of floats drawn 
    # from an exponential distribution
    # 1-rand() to avoid log(0), since rand() 
    # returns a float in [0, 1[
    collection = -np.log(1 - np.random.rand(n_items))

    # compute the normalization
    Z = sum(collection)

    #return the discrete distribution
    return collection/Z

For the Pareto case:

def create_discrete_dist_sqrt(n_items):
    # create a collection of floats 
    # distributed following a Pareto
    collection = 1/np.sqrt(1-np.random.rand(n_items))

    # compute the normalization
    Z = sum(collection)

    # return the distribution
    return collection/Z

The accept/reject method consists in taking the max of these $p_i$ and build a rectangle of width $n_{outcomes}$, height $\max_i p_i$. Then draw a random point in the rectangle and see if it is under the probability distribution, in that case accept, otherwise reject.

def accept_reject(prob, n):
height = np.amax(prob)
width = len(prob)

samples = []
for i in xrange(0, n):
    while(True):
        x = np.random.rand()*width
        y = np.random.rand()*height

        # cast to int to avoid warnings 
        # because vector index is a float
        ix = int(np.floor(x))

        if prob[ix] > y:
            samples.append(ix)
            break

return samples

The question is how does the running time increase with $n_{outcomes}$?

A clue is that the largest number obtained by sampling N numbers from an exp distribution is $O(\log N)$, so my reasoning is: the largest number is $O(\log N)$, the others will be small, since the exponential distribution drops pretty quickly, therefore I can approximate the difference between the highest rectangle and the others with $O(\log N)$, which is also my guess for the rejection probability, therefore the running time will be $\sim \log N$. But:

• I'm not completely convinced by this reasoning

• I don't know what's the estimate for the largest number sampled in N iterations from a Pareto.

EDIT

running time (RT) and rejection rate (rr) in function of $n_{outcomes}$

fit of $\max p(x)$ with $f(x)=a*\log x + b$

Answer 1

Let's consider two cases, the first where the probabilities are bounded above and the second where they are not. An example of the first case is the Exponential distribution, the second is a Gamma distribution with shape parameter $< 1$.

Let's denote the probabilities by $p$ and, in the first case, their $\sup$ by $p_{sup}$. It should be clear that in the first case the acceptance rate will approach $\text{E} [p/p_{sup}]$. Since this is evidently a constant, the runtime for your algorithm will be $O(1)$ in $N$.

In the second case, $p$ is unbounded above, and, intuitively, we expect the acceptance rate to go to zero / the runtime to go to $\infty$ as $N \to \infty$. It seems to me that how fast it does depends upon the shape of the distribution, but it may be that the constraint imposed on the functional form by the requirement that the distribution integrate to 1 is strong enough to provide an upper bound on the convergence rate.

I created some simulations in the second case, using a Gamma distribution with shape parameter 0.2 as my baseline, to investigate empirically the convergence rate. The R code, written for clarity rather than efficiency, is:

foo <- function(N) {
   sum <- 0
   for (j in 1:1000) {
      x <- rgamma(N, 0.2)
      sum <- sum + mean(x)/max(x)
   }
   sum <- sum / 1000
}
N <- c(10,100,1000,10000,100000)
a <- rep(0,length(N))
for (i in 1:length(N)) a[i] <- foo(N[i])

plot(1/a ~ log(N), lty='b', lwd = 2)

The plot indicates, but of course is not proof, that the increase in runtime was proportional to $\log N$:

Answer 2

Finaly, I figured it out.

Defining the rejection rate (RR) as the ratio between the rejection area and the total area of the rectangle, we have:

$$ RR = \left ( N \cdot p_{max} - \displaystyle \sum_{i=1}^N p_i \right )\frac{1}{N \cdot p_{max}} = 1-\frac{1}{N \cdot p_{max}} $$

since the probability distribution is normalized ($\sum_{i=1}^N p_i = 1$). The acceptance rate (AR) is just

$$ AR = \frac{1}{N \cdot p_{max}} $$

This is the probability of success, so the expected number of iterations (i.e. the running time of the algorithm) is just the inverse of this:

$$ RT = N \cdot p_{max} $$

We just need to find its asymptotic behavior with $N$, and it depends on the asymptotic behavior of $p_{max}$. Let's discuss $p_{max}$, by definition it is:

$$ p_{max} = \frac{\max_i \{ R_i\}}{\displaystyle \sum_{i=1} R_i} $$

where $R_i$ are the random numbers drawn in this case from a uniform distribution $U(0, 1)$. The denominator is extensive, so it's proportional to N, for the numerator we have to use limit theorems for extremes from probability theory.

Limit theorems for extremes

We're considering a sequence of $\{R_i\}$ i.i.d. random variables, we can define a random variable $Z_N = \max \{ R_1, R_2, \ldots, R_N \}$. Our problem is to find two coefficient $a_N$, $b_N$ such that:

$$ Z_N = a_N + b_N \Lambda $$

and

$$ P\{\Lambda < x\} \to H(x) \mbox{ as } N \to \infty $$

and H(x) is non degenerate (it does not take a single value).

The idea is that

$$ P\{Z_N \le x\} = P\{R_i \le x\, \forall i\} = F(x)^N $$

Therefore we need to find

$$ P\{Z_N \le a_N + b_N x\} = F(a_N + b_N x)^N \to H(x) \mbox{ as } N \to infty $$

This means that this condition must realize: $F(a_N + b_N x) \simeq 1$, so we can expand

$$ P\{Z_N \le a_N + b_N x\} = e^{-N\ln F} \simeq e^{-N(1-F)} $$

therefore $a_N$ and $b_N$ must be such that

$$ \lim_{N \to \infty} N(1-F) = \lim_{N \to \infty} N \int_{a_N + b_N x}^{\infty} dy f(y) = c(x) \mbox{ is finite.} $$

In the case of a uniform distribution $U(0,1)$, we have possible sequences:

$$ a_N = 1 $$

$$ b_n\mbox{: } N \int_{1-b_N}^{1} dy = 1 \to b_N = \frac{1}{N} $$

and the limiting distribution for the max is an exponential distribution. Using our results for $a_N$ and $b_N$ we can estimate the max:

$$ Z_N = 1 + \frac{1}{N} \Lambda = O(1) $$

where $\Lambda$ is distribution according to $P(x)$ which converges to an exponential distribution as $N \to \infty$.

Consequently the estimate for $p_{max}$ is $p_{max} = O(1/N)$, and therefore the running time (RT) of the algorithm can be estimated as the rejection rate

$$ RT = N \cdot p_{max} - 1 = O(1) $$

For the exponential distribution the possible sequences are:

$$ a_N \mbox{: } N\int_{a_N}^{\infty} dx e^{-x} = 1 \to a_N = \log N $$

$$ b_N\mbox{: } b_N = \frac{\int_{a_N}^\infty dy \int_{y}^\infty dx e^{-x}}{\int_{a_N}^\infty dx e^{-x}} = \frac{\log N}{N} $$

therefore the estimate for the max is:

$$ Z_N = \log N + \frac{\log N}{N} \Lambda = O(\log N) $$

where $\Lambda$ is a random variable with probability distribution which converges to a Gumbel distribution as $N \to \infty$.

$$p_{max} = O \left (\frac{\log N}{N} \right )$$

and the running time of the accept/reject algorithm is

$$ RT = O(\log N) $$

For a Pareto distribution with $x_m=1$ and $\alpha=2$, the sequences are:

$$ a_N = 0 $$

$$ b_N\mbox{: } N \int_{b_N}^{\infty} dx 2x^{-3} = 1 \to b_N = \sqrt{N} $$

therefore

$$ Z_N = \sqrt{N} \Lambda = O(\sqrt{N}) $$

where $\Lambda$ is a random variable with a distribution which converges to a Fréchet $P(x) = e^{-x^{-2}}$, which has mean $\mathbb{E}[\Lambda] = \sqrt{\pi}$ and infinite variance. From this formula it is also evident that the fluctuations of the max are of the same order of the fluctuations of the sum.

A naive analysis may conclude that

$$ p_{max} = O(\sqrt{N}/N) $$

and finally

$$ RT = O(\sqrt{N}) $$

but being the variance of the Fréchet divergent, there is a finite probability that $Z_N$ will strongly deviate from the mean and the behavior $O(\sqrt{N})$ will be hidden by a large amount of noise.