Eu tenho um problema de otimização que oNelder-Mead método vai resolver, mas que eu também gostaria de resolver usandoBFGS ou Newton-Raphson, ou algo que tenha uma função de gradiente, para mais velocidade e, esperançosamente, estimativas mais precisas. Eu escrevi essa função de gradiente seguindo (eu pensei) o exemplo nooptim / optimx documentação, mas quando eu usá-lo comBFGS meus valores iniciais ou não se movem (optim()), ou então a função completa não é executada (optimx(), que retornaError: Gradient function might be wrong - check it!). Me desculpe, há um pouco de código envolvido em reproduzir isso, mas aqui vai:

Esta é a função para a qual quero obter estimativas de parâmetro (isto é para suavizar as taxas de mortalidade de idosos, onde x é idade, começando aos 80 anos):

    KannistoMu <- function(pars, x = .5:30.5){
      a <- pars["a"]
      b <- pars["b"]
      (a * exp(b * x)) / (1 + a * exp(b * x))
    }

E aqui está uma função de verossimilhança para estimar a partir das taxas observadas (definidas como mortes,.Dx superexposição,.Exp):

    KannistoLik1 <- function(pars, .Dx, .Exp, .x. = .5:30.5){
      mu <- KannistoMu(exp(pars), x = .x.)
      # take negative and minimize it (default optimizer behavior)
      -sum(.Dx * log(mu) - .Exp * mu, na.rm = TRUE) 
    }

entendeexp(pars) lá porque eu doulog(pars) para otimizar ao longo, a fim de restringir o finala eb ser positivo.

Dados de exemplo (1962 fêmeas do Japão, se alguém estiver curioso):

    .Dx <- structure(c(10036.12, 9629.12, 8810.11, 8556.1, 7593.1, 6975.08, 
      6045.08, 4980.06, 4246.06, 3334.04, 2416.03, 1676.02, 1327.02, 
      980.02, 709, 432, 350, 217, 134, 56, 24, 21, 10, 8, 3, 1, 2, 
      1, 0, 0, 0), .Names = c("80", "81", "82", "83", "84", "85", "86", 
      "87", "88", "89", "90", "91", "92", "93", "94", "95", "96", "97", 
      "98", "99", "100", "101", "102", "103", "104", "105", "106", 
      "107", "108", "109", "110"))
    .Exp <- structure(c(85476.0333333333, 74002.0866666667, 63027.5183333333, 
      53756.8983333333, 44270.9, 36749.85, 29024.9333333333, 21811.07, 
      16912.315, 11917.9583333333, 7899.33833333333, 5417.67, 3743.67833333333, 
      2722.435, 1758.95, 1043.985, 705.49, 443.818333333333, 223.828333333333, 
      93.8233333333333, 53.1566666666667, 27.3333333333333, 16.1666666666667, 
      10.5, 4.33333333333333, 3.16666666666667, 3, 2.16666666666667, 
      1.5, 0, 1), .Names = c("80", "81", "82", "83", "84", "85", "86", 
      "87", "88", "89", "90", "91", "92", "93", "94", "95", "96", "97", 
      "98", "99", "100", "101", "102", "103", "104", "105", "106", 
      "107", "108", "109", "110"))

Os seguintes trabalhos para oNelder-Mead método:

    NMab <- optim(log(c(a = .1, b = .1)), 
      fn = KannistoLik1, method = "Nelder-Mead",
      .Dx = .Dx, .Exp = .Exp)
    exp(NMab$par) 
    # these are reasonable estimates
       a         b 
    0.1243144 0.1163926

Esta é a função gradiente que eu criei:

    Kannisto.gr <- function(pars, .Dx, .Exp, x = .5:30.5){
      a <- exp(pars["a"])
      b <- exp(pars["b"])
      d.a <- (a * exp(b * x) * .Exp + (-a * exp(b * x) - 1) * .Dx) /
        (a ^ 3 * exp(2 * b * x) + 2 * a ^ 2 * exp(b * x) + a)
      d.b <- (a * x * exp(b * x) * .Exp + (-a * x * exp(b * x) - x) * .Dx) /
        (a ^ 2 * exp(2 * b * x) + 2 * a * exp(b * x) + 1)
      -colSums(cbind(a = d.a, b = d.b), na.rm = TRUE)
    }

A saída é um vetor de comprimento 2, a mudança em relação aos parâmetrosa eb. Eu também tenho uma versão mais feia chegou explorando a saída dederiv(), que retorna a mesma resposta, e que eu não posto (apenas para confirmar que as derivadas estão certas).

Se eu fornecer paraoptim() como segue, comBFGS como o método, as estimativas não se movem dos valores iniciais:

    BFGSab <- optim(log(c(a = .1, b = .1)), 
      fn = KannistoLik1, gr = Kannisto.gr, method = "BFGS",
      .Dx = .Dx, .Exp = .Exp)
    # estimates do not change from starting values:
    exp(BFGSab$par) 
      a   b 
    0.1 0.1

