Entendiendo el sesgo de módulo | Julian a través del lente

Introducción

Suele decirse que este código:

unsigned int randomNumber = rand() % k;

es una mala idea, al menos si esperás obtener una distribución uniforme. Voy a intentar explorar este tema de manera más formal de lo que he visto hasta ahora.

La razón por la que es problemático es bastante elemental y fácil de entender: imaginá que tenés un generador aleatorio que produce valores entre $0$ y $9$ (es decir, RAND_MAX es $10$ ), y querés valores entre $0$ y $2$ (así que establecerías $k = 3$ ). Entonces, tenemos el siguiente mapeo:

rand()	Random Number
0, 3, 6, 9	0
1, 4, 7	1
2, 5, 8	2

Así, en la práctica, $P (randomNumber = 2) = \frac{3}{10}$ ; en cambio, $P (randomNumber = 0) = \frac{4}{10}$ . Esencialmente, el problema radica en que como $k$ no divide de manera exacta a RAND_MAX, esto hace que hasta $1$ valores adicionales caigan en los primeros números del rango. Para hacerte una idea, podés jugar con el gráfico de abajo, que muestra las probabilidades precisas para distintos valores de RAND_MAX y $k$ :

En la práctica, el efecto del parámetro $k$ se vuelve más pronunciado cuando el valor de $M % k$ crece; hasta que se hace demasiado grande (es decir, $k = M - 1$ ), momento en el que es casi como si el problema no existiera. De manera más formal, podemos analizar esto definiendo $X$ como una variable aleatoria discreta tal que $R_{X} = {0, ..., M - 1}$ , donde $M - 1 \in N$ es el número más alto alcanzable por el generador aleatorio y $P (X = i) = \frac{1}{M}$ .

Sea $Y = X % k$ . Lo que queremos es encontrar cuál es la función de distribución de probabilidad de $Y$ , ya que eso nos permitirá entender con precisión lo que está ocurriendo. En primer lugar, notemos que $R_{Y} = {0, ..., k - 1}$ . De ahí, la probabilidad:

