43
$\begingroup$

I'm looking for a way to generate random numbers that appear to be uniform distributed -- and every test will show them to be uniform -- except that they are more evenly distributed than true uniform data.

The problem I have with the "true" uniform randoms is that they will occasionally cluster. This effect is stronger at a low sample size. Roughly said: when I draw two Uniform randoms in U[0;1], chances are around 10% that they are within a range of 0.1, and 1% that they are within 0.01.

So I'm looking for a good way to generate random numbers that are more evenly distributed than uniform randoms.

Use case example: say I'm doing a computer game, and I want to place treasure randomly on a map (not caring about any other thing). I don't want the treasure to be all in one place, it should be all over the map. With uniform randoms, if I place, say, 10 objects, the chances are not that low that there are 5 or so really close to each other. This may give one player an advantage over another. Think of minesweeper, chances (albeit low, if there are enough mines) are that you are really lucky and win with a single click.

A very naive approach for my problem is to divide the data into a grid. As long as the number is large enough (and has factors), one can enforce extra uniformness this way. So instead of drawing 12 random variables from U[0;1], I can draw 6 from U[0;.5] and 6 from U[0.5;1], or 4 from U[0;1/3] + 4 from U[1/3;2/3] + 4 from U[2/3; 1].

Is there any better way to get this extra evenness into the uniform? It probably only works for batch randoms (when drawing a single random, I obviously have to consider the whole range). In particular, I can shuffle the records again afterwards (so it's not the first four from the first third).

How about doing it incrementally? So the first is on U[0;1], then two from each halves, one from each third, one from each fourth? Has this been investigated, and how good is it? I might have to be careful to use different generators for x and y to not get them correlated (the first xy would always be in the bottom half, the second in the left half and bottom third, the third in center third and top third... so at least some random bin permutation is also needed. and in the long run, it will be too even, I guess.

As a side node, is there a well-known test whether some distribution is too evenly distributed to be truly uniform? So testing "true uniform" vs. "someone messed with the data and distributed the items more evenly". If I recall correctly, Hopkins Statistic can measure this, but can it be used for testing, too? Also somewhat an inverse K-S-Test: if the largest deviation is below a certain expected threshold, the data is too evenly distributed?

$\endgroup$
  • 7
    $\begingroup$ Have you heard of Halton sequences? For "too evenly," people (beginning with Fisher's investigation of Mendel's pea experiment results) have referred to the (usual) chi-squared statistic to the lower tail of a chi-squared distribution. $\endgroup$ – whuber Oct 14 '12 at 16:17
  • $\begingroup$ One way of formalizing this would be to want a distribution $g(x_1, ..., x_n)$ such that (1) $g(\cdot)$ marginalizes to $1$ over $x_1, ..., x_{n - 1}$, (2) $g$ is symmetric, i.e. $X_1, ..., X_n$ are exchangeable, and (3) $g(x_1, ..., x_n)$ is large when $x_1, ..., x_n$ are dispersed. I think there is a real problem with (2) and (3) since infinite exchangeable sequences in $\mathbb R$ can't be negatively correlated, so the larger $n$ we want to use the less repulsion we can enforce; on the other hand, for large $n$, we should have good spread anyways. $\endgroup$ – guy Oct 14 '12 at 17:38
  • $\begingroup$ Halton sequences is pretty close to the approach that I was thinking of. Including skipping the first few entries to reduce the risk of correlation. I was also thinking of using a random permuation for each level. Thank you for this pointer, as this gives me a good point to search for related methods! $\endgroup$ – Anony-Mousse Oct 14 '12 at 19:38
  • $\begingroup$ wrt. Halton sequences again. I need to have them non-deterministic, at least except for an initial seed. I see two ways here. I can do a cyclic shift by a random offset + a random start offset + step size. The problem is that of course the "treasure" to remain to the game example should also not be in the same positions relative to each other each time. Or I could use this uniform-from-subinterval approach that I had in my question to add some amount of "random twist". So to say: Halton seems again too predictable and regular for my use. $\endgroup$ – Anony-Mousse Oct 14 '12 at 20:07
  • 3
    $\begingroup$ en.wikipedia.org/wiki/Low-discrepancy_sequence or mathworld.wolfram.com/QuasirandomSequence.html . Several of the common tests of uniform RNGs (such as those in the Diehard/Dieharder batteries of tests) are sensitive to such things; for example, there are too few 'small distances' between points. $\endgroup$ – Glen_b Oct 14 '12 at 22:44
59
$\begingroup$

Yes, there are many ways to produce a sequence of numbers that are more evenly distributed than random uniforms. In fact, there is a whole field dedicated to this question; it is the backbone of quasi-Monte Carlo (QMC). Below is a brief tour of the absolute basics.

Measuring uniformity

There are many ways to do this, but the most common way has a strong, intuitive, geometric flavor. Suppose we are concerned with generating $n$ points $x_1,x_2,\ldots,x_n$ in $[0,1]^d$ for some positive integer $d$. Define $$\newcommand{\I}{\mathbf 1} D_n := \sup_{R \in \mathcal R}\,\left|\frac{1}{n}\sum_{i=1}^n \I_{(x_i \in R)} - \mathrm{vol}(R)\right| \>, $$ where $R$ is a rectangle $[a_1, b_1] \times \cdots \times [a_d, b_d]$ in $[0,1]^d$ such that $0 \leq a_i \leq b_i \leq 1$ and $\mathcal R$ is the set of all such rectangles. The first term inside the modulus is the "observed" proportion of points inside $R$ and the second term is the volume of $R$, $\mathrm{vol}(R) = \prod_i (b_i - a_i)$.

The quantity $D_n$ is often called the discrepancy or extreme discrepancy of the set of points $(x_i)$. Intuitively, we find the "worst" rectangle $R$ where the proportion of points deviates the most from what we would expect under perfect uniformity.

This is unwieldy in practice and difficult to compute. For the most part, people prefer to work with the star discrepancy, $$ D_n^\star = \sup_{R \in \mathcal A} \,\left|\frac{1}{n}\sum_{i=1}^n \I_{(x_i \in R)} - \mathrm{vol}(R)\right| \>. $$ The only difference is the set $\mathcal A$ over which the supremum is taken. It is the set of anchored rectangles (at the origin), i.e., where $a_1 = a_2 = \cdots = a_d = 0$.

Lemma: $D_n^\star \leq D_n \leq 2^d D_n^\star$ for all $n$, $d$.
Proof. The left hand bound is obvious since $\mathcal A \subset \mathcal R$. The right-hand bound follows because every $R \in \mathcal R$ can be composed via unions, intersections and complements of no more than $2^d$ anchored rectangles (i.e., in $\mathcal A$).

Thus, we see that $D_n$ and $D_n^\star$ are equivalent in the sense that if one is small as $n$ grows, the other will be too. Here is a (cartoon) picture showing candidate rectangles for each discrepancy.

extremal and star discrepancy

Examples of "good" sequences

Sequences with verifiably low star discrepancy $D_n^\star$ are often called, unsurprisingly, low discrepancy sequences.

van der Corput. This is perhaps the simplest example. For $d=1$, the van der Corput sequences are formed by expanding the integer $i$ in binary and then "reflecting the digits" around the decimal point. More formally, this is done with the radical inverse function in base $b$, $$\newcommand{\rinv}{\phi} \rinv_b(i) = \sum_{k=0}^\infty a_k b^{-k-1} \>, $$ where $i = \sum_{k=0}^\infty a_k b^k$ and $a_k$ are the digits in the base $b$ expansion of $i$. This function forms the basis for many other sequences as well. For example, $41$ in binary is $101001$ and so $a_0 = 1$, $a_1 = 0$, $a_2 = 0$, $a_3 = 1$, $a_4 = 0$ and $a_5 = 1$. Hence, the 41st point in the van der Corput sequence is $x_{41} = \rinv_2(41) = 0.100101\,\text{(base 2)} = 37/64$.

Note that because the least significant bit of $i$ oscillates between $0$ and $1$, the points $x_i$ for odd $i$ are in $[1/2,1)$, whereas the points $x_i$ for even $i$ are in $(0,1/2)$.

Halton sequences. Among the most popular of classical low-discrepancy sequences, these are extensions of the van der Corput sequence to multiple dimensions. Let $p_j$ be the $j$th smallest prime. Then, the $i$th point $x_i$ of the $d$-dimensional Halton sequence is $$ x_i = (\rinv_{p_1}(i), \rinv_{p_2}(i),\ldots,\rinv_{p_d}(i)) \>. $$ For low $d$ these work quite well, but have problems in higher dimensions.

Halton sequences satisfy $D_n^\star = O(n^{-1} (\log n)^d)$. They are also nice because they are extensible in that the construction of the points does not depend on an a priori choice of the length of the sequence $n$.

Hammersley sequences. This is a very simple modification of the Halton sequence. We instead use $$ x_i = (i/n, \rinv_{p_1}(i), \rinv_{p_2}(i),\ldots,\rinv_{p_{d-1}}(i)) \>. $$ Perhaps surprisingly, the advantage is that they have better star discrepancy $D_n^\star = O(n^{-1}(\log n)^{d-1})$.

Here is an example of the Halton and Hammersley sequences in two dimensions.

Halton and Hammersley

Faure-permuted Halton sequences. A special set of permutations (fixed as a function of $i$) can be applied to the digit expansion $a_k$ for each $i$ when producing the Halton sequence. This helps remedy (to some degree) the problems alluded to in higher dimensions. Each of the permutations has the interesting property of keeping $0$ and $b-1$ as fixed points.

Lattice rules. Let $\beta_1, \ldots, \beta_{d-1}$ be integers. Take $$ x_i = (i/n, \{i \beta_1 / n\}, \ldots, \{i \beta_{d-1}/n\}) \>, $$ where $\{y\}$ denotes the fractional part of $y$. Judicious choice of the $\beta$ values yields good uniformity properties. Poor choices can lead to bad sequences. They are also not extensible. Here are two examples.

Good and bad lattices

$(t,m,s)$ nets. $(t,m,s)$ nets in base $b$ are sets of points such that every rectangle of volume $b^{t-m}$ in $[0,1]^s$ contains $b^t$ points. This is a strong form of uniformity. Small $t$ is your friend, in this case. Halton, Sobol' and Faure sequences are examples of $(t,m,s)$ nets. These lend themselves nicely to randomization via scrambling. Random scrambling (done right) of a $(t,m,s)$ net yields another $(t,m,s)$ net. The MinT project keeps a collection of such sequences.

Simple randomization: Cranley-Patterson rotations. Let $x_i \in [0,1]^d$ be a sequence of points. Let $U \sim \mathcal U(0,1)$. Then the points $\hat x_i = \{x_i + U\}$ are uniformly distributed in $[0,1]^d$.

Here is an example with the blue dots being the original points and the red dots being the rotated ones with lines connecting them (and shown wrapped around, where appropriate).

Cranley Patterson

Completely uniformly distributed sequences. This is an even stronger notion of uniformity that sometimes comes into play. Let $(u_i)$ be the sequence of points in $[0,1]$ and now form overlapping blocks of size $d$ to get the sequence $(x_i)$. So, if $s = 3$, we take $x_1 = (u_1,u_2,u_3)$ then $x_2 = (u_2,u_3,u_4)$, etc. If, for every $s \geq 1$, $D_n^\star(x_1,\ldots,x_n) \to 0$, then $(u_i)$ is said to be completely uniformly distributed. In other words, the sequence yields a set of points of any dimension that have desirable $D_n^\star$ properties.

As an example, the van der Corput sequence is not completely uniformly distributed since for $s = 2$, the points $x_{2i}$ are in the square $(0,1/2) \times [1/2,1)$ and the points $x_{2i-1}$ are in $[1/2,1) \times (0,1/2)$. Hence there are no points in the square $(0,1/2) \times (0,1/2)$ which implies that for $s=2$, $D_n^\star \geq 1/4$ for all $n$.

Standard references

The Niederreiter (1992) monograph and the Fang and Wang (1994) text are places to go for further exploration.

$\endgroup$
  • 4
    $\begingroup$ This answer is excellent, and I just wanted to appreciate the effort you put into it. Thank you! $\endgroup$ – Anony-Mousse Oct 16 '12 at 23:36
  • 1
    $\begingroup$ One small followup question. Halton sequences look good, because they also appear to be not too regular. The lattice stuff is much to regular for me, and also the Hammersley sequence seems to have a lot of objects on lines through the origin. What is a good way to control a balance between true uniform and fake uniform? Just take 80% contribution from Halton + 20% uniform random? $\endgroup$ – Anony-Mousse Oct 16 '12 at 23:48
  • 1
    $\begingroup$ +10k and definitely with a record low (87!!!!) answers! Oh, and I like this post very much. I bookmarked the question because of it, actually. Well done, @cardinal. $\endgroup$ – Macro Jan 10 '13 at 2:54
  • $\begingroup$ @Macro: Thank you for such a nice comment! You are very kind. I think this 10K thing may be temporary for me. I suspect I may fall well below 10K as soon as Procrastinator's votes are reverted. I'm surprised this hasn't happened, yet, actually. I believe they cast almost 3000 votes on this site. Thanks also for posting here; somehow I never saw Anony-Mousse's follow-up questions! $\endgroup$ – cardinal Jan 10 '13 at 3:06
  • $\begingroup$ @Anony-Mousse: Apologies for the terrible delay in responding. I must have overlooked these comments. I think creating a balance would depend on your goals. Theoretically speaking, introducing any random uniform points are bound to destroy the optimal properties of $D^\star$, for example. As a practical matter, it may be better to use a very small jitter of the QMC points where the jitter is chosen based on the $D^\star$ properties of the sequence. You could also introduce random rigid-body transformations on all the points, e.g., shifts and coordinate rotations. $\endgroup$ – cardinal Jan 18 '13 at 20:49
3
$\begingroup$

One way to do this would be to generate uniform random numbers, then test for "closeness" using any method you like and then delete random items that are too close to others and choose another set of random uniforms to make up for them.

Would such a distribution pass every test of uniformity? I sure hope not! It's no longer uniformly distributed, it is now some other distribution.

One uninuitive aspect of probability is that chance is clumpy. There are more runs in random data than people think there will be. I think Tversky did some research on this (he researched so much, though, that it's hard to remember).

$\endgroup$
  • 2
    $\begingroup$ One of the (many) problems with this approach is it's very hard to characterize the resulting distribution. $\endgroup$ – whuber Oct 14 '12 at 16:18
  • $\begingroup$ The OP seems most concerned with small sample sizes. This would suggest he doesn't need to care about the whole distribution. Suppose you have a set of coordinates, you generate another and then calculate the euclidean distance with respect to all the others. If the smallest distance is below some threshold, throw the number out and generate a new one. I think Peter's solution works fine. $\endgroup$ – John Oct 14 '12 at 17:21
  • $\begingroup$ @whuber He doesn't seem to be interested in that, although I might be wrong. $\endgroup$ – Peter Flom Oct 14 '12 at 18:18
  • 2
    $\begingroup$ Let me state my objection a little more clearly, Peter: when you remove and/or adjust pseudorandom values in an ad hoc way in order to approximate some desired property, such as lack of clustering, it is difficult to assure that the resulting sequences have any desirable properties. With your method, for instance, could you even tell us what the first moment of the resulting process would be? (That is, can you even assure us that the intensity is uniform?) What about the second moment? Usually these constitute the minimum information needed to use the sequences effectively for inference. $\endgroup$ – whuber Oct 15 '12 at 13:21
  • 2
    $\begingroup$ OK, but, in the example in the question, he wants to place treasure on a map in a game. That won't involve inference or moments or anything of the sort. I admit my method wouldn't be good for a lot of purposes, but I think it matches up with the example. Of course, maybe the example isn't really what he wants.... Maybe he wants something more formal, in which case all the other answers should be looked at. $\endgroup$ – Peter Flom Oct 15 '12 at 13:48
3
$\begingroup$

This is known as a "hard-core" poisson point process - so named by Brian Ripley in the 1970s; i.e. you want it to be random, but you don't want any points to be too close together. The "hard-core" can be imagined as a buffer zone around which other points cannot intrude.

Imagine you're recording the position of some cars in a city - but you're only recording the point at the nominal centre of the car. While they're on the streets no two point pairs can come close together because the points are protect by the "hard-core" of the bodywork - we'll ignore the potential super-position in multi-storey car parks :-)

There are procedures for generating such point processes - one way is just to generate points uniformly and then remove any that are too close together!

For some detail on such processes, refer for example to this

$\endgroup$
2
$\begingroup$

With respect to batch generation in advance, I would generate a large number of sets of pseudorandom variates, and then test them with a test such as the Kolmogorov-Smirnov test. You will want to select the set that has the highest p-value (i.e., $p \approx 1$ is ideal). Note that this will be slow, but as $N$ gets larger it probably becomes less necessary.

With respect to incremental generation, you essentially are looking for a series with a moderately negative autocorrelation. I'm not sure what the best way to do that would be, since I have very limited experience with time-series, but I suspect there are existing algorithms for this.

With respect to a test for "too even", any test for whether a sample follows a specific distribution (such as the KS noted above) will do, you just want to check if $p > (1-\alpha)$, rather than the standard approach. I wrote about an example of this alternative approach here: chi-squared always a one-sided test.

$\endgroup$
1
$\begingroup$

I would formalize your problem this way: You want a distribution over $[0,1]^n$ such that the density is $f(x) \propto e^{\left(\frac1k\sum_{ij}\lvert x_i-x_j \rvert^{k}\right)^{\frac1k}}$ for some $k<0$ quantifying the repulsion of points.

One easy way to generate such vectors is to do Gibbs sampling.

$\endgroup$
  • $\begingroup$ Can you elaborate on this? Gibbs sampling does not seem to help here, as conditional distribution = marginal distribuion = uniform? Or is your suggestion to use the previous samples to produce "holes" in the distribution to sample from? $\endgroup$ – Anony-Mousse Oct 14 '12 at 19:52
  • $\begingroup$ Pick a uniform random vector, and then repeatedly uniformly choose an index $i$ and resample $x_i$. Calculate the ratio $r$ of $f(x)$ before and after the resampling and reject your resampling with odds $r$. This is much faster than the other answers you've gotten when you have a very long vector because you are performing local rather than global rejections. $\endgroup$ – Neil G Oct 14 '12 at 21:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.