Quando eu olho para o$counts elemento da saída, diz queKannistoLik1() foi chamado 31 vezes eKannisto.gr() apenas 1 vez.$convergence é0, então eu acho que isso converge (se eu der início menos razoável, eles também ficam). Reduzi a tolerância, etc, e nada muda. Quando eu tento a mesma ligaçãooptimx() (não mostrado), recebo a advertência que mencionei acima e nenhum objeto é retornado. Eu obtenho os mesmos resultados ao especificargr = Kannisto.gr com o"CG". Com o"L-BFGS-B" Eu recebo os mesmos valores iniciais como estimativa, mas também é relatado que tanto a função quanto o gradiente foram chamados 21 vezes, e há uma mensagem de erro:"ERROR: BNORMAL_TERMINATION_IN_LNSRCH"

Eu estou esperando que haja algum pequeno detalhe na forma como a função de gradiente é escrita que irá resolver isso, como este aviso posterior eooptimx Comportamento são francamente insinuando que a função simplesmente não está certa (eu acho). Eu também tentei omaxNR() maximizador domaxLik pacote e observado comportamento semelhante (valores iniciais não se movem). Alguém pode me dar um ponteiro? Muito grato

[Edit] @Vincent sugeriu comparar com a saída de uma aproximação numérica:

    library(numDeriv)
    grad( function(u) KannistoLik1( c(a=u[1], b=u[2]), .Dx, .Exp ), log(c(.1,.1)) )
    [1] -14477.40  -7458.34
    Kannisto.gr(log(c(a=.1,b=.1)), .Dx, .Exp)
     a        b 
    144774.0  74583.4 

sinal tão diferente e desligado por um fator de 10? Eu mudo a função gradiente para seguir o exemplo:

    Kannisto.gr2 <- function(pars, .Dx, .Exp, x = .5:30.5){
      a <- exp(pars["a"])
      b <- exp(pars["b"])
      d.a <- (a * exp(b * x) * .Exp + (-a * exp(b * x) - 1) * .Dx) /
        (a ^ 3 * exp(2 * b * x) + 2 * a ^ 2 * exp(b * x) + a)
      d.b <- (a * x * exp(b * x) * .Exp + (-a * x * exp(b * x) - x) * .Dx) /
        (a ^ 2 * exp(2 * b * x) + 2 * a * exp(b * x) + 1)
      colSums(cbind(a=d.a,b=d.b), na.rm = TRUE) / 10
    }
    Kannisto.gr2(log(c(a=.1,b=.1)), .Dx, .Exp)
    # same as numerical:
      a         b 
    -14477.40  -7458.34 

Experimente no otimizador:

    BFGSab <- optim(log(c(a = .1, b = .1)), 
      fn = KannistoLik1, gr = Kannisto.gr2, method = "BFGS",
      .Dx = .Dx, .Exp = .Exp)
    # not reasonable results:
    exp(BFGSab$par) 
      a   b 
    Inf Inf 
    # and in fact, when not exp()'d, they look oddly familiar:
    BFGSab$par
      a         b 
    -14477.40  -7458.34 

Seguindo a resposta de Vincent, eu fiz o reescalonamento da função gradiente e useiabs() ao invés deexp() para manter os parâmetros positivos. As mais recentes e melhores funções de objetivo e gradiente:

    KannistoLik2 <- function(pars, .Dx, .Exp, .x. = .5:30.5){
      mu <- KannistoMu.c(abs(pars), x = .x.)
      # take negative and minimize it (default optimizer behavior)
      -sum(.Dx * log(mu) - .Exp * mu, na.rm = TRUE) 
    }

    # gradient, to be down-scaled in `optim()` call
    Kannisto.gr3 <- function(pars, .Dx, .Exp, x = .5:30.5){
      a <- abs(pars["a"])
      b <- abs(pars["b"])
      d.a <- (a * exp(b * x) * .Exp + (-a * exp(b * x) - 1) * .Dx) /
        (a ^ 3 * exp(2 * b * x) + 2 * a ^ 2 * exp(b * x) + a)
      d.b <- (a * x * exp(b * x) * .Exp + (-a * x * exp(b * x) - x) * .Dx) /
        (a ^ 2 * exp(2 * b * x) + 2 * a * exp(b * x) + 1)
      colSums(cbind(a = d.a, b = d.b), na.rm = TRUE) 
    }

    # try it out:
    BFGSab2 <- optim(
      c(a = .1, b = .1), 
      fn = KannistoLik2, 
      gr = function(...) Kannisto.gr3(...) * 1e-7, 
      method = "BFGS",
      .Dx = .Dx, .Exp = .Exp
    )
    # reasonable:
    BFGSab2$par
            a         b 
    0.1243249 0.1163924 

    # better:
    KannistoLik2(exp(NMab1$par),.Dx = .Dx, .Exp = .Exp) > KannistoLik2(BFGSab2$par,.Dx = .Dx, .Exp = .Exp)
    [1] TRUE

Isso foi resolvido muito mais rápido do que eu esperava, e aprendi mais que alguns truques. Obrigado Vincent!