P (Y = i) = P (\cup_{j \in J_{i}} (X = j)) = j \in J_{i} \sum P (X = j) = \frac{# J _{i}}{M}

Donde $J_{i} = {x \in R_{X} : x \equiv i (mod k)}$ . En concreto, la probabilidad de que $Y$ pertenezca a una clase de equivalencia particular es el número de elementos de esa clase dividido por el total de elementos de los que podemos elegir. Entonces, lo siguiente es calcular $# J_{i}$ .

Los elementos que pertenecen a $J_{i}$ tienen todos la misma forma: $q k + i$ , con $q \in Z$ . Sin embargo, están acotados superiormente por $M - 1$ e inferiormente por $i$ :

i \leq q k + i \leq M - 1 ⟹ 0 \leq q \leq ⌊ \frac{M - 1 - i}{k} ⌋

Por lo tanto, $# J_{i} = ⌊ \frac{M - 1 - i}{k} ⌋ + 1$ . Lo que significa que nuestra función de distribución de probabilidad es

P (Y = i) = \frac{⌊ \frac{M - 1 - i}{k} ⌋ + 1}{M}

Queremos verificar que esto suma $1$ , ya que eso confirmaría que es una distribución de probabilidad. Usaré las identidades $⌈ \frac{a}{b} ⌉ = ⌊ \frac{a - 1}{b} ⌋ + 1$ y $a = \sum_{i = 0}^{b - 1} ⌈ \frac{a - i}{b} ⌉$ ; ambas demostraciones pueden encontrarse en {1}:

i \in R_{Y} \sum P (Y = i) = i = 0 \sum k - 1 (\frac{⌊ \frac{M - 1 - i}{k} ⌋ + 1}{M}) = \frac{1}{M} i = 0 \sum k - 1 (⌊ \frac{M - 1 - i}{k} ⌋ + 1) = \frac{1}{M} i = 0 \sum k - 1 ⌈ \frac{M - i}{k} ⌉ = \frac{1}{M} M = 1

El hecho de que la distribución dependa del valor de $i$ ya es un problema: buscábamos una distribución uniforme y, por lo tanto, esperábamos que $P (Y = i) = \frac{1}{k}$ . Si graficás la función, verás exactamente el mismo gráfico con el que estuviste jugando arriba.

¿Qué tan grave es?

Qué tan grave sea depende de lo que estés haciendo. Una pregunta más interesante es cómo cuantificar esa gravedad. Una forma de hacerlo es medir qué tan diferente es la distribución de $Y$ respecto a la distribución de $Z$ ; hay una infinidad de formas de hacer esto , pero voy a usar la divergencia de Kullback-Leibler .

Entonces, supongamos que $Z$ es una variable aleatoria con distribución uniforme tal que $R_{Z} = {0, .., k - 1}$ (es decir, la distribución real que queríamos obtener). Dado algún $M$ , sabemos que $k$ puede tomar valores entre $1 \leq k \leq M - 1$ ; para cada uno de esos valores, tenemos una distribución de probabilidad, llamémosla $Y$ , y podemos calcular $KL (Z ∣∣ Y)$ :

Si explorás distintos valores de $M$ , quedará claro que $KL (Z ∣∣ Y)$ solo alcanza 0 cuando $k$ divide exactamente a $M$ ; y por lo tanto algunos valores de $M$ nunca llegan a un punto que alcance $0$ ; específicamente, estos son los números primos.

Además, el espacio entre los ceros se va haciendo cada vez mayor a medida que $M$ crece. Esto se debe básicamente a que los divisores pares de $M$ están más separados cuanto más grandes son; es fácil verlo mirando la descomposición prima de $M = \prod_{i = 1}^{b} p_{i}^{r_{i}}$ : todo divisor par puede verse exactamente como esa descomposición, solo que con un valor menor o igual para cada $r_{i}$ , lo que significa que los aumentos en el valor son multiplicativos.

La última clave para entender esto son las parábolas. ¿Por qué aparece una parábola invertida entre cada divisor par? Es simplemente por cómo funciona el módulo:

Al aumentar $k$ , un valor más queda sesgado. Esto significa que la distribución generada se aleja de la distribución objetivo $Z$ .
Cuando la cantidad de valores sesgados ( $M % k$ ) es aproximadamente $\frac{k}{2}$ , alcanzamos el punto más alto de la parábola.
Después de esto solo puede decrecer: a medida que más valores quedan sesgados, la distribución se parece cada vez más a una distribución uniforme nuevamente.

¿Cómo puedo corregirlo?

Quizás la solución más sencilla es hacer:

unsigned int randomNumber = k;
while (randomNumber >= k) randomNumber = rand();

Esto funciona exactamente como se espera: el número producido siempre está entre $0$ y $k - 1$ , y se distribuye de manera uniforme en el rango. Para verlo, imaginá que definimos $T_{i}$ como variables independientes e idénticamente distribuidas con distribución uniforme entre $0$ y $M - 1$ ; entonces podemos calcular la función de distribución de probabilidad de $T$ :

P (T = i) = r = 1 \sum \infty P (\land_{j = 1}^{r - 1} T_{j} \geq k \land T_{r} = i) = r = 1 \sum \infty P (T_{r} = i) j = 1 \prod r - 1 P (T_{j} \geq k) = r = 1 \sum \infty \frac{1}{M} (1 - \frac{k}{M})^{r - 1} = \frac{1}{M} r = 0 \sum \infty (1 - \frac{k}{M})^{r} = \frac{1}{M} \frac{1}{1 - ( 1 - \frac{k}{M} )} = \frac{1}{k}

¡Lo cual ahora luce exactamente como queríamos! Quizás la siguiente pregunta interesante es cuántas veces llamaremos a rand() hasta encontrar un número que satisfaga la restricción.

Esto tiene una respuesta sencilla: la probabilidad de obtener un número aleatorio en el rango es $\frac{k}{M}$ , y los números aleatorios generados son independientes entre sí; lo que significa que el número de iteraciones $N$ es una variable aleatoria geométrica . Por lo tanto, tenemos que $N \sim Geometric (\frac{k}{M})$ , y así podemos esperar tomar $E [N] = \frac{M}{k}$ iteraciones hasta generar un número aleatorio.

La esperanza señala un problema en nuestro código: si $M ≫ k$ , lo que ocurre frecuentemente en implementaciones prácticas donde $M$ es el número máximo representable, entonces el número de iteraciones va a ser alto.

¿Y ahora qué?

Bien, el problema con la última idea es que descartamos la mayor parte de los números que generamos, mientras que el problema con la primera idea es que existe un conjunto de números que sesgan nuestro generador. La siguiente idea es intentar descartar únicamente los números que producen el sesgo.

Podemos escribir $M$ como $⌊ \frac{M}{k} ⌋ k + (M % k)$ ; si aplicáramos la técnica del módulo sobre un generador con salida entre $0$ y $⌊ \frac{M}{k} ⌋ k - 1$ , obtendríamos una distribución uniforme, como vimos antes (porque el rango sería exactamente divisible). En efecto, todo lo que necesitamos hacer es descartar los $M % k$ valores al final o al comienzo del rango, y luego aplicar el módulo sobre el resultado; esto es prácticamente una «fusión» de las dos ideas anteriores. Quedaría algo así:

unsigned int threshold = M - (M % k);
unsigned int randomNumber;
do {
    randomNumber = rand();
} while (randomNumber >= threshold)
randomNumber = randomNumber % k;

En realidad no necesitamos demostrar que esto es correcto: es evidente que las primeras cuatro líneas muestrean de manera uniforme un número del rango ${0, ..., ⌊ \frac{M}{k} ⌋ k - 1}$ (a partir de la demostración del esquema anterior), y luego aplicar el módulo sobre ese valor lo convierte en una distribución uniforme con el rango deseado (gracias a la demostración del primer esquema).

El único cambio interesante aquí es la cantidad de llamadas a rand(). Ahora, la probabilidad de que un número esté en el rango deseado para el bucle es $\frac{M - ( M % k )}{M} = 1 - \frac{M % k}{M}$ , y por lo tanto el número de iteraciones es $N \sim Geometric (1 - \frac{M % k}{M})$ . Podemos ahora ver la cantidad esperada de iteraciones para distintos valores de $k$ y $M$ :

En efecto, esto es prácticamente lo que esperaríamos que ocurriera: cuanto más grande es $M % k$ , más valores desperdiciamos y, por lo tanto, más iteraciones se necesitan. Podemos demostrar, sin embargo, que el número esperado de iteraciones está acotado por $2$ :

E [N] = \frac{1}{1 - \frac{M % k}{M}} = \frac{1}{1 - \frac{M - ⌊ \frac{M}{k} ⌋ k}{M}} = \frac{M}{⌊ \frac{M}{k} ⌋ k} = 1 + \frac{M % k}{⌊ \frac{M}{k} ⌋ k}

Donde hemos usado repetidamente el hecho de que $M = ⌊ \frac{M}{k} ⌋ k + (M % k)$ . Notemos que $\frac{M % k}{⌊ \frac{M}{k} ⌋ k}$ es estrictamente menor que $1$ , porque si fuera $1$ o más, podríamos incrementar el valor de $⌊ \frac{M}{k} ⌋$ , lo que contradice la definición del propio operador de división. Por lo tanto, toda la expresión es menor que 2. Esto significa que nuestro algoritmo tiene un tiempo de ejecución esperado de $O (1)$ , asumiendo que la operación rand es $O (1)$ .

Quizás la última pregunta interesante respecto a este algoritmo es cuál es la probabilidad de hacer más de $t$ iteraciones antes de producir un número aleatorio. Trabajando con la función de distribución acumulada de la distribución geométrica se obtiene:

P (N > t) = (\frac{M % k}{M})^{t} = (1 - ⌊ \frac{M}{k} ⌋ \frac{k}{M})^{t}

Cuanto más cerca esté $⌊ \frac{M}{k} ⌋ \frac{k}{M}$ de $1$ , menos probable será necesitar más iteraciones. Por supuesto, este término se parece bastante a la cantidad esperada de iteraciones, solo que invertido y dividido por dos:

Es bastante fácil ver que este término es siempre mayor que $\frac{1}{2}$ . Sustituyendo esto en la fórmula anterior, obtenemos que $P (N > t) \leq (\frac{1}{2})^{t}$ .

{1} Knuth, Donald. 1994. Concrete Mathematics.

Introducción

¿Qué tan grave es?

¿Cómo puedo corregirlo?

¿Y ahora qué?

Relacionados