估计幂律分布中的指数截断

Question

我在做社交网络分析的时候，遇到了一个问题，就是如何给网络的度数拟合一个概率分布。

我有一个概率分布 P(X >= x)，通过观察发现，它遵循一种带有指数截断的幂律分布，而不是纯粹的幂律分布（也就是一条直线）。

那么，带有指数截断的幂律分布的公式是：

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

我想知道如何用Python来估计参数 alpha 和 beta。

我知道有一个叫做scipy.stats.powerlaw的包，并且里面有一个 .fit() 函数，但这个函数似乎不太管用，因为它只返回图形的位置和尺度，这对正态分布可能有用，但对我这个情况不太适合？而且关于这个包的教程也不多。

另外，我也知道 CLauset等人的实现，但他们似乎没有提供估计其他分布参数的方法。

Answer 1

这里有一种方法，可以通过在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)
}

想了解更多信息，可以查看 https://github.com/jeffalstott/powerlaw/

Answer 2

函数 scipy.stats.powerlaw.fit 可能还是能满足你的需求。关于 scipy.stats 中的分布，确实有点让人困惑（每个分布的文档都提到可选参数 loc 和 scale，尽管并不是所有的分布都使用这些参数，而且每个分布对它们的使用方式也不一样）。如果你查看文档：

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

你会发现还有一个非可选参数 "a"，它是“形状参数”。在 powerlaw 的情况下，这个参数只有一个。你不用担心 "loc" 和 "scale"。

补充：抱歉，我忘了你还想要 beta 参数。你最好的选择可能是自己定义你想要的 powerlaw 函数，然后使用 scipy 的通用拟合算法来学习参数。例如：

http://www.scipy.org/Cookbook/FittingData#head-5eba0779a34c07f5a596bbcf99dbc7886eac18e5

Answer 3

Powerlaw库可以直接用来估算参数，步骤如下：

1. 安装所有需要的Python依赖：

   pip install powerlaw mpmath scipy
   

2. 在Python环境中运行powerlaw包的拟合：

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

3. 从结果中获取参数：

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