Python scipy.stats.powerlaw exponente negativo

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 una a negativa para generar las muestras aleatorias y cómo puedo proporcionar un coeficiente / transformada 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!

EDITAR: Debo agregar que estoy tratando de replicar la función RANDOMP de IDL:

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

Un PDF, integrado sobre su dominio, debe ser igual a uno. En otras palabras, el área bajo la curva de la 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 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 . Por lo tanto, 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 ' t llegar a elegir esta constante independiente de a .

Ahora, 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 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.

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

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

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

Si desea generar una distribución de ley de potencia, puede usar una desviación aleatoria. Solo tienes que 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)

e 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 corregir el promedio de distribución, tiene que estudiar Redes complejas porque el k_min depende de k_max y la conectividad promedio.

El paquete powerlaw de Python puede hacer esto. Considere para a>1 una distribución de ley de potencia con función de densidad de probabilidad

 f(x) = c * x^(-a) 

para x > x_min y f(x) = 0 contrario. Aquí c es un factor de normalización y se determina como

 c = (a-1) * x_min^(a-1). 

En el siguiente ejemplo, es a = 1.5 y x_min = 1.0 y al comparar la función de densidad de probabilidad estimada de la muestra aleatoria con el PDF de la expresión anterior se obtiene el resultado esperado.

 import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as pl import numpy as np import powerlaw a, xmin = 1.5, 1.0 N = 10000 # generates random variates of power law distribution vrs = powerlaw.Power_Law(xmin=xmin, parameters=[a]).generate_random(N) # plotting the PDF estimated from variates bin_min, bin_max = np.min(vrs), np.max(vrs) bins = 10**(np.linspace(np.log10(bin_min), np.log10(bin_max), 100)) counts, edges = np.histogram(vrs, bins, density=True) centers = (edges[1:] + edges[:-1])/2. # plotting the expected PDF xs = np.linspace(bin_min, bin_max, 100000) pl.plot(xs, [(a-1)*xmin**(a-1)*x**(-a) for x in xs], color='red') pl.plot(centers, counts, '.') pl.xscale('log') pl.yscale('log') pl.savefig('powerlaw_variates.png') 

devoluciones

ley de potencia