# Sampling from discrete distribution - running time of accept reject

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$ • Petersen JL. Estimating the parameters of a Pareto distribution. Univ Montana, Missoula, MT, USA Tech Rep. 2000 May give you some ideas on how to regress Pareto in general.
– Carl
Oct 1, 2016 at 18:55
• yeah thanks, but my question is on the link between running time and type of distribution actually Oct 1, 2016 at 19:31
– Carl
Oct 1, 2016 at 20:12
• I'm not sampling an exponential distribution, I'm sampling a discrete distribution built by drawing random numbers from an exponential and then normalizing (same with Pareto). Btw edited the question, added also the procedure, hope that now it is more clear :) Oct 1, 2016 at 20:18
• Yes, clearer. However, not my area.
– Carl
Oct 1, 2016 at 20:31

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$: • I'm not convinced at 100%. $p_{sup}$ depends on N, at least in the case of exp distribution, since "the largest number obtained by sampling N numbers from an exp distribution is $O(\log N)$", therefore I should expect a decrease in the acceptance rate with $\log N$ and consequently a increase of running time with $\log N$? Do you think it is correct? And in the case of Pareto? I don't have an estimate for that. Oct 1, 2016 at 21:53
• You are right, but you're thinking about it the wrong way (or I am!) The largest number you care about is from $p(x)$, not from $x$ itself. Although you're right about $\max x \to \infty$ being $O(\log N)$ of course, $\max p(x)$ from a standard exponential goes to 1, and rather quickly at that - 10 draws will give you an average $\max p(x)$ of about 0.91, Oct 1, 2016 at 22:13
• I'm quite convinced exponential is $O(1)$ after running some simulations, even though I don't get mathematically how $\max p(x) \to 1$ Oct 1, 2016 at 22:35
• The maximum probability from an exponential distribution is at $x=0$, where $p(x) = 1$ $(\exp(-0))$. As you sample more and more, you're going to get draws closer and closer to 0, which means the maximum $p(x_{sample})$ will approach 1. Oct 1, 2016 at 23:00
• To clarify a little, your acceptance probability for a particular proposal $y$ drawn from your initial sample $x$ will be $p(y) / \max(p(x))$, So what you care about, from the perspective of acceptance probabilities, is $p(x)$, not $x$ itself. Oct 1, 2016 at 23:28

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\}}\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.