math - sorteos - Generación aleatoria de números racionales

numeros aleatorios google

Question

Los racionales son enumerables. Por ejemplo, este código encuentra k-th racional en el intervalo abierto 0..1, ordenando que {n1, d1} sea ​​anterior a {n2, d2} if (d1<d2 || (d1==d2 && n1<n2)) suponiendo que {n,d} sea ​​coprime.

RankedRational[i_Integer?Positive] := 
 Module[{sum = 0, eph = 1, den = 1},
  While[sum < i, sum += (eph = EulerPhi[++den])];
  Select[Range[den - 1], CoprimeQ[#, den] &][[i - (sum - eph)]]/den
  ]

In[118]:= Table[RankedRational[i], {i, 1, 11}]

Out[118]= {1/2, 1/3, 2/3, 1/4, 3/4, 1/5, 2/5, 3/5, 4/5, 1/6, 5/6}

Ahora me gustaría generar racionales aleatorios, dado un límite superior en el orden del denominador uniformemente, de modo que, para un denominador suficientemente grande, los racionales se distribuyan uniformemente a lo largo del intervalo de la unidad.

Intuitivamente, uno podría elegir entre todos los racionales con pequeños denominadores con pesos iguales:

RandomRational1[maxden_, len_] := 
 RandomChoice[(Table[
     i/j, {j, 2, maxden}, {i, 
      Select[Range[j - 1], CoprimeQ[#, j] &]}] // Flatten), len]

¿Se pueden generar racionales aleatorios con esta distribución de manera más eficiente, sin construirlos todos? No se necesita mucho para que esta mesa se vuelva enorme.

In[197]:= Table[RankedRational[10^k] // Denominator, {k, 2, 10}]

Out[197]= {18, 58, 181, 573, 1814, 5736, 18138, 57357, 181380}

¿O tal vez es posible generar racionalmente con un denominador limitado que tenga una distribución diferente de "sentir como uniforme"?

EDITAR Este es el código de Mathematica que ejecuta la generación de aceptación / rechazo sugerida por btilly.

Clear[RandomFarey];
RandomFarey[n_, len_] := Module[{pairs, dim = 0, res, gcds},
  Join @@ Reap[While[dim < len,
      gcds = cfGCD[pairs = cfPairs[n, len - dim]];
      pairs = Pick[pairs, gcds, 1];
      If[pairs =!= {}, 
       dim += Length@Sow[res = pairs[[All, 1]]/pairs[[All, 2]]]];
      ]][[2, -1]]
  ]

La siguiente función compilada genera pares de números enteros {i,j} tales que 1<=i < j<=n :

cfPairs = 
  Compile[{{n, _Integer}, {len, _Integer}}, 
   Table[{i, RandomInteger[{i + 1, n}]}, {i, 
     RandomChoice[2 (n - Range[n - 1])/(n (n - 1.0)) -> Range[n - 1], 
      len]}]];

y la siguiente función compilada calcula gcd. Supone que la entrada es un par de enteros positivos.

cfGCD = Compile[{{prs, _Integer, 1}}, Module[{a, b, p, q, mod},
    a = prs[[1]]; b = prs[[2]]; p = Max[a, b]; q = Min[a, b]; 
    While[q > 0, mod = Mod[p, q]; p = q; q = mod]; p], 
   RuntimeAttributes -> Listable];

Entonces

In[151]:= data = RandomFarey[12, 10^6]; // AbsoluteTiming

Out[151]= {1.5423084, Null}

In[152]:= cdf = CDF[EmpiricalDistribution[data], x];

In[153]:= Plot[{cdf, x}, {x, 0, 1}, ImageSize -> 300]

Answer 1

Aquí hay algunos pensamientos al azar sobre el problema que plantea. No he revisado cuidadosamente los cálculos, así que podría estar apagado por 1 aquí o allí. Pero representa el tipo de razonamiento que seguiría.

Consideremos solo las fracciones en el intervalo (0,1). Es mucho más fácil de esa manera. Podemos tratar más adelante con 1/1 y fracciones impropias.

El árbol Stern-Brocot enumera de manera única cada fracción común positiva reducida (y por lo tanto cada número racional positivo menor o igual a uno) una vez, en orden y en forma reducida, como un nodo en el árbol. En este árbol binario, cualquier nodo y, por lo tanto, cualquier fracción se puede alcanzar mediante una secuencia finita de giros de izquierda a derecha comenzando desde el nivel superior (por conveniencia, llamémoslo nivel -1), que contiene 0/1 y 1/0. [Sí, 1/0. ¡Eso no es un error de impresión!]

Dado un denominador, k, necesitarías tomar como máximo k giros para alcanzar cualquier fracción reducida j / k, donde j es menor que k. Por ejemplo, si el denominador fuera 101, todas las fracciones posibles con un denominador de 101 o menos estarán en el árbol en algún lugar entre el Nivel 1 (que contiene 1/1) y el Nivel 101 (que contiene 1/101 en la posición más a la izquierda).

Supongamos que tenemos un generador de números que genera 0 y 1. (No me pregunten cómo hacerlo, no tengo ni idea). Lef decide arbitrariamente que Izquierda = 0 y Derecha = 1.

Supongamos que tenemos otro generador de números que puede generar números enteros aleatoriamente entre 1 y n. Supongamos además que el primer número generado es 0, es decir. gire a la izquierda: esto garantiza que la fracción caerá en el intervalo (0,1).

Seleccione el máximo denominador, k. Generar aleatoriamente un número, m, entre 1 y k. Luego genere una lista aleatoria de R y L. Atraviesa (es decir, desciende) el árbol de Stern-Brocot, siguiendo la lista de giros. Deténgase cuando llegue a la fracción de destino.

Si esa fracción tiene un denominador igual o menor que k, deténgase, ese es su número.

Si el denominador es mayor que k, ascienda el árbol (a lo largo del mismo camino que descendió) hasta que llegue a una fracción con un denominador no mayor que k.

No sé que la generación de números sea verdaderamente aleatoria. Ni siquiera sabría cómo decirlo. Pero por lo que es worthe, no detecto ninguna fuente obvia de sesgo.

Answer 2

Con un límite en el denominador, los racionales no están uniformemente distribuidos (1/2 se separa de todo lo demás por un buen hueco, por ejemplo.

Dicho esto, sería algo así como

In[300]:= Rationalize[RandomReal[1, 10], 0.001]

Out[300]= {17/59, 45/68, 11/31, 9/16, 1/17, 13/22, 7/10, 1/17, 5/21, 8/39}

¿trabajo para ti?

Answer 3

Recomiendo mirar el juego "adivinar el número" para obtener números racionales arbitrarios. para obtener algo de inspiración sobre su problema subyacente.

Si su objetivo es ser aproximadamente uniforme lo antes posible, y no le importa elegir diferentes racionales con diferentes probabilidades, el siguiente algoritmo debería ser eficiente.

lower = fractions.Fraction(0)
upper = fractions.Fraction(1)

while lower < upper:
    mid = (upper + lower)/2
    if 0 == random_bit():
        upper = largest_rational_under(mid, denominator_bound)
    else:
        lower = smallest_rational_over_or_equal(mid, denominator_bound)

Tenga en cuenta que ambas funciones auxiliares se pueden calcular caminando el Árbol Stern-Brocot hacia la mitad. También tenga en cuenta que, con algunas modificaciones menores, puede transformar fácilmente esto en un algoritmo iterativo que escupe una secuencia de números racionales, y eventualmente convergerá con la misma probabilidad en cualquier parte del intervalo. Considero que esa propiedad es amable.

Si quiere la distribución exacta que especificó originalmente, y rand(n) le da un entero aleatorio de 1 a n , entonces el siguiente pseudocódigo funcionará para denominador enlazado n :

Try:
    k = rand(n * (n+1) / 2)
    do binary search for largest j with j * (j-1) / 2 < k
    i = k - (j * (j-1) / 2)
    if (i, j) are not relatively prime:
        redo Try
answer = i/j

En promedio, para n grande, deberá Try aproximadamente 2.55 veces. Entonces en la práctica esto debería ser bastante eficiente.