c++ - rango - metodo random java

Forma rápida de generar bits pseudoaleatorios con una probabilidad dada de 0 o 1 para cada bit (10)

Normalmente, un generador de números aleatorios devuelve una secuencia de bits para los cuales la probabilidad de observar un 0 o un 1 en cada posición es igual (es decir, 50%). Llamemos a esto un PRNG imparcial.

Necesito generar una cadena de bits pseudoaleatorios con la siguiente propiedad: la probabilidad de ver un 1 en cada posición es p (es decir, la probabilidad de ver un 0 es 1-p). El parámetro p es un número real entre 0 y 1; en mi problema sucede que tiene una resolución de 0.5%, es decir, puede tomar los valores 0%, 0.5%, 1%, 1.5%, ..., 99.5%, 100%.

Tenga en cuenta que p es una probabilidad y no una fracción exacta. El número real de bits establecido en 1 en una secuencia de n bits debe seguir la distribución binomial B (n, p).

Hay un método ingenuo que puede usar un PRNG imparcial para generar el valor de cada bit (pseudocódigo):

generate_biased_stream(n, p): result = [] for i in 1 to n: if random_uniform(0, 1) < p: result.append(1) else: result.append(0) return result

Tal implementación es mucho más lenta que la que genera un flujo imparcial, ya que llama a la función de generador de números aleatorios una vez por cada bit; mientras que un generador de flujo imparcial lo llama una vez por tamaño de palabra (por ejemplo, puede generar 32 o 64 bits aleatorios con una sola llamada).

Quiero una implementación más rápida, incluso sacrifica un poco la aleatoriedad. Una idea que se me ocurre es calcular previamente una tabla de búsqueda: para cada uno de los 200 valores posibles de p, calcule los valores de C de 8 bits utilizando el algoritmo más lento y guárdelos en una tabla. Luego, el algoritmo rápido elegiría uno de estos al azar para generar 8 bits sesgados.

Una parte posterior del cálculo del sobre para ver cuánta memoria se necesita: C debe ser al menos 256 (el número de valores posibles de 8 bits), probablemente más para evitar efectos de muestreo; digamos 1024. Quizás el número debería variar dependiendo de p, pero seamos simples y digamos que el promedio es 1024. Dado que hay 200 valores de p => el uso total de memoria es 200 KB. Esto no está mal y podría caber en la caché L2 (256 KB). Todavía necesito evaluarlo para ver si hay efectos de muestreo que introducen sesgos, en cuyo caso habrá que aumentar C.

Una deficiencia de esta solución es que puede generar solo 8 bits a la vez, incluso con mucho trabajo, mientras que un PRNG imparcial puede generar 64 a la vez con solo unas pocas instrucciones aritméticas.

Me gustaría saber si hay un método más rápido, basado en operaciones de bits en lugar de tablas de búsqueda. Por ejemplo, modificar el código de generación de números aleatorios directamente para introducir un sesgo para cada bit. Esto lograría el mismo rendimiento que un PRNG imparcial.

Editar 5 de marzo

Gracias a todos por sus sugerencias, tengo muchas ideas y sugerencias interesantes. Aquí están los mejores:

Cambie los requisitos del problema para que p tenga una resolución de 1/256 en lugar de 1/200. Esto permite usar bits de manera más eficiente y también brinda más oportunidades de optimización. Creo que puedo hacer este cambio.
Utilice la codificación aritmética para consumir eficientemente bits de un generador imparcial. Con el cambio de resolución anterior, esto se vuelve mucho más fácil.
Algunas personas sugirieron que los PRNG son muy rápidos, por lo que el uso de la codificación aritmética podría hacer que el código sea más lento debido a la sobrecarga introducida. En cambio, siempre debería consumir el peor número de bits y optimizar ese código. Ver los puntos de referencia a continuación.
@rici sugirió usar SIMD. Esta es una buena idea, que funciona solo si siempre consumimos un número fijo de bits.

Puntos de referencia (sin decodificación aritmética)

Nota: como muchos de ustedes han sugerido, cambié la resolución de 1/200 a 1/256.

Escribí varias implementaciones del método ingenuo que simplemente toma 8 bits imparciales aleatorios y genera 1 bit sesgado:

