Mínimos Cuadrados No Lineales

Para ajustar un modelo m\left( x,t\right) de n parámetros a un conjunto de m datos \left( t_{i},y_{i}\right) , definimos f\left( x\right) =\frac{1}{2}R\left( x\right) ^{T}R\left( x\right) con R:\mathbb{R}^{n}\rightarrow \mathbb{R}^{m} y R\left( x\right) no lineal con r_{i}\left( x\right) =m\left( x,t_{i}\right) -y_{i} cómo i-ésima componente.

El problema finalmente es:

Minimizar f\left( x\right) =\frac{1}{2}R\left( x\right) ^{T}R\left( x\right) =\frac{1}{2}\sum\limits_{i=1}^{m}r_{i}^{2}\left( x\right)

Esto es equivalente a minimizar la suma de los cuadrados de las diferencias entre el modelo y los datos.

f\left( x\right) es una función semi-paraboidal en el sentido de que es la suma de parábolas y por tanto es minimizable.

Para hallar x tal que f\left( x\right) se minimiza, hallamos las soluciones del gradiente de f\left( x\right) , de este modo, definiendo J\left( x\right) \in \mathbb{R} ^{m\times n} como el jacobiano de R\left( x\right) ,J\left( x\right) _{ij}=\frac{\partial r_{i}\left( x\right) }{\partial x_{j}} tenemos las 2 primeras derivadas de f\left( x\right) :

\nabla f\left( x\right) =\sum\limits_{i=1}^{m}r_{i}\left( x\right) \cdot \nabla r_{i}\left( x\right) =J\left( x\right) ^{T}R\left( x\right)

\nabla ^{2}f\left( x\right) =\sum\limits_{i=1}^{m}\left( \nabla r_{i}\left( x\right) \cdot \nabla r_{i}\left( x\right) ^{T}+r_{i}\left( x\right) \cdot \nabla ^{2}r_{i}\left( x\right) \right) =J\left( x\right) ^{T}J\left( x\right) +S\left( x\right)

Siendo S\left( x\right) =\sum\limits_{i=1}^{m}r_{i}\left( x\right) \cdot \nabla ^{2}r_{i}\left( x\right)

Aplicando el método de Newton para resolver \nabla f\left( x\right) tenemos:

x_{n+1}=x_{n}-\left( \nabla ^{2}f\left( x_{n}\right) \right) ^{-1}\nabla f\left( x_{n}\right)

x_{n+1}=x_{n}-\left( J\left( x_{n}\right) ^{T}J\left( x_{n}\right) +S\left( x_{n}\right) \right) ^{-1}J\left( x_{n}\right) ^{T}R\left( x_{n}\right)


EJEMPLO

Ajustar el modelo y\left( t\right) =e^{tx} a los siguientes datos:

\begin{tabular}{|l|l|l|l|} \hline t & 1 & 2 & 3 \\ \hline y & 2 & 4 & 8 \\ \hline \end{tabular}

De aquí, m\left( x,t_{i}\right) =e^{t_{i}x} entonces:

R=r_{i}\left( x\right) =m\left( x,t_{i}\right) -y_{0}\left( t_{i}\right) =\left( \begin{array}{c} e^{x}-2 \\ e^{2x}-4 \\ e^{3x}-8\end{array}\right)

Ahora f\left( x\right) =\frac{1}{2}R\left( x\right) ^{T}R\left( x\right) =\frac{1}{2}\left[ \left( e^{x}-2\right) ^{2}+\left( e^{2x}-4\right) ^{2}+\left( e^{3x}-8\right) ^{2}\right]

De este modo:

\nabla f\left( x\right) =e^{x}\left( e^{x}-2\right) +2e^{2x}\left( e^{2x}-4\right) +3e^{3x}\left( e^{3x}-8\right)

\nabla ^{2}f\left( x\right) =\left( e^{x}\right) ^{2}+\left( 2e^{2x}\right) ^{2}+\left( 3e^{3x}\right) ^{2}+e^{x}\left( e^{x}-2\right) +4e^{2x}\left( e^{2x}-4\right) +9e^{3x}\left( e^{3x}-8\right)

Finalmente aplicando el método de Newton con:

x_{n+1}=x_{n}-\left( \nabla ^{2}f\left( x_{n}\right) \right) ^{-1}\nabla f\left( x_{n}\right)

Con x_{0}=1 y 8 decimales tenemos:

