Comment générer efficacement des valeurs triées uniformément réparties dans un intervalle?

Disons que je veux générer un ensemble de nombres aléatoires à partir de l'intervalle (a, b). La séquence générée doit également avoir la propriété d'être triée. Je peux penser à deux façons d'y parvenir.

Soit nla longueur de la séquence à générer.

1er algorithme:

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

2ème algorithme:

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

Ma question est la suivante: l'algorithme 1 produit-il des séquences aussi bonnes que celles générées par l'algorithme 2?

Le premier algorithme échoue gravement pour deux raisons:

1. Prendre la parole de peut le réduire considérablement. En effet, quand , il sera nul, vous donnant un ensemble dont les valeurs sont toutes les mêmes!b - a < $(a-b)/n$$b-a \lt n$

2. Lorsque vous ne prenez pas la parole, les valeurs résultantes sont trop uniformément réparties. Par exemple, dans tout échantillon aléatoire simple de iid variables uniformes (disons entre et ), il y a chance que le le plus grand ne sera pas dans l'intervalle supérieur de à . Avec l'algorithme 1, il y a chances que le maximum soit dans cet intervalle. À certaines fins, cette super-uniformité est bonne, mais en général, c'est une erreur terrible car (a) de nombreuses statistiques seront ruinées mais (b) il peut être très difficile de déterminer pourquoi.a = 0 b = 1 ( 1 - 1 / n ) n ≈ 1 / e ≈ 37 % 1 - 1 / n 1 100 $n$$a=0$$b=1$$(1-1/n)^n\approx 1/e\approx 37\%$$1-1/n$$1$$100\%$

3. Si vous voulez éviter le tri, générez plutôt variables indépendantes exponentiellement distribuées. Normaliser leur somme cumulée à la plage en divisant par la somme. Supprimez la plus grande valeur (qui sera toujours ). Redimensionner à la plage .( 0 , 1 ) 1 ( a , b $n+1$$(0,1)$$1$$(a,b)$

Les histogrammes des trois algorithmes sont affichés. (Chacun représente les résultats cumulatifs de ensembles indépendants de valeurs chacun.) L'absence de toute variation visible dans l'histogramme de l'algorithme 1 montre le problème là-bas. La variation dans les deux autres algorithmes est exactement ce à quoi s'attendre - et ce dont vous avez besoin d'un générateur de nombres aléatoires.n = $1000$$n=100$

Pour de nombreuses autres façons (amusantes) de simuler des variations uniformes indépendantes, voir Simulation de tirages à partir d'une distribution uniforme à l'aide de tirages à partir d'une distribution normale .

Voici le Rcode qui a produit la 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")

Le premier algorithme produit des nombres trop espacés

Voir également les séries à faible écart .

En supposant que vous vouliez 2 nombres aléatoires dans . Avec de vraies données uniformes, la probabilité est de 50:50, elles sont à la fois supérieures ou inférieures à 0,5 en même temps. Avec votre approche, la chance est 0. Vos données ne sont donc pas uniformes.$[0;1]$

(Comme l' a souligné, cela peut être une propriété souhaitée , par exemple pour la stratification. Série à faible discrépance comme Halton et Sobel n'ont leurs cas d'utilisation.)

Une approche appropriée mais coûteuse (pour les valeurs réelles)

... consiste à utiliser des nombres aléatoires distribués bêta. La statistique d'ordre de rang de la distribution uniforme est distribuée bêta. Vous pouvez l'utiliser pour dessiner au hasard le plus petit , puis le deuxième plus petit, ... répétez.

$[0;1]$$\text{Beta}[1,n]$$n$$1-X\sim \text{Beta}[n, 1]$$-\ln (1-X)\sim \text{Exponential}[n]$$\frac{-\ln(U[0;1])}{n}$

Ce qui donne l'algorithme suivant:

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

Il peut y avoir des instabilités numériques impliquées, et le calcul powet une division pour chaque objet peuvent s'avérer plus lents que le tri.

Pour les valeurs entières, vous devrez peut-être utiliser une distribution différente.

Le tri est incroyablement bon marché, alors utilisez-le

$O(n \log n)$

Cela dépend également de ce que vous faites avec les nombres aléatoires. Pour les problèmes d'intégration numérique, la première méthode (lorsqu'elle est corrigée en supprimant l'opérateur de plancher) produira un ensemble de points supérieur. Ce que vous faites est une forme d’échantillonnage stratifié et il a l’avantage d’éviter l’agglutination. il est impossible d'obtenir toutes vos valeurs dans la plage 0- (ba) / n par exemple. Cela dit, pour d'autres applications, cela pourrait être très mauvais, cela dépend de ce que vous voulez en faire.