Sin SIMD
Con SIMD usando la biblioteca de clases de vectores de Agner Fog, como lo sugiere @rici
Con SIMD usando intrínsecos

Utilizo dos generadores de números seudoaleatorios imparciales:

xorshift128plus
Ranvec1 (Mersenne Twister-like) de la biblioteca de Agner Fog.

También mido la velocidad del PRNG imparcial para comparar. Aquí están los resultados:

RNG: Ranvec1(Mersenne Twister for Graphics Processors + Multiply with Carry) Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline) Gbps/s: 16.081 16.125 16.093 [Gb/s] Number of ones: 536,875,204 536,875,204 536,875,204 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency Gbps/s: 0.778 0.783 0.812 [Gb/s] Number of ones: 104,867,269 104,867,269 104,867,269 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=vectorclass Gbps/s: 2.176 2.184 2.145 [Gb/s] Number of ones: 104,859,067 104,859,067 104,859,067 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=intrinsics Gbps/s: 2.129 2.151 2.183 [Gb/s] Number of ones: 104,859,067 104,859,067 104,859,067 Theoretical : 104,857,600

SIMD aumenta el rendimiento en un factor de 3 en comparación con el método escalar. Es 8 veces más lento que el generador imparcial, como se esperaba.

El generador polarizado más rápido alcanza 2.1 Gb / s.

RNG: xorshift128plus Method: Unbiased with 1/1 efficiency (incorrect, baseline) Gbps/s: 18.300 21.486 21.483 [Gb/s] Number of ones: 536,867,655 536,867,655 536,867,655 Theoretical : 104,857,600 Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline) Gbps/s: 22.660 22.661 24.662 [Gb/s] Number of ones: 536,867,655 536,867,655 536,867,655 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency Gbps/s: 1.065 1.102 1.078 [Gb/s] Number of ones: 104,868,930 104,868,930 104,868,930 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=vectorclass Gbps/s: 4.972 4.971 4.970 [Gb/s] Number of ones: 104,869,407 104,869,407 104,869,407 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=intrinsics Gbps/s: 4.955 4.971 4.971 [Gb/s] Number of ones: 104,869,407 104,869,407 104,869,407 Theoretical : 104,857,600

Para xorshift, SIMD aumenta el rendimiento en un factor de 5 en comparación con el método escalar. Es 4 veces más lento que el generador imparcial. Tenga en cuenta que esta es una implementación escalar de xorshift.

El generador polarizado más rápido alcanza 4.9 Gb / s.

RNG: xorshift128plus_avx2 Method: Unbiased with 1/1 efficiency (incorrect, baseline) Gbps/s: 18.754 21.494 21.878 [Gb/s] Number of ones: 536,867,655 536,867,655 536,867,655 Theoretical : 104,857,600 Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline) Gbps/s: 54.126 54.071 54.145 [Gb/s] Number of ones: 536,874,540 536,880,718 536,891,316 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency Gbps/s: 1.093 1.103 1.063 [Gb/s] Number of ones: 104,868,930 104,868,930 104,868,930 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=vectorclass Gbps/s: 19.567 19.578 19.555 [Gb/s] Number of ones: 104,836,115 104,846,215 104,835,129 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=intrinsics Gbps/s: 19.551 19.589 19.557 [Gb/s] Number of ones: 104,831,396 104,837,429 104,851,100 Theoretical : 104,857,600

Esta implementación utiliza AVX2 para ejecutar 4 generadores xorshift imparciales en paralelo.

El generador polarizado más rápido alcanza 19.5 Gb / s.

Puntos de referencia para la decodificación aritmética

Las pruebas simples muestran que el código de decodificación aritmética es el cuello de botella, no el PRNG. Así que solo estoy comparando el PRNG más caro.

