How to generate sorted uniformly distributed values in an interval efficiently?

Question

Let's say I want to generate a set of random numbers from the interval (a, b). The generated sequence should also have the property that it is sorted. I can think of two ways to achieve this.

Let n be the length of the sequence to be generated.

1st Algorithm:

Let `offset = floor((b - a) / n)`
for i = 1 up to n:
   generate a random number r_i from (a, a+offset)
   a = a + offset
   add r_i to the sequence r

2nd Algorithm:

for i = 1 up to n:
    generate a random number s_i from (a, b)
    add s_i to the sequence s
sort(r)

My question is, does algorithm 1 produce sequences that is as good as the ones generated by algorithm 2?

Answer 1

The first algorithm fails badly for two reasons:

1. Taking the floor of $(a-b)/n$ can reduce it drastically. Indeed, when $b-a \lt n$, it will be zero, giving you a set whose values are all the same!

2. When you don't take the floor, the resulting values are too evenly distributed. For instance, in any simple random sample of $n$ iid uniform variates (say between $a=0$ and $b=1$), there is a $(1-1/n)^n\approx 1/e\approx 37\%$ chance that the largest will not be in the upper interval from $1-1/n$ to $1$. With algorithm 1, there is a $100\%$ chance that the maximum will be in that interval. For some purposes this super-uniformity is good, but in general it is a terrible error because (a) many statistics will be ruined but (b) it can be very difficult to determine why.

3. If you want to avoid sorting, instead generate $n+1$ independent exponentially-distributed variates. Normalize their cumulative sum to the range $(0,1)$ by dividing by the sum. Drop the largest value (which will always be $1$). Rescale to the range $(a,b)$.

Histograms of all three algorithms are shown. (Each depicts the cumulative results of $1000$ independent sets of $n=100$ values each.) The lack of any visible variation in the histogram for Algorithm 1 shows the problem there. The variation in the other two algorithms is exactly what is to be expected--and what you need from a random number generator.

For many more (amusing) ways to simulate independent uniform variates, see Simulating draws from a Uniform Distribution using draws from a Normal Distribution.

Here is the R code that produced the figure.

b <- 1
a <- 0
n <- 100
n.iter <- 1e3

offset <- (b-a)/n
as <- seq(a, by=offset, length.out=n)
sim.1 <- matrix(runif(n.iter*n, as, as+offset), nrow=n)
sim.2 <- apply(matrix(runif(n.iter*n, a, b), nrow=n), 2, sort)
sim.3 <- apply(matrix(rexp(n.iter*(n+1)), nrow=n+1), 2, function(x) {
  a + (b-a) * cumsum(x)[-(n+1)] / sum(x)
})

par(mfrow=c(1,3))
hist(sim.1, main="Algorithm 1")
hist(sim.2, main="Algorithm 2")
hist(sim.3, main="Exponential")

Answer 2

The first algorithm producess too evenly spaced numbers

See also low discrepancy series.

Assuming you want 2 random numbers in $[0;1]$. With real uniform data, the chance is 50:50 they are both larger or smaller than 0.5 at the same time. With your approach, the chance is 0. So your data is not uniform.

(As pointed out, this may be a desired property e.g. for stratification. Low-discrepancy series like Halton and Sobel do have their use cases.)

A proper but expensive approach (for real values)

... is to use beta-distributed random numbers. The rank order statistic of the uniform distribution is beta distributed. You can use this to randomly draw the smallest, then the second smallest, ... repeat.

Assuming the data is to be generated in $[0;1]$. The smallest value is $\text{Beta}[1,n]$ distributed. (For subsequent cases, reduce $n$ and rescale to the remaining interval). To generate a general beta random, we would need to generate two Gamma distributed random values. But $1-X\sim \text{Beta}[n, 1]$. Then $-\ln (1-X)\sim \text{Exponential}[n]$. We can sample random numbers from this distribution as $\frac{-\ln(U[0;1])}{n}$ for this.

\begin{align*} -\ln (1-x) &= \frac{-\ln(1-u)}{n} \\ 1-x &= u^\frac{1}{n} \\ x &= 1 - u^\frac{1}{n} \end{align*}

Which yields the following algorithm:

x = a
for i in range(n, 0, -1):
    x += (b-x) * (1 - pow(rand(), 1. / i))
    result.append(x) 

There may be numerical instabilities involved, and computing pow and a division for every object may turn out to be slower than sorting.

For integer values you may need to use a different distribution.

Sorting is incredibly cheap, so just use it

But don't bother. Sorting is so ridiculously cheap, so just sort. Over the years, we have well understood how to implement sorting algorithms that sorting doubles is not worth avoiding. Theoretically it's $O(n \log n)$ but the constant term is so ridiculously small in a good implementation that this is the perfect example how useless theoretical complexity results can be. Run a benchmark. Generate 1 million randoms with and without sorting. Run it a few times, and I wouldn't be surprised if quite often the sorting beats the non-sorting, because the cost of sorting will still be much less than your measurement error.

Answer 3

It also depends on what you are doing with the random numbers. For numerical integration problems method one(when corrected by removing the floor operator) would produce superior point set. What you are doing is a form of stratified sampling and it has the advantage that it avoids clumping. it's impossible to get all your values in 0-(b-a)/n range for example. That said for other applications this could be very bad, it depends on what you want to do with it.