python scipy.stats.powerlaw exponente negativo

statistics idl-programming-language

Question

Quiero proporcionar un exponente negativo para la rutina scipy.stats.powerlaw, por ejemplo, a = -1.5, para dibujar muestras aleatorias:

"""
powerlaw.pdf(x, a) = a * x**(a-1)
"""

from scipy.stats import powerlaw
R = powerlaw.rvs(a, size=100)

¿Por qué se requiere un> 0, cómo puedo suministrar un negativo a para generar las muestras aleatorias, y cómo puedo suministrar un coeficiente / transformación de normalización, es decir,

PDF(x,C,a) = C * x**a

La documentación está aquí

http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.powerlaw.html

¡Gracias!

EDIT: debería agregar que estoy tratando de replicar la función RANDOMP de IDL:

http://idlastro.gsfc.nasa.gov/ftp/pro/math/randomp.pro

Answer 1

Si desea generar una distribución de ley de potencia, puede usar una desviación aleatoria. Solo debe generar un número aleatorio entre [0,1] y aplicar el método inverso ( Wolfran ). En este caso, la función de densidad de probabilidad es:

p (k) = k ^ (- gamma)

y y es la variable uniforme entre 0 y 1.

y ~ U (0,1)

import numpy as np

def power_law(k_min, k_max, y, gamma):
    return ((k_max**(-gamma+1) - k_min**(-gamma+1))*y  + k_min**(-gamma+1.0))**1.0/(-gamma + 1.0)

Ahora para generar una distribución, solo tienes que crear una matriz

nodes = 1000
scale_free_distribution = np.zeros(nodes, float)
k_min = 1.0
k_max = 100*k_min
gamma = 3.0

for n in range(nodes):
    scale_free_distribution[n] = power_law(k_min, k_max,np.random.uniform(0,1), gamma)

Esto funcionará para generar una distribución de ley de potencia con gamma = 3.0, si desea arreglar el promedio de distribución, debe estudiar las redes complejas porque el k_min depende de k_max y la conectividad promedio.

Answer 2

Si r es una desviación aleatoria uniforme de U (0,1), entonces x en la siguiente expresión es una desviación aleatoria distribuida de la ley de potencias:

x = xmin * (1-r) ** (-1/(alpha-1))

donde xmin es el valor más pequeño (positivo) por encima del cual se cumple la distribución de la ley de poder, y alfa es el exponente de la distribución.

Answer 3

Un PDF, integrado en su dominio, debe ser igual a uno. En otras palabras, el área bajo una curva de función de densidad de probabilidad debe ser igual a uno.

In [36]: import scipy.integrate as integrate
In [40]: y, err = integrate.quad(lambda x: 0.5*x**(-0.5), 0, 1)

In [41]: y
Out[41]: 0.9999999999999998  # The integral is close to 1

La función de densidad de powerlaw tiene un dominio de 0 <= x <= 1. En este dominio, la integral de x**b es finita para cualquier b > -1. Cuando b es más pequeño, x**b explota demasiado rápido cerca de x = 0 . Entonces, no es una función de densidad de probabilidad válida cuando b <= -1 .

In [38]: integrate.quad(lambda x: x**(-1), 0, 1)
UserWarning: The maximum number of subdivisions (50) has been achieved...
# The integral blows up

Por lo tanto, para x**(a-1) , a debe satisfacer a-1 > -1 o equivalentemente, a > 0 .

La primera constante a en a * x**(a-1) es la constante de normalización que hace que la integral de a * x**(a-1) sobre el dominio [0,1] sea igual a 1. Así que no lo hace t puede elegir esta constante independientemente de a .

Ahora bien, si cambia el dominio para que esté a una distancia medible de 0, entonces sí, podría definir un PDF de la forma C * x**a para a negativo a . Pero tendría que indicar qué dominio desea, y no creo que haya (todavía) un PDF disponible en scipy.stats para esto.