RNG: Ranvec1(Mersenne Twister for Graphics Processors + Multiply with Carry) Method: Arithmetic decoding (floating point) Gbps/s: 0.068 0.068 0.069 [Gb/s] Number of ones: 10,235,580 10,235,580 10,235,580 Theoretical : 10,240,000 Method: Arithmetic decoding (fixed point) Gbps/s: 0.263 0.263 0.263 [Gb/s] Number of ones: 10,239,367 10,239,367 10,239,367 Theoretical : 10,240,000 Method: Unbiased with 1/1 efficiency (incorrect, baseline) Gbps/s: 12.687 12.686 12.684 [Gb/s] Number of ones: 536,875,204 536,875,204 536,875,204 Theoretical : 104,857,600 Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline) Gbps/s: 14.536 14.536 14.536 [Gb/s] Number of ones: 536,875,204 536,875,204 536,875,204 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency Gbps/s: 0.754 0.754 0.754 [Gb/s] Number of ones: 104,867,269 104,867,269 104,867,269 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=vectorclass Gbps/s: 2.094 2.095 2.094 [Gb/s] Number of ones: 104,859,067 104,859,067 104,859,067 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=intrinsics Gbps/s: 2.094 2.094 2.095 [Gb/s] Number of ones: 104,859,067 104,859,067 104,859,067 Theoretical : 104,857,600

El método simple de punto fijo alcanza 0.25 Gb / s, mientras que el método escalar ingenuo es 3 veces más rápido y el método SIMD ingenuo es 8 veces más rápido. Puede haber formas de optimizar y / o paralelizar aún más el método de decodificación aritmética, pero debido a su complejidad, he decidido detenerme aquí y elegir la implementación SIMD ingenua.

Gracias a todos por la ayuda.

Que hace

Esta implementación realiza una sola llamada al módulo de kernel de dispositivo aleatorio a través de la interfaz del archivo de caracteres especiales "/ dev / urandom" para obtener la cantidad de datos aleatorios necesarios para representar todos los valores en una resolución dada. La resolución máxima posible es 1/256 ^ 2, de modo que 0.005 se puede representar mediante:

328/256 ^ 2,

es decir:

resolución: 256 * 256

x: 328

con error 0.000004883.

Como hace eso

La implementación calcula el número de bits bits_per_byte que es el número de bits distribuidos uniformemente necesarios para manejar la resolución dada, es decir, representan todos los @resolution valores. Luego realiza una sola llamada al dispositivo de aleatorización ("/ dev / urandom" si URANDOM_DEVICE está definido; de lo contrario, utilizará ruido adicional de los controladores del dispositivo a través de la llamada a "/ dev / random" que puede bloquearse si no hay suficiente entropía en bits) para obtener el número requerido de bytes distribuidos uniformemente y rellena una matriz rnd_bytes de bytes aleatorios. Finalmente, lee el número de bits necesarios por cada muestra de Bernoulli de cada bytes bytes_por_byte bytes de la matriz rnd_bytes y compara el valor entero de estos bits con la probabilidad de éxito en un único resultado de Bernoulli dado por x/resolution . Si el valor alcanza, es decir, cae en el segmento de x/resolution longitud que elegimos arbitrariamente como segmento [0, x / resolución), luego notamos el éxito e insertamos 1 en la matriz resultante.

Leer desde un dispositivo aleatorio:

/* if defined use /dev/urandom (will not block), * if not defined use /dev/random (may block)*/ #define URANDOM_DEVICE 1 /* * @brief Read @outlen bytes from random device * to array @out. */ int get_random_samples(char *out, size_t outlen) { ssize_t res; #ifdef URANDOM_DEVICE int fd = open("/dev/urandom", O_RDONLY); if (fd == -1) return -1; res = read(fd, out, outlen); if (res < 0) { close(fd); return -2; } #else size_t read_n; int fd = open("/dev/random", O_RDONLY); if (fd == -1) return -1; read_n = 0; while (read_n < outlen) { res = read(fd, out + read_n, outlen - read_n); if (res < 0) { close(fd); return -3; } read_n += res; } #endif /* URANDOM_DEVICE */ close(fd); return 0; }

Complete el vector de muestras de Bernoulli:

/* * @brief Draw vector of Bernoulli samples. * @details @x and @resolution determines probability * of success in Bernoulli distribution * and accuracy of results: p = x/resolution. * @param resolution: number of segments per sample of output array * as power of 2: max resolution supported is 2^24=16777216 * @param x: determines used probability, x = [0, resolution - 1] * @param n: number of samples in result vector */ int get_bernoulli_samples(char *out, uint32_t n, uint32_t resolution, uint32_t x) { int res; size_t i, j; uint32_t bytes_per_byte, word; unsigned char *rnd_bytes; uint32_t uniform_byte; uint8_t bits_per_byte; if (out == NULL || n == 0 || resolution == 0 || x > (resolution - 1)) return -1; bits_per_byte = log_int(resolution); bytes_per_byte = bits_per_byte / BITS_PER_BYTE + (bits_per_byte % BITS_PER_BYTE ? 1 : 0); rnd_bytes = malloc(n * bytes_per_byte); if (rnd_bytes == NULL) return -2; res = get_random_samples(rnd_bytes, n * bytes_per_byte); if (res < 0) { free(rnd_bytes); return -3; } i = 0; while (i < n) { /* get Bernoulli sample */ /* read byte */ j = 0; word = 0; while (j < bytes_per_byte) { word |= (rnd_bytes[i * bytes_per_byte + j] << (BITS_PER_BYTE * j)); ++j; } uniform_byte = word & ((1u << bits_per_byte) - 1); /* decision */ if (uniform_byte < x) out[i] = 1; else out[i] = 0; ++i; } free(rnd_bytes); return 0; }

Uso:

int main(void) { int res; char c[256]; res = get_bernoulli_samples(c, sizeof(c), 256*256, 328); /* 328/(256^2) = 0.0050 */ if (res < 0) return -1; return 0; }

Código completo , results .

Desde un punto de vista teórico de la información, una secuencia de bits sesgada (con p != 0.5 ) contiene menos información que una secuencia imparcial, por lo que en teoría debería tomar (en promedio) menos de 1 bit de la entrada imparcial para producir un solo bit de la secuencia de salida sesgada. Por ejemplo, la entropy de una variable aleatoria de Bernoulli con p = 0.1 es -0.1 * log2(0.1) - 0.9 * log2(0.9) bits, que es alrededor de 0.469 bits. Eso sugiere que para el caso p = 0.1 deberíamos poder producir un poco más de dos bits del flujo de salida por bit de entrada imparcial.

A continuación, doy dos métodos para producir los bits sesgados. Ambos logran una eficiencia cercana a la óptima, en el sentido de que requieren la menor cantidad posible de bits insesgados de entrada.

Método 1: codificación (des) aritmética

Un método práctico es decodificar su flujo de entrada imparcial utilizando (des) codificación aritmética , como ya se describió en la respuesta de alexis . Para este caso simple, no es difícil codificar algo. Aquí hay un pseudocódigo no optimizado ( tos, Python ) que hace esto:

import random def random_bits(): """ Infinite generator generating a stream of random bits, with 0 and 1 having equal probability. """ global bit_count # keep track of how many bits were produced while True: bit_count += 1 yield random.choice([0, 1]) def bernoulli(p): """ Infinite generator generating 1-bits with probability p and 0-bits with probability 1 - p. """ bits = random_bits() low, high = 0.0, 1.0 while True: if high <= p: # Generate 1, rescale to map [0, p) to [0, 1) yield 1 low, high = low / p, high / p elif low >= p: # Generate 0, rescale to map [p, 1) to [0, 1) yield 0 low, high = (low - p) / (1 - p), (high - p) / (1 - p) else: # Use the next random bit to halve the current interval. mid = 0.5 * (low + high) if next(bits): low = mid else: high = mid

Aquí hay un ejemplo de uso:

import itertools bit_count = 0 # Generate a million deviates. results = list(itertools.islice(bernoulli(0.1), 10**6)) print("First 50:", ''''.join(map(str, results[:50]))) print("Biased bits generated:", len(results)) print("Unbiased bits used:", bit_count) print("mean:", sum(results) / len(results))

Lo anterior da el siguiente resultado de muestra:

First 50: 00000000000001000000000110010000001000000100010000 Biased bits generated: 1000000 Unbiased bits used: 469036 mean: 0.100012

Según lo prometido, hemos generado 1 millón de bits de nuestro flujo sesgado de salida utilizando menos de quinientos mil del flujo imparcial de origen.

Para fines de optimización, al traducir esto a C / C ++ puede tener sentido codificar esto usando aritmética de punto fijo basada en enteros en lugar de punto flotante.

Método 2: algoritmo basado en enteros

