由于我一直在进行社交网络分析，我遇到了拟合网络度数概率分布的问题。
因此，我有一个概率分布 P(X >= x)，从可视化检查中可以得知它遵循一个具有指数截尾的幂律而不是纯粹的幂律（一条直线）。
那么，给定具有指数截尾的幂律分布方程：
f(x) = x**alpha * exp(beta*x)
如何使用Python估计参数 alpha 和 beta？
我知道scipy.stats.powerlaw包存在，并且他们有一个.fit()函数，但看起来它并不能胜任这项工作，因为它只返回图形的位置和比例尺，这似乎仅对正态分布有用？此外，关于该包的教程也不足够。
P.S. 我很清楚 CLauset等人 的实现，但他们似乎没有提供估计其他分布参数的方法。
Powerlaw库可以直接用于以下参数的估计:
Install all the pythons dependencies:
pip install powerlaw mpmath scipy
Run the powerlaw package fit in a python environment:
import powerlaw
data = [5, 4, ... ]
results = powerlaw.Fit(data)
get the parameters from the results
results.truncated_power_law.parameter1 # power law parameter (alpha)
results.truncated_power_law.parameter2 # exponential cut-off parameter (beta)
函数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
这里提供一种在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)
}