jak określić konkretną funkcję gradientu do użycia w optym() lub innym optymalizatorze

Question

Mam problem z optymalizacją, który rozwiąże metoda Nelder-Mead, ale chciałbym również rozwiązać za pomocą BFGS lub Newtona-Raphsona, lub coś, co zajmuje funkcja gradientu, dla większej prędkości i, miejmy nadzieję, dokładniejszych oszacowań. Napisałem taką funkcję gradientu po (pomyślałem) przykładzie z dokumentacji optim/optimx, ale kiedy używam go z BFGS, moje wartości początkowe albo się nie poruszają (optim()), albo też funkcja nie działa (optimx() , która zwraca Error: Gradient function might be wrong - check it!). Przykro mi, że jest trochę kodu zaangażowanego w odtworzenie tego, ale tutaj:

Jest to funkcja, dla której chcę uzyskać prognozy parametrów (służy to do wygładzania współczynników umieralności w podeszłym wieku, gdzie x to wiek, począwszy od lat 80):

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

a tu funkcję dziennika wiarygodności oszacowania go z zaobserwowanych stóp (definiowane jako zgonów .Dx nad ekspozycją, .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) 
    } 

widać exp(pars) w tam beca użyj I dać log(pars), aby zoptymalizować, aby ograniczyć ostateczne a i b, aby były pozytywne.

przykład dane (1962 Japonia samice, jeśli ktoś jest ciekawy):

.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")) 

następujące prace dotyczące sposobu Nelder-Mead:

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 

Ta funkcja gradient wpadł:

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

Wyjście to wektor o długości 2, zmiana w stosunku do p Arameters a i b. Mam też brzydszą wersję, wykorzystującą wyjście deriv(), które zwraca tę samą odpowiedź, a której nie publikuję (tylko po to, aby potwierdzić, że instrumenty pochodne mają rację).

Gdybym dostarczyć go do optim() następująco, z BFGS jako metody oszacowania nie poruszają się od wartości wyjściowych:

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 

Kiedy patrzę na elemencie $counts wyjścia, to mówi, że KannistoLik1() został wywołany 31 razy i Kannisto.gr() tylko raz. $convergence jest 0, więc domyślam się, że to się łączy (jeśli daję mniej rozsądnych początków, oni też zostaną). Zmniejszyłem tolerancję, itd. I nic się nie zmieniło. Kiedy próbuję wykonać to samo połączenie w optimx() (nie pokazano), otrzymuję warowanie, o którym wspomniałem powyżej, i żaden obiekt nie jest zwracany. Te same wyniki uzyskuję, określając gr = Kannisto.gr przy pomocy.Dzięki metodzie "L-BFGS-B" dostaję te same wartości wyjściowych z powrotem jako oszacowania, ale jest również, że zarówno funkcja i jej nachylenia, nazywane były 21 razy, a tam jest komunikat o błędzie: "ERROR: BNORMAL_TERMINATION_IN_LNSRCH"

Mam nadzieję, że nie jest pewne pomniejsze szczegóły w sposobie, w jaki napisana jest funkcja gradientu, która rozwiąże to, ponieważ to późniejsze ostrzeżenie i zachowanie optimx wskazują wprost, że funkcja po prostu nie jest odpowiednia (chyba). Próbowałem również maksymalizator maxNR() z pakietu maxLik i zaobserwowałem podobne zachowanie (wartości początkowe się nie poruszają). Czy ktoś może mi wskazać? Bardzo zobowiązany

[Edytuj] @Vincent zaproponował mi porównać z wyjściem ze zbliżenia liczbowym:

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 

więc inny znak, i wyłącza o czynnik 10? I zmienić funkcję gradientu do naśladownictwa:

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 

Spróbuj w optymalizator:

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 

następującej odpowiedzi Vincenta, to przeskalowane funkcji gradientu i używane abs() zamiast exp() aby zachować parametry pozytywne. Najnowsze i lepsze wykonywanie obiektywne i gradientu funkcje:

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 

ten został rozwiązany znacznie szybciej niż się spodziewałem, a ja nauczyłem się więcej niż kilka sztuczek. Dzięki Vincent!

Answer 1

Aby sprawdzić, czy gradient jest poprawna, można porównać go z zbliżenia numerycznej:

library(numDeriv); 
grad(function(u) KannistoLik1(c(a=u[1], b=u[2]), .Dx, .Exp), c(1,1)); 
Kannisto.gr(c(a=1,b=1), .Dx, .Exp) 

Znaki są źle: algorytm nie widzi żadnej poprawy gdy porusza się w tym kierunku, dlatego się nie porusza.

Można użyć jakiegoś systemu algebry komputerowej (tutaj maxima) wykonać obliczenia dla Ciebie:

display2d: false; 
f(a,b,x) := a * exp(b*x)/(1 + a * exp(b*x)); 
l(a,b,d,e,x) := - d * log(f(a,b,x)) + e * f(a,b,x); 
factor(diff(l(exp(a),exp(b),d,e,x),a)); 
factor(diff(l(exp(a),exp(b),d,e,x),b)); 

Po prostu skopiuj i wklej wynik do R:

f_gradient <- function(u, .Dx, .Exp, .x.=.5:30.5) { 
    a <- u[1] 
    b <- u[1] 
    x <- .x. 
    d <- .Dx 
    e <- .Exp 
    c(
    sum((e*exp(exp(b)*x+a)-d*exp(exp(b)*x+a)-d)/(exp(exp(b)*x+a)+1)^2), 
    sum(exp(b)*x*(e*exp(exp(b)*x+a)-d*exp(exp(b)*x+a)-d)/(exp(exp(b)*x+a)+1)^2) 
) 
} 

library(numDeriv) 
grad(function(u) KannistoLik1(c(a=u[1], b=u[2]), .Dx, .Exp), c(1,1)) 
f_gradient(c(a=1,b=1), .Dx, .Exp) # Identical 

Jeśli na ślepo umieścić gradient w optymalizacji, jest numeryczny problem niestabilności: podane rozwiązanie to (Inf,Inf) ... Aby temu zapobiec, można zmienić skalę gradientu (lepszym rozwiązaniem byłoby użycie mniej wybuchowej transformacji niż wykładnicza, aby zapewnić, że parametry pozostaną dodatnie).

BFGSab <- optim(
    log(c(a = .1, b = .1)), 
    fn = KannistoLik1, 
    gr = function(...) f_gradient(...) * 1e-3, 
    method = "BFGS", 
    .Dx = .Dx, .Exp = .Exp 
) 
exp(BFGSab$par) # Less precise than Nelder-Mead