Comment insérer un gradient numérique dans constrOptim

Question

Je faisais une estimation du maximum de vraisemblance en utilisant optim() et c'était assez facile. Il est une distribution logistique généralisée avec 4 paramètres et quelques restrictions, tous énumérés dans la fonction de vraisemblance:

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

Je l'ai fait est fonction de vraisemblance et travaillé flawessly ainsi:

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) 

Comme cette fonction ne le gradient numériquement je n'avais pas beaucoup de problème.

Puis je voulais changer à constrOptim() parce que les limites sont en réalité 0 et pas un petit nombre sur les 3 premiers paramètres. Mais, le problème que je rencontre est que l'argument grad doit être spécifié et je ne peux pas dériver cette fonction pour donner une fonction de gradient, donc je dois le faire numériquement comme dans optim(), ça marche si je mets grad = NULL mais je ne fais pas Je veux la méthode de Nelder-Mead mais BFGS.

J'ai essayé de cette façon, mais pas de beaucoup sucess:

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

Votre notation est un peu compliqué, peut-être que vous confondre.

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)