En lugar de intentar convertir el método de decodificación aritmética para usar enteros directamente, aquí hay un enfoque más simple. Ya no es una decodificación aritmética, pero no es totalmente ajena, y logra casi la misma relación de bits de salida / sesgo de entrada-imparcial de entrada que la versión de punto flotante anterior. Está organizado de modo que todas las cantidades quepan en un entero de 32 bits sin signo, por lo que debería ser fácil de traducir a C / C ++. El código está especializado para el caso donde p es un múltiplo exacto de 1/200 , pero este enfoque funcionaría para cualquier p que pueda expresarse como un número racional con un denominador razonablemente pequeño.

def bernoulli_int(p): """ Infinite generator generating 1-bits with probability p and 0-bits with probability 1 - p. p should be an integer multiple of 1/200. """ bits = random_bits() # Assuming that p has a resolution of 0.05, find p / 0.05. p_int = int(round(200*p)) value, high = 0, 1 while True: if high < 2**31: high = 2 * high value = 2 * value + next(bits) else: # Throw out everything beyond the last multiple of 200, to # avoid introducing a bias. discard = high - high % 200 split = high // 200 * p_int if value >= discard: # rarer than 1 time in 10 million value -= discard high -= discard elif value >= split: yield 0 value -= split high = discard - split else: yield 1 high = split

La observación clave es que cada vez que alcanzamos el comienzo del ciclo while, el value se distribuye uniformemente entre todos los enteros en [0, high) , y es independiente de todos los bits de salida anteriores. Si le importa la velocidad más que la corrección perfecta, puede deshacerse del discard y el value >= discard rama de value >= discard : eso es justo allí para garantizar que produzcamos 0 y 1 con las probabilidades correctas. Deje esa complicación fuera y obtendrá las probabilidades correctas en su lugar. Además, si hace que la resolución para p igual a 1/256 lugar de 1/200 , las operaciones de división y módulo que pueden llevar mucho tiempo pueden reemplazarse con operaciones de bit.

Con el mismo código de prueba que antes, pero usando bernoulli_int en lugar de bernoulli , obtengo los siguientes resultados para p=0.1 :

First 50: 00000010000000000100000000000000000000000110000100 Biased bits generated: 1000000 Unbiased bits used: 467997 mean: 0.099675

Digamos que la probabilidad de que aparezca 1 es 6,25% (1/16). Hay 16 patrones de bits posibles para un número de 4 bits: 0000,0001, ..., 1110,1111 .

Ahora, solo genere un número aleatorio como solía hacerlo y reemplace cada 1111 en un límite de mordisco con un 1 y cambie todo lo demás a 0 .

Ajuste en consecuencia para otras probabilidades.

Obtendrá un comportamiento teóricamente óptimo, es decir, hará un uso realmente mínimo del generador de números aleatorios y podrá modelar cualquier probabilidad p exactamente, si aborda esto utilizando la codificación aritmética .

La codificación aritmética es una forma de compresión de datos que representa el mensaje como un subintervalo de un rango de números. Proporciona una codificación teóricamente óptima y puede usar un número fraccionario de bits para cada símbolo de entrada.

La idea es la siguiente: imagine que tiene una secuencia de bits aleatorios, que son 1 con probabilidad p . Por conveniencia, usaré q para la probabilidad de que el bit sea cero. ( q = 1-p ). La codificación aritmética se asigna a cada parte de bit del rango de números. Para el primer bit, asigne el intervalo [0, q) si la entrada es 0, y el intervalo [q, 1) si la entrada es 1. Los bits subsiguientes asignan subintervalos proporcionales del rango actual. Por ejemplo, suponga que q = 1/3 La entrada 1 0 0 se codificará así:

Initially [0, 1), range = 1 After 1 [0.333, 1), range = 0.6666 After 0 [0.333, 0.5555), range = 0.2222 After 0 [0.333, 0.407407), range = 0.074074

El primer dígito, 1, selecciona los dos tercios superiores (1-q) del rango; el segundo dígito, 0 , selecciona el tercio inferior de eso , y así sucesivamente. Después del primer y segundo paso, el intervalo se extiende sobre el punto medio; pero después del tercer paso está completamente por debajo del punto medio, por lo que se puede generar el primer dígito comprimido: 0 . El proceso continúa y se agrega un símbolo EOF especial como terminador.

