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

Question

Estoy escribiendo algunas pruebas para una aplicación de línea de comandos C++ de Linux. 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.

Lo ideal sería que hubiera 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

Este page at Wolfram MathWorld explica cómo obtener una distribución de ley de potencia 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 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 es la variable distribuida de la ley de poder.

Answer 2

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) a obtenga la transformación que necesita de la distribución uniforme [0,1] a su gusto.

Así que empiezas definiendo la distribución que deseas.

P = F(x) 

(para x en [0,1]) después se integra 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 como C en la última línea y usa y.

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

Alternativamente puede utilizar 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 ineficaz 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 un reverso búsqueda para obtener el resultado.

---

El verdadero stinker aquí es que simples x^-n distribuciones no son normalizable en la gama [0,1], por lo que no se puede usar el teorema de muestreo. Pruebe (x + 1)^- n en su lugar ...

Answer 3

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

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

Al definir su distribución de ley de potencia, debe ser capaz de conectarlo con los generadores y motores existentes. El libro Las extensiones de biblioteca estándar de C++ por Pete Becker tiene un gran capítulo en <random>.

Here is an article acerca de cómo crear otras distribuciones (con ejemplos de Cauchy, Chi-cuadrado, 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 ser (pseudo) -seudocódigo.

Una pequeña diferencia entre el Wolfram MathWorld formula en la respuesta aceptada y otra, quizás más común, ecuaciones es el hecho de que la ley de potencia exponenten (que por lo general se denota como alfa) no lleva un signo negativo explícito. Por lo tanto, 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.

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 representan 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