Estimación del límite exponencial en una distribución de ley de potencia

Question

Como he estado haciendo algunos análisis de redes sociales, he tropezado 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, por inspección visual, sigue una ley de potencia con un límite exponencial en lugar de una ley pura de potencia (una línea recta).

Por lo tanto, dado que la ecuación de la ley de distribución de potencia con corte exponencial es:

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

¿Cómo puede ocurrir 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 hacer el trabajo, ya que solo devuelve la ubicación y la escala de la gráfica, que parece ser útil solo para la distribución normal ? Tampoco hay suficientes tutoriales en este paquete.

P.S. Soy muy consciente de la implementación de CLauset et al, pero no parecen proporcionar formas de estimar los parámetros de las distribuciones alternativas.

Answer 1

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 los parámetros opcionales loc y scale, aunque no todos 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 "a" no opcional, que son los "parámetros de forma". En el caso de Powerlaw, este contiene un solo parámetro. No te preocupes por "loc" y "scale".

Editar: Lo sentimos, se olvidó de que también quería el parámetro beta. Lo mejor que puede hacer es definir la función powerlaw que usted quiere, y luego usar los algoritmos genéricos de adaptación de scipy para aprender los parámetros.Por ejemplo: http://www.scipy.org/Cookbook/FittingData#head-5eba0779a34c07f5a596bbcf99dbc7886eac18e5

Answer 2

Aquí es un medio para estimar el exponente de escala y la tasa exponencial de la ley de potencia con exponencial de corte mediante la maximización de 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) 
} 

Salida https://github.com/jeffalstott/powerlaw/ para obtener más información

Answer 3

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

1. instalar todos los pitones dependencias:

   pip install powerlaw mpmath scipy 
   

2. ejecutar el paquete powerlaw ajuste en un entorno Python:

   import powerlaw 
   data = [5, 4, ... ] 
   results = powerlaw.Fit(data) 
   

3. obtener los parámetros a partir de los resultados

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