Wie man einen numerischen Gradienten in constrOptim einfügt

Question

Ich habe Maximum-Likelihood-Schätzung mit optim() gemacht und es war ziemlich einfach. Es ist eine verallgemeinerte logistische Verteilung mit 4 Parametern und ein paar Einschränkungen, die alle in der Likelihood-Funktion aufgeführt:

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

Ich habe es Likelihood-Funktion ist und arbeitete flawessly auf diese Weise:

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) 

Da diese Funktion tut Den Gradienten hatte ich numerisch nicht viel problem.

Dann wollte ich auf constrOptim() ändern, weil die Grenzen tatsächlich 0 sind und nicht eine kleine Zahl auf den ersten 3 Parametern. Aber, das Problem, das ich gegenüberstelle, ist, dass das Argument grad angegeben werden muss und ich kann diese Funktion nicht ableiten, um eine Gradientenfunktion zu geben, also muss ich es numerisch wie in optim() tun, es funktioniert, wenn ich grad = NULL setze, aber ich don ' Ich möchte Nelder-Mead-Methode, aber BFGS.

Ich habe auf diese Weise versucht, aber nicht viel Erfolg:

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

Ihre Schreibweise ist ein bisschen kompliziert, vielleicht, dass Sie verwirrt.

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)