constrOptim 내부에 숫자 그라디언트를 삽입하는 방법

Question

optim()을 사용하여 최대 우도 추정을하고 있었고 매우 쉽습니다. 나는 그것이 우도 함수의 제작 및 flawessly 이런 식으로 일을

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) 
} 

:이 함수가 있기 때문에

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) 

을 그것은 모든 가능성 기능에 나열된 4 개 매개 변수로 일반화 물류 유통 및 제한의 커플입니다 그라데이션 수치는별로 문제가 없었습니다.

그런 다음 처음 세 개의 매개 변수에서 경계가 실제로 0이고 작은 수가 아니기 때문에 constrOptim()으로 변경하려고했습니다. 하지만, 내가 직면하는 문제는 인수가 grad이 지정되어야하고 그라디언트 함수를 제공하는 함수를 파생 할 수 없다는 것입니다. 따라서 optim()으로 수치 적으로해야합니다. grad = NULL을 넣으면 작동하지만, 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)