cách xác định hàm gradient để sử dụng trong optim() hoặc trình tối ưu hóa khác

Question

Tôi có vấn đề về tối ưu hóa mà phương pháp Nelder-Mead sẽ giải quyết, nhưng tôi cũng muốn giải quyết bằng cách sử dụng BFGS hoặc Newton-Raphson, hoặc thứ gì đó một hàm gradient, cho tốc độ cao hơn và hy vọng các ước tính chính xác hơn. Tôi đã viết một hàm gradient như sau (tôi nghĩ) ví dụ trong tài liệu optim/optimx, nhưng khi tôi sử dụng nó với BFGS giá trị bắt đầu của tôi không di chuyển (optim()), hoặc chức năng khác hoàn toàn không chạy (optimx()) , trả về Error: Gradient function might be wrong - check it!). Tôi xin lỗi có một chút mã liên quan đến việc tái tạo điều này, nhưng ở đây đi:

Đây là hàm tôi muốn lấy ước tính tham số (điều này là để làm mịn tỷ lệ tử vong ở độ tuổi già, trong đó x là tuổi, bắt đầu từ 80 tuổi):

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

và đây là một hàm log likelihood cho việc ước tính nó từ tỷ lệ quan sát (định nghĩa là trường hợp tử vong, .Dx qua tiếp xúc, .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) 
    } 

bạn thấy exp(pars) trong đó BECA sử dụng tôi cung cấp cho log(pars) để tối ưu hóa, để hạn chế a và b cuối cùng là dương.

dữ liệu Ví dụ (1962 Nhật Bản nữ, nếu có ai là tò mò):

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

Các công trình sau đây đối với phương pháp 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 

Đây là chức năng Gradient tôi đến với:

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

Kết quả là một vectơ có chiều dài 2, thay đổi liên quan đến p đường kính a và b. Tôi cũng có một phiên bản xấu hơn đến bằng cách khai thác đầu ra của deriv(), trả về cùng một câu trả lời và tôi không đăng (chỉ để xác nhận rằng các dẫn xuất là đúng).

Nếu tôi cung cấp nó để optim() như sau, với BFGS như phương pháp, dự toán không di chuyển từ các giá trị khởi đầu:

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 

Khi tôi nhìn vào các yếu tố $counts của đầu ra, nó nói rằng KannistoLik1() được gọi là 31 lần và Kannisto.gr() chỉ 1 lần. $convergence là 0, vì vậy tôi đoán nó nghĩ rằng nó hội tụ (nếu tôi đưa ra khởi đầu ít hợp lý hơn, họ cũng sẽ đặt). Tôi giảm dung sai, vv, và không có gì thay đổi. Khi tôi thử cùng một cuộc gọi trong số optimx() (không được hiển thị), tôi nhận được cảnh báo tôi đã đề cập ở trên và không có đối tượng nào được trả lại. Tôi nhận được kết quả tương tự khi chỉ định gr = Kannisto.gr với số "CG".Với phương pháp "L-BFGS-B" tôi nhận được các giá trị khởi đầu cùng trở lại như ước tính, nhưng nó cũng được thông báo rằng cả hai chức năng và độ dốc được gọi là 21 lần, và có một thông báo lỗi: "ERROR: BNORMAL_TERMINATION_IN_LNSRCH"

Tôi hy vọng rằng có một số chi tiết nhỏ trong cách chức năng gradient được viết mà sẽ giải quyết điều này, như cảnh báo sau này và hành vi optimx được thẳng thừng gợi ý rằng chức năng đơn giản là không đúng (tôi nghĩ). Tôi cũng đã thử trình tối ưu hóa maxNR() từ gói maxLik và đã quan sát hành vi tương tự (giá trị bắt đầu không di chuyển). Bất cứ ai có thể cho tôi một con trỏ? Nhiều nghĩa vụ

[Chỉnh sửa] @Vincent đề nghị tôi so sánh với sản lượng từ một xấp xỉ số:

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 

dấu hiệu rất khác nhau, và tắt bằng một yếu tố của 10? Tôi thay đổi chức năng gradient để làm theo phù hợp:

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 

Hãy thử nó trong tôi ưu hoa:

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 

Tiếp theo câu trả lời của Vincent, tôi rescaled chức năng gradient, và sử dụng abs() thay vì exp() để giữ các thông số tích cực. Các chức năng mục tiêu và gradient mới nhất và hiệu quả nhất:

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 

Điều này đã được giải quyết nhanh hơn nhiều so với mong đợi và tôi đã học được nhiều hơn một vài thủ thuật. Cảm ơn Vincent!

Answer 1

Để kiểm tra xem gradient là đúng, bạn có thể so sánh nó với một xấp xỉ số:

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) 

Các dấu hiệu là sai: các thuật toán không thấy bất kỳ cải thiện khi nó di chuyển theo hướng này, và do đó không di chuyển.

Bạn có thể sử dụng một số hệ thống đại số máy tính (ở đây, Maxima) để làm phép tính dành cho bạn:

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

tôi chỉ cần sao chép và dán kết quả vào 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 

Nếu bạn mù quáng đặt gradient trong tối ưu hóa, có một vấn đề không ổn định số: giải pháp được đưa ra là (Inf,Inf) ... Để ngăn chặn nó, bạn có thể rescale gradient (cách giải quyết tốt hơn là sử dụng phép biến đổi ít nổ hơn số mũ, để đảm bảo rằng các thông số vẫn dương).

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