constrOptimの中に数値的なグラデーションを挿入するには

Question

optim（）を使って最尤推定を行っていて、とても簡単でした。これは、尤度関数に記載されているすべての、4つのパラメータと制限のカップルと一般的なロジスティック分布です：

genlogis.loglikelihood <- function(param = c(sqrt(2/pi),0.5, 2, 0), x){ 

    if(length(param) < 3 | length(param) > 4){ 
    stop('Incorrect number of parameters: param = c(a,b,p,location)') 
    } 

    if(length(param) == 3){ 
    #warning('Location parameter is set to 0') 
    location = 0 
    } 

    if(length(param) == 4){ 
    location = param[4] 
    } 

    a = param[1] 
    b = param[2] 
    p = param[3] 

    if(!missing(a)){ 
    if(a < 0){ 
     stop('The argument "a" must be positive.') 
    } 
    } 
    if(!missing(b)){ 
    if(b < 0){ 

     stop('The argument "b" must be positive.') 
    } 
    } 
    if(!missing(p)){ 
    if(p < 0){ 
     stop('The argument "p" must be positive.') 
    } 
    } 

    if(p == 0 && b > 0 && a > 0){ 
    stop('If "p" equals to 0, "b" or "a" must be 
     0 otherwise there is identifiability problem.') 
    } 
    if(b == 0 && a == 0){ 
    stop('The distribution is not defined for "a" 
     and "b" equal to 0 simultaneously.') 
    } 

    z <- sum(log((a+b*(1+p)*abs((x-location))^p) * exp(-((x-location)*(a+b*abs((x-location))^p))))) - 
      sum(2*log(exp(-((x-location)*(a+b*abs((x-location))^p))) + 1)) 
    if(!is.finite(z)){ 
    z <- 1e+20 
    } 

    return(-z) 
} 

私はそれが尤度関数だ作られたとflawesslyこの方法を働いた：

opt <- function(parameters, data){ 
      optim(par = parameters, fn = genlogis.loglikelihood, x=data, 
        lower = c(0.00001,0.00001,0.00001, -Inf), 
        upper = c(Inf,Inf,Inf,Inf), method = 'L-BFGS-B') 
     } 
opt(c(0.3, 1.01, 2.11, 3.5), faithful$eruptions) 

この機能はありませんので、勾配数値的にはそれほど問題はありませんでした。境界は実際には0ではなく最初の3つのパラメータの数が少ないですので

は、それから私は、constrOptim()に変更するようでした。しかし、私が直面する問題は、「私はgrad = NULLを置くが、私はドン場合、それは動作しますが、引数gradを指定する必要があり、私は勾配関数を与えるために、その関数を導出することはできませんので、私は数値的optim()のようにそれをしなければならないということですNelder-Mead法ではなくBFGS法が必要です。

私はあまりない成功事例のこの方法を試してみた：

opt2 <- function(initial, data){ 
    ui <- rbind(c(1, 0, 0, 0), c(0,1,0,0), c(0,0,1,0)) 
    ci <- c(0,0,0)  
      constrOptim(theta = initial, f = genlogis.loglikelihood(param, x), 
         grad = numDeriv::grad(func = function(x, param) genlogis.loglikelihood(param, x), param = theta, x = data) 
         , x = data, ui = ui, ci = ci) 
     } 

Answer 1

あなたの表記は少し複雑です、多分あなたを混同しています。

opt2 <- function(parameters, data){ 
    fn = function(p) genlogis.loglikelihood(p, x = data) 
    gr = function(p) numDeriv::grad(fn, p) 
    ui <- rbind(c(1, 0, 0, 0), c(0,1,0,0), c(0,0,1,0)) 
    ci <- c(0,0,0)  
    constrOptim(theta = parameters, f = fn, grad = gr, 
       ui = ui, ci = ci, method="BFGS") 
} 
opt2(c(0.3, 1.01, 2.11, 3.5), faithful$eruptions)