Di chuyển ra ngoài chức năng tối ưu của R

Question

Tôi đang cố gắng sử dụng R để ước tính mô hình logit đa thức với thông số kỹ thuật thủ công. Tôi đã tìm thấy một vài gói cho phép bạn ước tính các mô hình MNL here hoặc here.

Tôi đã tìm thấy một số tác phẩm khác về "lăn" riêng MLE chức năng của bạn here. Tuy nhiên, từ việc đào bới của tôi - tất cả các chức năng và gói này đều dựa trên hàm nội bộ optim.

Trong các thử nghiệm benchmark của tôi, optim là nút cổ chai. Sử dụng bộ dữ liệu mô phỏng với ~ 16000 quan sát và 7 tham số, R mất khoảng 90 giây trên máy tính của tôi. Mô hình tương đương trong Biogeme mất ~ 10 giây. Một đồng nghiệp viết mã của riêng mình trong các báo cáo Ox khoảng 4 giây cho cùng một kiểu máy này.

Có ai có kinh nghiệm viết chức năng MLE của riêng mình hay có thể chỉ cho tôi theo hướng được tối ưu hóa vượt quá chức năng mặc định optim (không có ý định chơi chữ)?

Nếu bất cứ ai muốn mã R để tái tạo mô hình, cho tôi biết - tôi glady sẽ cung cấp nó. Tôi chưa cung cấp nó vì nó không liên quan trực tiếp đến vấn đề tối ưu hóa chức năng optim và để bảo toàn không gian ...

EDIT: Cảm ơn mọi người đã suy nghĩ của bạn. Dựa trên vô số bình luận dưới đây, chúng tôi có thể nhận được R trong cùng một ballpark như Biogeme cho các mô hình phức tạp hơn, và R thực sự nhanh hơn cho một số mô hình nhỏ hơn/đơn giản hơn mà chúng tôi chạy. Tôi nghĩ rằng giải pháp lâu dài cho vấn đề này sẽ liên quan đến việc viết một hàm tối đa riêng biệt dựa trên thư viện fortran hoặc C, nhưng chắc chắn mở ra các cách tiếp cận khác.

Answer 1

Cố gắng với nlm() chức năng chưa? Không biết nếu nó nhanh hơn nhiều, nhưng nó cải thiện tốc độ. Ngoài ra kiểm tra các tùy chọn. tối ưu sử dụng thuật toán chậm làm mặc định. Bạn có thể tăng tốc độ gấp 5 lần bằng cách sử dụng thuật toán Quasi-Newton (phương pháp = "BFGS") thay vì mặc định. Nếu bạn không quan tâm quá nhiều về các chữ số cuối cùng, bạn cũng có thể đặt mức dung sai cao hơn nlm() để đạt được tốc độ thêm.

f <- function(x) sum((x-1:length(x))^2) 

a <- 1:5 

system.time(replicate(500, 
    optim(a,f) 
)) 
    user system elapsed 
    0.78 0.00 0.79 

system.time(replicate(500, 
    optim(a,f,method="BFGS") 
)) 
    user system elapsed 
    0.11 0.00 0.11 

system.time(replicate(500, 
    nlm(f,a) 
)) 
    user system elapsed 
    0.10 0.00 0.09 

system.time(replicate(500, 
     nlm(f,a,steptol=1e-4,gradtol=1e-4) 
)) 
    user system elapsed 
    0.03 0.00 0.03 

Answer 2

Bạn có xem tài liệu trên CRAN Task View for Optimization không?

Answer 3

FWIW, tôi đã làm điều này trong C-ish, sử dụng OPTIF9. Bạn sẽ khó ép nhanh hơn thế. Có rất nhiều cách để điều gì đó chạy chậm hơn, chẳng hạn như bằng cách chạy một trình thông dịch như R.

Nhập: Từ nhận xét, rõ ràng OPTIF9 được sử dụng làm công cụ tối ưu hóa. Điều đó có nghĩa là phần lớn thời gian được dùng để đánh giá hàm mục tiêu trong R. Trong khi có thể các hàm C đang được sử dụng bên dưới cho một số thao tác, thì vẫn còn có phí thông dịch. Có một cách nhanh chóng để xác định những dòng mã và các cuộc gọi hàm trong R chịu trách nhiệm cho phần lớn thời gian, và đó là tạm dừng các lệnh đó bằng phím Escape và kiểm tra ngăn xếp. Nếu một tuyên bố tốn X% thời gian, thì nó nằm trên đống X% thời gian. Bạn có thể thấy rằng có các hoạt động sẽ không đi đến C và phải như vậy. Bất kỳ yếu tố tăng tốc nào bạn nhận được theo cách này sẽ được giữ nguyên khi bạn tìm cách để song song với việc thực thi R.