c++ - sorteos - tabla numeros aleatorios

¿Generador de números aleatorios que produce una distribución de ley de potencia?

Question

Estoy escribiendo algunas pruebas para una aplicación Linux de línea de comandos C ++. Me gustaría generar un montón de enteros con una ley de poder / distribución de cola larga. Es decir, obtengo algunos números con mucha frecuencia, pero la mayoría son relativamente poco frecuentes.

Idealmente, solo habría algunas ecuaciones mágicas que podría usar con rand () o una de las funciones aleatorias stdlib. Si no, un pedazo fácil de usar de C / C ++ sería genial.

¡Gracias!

Answer 1

Esta página en Wolfram MathWorld discute cómo obtener una distribución de ley de poder a partir de una distribución uniforme (que es lo que proporcionan la mayoría de los generadores de números aleatorios).

La respuesta corta (derivación en el enlace de arriba):

x = [(x1^(n+1) - x0^(n+1))*y + x0^(n+1)]^(1/(n+1))

donde y es una variable uniforme, n es la potencia de distribución, x0 y x1 definen el rango de la distribución, y x es la variable distribuida de la ley de potencia.

Answer 2

No puedo comentar las matemáticas necesarias para producir una distribución de la ley de poder (las otras publicaciones tienen sugerencias) pero sugiero que se familiarice con las instalaciones de números aleatorios de la Biblioteca de estándares TR1 C ++ en <random> . Estos proporcionan más funcionalidad que std::rand y std::srand . El nuevo sistema especifica una API modular para generadores, motores y distribuciones y suministra muchos presets.

Los ajustes preestablecidos de distribución incluidos son:

• uniform_int
• bernoulli_distribution
• geometric_distribution
• poisson_distribution
• binomial_distribution
• uniform_real
• exponential_distribution
• normal_distribution
• gamma_distribution

Cuando defina la distribución de la ley de potencia, debería poder conectarla con los generadores y motores existentes. El libro The C ++ Standard Library Extensions de Pete Becker tiene un gran capítulo sobre <random> .

Aquí hay un artículo sobre cómo crear otras distribuciones (con ejemplos para Cauchy, Chi-squared, Student t y Snedecor F)

Answer 3

Si conoce la distribución que desea (llamada Función de distribución de probabilidad (PDF)) y la ha normalizado correctamente, puede integrarla para obtener la Función de distribución acumulativa (CDF), luego invierta la CDF (si es posible) para obtener la transformación que necesita una distribución uniforme [0,1] a la deseada.

Entonces comienzas por definir la distribución que deseas.

P = F(x)

(para x en [0,1]) luego integrado para dar

C(y) = /int_0^y F(x) dx

Si esto puede invertirse obtienes

y = F^{-1}(C)

Entonces llame a rand() y conecte el resultado como C en la última línea y use y.

Este resultado se llama Teorema Fundamental de Muestreo. Esto es una molestia debido a los requisitos de normalización y la necesidad de invertir la función analíticamente.

Alternativamente puede usar una técnica de rechazo: arroje un número uniformemente en el rango deseado, luego arroje otro número y compare con el PDF en el lugar indeseado por su primer lanzamiento. Rechazar si el segundo lanzamiento excede el PDF. Tiende a ser ineficiente para PDF con mucha región de baja probabilidad, como aquellos con colas largas ...

Un enfoque intermedio implica invertir el CDF por fuerza bruta: almacena el CDF como una tabla de búsqueda y realiza una búsqueda inversa para obtener el resultado.

El verdadero inconveniente aquí es que las distribuciones simples x^-n no son normalizables en el rango [0,1] , por lo que no puede usar el teorema de muestreo. Pruebe (x + 1) ^ - n en su lugar ...

Answer 4

Solo quería llevar a cabo una simulación real como complemento a la respuesta (legítimamente) aceptada. Aunque en R, el código es tan simple como ser (pseudo) -seudocódigo.

Una pequeña diferencia entre la fórmula de Wolfram MathWorld en la respuesta aceptada y otras ecuaciones, quizás más comunes, es el hecho de que el exponente de la ley de potencia n (que típicamente se denota como alfa) no lleva un signo negativo explícito. Entonces, el valor alfa elegido debe ser negativo, y típicamente entre 2 y 3.

x0 y x1 representan los límites inferior y superior de la distribución.

Asi que aqui esta:

x1 = 5 # Maximum value
x0 = 0.1 # It can''t be zero; otherwise X^0^(neg) is 1/0.
alpha = -2.5 # It has to be negative.
y = runif(1e5) # Number of samples
x = ((x1^(alpha+1) - x0^(alpha+1))*y + x0^(alpha+1))^(1/(alpha+1))
hist(x, prob = T, breaks=40, ylim=c(0,10), xlim=c(0,1.2), border=F,
col="yellowgreen", main="Power law density")
lines(density(x), col="chocolate", lwd=1)
lines(density(x, adjust=2), lty="dotted", col="darkblue", lwd=2)

o trazado en escala logarítmica:

h = hist(x, prob=T, breaks=40, plot=F)
plot(h$count, log="xy", type=''l'', lwd=1, lend=2,
xlab="", ylab="", main="Density in logarithmic scale")

Aquí está el resumen de los datos:

> summary(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.1000 0.1208 0.1584 0.2590 0.2511 4.9388