Estimación de corte exponencial en una distribución de ley de potencia.

Como he estado haciendo algunos análisis de redes sociales, me he topado con el problema de ajustar una distribución de probabilidad en el grado de la red.

Entonces, tengo una distribución de probabilidad P(X >= x) que, de la inspección visual, sigue una ley de potencia con un corte exponencial en lugar de una ley de potencia pura (una línea recta).

Entonces, dado que la ecuación para la distribución de la ley de potencia con corte exponencial es:

f (x) = x ** alfa * exp (beta * x)

¿Cómo puedo estimar los parámetros alpha y beta usando Python?

Sé que existe el paquete scipy.stats.powerlaw y tienen una función .fit() pero eso no parece funcionar, ya que solo devuelve la ubicación y la escala del gráfico, lo que parece ser útil solo para la distribución normal. Tampoco hay suficientes tutoriales en este paquete.

PD: Soy muy consciente de la implementación de CLauset et al, pero parece que no proporcionan formas de estimar los parámetros de distribuciones alternativas.

La función scipy.stats.powerlaw.fit aún puede funcionar para sus propósitos. Es un poco confuso cómo funcionan las distribuciones en scipy.stats (la documentación de cada uno se refiere a la ubicación y escala de los parámetros opcionales, aunque no todos ellos usan estos parámetros, y cada uno los usa de manera diferente). Si nos fijamos en los documentos:

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

También hay un segundo parámetro no opcional “a”, que es el de “parámetros de forma”. En el caso de powerlaw, este contiene un solo parámetro. No te preocupes por “loc” y “escala”.

Edición: Lo sentimos, olvidé que también querías el parámetro beta. Lo mejor que puedes hacer es definir la función de powerlaw que deseas, y luego usar los algoritmos de adaptación generics de scipy para aprender los parámetros. Por ejemplo: http://www.scipy.org/Cookbook/FittingData#head-5eba0779a34c07f5a596bbcf99dbc7886eac18e5

Aquí hay un medio para estimar el exponente de escala y la tasa exponencial de la ley de potencia con corte exponencial maximizando la probabilidad en R:

 # Input: Data vector, lower threshold # Output: List, giving type ("powerexp"), scaling exponent, exponential rate, lower threshold, log-likelihood powerexp.fit <- function(data,threshold=1,method="constrOptim",initial_rate=-1) { x <- data[data>=threshold] negloglike <- function(theta) { -powerexp.loglike(x,threshold,exponent=theta[1],rate=theta[2]) } # Fit a pure power-law distribution pure_powerlaw <- pareto.fit(data,threshold) # Use this as a first guess at the exponent initial_exponent <- pure_powerlaw$exponent if (initial_rate < 0) { initial_rate <- exp.fit(data,threshold)$rate } minute_rate <- 1e-6 theta_0 <- as.vector(c(initial_exponent,initial_rate)) theta_1 <- as.vector(c(initial_exponent,minute_rate)) switch(method, constrOptim = { # Impose the constraint that rate >= 0 # and that exponent >= -1 ui <- rbind(c(1,0),c(0,1)) ci <- c(-1,0) # Can't start with values on the boundary of the feasible set so add # tiny amounts just in case if (theta_0[1] == -1) {theta_0[1] <- theta_0[1] + minute_rate} if (theta_0[2] == 0) {theta_0[2] <- theta_0[2] + minute_rate} est <- constrOptim(theta=theta_0,f=negloglike,grad=NULL,ui=ui,ci=ci) alpha <- est$par[1] lambda <- est$par[2] loglike <- -est$value}, optim = { est <- optim(par=theta_0,fn=negloglike) alpha <- est$par[1] lambda <- est$par[2] loglike <- -est$value}, nlm = { est.0 <- nlm(f=negloglike,p=theta_0) est.1 <- nlm(f=negloglike,p=theta_1) est <- est.0 if (-est.1$minimum > -est.0$minimum) { est <- est.1;cat("NLM had to switch\n") } alpha <- est$estimate[1] lambda <- est$estimate[2] loglike <- -est$minimum}, {cat("Unknown method",method,"\n"); alpha<-NA; lambda<-NA; loglike<-NA} ) fit <- list(type="powerexp", exponent=alpha, rate=lambda, xmin=threshold, loglike=loglike, samples.over.threshold=length(x)) return(fit) } 

Echa un vistazo a https://github.com/jeffalstott/powerlaw/ para más información

La biblioteca Powerlaw se puede usar directamente para estimar los parámetros de la siguiente manera:

  1. Instale todas las dependencias de los pitones:

     pip install powerlaw mpmath scipy 
  2. Ejecute el paquete powerlaw en un entorno python:

     import powerlaw data = [5, 4, ... ] results = powerlaw.Fit(data) 
  3. Obtener los parámetros de los resultados.

     results.truncated_power_law.parameter1 # power law parameter (alpha) results.truncated_power_law.parameter2 # exponential cut-off parameter (beta) 

También estoy trabajando en el campo de las redes y tuve que hacer un ajuste muy similar al tuyo. Encontré una solución muy fácil y rápida aquí , y la mejor parte es que no necesita instalar ningún paquete aparte de Scipy (que estoy seguro de que ya tiene).

La distribución que quería ajustar es una ley de corte de energía con un cambio, como la que se describe en este documento . Usando tu misma notación, mi ajuste es

 f(x) = (x + x0)**alpha * exp(-beta*x) 

así que solo agrega el tercer parámetro x0 a tu distribución. Tenga en cuenta que asumo que la beta es positiva y que simplemente tomo el signo de afuera (creo que esto aclara que su exponencial está disminuyendo).

La implementación es así:

 import numpy as np. import scipy.optimize as opt def distribution(x, alpha, beta, x0): return (x + x0)**alpha * np.exp(-beta *x) # ... I prepare my data here fit = opt.curve_fit(distribution, x_data, y_data) # you can pass guess for the parameters/errors alpha, beta, x0 = fit[0] 

Este es el resultado:

ajuste