冪律（）`的`包功能fitdistrplus`

Question

我生成使用包rplcon()功能poweRlaw

data <- rplcon(1000,10,2)

現在，我想知道哪些已知分佈擬合數據的一些隨機變量最好。 Lognorm？ EXP？伽瑪？冪律？指數截斷的冪律？

於是我就用功能fitdist()封裝fitdistrplus：

fit.lnormdl <- fitdist(data,"lnorm") 
fit.gammadl <- fitdist(data, "gamma", lower = c(0, 0)) 
fit.expdl <- fitdist(data,"exp") 

由於冪律分佈和冪律指數截止並非根據CRAN Task View: Probability Distributions基礎概率函數，所以我寫的d，P，基於對?fitdist

dplcon <- function (x, xmin, alpha, log = FALSE) 
{ 
    if (log) { 
     pdf = log(alpha - 1) - log(xmin) - alpha * (log(x/xmin)) 
     pdf[x < xmin] = -Inf 
    } 
    else { 
     pdf = (alpha - 1)/xmin * (x/xmin)^(-alpha) 
     pdf[x < xmin] = 0 
    } 
    pdf 
} 
pplcon <- function (q, xmin, alpha, lower.tail = TRUE) 
{ 
    cdf = 1 - (q/xmin)^(-alpha + 1) 
    if (!lower.tail) 
     cdf = 1 - cdf 
    cdf[q < round(xmin)] = 0 
    cdf 
} 
qplcon <- function(p,xmin,alpha) alpha*p^(1/(1-xmin)) 

最後的例子4功法的q函數，我用下面的代碼來獲取參數xmin和p的alpha奧爾法：

fitpl <- fitdist(data,"plcon",start = list(xmin=1,alpha=1)) 

但它拋出一個錯誤：

<simpleError in optim(par = vstart, fn = fnobj, fix.arg = fix.arg, obs = data,  ddistnam = ddistname, hessian = TRUE, method = meth, lower = lower,  upper = upper, ...): function cannot be evaluated at initial parameters> 
Error in fitdist(data, "plcon", start = list(xmin = 1, alpha = 1)) : 
    the function mle failed to estimate the parameters, 
       with the error code 100 

我嘗試在谷歌和計算器進行搜索，所以很多類似錯誤的問題出現，但閱讀和嘗試後，沒有解決方案，在工作我的問題，我應該怎麼做才能正確地完成它以獲取參數？ 謝謝大家幫我個忙！

Answer 1

這是一個有趣的發現，我並不完全滿意這個發現，但我會告訴你我發現了什麼，看看它是否有幫助。

在調用fitdist函數時，默認情況下它想要使用來自同一包中的mledist。這本身導致對stats::optim的調用，這是一種通用優化功能。在它的返回值中，它提供了一個收斂錯誤代碼，詳情請參閱?optim。您看到的100不是由optim返回的那個。因此，我拆開了代碼mledist和fitdist以查找錯誤代碼的來源。不幸的是，它被定義在多個案例中，並且是一般陷阱錯誤代碼。如果你分解了所有的代碼，那麼fitdist在這裏試圖做的是如下，預先進行各種檢查等。

fnobj <- function(par, fix.arg, obs, ddistnam) { 
    -sum(do.call(ddistnam, c(list(obs), as.list(par), 
          as.list(fix.arg), log = TRUE))) 
} 

vstart = list(xmin=5,alpha=5) 
fnobj <- function(par, fix.arg obs, ddistnam) { 
    -sum(do.call(ddistnam, c(list(obs), as.list(par), 
          as.list(fix.arg), log = TRUE))) 
} 
ddistname=dplcon 
fix.arg = NULL 
meth = "Nelder-Mead" 
lower = -Inf 
upper = Inf 
optim(par = vstart, fn = fnobj, 
     fix.arg = fix.arg, obs = data, ddistnam = ddistname, 
     hessian = TRUE, method = meth, lower = lower, 
     upper = upper) 

如果我們運行這段代碼，我們找到更多有用的錯誤「功能不能在初始參數進行評估」。如果我們看看函數定義，這是有道理的。有xmin=0或alpha=1將產生-Inf的對數似然。好吧，想想嘗試不同的初始值，我嘗試了幾個隨機選擇，但都返回了一個新的錯誤，「非有限差分值1」。

搜索optim來源進一步爲這兩個錯誤的來源，他們不是R源本身的一部分，但是有一個.External2調用，所以我只能假設錯誤來自那裏。非有限錯誤意味着某個函數評估的某個地方給出了非數字結果。當alpha <= 1或xmin <= 0時，功能dplcon將這樣做。fitdist可讓您指定傳遞給mledist或其他參數的其他參數（取決於您選擇的方法，mle是默認值），其中lower是用於控制要優化的參數的下限的參數。所以，我想強加這些限制，並再次嘗試：

fitpl <- fitdist(data,"plcon",start = list(xmin=1,alpha=2), lower = c(xmin = 0, alpha = 1)) 

煩人，這仍然給出了一個錯誤代碼100跟蹤下來產生錯誤「L-BFGS-B需要的‘FN’有限值」。優化方法已從默認的Nelder-Mead更改爲您指定的邊界，並在外部C代碼調用此錯誤出現的某處，可能接近xmin或alpha的限制，其中數值計算的穩定性接近無限時爲重要。

，我決定做位數的匹配，而不是最大的可能性，試圖找出更多

fitpl <- fitdist(data,"plcon",start = list(xmin=1,alpha=2), 
    method= "qme",probs = c(1/3,2/3)) 
fitpl 
## Fitting of the distribution ' plcon ' by matching quantiles 
## Parameters: 
##   estimate 
## xmin 0.02135157 
## alpha 46.65914353 

這表明的xmin最佳值接近於0，這是極限。我不滿意的原因是我無法使用fitdist獲得分佈的最大似然擬合，但希望這個解釋有幫助，而分位數匹配提供了另一種選擇。

編輯：

學習多一點有關電源律分佈在總體上是有道理的，因爲你認爲這種不工作後。參數功率參數具有似然函數，其可以在給定xmin的條件下最大化。然而，由於似然函數在xmin中增加，所以xmin不存在這樣的表達式。通常估計xmin來自Kolmogorov - Smirnov統計量，請參閱this mathoverflow問題和poweRlaw包的d_jss_paper vignette以獲取更多信息和相關參考。

有功能可以估計poweRlaw程序包本身的冪律分佈參數。

m = conpl$new(data) 
xminhat = estimate_xmin(m)$xmin 
m$setXmin(xminhat) 
alphahat = estimate_pars(m)$pars 
c(xmin = xminhat, alpha = alphahat)