Cómo manejar restricciones de límites cuando se usa `nls.lm` en R

mathematical-optimization (2)

En primer lugar, no soy un experto en Matlab y Optimización y nunca he usado R.

No estoy seguro de ver cuál es tu pregunta real, pero tal vez pueda arrojar algo de luz sobre tu perplejidad:

LM es un enfoque ligeramente mejorado de Gauss-Newton: para problemas con varios mínimos locales es muy sensible a los estados iniciales. La inclusión de límites típicamente genera más de esos mínimos.

TRR es similar a LM, pero más robusto. Tiene mejores capacidades para "saltar" de malos mínimos locales. Es bastante factible que se comporte mejor, pero que se comporte peor, que un LM. En realidad explicar por qué es muy difícil. Necesitaría estudiar los algoritmos en detalle y observar cómo se comportan en esta situación.

No puedo explicar la diferencia entre la implementación de Matlab y la de R, pero hay varias extensiones de TRR que quizás Matlab usa y R no. ¿Su enfoque de usar LM y TRR converge alternativamente mejor que TRR solo?

Hice esta pregunta hace un tiempo. No estoy seguro de si debo publicar esto como una respuesta o una nueva pregunta. No tengo una respuesta, pero "resolví" el problema aplicando el algoritmo de Levenberg-Marquardt usando nls.lm en R y cuando la solución está en el límite, ejecuto el algoritmo reflectivo de la región de confianza (TRR, implementado en R ) alejarse de ello. Ahora tengo nuevas preguntas.

Desde mi experiencia, de esta manera, el programa alcanza el nivel óptimo y no es tan sensible a los valores iniciales. Pero este es solo un método práctico para alejarse de los problemas que encuentro usando nls.lm y también otras funciones de optimización en R. Me gustaría saber por qué nls.lm comporta de esta manera para problemas de optimización con restricciones de límites y cómo manejar el restricciones de límites al usar nls.lm en la práctica.

A continuación, proporcioné un ejemplo que ilustra los dos problemas usando nls.lm

Es sensible a los valores iniciales.
Se detiene cuando algún parámetro alcanza el límite.

Un ejemplo reproducible: Dataset de enfoque D

library(devtools) install_github("KineticEval","zhenglei-gao") library(KineticEval) data(FOCUS2006D) km <- mkinmod.full(parent=list(type="SFO",M0 = list(ini = 0.1,fixed = 0,lower = 0.0,upper =Inf),to="m1"),m1=list(type="SFO"),data=FOCUS2006D) system.time(Fit.TRR <- KinEval(km,evalMethod = ''NLLS'',optimMethod = ''TRR'')) system.time(Fit.LM <- KinEval(km,evalMethod = ''NLLS'',optimMethod = ''LM'',ctr=kingui.control(runTRR=FALSE))) compare_multi_kinmod(km,rbind(Fit.TRR$par,Fit.LM$par)) dev.print(jpeg,"LMvsTRR.jpeg",width=480)

Las ecuaciones diferenciales que describen el modelo / sistema son:

"d_parent = - k_parent * parent" "d_m1 = - k_m1 * m1 + k_parent * f_parent_to_m1 * parent"

En la gráfica de la izquierda está el modelo con valores iniciales, y en el centro está el modelo ajustado que usa "TRR" (similar al algoritmo en la función lsqnonlin ), a la derecha está el modelo ajustado que usa "LM" con nls.lm Mirando los parámetros ajustados ( Fit.LM$par ) encontrará que un parámetro ajustado ( f_parent_to_m1 ) está en el límite 1 . Si cambio el valor de inicio para un parámetro M0_parent de 0.1 a 100, nls.lm los mismos resultados usando nls.lm y lsqnonlin . Tengo muchos casos como este.

newpars <- rbind(Fit.TRR$par,Fit.LM$par) rownames(newpars)<- c("TRR(lsqnonlin)","LM(nls.lm)") newpars M0_parent k_parent k_m1 f_parent_to_m1 TRR(lsqnonlin) 99.59848 0.09869773 0.005260654 0.514476 LM(nls.lm) 84.79150 0.06352110 0.014783294 1.000000

Excepto por los problemas anteriores, a menudo sucede que el Hessian devuelto por nls.lm no es invertible (especialmente cuando algunos parámetros están en el límite), por lo que no puedo obtener una estimación de la matriz de covarianza. Por otro lado, el algoritmo "TRR" (en Matlab) casi siempre da una estimación al calcular el jacobiano en el punto de solución. Creo que esto es útil, pero también estoy seguro de que los algoritmos de optimización R (los que he probado) no hicieron esto por una razón. Me gustaría saber si estoy equivocado al utilizar la forma Matlab de calcular la matriz de covarianza para obtener un error estándar para las estimaciones de los parámetros.

Una última nota, lsqnonlin en mi publicación anterior que lsqnonlin Matlab supera las funciones de optimización de R en casi todos los casos. Estaba equivocado. El algoritmo "Trust-Region-Reflective" utilizado en Matlab es, de hecho, más lento (a veces mucho más lento) si también se implementa en R como se puede ver en el ejemplo anterior. Sin embargo, aún es más estable y alcanza una solución mejor que los algoritmos básicos de optimización de la R.

Usando el paquete mkin, puede encontrar los parámetros usando el algoritmo "Puerto" (que también es una especie de algoritmo TRR, por lo que pude ver en su documentación), o el algoritmo "Marq", que usa nls.lm en el fondo. Luego puede usar valores iniciales "normales" o valores iniciales "malos".

library(mkin) packageVersion("mkin")

La versión mkin reciente puede acelerar el proceso considerablemente, ya que compila los modelos a partir del código C generado automáticamente si el compilador está disponible en su sistema (por ejemplo, tiene r-base-dev instalado en Debian / Ubuntu, o Rtools en Windows).

Esto define el modelo:

m <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"), use_of_ff = "max")

Puedes comprobar que las ecuaciones diferenciales son correctas:

cat(m$diffs, sep = "/n")

Luego encajamos en cuatro variantes, Puerto y LM, con o sin M0 fijado a 0.1:

f.Port = mkinfit(m, FOCUS_2006_D) f.Port.M0 = mkinfit(m, FOCUS_2006_D, state.ini = c(parent = 0.1, m1 = 0)) f.LM = mkinfit(m, FOCUS_2006_D, method.modFit = "Marq") f.LM.M0 = mkinfit(m, FOCUS_2006_D, state.ini = c(parent = 0.1, m1 = 0), method.modFit = "Marq")

Entonces nos fijamos en los resultados:

results <- sapply(list(Port = f.Port, Port.M0 = f.Port.M0, LM = f.LM, LM.M0 = f.LM.M0), function(x) round(summary(x)$bpar[, "Estimate"], 5))

cuales son

Port Port.M0 LM LM.M0 parent_0 99.59848 99.59848 99.59848 39.52278 k_parent 0.09870 0.09870 0.09870 0.00000 k_m1 0.00526 0.00526 0.00526 0.00000 f_parent_to_m1 0.51448 0.51448 0.51448 1.00000

Así que podemos ver que el algoritmo de puerto encuentra la mejor solución (según mi conocimiento) incluso con valores de inicio incorrectos. El problema de la velocidad que uno puede tener con modelos más complicados se alivia con la generación automática de código C.