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

Question

由于我一直在做一些社会网络分析，我偶然发现了在网络程度上拟合概率分布的问题。

所以，我有一个概率分布的P(X >= x)，从视觉上看，它遵循指数截止的幂律，而不是纯幂律(直线)。

因此，假设幂律分布的指数截止方程是：

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

如何使用Python估计参数alpha和beta？

我知道scipy.stats.powerlaw包是存在的，他们有一个.fit()函数，但这似乎不起作用，因为它只返回绘图的位置和规模，这似乎只对正态分布有用吗？在这个软件包上也没有足够的教程。

我很清楚CLauset等人的实现，但它们似乎没有提供估计备用发行版参数的方法。

Answer 1

强力法库可以直接用于估计参数如下：

1. 安装所有pythons依赖项： pip安装powerlaw mpmath
2. 在python环境中运行powerlaw包： 进口动力数据= 5，4，.结果=powerlaw.Fit(数据)
3. 从结果中获取参数 results.truncated_power_law.parameter1 #幂律参数(α) results.truncated_power_law.parameter2 #指数截断参数(β)

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函数，然后使用and的通用拟合算法来学习参数。例如：http://www.scipy.org/Cookbook/FittingData#head-5eba0779a34c07f5a596bbcf99dbc7886eac18e5

Answer 3

本文通过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/以获得更多信息