Ir más allá de la función de optimización de R

Question

Estoy tratando de usar R para estimar un modelo logit multinomial con una especificación manual. He encontrado algunos paquetes que le permiten estimar los modelos MNL here o here.

He encontrado algunas otras escrituras en "rodar" tu propia función MLE here. Sin embargo, desde mi búsqueda, todas estas funciones y paquetes se basan en la función interna optim.

En mis pruebas de referencia, optim es el cuello de botella. Usando un conjunto de datos simulado con ~ 16000 observaciones y 7 parámetros, R toma alrededor de 90 segundos en mi máquina. El modelo equivalente en Biogeme demora ~ 10 segundos. Un colega que escribe su propio código en Ox informa alrededor de 4 segundos para este mismo modelo.

¿Alguien tiene experiencia en escribir su propia función MLE o puede dirigirme hacia algo que está optimizado más allá de la función predeterminada optim (sin juego de palabras)?

Si alguien quiere que el código R vuelva a crear el modelo, avíseme: se lo facilitaré glady. No lo he proporcionado, ya que no es directamente relevante para el problema de optimizar la función optim y para preservar el espacio ...

EDITAR: Gracias a todos por su opinión. Basándonos en una miríada de comentarios a continuación, pudimos obtener R en el mismo estadio que Biogeme para modelos más complicados, y R fue en realidad más rápido para varios modelos más pequeños/simples que ejecutamos. Creo que la solución a largo plazo a este problema implicará escribir una función de maximización separada que dependa de una biblioteca fortran o C, pero ciertamente estoy abierto a otros enfoques.

Answer 1

¿Ya funcionó con la función nlm()? No sé si es mucho más rápido, pero mejora la velocidad. También verifica las opciones. optim usa un algoritmo lento como el predeterminado. Puede ganar una aceleración> 5 veces usando el algoritmo Quasi-Newton (method = "BFGS") en lugar del predeterminado. Si no le preocupan demasiado los últimos dígitos, también puede establecer los niveles de tolerancia más altos de nlm() para obtener más velocidad.

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

¿Consideró el material en el CRAN Task View for Optimization?

Answer 3

FWIW, he hecho esto en C-ish, usando OPTIF9. Estaría en apuros para ir más rápido que eso. Hay muchas maneras de que algo vaya más lento, como ejecutar un intérprete como R.

Agregado: De los comentarios, está claro que OPTIF9 se usa como motor de optimización. Eso significa que lo más probable es que se dedique la mayor parte del tiempo a evaluar la función objetivo en R. Si bien es posible que las funciones C se utilicen debajo para algunas de las operaciones, todavía hay una sobrecarga del intérprete. Hay una manera rápida de determinar qué líneas de código y llamadas de función en R son responsables la mayor parte del tiempo, y eso es pausarlo con la tecla Escape y examinar la pila. Si una declaración cuesta X% de tiempo, está en la pila X% del tiempo. Puede encontrar que hay operaciones que no van a C y deberían serlo. Cualquier factor de aceleración que obtienes de esta manera se conservará cuando encuentres una forma de paralelizar la ejecución de R.