\begin{tabular}{|l|l|} \hline n & x \\ \hline 0 & 1.00000000 \\ \hline 1 & 0.87299182 \\ \hline 2 & 0.77360001 \\ \hline 3 & 0.71421136 \\ \hline 4 & 0.69491616 \\ \hline 5 & 0.69316064 \\ \hline 6 & 0.69314718 \\ \hline 7 & 0.69314718 \\ \hline \end{tabular}

Obtenemos una convergencia en 7 iteraciones con x_{7}=\ln 2\approx 0.69314718

De modo que y\left( t\right) =e^{t\ln 2}=2^{t} la cual se ajusta a los datos.


Post #10mycomplexsoul

7 pensamientos en “Mínimos Cuadrados No Lineales

  1. Perdona si me atrevo a preguntarte. Estoy estimando un modelo a través de una regresión de tipo logística. El caso es que mi variable dependiente son directamente las probabilidades, por lo que la estimación de los parametros se me hace muy costosa en el spss. Porque la regresión logística la hace a partir de la categórica y no de las probabilidades. Podrías echarme una mano? Gracias

  2. HOla el ejemplo me aclara muy bien todo el rollo teorico del inicio. Sin embargo, tengo mis dudas sobre como seria para dos parametros. Por ejemplo con un modelo como:

    y= exp(a(t-b))

    donde los parametros a estimar son a y b?
    se puede aplicar el metodo?
    como quedaria R?

    racias un saludo.

    • Respondiendo a tus dudas: si, se puede aplicar y en este caso R quedaría como sigue:

      R(t)=\left( \begin{array}{c} e^{a_{1}\left(t-b_{1}\right)}-y_{1} \\ e^{a_{2}\left(t-b_{2}\right)}-y_{2} \\ ...\end{array}\right)

      Y así dependiendo cuantos datos tengas que ajustar.
      Desde luego operar con todas esas variables puede complicarse a nivel del Jacobiano, mi recomendación sería “linealizar” el modelo que sugieres, es decir hacer un cambio de variables como sigue:

      y = e^{a\left(t-b\right)}
      \ln y = a*t - a*b

      Renombrando c= -a*b y z = \ln y tenemos z = a*t + c un modelo lineal, que evidentemente es más fácil de manejar.

      Con cualquier camino deberías llegar a los mismos resultados. Saludos.

      • Hola,
        Si bien el modelo linealizado es más fácil de tratar los resultados obtenidos pueden no ser correctos, en efecto los errores de los datos medidos pueden hacer que uno de los estimadores de los parámetros buscados no converga al valor verdadero. Esto es, particularmente notable cuando en los valores medidos de la variable t tienes errores sistemáticos.

  3. Hola,

    Disculpa. Estoy programando en MatLab este método, pero tengo una duda. Necesito optimizar una función con dos parámetros. El problema que encuentro es que no sé como obtener S(x). En este caso debería resultarme una matriz de 2×2, pues JT*J es de 2×2 así que la suma JT*J+S debe ser de las mismas dimensiones. Además que debe ser cuadrada para poder obtener la inversa cuando se obtenga el valor de x iterativamente. También debo realizar la obtención para un caso de 3 variables, en este caso asumo que S(x) debe ser de 3×3. Pero igualmente no sé como obtener este S(x).

    Puedes ayudarme, por favor.

    • Hola, supongo que el problema real está en la segunda derivada de r_{i}\left ( x \right ) que resulta ser una matriz nxn (2×2 en el caso que describes), tomando en cuenta que:

      \triangledown r_{i} \left ( x \right )_{j} = \left [ \frac{\partial }{\partial x_{j}} \left ( r_{i}\left ( x \right ) \right ) \right ]_{j}

      Tenemos entonces:

      \triangledown ^{2} r_{i} \left ( x \right )_{jk} = \left [ \frac{\partial }{\partial x_{k}} \frac{\partial }{\partial x_{j}} \left ( r_{i}\left ( x \right ) \right ) \right ]_{jk}

      Es decir, cada fila de la matriz es la derivada con respecto a la k-ésima variable del gradiente de r_{i}\left ( x \right ). Con esto y la definición de S(x) debería ser posible calcularlo.

      Espero sea de utilidad. Saludos.

      • Muchas Gracias,

        Si ayer estuve muy cansado, hoy desperté y recorde que existen las derivadas cruzadas (jajaja). Estoy realizando la implementación del método en MatLab. Cuando lo termine quisiera compartirlo, ahí te aviso. Es para mi tesis que es sobre la aplicación de métodos de optimización en el modelamiento del comportamiento de material viscoelástico (hormigón asfáltico). Sígueme en twitter en @topicster. Sería genial poder intercambiar conocimientos.

        Adiós.

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s