¿Qué tiene esto que ver con tu problema? La salida comprimida tendrá ceros aleatorios y unos con igual probabilidad. Entonces, para obtener bits con probabilidad p , solo pretenda que la salida de su RNG es el resultado de la codificación aritmética como se indicó anteriormente, y aplique el proceso del decodificador. Es decir, lea los bits como si subdividieran el intervalo de línea en piezas cada vez más pequeñas. Por ejemplo, después de leer 01 del RNG, estaremos en el rango [0.25, 0.5) . Siga leyendo bits hasta que se "decodifique" suficiente salida. Como está imitando la descompresión, obtendrá más bits aleatorios de los que ingresó. Debido a que la codificación aritmética es teóricamente óptima, no hay forma posible de convertir la salida RNG en bits más sesgados sin sacrificar la aleatoriedad: está obteniendo el verdadero máximo.

El problema es que no puedes hacer esto en un par de líneas de código, y no conozco una biblioteca a la que pueda señalarte (aunque debe haber alguna que puedas usar). Aún así, es bastante simple. El artículo anterior proporciona código para un codificador y decodificador de uso general, en C. Es bastante sencillo y admite múltiples símbolos de entrada con probabilidades arbitrarias; en su caso, es posible una implementación mucho más simple (como muestra ahora la respuesta de Mark Dickinson ), ya que el modelo de probabilidad es trivial. Para un uso prolongado, se necesitaría un poco más de trabajo para producir una implementación robusta que no haga muchos cálculos de punto flotante para cada bit.

Wikipedia también tiene una discusión interesante sobre la codificación aritmética considerada como cambio de radix, que es otra forma de ver su tarea.

Si está preparado para aproximar p basado en 256 valores posibles, y tiene un PRNG que puede generar valores uniformes en los que los bits individuales son independientes entre sí, entonces puede usar la comparación vectorizada para producir múltiples bits sesgados a partir de un único azar número.

Solo vale la pena hacerlo si (1) te preocupa la calidad de los números aleatorios y (2) es probable que necesites una gran cantidad de bits con el mismo sesgo. El segundo requisito parece estar implícito en la pregunta original, que critica una solución propuesta, como sigue: "Una deficiencia de esta solución es que puede generar solo 8 bits a la vez, incluso con mucho trabajo, mientras que un PRNG imparcial puede generar 64 a la vez con solo unas pocas instrucciones aritméticas ". Aquí, la implicación parece ser que es útil generar un gran bloque de bits sesgados en una sola llamada.

La calidad de los números aleatorios es un tema difícil. Es difícil, si no imposible, medir, y por lo tanto, diferentes personas propondrán diferentes métricas que enfatizan y / o devalúan diferentes aspectos de la "aleatoriedad". En general, es posible cambiar la velocidad de generación de números aleatorios por una "calidad" más baja; si vale la pena hacerlo depende de su aplicación precisa.

Las pruebas más simples posibles de calidad de números aleatorios implican la distribución de valores individuales y la duración del ciclo del generador. Las implementaciones estándar de las funciones random rand y Posix de la biblioteca C generalmente pasarán la prueba de distribución, pero las duraciones de los ciclos no son adecuadas para aplicaciones de larga duración.

Sin embargo, estos generadores suelen ser extremadamente rápidos: la implementación glibc de random requiere solo unos pocos ciclos, mientras que el clásico generador congruencial lineal (LCG) requiere una multiplicación y una suma. (O, en el caso de la implementación de glibc, tres de los anteriores para generar 31 bits). Si eso es suficiente para sus requisitos de calidad, entonces no tiene sentido tratar de optimizar, particularmente si la probabilidad de sesgo cambia con frecuencia.

Tenga en cuenta que la duración del ciclo debe ser mucho más larga que la cantidad esperada de muestras; idealmente, debería ser mayor que el cuadrado de ese número, por lo que un generador lineal congruencial (LCG) con una longitud de ciclo de 2 ³¹ no es apropiado si espera generar gigabytes de datos aleatorios. Incluso el generador de retroalimentación aditiva no lineal trinomial de Gnu, cuya duración del ciclo se afirma que es aproximadamente 2 ³⁵ , no debe usarse en aplicaciones que requerirán millones de muestras.

Otro problema de calidad, que es mucho más difícil de probar, se relaciona con la independencia en muestras consecutivas. Las longitudes de ciclo corto fallan completamente en esta métrica, porque una vez que comienza la repetición, los números aleatorios generados se correlacionan con precisión con los valores históricos. El algoritmo trinomial de Gnu, aunque su ciclo es más largo, tiene una correlación clara como resultado del hecho de que el i ^número aleatorio generado, r _i , es siempre uno de los dos valores r _i

