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

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?

• BTW it's remarkably easy to generate a list of sorted random numbers in R. In order to generate an array of $k$ sets of $n$ random numbers over an uniform interval $[a, b]$, the following code works: rand_array <- replicate(k, sort(runif(n, a, b)). – RobertF Jul 8 '16 at 17:57

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")
• What do you think of the algorithm (based on rank order statistics) in my answer? ;-) – Anony-Mousse Jan 21 '15 at 12:03
• @Anony It is a less efficient version of my algorithm 3. (Yours seems to involve a lot of unnecessary rescaling.) You generate the exponential variates by taking logs of uniforms, which is standard. – whuber Jan 21 '15 at 15:34

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.

• +1 I think this is a useful contribution to the question, especially by characterizing Algorithm 1 in terms of stratification. – whuber Jan 20 '15 at 22:35

The first algorithm producess too evenly spaced numbers

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.

• There can be reasons to avoid sorting. One is when you want to generate a huge number of random variates, so many that a standard sort routine cannot handle them. – whuber Jan 21 '15 at 15:37
• I think the numerical issues with sums using floating point math become a problem much earlier. (And the problems with cyclic patterns in pseudo random numbers!) It's fairly easy to scale the sorting approach to terabytes, and to exabytes on distributed systems. – Anony-Mousse Jan 22 '15 at 7:09
• With scaling that large, the log term starts becoming more ... interesting. Although it's good to be concerned about floating point errors, they are not going to be of any consequence until you are summing more than about $10^{12}$ values and the problem is easily solved (albeit by more programming, I admit) by breaking the sums into subgroups. My point is that when you are performing a calculation that needs to step in sequence through a set of uniform variates, the non-sorting methods completely avoid having to generate, store, and sort all of them initially. – whuber Jan 22 '15 at 16:23
• Ok, not having to store them is an argument. But then you'll need my approach, your variant 3 using the cumulative sum won't work. – Anony-Mousse Jan 22 '15 at 17:56
• That is an excellent point. Now I see the virtue of the extra calculations! (+1) – whuber Jan 22 '15 at 17:59