¿Cómo generar eficientemente valores ordenados distribuidos uniformemente en un intervalo?

12

Digamos que quiero generar un conjunto de números aleatorios a partir del intervalo (a, b). La secuencia generada también debe tener la propiedad de que está ordenada. Puedo pensar en dos formas de lograr esto.

Deje nser la longitud de la secuencia que se generará.

1er algoritmo:

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

2do algoritmo:

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

Mi pregunta es, ¿el algoritmo 1 produce secuencias que son tan buenas como las generadas por el algoritmo 2?

ultrajohn
fuente
Rkn[a,b]rand_array <- replicate(k, sort(runif(n, a, b))

Respuestas:

18

El primer algoritmo falla gravemente por dos razones:

  1. (ab)/nba<n

  2. na=0b=1(11/n)n1/e37%11/n1100%posibilidad de que el máximo sea en ese intervalo. Para algunos propósitos, esta súper uniformidad es buena, pero en general es un error terrible porque (a) muchas estadísticas se arruinarán pero (b) puede ser muy difícil determinar por qué.

  3. n+1(0,1)1(a,b)

1000n=100

Para conocer muchas más formas (divertidas) de simular variaciones uniformes independientes, consulte Simulación de sorteos de una distribución uniforme mediante sorteos de una distribución normal .

Figura: histogramas

Aquí está el Rcódigo que produjo la figura.

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")
whuber
fuente
¿Qué opinas del algoritmo (basado en estadísticas de orden de rango) en mi respuesta? ;-)
HA SALIDO - Anony-Mousse
@Anony Es una versión menos eficiente de mi algoritmo 3. (El tuyo parece implicar una gran cantidad de cambios de escala innecesarios). Generas las variaciones exponenciales tomando registros de uniformes, lo cual es estándar.
whuber
6

El primer algoritmo produce números demasiado espaciados uniformemente

Ver también series de baja discrepancia .

[0;1]

(Como se ha señalado, esto puede ser por ejemplo una propiedad deseada para la estratificación. Serie de baja discrepancia como Halton y Sobel no tienen sus casos de uso.)

Un enfoque adecuado pero costoso (para valores reales)

... es usar números aleatorios distribuidos en beta. La estadística de orden de rango de la distribución uniforme es beta distribuida. Puede usar esto para dibujar aleatoriamente el más pequeño , luego el segundo más pequeño, ... repita.

[0;1]Beta[1,n]n1XBeta[n,1]ln(1X)Exponential[n]ln(U[0;1])n

ln(1x)=ln(1u)n1x=u1nx=1u1n

Lo que produce el siguiente algoritmo:

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

Puede haber inestabilidades numéricas involucradas, y la computación powy una división para cada objeto pueden resultar más lentas que la clasificación.

Para valores enteros puede que necesite usar una distribución diferente.

Ordenar es increíblemente barato, así que solo úsalo

O(nlogn)

HA SALIDO - Anony-Mousse
fuente
1
Puede haber razones para evitar la clasificación. Una es cuando desea generar una gran cantidad de variantes aleatorias, tantas que una rutina de clasificación estándar no puede manejarlas.
whuber
Creo que los problemas numéricos con las sumas que usan las matemáticas de punto flotante se convirtieron en un problema mucho antes. (¡Y los problemas con los patrones cíclicos en números pseudoaleatorios!) Es bastante fácil escalar el enfoque de clasificación a terabytes y exabytes en sistemas distribuidos.
HA SALIDO - Anony-Mousse
1012
Ok, no tener que almacenarlos es un argumento. Pero entonces necesitarás mi enfoque, tu variante 3 usando la suma acumulativa no funcionará.
HA SALIDO - Anony-Mousse
Ese es un excelente punto. ¡Ahora veo la virtud de los cálculos adicionales! (+1)
whuber
5

También depende de lo que esté haciendo con los números aleatorios. Para los problemas de integración numérica, el método uno (cuando se corrige quitando el operador de piso) produciría un conjunto de puntos superior. Lo que está haciendo es una forma de muestreo estratificado y tiene la ventaja de que evita la aglomeración. Es imposible obtener todos sus valores en el rango 0- (ba) / n, por ejemplo. Dicho esto para otras aplicaciones, esto podría ser muy malo, depende de lo que quieras hacer con él.

usuario67054
fuente
2
+1 Creo que esta es una contribución útil a la pregunta, especialmente al caracterizar el Algoritmo 1 en términos de estratificación.
whuber