−3 + r _i

−31 or r _i

−3 + r _i

−31 +1. Esto puede tener consecuencias surprising o al menos puzzling , particularmente con los experimentos de Bernoulli.

Aquí hay una implementación que usa la útil biblioteca de clases de vectores de Agner Fog, que abstrae muchos de los detalles molestos en intrínsecos SSE, y también viene útil con un generador de números aleatorios vectorizados rápidos (que se encuentra en special.zip dentro del archivo vectorclass.zip ), que nos permite generar 256 bits a partir de ocho llamadas al PRNG de 256 bits. Puede leer la explicación del Dr. Fog de por qué considera que incluso el tornado Mersenne tiene problemas de calidad, y su solución propuesta; No estoy calificado para comentar, realmente, pero al menos parece dar los resultados esperados en los experimentos de Bernoulli que he intentado con él.

#include "vectorclass/vectorclass.h" #include "vectorclass/ranvec1.h" class BiasedBits { public: // Default constructor, seeded with fixed values BiasedBits() : BiasedBits(1) {} // Seed with a single seed; other possibilities exist. BiasedBits(int seed) : rng(3) { rng.init(seed); } // Generate 256 random bits, each with probability `p/256` of being 1. Vec8ui random256(unsigned p) { if (p >= 256) return Vec8ui{ 0xFFFFFFFF }; Vec32c output{ 0 }; Vec32c threshold{ 127 - p }; for (int i = 0; i < 8; ++i) { output += output; output -= Vec32c(Vec32c(rng.uniform256()) > threshold); } return Vec8ui(output); } private: Ranvec1 rng; };

En mi prueba, eso produjo y contó 268435456 bits en 260 ms, o un bit por nanosegundo. La máquina de prueba es un i5, por lo que no tiene AVX2; YMMV.

En el caso de uso real, con 201 valores posibles para p , el cálculo de los valores umbral de 8 bits será molestamente impreciso. Si esa imprecisión no es deseada, puede adaptar lo anterior para usar umbrales de 16 bits, a costa de generar el doble de números aleatorios.

Alternativamente, puede realizar una vectorización manual basada en umbrales de 10 bits, lo que le daría una muy buena aproximación a incrementos de 0.5%, utilizando el truco estándar de manipulación de bits para hacer la comparación de umbral vectorizada al verificar el préstamo en cada 10 bits de la resta del vector de valores y el umbral repetido. Combinado con, por ejemplo, std::mt19937_64 , eso le daría un promedio de seis bits cada número aleatorio de 64 bits.

Si p está cerca de 0, puede calcular la probabilidad de que el enésimo bit sea el primer bit que sea 1; luego calcula un número aleatorio entre 0 y 1 y elige n en consecuencia. Por ejemplo, si p = 0.005 (0.5%), y el número aleatorio es 0.638128, puede calcular (supongo que aquí) n = 321, por lo que debe llenar con 321 0 bits y un conjunto de bits.

Si p está cerca de 1, use 1-p en lugar de p, y establezca 1 bits más un bit 0.

Si p no está cerca de 1 o 0, haga una tabla de las 256 secuencias de 8 bits, calcule sus probabilidades acumulativas, luego obtenga un número aleatorio, haga una búsqueda binaria en la matriz de probabilidades acumulativas y puede establecer 8 bits .

Suponiendo que tiene acceso a un generador de bits aleatorios, puede generar un valor para comparar con p bit por bit y abortar tan pronto como pueda probar que el valor generado es menor o mayor o igual que p .

Proceda de la siguiente manera para crear un elemento en una secuencia con la probabilidad dada p :

Comience con 0. en binario
Añadir un bit aleatorio; suponiendo que se ha dibujado un 1 , obtendrá 0.1
Si el resultado (en notación binaria) es probablemente más pequeño que p , genera un 1
Si el resultado es probablemente mayor o igual a p , genera un 0
De lo contrario (si no se puede descartar ninguno), continúe con el paso 2.

Supongamos que p en notación binaria es 0.1001101... ; si este proceso genera cualquiera de 0.0 , 0.1000 , 0.10010 , ..., el valor no puede ser más grande o igual que p más; si alguno de 0.11 , 0.101 , 0.100111 , ... se genera, el valor no puede volverse más pequeño que p .

