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

Question

Estoy escribiendo algunas pruebas para un C ++ línea de comandos de Linux aplicación. Me gustaría generar un montón de números enteros con una distribución de cola larga de ley de potencia /. Es decir, consigo un algunos números muy frecuentemente, pero la mayoría de ellos con poca frecuencia.

Lo ideal sería sólo algunas ecuaciones mágicas que podría utilizar con rand () o una de las funciones aleatorias stdlib. Si no es así, un fácil utilizar el trozo de C / C ++ sería grande.

Gracias!

Answer 1

Esta página rel="noreferrer"> href="http://mathworld.wolfram.com/RandomNumber.html" discute cómo conseguir una distribución de ley potencial de una distribución uniforme (que es lo que la mayoría de los generadores de números aleatorios proporcionan).

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

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

donde y es una variable aleatoria uniforme, n es el poder de distribución, x0 y x1 definir el rango de la distribución, y x variate se distribuye el poder de la ley.

Answer 2

Si conoce la distribución que desee (llamada de la función de distribución de probabilidad (PDF)) y tiene correctamente normalizada, se puede integrar para obtener la función de distribución acumulativa (CDF), y luego invertir la CDF (si es posible) para obtener la transformación que necesita de la distribución uniforme [0,1] a su deseada.

Así se inicia mediante la definición de la distribución que desee.

P = F(x)

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

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

Si esto se puede invertir a obtener

y = F^{-1}(C)

Así que llame rand() y enchufe el resultado en C como en la última línea y el uso de y.

Este resultado se llama el teorema fundamental del muestreo. Esta es una molestia debido al requisito de la normalización y la necesidad de invertir la función analítica.

Como alternativa se puede utilizar una técnica de rechazo: lanzar una serie de manera uniforme en el rango deseado, a continuación, lanzar otro número y comparar con el PDF en la ubicación indeicated por su primer tiro. Rechazar si la segunda tirada supera el PDF. Tiende a ser ineficiente para archivos PDF con un montón de baja probabilidad región, como los que tienen colas largas ...

Un enfoque intermedio implica la inversión de la CDF por la fuerza bruta: almacenar la CDF como una tabla de búsqueda, y hacer una búsqueda inversa para obtener el resultado

.

---

El stinker real aquí es que las distribuciones x^-n simples no son normalizable en el [0,1] rango, por lo que no puede usar el teorema de muestreo. Pruebe (x + 1) ^ - N en vez ...

Answer 3

No puedo comentar sobre la matemática necesaria para producir una distribución de ley de potencia (los otros mensajes tienen sugerencias) pero sugeriría que se familiarice con las instalaciones de números aleatorios TR1 biblioteca de C ++ estándar 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 un montón de ajustes preestablecidos.

Los preajustes de distribución incluidos son:

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

Cuando se define la distribución de ley de potencia, debe ser capaz de conectarlo con los generadores y motores existentes. El libro El C ++ Standard Library extensiones por Pete Becker tiene un gran capítulo sobre <random>.

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

Answer 4

Sólo quería llevar a cabo una simulación real como complemento a la respuesta (con razón) aceptado. Aunque en R, el código es tan simple como sea (pseudo) -pseudo-código.

Una pequeña diferencia entre el Wolfram MathWorld fórmula en la respuesta aceptada y otra, quizás más común, ecuaciones es el hecho de que el ley de potencias exponente n (que por lo general se denota como alfa) no lleva un signo negativo explícito. Así que el valor alfa elegido tiene que ser negativa, y por lo general entre 2 y 3.

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

Así que aquí está:

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