Diferencia entre el algoritmo de Levenberg-Marquardt y la ODR

Pude ajustar curvas a un conjunto de datos x / y usandopico-o-mat, Como se muestra abajo. Eso es un fondo lineal y 10 curvas lorentzianas.

Ya que necesito encajar muchas curvas similares, escribí una rutina de guión ajustada, usandompfit.py, que es un algoritmo de Levenberg-Marquardt. Sin embargo, el ajuste lleva más tiempo y, en mi opinión, es menos preciso que el resultado de peak-o-mat:

Valores iniciales

Resultado ajustado con fondo lineal fijo (valores para el fondo lineal tomado del resultado peak-o-mat)

Ajuste el resultado con todas las variables gratis.

Creo que los valores iniciales ya están muy cerca, pero incluso con el fondo lineal fijo, el lorentziano izquierdo es claramente una degradación del ajuste.

El resultado es aún peor para el ajuste libre total.

Peak-o-mat parece utilizarscipy.odr.odrpack. Ahora que es más probable:

¿Hice algún error de implementación?¿Odrpack es más adecuado para este problema en particular?

Ajustarse a un problema más simple (datos lineales con un pico en el medio) muestra una muy buena correlación entre peak-o-mat y mi script. Tampoco encontré mucho sobre ordpack.

Editar: Parece que podría responder la pregunta yo solo, sin embargo, la respuesta es un poco inquietante. El uso de scipy.odr (que permite el ajuste con el método odr o lesssq) da el resultado como peak-o-mat, incluso sin restricciones.

La imagen de abajo muestra nuevamente los datos, los valores de inicio (casi perfectos) y luego los ajustes Odr y Lesssq se ajustan. Las curvas de componentes son para el odr one.

Cambiaré a ODR, pero esto todavía me deja molesto. Los métodos (mpfit.py, scipy.optimize.leastsq, scipy.odr en modo lesssq) 'deberían' dar los mismos resultados.

Y para las personas que se encuentran en esta publicación: para hacer el ajuste de odr, se debe especificar un error para los valores x e y. Si no hay ningún error, use valores pequeños con sx << sy.

linear = odr.Model(f)
mydata = odr.RealData(x, y, sx = 1e-99, sy = 0.01)
myodr = odr.ODR(mydata, linear, beta0 = beta0, maxit = 2000)
myoutput1 = myodr.run() 

Respuestas a la pregunta(1)

Su respuesta a la pregunta