Para mí, parece que este método usa aproximadamente dos bits aleatorios en la expectativa. La codificación aritmética (como se muestra en la respuesta de Mark Dickinson) consume como máximo un bit aleatorio por bit sesgado (en promedio) para fijo p ; El costo de la modificación p no está claro.

Uh, los generadores de números pseudoaleatorios son generalmente bastante rápidos. No estoy seguro de qué idioma es este (Python, tal vez), pero "result.append" (que casi con certeza contiene asignación de memoria) es probablemente más lento que "random_uniform" (que solo hace un poco de matemática).

Si desea optimizar el rendimiento de este código:

Verifica que sea un problema. Las optimizaciones son un poco de trabajo y hacen que el código sea más difícil de mantener. No los hagas a menos que sea necesario.
Perfílalo. Ejecute algunas pruebas para determinar qué partes del código son realmente las más lentas. Esas son las partes que necesita para acelerar.
Realice sus cambios y verifique que en realidad sean más rápidos. Los compiladores son bastante inteligentes; A menudo, el código claro se compilará en un código mejor que algo complejo que podría aparecer más rápido.

Si está trabajando en un lenguaje compilado (incluso compilado JIT), recibe un impacto de rendimiento por cada transferencia de control (si, mientras, llamada de función, etc.). Elimina lo que puedas. La asignación de memoria también es (generalmente) bastante costosa.

Si está trabajando en un idioma interpretado, todas las apuestas están canceladas. Es muy probable que el código más simple sea el mejor. La sobrecarga del intérprete empequeñecerá todo lo que esté haciendo, por lo tanto, reduzca su trabajo tanto como sea posible.

Solo puedo adivinar dónde están sus problemas de rendimiento:

Asignación de memoria. Preasigne la matriz en su tamaño completo y complete las entradas más adelante. Esto asegura que no será necesario reasignar la memoria mientras agrega las entradas.
Ramas Es posible que pueda evitar el "si" lanzando el resultado o algo similar. Esto dependerá mucho del compilador. Verifique el ensamblaje (o perfil) para verificar que haga lo que desea.
Tipos numéricos Descubra el tipo que su generador de números aleatorios usa de forma nativa y haga su aritmética en ese tipo. Por ejemplo, si el generador devuelve números enteros sin signo de 32 bits, escale primero "p" a ese rango y luego úselo para la comparación.

Por cierto, si realmente desea utilizar la menor cantidad de aleatoriedad posible, utilice la "codificación aritmética" para decodificar su secuencia aleatoria. No será rápido

Una cosa que puede hacer es muestrear varias veces desde el generador imparcial subyacente, obteniendo varias palabras de 32 bits o 64 bits, y luego realizando aritmética booleana bit a bit. Como ejemplo, para 4 palabras b1,b2,b3,b4 , puede obtener las siguientes distribuciones:

expression | p(bit is 1) -----------------------+------------- b1 & b2 & b3 & b4 | 6.25% b1 & b2 & b3 | 12.50% b1 & b2 & (b3 | b4) | 18.75% b1 & b2 | 25.00% b1 | (b2 & (b3 | b4)) | 31.25% b1 & (b2 | b3) | 37.50% b1 & (b2 | b3 | b4)) | 43.75% b1 | 50.00%

Se pueden hacer construcciones similares para resoluciones más finas. Se vuelve un poco tedioso y aún requiere más llamadas de generador, pero al menos no una por bit. Esto es similar a la respuesta de a3f, pero es probable que sea más fácil de implementar y, sospecho, más rápido que escanear palabras para 0xF nybbles 0xF.

Tenga en cuenta que para su resolución deseada del 0.5%, necesitaría 8 palabras imparciales para una palabra sesgada, lo que le daría una resolución de (0.5 ^ 8) = 0.390625%.

Una forma que daría un resultado preciso es generar primero aleatoriamente para un bloque de k bits el número de 1 bits que sigue a la distribución binomial, y luego generar una palabra de k bits con exactamente esa cantidad de bits utilizando uno de los métodos here . Por ejemplo, el método por mic006 solo necesita sobre números aleatorios log k k-bit, y el mío solo necesita uno.