3

rplcon()パッケージ内の関数を使用していくつかの確率変数を生成しますpoweRlaw

data <- rplcon(1000,10,2)

ここで、データに最もよく適合する既知の分布を知りたいと思います。ログノルム？経験？ガンマ？べき法則？指数カットオフのべき乗則?

だから私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によると、べき乗則分布と指数カットオフを伴うべき乗則は基本確率関数ではないため、次の例 4 に基づいてべき乗則の d,p,q 関数を記述します。?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))

最後に、以下のコードを使用してパラメーターxminと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

Google と stackoverflow で検索しようとすると、同様のエラーの質問がたくさん表示されますが、読んで試しても解決策がありません。パラメーターを取得するために正しく完了するにはどうすればよいですか? 私に好意を寄せてくれたみんなに感謝します！

4

1 に答える 1

5

これは興味深い発見であり、私はこの発見に完全に満足しているわけではありませんが、私が発見したことをお話しし、それが役立つかどうかを確認します.

関数を呼び出すと、fitdistデフォルトでmledistは同じパッケージから使用したいと考えています。stats::optimこれ自体が、一般的な最適化関数である呼び出しになります。戻り値では、収束エラー コードが返されます。詳細については、 を参照?optimしてください。表示される100は、 によって返されるものの 1 つではありません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これら2つのエラーのソースをさらに検索すると、Rソース自体の一部ではありませんが、.External2呼び出しがあるため、そこからエラーが発生したとしか思えません。非有限エラーは、どこかの関数評価の 1 つが数値以外の結果を与えることを意味します。関数は、またはのdplconときにそうします。最適化するパラメーターの下限を制御するための引数またはその他 (選択したメソッドに応じて、mle がデフォルト) に渡される追加の引数を指定できます。だから私はこれらの制限を課して再試行しようとしました:alpha <= 1xmin <= 0fitdistmledistlower

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

厄介なことに、これでもエラー コード 100 が返されます。これを追跡すると、「L-BFGS-B には 'fn' の有限値が必要です」というエラーが表示されます。xmin境界を指定すると、最適化方法がデフォルトの Nelder-Mead から変更され、外部 C コード呼び出しのどこかでこのエラーが発生します。おそらく、またはの限界に近い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 統計から得られます。詳細と関連する参考文献については、このmathoverflow の質問と poweRlaw パッケージの d_jss_paper ビネットを参照してください。

poweRlawパッケージ自体にべき乗分布のパラメータを推定する機能があります。

m = conpl$new(data)
xminhat = estimate_xmin(m)$xmin
m$setXmin(xminhat)
alphahat = estimate_pars(m)$pars
c(xmin = xminhat, alpha = alphahat)
于 2016-05-11T08:40:02.980 に答える