Generando con precisión variantes de distribución discreta de la ley de potencia

¿Cuáles son los mejores métodos para generar con precisión enteros aleatorios distribuidos de acuerdo con una ley de potencia? La probabilidad de obtener ( ) debería ser igual a y el método debería funcionar bien para cualquier .$k$$k=1,2,\ldots$$p_k = k^{-\gamma} / \zeta(\gamma)$$\gamma > 1$

Puedo ver dos enfoques ingenuos:

1. Calcule hasta un gran para que esté "lo suficientemente cerca" de 1, luego genere números enteros de acuerdo con estas probabilidades. Esto simplemente no funcionará si está cerca de 1 ya que necesitaría ser enorme.$p_k$$k_\text{max}$$\sum_{k=1}^{k_\text{max}}$$\gamma$$k_\text{max}$

2. Dibuje números reales de una distribución continua de la ley de potencia (un problema más fácil que sé cómo resolver) y redondee a los enteros de alguna manera. Es posible calcular analíticamente la probabilidad precisa de obtener cada número entero con el método anterior. Podría usar el rechazo para corregirlos en (que también se puede calcular siempre que pueda evaluar la ). (Esto sería un poco complicado ya que tendría que redondear de manera que obtenga enteros con mayor probabilidad que para mayor que algún valor pequeño, y manejar menos que eso por separado).$p_k$$\zeta$$p_k$$k$$k$

¿Existe un método mejor que también sea preciso (no aproximado)?

Creo que (una versión ligeramente modificada de) el método 2 es bastante sencillo, en realidad

Usando la definición de la función de distribución de Pareto dada en Wikipedia

si toma y entonces la relación de a se maximiza en , lo que significa que puede escalar según la relación en y usar un muestreo de rechazo directo. Parece ser razonablemente eficiente.$x_m=\frac{1}{2}$$\alpha=\gamma$$p_x$$q_x=F_X(x+\frac{1}{2})-F_X(x-\frac{1}{2})$$x=1$$x=1$

Para ser más explícito: si genera a partir de un Pareto con y y redondea al entero más cercano (en lugar de truncar), entonces parece posible utilizar el muestreo de rechazo con : cada valor generado de de ese proceso se acepta con probabilidad .$x_m=\frac{1}{2}$$\alpha=\gamma$$M = p_1/q_1$$x$$\frac{p_x}{Mq_x}$

( aquí fue ligeramente redondeado ya que soy flojo; en realidad, el ajuste para este caso sería un poco diferente, pero no lo suficiente como para verse diferente en la trama; de hecho, la imagen pequeña hace que parezca un poco demasiado pequeño cuando en realidad es una fracción demasiado grande)$M$

Un ajuste más cuidadoso de y ( para entre 0 y 1) probablemente aumentaría aún más la eficiencia, pero este enfoque funciona razonablemente bien en los casos con los que he jugado.$x_m$$\alpha$$\alpha=\gamma-a$$a$

Si puede dar una idea del rango típico de valores de , puedo echar un vistazo más de cerca a la eficiencia allí.$\gamma$

El método 1 se puede adaptar para ser exacto, también, realizando el método 1 casi siempre, y luego aplicando otro método para lidiar con la cola. Esto se puede hacer de maneras muy rápidas.

Por ejemplo, si toma un vector entero de longitud 256 y llena los primeros valores de con , los siguientes valores de con y así sucesivamente hasta , eso será casi usa toda la matriz. Las pocas celdas restantes indican luego pasar a un segundo método que combina el manejo de la cola derecha y también los pequeños bits de probabilidad 'sobrantes' de la parte izquierda.$\lfloor 256 p_1\rfloor$1$\lfloor 256 p_2\rfloor$2$256 p_i <1$

El remanente izquierdo se puede hacer mediante una serie de enfoques (incluso con, digamos 'cuadrar el histograma' si está automatizado, pero no tiene que ser tan eficiente como eso), y la cola derecha se puede hacer usando algo como el enfoque de aceptar-rechazar anterior.

El algoritmo básico consiste en generar un número entero de 1 a 256 (que requiere solo 8 bits del rng; si la eficiencia es primordial, las operaciones de bits pueden sacarlos de la parte superior, dejando el resto del número uniforme (lo mejor sería dejado como un valor entero no normalizado hasta este punto) que puede usarse para tratar el remanente izquierdo y la cola derecha si es necesario.

Cuidadosamente implementado, este tipo de cosas puede ser muy rápido. Puede usar diferentes valores de que 256 (por ejemplo, podría ser una posibilidad), pero no todo es igual. Sin embargo, si toma una tabla muy grande, puede que no queden suficientes bits en el uniforme para que sea adecuada para generar la cola y necesita un segundo valor uniforme allí (pero rara vez se necesita, por lo que no es demasiado un problema)$2^k$$2^{16}$

En el mismo ejemplo de zeta (2) que el anterior, tendría 212 1's, 26 2' s, 7 3's, 3 4' s, uno 5y los valores de 250-256 tratarían con el remanente. Más del 97% del tiempo genera uno de los valores de la tabla (1-5).

Hasta donde yo sé, el estado del arte en leyes de poder es el documento de Clauset, Shalizi y Newman que discute su problema en el Apéndice D. Note en particular (donde es un dibujo de una ley de poder continua) dicen:$y$

Otros enfoques aproximados para generar enteros, como redondear hacia abajo (truncar) el valor de y, dan resultados sustancialmente más pobres y no deben usarse.

Como alternativa a la respuesta aceptada, Clauset et al. El método para obtener dibujos precisos de la distribución de la ley de potencia discreta es dibujar un aleatorio uniforme y luego hacer donde es el cdf complementario de la ley de potencia discreta. Necesita la función zeta para calcular pero solo debe calcularse con cierta precisión, por lo que es posible generar sorteos que tengan la distribución discreta de la ley de potencia de esta manera. Necesitas usar el método de bisección para resolver la ecuación .$r \in [0,1)$$x= P^{-1}(1-r)$$P(x) = \sum_{a=x}^\infty P(X=a)$$P(x)$$P(x) = 1-r$

Debido a que el cálculo exacto es costoso, también se proporciona un método aproximado, que es definir que no es lo mismo que solo redondear valores de la ley de potencia continua. El error de esta aproximación se da en la ecuación (D.7) de Clauset et al. y depende de .