# Sampling from discrete probability distribution from first principles

I have a set $S=\{a_1,a_2,\dots,a_n\}$. The probability with which each of the element is selected is $\{p_1,p_2,\dots,p_n\}$ respectively (where of course $p_1+p_2+\cdots+p_n=1$).

I want to simulate an experiment which does that. However I wish to do that without any libraries (i.e. from first principles).

I'm using the following method:

1) I map the elements on the real number line as follows $X(a_1)=1$; $X(a_2)=2$;$~\dots$;$X(a_n)=n$.

2) Then I calculate the cumulative probability distribution function for each coordinate (i.e $P(x < X)$) as follows: $\mathrm{cdf}(x)= P(a_1) + P(a_2) + \cdots + P(a_i)$ such that $X(a_i) \le x < X(a_{i+1})$ (thus the cdf is a step function).

3) I randomly select a real number $q \in (0,1)$ and calculate the $x$-coordinate where the line $y = q$ intersects the cdf. Since the cdf is a step function with jumps at $1,2,\dots,n$ the point would have an integer $x$-coordinate between $1$ and $n$. Let the $x$-coordinate be $m$.

4) I select that $a_i$ such that $X(a_i) = m$.

My question is does this method simulate the experiment without any bias?

I'm not getting the required results, which is why I'm a bit skeptical.

Any help will be greatly appreciated! Thanks!

-
A simple way of getting the result you want is to choose $a_i$ if $q$ is in the $i$-th interval in the following list of $n$ intervals: $(0, p_1], (p_1, p_1+p_2], (p_1+p_2, p_1+p_2+p_3], \ldots, (1-p_n, 1]$ where that last $1-p_n$ is shorthand for $p_1+p_2+\cdots+p_{n-1}$ – Dilip Sarwate May 27 '12 at 19:43

I haven't completely looked at the details of what you are doing, but it looks like you are in the neighborhood of what is typically done. Step 1 is to create a table with $T_i = P(X \le a_i)$ for each $a_i$, which it looks like you have done. Next, generate a uniform random variable, $U \in (0, 1)$. Then, iterate through the table and take $X = a_i$ for the smallest value of $i$ such that $U < T_i$. So, for example, imagine I have uniform mass on $\{1, 2, 3\}$. I draw $U = .4$. $P(X \le 1) = \frac 1 3$, so I don't take $X = 1$. But $P(X \le 2) = \frac 2 3$ and $U < \frac 2 3$ so I take $X = 2$.

There are various ways to optimize this process if speed is an issue. In particular, we aren't using the ordering of the $a_i$ in any essential way, so we might remap the $a_i$ to the integers so that the most likely values are the smaller ones before creating a new table $T_i = P(X_i \le i)$. This way the search terminates faster since we start by checking values which have more mass instead of starting at the lower tail of the distribution.

For countable distributions we can do something similar to this; the trick here is to not have our table be exhaustive and add to it as needed when we generate a particularly large value of $U$.

EDIT: Made the inequalities not strict.

-