# Sampling from a distribution with given probabilities

I have found myself needing to sample numbers from a given distribution and was wondering what the most efficient method of doing so is.

I have a set of N elements 0 to N-1 and an array mapping those elements to the probability of their selection, pp.

My approach is to create a new array, cdf, and accumulate the entries of pp so that the first element of cdf is pp[0] and the last element is the sum of all the elements in pp, namely 1. Then, since this array is necessarily sorted, I choose a random number between 0 and 1 and perform a binary search in the array.

This is my code so far (in Java):

public class Sample {
/**
* Selects q elements chosen at random according to the probabilities
* given in pp.
*
* @param pp  The probabilites corresponding to each index
*            (item 0 has probability pp[0] and so forth)
* @param q   The number of elements to sample
* @return    An array containing the elements sampled
*/
public static int [] sample(double [] pp, int q) {
double [] cdf = pp.clone();
for (int i = 1; i < cdf.length; i++)
cdf[i] += cdf[i - 1];

int [] values = new int[q];
for (int i = 0; i < q; i++) {
// binarySearch returns a negative number at the insertion point
// if the exact number isn't found. We don't expect to hit an
// exact match with Math.random(), so the Math.abs() call
// can't be too expensive.
values[i] = Math.abs(Arrays.binarySearch(cdf, Math.random()));
}

return values;
}
}


I'm mostly unsure about the behavior of the Arrays.binarySearch function. Am I at risk of missing the first/last element, or going out of bounds? The running time of this algorithm is $O(N)$ for the cdf construction, plus $O(q \log N)$ for the actual sampling. The second operation dominates ($q$ could be greater than $N$), so this is an $O(q \log N)$ algorithm, which is slower than sorting! Is there any way to do this more quickly?

Boundary issues.

The use of double brings some precision issues. For instance, if you have seven elements, each with a probability of 1/7, then your cdf array will contain :

0 : 0.14285714285714285d
1 : 0.2857142857142857d
2 : 0.42857142857142855d
3 : 0.5714285714285714d
4 : 0.7142857142857142d
5 : 0.857142857142857d
6 : 0.9999999999999998d


If we feed this to your code, and manipulate the random value, we can test what happens at border cases.

• If random returns 0d, your method returns 1.
• If random returns 0.15d, your method returns 2.
• If random returns 0.9999999999999997d your method returns 7.
• If random returns 0.9999999999999998d your method returns 6.
• If random returns 0.9999999999999999d your method returns 8.

So yes, there are boundary issues.

1. You treat the output of binarySearch the same when the value is found as when it is not found. This explains the 6 for 0.9999999999999998d. Your documentation says you're gambling on this not occurring, but you're motivation to do so is not good. You're simply guessing it will aversely affect performance, to check the output. But very likely this is not going to be a performance bottleneck at all.
2. That 6 is actually correct, since binarySearch, when not finding the value returns -i-1 where i is the insertion index. By interpreting this value as you do, you end up with 1-based indexes for most runs.
3. Due to imprecision, you can get return values that are out of bounds (the 8). While this is the same sort of gamble you could make as in point 1, in this instance you are not just sacrificing some precision, but algorithm correctness. In the field this may lead to a bug that will be very hard, if not impossible to find.

Solution to all these is correctly handling the result of binarySearch(). If it is negative convert it to the insertion index. If it is positive, use as is. If the insertion index is at cdf.length, ignore this and do a reroll.

Other issues.

• As others have pointed out : naming can be improved. use names that explain what you are doing.
• cloning seems unnecessary, you end up still replacing all but one value in the cloned array. Simply create a new array of the appropriate size and populate it.

General performance.

I don't think there is a faster way to do this (in terms of O() ). If q == 1, you could just calculate the cumulative probabilities until you go over the random() value. But for q > 1, it simply makes sense to reuse the already calculated cumulative probabilities and apply binarySearch.

• Thanks for your clarification of the binarySearch function. That +1 was a bit confusing. It turns out that there actually is a faster way to do this, based on the Alias Method -- See here: jstatsoft.org/v11/i03/paper – Alex Reinking Jun 30 '14 at 18:50

Just a suggestion:

I think you can use better variable names than pp and q.

cdf can be clonedArray, and the same for other variable names.

This doesn't matter so much but increases the readability of the program.

• You're probably right about that. I was sticking a bit too closely to the mathematical notation... – Alex Reinking Jun 27 '14 at 5:12

You return indices 1 to N, but I think it's better from 0 to (N - 1). Unless you generate billions of values, I would stick with collections instead of arrays (there is a boxing/unboxing cost when using collections).

I made the code OO. And split the methods such that each method only does one think (separation of concerns).

public class ProbabilitySampler {
private Random random = new Random();
private List<Double> probabilities;
private double[] cdf;

public ProbabilitySampler(List<Double> probabilities) {
super();
// TODO check sum of probabilities is very close to 1.0
this.probabilities = Collections.unmodifiableList(probabilities);
this.cdf = buildCDF(this.probabilities);
}

public static double[] buildCDF(List<Double> probabilities) {
double[] cdf = new double[probabilities.size()];
cdf[0] = probabilities.get(0);
for (int i = 1; i < probabilities.size(); i++)
cdf[i] = cdf[i - 1] + probabilities.get(i);
return cdf;
}

public Integer sample() {
int index = Arrays.binarySearch(cdf, random.nextDouble());
return (index >= 0) ? index : (-index - 1);
}

public static void main(String[] args) {
List<Double> probabilities = Arrays.asList(0.32, 0.68);
ProbabilitySampler probabilitySampler = new ProbabilitySampler(probabilities);

int nSamples = 100000;
final List<Integer> distribution = new ArrayList<>(Collections.nCopies(probabilities.size(), 0));
IntStream
.range(0, nSamples)
.map(i -> probabilitySampler.sample())
.forEach(randomItem -> distribution.set(randomItem, distribution.get(randomItem) + 1));
System.out.println(distribution);
}
}


[EDIT: I had originally defined cdf as a List, but I changed it to an array after noticing that Arrays.binarySearch is significantly faster than Collections.binarySearch.]

There is some trick you can use. It's a bit far fetched however, and I don't recommend using it. Also, the complexity remains $log(N)$, but it could provide a big speed up, if your array of probabilities is really long, which I doubt is the case.

You basically restrict the initial binary search range by using a "discretized" CDF array which contains indices to the real CDF array. For example, lets say you have only 2 values: the first one with probability 0.32 and the second one with probability (1 - 0.32).

int[] discretizedCDF = new int[]{0, 0, 0, 1, 1, 1, 1, 1, 1, 1};
double rand = Math.random();
double discretizedRand = rand * discretizedCDF.length;


And then from this discretizedRand and discretizedCDF you can get lower and upper indices for the binary search in the real CDF array. There is a version of Arrays.binarySearch that takes a starting range. So you still have to carry the binary search, but the search span can be greatly reduced. The longer discretizedCDF is, the shorter the search span.

I was lazy and did not write the full algorithm, but I hope you understand what I meant. Also, you would write a method to generate discretizedCDF from the true CDF instead of writing by hand as I did in the example.

Again, if you have few elements in your probabilities array, it is not worth it. And this trick does not provide any speed up if all the elements are clustered in a small range in the CDF.