1

極値分析を行っています。さまざまな理由から、fevd パッケージを使用したくありません (最初の理由は、他の方法ではできないことを微調整できるようにしたいからです)。私は自分のコードを書きました。それはほとんど非常に単純で、すべてを解決したと思っていました。しかし、一部のパラメーターの組み合わせでは、対数尤度分析 ( optimに基づく) から得られるヘッシアンは正しくありません。

一歩ずつ進んでいきます。私のコード - またはその選択した部分 - は次のようになります。

# routines for non stationary
Log_lik_GEV <- function(dataIN,scaleIN,shapeIN,locationIN){
    # simply calculate the negative log likelihood value for a set of X and parameters, for the GPD
    #xi, mu, sigma  - xi is the shape parameter, mu the location parameter, and sigma is the scale parameter.
    # shape = xi
    # location = mu
    # scale = beta
    library(fExtremes)
    #dgev   Density of the GEV Distribution, dgev(x, xi = 1, mu = 0, sigma = 1)

    LLvalues <- dgev(dataIN, xi = shapeIN, mu = locationIN, beta = scaleIN) 
    NLL <- -sum(log(LLvalues[is.finite(LLvalues)]))
    return(NLL)
}


function_MLE <- function(par , dataIN){
    scoreLL <- 0
    shape_param <- par[1]
    scale_param <- par[2]
    location_param <- par[3]
    scoreLL <- Log_lik_GEV(dataIN, scale_param, shape_param, location_param)
    if (abs(shape_param) > 0.3) scoreLL <- scoreLL*10000000
    if ((scale_param) <= 0) {
        scale_param <- abs(scale_param)
        par[2] <- abs(scale_param)
        scoreLL <- scoreLL*1000000000
    }
    sum(scoreLL) 
}

kernel_estimation <- function(dati_AM, shape_o, scale_o, location_o) {

     paramOUT <- optim(par = c(shape_o, scale_o, location_o), fn = function_MLE, dataIN = dati_AM, control = list(maxit = 3000, reltol = 0.00000001), hessian = TRUE)

     # calculation std errors
     covmat <- solve(paramOUT$hessian)
     stde <- sqrt(diag(covmat))
     print(covmat)

     print('')

     result <- list(shape_gev =paramOUT$par[1], scale_gev = paramOUT$par[2],location_gev =paramOUT$par[3], var_covar = covmat)

     return(result)
}

場合によっては、すべてがうまく機能します。私のルーチンと fevd ルーチンを実行すると、まったく同じ結果が得られます。場合によっては (shape=-0.29 が非常に強い負/ワイブルの場合の特定のケース)、私のルーチンは負の分散とファンキーなヘシアンを与えます。常に間違っているわけではありませんが、一部のパラメーターの組み合わせは明らかに有効なヘッシアンを与えていません (注: パラメーターはまだ正しく推定されています。つまり、fevd の結果と同じですが、共分散行列は完全にオフになっています)。

2 つの手順のヘッセ行列を比較したこの投稿を見つけましたが、実際、最適化は不安定なようです。ただし、ルーチンで maxLik を単純に置き換えると、まったく収束しません (収束が発生した場合でも)。

 paramOUT = maxLik(function_MLE, start =c(shape_o, scale_o, location_o),   
                          dataIN=dati_AM, method ='NR' )

正しい初期値であっても、異なる初期値を与えようとしましたが、収束しません。

私の例では最適なルーチンが正しく使用されていると思うので、データを提供していません。簡単に言うと、パラメータの組み合わせによっては、数値結果が安定しません。私の質問は：

1) maxLik の使用方法に何か不足がありますか?

2) ヘシアンを抽出できる maxLik 以外の最適化ルーチンはありますか?

ありがとう

4

0 に答える 0